1 |
mocchiut |
1.1 |
/** |
2 |
|
|
* \file TrkLevel2.cpp |
3 |
|
|
* \author Elena Vannuccini |
4 |
|
|
*/ |
5 |
|
|
#include <TrkLevel2.h> |
6 |
|
|
#include <iostream> |
7 |
pam-fi |
1.14 |
#include <math.h> |
8 |
mocchiut |
1.1 |
using namespace std; |
9 |
|
|
//...................................... |
10 |
|
|
// F77 routines |
11 |
|
|
//...................................... |
12 |
|
|
extern "C" { |
13 |
|
|
void dotrack_(int*, double*, double*, double*, double*, int*); |
14 |
pam-fi |
1.2 |
void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*); |
15 |
pam-fi |
1.16 |
// int readb_(const char*); |
16 |
pam-fi |
1.26 |
// int readb_(); |
17 |
|
|
void mini2_(int*,int*,int*); |
18 |
|
|
void guess_(); |
19 |
|
|
void gufld_(float*, float*); |
20 |
mocchiut |
1.1 |
} |
21 |
pam-fi |
1.26 |
|
22 |
mocchiut |
1.1 |
//-------------------------------------- |
23 |
|
|
// |
24 |
|
|
// |
25 |
|
|
//-------------------------------------- |
26 |
|
|
TrkTrack::TrkTrack(){ |
27 |
pam-fi |
1.21 |
// cout << "TrkTrack::TrkTrack()" << endl; |
28 |
pam-fi |
1.3 |
seqno = -1; |
29 |
|
|
image = -1; |
30 |
mocchiut |
1.1 |
chi2 = 0; |
31 |
pam-fi |
1.14 |
nstep = 0; |
32 |
|
|
for(int it1=0;it1<5;it1++){ |
33 |
|
|
al[it1] = 0; |
34 |
|
|
for(int it2=0;it2<5;it2++)coval[it1][it2] = 0; |
35 |
mocchiut |
1.1 |
}; |
36 |
|
|
for(int ip=0;ip<6;ip++){ |
37 |
pam-fi |
1.21 |
xgood[ip] = 0; |
38 |
|
|
ygood[ip] = 0; |
39 |
|
|
xm[ip] = 0; |
40 |
|
|
ym[ip] = 0; |
41 |
|
|
zm[ip] = 0; |
42 |
|
|
resx[ip] = 0; |
43 |
|
|
resy[ip] = 0; |
44 |
|
|
xv[ip] = 0; |
45 |
|
|
yv[ip] = 0; |
46 |
|
|
zv[ip] = 0; |
47 |
|
|
axv[ip] = 0; |
48 |
|
|
ayv[ip] = 0; |
49 |
|
|
dedx_x[ip] = 0; |
50 |
|
|
dedx_y[ip] = 0; |
51 |
|
|
}; |
52 |
|
|
clx = 0; |
53 |
|
|
cly = 0; |
54 |
pam-fi |
1.29 |
// clx = new TRefArray(6,0); //forse causa memory leak??? |
55 |
|
|
// cly = new TRefArray(6,0); //forse causa memory leak??? |
56 |
mocchiut |
1.1 |
}; |
57 |
|
|
//-------------------------------------- |
58 |
|
|
// |
59 |
|
|
// |
60 |
|
|
//-------------------------------------- |
61 |
|
|
TrkTrack::TrkTrack(const TrkTrack& t){ |
62 |
pam-fi |
1.3 |
seqno = t.seqno; |
63 |
mocchiut |
1.1 |
image = t.image; |
64 |
|
|
chi2 = t.chi2; |
65 |
pam-fi |
1.14 |
nstep = t.nstep; |
66 |
|
|
for(int it1=0;it1<5;it1++){ |
67 |
|
|
al[it1] = t.al[it1]; |
68 |
|
|
for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2]; |
69 |
mocchiut |
1.1 |
}; |
70 |
|
|
for(int ip=0;ip<6;ip++){ |
71 |
pam-fi |
1.21 |
xgood[ip] = t.xgood[ip]; |
72 |
|
|
ygood[ip] = t.ygood[ip]; |
73 |
|
|
xm[ip] = t.xm[ip]; |
74 |
|
|
ym[ip] = t.ym[ip]; |
75 |
|
|
zm[ip] = t.zm[ip]; |
76 |
|
|
resx[ip] = t.resx[ip]; |
77 |
|
|
resy[ip] = t.resy[ip]; |
78 |
|
|
xv[ip] = t.xv[ip]; |
79 |
|
|
yv[ip] = t.yv[ip]; |
80 |
|
|
zv[ip] = t.zv[ip]; |
81 |
|
|
axv[ip] = t.axv[ip]; |
82 |
|
|
ayv[ip] = t.ayv[ip]; |
83 |
|
|
dedx_x[ip] = t.dedx_x[ip]; |
84 |
|
|
dedx_y[ip] = t.dedx_y[ip]; |
85 |
|
|
}; |
86 |
pam-fi |
1.22 |
clx = 0; |
87 |
|
|
cly = 0; |
88 |
|
|
if(t.clx)clx = new TRefArray(*(t.clx)); |
89 |
|
|
if(t.cly)cly = new TRefArray(*(t.cly)); |
90 |
pam-fi |
1.9 |
|
91 |
mocchiut |
1.1 |
}; |
92 |
|
|
//-------------------------------------- |
93 |
|
|
// |
94 |
|
|
// |
95 |
|
|
//-------------------------------------- |
96 |
pam-fi |
1.21 |
void TrkTrack::Copy(TrkTrack& t){ |
97 |
|
|
|
98 |
|
|
t.seqno = seqno; |
99 |
|
|
t.image = image; |
100 |
|
|
t.chi2 = chi2; |
101 |
|
|
t.nstep = nstep; |
102 |
|
|
for(int it1=0;it1<5;it1++){ |
103 |
|
|
t.al[it1] = al[it1]; |
104 |
|
|
for(int it2=0;it2<5;it2++)t.coval[it1][it2] = coval[it1][it2]; |
105 |
|
|
}; |
106 |
|
|
for(int ip=0;ip<6;ip++){ |
107 |
|
|
t.xgood[ip] = xgood[ip]; |
108 |
|
|
t.ygood[ip] = ygood[ip]; |
109 |
|
|
t.xm[ip] = xm[ip]; |
110 |
|
|
t.ym[ip] = ym[ip]; |
111 |
|
|
t.zm[ip] = zm[ip]; |
112 |
|
|
t.resx[ip] = resx[ip]; |
113 |
|
|
t.resy[ip] = resy[ip]; |
114 |
|
|
t.xv[ip] = xv[ip]; |
115 |
|
|
t.yv[ip] = yv[ip]; |
116 |
|
|
t.zv[ip] = zv[ip]; |
117 |
|
|
t.axv[ip] = axv[ip]; |
118 |
|
|
t.ayv[ip] = ayv[ip]; |
119 |
|
|
t.dedx_x[ip] = dedx_x[ip]; |
120 |
|
|
t.dedx_y[ip] = dedx_y[ip]; |
121 |
|
|
|
122 |
|
|
}; |
123 |
|
|
|
124 |
|
|
}; |
125 |
|
|
//-------------------------------------- |
126 |
|
|
// |
127 |
|
|
// |
128 |
|
|
//-------------------------------------- |
129 |
mocchiut |
1.1 |
/** |
130 |
|
|
* Evaluates the trajectory in the apparatus associated to the track. |
131 |
|
|
* It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling TrkLevel2::LoadField() ), otherwise an error message is returned. |
132 |
|
|
* @param t pointer to an object of the class Trajectory, |
133 |
|
|
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
134 |
|
|
* @return error flag. |
135 |
|
|
*/ |
136 |
|
|
int TrkTrack::DoTrack(Trajectory* t){ |
137 |
|
|
|
138 |
|
|
double *dxout = new double[t->npoint]; |
139 |
|
|
double *dyout = new double[t->npoint]; |
140 |
|
|
double *dzin = new double[t->npoint]; |
141 |
|
|
double dal[5]; |
142 |
|
|
|
143 |
|
|
int ifail = 0; |
144 |
|
|
|
145 |
|
|
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
146 |
|
|
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; |
147 |
|
|
|
148 |
pam-fi |
1.26 |
TrkParams::Load(1); |
149 |
|
|
if( !TrkParams::IsLoaded(1) ){ |
150 |
|
|
cout << "int TrkTrack::DoTrack(Trajectory* t) --- ERROR --- m.field not loaded"<<endl; |
151 |
|
|
return 0; |
152 |
|
|
} |
153 |
mocchiut |
1.1 |
dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail); |
154 |
|
|
|
155 |
|
|
for (int i=0; i<t->npoint; i++){ |
156 |
|
|
t->x[i] = (float)*dxout++; |
157 |
|
|
t->y[i] = (float)*dyout++; |
158 |
|
|
} |
159 |
|
|
|
160 |
|
|
// delete [] dxout; |
161 |
|
|
// delete [] dyout; |
162 |
|
|
// delete [] dzin; |
163 |
|
|
|
164 |
|
|
return ifail; |
165 |
|
|
}; |
166 |
|
|
//-------------------------------------- |
167 |
|
|
// |
168 |
|
|
// |
169 |
|
|
//-------------------------------------- |
170 |
pam-fi |
1.2 |
/** |
171 |
|
|
* Evaluates the trajectory in the apparatus associated to the track. |
172 |
|
|
* It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling TrkLevel2::LoadField() ), otherwise an error message is returned. |
173 |
|
|
* @param t pointer to an object of the class Trajectory, |
174 |
|
|
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
175 |
|
|
* @return error flag. |
176 |
|
|
*/ |
177 |
|
|
int TrkTrack::DoTrack2(Trajectory* t){ |
178 |
|
|
|
179 |
|
|
double *dxout = new double[t->npoint]; |
180 |
|
|
double *dyout = new double[t->npoint]; |
181 |
|
|
double *dthxout = new double[t->npoint]; |
182 |
|
|
double *dthyout = new double[t->npoint]; |
183 |
|
|
double *dtlout = new double[t->npoint]; |
184 |
|
|
double *dzin = new double[t->npoint]; |
185 |
|
|
double dal[5]; |
186 |
|
|
|
187 |
|
|
int ifail = 0; |
188 |
|
|
|
189 |
|
|
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
190 |
|
|
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; |
191 |
|
|
|
192 |
pam-fi |
1.26 |
TrkParams::Load(1); |
193 |
|
|
if( !TrkParams::IsLoaded(1) ){ |
194 |
|
|
cout << "int TrkTrack::DoTrack2(Trajectory* t) --- ERROR --- m.field not loaded"<<endl; |
195 |
|
|
return 0; |
196 |
|
|
} |
197 |
pam-fi |
1.2 |
dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
198 |
|
|
|
199 |
|
|
for (int i=0; i<t->npoint; i++){ |
200 |
|
|
t->x[i] = (float)*dxout++; |
201 |
|
|
t->y[i] = (float)*dyout++; |
202 |
|
|
t->thx[i] = (float)*dthxout++; |
203 |
|
|
t->thy[i] = (float)*dthyout++; |
204 |
|
|
t->tl[i] = (float)*dtlout++; |
205 |
|
|
} |
206 |
|
|
|
207 |
|
|
// delete [] dxout; |
208 |
|
|
// delete [] dyout; |
209 |
|
|
// delete [] dzin; |
210 |
|
|
|
211 |
|
|
return ifail; |
212 |
|
|
}; |
213 |
|
|
//-------------------------------------- |
214 |
|
|
// |
215 |
|
|
// |
216 |
|
|
//-------------------------------------- |
217 |
mocchiut |
1.1 |
//float TrkTrack::BdL(){ |
218 |
|
|
//}; |
219 |
|
|
//-------------------------------------- |
220 |
|
|
// |
221 |
|
|
// |
222 |
|
|
//-------------------------------------- |
223 |
|
|
Float_t TrkTrack::GetRigidity(){ |
224 |
|
|
Float_t rig=0; |
225 |
|
|
if(chi2>0)rig=1./al[4]; |
226 |
|
|
if(rig<0) rig=-rig; |
227 |
|
|
return rig; |
228 |
|
|
}; |
229 |
|
|
// |
230 |
|
|
Float_t TrkTrack::GetDeflection(){ |
231 |
|
|
Float_t def=0; |
232 |
|
|
if(chi2>0)def=al[4]; |
233 |
|
|
return def; |
234 |
|
|
}; |
235 |
|
|
// |
236 |
|
|
Float_t TrkTrack::GetDEDX(){ |
237 |
|
|
Float_t dedx=0; |
238 |
pam-fi |
1.30 |
// for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*xgood[i]+dedx_y[i]*ygood[i]; |
239 |
|
|
for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*XGood(i)+dedx_y[i]*YGood(i); |
240 |
mocchiut |
1.1 |
dedx = dedx/(this->GetNX()+this->GetNY()); |
241 |
|
|
return dedx; |
242 |
|
|
}; |
243 |
|
|
|
244 |
|
|
//-------------------------------------- |
245 |
|
|
// |
246 |
|
|
// |
247 |
|
|
//-------------------------------------- |
248 |
|
|
void TrkTrack::Dump(){ |
249 |
|
|
cout << endl << "========== Track " ; |
250 |
pam-fi |
1.13 |
cout << endl << "seq. n. : "<< seqno; |
251 |
|
|
cout << endl << "image n. : "<< image; |
252 |
|
|
cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " "; |
253 |
mocchiut |
1.1 |
cout << endl << "chi^2 : "<< chi2; |
254 |
pam-fi |
1.13 |
cout << endl << "n.step : "<< nstep; |
255 |
pam-fi |
1.30 |
cout << endl << "xgood : "; for(int i=0; i<6; i++)cout << XGood(i) ; |
256 |
|
|
cout << endl << "ygood : "; for(int i=0; i<6; i++)cout << YGood(i) ; |
257 |
mocchiut |
1.1 |
cout << endl << "xm : "; for(int i=0; i<6; i++)cout << xm[i] << " "; |
258 |
|
|
cout << endl << "ym : "; for(int i=0; i<6; i++)cout << ym[i] << " "; |
259 |
|
|
cout << endl << "zm : "; for(int i=0; i<6; i++)cout << zm[i] << " "; |
260 |
pam-fi |
1.13 |
cout << endl << "xv : "; for(int i=0; i<6; i++)cout << xv[i] << " "; |
261 |
|
|
cout << endl << "yv : "; for(int i=0; i<6; i++)cout << yv[i] << " "; |
262 |
|
|
cout << endl << "zv : "; for(int i=0; i<6; i++)cout << zv[i] << " "; |
263 |
|
|
cout << endl << "resx : "; for(int i=0; i<6; i++)cout << resx[i] << " "; |
264 |
|
|
cout << endl << "resy : "; for(int i=0; i<6; i++)cout << resy[i] << " "; |
265 |
|
|
cout << endl << "coval : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" "; |
266 |
|
|
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[1][i]<<" "; |
267 |
|
|
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[2][i]<<" "; |
268 |
|
|
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[3][i]<<" "; |
269 |
|
|
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[4][i]<<" "; |
270 |
mocchiut |
1.1 |
cout << endl << "dedx_x : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " "; |
271 |
|
|
cout << endl << "dedx_y : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " "; |
272 |
pam-fi |
1.14 |
cout << endl; |
273 |
mocchiut |
1.1 |
} |
274 |
pam-fi |
1.13 |
/** |
275 |
|
|
* Set the TrkTrack position measurements |
276 |
|
|
*/ |
277 |
|
|
void TrkTrack::SetMeasure(double *xmeas, double *ymeas, double *zmeas){ |
278 |
|
|
for(int i=0; i<6; i++) xm[i]=*xmeas++; |
279 |
|
|
for(int i=0; i<6; i++) ym[i]=*ymeas++; |
280 |
|
|
for(int i=0; i<6; i++) zm[i]=*zmeas++; |
281 |
|
|
} |
282 |
|
|
/** |
283 |
|
|
* Set the TrkTrack position resolution |
284 |
|
|
*/ |
285 |
|
|
void TrkTrack::SetResolution(double *rx, double *ry){ |
286 |
|
|
for(int i=0; i<6; i++) resx[i]=*rx++; |
287 |
|
|
for(int i=0; i<6; i++) resy[i]=*ry++; |
288 |
|
|
} |
289 |
|
|
/** |
290 |
|
|
* Set the TrkTrack good measurement |
291 |
|
|
*/ |
292 |
|
|
void TrkTrack::SetGood(int *xg, int *yg){ |
293 |
pam-fi |
1.30 |
// NB! si perdera` l'informazione sul numero del cluster |
294 |
pam-fi |
1.13 |
for(int i=0; i<6; i++) xgood[i]=*xg++; |
295 |
|
|
for(int i=0; i<6; i++) ygood[i]=*yg++; |
296 |
|
|
} |
297 |
|
|
|
298 |
|
|
/** |
299 |
|
|
* Load the magnetic field |
300 |
|
|
*/ |
301 |
pam-fi |
1.16 |
void TrkTrack::LoadField(TString path){ |
302 |
|
|
|
303 |
pam-fi |
1.26 |
// strcpy(path_.path,path.Data()); |
304 |
|
|
// path_.pathlen = path.Length(); |
305 |
|
|
// path_.error = 0; |
306 |
|
|
// readb_(); |
307 |
|
|
|
308 |
|
|
TrkParams::Set(path,1); |
309 |
pam-fi |
1.28 |
TrkParams::Load(1); |
310 |
pam-fi |
1.16 |
|
311 |
pam-fi |
1.13 |
}; |
312 |
pam-fi |
1.19 |
|
313 |
pam-fi |
1.25 |
|
314 |
pam-fi |
1.19 |
/** |
315 |
|
|
* Method to fill minimization-routine common |
316 |
|
|
*/ |
317 |
|
|
void TrkTrack::FillMiniStruct(cMini2track& track){ |
318 |
|
|
|
319 |
|
|
for(int i=0; i<6; i++){ |
320 |
|
|
|
321 |
pam-fi |
1.30 |
// track.xgood[i]=xgood[i]; |
322 |
|
|
// track.ygood[i]=ygood[i]; |
323 |
|
|
track.xgood[i]=XGood(i); |
324 |
|
|
track.ygood[i]=YGood(i); |
325 |
pam-fi |
1.19 |
|
326 |
|
|
track.xm[i]=xm[i]; |
327 |
|
|
track.ym[i]=ym[i]; |
328 |
|
|
track.zm[i]=zm[i]; |
329 |
|
|
|
330 |
|
|
// --- temporaneo ---------------------------- |
331 |
|
|
// andrebbe inserita la dimensione del sensore |
332 |
|
|
float segment = 100.; |
333 |
|
|
track.xm_a[i]=xm[i]; |
334 |
|
|
track.xm_b[i]=xm[i]; |
335 |
|
|
track.ym_a[i]=ym[i]; |
336 |
|
|
track.ym_b[i]=ym[i]; |
337 |
pam-fi |
1.30 |
if( XGood(i) && !YGood(i) ){ |
338 |
pam-fi |
1.19 |
track.ym_a[i] = track.ym_a[i]+segment; |
339 |
|
|
track.ym_b[i] = track.ym_b[i]-segment; |
340 |
pam-fi |
1.30 |
}else if( !XGood(i) && YGood(i)){ |
341 |
pam-fi |
1.19 |
track.xm_a[i] = track.xm_a[i]+segment; |
342 |
|
|
track.xm_b[i] = track.xm_b[i]-segment; |
343 |
|
|
} |
344 |
|
|
// --- temporaneo ---------------------------- |
345 |
|
|
|
346 |
|
|
track.resx[i]=resx[i]; |
347 |
|
|
track.resy[i]=resy[i]; |
348 |
|
|
} |
349 |
|
|
|
350 |
|
|
for(int i=0; i<5; i++) track.al[i]=al[i]; |
351 |
|
|
track.zini = 23.5; |
352 |
|
|
// ZINI = 23.5 !!! it should be the same parameter in all codes |
353 |
|
|
|
354 |
|
|
} |
355 |
|
|
/** |
356 |
|
|
* Method to set values from minimization-routine common |
357 |
|
|
*/ |
358 |
|
|
void TrkTrack::SetFromMiniStruct(cMini2track *track){ |
359 |
|
|
|
360 |
|
|
for(int i=0; i<5; i++) { |
361 |
|
|
al[i]=track->al[i]; |
362 |
|
|
for(int j=0; j<5; j++) coval[i][j]=track->cov[i][j]; |
363 |
|
|
} |
364 |
|
|
chi2 = track->chi2; |
365 |
|
|
nstep = track->nstep; |
366 |
|
|
for(int i=0; i<6; i++){ |
367 |
|
|
xv[i] = track->xv[i]; |
368 |
|
|
yv[i] = track->yv[i]; |
369 |
|
|
zv[i] = track->zv[i]; |
370 |
|
|
xm[i] = track->xm[i]; |
371 |
|
|
ym[i] = track->ym[i]; |
372 |
|
|
zm[i] = track->zm[i]; |
373 |
|
|
axv[i] = track->axv[i]; |
374 |
|
|
ayv[i] = track->ayv[i]; |
375 |
|
|
} |
376 |
|
|
|
377 |
|
|
} |
378 |
pam-fi |
1.13 |
/** |
379 |
|
|
* Tracking method. It calls F77 mini routine. |
380 |
|
|
*/ |
381 |
|
|
void TrkTrack::Fit(double pfixed, int& fail, int iprint){ |
382 |
pam-fi |
1.15 |
|
383 |
|
|
float al_ini[] = {0.,0.,0.,0.,0.}; |
384 |
|
|
|
385 |
pam-fi |
1.13 |
extern cMini2track track_; |
386 |
|
|
fail = 0; |
387 |
pam-fi |
1.19 |
FillMiniStruct(track_); |
388 |
pam-fi |
1.15 |
|
389 |
pam-fi |
1.19 |
// if fit variables have been reset, evaluate the initial guess |
390 |
pam-fi |
1.15 |
if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_(); |
391 |
|
|
|
392 |
pam-fi |
1.16 |
// --------------------- free momentum |
393 |
pam-fi |
1.14 |
if(pfixed==0.) { |
394 |
pam-fi |
1.19 |
track_.pfixed=0.; |
395 |
pam-fi |
1.14 |
} |
396 |
pam-fi |
1.16 |
// --------------------- fixed momentum |
397 |
pam-fi |
1.13 |
if(pfixed!=0.) { |
398 |
pam-fi |
1.19 |
al[4]=1./pfixed; |
399 |
|
|
track_.pfixed=pfixed; |
400 |
pam-fi |
1.13 |
} |
401 |
pam-fi |
1.15 |
|
402 |
pam-fi |
1.19 |
// store temporarily the initial guess |
403 |
pam-fi |
1.15 |
for(int i=0; i<5; i++) al_ini[i]=track_.al[i]; |
404 |
|
|
|
405 |
pam-fi |
1.19 |
// ------------------------------------------ |
406 |
|
|
// call mini routine |
407 |
pam-fi |
1.26 |
TrkParams::Load(1); |
408 |
|
|
if( !TrkParams::IsLoaded(1) ){ |
409 |
|
|
cout << "void TrkTrack::Fit(double pfixed, int& fail, int iprint) --- ERROR --- m.field not loaded"<<endl; |
410 |
|
|
return; |
411 |
|
|
} |
412 |
pam-fi |
1.13 |
int istep=0; |
413 |
|
|
int ifail=0; |
414 |
|
|
mini2_(&istep,&ifail, &iprint); |
415 |
|
|
if(ifail!=0) { |
416 |
pam-fi |
1.19 |
if(iprint)cout << "ERROR: ifail= " << ifail << endl; |
417 |
pam-fi |
1.13 |
fail = 1; |
418 |
|
|
} |
419 |
pam-fi |
1.19 |
// ------------------------------------------ |
420 |
pam-fi |
1.15 |
|
421 |
pam-fi |
1.19 |
SetFromMiniStruct(&track_); |
422 |
pam-fi |
1.13 |
// cout << endl << "eta ===> " << track_.al[4] << endl; |
423 |
|
|
|
424 |
pam-fi |
1.19 |
// for(int i=0; i<5; i++) al[i]=track_.al[i]; |
425 |
|
|
// chi2=track_.chi2; |
426 |
|
|
// nstep=track_.nstep; |
427 |
|
|
// for(int i=0; i<6; i++) xv[i]=track_.xv[i]; |
428 |
|
|
// for(int i=0; i<6; i++) yv[i]=track_.yv[i]; |
429 |
|
|
// for(int i=0; i<6; i++) zv[i]=track_.zv[i]; |
430 |
|
|
// for(int i=0; i<6; i++) axv[i]=track_.axv[i]; |
431 |
|
|
// for(int i=0; i<6; i++) ayv[i]=track_.ayv[i]; |
432 |
|
|
// for(int i=0; i<5; i++) { |
433 |
|
|
// for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j]; |
434 |
|
|
// } |
435 |
pam-fi |
1.15 |
|
436 |
pam-fi |
1.20 |
if(fail){ |
437 |
|
|
if(iprint)cout << " >>>> fit failed >>>> drawing initial par"<<endl; |
438 |
|
|
for(int i=0; i<5; i++) al[i]=al_ini[i]; |
439 |
|
|
} |
440 |
pam-fi |
1.15 |
|
441 |
pam-fi |
1.13 |
}; |
442 |
|
|
/* |
443 |
pam-fi |
1.15 |
* Reset the fit parameters |
444 |
pam-fi |
1.13 |
*/ |
445 |
|
|
void TrkTrack::FitReset(){ |
446 |
|
|
for(int i=0; i<5; i++) al[i]=-9999.; |
447 |
|
|
chi2=0.; |
448 |
|
|
nstep=0; |
449 |
|
|
for(int i=0; i<6; i++) xv[i]=0.; |
450 |
|
|
for(int i=0; i<6; i++) yv[i]=0.; |
451 |
|
|
for(int i=0; i<6; i++) zv[i]=0.; |
452 |
|
|
for(int i=0; i<6; i++) axv[i]=0.; |
453 |
|
|
for(int i=0; i<6; i++) ayv[i]=0.; |
454 |
|
|
for(int i=0; i<5; i++) { |
455 |
|
|
for(int j=0; j<5; j++) coval[i][j]=0.; |
456 |
|
|
} |
457 |
|
|
} |
458 |
pam-fi |
1.31 |
/* |
459 |
|
|
* Set the tracking mode |
460 |
|
|
*/ |
461 |
pam-fi |
1.27 |
void TrkTrack::SetTrackingMode(int trackmode){ |
462 |
|
|
extern cMini2track track_; |
463 |
|
|
track_.trackmode = trackmode; |
464 |
|
|
} |
465 |
pam-fi |
1.31 |
/* |
466 |
|
|
* Set the factor scale for tracking precision |
467 |
|
|
*/ |
468 |
|
|
void TrkTrack::SetPrecisionFactor(double fact){ |
469 |
|
|
extern cMini2track track_; |
470 |
|
|
track_.fact = fact; |
471 |
|
|
} |
472 |
|
|
/* |
473 |
|
|
* Set the factor scale for tracking precision |
474 |
|
|
*/ |
475 |
|
|
void TrkTrack::SetStepMin(int istepmin){ |
476 |
|
|
extern cMini2track track_; |
477 |
|
|
track_.istepmin = istepmin; |
478 |
|
|
} |
479 |
pam-fi |
1.13 |
|
480 |
pam-fi |
1.29 |
|
481 |
|
|
/* |
482 |
|
|
* Method to retrieve the X-view clusters associated to the track. |
483 |
|
|
* @param ip Tracker plane (0-5) |
484 |
|
|
*/ |
485 |
|
|
TrkCluster *TrkTrack::GetClusterX(int ip){ |
486 |
|
|
cout << " TrkCluster *TrkTrack::GetClusterX(int ip) -- momentaneamente fuori servizio --"<< endl; |
487 |
|
|
if(!clx)return NULL; |
488 |
|
|
TrkCluster *pt = (TrkCluster*)(clx->At(ip)); |
489 |
|
|
return pt; |
490 |
|
|
}; |
491 |
|
|
/* |
492 |
|
|
* Method to retrieve the Y-view clusters associated to the track. |
493 |
|
|
* @param ip Tracker plane (0-5) |
494 |
|
|
*/ |
495 |
|
|
TrkCluster *TrkTrack::GetClusterY(int ip){ |
496 |
|
|
cout << " TrkCluster *TrkTrack::GetClusterY(int ip) -- momentaneamente fuori servizio --"<< endl; |
497 |
|
|
if(!cly)return NULL; |
498 |
|
|
TrkCluster *pt = (TrkCluster*)(cly->At(ip)); |
499 |
|
|
return pt; |
500 |
|
|
}; |
501 |
|
|
|
502 |
|
|
|
503 |
mocchiut |
1.1 |
//-------------------------------------- |
504 |
|
|
// |
505 |
|
|
// |
506 |
|
|
//-------------------------------------- |
507 |
pam-fi |
1.10 |
void TrkTrack::Clear(){ |
508 |
pam-fi |
1.21 |
// cout << "TrkTrack::Clear()"<<endl; |
509 |
|
|
seqno = -1; |
510 |
|
|
image = -1; |
511 |
|
|
chi2 = 0; |
512 |
|
|
nstep = 0; |
513 |
|
|
for(int it1=0;it1<5;it1++){ |
514 |
|
|
al[it1] = 0; |
515 |
|
|
for(int it2=0;it2<5;it2++)coval[it1][it2] = 0; |
516 |
|
|
}; |
517 |
|
|
for(int ip=0;ip<6;ip++){ |
518 |
|
|
xgood[ip] = 0; |
519 |
|
|
ygood[ip] = 0; |
520 |
|
|
xm[ip] = 0; |
521 |
|
|
ym[ip] = 0; |
522 |
|
|
zm[ip] = 0; |
523 |
|
|
resx[ip] = 0; |
524 |
|
|
resy[ip] = 0; |
525 |
|
|
xv[ip] = 0; |
526 |
|
|
yv[ip] = 0; |
527 |
|
|
zv[ip] = 0; |
528 |
|
|
axv[ip] = 0; |
529 |
|
|
ayv[ip] = 0; |
530 |
|
|
dedx_x[ip] = 0; |
531 |
|
|
dedx_y[ip] = 0; |
532 |
|
|
|
533 |
|
|
}; |
534 |
|
|
if(clx)clx->Clear(); |
535 |
|
|
if(cly)cly->Clear(); |
536 |
pam-fi |
1.10 |
}; |
537 |
|
|
//-------------------------------------- |
538 |
|
|
// |
539 |
|
|
// |
540 |
|
|
//-------------------------------------- |
541 |
pam-fi |
1.11 |
void TrkTrack::Delete(){ |
542 |
pam-fi |
1.21 |
// cout << "TrkTrack::Delete()"<<endl; |
543 |
|
|
// Clear(); |
544 |
|
|
if(clx)delete clx; |
545 |
|
|
if(cly)delete cly; |
546 |
pam-fi |
1.11 |
}; |
547 |
pam-fi |
1.21 |
//-------------------------------------- |
548 |
pam-fi |
1.11 |
// |
549 |
|
|
// |
550 |
|
|
//-------------------------------------- |
551 |
pam-fi |
1.10 |
|
552 |
|
|
//-------------------------------------- |
553 |
|
|
// |
554 |
|
|
// |
555 |
|
|
//-------------------------------------- |
556 |
mocchiut |
1.1 |
TrkSinglet::TrkSinglet(){ |
557 |
pam-fi |
1.21 |
// cout << "TrkSinglet::TrkSinglet() " << GetUniqueID()<<endl; |
558 |
mocchiut |
1.1 |
plane = 0; |
559 |
|
|
coord[0] = 0; |
560 |
|
|
coord[1] = 0; |
561 |
|
|
sgnl = 0; |
562 |
pam-fi |
1.21 |
// cls = 0; |
563 |
mocchiut |
1.1 |
}; |
564 |
|
|
//-------------------------------------- |
565 |
|
|
// |
566 |
|
|
// |
567 |
|
|
//-------------------------------------- |
568 |
|
|
TrkSinglet::TrkSinglet(const TrkSinglet& s){ |
569 |
pam-fi |
1.21 |
// cout << "TrkSinglet::TrkSinglet(const TrkSinglet& s) " << GetUniqueID()<<endl; |
570 |
mocchiut |
1.1 |
plane = s.plane; |
571 |
|
|
coord[0] = s.coord[0]; |
572 |
|
|
coord[1] = s.coord[1]; |
573 |
|
|
sgnl = s.sgnl; |
574 |
pam-fi |
1.9 |
// cls = 0;//<<<<pointer |
575 |
pam-fi |
1.21 |
cls = TRef(s.cls); |
576 |
mocchiut |
1.1 |
}; |
577 |
|
|
//-------------------------------------- |
578 |
|
|
// |
579 |
|
|
// |
580 |
|
|
//-------------------------------------- |
581 |
|
|
void TrkSinglet::Dump(){ |
582 |
|
|
int i=0; |
583 |
|
|
cout << endl << "========== Singlet " ; |
584 |
|
|
cout << endl << "plane : " << plane; |
585 |
|
|
cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++; |
586 |
|
|
cout << endl << "sgnl : " << sgnl; |
587 |
|
|
} |
588 |
|
|
//-------------------------------------- |
589 |
|
|
// |
590 |
|
|
// |
591 |
|
|
//-------------------------------------- |
592 |
pam-fi |
1.21 |
void TrkSinglet::Clear(){ |
593 |
|
|
// cout << "TrkSinglet::Clear() " << GetUniqueID()<<endl; |
594 |
|
|
// cls=0; |
595 |
|
|
plane=-1; |
596 |
|
|
coord[0]=-999; |
597 |
|
|
coord[1]=-999; |
598 |
|
|
sgnl=0; |
599 |
|
|
|
600 |
|
|
} |
601 |
|
|
//-------------------------------------- |
602 |
|
|
// |
603 |
|
|
// |
604 |
|
|
//-------------------------------------- |
605 |
mocchiut |
1.1 |
TrkLevel2::TrkLevel2(){ |
606 |
mocchiut |
1.24 |
// cout <<"TrkLevel2::TrkLevel2()"<<endl; |
607 |
mocchiut |
1.1 |
for(Int_t i=0; i<12 ; i++){ |
608 |
pam-fi |
1.10 |
good[i] = -1; |
609 |
|
|
}; |
610 |
pam-fi |
1.23 |
// okkio!! memory-leak |
611 |
|
|
// Track = new TClonesArray("TrkTrack"); |
612 |
|
|
// SingletX = new TClonesArray("TrkSinglet"); |
613 |
|
|
// SingletY = new TClonesArray("TrkSinglet"); |
614 |
pam-fi |
1.21 |
Track = 0; |
615 |
|
|
SingletX = 0; |
616 |
|
|
SingletY = 0; |
617 |
pam-fi |
1.6 |
|
618 |
mocchiut |
1.1 |
} |
619 |
|
|
//-------------------------------------- |
620 |
|
|
// |
621 |
|
|
// |
622 |
|
|
//-------------------------------------- |
623 |
pam-fi |
1.23 |
void TrkLevel2::Set(){ |
624 |
|
|
if(!Track)Track = new TClonesArray("TrkTrack"); |
625 |
|
|
if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); |
626 |
|
|
if(!SingletY)SingletY = new TClonesArray("TrkSinglet"); |
627 |
|
|
} |
628 |
|
|
//-------------------------------------- |
629 |
|
|
// |
630 |
|
|
// |
631 |
|
|
//-------------------------------------- |
632 |
mocchiut |
1.1 |
void TrkLevel2::Dump(){ |
633 |
pam-fi |
1.10 |
|
634 |
|
|
// |
635 |
mocchiut |
1.1 |
cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-"; |
636 |
pam-fi |
1.10 |
cout << endl << "good : "; for(int i=0; i<12; i++) cout << good[i]<<" "; |
637 |
mocchiut |
1.1 |
cout << endl << "ntrk() : " << this->ntrk() ; |
638 |
|
|
cout << endl << "nclsx() : " << this->nclsx(); |
639 |
|
|
cout << endl << "nclsy() : " << this->nclsy(); |
640 |
pam-fi |
1.21 |
// TClonesArray &t = *Track; |
641 |
|
|
// TClonesArray &sx = *SingletX; |
642 |
|
|
// TClonesArray &sy = *SingletY; |
643 |
|
|
// for(int i=0; i<ntrk(); i++) ((TrkTrack *)t[i])->Dump(); |
644 |
|
|
// for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); |
645 |
|
|
// for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); |
646 |
|
|
if(Track){ |
647 |
|
|
TClonesArray &t = *Track; |
648 |
|
|
for(int i=0; i<ntrk(); i++) ((TrkTrack *)t[i])->Dump(); |
649 |
|
|
} |
650 |
|
|
if(SingletX){ |
651 |
|
|
TClonesArray &sx = *SingletX; |
652 |
|
|
for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); |
653 |
|
|
} |
654 |
|
|
if(SingletY){ |
655 |
|
|
TClonesArray &sy = *SingletY; |
656 |
|
|
for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); |
657 |
|
|
} |
658 |
mocchiut |
1.1 |
} |
659 |
|
|
//-------------------------------------- |
660 |
|
|
// |
661 |
|
|
// |
662 |
|
|
//-------------------------------------- |
663 |
|
|
/** |
664 |
|
|
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). |
665 |
|
|
*/ |
666 |
pam-fi |
1.26 |
// void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){ |
667 |
pam-fi |
1.10 |
|
668 |
pam-fi |
1.26 |
// // temporary objects: |
669 |
|
|
// TrkSinglet* t_singlet = new TrkSinglet(); |
670 |
|
|
// TrkTrack* t_track = new TrkTrack(); |
671 |
|
|
|
672 |
|
|
// // **** general variables **** |
673 |
|
|
// // good2 = l2->good2; |
674 |
|
|
// for(Int_t i=0; i<12 ; i++){ |
675 |
|
|
// // crc[i] = l2->crc[i]; |
676 |
|
|
// good[i] = l2->good[i]; |
677 |
|
|
// }; |
678 |
|
|
// // *** TRACKS *** |
679 |
|
|
// if(!Track) Track = new TClonesArray("TrkTrack"); |
680 |
|
|
// TClonesArray &t = *Track; |
681 |
|
|
// for(int i=0; i<l2->ntrk; i++){ |
682 |
|
|
// t_track->seqno = i;// NBNBNBNB deve sempre essere = i |
683 |
|
|
// t_track->image = l2->image[i]-1; |
684 |
|
|
// // cout << "track "<<i<<t_track->seqno << t_track->image<<endl; |
685 |
|
|
// t_track->chi2 = l2->chi2_nt[i]; |
686 |
|
|
// t_track->nstep = l2->nstep_nt[i]; |
687 |
|
|
// for(int it1=0;it1<5;it1++){ |
688 |
|
|
// t_track->al[it1] = l2->al_nt[i][it1]; |
689 |
|
|
// for(int it2=0;it2<5;it2++) |
690 |
|
|
// t_track->coval[it1][it2] = l2->coval[i][it2][it1]; |
691 |
|
|
// }; |
692 |
|
|
// for(int ip=0;ip<6;ip++){ |
693 |
|
|
// t_track->xgood[ip] = l2->xgood_nt[i][ip]; |
694 |
|
|
// t_track->ygood[ip] = l2->ygood_nt[i][ip]; |
695 |
|
|
// t_track->xm[ip] = l2->xm_nt[i][ip]; |
696 |
|
|
// t_track->ym[ip] = l2->ym_nt[i][ip]; |
697 |
|
|
// t_track->zm[ip] = l2->zm_nt[i][ip]; |
698 |
|
|
// t_track->resx[ip] = l2->resx_nt[i][ip]; |
699 |
|
|
// t_track->resy[ip] = l2->resy_nt[i][ip]; |
700 |
|
|
// t_track->xv[ip] = l2->xv_nt[i][ip]; |
701 |
|
|
// t_track->yv[ip] = l2->yv_nt[i][ip]; |
702 |
|
|
// t_track->zv[ip] = l2->zv_nt[i][ip]; |
703 |
|
|
// t_track->axv[ip] = l2->axv_nt[i][ip]; |
704 |
|
|
// t_track->ayv[ip] = l2->ayv_nt[i][ip]; |
705 |
|
|
// t_track->dedx_x[ip] = l2->dedx_x[i][ip]; |
706 |
|
|
// t_track->dedx_y[ip] = l2->dedx_y[i][ip]; |
707 |
|
|
// // t_track->clx[ip] = 0; |
708 |
|
|
// // t_track->cly[ip] = 0; |
709 |
|
|
// }; |
710 |
|
|
// new(t[i]) TrkTrack(*t_track); |
711 |
|
|
// t_track->Clear(); |
712 |
|
|
// }; |
713 |
|
|
// // *** SINGLETS *** |
714 |
|
|
// if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); |
715 |
|
|
// TClonesArray &sx = *SingletX; |
716 |
|
|
// for(int i=0; i<l2->nclsx; i++){ |
717 |
|
|
// t_singlet->plane = l2->planex[i]; |
718 |
|
|
// t_singlet->coord[0] = l2->xs[i][0]; |
719 |
|
|
// t_singlet->coord[1] = l2->xs[i][1]; |
720 |
|
|
// t_singlet->sgnl = l2->signlxs[i]; |
721 |
|
|
// // t_singlet->cls = 0; |
722 |
|
|
// new(sx[i]) TrkSinglet(*t_singlet); |
723 |
|
|
// t_singlet->Clear(); |
724 |
|
|
// } |
725 |
|
|
// if(!SingletY)SingletY = new TClonesArray("TrkSinglet"); |
726 |
|
|
// TClonesArray &sy = *SingletY; |
727 |
|
|
// for(int i=0; i<l2->nclsy; i++){ |
728 |
|
|
// t_singlet->plane = l2->planey[i]; |
729 |
|
|
// t_singlet->coord[0] = l2->ys[i][0]; |
730 |
|
|
// t_singlet->coord[1] = l2->ys[i][1]; |
731 |
|
|
// t_singlet->sgnl = l2->signlys[i]; |
732 |
|
|
// // t_singlet->cls = 0; |
733 |
|
|
// new(sy[i]) TrkSinglet(*t_singlet); |
734 |
|
|
// t_singlet->Clear(); |
735 |
|
|
// }; |
736 |
pam-fi |
1.9 |
|
737 |
pam-fi |
1.26 |
// delete t_track; |
738 |
|
|
// delete t_singlet; |
739 |
|
|
// } |
740 |
pam-fi |
1.9 |
//-------------------------------------- |
741 |
|
|
// |
742 |
|
|
// |
743 |
|
|
//-------------------------------------- |
744 |
|
|
/** |
745 |
|
|
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). |
746 |
pam-fi |
1.29 |
* Ref to Level1 data (clusters) is also set. If l1==NULL no references are set. |
747 |
|
|
* (NB It make sense to set references only if events are stored in a tree that contains also the Level1 branch) |
748 |
pam-fi |
1.9 |
*/ |
749 |
|
|
void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){ |
750 |
|
|
|
751 |
pam-fi |
1.29 |
// cout << "void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1)"<<endl; |
752 |
|
|
Clear(); |
753 |
pam-fi |
1.9 |
// temporary objects: |
754 |
pam-fi |
1.29 |
TrkSinglet* t_singlet = new TrkSinglet(); |
755 |
|
|
TrkTrack* t_track = new TrkTrack(); |
756 |
|
|
|
757 |
|
|
// ----------------- |
758 |
|
|
// general variables |
759 |
|
|
// ----------------- |
760 |
|
|
for(Int_t i=0; i<12 ; i++){ |
761 |
|
|
good[i] = l2->good[i]; |
762 |
|
|
}; |
763 |
|
|
// -------------- |
764 |
|
|
// *** TRACKS *** |
765 |
|
|
// -------------- |
766 |
|
|
if(!Track) Track = new TClonesArray("TrkTrack"); |
767 |
|
|
TClonesArray &t = *Track; |
768 |
|
|
//----------------------------------------------------- |
769 |
|
|
if( l1 && !t_track->clx )t_track->clx = new TRefArray(6,0); |
770 |
|
|
if( l1 && !t_track->cly )t_track->cly = new TRefArray(6,0); |
771 |
|
|
//----------------------------------------------------- |
772 |
|
|
for(int i=0; i<l2->ntrk; i++){ |
773 |
|
|
// cout <<" TRACK "<<i<<" ------------------ "<<endl; |
774 |
|
|
t_track->seqno = i;// NBNBNBNB deve sempre essere = i |
775 |
|
|
t_track->image = l2->image[i]-1; |
776 |
|
|
t_track->chi2 = l2->chi2_nt[i]; |
777 |
|
|
t_track->nstep = l2->nstep_nt[i]; |
778 |
|
|
for(int it1=0;it1<5;it1++){ |
779 |
|
|
t_track->al[it1] = l2->al_nt[i][it1]; |
780 |
|
|
for(int it2=0;it2<5;it2++) |
781 |
|
|
t_track->coval[it1][it2] = l2->coval[i][it2][it1]; |
782 |
pam-fi |
1.9 |
}; |
783 |
pam-fi |
1.29 |
for(int ip=0;ip<6;ip++){ |
784 |
|
|
t_track->xgood[ip] = l2->cltrx[i][ip];//l2->xgood_nt[i][ip]; |
785 |
|
|
t_track->ygood[ip] = l2->cltry[i][ip];//l2->ygood_nt[i][ip]; |
786 |
|
|
t_track->xm[ip] = l2->xm_nt[i][ip]; |
787 |
|
|
t_track->ym[ip] = l2->ym_nt[i][ip]; |
788 |
|
|
t_track->zm[ip] = l2->zm_nt[i][ip]; |
789 |
|
|
t_track->resx[ip] = l2->resx_nt[i][ip]; |
790 |
|
|
t_track->resy[ip] = l2->resy_nt[i][ip]; |
791 |
|
|
t_track->xv[ip] = l2->xv_nt[i][ip]; |
792 |
|
|
t_track->yv[ip] = l2->yv_nt[i][ip]; |
793 |
|
|
t_track->zv[ip] = l2->zv_nt[i][ip]; |
794 |
|
|
t_track->axv[ip] = l2->axv_nt[i][ip]; |
795 |
|
|
t_track->ayv[ip] = l2->ayv_nt[i][ip]; |
796 |
|
|
t_track->dedx_x[ip] = l2->dedx_x[i][ip]; |
797 |
|
|
t_track->dedx_y[ip] = l2->dedx_y[i][ip]; |
798 |
|
|
// cout << " ip "<<ip<<" l2->cltrx[i][ip] "<< l2->cltrx[i][ip]<<endl; |
799 |
|
|
// cout << " ip "<<ip<<" l2->cltry[i][ip] "<< l2->cltry[i][ip]<<endl; |
800 |
|
|
//----------------------------------------------------- |
801 |
|
|
//----------------------------------------------------- |
802 |
|
|
// if(l1 && t_track->xgood[ip])t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip); |
803 |
|
|
// if(l1 && t_track->ygood[ip])t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip); |
804 |
|
|
if(l2->xgood_nt[i][ip]){ |
805 |
|
|
// cout << " ip "<<ip<<" l2->cltrx[i][ip] "<< l2->cltrx[i][ip]<<endl; |
806 |
|
|
// cout << " ip "<<ip<<" l2->cltrx[i][ip] "<< l2->cltrx[i][ip]<<" "; |
807 |
|
|
// if( l1->GetCluster(l2->cltrx[i][ip]-1)->TestBit(l1->GetCluster(l2->cltrx[i][ip]-1)->kIsReferenced) )cout << ">> is referenced "; |
808 |
|
|
if(l1)t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip); |
809 |
|
|
// cout << " --- "<<l1->GetCluster(l2->cltrx[i][ip]-1)->GetUniqueID()<<endl; |
810 |
|
|
// t_track->xgood[ip] = l2->cltrx[i][ip]; // WORK-AROUND ***** |
811 |
|
|
}else{ |
812 |
|
|
if(l1)t_track->clx->RemoveAt(ip); |
813 |
|
|
} |
814 |
|
|
if(l2->ygood_nt[i][ip]){ |
815 |
|
|
// cout << " ip "<<ip<<" l2->cltry[i][ip] "<< l2->cltry[i][ip]<<endl; |
816 |
|
|
// cout << " ip "<<ip<<" l2->cltry[i][ip] "<< l2->cltry[i][ip]<<" "; |
817 |
|
|
// if( l1->GetCluster(l2->cltry[i][ip]-1)->TestBit(l1->GetCluster(l2->cltry[i][ip]-1)->kIsReferenced) )cout << ">> is referenced "; |
818 |
|
|
if(l1)t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip); |
819 |
|
|
// cout << " --- "<<l1->GetCluster(l2->cltry[i][ip]-1)->GetUniqueID()<<endl; |
820 |
|
|
// t_track->ygood[ip] = l2->cltry[i][ip]; // WORK-AROUND ***** |
821 |
|
|
}else{ |
822 |
|
|
if(l1)t_track->cly->RemoveAt(ip); |
823 |
|
|
} |
824 |
|
|
//----------------------------------------------------- |
825 |
|
|
//----------------------------------------------------- |
826 |
|
|
}; |
827 |
|
|
new(t[i]) TrkTrack(*t_track); |
828 |
|
|
t_track->Clear(); |
829 |
|
|
}; |
830 |
|
|
// ---------------- |
831 |
|
|
// *** SINGLETS *** |
832 |
|
|
// ---------------- |
833 |
|
|
if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); |
834 |
|
|
TClonesArray &sx = *SingletX; |
835 |
|
|
for(int i=0; i<l2->nclsx; i++){ |
836 |
|
|
t_singlet->plane = l2->planex[i]; |
837 |
|
|
t_singlet->coord[0] = l2->xs[i][0]; |
838 |
|
|
t_singlet->coord[1] = l2->xs[i][1]; |
839 |
|
|
t_singlet->sgnl = l2->signlxs[i]; |
840 |
|
|
//----------------------------------------------------- |
841 |
|
|
if(l1) t_singlet->cls = l1->GetCluster(l2->clsx[i]-1); |
842 |
|
|
//----------------------------------------------------- |
843 |
|
|
new(sx[i]) TrkSinglet(*t_singlet); |
844 |
|
|
t_singlet->Clear(); |
845 |
|
|
} |
846 |
|
|
if(!SingletY)SingletY = new TClonesArray("TrkSinglet"); |
847 |
|
|
TClonesArray &sy = *SingletY; |
848 |
|
|
for(int i=0; i<l2->nclsy; i++){ |
849 |
|
|
t_singlet->plane = l2->planey[i]; |
850 |
|
|
t_singlet->coord[0] = l2->ys[i][0]; |
851 |
|
|
t_singlet->coord[1] = l2->ys[i][1]; |
852 |
|
|
t_singlet->sgnl = l2->signlys[i]; |
853 |
pam-fi |
1.26 |
//----------------------------------------------------- |
854 |
pam-fi |
1.29 |
if(l1) t_singlet->cls = l1->GetCluster(l2->clsy[i]-1); |
855 |
pam-fi |
1.26 |
//----------------------------------------------------- |
856 |
pam-fi |
1.29 |
new(sy[i]) TrkSinglet(*t_singlet); |
857 |
|
|
t_singlet->Clear(); |
858 |
|
|
}; |
859 |
pam-fi |
1.5 |
|
860 |
pam-fi |
1.29 |
delete t_track; |
861 |
|
|
delete t_singlet; |
862 |
mocchiut |
1.1 |
} |
863 |
pam-fi |
1.7 |
/** |
864 |
|
|
* Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common). |
865 |
|
|
*/ |
866 |
|
|
|
867 |
|
|
void TrkLevel2::GetLevel2Struct(cTrkLevel2 *l2) const { |
868 |
|
|
|
869 |
|
|
// general variables |
870 |
pam-fi |
1.10 |
// l2->good2 = good2 ; |
871 |
pam-fi |
1.7 |
for(Int_t i=0; i<12 ; i++){ |
872 |
pam-fi |
1.10 |
// l2->crc[i] = crc[i]; |
873 |
|
|
l2->good[i] = good[i]; |
874 |
pam-fi |
1.7 |
}; |
875 |
|
|
// *** TRACKS *** |
876 |
|
|
|
877 |
pam-fi |
1.21 |
if(Track){ |
878 |
|
|
l2->ntrk = Track->GetEntries(); |
879 |
|
|
for(Int_t i=0;i<l2->ntrk;i++){ |
880 |
|
|
l2->image[i] = 1 + ((TrkTrack *)Track->At(i))->image; |
881 |
|
|
l2->chi2_nt[i] = ((TrkTrack *)Track->At(i))->chi2; |
882 |
|
|
l2->nstep_nt[i] = ((TrkTrack *)Track->At(i))->nstep; |
883 |
|
|
for(int it1=0;it1<5;it1++){ |
884 |
|
|
l2->al_nt[i][it1] = ((TrkTrack *)Track->At(i))->al[it1]; |
885 |
|
|
for(int it2=0;it2<5;it2++) |
886 |
|
|
l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2]; |
887 |
|
|
}; |
888 |
|
|
for(int ip=0;ip<6;ip++){ |
889 |
pam-fi |
1.30 |
l2->xgood_nt[i][ip] = ((TrkTrack *)Track->At(i))->XGood(ip); |
890 |
|
|
l2->ygood_nt[i][ip] = ((TrkTrack *)Track->At(i))->YGood(ip); |
891 |
pam-fi |
1.21 |
l2->xm_nt[i][ip] = ((TrkTrack *)Track->At(i))->xm[ip]; |
892 |
|
|
l2->ym_nt[i][ip] = ((TrkTrack *)Track->At(i))->ym[ip]; |
893 |
|
|
l2->zm_nt[i][ip] = ((TrkTrack *)Track->At(i))->zm[ip]; |
894 |
|
|
l2->resx_nt[i][ip] = ((TrkTrack *)Track->At(i))->resx[ip]; |
895 |
|
|
l2->resy_nt[i][ip] = ((TrkTrack *)Track->At(i))->resy[ip]; |
896 |
|
|
l2->xv_nt[i][ip] = ((TrkTrack *)Track->At(i))->xv[ip]; |
897 |
|
|
l2->yv_nt[i][ip] = ((TrkTrack *)Track->At(i))->yv[ip]; |
898 |
|
|
l2->zv_nt[i][ip] = ((TrkTrack *)Track->At(i))->zv[ip]; |
899 |
|
|
l2->axv_nt[i][ip] = ((TrkTrack *)Track->At(i))->axv[ip]; |
900 |
|
|
l2->ayv_nt[i][ip] = ((TrkTrack *)Track->At(i))->ayv[ip]; |
901 |
|
|
l2->dedx_x[i][ip] = ((TrkTrack *)Track->At(i))->dedx_x[ip]; |
902 |
|
|
l2->dedx_y[i][ip] = ((TrkTrack *)Track->At(i))->dedx_y[ip]; |
903 |
|
|
}; |
904 |
|
|
} |
905 |
|
|
} |
906 |
|
|
// *** SINGLETS *** |
907 |
|
|
if(SingletX){ |
908 |
|
|
l2->nclsx = SingletX->GetEntries(); |
909 |
|
|
for(Int_t i=0;i<l2->nclsx;i++){ |
910 |
|
|
l2->planex[i] = ((TrkSinglet *)SingletX->At(i))->plane; |
911 |
|
|
l2->xs[i][0] = ((TrkSinglet *)SingletX->At(i))->coord[0]; |
912 |
|
|
l2->xs[i][1] = ((TrkSinglet *)SingletX->At(i))->coord[1]; |
913 |
|
|
l2->signlxs[i] = ((TrkSinglet *)SingletX->At(i))->sgnl; |
914 |
|
|
} |
915 |
pam-fi |
1.7 |
} |
916 |
|
|
|
917 |
pam-fi |
1.21 |
if(SingletY){ |
918 |
|
|
l2->nclsy = SingletY->GetEntries(); |
919 |
|
|
for(Int_t i=0;i<l2->nclsy;i++){ |
920 |
|
|
l2->planey[i] = ((TrkSinglet *)SingletY->At(i))->plane; |
921 |
|
|
l2->ys[i][0] = ((TrkSinglet *)SingletY->At(i))->coord[0]; |
922 |
|
|
l2->ys[i][1] = ((TrkSinglet *)SingletY->At(i))->coord[1]; |
923 |
|
|
l2->signlys[i] = ((TrkSinglet *)SingletY->At(i))->sgnl; |
924 |
|
|
} |
925 |
pam-fi |
1.7 |
} |
926 |
|
|
} |
927 |
mocchiut |
1.1 |
//-------------------------------------- |
928 |
|
|
// |
929 |
|
|
// |
930 |
|
|
//-------------------------------------- |
931 |
|
|
void TrkLevel2::Clear(){ |
932 |
|
|
for(Int_t i=0; i<12 ; i++){ |
933 |
pam-fi |
1.21 |
good[i] = -1; |
934 |
|
|
}; |
935 |
|
|
// if(Track)Track->Clear("C"); |
936 |
|
|
// if(SingletX)SingletX->Clear("C"); |
937 |
|
|
// if(SingletY)SingletY->Clear("C"); |
938 |
|
|
if(Track)Track->Delete(); |
939 |
|
|
if(SingletX)SingletX->Delete(); |
940 |
|
|
if(SingletY)SingletY->Delete(); |
941 |
|
|
} |
942 |
|
|
// //-------------------------------------- |
943 |
|
|
// // |
944 |
|
|
// // |
945 |
|
|
// //-------------------------------------- |
946 |
pam-fi |
1.11 |
void TrkLevel2::Delete(){ |
947 |
|
|
|
948 |
pam-fi |
1.21 |
// cout << "void TrkLevel2::Delete()"<<endl; |
949 |
|
|
Clear(); |
950 |
|
|
if(Track)delete Track; |
951 |
|
|
if(SingletX)delete SingletX; |
952 |
|
|
if(SingletY)delete SingletY; |
953 |
|
|
|
954 |
pam-fi |
1.11 |
} |
955 |
|
|
//-------------------------------------- |
956 |
|
|
// |
957 |
|
|
// |
958 |
|
|
//-------------------------------------- |
959 |
mocchiut |
1.1 |
/** |
960 |
|
|
* Sort physical tracks and stores them in a TObjectArray, ordering by increasing chi**2 value (in case of track image, it selects the one with lower chi**2). The total number of physical tracks is given by GetNTracks() and the it-th physical track can be retrieved by means of the method GetTrack(int it). |
961 |
|
|
* This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used. |
962 |
|
|
*/ |
963 |
pam-fi |
1.8 |
TRefArray *TrkLevel2::GetTracks_NFitSorted(){ |
964 |
pam-fi |
1.3 |
|
965 |
pam-fi |
1.21 |
if(!Track)return 0; |
966 |
|
|
|
967 |
|
|
TRefArray *sorted = new TRefArray(); |
968 |
pam-fi |
1.8 |
|
969 |
pam-fi |
1.21 |
TClonesArray &t = *Track; |
970 |
pam-fi |
1.8 |
// TClonesArray &ts = *PhysicalTrack; |
971 |
pam-fi |
1.21 |
int N = ntrk(); |
972 |
|
|
vector<int> m(N); for(int i=0; i<N; i++)m[i]=1; |
973 |
pam-fi |
1.8 |
// int m[50]; for(int i=0; i<N; i++)m[i]=1; |
974 |
|
|
|
975 |
pam-fi |
1.21 |
int indo=0; |
976 |
|
|
int indi=0; |
977 |
pam-fi |
1.26 |
while(N > 0){ |
978 |
|
|
// while(N != 0){ |
979 |
pam-fi |
1.21 |
int nfit =0; |
980 |
|
|
float chi2ref = numeric_limits<float>::max(); |
981 |
pam-fi |
1.9 |
|
982 |
pam-fi |
1.21 |
// first loop to search maximum num. of fit points |
983 |
|
|
for(int i=0; i < ntrk(); i++){ |
984 |
|
|
if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){ |
985 |
|
|
nfit = ((TrkTrack *)t[i])->GetNtot(); |
986 |
|
|
} |
987 |
|
|
} |
988 |
|
|
//second loop to search minimum chi2 among selected |
989 |
pam-fi |
1.26 |
for(int i=0; i<ntrk(); i++){ |
990 |
pam-fi |
1.21 |
Float_t chi2 = ((TrkTrack *)t[i])->chi2; |
991 |
pam-fi |
1.26 |
if(chi2 < 0) chi2 = -chi2*1000; |
992 |
pam-fi |
1.21 |
if( chi2 < chi2ref |
993 |
|
|
&& ((TrkTrack *)t[i])->GetNtot() == nfit |
994 |
|
|
&& m[i]==1){ |
995 |
|
|
chi2ref = ((TrkTrack *)t[i])->chi2; |
996 |
|
|
indi = i; |
997 |
|
|
}; |
998 |
|
|
}; |
999 |
|
|
if( ((TrkTrack *)t[indi])->HasImage() ){ |
1000 |
|
|
m[((TrkTrack *)t[indi])->image] = 0; |
1001 |
|
|
N--; |
1002 |
pam-fi |
1.8 |
|
1003 |
pam-fi |
1.26 |
// cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl; |
1004 |
pam-fi |
1.21 |
}; |
1005 |
|
|
sorted->Add( (TrkTrack*)t[indi] ); |
1006 |
pam-fi |
1.3 |
|
1007 |
pam-fi |
1.21 |
m[indi] = 0; |
1008 |
pam-fi |
1.26 |
// cout << "SORTED "<< indo << " "<< indi << " "<< N << " "<<((TrkTrack *)t[indi])->image<<" "<<chi2ref<<endl; |
1009 |
pam-fi |
1.21 |
N--; |
1010 |
|
|
indo++; |
1011 |
|
|
} |
1012 |
|
|
m.clear(); |
1013 |
pam-fi |
1.26 |
// cout << "GetTracks_NFitSorted(it): Done"<< endl; |
1014 |
pam-fi |
1.8 |
|
1015 |
pam-fi |
1.21 |
return sorted; |
1016 |
pam-fi |
1.6 |
// return PhysicalTrack; |
1017 |
pam-fi |
1.3 |
} |
1018 |
mocchiut |
1.1 |
//-------------------------------------- |
1019 |
|
|
// |
1020 |
|
|
// |
1021 |
|
|
//-------------------------------------- |
1022 |
|
|
/** |
1023 |
|
|
* Retrieves the is-th stored track. |
1024 |
|
|
* @param it Track number, ranging from 0 to ntrk(). |
1025 |
|
|
* Fitted tracks ( images included ) are stored in a TObjectArray ( TrkLevel2::Track ) in the same order they are returned by the F77 fitting routine. |
1026 |
|
|
*/ |
1027 |
|
|
TrkTrack *TrkLevel2::GetStoredTrack(int is){ |
1028 |
|
|
|
1029 |
|
|
if(is >= this->ntrk()){ |
1030 |
|
|
cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl; |
1031 |
|
|
cout << " Stored tracks ntrk() = "<< this->ntrk() << endl; |
1032 |
|
|
return 0; |
1033 |
|
|
} |
1034 |
pam-fi |
1.21 |
if(!Track){ |
1035 |
|
|
cout << "TrkTrack *TrkLevel2::GetStoredTrack(int is) >> (TClonesArray*) Track ==0 "<<endl; |
1036 |
|
|
}; |
1037 |
mocchiut |
1.1 |
TClonesArray &t = *(Track); |
1038 |
|
|
TrkTrack *track = (TrkTrack*)t[is]; |
1039 |
|
|
return track; |
1040 |
|
|
} |
1041 |
|
|
//-------------------------------------- |
1042 |
|
|
// |
1043 |
|
|
// |
1044 |
|
|
//-------------------------------------- |
1045 |
|
|
/** |
1046 |
pam-fi |
1.6 |
* Retrieves the is-th stored X singlet. |
1047 |
|
|
* @param it Singlet number, ranging from 0 to nclsx(). |
1048 |
|
|
*/ |
1049 |
|
|
TrkSinglet *TrkLevel2::GetSingletX(int is){ |
1050 |
|
|
|
1051 |
|
|
if(is >= this->nclsx()){ |
1052 |
|
|
cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl; |
1053 |
|
|
cout << " Stored x-singlets nclsx() = "<< this->nclsx() << endl; |
1054 |
|
|
return 0; |
1055 |
|
|
} |
1056 |
pam-fi |
1.21 |
if(!SingletX)return 0; |
1057 |
pam-fi |
1.6 |
TClonesArray &t = *(SingletX); |
1058 |
|
|
TrkSinglet *singlet = (TrkSinglet*)t[is]; |
1059 |
|
|
return singlet; |
1060 |
|
|
} |
1061 |
|
|
//-------------------------------------- |
1062 |
|
|
// |
1063 |
|
|
// |
1064 |
|
|
//-------------------------------------- |
1065 |
|
|
/** |
1066 |
|
|
* Retrieves the is-th stored Y singlet. |
1067 |
|
|
* @param it Singlet number, ranging from 0 to nclsx(). |
1068 |
|
|
*/ |
1069 |
|
|
TrkSinglet *TrkLevel2::GetSingletY(int is){ |
1070 |
|
|
|
1071 |
|
|
if(is >= this->nclsy()){ |
1072 |
|
|
cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl; |
1073 |
|
|
cout << " Stored y-singlets nclsy() = "<< this->nclsx() << endl; |
1074 |
|
|
return 0; |
1075 |
|
|
} |
1076 |
pam-fi |
1.21 |
if(!SingletY)return 0; |
1077 |
pam-fi |
1.6 |
TClonesArray &t = *(SingletY); |
1078 |
|
|
TrkSinglet *singlet = (TrkSinglet*)t[is]; |
1079 |
|
|
return singlet; |
1080 |
|
|
} |
1081 |
|
|
//-------------------------------------- |
1082 |
|
|
// |
1083 |
|
|
// |
1084 |
|
|
//-------------------------------------- |
1085 |
|
|
/** |
1086 |
mocchiut |
1.1 |
* Retrieves the it-th "physical" track, sorted by the method GetNTracks(). |
1087 |
|
|
* @param it Track number, ranging from 0 to GetNTracks(). |
1088 |
|
|
*/ |
1089 |
pam-fi |
1.10 |
|
1090 |
pam-fi |
1.8 |
TrkTrack *TrkLevel2::GetTrack(int it){ |
1091 |
|
|
|
1092 |
|
|
if(it >= this->GetNTracks()){ |
1093 |
|
|
cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; |
1094 |
|
|
cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl; |
1095 |
|
|
return 0; |
1096 |
|
|
} |
1097 |
|
|
|
1098 |
|
|
TRefArray *sorted = GetTracks(); //TEMPORANEO |
1099 |
pam-fi |
1.21 |
if(!sorted)return 0; |
1100 |
pam-fi |
1.8 |
TrkTrack *track = (TrkTrack*)sorted->At(it); |
1101 |
pam-fi |
1.21 |
sorted->Clear(); |
1102 |
|
|
delete sorted; |
1103 |
pam-fi |
1.8 |
return track; |
1104 |
mocchiut |
1.1 |
} |
1105 |
pam-fi |
1.6 |
/** |
1106 |
|
|
* Give the number of "physical" tracks, sorted by the method GetTracks(). |
1107 |
|
|
*/ |
1108 |
pam-fi |
1.5 |
Int_t TrkLevel2::GetNTracks(){ |
1109 |
pam-fi |
1.8 |
|
1110 |
|
|
Float_t ntot=0; |
1111 |
pam-fi |
1.21 |
if(!Track)return 0; |
1112 |
pam-fi |
1.8 |
TClonesArray &t = *Track; |
1113 |
mocchiut |
1.12 |
for(int i=0; i<ntrk(); i++) { |
1114 |
pam-fi |
1.8 |
if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.; |
1115 |
|
|
else ntot+=0.5; |
1116 |
|
|
} |
1117 |
|
|
return (Int_t)ntot; |
1118 |
|
|
|
1119 |
pam-fi |
1.5 |
}; |
1120 |
mocchiut |
1.1 |
//-------------------------------------- |
1121 |
|
|
// |
1122 |
|
|
// |
1123 |
|
|
//-------------------------------------- |
1124 |
|
|
/** |
1125 |
|
|
* Retrieves (if present) the image of the it-th "physical" track, sorted by the method GetNTracks(). |
1126 |
|
|
* @param it Track number, ranging from 0 to GetNTracks(). |
1127 |
|
|
*/ |
1128 |
pam-fi |
1.8 |
TrkTrack *TrkLevel2::GetTrackImage(int it){ |
1129 |
|
|
|
1130 |
pam-fi |
1.21 |
if(it >= this->GetNTracks()){ |
1131 |
|
|
cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; |
1132 |
|
|
cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl; |
1133 |
|
|
return 0; |
1134 |
|
|
} |
1135 |
pam-fi |
1.8 |
|
1136 |
pam-fi |
1.21 |
TRefArray* sorted = GetTracks(); //TEMPORANEO |
1137 |
|
|
if(!sorted)return 0; |
1138 |
|
|
TrkTrack *track = (TrkTrack*)sorted->At(it); |
1139 |
pam-fi |
1.8 |
|
1140 |
pam-fi |
1.21 |
if(!track->HasImage()){ |
1141 |
|
|
cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl; |
1142 |
|
|
return 0; |
1143 |
|
|
} |
1144 |
|
|
if(!Track)return 0; |
1145 |
|
|
TrkTrack *image = (TrkTrack*)(*Track)[track->image]; |
1146 |
|
|
|
1147 |
|
|
sorted->Delete(); |
1148 |
|
|
delete sorted; |
1149 |
pam-fi |
1.8 |
|
1150 |
pam-fi |
1.21 |
return image; |
1151 |
pam-fi |
1.8 |
|
1152 |
mocchiut |
1.1 |
} |
1153 |
|
|
//-------------------------------------- |
1154 |
|
|
// |
1155 |
|
|
// |
1156 |
|
|
//-------------------------------------- |
1157 |
|
|
/** |
1158 |
|
|
* Loads the magnetic field. |
1159 |
|
|
* @param s Path of the magnetic-field files. |
1160 |
|
|
*/ |
1161 |
pam-fi |
1.16 |
void TrkLevel2::LoadField(TString path){ |
1162 |
|
|
// |
1163 |
pam-fi |
1.26 |
// strcpy(path_.path,path.Data()); |
1164 |
|
|
// path_.pathlen = path.Length(); |
1165 |
|
|
// path_.error = 0; |
1166 |
|
|
// readb_(); |
1167 |
|
|
|
1168 |
|
|
TrkParams::Set(path,1); |
1169 |
pam-fi |
1.28 |
TrkParams::Load(1); |
1170 |
pam-fi |
1.26 |
|
1171 |
pam-fi |
1.16 |
// |
1172 |
mocchiut |
1.1 |
}; |
1173 |
pam-fi |
1.25 |
/** |
1174 |
|
|
* Get BY (kGauss) |
1175 |
|
|
* @param v (x,y,z) coordinates in cm |
1176 |
|
|
*/ |
1177 |
|
|
float TrkLevel2::GetBX(float* v){ |
1178 |
|
|
float b[3]; |
1179 |
|
|
gufld_(v,b); |
1180 |
|
|
return b[0]/10.; |
1181 |
|
|
} |
1182 |
|
|
/** |
1183 |
|
|
* Get BY (kGauss) |
1184 |
|
|
* @param v (x,y,z) coordinates in cm |
1185 |
|
|
*/ |
1186 |
|
|
float TrkLevel2::GetBY(float* v){ |
1187 |
|
|
float b[3]; |
1188 |
|
|
gufld_(v,b); |
1189 |
|
|
return b[1]/10.; |
1190 |
|
|
} |
1191 |
|
|
/** |
1192 |
|
|
* Get BY (kGauss) |
1193 |
|
|
* @param v (x,y,z) coordinates in cm |
1194 |
|
|
*/ |
1195 |
|
|
float TrkLevel2::GetBZ(float* v){ |
1196 |
|
|
float b[3]; |
1197 |
|
|
gufld_(v,b); |
1198 |
|
|
return b[2]/10.; |
1199 |
|
|
} |
1200 |
mocchiut |
1.1 |
//-------------------------------------- |
1201 |
|
|
// |
1202 |
|
|
// |
1203 |
|
|
//-------------------------------------- |
1204 |
|
|
/** |
1205 |
pam-fi |
1.6 |
* Get tracker-plane (mechanical) z-coordinate |
1206 |
|
|
* @param plane_id plane index (1=TOP,2,3,4,5,6=BOTTOM) |
1207 |
|
|
*/ |
1208 |
|
|
Float_t TrkLevel2::GetZTrk(Int_t plane_id){ |
1209 |
|
|
switch(plane_id){ |
1210 |
|
|
case 1: return ZTRK1; |
1211 |
|
|
case 2: return ZTRK2; |
1212 |
|
|
case 3: return ZTRK3; |
1213 |
|
|
case 4: return ZTRK4; |
1214 |
|
|
case 5: return ZTRK5; |
1215 |
|
|
case 6: return ZTRK6; |
1216 |
|
|
default: return 0.; |
1217 |
|
|
}; |
1218 |
|
|
}; |
1219 |
|
|
//-------------------------------------- |
1220 |
|
|
// |
1221 |
|
|
// |
1222 |
|
|
//-------------------------------------- |
1223 |
|
|
/** |
1224 |
pam-fi |
1.2 |
* Trajectory default constructor. |
1225 |
|
|
* (By default is created with z-coordinates inside the tracking volume) |
1226 |
|
|
*/ |
1227 |
|
|
Trajectory::Trajectory(){ |
1228 |
|
|
npoint = 10; |
1229 |
|
|
x = new float[npoint]; |
1230 |
|
|
y = new float[npoint]; |
1231 |
|
|
z = new float[npoint]; |
1232 |
|
|
thx = new float[npoint]; |
1233 |
|
|
thy = new float[npoint]; |
1234 |
|
|
tl = new float[npoint]; |
1235 |
pam-fi |
1.6 |
float dz = ((ZTRK1)-(ZTRK6))/(npoint-1); |
1236 |
pam-fi |
1.2 |
for(int i=0; i<npoint; i++){ |
1237 |
|
|
x[i] = 0; |
1238 |
|
|
y[i] = 0; |
1239 |
pam-fi |
1.6 |
z[i] = (ZTRK1) - i*dz; |
1240 |
pam-fi |
1.2 |
thx[i] = 0; |
1241 |
|
|
thy[i] = 0; |
1242 |
|
|
tl[i] = 0; |
1243 |
|
|
} |
1244 |
|
|
} |
1245 |
|
|
//-------------------------------------- |
1246 |
|
|
// |
1247 |
|
|
// |
1248 |
|
|
//-------------------------------------- |
1249 |
|
|
/** |
1250 |
mocchiut |
1.1 |
* Trajectory constructor. |
1251 |
pam-fi |
1.2 |
* (By default is created with z-coordinates inside the tracking volume) |
1252 |
mocchiut |
1.1 |
* \param n Number of points |
1253 |
|
|
*/ |
1254 |
|
|
Trajectory::Trajectory(int n){ |
1255 |
pam-fi |
1.2 |
if(n<=0){ |
1256 |
|
|
cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl; |
1257 |
|
|
n=10; |
1258 |
|
|
} |
1259 |
mocchiut |
1.1 |
npoint = n; |
1260 |
|
|
x = new float[npoint]; |
1261 |
|
|
y = new float[npoint]; |
1262 |
|
|
z = new float[npoint]; |
1263 |
pam-fi |
1.2 |
thx = new float[npoint]; |
1264 |
|
|
thy = new float[npoint]; |
1265 |
|
|
tl = new float[npoint]; |
1266 |
pam-fi |
1.6 |
float dz = ((ZTRK1)-(ZTRK6))/(npoint-1); |
1267 |
mocchiut |
1.1 |
for(int i=0; i<npoint; i++){ |
1268 |
pam-fi |
1.2 |
x[i] = 0; |
1269 |
mocchiut |
1.1 |
y[i] = 0; |
1270 |
pam-fi |
1.6 |
z[i] = (ZTRK1) - i*dz; |
1271 |
pam-fi |
1.2 |
thx[i] = 0; |
1272 |
|
|
thy[i] = 0; |
1273 |
|
|
tl[i] = 0; |
1274 |
mocchiut |
1.1 |
} |
1275 |
|
|
} |
1276 |
|
|
//-------------------------------------- |
1277 |
|
|
// |
1278 |
|
|
// |
1279 |
|
|
//-------------------------------------- |
1280 |
|
|
/** |
1281 |
|
|
* Trajectory constructor. |
1282 |
|
|
* \param n Number of points |
1283 |
|
|
* \param pz Pointer to float array, defining z coordinates |
1284 |
|
|
*/ |
1285 |
|
|
Trajectory::Trajectory(int n, float* zin){ |
1286 |
pam-fi |
1.2 |
npoint = 10; |
1287 |
|
|
if(n>0)npoint = n; |
1288 |
mocchiut |
1.1 |
x = new float[npoint]; |
1289 |
|
|
y = new float[npoint]; |
1290 |
|
|
z = new float[npoint]; |
1291 |
pam-fi |
1.2 |
thx = new float[npoint]; |
1292 |
|
|
thy = new float[npoint]; |
1293 |
|
|
tl = new float[npoint]; |
1294 |
|
|
int i=0; |
1295 |
|
|
do{ |
1296 |
pam-fi |
1.21 |
x[i] = 0; |
1297 |
|
|
y[i] = 0; |
1298 |
|
|
z[i] = zin[i]; |
1299 |
|
|
thx[i] = 0; |
1300 |
|
|
thy[i] = 0; |
1301 |
|
|
tl[i] = 0; |
1302 |
|
|
i++; |
1303 |
pam-fi |
1.2 |
}while(zin[i-1] > zin[i] && i < npoint); |
1304 |
|
|
npoint=i; |
1305 |
|
|
if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl; |
1306 |
mocchiut |
1.1 |
} |
1307 |
pam-fi |
1.21 |
void Trajectory::Delete(){ |
1308 |
|
|
|
1309 |
|
|
if(x) delete [] x; |
1310 |
|
|
if(y) delete [] y; |
1311 |
|
|
if(z) delete [] z; |
1312 |
|
|
if(thx) delete [] thx; |
1313 |
|
|
if(thy) delete [] thy; |
1314 |
|
|
if(tl) delete [] tl; |
1315 |
|
|
|
1316 |
|
|
} |
1317 |
mocchiut |
1.1 |
//-------------------------------------- |
1318 |
|
|
// |
1319 |
|
|
// |
1320 |
|
|
//-------------------------------------- |
1321 |
|
|
/** |
1322 |
|
|
* Dump the trajectory coordinates. |
1323 |
|
|
*/ |
1324 |
|
|
void Trajectory::Dump(){ |
1325 |
|
|
cout <<endl<< "Trajectory ========== "<<endl; |
1326 |
|
|
for (int i=0; i<npoint; i++){ |
1327 |
pam-fi |
1.2 |
cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ; |
1328 |
|
|
cout <<" -- " << thx[i] <<" "<< thy[i] ; |
1329 |
|
|
cout <<" -- " << tl[i] << endl; |
1330 |
mocchiut |
1.1 |
}; |
1331 |
|
|
} |
1332 |
pam-fi |
1.2 |
//-------------------------------------- |
1333 |
|
|
// |
1334 |
|
|
// |
1335 |
|
|
//-------------------------------------- |
1336 |
|
|
/** |
1337 |
|
|
* Get trajectory length between two points |
1338 |
|
|
* @param ifirst first point (default 0) |
1339 |
|
|
* @param ilast last point (default npoint) |
1340 |
|
|
*/ |
1341 |
|
|
float Trajectory::GetLength(int ifirst, int ilast){ |
1342 |
|
|
if( ifirst<0 ) ifirst = 0; |
1343 |
|
|
if( ilast>=npoint) ilast = npoint-1; |
1344 |
|
|
float l=0; |
1345 |
|
|
for(int i=ifirst;i<=ilast;i++){ |
1346 |
|
|
l=l+tl[i]; |
1347 |
|
|
}; |
1348 |
|
|
if(z[ilast] > ZINI)l=l-tl[ilast]; |
1349 |
|
|
if(z[ifirst] < ZINI) l=l-tl[ifirst]; |
1350 |
|
|
|
1351 |
|
|
return l; |
1352 |
mocchiut |
1.1 |
|
1353 |
pam-fi |
1.2 |
} |
1354 |
pam-fi |
1.6 |
|
1355 |
pam-fi |
1.19 |
/** |
1356 |
|
|
* Evaluates the trajectory in the apparatus associated to the track. |
1357 |
|
|
* It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling TrkLevel2::LoadField() ), otherwise an error message is returned. |
1358 |
|
|
* @param t pointer to an object of the class Trajectory, |
1359 |
|
|
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
1360 |
|
|
* @return error flag. |
1361 |
|
|
*/ |
1362 |
|
|
int Trajectory::DoTrack2(float* al){ |
1363 |
|
|
|
1364 |
|
|
double *dxout = new double[npoint]; |
1365 |
|
|
double *dyout = new double[npoint]; |
1366 |
|
|
double *dthxout = new double[npoint]; |
1367 |
|
|
double *dthyout = new double[npoint]; |
1368 |
|
|
double *dtlout = new double[npoint]; |
1369 |
|
|
double *dzin = new double[npoint]; |
1370 |
|
|
double dal[5]; |
1371 |
|
|
|
1372 |
|
|
int ifail = 0; |
1373 |
|
|
|
1374 |
|
|
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
1375 |
|
|
for (int i=0; i<npoint; i++) dzin[i] = (double)z[i]; |
1376 |
|
|
|
1377 |
pam-fi |
1.26 |
TrkParams::Load(1); |
1378 |
|
|
if( !TrkParams::IsLoaded(1) ){ |
1379 |
|
|
cout << "int Trajectory::DoTrack2(float* al) --- ERROR --- m.field not loaded"<<endl; |
1380 |
|
|
return 0; |
1381 |
|
|
} |
1382 |
pam-fi |
1.19 |
dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
1383 |
|
|
|
1384 |
|
|
for (int i=0; i<npoint; i++){ |
1385 |
|
|
x[i] = (float)*dxout++; |
1386 |
|
|
y[i] = (float)*dyout++; |
1387 |
|
|
thx[i] = (float)*dthxout++; |
1388 |
|
|
thy[i] = (float)*dthyout++; |
1389 |
|
|
tl[i] = (float)*dtlout++; |
1390 |
|
|
} |
1391 |
|
|
|
1392 |
|
|
return ifail; |
1393 |
|
|
}; |
1394 |
pam-fi |
1.10 |
|
1395 |
mocchiut |
1.1 |
ClassImp(TrkLevel2); |
1396 |
|
|
ClassImp(TrkSinglet); |
1397 |
|
|
ClassImp(TrkTrack); |
1398 |
|
|
ClassImp(Trajectory); |