4 |
*/ |
*/ |
5 |
#include <TrkLevel2.h> |
#include <TrkLevel2.h> |
6 |
#include <iostream> |
#include <iostream> |
7 |
|
#include <math.h> |
8 |
using namespace std; |
using namespace std; |
9 |
//...................................... |
//...................................... |
10 |
// F77 routines |
// F77 routines |
12 |
extern "C" { |
extern "C" { |
13 |
void dotrack_(int*, double*, double*, double*, double*, int*); |
void dotrack_(int*, double*, double*, double*, double*, int*); |
14 |
void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*); |
void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*); |
15 |
int readb_(const char*); |
void dotrack3_(int*, double*, double*, double*, double*,double*, double*, double*,double*,int*); |
16 |
|
void mini2_(int*,int*,int*); |
17 |
|
void guess_(); |
18 |
|
void gufld_(float*, float*); |
19 |
|
float risxeta2_(float *); |
20 |
|
float risxeta3_(float *); |
21 |
|
float risxeta4_(float *); |
22 |
|
float risyeta2_(float *); |
23 |
} |
} |
24 |
|
|
25 |
//-------------------------------------- |
//-------------------------------------- |
26 |
// |
// |
27 |
// |
// |
28 |
//-------------------------------------- |
//-------------------------------------- |
29 |
TrkTrack::TrkTrack(){ |
TrkTrack::TrkTrack(){ |
30 |
|
// cout << "TrkTrack::TrkTrack()" << endl; |
31 |
seqno = -1; |
seqno = -1; |
32 |
image = -1; |
image = -1; |
33 |
chi2 = 0; |
chi2 = 0; |
34 |
|
nstep = 0; |
35 |
for(int it1=0;it1<5;it1++){ |
for(int it1=0;it1<5;it1++){ |
36 |
al[it1] = 0; |
al[it1] = 0; |
37 |
for(int it2=0;it2<5;it2++) |
for(int it2=0;it2<5;it2++)coval[it1][it2] = 0; |
|
coval[it1][it2] = 0; |
|
38 |
}; |
}; |
39 |
for(int ip=0;ip<6;ip++){ |
for(int ip=0;ip<6;ip++){ |
40 |
xgood[ip] = 0; |
xgood[ip] = 0; |
44 |
zm[ip] = 0; |
zm[ip] = 0; |
45 |
resx[ip] = 0; |
resx[ip] = 0; |
46 |
resy[ip] = 0; |
resy[ip] = 0; |
47 |
|
tailx[ip] = 0; |
48 |
|
taily[ip] = 0; |
49 |
xv[ip] = 0; |
xv[ip] = 0; |
50 |
yv[ip] = 0; |
yv[ip] = 0; |
51 |
zv[ip] = 0; |
zv[ip] = 0; |
53 |
ayv[ip] = 0; |
ayv[ip] = 0; |
54 |
dedx_x[ip] = 0; |
dedx_x[ip] = 0; |
55 |
dedx_y[ip] = 0; |
dedx_y[ip] = 0; |
56 |
}; |
multmaxx[ip] = 0; |
57 |
|
multmaxy[ip] = 0; |
58 |
|
seedx[ip] = 0; |
59 |
|
seedy[ip] = 0; |
60 |
|
xpu[ip] = 0; |
61 |
|
ypu[ip] = 0; |
62 |
|
|
63 |
|
}; |
64 |
|
|
65 |
|
// TrkParams::SetTrackingMode(); |
66 |
|
// TrkParams::SetPrecisionFactor(); |
67 |
|
// TrkParams::SetStepMin(); |
68 |
|
TrkParams::SetMiniDefault(); |
69 |
|
TrkParams::SetPFA(); |
70 |
|
|
71 |
|
int ngf = TrkParams::nGF; |
72 |
|
for(int i=0; i<ngf; i++){ |
73 |
|
xGF[i] = 0.; |
74 |
|
yGF[i] = 0.; |
75 |
|
} |
76 |
|
|
77 |
|
|
78 |
}; |
}; |
79 |
//-------------------------------------- |
//-------------------------------------- |
80 |
// |
// |
84 |
seqno = t.seqno; |
seqno = t.seqno; |
85 |
image = t.image; |
image = t.image; |
86 |
chi2 = t.chi2; |
chi2 = t.chi2; |
87 |
|
nstep = t.nstep; |
88 |
for(int it1=0;it1<5;it1++){ |
for(int it1=0;it1<5;it1++){ |
89 |
al[it1] = t.al[it1]; |
al[it1] = t.al[it1]; |
90 |
for(int it2=0;it2<5;it2++) |
for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2]; |
|
coval[it1][it2] = t.coval[it1][it2]; |
|
91 |
}; |
}; |
92 |
for(int ip=0;ip<6;ip++){ |
for(int ip=0;ip<6;ip++){ |
93 |
xgood[ip] = t.xgood[ip]; |
xgood[ip] = t.xgood[ip]; |
97 |
zm[ip] = t.zm[ip]; |
zm[ip] = t.zm[ip]; |
98 |
resx[ip] = t.resx[ip]; |
resx[ip] = t.resx[ip]; |
99 |
resy[ip] = t.resy[ip]; |
resy[ip] = t.resy[ip]; |
100 |
|
tailx[ip] = t.tailx[ip]; |
101 |
|
taily[ip] = t.taily[ip]; |
102 |
xv[ip] = t.xv[ip]; |
xv[ip] = t.xv[ip]; |
103 |
yv[ip] = t.yv[ip]; |
yv[ip] = t.yv[ip]; |
104 |
zv[ip] = t.zv[ip]; |
zv[ip] = t.zv[ip]; |
106 |
ayv[ip] = t.ayv[ip]; |
ayv[ip] = t.ayv[ip]; |
107 |
dedx_x[ip] = t.dedx_x[ip]; |
dedx_x[ip] = t.dedx_x[ip]; |
108 |
dedx_y[ip] = t.dedx_y[ip]; |
dedx_y[ip] = t.dedx_y[ip]; |
109 |
}; |
multmaxx[ip] = t.multmaxx[ip]; |
110 |
|
multmaxy[ip] = t.multmaxy[ip]; |
111 |
|
seedx[ip] = t.seedx[ip]; |
112 |
|
seedy[ip] = t.seedy[ip]; |
113 |
|
xpu[ip] = t.xpu[ip]; |
114 |
|
ypu[ip] = t.ypu[ip]; |
115 |
|
}; |
116 |
|
|
117 |
|
// TrkParams::SetTrackingMode(); |
118 |
|
// TrkParams::SetPrecisionFactor(); |
119 |
|
// TrkParams::SetStepMin(); |
120 |
|
TrkParams::SetMiniDefault(); |
121 |
|
TrkParams::SetPFA(); |
122 |
|
|
123 |
|
int ngf = TrkParams::nGF; |
124 |
|
for(int i=0; i<ngf; i++){ |
125 |
|
xGF[i] = t.xGF[i]; |
126 |
|
yGF[i] = t.yGF[i]; |
127 |
|
} |
128 |
}; |
}; |
129 |
//-------------------------------------- |
//-------------------------------------- |
130 |
// |
// |
131 |
// |
// |
132 |
//-------------------------------------- |
//-------------------------------------- |
133 |
/** |
void TrkTrack::Copy(TrkTrack& t){ |
|
* Evaluates the trajectory in the apparatus associated to the track. |
|
|
* 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. |
|
|
* @param t pointer to an object of the class Trajectory, |
|
|
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
|
|
* @return error flag. |
|
|
*/ |
|
|
int TrkTrack::DoTrack(Trajectory* t){ |
|
|
|
|
|
double *dxout = new double[t->npoint]; |
|
|
double *dyout = new double[t->npoint]; |
|
|
double *dzin = new double[t->npoint]; |
|
|
double dal[5]; |
|
134 |
|
|
135 |
int ifail = 0; |
t.seqno = seqno; |
136 |
|
t.image = image; |
137 |
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
t.chi2 = chi2; |
138 |
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; |
t.nstep = nstep; |
139 |
|
for(int it1=0;it1<5;it1++){ |
140 |
|
t.al[it1] = al[it1]; |
141 |
|
for(int it2=0;it2<5;it2++)t.coval[it1][it2] = coval[it1][it2]; |
142 |
|
}; |
143 |
|
for(int ip=0;ip<6;ip++){ |
144 |
|
t.xgood[ip] = xgood[ip]; |
145 |
|
t.ygood[ip] = ygood[ip]; |
146 |
|
t.xm[ip] = xm[ip]; |
147 |
|
t.ym[ip] = ym[ip]; |
148 |
|
t.zm[ip] = zm[ip]; |
149 |
|
t.resx[ip] = resx[ip]; |
150 |
|
t.resy[ip] = resy[ip]; |
151 |
|
t.tailx[ip] = tailx[ip]; |
152 |
|
t.taily[ip] = taily[ip]; |
153 |
|
t.xv[ip] = xv[ip]; |
154 |
|
t.yv[ip] = yv[ip]; |
155 |
|
t.zv[ip] = zv[ip]; |
156 |
|
t.axv[ip] = axv[ip]; |
157 |
|
t.ayv[ip] = ayv[ip]; |
158 |
|
t.dedx_x[ip] = dedx_x[ip]; |
159 |
|
t.dedx_y[ip] = dedx_y[ip]; |
160 |
|
t.multmaxx[ip] = multmaxx[ip]; |
161 |
|
t.multmaxy[ip] = multmaxy[ip]; |
162 |
|
t.seedx[ip] = seedx[ip]; |
163 |
|
t.seedy[ip] = seedy[ip]; |
164 |
|
t.xpu[ip] = xpu[ip]; |
165 |
|
t.ypu[ip] = ypu[ip]; |
166 |
|
|
167 |
|
}; |
168 |
|
int ngf = TrkParams::nGF; |
169 |
|
for(int i=0; i<ngf; i++){ |
170 |
|
t.xGF[i] = xGF[i]; |
171 |
|
t.yGF[i] = yGF[i]; |
172 |
|
} |
173 |
|
|
|
dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail); |
|
174 |
|
|
175 |
for (int i=0; i<t->npoint; i++){ |
}; |
176 |
t->x[i] = (float)*dxout++; |
//-------------------------------------- |
177 |
t->y[i] = (float)*dyout++; |
// |
178 |
} |
// |
179 |
|
//-------------------------------------- |
180 |
|
/** |
181 |
|
* |
182 |
|
* >>> OBSOLETE !!! use TrkTrack::DoTrack(Trajectory* t) instead |
183 |
|
* |
184 |
|
*/ |
185 |
|
int TrkTrack::DoTrack2(Trajectory* t){ |
186 |
|
|
187 |
// delete [] dxout; |
cout << endl; |
188 |
// delete [] dyout; |
cout << " int TrkTrack::DoTrack2(Trajectory* t) --->> NB NB !! this method is going to be eliminated !!! "<<endl; |
189 |
// delete [] dzin; |
cout << " >>>> replace it with TrkTrack::DoTrack(Trajectory* t) <<<<"<<endl; |
190 |
|
cout << " (Sorry Wolfgang!! Don't be totally confused!! By Elena)"<<endl; |
191 |
|
cout << endl; |
192 |
|
|
193 |
|
return DoTrack(t); |
194 |
|
|
|
return ifail; |
|
195 |
}; |
}; |
196 |
//-------------------------------------- |
//-------------------------------------- |
197 |
// |
// |
198 |
// |
// |
199 |
//-------------------------------------- |
//-------------------------------------- |
200 |
/** |
/** |
201 |
* Evaluates the trajectory in the apparatus associated to the track. |
* Evaluates the trajectory in the apparatus associated to the track state-vector. |
202 |
* 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. |
* It integrates the equations of motion in the magnetic field. |
203 |
* @param t pointer to an object of the class Trajectory, |
* @param t pointer to an object of the class Trajectory, |
204 |
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
* which z coordinates should be previously assigned. |
205 |
* @return error flag. |
* @return error flag. |
206 |
*/ |
*/ |
207 |
int TrkTrack::DoTrack2(Trajectory* t){ |
int TrkTrack::DoTrack(Trajectory* t){ |
208 |
|
|
209 |
double *dxout = new double[t->npoint]; |
double *dxout = new double[t->npoint]; |
210 |
double *dyout = new double[t->npoint]; |
double *dyout = new double[t->npoint]; |
219 |
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
220 |
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; |
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; |
221 |
|
|
222 |
|
TrkParams::Load(1); |
223 |
|
if( !TrkParams::IsLoaded(1) ){ |
224 |
|
cout << "int TrkTrack::DoTrack(Trajectory* t) --- ERROR --- m.field not loaded"<<endl; |
225 |
|
return 0; |
226 |
|
} |
227 |
dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
228 |
|
|
229 |
for (int i=0; i<t->npoint; i++){ |
for (int i=0; i<t->npoint; i++){ |
230 |
t->x[i] = (float)*dxout++; |
t->x[i] = (float)*(dxout+i); |
231 |
t->y[i] = (float)*dyout++; |
t->y[i] = (float)*(dyout+i); |
232 |
t->thx[i] = (float)*dthxout++; |
t->thx[i] = (float)*(dthxout+i); |
233 |
t->thy[i] = (float)*dthyout++; |
t->thy[i] = (float)*(dthyout+i); |
234 |
t->tl[i] = (float)*dtlout++; |
t->tl[i] = (float)*(dtlout+i); |
235 |
|
} |
236 |
|
|
237 |
|
delete [] dxout; |
238 |
|
delete [] dyout; |
239 |
|
delete [] dzin; |
240 |
|
delete [] dthxout; |
241 |
|
delete [] dthyout; |
242 |
|
delete [] dtlout; |
243 |
|
|
244 |
|
return ifail; |
245 |
|
}; |
246 |
|
//-------------------------------------- |
247 |
|
// |
248 |
|
// |
249 |
|
//-------------------------------------- |
250 |
|
//float TrkTrack::BdL(){ |
251 |
|
//}; |
252 |
|
//-------------------------------------- |
253 |
|
// |
254 |
|
// |
255 |
|
//-------------------------------------- |
256 |
|
Float_t TrkTrack::GetRigidity(){ |
257 |
|
Float_t rig=0; |
258 |
|
if(chi2>0)rig=1./al[4]; |
259 |
|
if(rig<0) rig=-rig; |
260 |
|
return rig; |
261 |
|
}; |
262 |
|
// |
263 |
|
Float_t TrkTrack::GetDeflection(){ |
264 |
|
Float_t def=0; |
265 |
|
if(chi2>0)def=al[4]; |
266 |
|
return def; |
267 |
|
}; |
268 |
|
// |
269 |
|
/** |
270 |
|
* Method to retrieve the dE/dx measured on a tracker view. |
271 |
|
* @param ip plane (0-5) |
272 |
|
* @param iv view (0=x 1=y) |
273 |
|
*/ |
274 |
|
Float_t TrkTrack::GetDEDX(int ip, int iv){ |
275 |
|
if(iv==0 && ip>=0 && ip<6)return fabs(dedx_x[ip]); |
276 |
|
else if(iv==1 && ip>=0 && ip<6)return fabs(dedx_y[ip]); |
277 |
|
else { |
278 |
|
cout << "TrkTrack::GetDEDX(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl; |
279 |
|
return 0.; |
280 |
|
} |
281 |
|
} |
282 |
|
/** |
283 |
|
* Method to evaluate the dE/dx measured on a tracker plane. |
284 |
|
* The two measurements on x- and y-view are averaged. |
285 |
|
* @param ip plane (0-5) |
286 |
|
*/ |
287 |
|
Float_t TrkTrack::GetDEDX(int ip){ |
288 |
|
if( (Int_t)XGood(ip)+(Int_t)YGood(ip) == 0 ) return 0; |
289 |
|
return (GetDEDX(ip,0)+GetDEDX(ip,1))/((Int_t)XGood(ip)+(Int_t)YGood(ip)); |
290 |
|
}; |
291 |
|
|
292 |
|
/** |
293 |
|
* Method to evaluate the dE/dx averaged over all planes. |
294 |
|
*/ |
295 |
|
Float_t TrkTrack::GetDEDX(){ |
296 |
|
Float_t dedx=0; |
297 |
|
for(Int_t ip=0; ip<6; ip++)dedx+=GetDEDX(ip,0)*XGood(ip)+GetDEDX(ip,1)*YGood(ip); |
298 |
|
dedx = dedx/(GetNX()+GetNY()); |
299 |
|
return dedx; |
300 |
|
}; |
301 |
|
/** |
302 |
|
* Returns 1 if the cluster on a tracker view includes bad strips |
303 |
|
* (at least one bad strip among the four strip used by p.f.a.) |
304 |
|
* @param ip plane (0-5) |
305 |
|
* @param iv view (0=x 1=y) |
306 |
|
*/ |
307 |
|
Bool_t TrkTrack::IsBad(int ip,int iv){ |
308 |
|
if(iv==0 && ip>=0 && ip<6)return (xgood[ip]<0) ; |
309 |
|
else if(iv==1 && ip>=0 && ip<6)return (ygood[ip]<0) ; |
310 |
|
else { |
311 |
|
cout << "TrkTrack::IsBad(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl; |
312 |
|
return 0.; |
313 |
|
} |
314 |
|
}; |
315 |
|
/** |
316 |
|
* Returns 1 if the signal on a tracker view is saturated. |
317 |
|
* @param ip plane (0-5) |
318 |
|
* @param iv view (0=x 1=y) |
319 |
|
*/ |
320 |
|
Bool_t TrkTrack::IsSaturated(int ip,int iv){ |
321 |
|
if(iv==0 && ip>=0 && ip<6)return (dedx_x[ip]<0) ; |
322 |
|
else if(iv==1 && ip>=0 && ip<6)return (dedx_y[ip]<0) ; |
323 |
|
else { |
324 |
|
cout << "TrkTrack::IsSaturated(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl; |
325 |
|
return 0.; |
326 |
|
} |
327 |
|
}; |
328 |
|
/** |
329 |
|
* Returns 1 if either the x or the y signal on a tracker plane is saturated. |
330 |
|
* @param ip plane (0-5) |
331 |
|
*/ |
332 |
|
Bool_t TrkTrack::IsSaturated(int ip){ |
333 |
|
return (IsSaturated(ip,0)||IsSaturated(ip,1)); |
334 |
|
}; |
335 |
|
/** |
336 |
|
* Returns 1 if there is at least a saturated signal along the track. |
337 |
|
*/ |
338 |
|
Bool_t TrkTrack::IsSaturated(){ |
339 |
|
for(int ip=0; ip<6; ip++)for(int iv=0; iv<2; iv++)if(IsSaturated(ip,iv))return true; |
340 |
|
return false; |
341 |
|
} |
342 |
|
/** |
343 |
|
* Returns the track "lever-arm" on the x view, defined as the distance (in planes) between |
344 |
|
* the upper and lower x measurements (the maximum value of lever-arm is 6). |
345 |
|
*/ |
346 |
|
Int_t TrkTrack::GetLeverArmX(){ |
347 |
|
int first_plane = -1; |
348 |
|
int last_plane = -1; |
349 |
|
for(Int_t ip=0; ip<6; ip++){ |
350 |
|
if( XGood(ip) && first_plane == -1 )first_plane = ip; |
351 |
|
if( XGood(ip) && first_plane != -1 )last_plane = ip; |
352 |
|
} |
353 |
|
if( first_plane == -1 || last_plane == -1){ |
354 |
|
cout<< "Int_t TrkTrack::GetLeverArmX() -- XGood(ip) always false ??? "<<endl; |
355 |
|
return 0; |
356 |
|
} |
357 |
|
return (last_plane-first_plane+1); |
358 |
|
} |
359 |
|
/** |
360 |
|
* Returns the track "lever-arm" on the y view, defined as the distance (in planes) between |
361 |
|
* the upper and lower y measurements (the maximum value of lever-arm is 6). |
362 |
|
*/ |
363 |
|
Int_t TrkTrack::GetLeverArmY(){ |
364 |
|
int first_plane = -1; |
365 |
|
int last_plane = -1; |
366 |
|
for(Int_t ip=0; ip<6; ip++){ |
367 |
|
if( YGood(ip) && first_plane == -1 )first_plane = ip; |
368 |
|
if( YGood(ip) && first_plane != -1 )last_plane = ip; |
369 |
|
} |
370 |
|
if( first_plane == -1 || last_plane == -1){ |
371 |
|
cout<< "Int_t TrkTrack::GetLeverArmY() -- YGood(ip) always false ??? "<<endl; |
372 |
|
return 0; |
373 |
|
} |
374 |
|
return (last_plane-first_plane+1); |
375 |
|
} |
376 |
|
/** |
377 |
|
* Returns the track "lever-arm" on the x+y view, defined as the distance (in planes) between |
378 |
|
* the upper and lower x,y (couple) measurements (the maximum value of lever-arm is 6). |
379 |
|
*/ |
380 |
|
Int_t TrkTrack::GetLeverArmXY(){ |
381 |
|
int first_plane = -1; |
382 |
|
int last_plane = -1; |
383 |
|
for(Int_t ip=0; ip<6; ip++){ |
384 |
|
if( XGood(ip) && YGood(ip) && first_plane == -1 )first_plane = ip; |
385 |
|
if( XGood(ip) && YGood(ip) && first_plane != -1 )last_plane = ip; |
386 |
|
} |
387 |
|
if( first_plane == -1 || last_plane == -1){ |
388 |
|
cout<< "Int_t TrkTrack::GetLeverArmXY() -- XGood(ip)*YGood(ip) always false ??? "<<endl; |
389 |
|
return 0; |
390 |
|
} |
391 |
|
return (last_plane-first_plane+1); |
392 |
|
} |
393 |
|
/** |
394 |
|
* Returns the number of hit planes |
395 |
|
*/ |
396 |
|
Int_t TrkTrack::GetNhit() { |
397 |
|
int np=0; |
398 |
|
for(Int_t ip=0; ip<6; ip++) np += (XGood(ip)||YGood(ip)) ; |
399 |
|
return np; |
400 |
|
}; |
401 |
|
/** |
402 |
|
* Returns the reduced chi-square of track x-projection |
403 |
|
*/ |
404 |
|
Float_t TrkTrack::GetChi2X(){ |
405 |
|
float chiq=0; |
406 |
|
for(int ip=0; ip<6; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.); |
407 |
|
if(GetNX()>3)chiq=chiq/(GetNX()-3); |
408 |
|
else chiq=0; |
409 |
|
// if(chiq==0)cout << " Float_t TrkTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl; |
410 |
|
return chiq; |
411 |
|
} |
412 |
|
/** |
413 |
|
* Returns the reduced chi-square of track y-projection |
414 |
|
*/ |
415 |
|
Float_t TrkTrack::GetChi2Y(){ |
416 |
|
float chiq=0; |
417 |
|
for(int ip=0; ip<6; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.); |
418 |
|
if(GetNY()>2)chiq=chiq/(GetNY()-2); |
419 |
|
else chiq=0; |
420 |
|
// if(chiq==0)cout << " Float_t TrkTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl; |
421 |
|
return chiq; |
422 |
|
} |
423 |
|
/** |
424 |
|
* Returns the logarythm of the likeliwood-function of track x-projection |
425 |
|
*/ |
426 |
|
Float_t TrkTrack::GetLnLX(){ |
427 |
|
float lnl=0; |
428 |
|
for(int ip=0; ip<6; ip++) |
429 |
|
if( XGood(ip) && tailx[ip]!=0 ) |
430 |
|
lnl += (tailx[ip]+1.) * log( (tailx[ip]*pow(resx[ip],2.) + pow(xv[ip]-xm[ip],2.)) / (tailx[ip]*pow(resx[ip],2)) ); |
431 |
|
if(GetNX()>3)lnl=lnl/(GetNX()-3); |
432 |
|
else lnl=0; |
433 |
|
if(lnl==0){ |
434 |
|
cout << " Float_t TrkTrack::GetLnLX() -- WARNING -- value not defined "<<lnl<<endl; |
435 |
|
Dump(); |
436 |
|
} |
437 |
|
return lnl; |
438 |
|
|
439 |
|
} |
440 |
|
/** |
441 |
|
* Returns the logarythm of the likeliwood-function of track y-projection |
442 |
|
*/ |
443 |
|
Float_t TrkTrack::GetLnLY(){ |
444 |
|
float lnl=0; |
445 |
|
for(int ip=0; ip<6; ip++) |
446 |
|
if( YGood(ip) && taily[ip]!=0 ) |
447 |
|
lnl += (taily[ip]+1.) * log( (taily[ip]*pow(resy[ip],2.) + pow(yv[ip]-ym[ip],2.)) / (taily[ip]*pow(resy[ip],2)) ); |
448 |
|
if(GetNY()>2)lnl=lnl/(GetNY()-2); |
449 |
|
else lnl=0; |
450 |
|
if(lnl==0){ |
451 |
|
cout << " Float_t TrkTrack::GetLnLY() -- WARNING -- value not defined "<<lnl<<endl; |
452 |
|
Dump(); |
453 |
|
} |
454 |
|
return lnl; |
455 |
|
|
456 |
|
} |
457 |
|
/** |
458 |
|
* Returns the effective angle, relative to the sensor, on each plane. |
459 |
|
* @param ip plane (0-5) |
460 |
|
* @param iv view (0=x 1=y) |
461 |
|
*/ |
462 |
|
Float_t TrkTrack::GetEffectiveAngle(int ip, int iv){ |
463 |
|
|
464 |
|
if(ip<0 || ip>5){ |
465 |
|
cout << "Float_t TrkTrack::GetEffectiveAngle(int "<<ip<<", int "<<iv<<") ==> wrong input"<<endl; |
466 |
|
return 0.; |
467 |
|
} |
468 |
|
|
469 |
|
float v[3]={xv[ip],yv[ip],zv[ip]}; |
470 |
|
//----------------------------------------- |
471 |
|
// effective angle (relative to the sensor) |
472 |
|
//----------------------------------------- |
473 |
|
float axv_geo = axv[ip]; |
474 |
|
float muhall_h = 297.61; //cm**2/Vs |
475 |
|
float BY = TrkParams::GetBY(v); |
476 |
|
float axv_eff = 0; |
477 |
|
if(ip==5) axv_geo = -1*axv_geo; |
478 |
|
if(ip==5) BY = -1*BY; |
479 |
|
axv_eff = 180.*atan( tan(axv_geo*acos(-1.)/180.) + muhall_h * BY * 0.0001)/acos(-1.); |
480 |
|
//----------------------------------------- |
481 |
|
// effective angle (relative to the sensor) |
482 |
|
//----------------------------------------- |
483 |
|
float ayv_geo = ayv[ip]; |
484 |
|
float muhall_e = 1258.18; //cm**2/Vs |
485 |
|
float BX = TrkParams::GetBX(v); |
486 |
|
float ayv_eff = 0; |
487 |
|
ayv_eff = 180.*atan( tan(ayv_geo*acos(-1.)/180.) + muhall_e * BX * 0.0001)/acos(-1.); |
488 |
|
|
489 |
|
if (iv==0)return axv_eff; |
490 |
|
else if(iv==1)return ayv_eff; |
491 |
|
else{ |
492 |
|
cout << "Float_t TrkTrack::GetEffectiveAngle(int "<<ip<<", int "<<iv<<") ==> wrong input"<<endl; |
493 |
|
return 0.; |
494 |
|
} |
495 |
|
|
496 |
|
}; |
497 |
|
|
498 |
|
//-------------------------------------- |
499 |
|
// |
500 |
|
// |
501 |
|
//-------------------------------------- |
502 |
|
void TrkTrack::Dump(){ |
503 |
|
cout << endl << "========== Track " ; |
504 |
|
cout << endl << "seq. n. : "<< seqno; |
505 |
|
cout << endl << "image n. : "<< image; |
506 |
|
cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " "; |
507 |
|
cout << endl << "chi^2 : "<< chi2; |
508 |
|
cout << endl << "n.step : "<< nstep; |
509 |
|
cout << endl << "xgood : "; for(int i=0; i<6; i++)cout << XGood(i) ; |
510 |
|
cout << endl << "ygood : "; for(int i=0; i<6; i++)cout << YGood(i) ; |
511 |
|
cout << endl << "xm : "; for(int i=0; i<6; i++)cout << xm[i] << " "; |
512 |
|
cout << endl << "ym : "; for(int i=0; i<6; i++)cout << ym[i] << " "; |
513 |
|
cout << endl << "zm : "; for(int i=0; i<6; i++)cout << zm[i] << " "; |
514 |
|
cout << endl << "xv : "; for(int i=0; i<6; i++)cout << xv[i] << " "; |
515 |
|
cout << endl << "yv : "; for(int i=0; i<6; i++)cout << yv[i] << " "; |
516 |
|
cout << endl << "zv : "; for(int i=0; i<6; i++)cout << zv[i] << " "; |
517 |
|
cout << endl << "resx : "; for(int i=0; i<6; i++)cout << resx[i] << " "; |
518 |
|
cout << endl << "resy : "; for(int i=0; i<6; i++)cout << resy[i] << " "; |
519 |
|
cout << endl << "tailx : "; for(int i=0; i<6; i++)cout << tailx[i] << " "; |
520 |
|
cout << endl << "taily : "; for(int i=0; i<6; i++)cout << taily[i] << " "; |
521 |
|
cout << endl << "coval : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" "; |
522 |
|
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[1][i]<<" "; |
523 |
|
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[2][i]<<" "; |
524 |
|
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[3][i]<<" "; |
525 |
|
cout << endl << " "; for(int i=0; i<5; i++)cout << coval[4][i]<<" "; |
526 |
|
cout << endl << "dedx_x : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " "; |
527 |
|
cout << endl << "dedx_y : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " "; |
528 |
|
cout << endl << "maxs x : "; for(int i=0; i<6; i++)cout << GetClusterX_MaxStrip(i) << " "; |
529 |
|
cout << endl << "maxs y : "; for(int i=0; i<6; i++)cout << GetClusterY_MaxStrip(i) << " "; |
530 |
|
cout << endl << "mult x : "; for(int i=0; i<6; i++)cout << GetClusterX_Multiplicity(i) << " "; |
531 |
|
cout << endl << "mult y : "; for(int i=0; i<6; i++)cout << GetClusterY_Multiplicity(i) << " "; |
532 |
|
cout << endl << "seed x : "; for(int i=0; i<6; i++)cout << GetClusterX_Seed(i) << " "; |
533 |
|
cout << endl << "seed y : "; for(int i=0; i<6; i++)cout << GetClusterY_Seed(i) << " "; |
534 |
|
cout << endl << "xpu : "; for(int i=0; i<6; i++)cout << xpu[i] << " "; |
535 |
|
cout << endl << "ypu : "; for(int i=0; i<6; i++)cout << ypu[i] << " "; |
536 |
|
|
537 |
|
cout << endl; |
538 |
|
} |
539 |
|
/** |
540 |
|
* Set the TrkTrack position measurements |
541 |
|
*/ |
542 |
|
void TrkTrack::SetMeasure(double *xmeas, double *ymeas, double *zmeas){ |
543 |
|
for(int i=0; i<6; i++) xm[i]=*xmeas++; |
544 |
|
for(int i=0; i<6; i++) ym[i]=*ymeas++; |
545 |
|
for(int i=0; i<6; i++) zm[i]=*zmeas++; |
546 |
|
} |
547 |
|
/** |
548 |
|
* Set the TrkTrack position resolution |
549 |
|
*/ |
550 |
|
void TrkTrack::SetResolution(double *rx, double *ry){ |
551 |
|
for(int i=0; i<6; i++) resx[i]=*rx++; |
552 |
|
for(int i=0; i<6; i++) resy[i]=*ry++; |
553 |
|
} |
554 |
|
/** |
555 |
|
* Set the TrkTrack tails position resolution |
556 |
|
*/ |
557 |
|
void TrkTrack::SetTail(double *tx, double *ty, double factor){ |
558 |
|
for(int i=0; i<6; i++) tailx[i]=factor*(*tx++); |
559 |
|
for(int i=0; i<6; i++) taily[i]=factor*(*ty++); |
560 |
|
} |
561 |
|
/** |
562 |
|
* Set the TrkTrack Student parameter (resx,resy,tailx,taily) |
563 |
|
* from previous gausian fit |
564 |
|
*@param flag =0 standard, =1 with noise correction |
565 |
|
*/ |
566 |
|
void TrkTrack::SetStudentParam(int flag){ |
567 |
|
float sx[11]={0.000128242, |
568 |
|
0.000136942, |
569 |
|
0.000162718, |
570 |
|
0.000202644, |
571 |
|
0.00025597, |
572 |
|
0.000317456, |
573 |
|
0.000349048, |
574 |
|
0.000384638, |
575 |
|
0.000457295, |
576 |
|
0.000512319, |
577 |
|
0.000538573}; |
578 |
|
float tx[11]={1.79402, |
579 |
|
2.04876, |
580 |
|
2.88376, |
581 |
|
3.3, |
582 |
|
3.14084, |
583 |
|
4.07686, |
584 |
|
4.44736, |
585 |
|
3.5179, |
586 |
|
3.38697, |
587 |
|
3.45739, |
588 |
|
3.18627}; |
589 |
|
float sy[11]={0.000483075, |
590 |
|
0.000466925, |
591 |
|
0.000431658, |
592 |
|
0.000428317, |
593 |
|
0.000433854, |
594 |
|
0.000444044, |
595 |
|
0.000482098, |
596 |
|
0.000537579, |
597 |
|
0.000636279, |
598 |
|
0.000741998, |
599 |
|
0.000864261}; |
600 |
|
float ty[11]={0.997032, |
601 |
|
1.11147, |
602 |
|
1.18526, |
603 |
|
1.61404, |
604 |
|
2.21908, |
605 |
|
3.08959, |
606 |
|
4.48833, |
607 |
|
4.42687, |
608 |
|
4.65253, |
609 |
|
4.52043, |
610 |
|
4.29926}; |
611 |
|
int index; |
612 |
|
float fact=0.; |
613 |
|
for(int i=0; i<6; i++) { |
614 |
|
index = int((fabs(axv[i])+1.)/2.); |
615 |
|
if(index>10) index=10; |
616 |
|
tailx[i]=tx[index]; |
617 |
|
if(flag==1) { |
618 |
|
if(fabs(axv[i])<=10.) fact = resx[i]/risxeta2_(&(axv[i])); |
619 |
|
if(fabs(axv[i])>10.&&fabs(axv[i])<=15.) fact = resx[i]/risxeta3_(&(axv[i])); |
620 |
|
if(fabs(axv[i])>15.) fact = resx[i]/risxeta4_(&(axv[i])); |
621 |
|
} else fact = 1.; |
622 |
|
resx[i] = sx[index]*fact; |
623 |
|
} |
624 |
|
for(int i=0; i<6; i++) { |
625 |
|
index = int((fabs(ayv[i])+1.)/2.); |
626 |
|
if(index>10) index=10; |
627 |
|
taily[i]=ty[index]; |
628 |
|
if(flag==1) fact = resy[i]/risyeta2_(&(ayv[i])); |
629 |
|
else fact = 1.; |
630 |
|
resy[i] = sy[index]*fact; |
631 |
|
} |
632 |
|
} |
633 |
|
/** |
634 |
|
* Set the TrkTrack good measurement |
635 |
|
*/ |
636 |
|
void TrkTrack::SetGood(int *xg, int *yg){ |
637 |
|
|
638 |
|
for(int i=0; i<6; i++) xgood[i]=*xg++; |
639 |
|
for(int i=0; i<6; i++) ygood[i]=*yg++; |
640 |
|
} |
641 |
|
|
642 |
|
/** |
643 |
|
* Load the magnetic field |
644 |
|
*/ |
645 |
|
void TrkTrack::LoadField(TString path){ |
646 |
|
|
647 |
|
// strcpy(path_.path,path.Data()); |
648 |
|
// path_.pathlen = path.Length(); |
649 |
|
// path_.error = 0; |
650 |
|
// readb_(); |
651 |
|
|
652 |
|
// TrkParams::SetTrackingMode(); |
653 |
|
// TrkParams::SetPrecisionFactor(); |
654 |
|
// TrkParams::SetStepMin(); |
655 |
|
TrkParams::SetMiniDefault(); |
656 |
|
|
657 |
|
TrkParams::Set(path,1); |
658 |
|
TrkParams::Load(1); |
659 |
|
if( !TrkParams::IsLoaded(1) ){ |
660 |
|
cout << "void TrkTrack::LoadField(TString path) --- ERROR --- m.field not loaded"<<endl; |
661 |
|
} |
662 |
|
|
663 |
|
}; |
664 |
|
|
665 |
|
|
666 |
|
/** |
667 |
|
* Method to fill minimization-routine common |
668 |
|
*/ |
669 |
|
void TrkTrack::FillMiniStruct(cMini2track& track){ |
670 |
|
|
671 |
|
for(int i=0; i<6; i++){ |
672 |
|
|
673 |
|
// cout << i<<" - "<<xgood[i]<<" "<<XGood(i)<<endl; |
674 |
|
// cout << i<<" - "<<ygood[i]<<" "<<YGood(i)<<endl; |
675 |
|
track.xgood[i] = (double) XGood(i); |
676 |
|
track.ygood[i] = (double) YGood(i); |
677 |
|
|
678 |
|
track.xm[i]= (double) xm[i]; |
679 |
|
track.ym[i]= (double) ym[i]; |
680 |
|
track.zm[i]= (double) zm[i]; |
681 |
|
|
682 |
|
// --- temporaneo ---------------------------- |
683 |
|
// float segment = 100.; |
684 |
|
// track.xm_a[i]=xm[i]; |
685 |
|
// track.xm_b[i]=xm[i]; |
686 |
|
// track.ym_a[i]=ym[i]; |
687 |
|
// track.ym_b[i]=ym[i]; |
688 |
|
// if( XGood(i) && !YGood(i) ){ |
689 |
|
// track.ym_a[i] = track.ym_a[i]+segment; |
690 |
|
// track.ym_b[i] = track.ym_b[i]-segment; |
691 |
|
// }else if( !XGood(i) && YGood(i)){ |
692 |
|
// track.xm_a[i] = track.xm_a[i]+segment; |
693 |
|
// track.xm_b[i] = track.xm_b[i]-segment; |
694 |
|
// } |
695 |
|
// --- temporaneo ---------------------------- |
696 |
|
|
697 |
|
if( XGood(i) || YGood(i) ){ |
698 |
|
//NB!! the length of the sensor is not exactely taken into account |
699 |
|
double segment = 7.;// 2.;//cm //Elena 10th |
700 |
|
// NB: i parametri di allineamento hanno una notazione particolare!!! |
701 |
|
// sensor = 0 (hybrid side), 1 |
702 |
|
// ladder = 0-2 (increasing x) |
703 |
|
// plane = 0-5 (from bottom to top!!!) |
704 |
|
int is = (int)GetSensor(i); if(i==5)is=1-is; |
705 |
|
int ip = 5-i; |
706 |
|
int il = (int)GetLadder(i); |
707 |
|
|
708 |
|
double omega = 0.; |
709 |
|
// double beta = 0.;// EM GCC 4.7 |
710 |
|
// double gamma = 0.; |
711 |
|
if( |
712 |
|
(is < 0 || is > 1 || ip < 0 || ip > 5 || il < 0 || il > 2) && |
713 |
|
true){ |
714 |
|
// se il piano risulta colpito, ladder e sensore devono essere |
715 |
|
// assegnati correttamente |
716 |
|
cout << " void TrkTrack::FillMiniStruct(cMini2track&) --- WARNING --- sensor not defined, cannot read alignment parameters "<<endl; |
717 |
|
cout << " is ip il = "<<is<<" "<<ip<<" "<<il<<endl; |
718 |
|
}else{ |
719 |
|
omega = alignparameters_.omega[is][il][ip]; |
720 |
|
// beta = alignparameters_.beta[is][il][ip];// EM GCC 4.7 unused |
721 |
|
// gamma = alignparameters_.gamma[is][il][ip];// EM GCC 4.7 unused |
722 |
|
} |
723 |
|
|
724 |
|
if( XGood(i) && !YGood(i) ){ |
725 |
|
track.xm_a[i] = (double) xm[i] - omega * segment; |
726 |
|
track.ym_a[i] = (double) ym[i] + segment; |
727 |
|
// track.zm_a[i] = zm[i] + beta * segment;//not used yet |
728 |
|
track.xm_b[i] = (double) xm[i] + omega * segment; |
729 |
|
track.ym_b[i] = (double) ym[i] - segment; |
730 |
|
// track.zm_b[i] = zm[i] - beta * segment;//not used yet |
731 |
|
}else if( !XGood(i) && YGood(i) ){ |
732 |
|
track.xm_a[i] = (double) xm[i] + segment; |
733 |
|
track.ym_a[i] = (double) ym[i] + omega * segment; |
734 |
|
// track.zm_a[i] = zm[i] - gamma * segment;//not used yet |
735 |
|
track.xm_b[i] = (double) xm[i] - segment; |
736 |
|
track.ym_b[i] = (double) ym[i] - omega * segment; |
737 |
|
// track.zm_b[i] = zm[i] + gamma * segment;//not used yet |
738 |
|
} |
739 |
|
} |
740 |
|
|
741 |
|
track.resx[i]= (double) resx[i]; |
742 |
|
track.resy[i]= (double) resy[i]; |
743 |
|
track.tailx[i]= (double) tailx[i]; |
744 |
|
track.taily[i]= (double) taily[i]; |
745 |
|
} |
746 |
|
|
747 |
|
for(int i=0; i<5; i++) track.al[i]= (double) al[i]; |
748 |
|
track.zini = 23.5; |
749 |
|
// ZINI = 23.5 !!! it should be the same parameter in all codes |
750 |
|
|
751 |
|
} |
752 |
|
/** |
753 |
|
* Method to set values from minimization-routine common |
754 |
|
*/ |
755 |
|
void TrkTrack::SetFromMiniStruct(cMini2track *track){ |
756 |
|
|
757 |
|
for(int i=0; i<5; i++) { |
758 |
|
al[i]= (float) (track->al[i]); |
759 |
|
for(int j=0; j<5; j++) coval[i][j]= (float) (track->cov[i][j]); |
760 |
|
} |
761 |
|
chi2 = (float) (track->chi2); |
762 |
|
nstep = (float) (track->nstep); |
763 |
|
for(int i=0; i<6; i++){ |
764 |
|
xv[i] = (float) (track->xv[i]); |
765 |
|
yv[i] = (float) (track->yv[i]); |
766 |
|
zv[i] = (float) (track->zv[i]); |
767 |
|
xm[i] = (float) (track->xm[i]); |
768 |
|
ym[i] = (float) (track->ym[i]); |
769 |
|
zm[i] = (float) (track->zm[i]); |
770 |
|
axv[i] = (float) (track->axv[i]); |
771 |
|
ayv[i] = (float) (track->ayv[i]); |
772 |
|
resx[i] = (float) (track->resx[i]); //Elena 10th |
773 |
|
resy[i] = (float) (track->resy[i]); |
774 |
|
} |
775 |
|
|
776 |
|
} |
777 |
|
/** |
778 |
|
* \brief Method to re-evaluate coordinates of clusters associated with a track. |
779 |
|
* |
780 |
|
* The method can be applied only after recovering level1 information |
781 |
|
* (either by reprocessing single events from level0 or from |
782 |
|
* the TrkLevel1 branch, if present); it calls F77 subroutines that |
783 |
|
* read the level1 common and fill the minimization-routine common. |
784 |
|
* Some clusters can be excluded or added by means of the methods: |
785 |
|
* |
786 |
|
* TrkTrack::ResetXGood(int ip) |
787 |
|
* TrkTrack::ResetYGood(int ip) |
788 |
|
* TrkTrack::SetXGood(int ip, int cid, int is) |
789 |
|
* TrkTrack::SetYGood(int ip, int cid, int is) |
790 |
|
* |
791 |
|
* NB! The method TrkTrack::SetGood(int *xg, int *yg) set the plane-mask (0-1) |
792 |
|
* for the minimization-routine common. It deletes the cluster information |
793 |
|
* (at least for the moment...) thus cannot be applied before |
794 |
|
* TrkTrack::EvaluateClusterPositions(). |
795 |
|
* |
796 |
|
* Different p.f.a. can be applied by calling (once) the method: |
797 |
|
* |
798 |
|
* TrkParams::SetPFA(0); //Set ETA p.f.a. |
799 |
|
* |
800 |
|
* @see TrkParams::SetPFA(int) |
801 |
|
*/ |
802 |
|
Bool_t TrkTrack::EvaluateClusterPositions(){ |
803 |
|
|
804 |
|
// cout << "void TrkTrack::GetClusterositions() "<<endl; |
805 |
|
|
806 |
|
TrkParams::Load(1); |
807 |
|
if( !TrkParams::IsLoaded(1) ){ |
808 |
|
cout << "Bool_t TrkTrack::EvaluateClusterPositions() ---ERROR--- m.field not loaded "<<endl; |
809 |
|
return false; |
810 |
|
} |
811 |
|
TrkParams::Load(4); |
812 |
|
if( !TrkParams::IsLoaded(4) ){ |
813 |
|
cout << "Bool_t TrkTrack::EvaluateClusterPositions() ---ERROR--- p.f.a. par. not loaded "<<endl; |
814 |
|
return false; |
815 |
|
} |
816 |
|
TrkParams::Load(5); |
817 |
|
if( !TrkParams::IsLoaded(5) ){ |
818 |
|
cout << "Bool_t TrkTrack::EvaluateClusterPositions() ---ERROR--- alignment par. not loaded "<<endl; |
819 |
|
return false; |
820 |
|
} |
821 |
|
|
822 |
|
for(int ip=0; ip<6; ip++){ |
823 |
|
// cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;; |
824 |
|
int icx = GetClusterX_ID(ip)+1;//0=no-cluster,1-N |
825 |
|
int icy = GetClusterY_ID(ip)+1;//0=no-cluster,1-N |
826 |
|
int sensor = GetSensor(ip)+1;//<< convenzione "Paolo" |
827 |
|
if(ip==5 && sensor!=0)sensor=3-sensor;//<< convenzione "Elena" |
828 |
|
int ladder = GetLadder(ip)+1; |
829 |
|
float ax = axv[ip]; |
830 |
|
float ay = ayv[ip]; |
831 |
|
float v[3]; |
832 |
|
v[0]=xv[ip]; |
833 |
|
v[1]=yv[ip]; |
834 |
|
v[2]=zv[ip]; |
835 |
|
float bfx = 10*TrkParams::GetBX(v);//Tesla |
836 |
|
float bfy = 10*TrkParams::GetBY(v);//Tesla |
837 |
|
int ipp=ip+1; |
838 |
|
xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy); |
839 |
|
// if(icx<0 || icy<0)return false; |
840 |
|
} |
841 |
|
return true; |
842 |
|
} |
843 |
|
/** |
844 |
|
* \brief Tracking method. It calls F77 mini routine. |
845 |
|
* |
846 |
|
* @param pfixed Particle momentum. If pfixed=0 the momentum |
847 |
|
* is left as a free parameter, otherwise it is fixed to the input value. |
848 |
|
* @param fail Output flag (!=0 if the fit failed). |
849 |
|
* @param iprint Flag to set debug mode ( 0 = no output; 1 = verbose; 2 = debug). |
850 |
|
* @param froml1 Flag to re-evaluate positions (see TrkTrack::GetClusterPositions()). |
851 |
|
* |
852 |
|
* The option to re-evaluate positions can be used only after recovering |
853 |
|
* level1 information, eg. by reprocessing the single event. |
854 |
|
* |
855 |
|
* Example: |
856 |
|
* |
857 |
|
* if( !event->GetTrkLevel0() )return false; |
858 |
|
* event->GetTrkLevel0()->ProcessEvent(); // re-processing level0->level1 |
859 |
|
* int fail=0; |
860 |
|
* event->GetTrkLevel2()->GetTrack(0)->Fit(0.,fail,0,1); |
861 |
|
* |
862 |
|
* @see EvaluateClusterPositions() |
863 |
|
* |
864 |
|
* The fitting procedure can be varied by changing the tracking mode, |
865 |
|
* the fit-precision factor, the minimum number of step, etc. |
866 |
|
* @see SetTrackingMode(int) |
867 |
|
* @see SetPrecisionFactor(double) |
868 |
|
* @see SetStepMin(int) |
869 |
|
* @see SetDeltaB(int,double) |
870 |
|
*/ |
871 |
|
void TrkTrack::Fit(double pfixed, int& fail, int iprint, int froml1){ |
872 |
|
|
873 |
|
TrkParams::Load(1); |
874 |
|
if( !TrkParams::IsLoaded(1) ){ |
875 |
|
cout << "void TrkTrack::Fit(double,int&,int,int) ---ERROR--- m.field not loaded "<<endl; |
876 |
|
return; |
877 |
|
} |
878 |
|
TrkParams::Load(5); |
879 |
|
if( !TrkParams::IsLoaded(5) ){ |
880 |
|
cout << "void TrkTrack::Fit(double,int&,int,int) ---ERROR--- align.param. not loaded "<<endl; |
881 |
|
return; |
882 |
|
} |
883 |
|
|
884 |
|
float al_ini[] = {0.,0.,0.,0.,0.}; |
885 |
|
|
886 |
|
extern cMini2track track_; |
887 |
|
fail = 0; |
888 |
|
|
889 |
|
// FillMiniStruct(track_); |
890 |
|
|
891 |
|
if(froml1!=0){ |
892 |
|
if( !EvaluateClusterPositions() ){ |
893 |
|
cout << "void TrkTrack::Fit("<<pfixed<<","<<fail<<","<<iprint<<","<<froml1<<") --- ERROR evaluating cluster positions "<<endl; |
894 |
|
FillMiniStruct(track_) ; |
895 |
|
fail = 1; |
896 |
|
return; |
897 |
|
} |
898 |
|
}else{ |
899 |
|
FillMiniStruct(track_); |
900 |
|
} |
901 |
|
|
902 |
|
// if fit variables have been reset, evaluate the initial guess |
903 |
|
if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_(); |
904 |
|
|
905 |
|
// --------------------- free momentum |
906 |
|
if(pfixed==0.) { |
907 |
|
track_.pfixed=0.; |
908 |
|
} |
909 |
|
// --------------------- fixed momentum |
910 |
|
if(pfixed!=0.) { |
911 |
|
al[4]=1./pfixed; |
912 |
|
track_.pfixed=pfixed; |
913 |
|
} |
914 |
|
|
915 |
|
// store temporarily the initial guess |
916 |
|
for(int i=0; i<5; i++) al_ini[i]=track_.al[i]; |
917 |
|
|
918 |
|
// ------------------------------------------ |
919 |
|
// call mini routine |
920 |
|
// ------------------------------------------ |
921 |
|
int istep=0; |
922 |
|
int ifail=0; |
923 |
|
mini2_(&istep,&ifail, &iprint); |
924 |
|
if(ifail!=0) { |
925 |
|
if(iprint)cout << "ERROR: ifail= " << ifail << endl; |
926 |
|
fail = 1; |
927 |
|
} |
928 |
|
if(chi2!=chi2){ |
929 |
|
if(iprint)cout << "ERROR: chi2= " << chi2 << endl; |
930 |
|
FitReset(); |
931 |
|
fail = 1; |
932 |
|
} |
933 |
|
// ------------------------------------------ |
934 |
|
|
935 |
|
SetFromMiniStruct(&track_); |
936 |
|
|
937 |
|
if(fail){ |
938 |
|
if(iprint)cout << " >>>> fit failed "<<endl; |
939 |
|
for(int i=0; i<5; i++) al[i]=al_ini[i]; |
940 |
|
} |
941 |
|
|
942 |
|
|
943 |
|
//ELENA 2015 |
944 |
|
int ngf = TrkParams::nGF; |
945 |
|
float* zgf = TrkParams::zGF; |
946 |
|
Trajectory tj = Trajectory(ngf,zgf); |
947 |
|
tj.DoTrack(al,ZINI); |
948 |
|
for(int i=0; i<14; i++){ |
949 |
|
xGF[i] = tj.x[i]; |
950 |
|
yGF[i] = tj.y[i]; |
951 |
|
} |
952 |
|
if(ngf!=14){ |
953 |
|
cout << "TrkTrack::Fit("<<pfixed<<","<<fail<<","<<iprint<<","<<froml1<<"<<) "<<endl;; |
954 |
|
cout << "TrkParams::nGF = "<<TrkParams::nGF<<" != 14 "<<endl; |
955 |
|
cout << "back-incompatibility !?!?! acceptance check not reliable"<<endl; |
956 |
|
} |
957 |
|
|
958 |
|
|
959 |
|
|
960 |
|
}; |
961 |
|
/** |
962 |
|
* Reset the fit parameters |
963 |
|
*/ |
964 |
|
void TrkTrack::FitReset(){ |
965 |
|
for(int i=0; i<5; i++) al[i]=-9999.; |
966 |
|
chi2=0.; |
967 |
|
nstep=0; |
968 |
|
// for(int i=0; i<6; i++) xv[i]=0.; |
969 |
|
// for(int i=0; i<6; i++) yv[i]=0.; |
970 |
|
// for(int i=0; i<6; i++) zv[i]=0.; |
971 |
|
// for(int i=0; i<6; i++) axv[i]=0.; |
972 |
|
// for(int i=0; i<6; i++) ayv[i]=0.; |
973 |
|
for(int i=0; i<5; i++) { |
974 |
|
for(int j=0; j<5; j++) coval[i][j]=0.; |
975 |
} |
} |
976 |
|
} |
977 |
|
/** |
978 |
|
* Set the tracking mode |
979 |
|
*/ |
980 |
|
void TrkTrack::SetTrackingMode(int trackmode){ |
981 |
|
extern cMini2track track_; |
982 |
|
track_.trackmode = trackmode; |
983 |
|
} |
984 |
|
/** |
985 |
|
* Set the factor scale for tracking precision |
986 |
|
*/ |
987 |
|
void TrkTrack::SetPrecisionFactor(double fact){ |
988 |
|
extern cMini2track track_; |
989 |
|
track_.fact = fact; |
990 |
|
} |
991 |
|
/** |
992 |
|
* Set the minimum number of steps for tracking precision |
993 |
|
*/ |
994 |
|
void TrkTrack::SetStepMin(int istepmin){ |
995 |
|
extern cMini2track track_; |
996 |
|
track_.istepmin = istepmin; |
997 |
|
} |
998 |
|
/** |
999 |
|
* Set deltaB parameters (id=0,1). By default they are set to zero. |
1000 |
|
*/ |
1001 |
|
void TrkTrack::SetDeltaB(int id, double db){ |
1002 |
|
if(id!=0 && id!=1)cout << "void TrkTrack::SetDeltaB(int id,double db) -- wrong input parameters: "<<id<<" "<<db<<endl; |
1003 |
|
TrkParams::SetDeltaB(id,db); |
1004 |
|
} |
1005 |
|
|
1006 |
// delete [] dxout; |
/** |
1007 |
// delete [] dyout; |
* Returns true if the track is inside the magnet cavity. |
1008 |
// delete [] dzin; |
* @param toll Tolerance around the nominal volume (toll>0 define an inner fiducial volume) |
1009 |
|
*/ |
1010 |
|
Bool_t TrkTrack::IsInsideCavity(float toll){ |
1011 |
|
|
1012 |
|
// float xmagntop, ymagntop, xmagnbottom, ymagnbottom; |
1013 |
|
// xmagntop = xv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*axv[0]/180.); |
1014 |
|
// ymagntop = yv[0] + (ZMAGNHIGH-zv[0])*tan(acos(-1.0)*ayv[0]/180.); |
1015 |
|
// xmagnbottom = xv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*axv[5]/180.); |
1016 |
|
// ymagnbottom = yv[5] + (ZMAGNLOW-zv[5])*tan(acos(-1.0)*ayv[5]/180.); |
1017 |
|
// if( xmagntop>XMAGNLOW && xmagntop<XMAGNHIGH && |
1018 |
|
// ymagntop>YMAGNLOW && ymagntop<YMAGNHIGH && |
1019 |
|
// xmagnbottom>XMAGNLOW && xmagnbottom<XMAGNHIGH && |
1020 |
|
// ymagnbottom>YMAGNLOW && ymagnbottom<YMAGNHIGH ) return(true); |
1021 |
|
// else return(false); |
1022 |
|
|
1023 |
|
int ngf = TrkParams::nGF; |
1024 |
|
for(int i=0; i<ngf; i++){ |
1025 |
|
// |
1026 |
|
// cout << endl << TrkParams::GF_element[i]; |
1027 |
|
if( |
1028 |
|
TrkParams::GF_element[i].CompareTo("CUF") && |
1029 |
|
TrkParams::GF_element[i].CompareTo("T2") && |
1030 |
|
TrkParams::GF_element[i].CompareTo("T3") && |
1031 |
|
TrkParams::GF_element[i].CompareTo("T4") && |
1032 |
|
TrkParams::GF_element[i].CompareTo("T5") && |
1033 |
|
TrkParams::GF_element[i].CompareTo("CLF") && |
1034 |
|
true)continue; |
1035 |
|
// apply condition only within the cavity |
1036 |
|
// cout << " -- "<<xGF[i]<<" "<<yGF[i]; |
1037 |
|
if( |
1038 |
|
xGF[i] <= TrkParams::xGF_min[i] + toll || |
1039 |
|
xGF[i] >= TrkParams::xGF_max[i] - toll || |
1040 |
|
yGF[i] <= TrkParams::yGF_min[i] + toll || |
1041 |
|
yGF[i] >= TrkParams::yGF_max[i] - toll || |
1042 |
|
false){ |
1043 |
|
|
1044 |
|
return false; |
1045 |
|
} |
1046 |
|
} |
1047 |
|
return true; |
1048 |
|
|
1049 |
|
|
1050 |
|
} |
1051 |
|
/** |
1052 |
|
* Returns true if the track is inside the nominal acceptance, which is defined |
1053 |
|
* by the intersection among magnet cavity, silicon-plane sensitive area and |
1054 |
|
* ToF sensitive area (nominal values from the official document used to |
1055 |
|
* calculate the geometrical factor) |
1056 |
|
* @param toll Tolerance around the nominal volume (toll>0 define an inner fiducial volume) |
1057 |
|
*/ |
1058 |
|
// Bool_t TrkTrack::IsInsideAcceptance(){ |
1059 |
|
|
1060 |
|
// int ngf = TrkParams::nGF; |
1061 |
|
// for(int i=0; i<ngf; i++){ |
1062 |
|
// if( |
1063 |
|
// xGF[i] <= TrkParams::xGF_min[i] || |
1064 |
|
// xGF[i] >= TrkParams::xGF_max[i] || |
1065 |
|
// yGF[i] <= TrkParams::yGF_min[i] || |
1066 |
|
// yGF[i] >= TrkParams::yGF_max[i] || |
1067 |
|
// false)return false; |
1068 |
|
// } |
1069 |
|
// return true; |
1070 |
|
|
1071 |
|
// } |
1072 |
|
Bool_t TrkTrack::IsInsideAcceptance(float toll){ |
1073 |
|
|
1074 |
|
|
1075 |
|
int ngf = TrkParams::nGF; |
1076 |
|
for(int i=0; i<ngf; i++){ |
1077 |
|
// |
1078 |
|
// cout << endl << TrkParams::GF_element[i]; |
1079 |
|
if( |
1080 |
|
TrkParams::GF_element[i].CompareTo("S11") && |
1081 |
|
TrkParams::GF_element[i].CompareTo("S12") && |
1082 |
|
TrkParams::GF_element[i].CompareTo("S21") && |
1083 |
|
TrkParams::GF_element[i].CompareTo("S22") && |
1084 |
|
TrkParams::GF_element[i].CompareTo("T1") && |
1085 |
|
TrkParams::GF_element[i].CompareTo("CUF") && |
1086 |
|
TrkParams::GF_element[i].CompareTo("T2") && |
1087 |
|
TrkParams::GF_element[i].CompareTo("T3") && |
1088 |
|
TrkParams::GF_element[i].CompareTo("T4") && |
1089 |
|
TrkParams::GF_element[i].CompareTo("T5") && |
1090 |
|
TrkParams::GF_element[i].CompareTo("CLF") && |
1091 |
|
TrkParams::GF_element[i].CompareTo("T6") && |
1092 |
|
TrkParams::GF_element[i].CompareTo("S31") && |
1093 |
|
TrkParams::GF_element[i].CompareTo("S32") && |
1094 |
|
true)continue; |
1095 |
|
// apply condition only within the cavity |
1096 |
|
// cout << " -- "<<xGF[i]<<" "<<yGF[i]; |
1097 |
|
if( |
1098 |
|
xGF[i] <= TrkParams::xGF_min[i] + toll || |
1099 |
|
xGF[i] >= TrkParams::xGF_max[i] - toll || |
1100 |
|
yGF[i] <= TrkParams::yGF_min[i] + toll || |
1101 |
|
yGF[i] >= TrkParams::yGF_max[i] - toll || |
1102 |
|
false){ |
1103 |
|
|
1104 |
|
return false; |
1105 |
|
} |
1106 |
|
} |
1107 |
|
return true; |
1108 |
|
} |
1109 |
|
|
1110 |
|
/** |
1111 |
|
* Returns true if the track is inside one of the surfaces which define the |
1112 |
|
* geometrical acceptance. |
1113 |
|
* @param surf tag of the surface (possible values are: S11 S12 S21 S22 T1 |
1114 |
|
* CUF T2 T3 T4 T5 CLF T6 S31 S32). |
1115 |
|
* @param toll Tolerance around the nominal surface (toll>0 define an inner |
1116 |
|
* fiducial surface) |
1117 |
|
*/ |
1118 |
|
Bool_t TrkTrack::IsInsideGFSurface(const char* surf, float toll){ |
1119 |
|
|
1120 |
|
|
1121 |
|
int ngf = TrkParams::nGF; |
1122 |
|
bool SURFOK = false; |
1123 |
|
for(int i=0; i<ngf; i++){ |
1124 |
|
if( !TrkParams::GF_element[i].CompareTo(surf) ){ |
1125 |
|
SURFOK=true; |
1126 |
|
if( |
1127 |
|
xGF[i] > TrkParams::xGF_min[i] + toll && |
1128 |
|
xGF[i] < TrkParams::xGF_max[i] - toll && |
1129 |
|
yGF[i] > TrkParams::yGF_min[i] + toll && |
1130 |
|
yGF[i] < TrkParams::yGF_max[i] - toll && |
1131 |
|
true)return true; |
1132 |
|
} |
1133 |
|
} |
1134 |
|
if( !SURFOK )cout << " Bool_t TrkTrack::IsInsideGFSurface(char* surf, float toll) --> suface "<<surf<<" not defined "<<endl; |
1135 |
|
return false; |
1136 |
|
|
1137 |
|
} |
1138 |
|
|
1139 |
|
/** |
1140 |
|
* Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track. |
1141 |
|
* If no cluster is associated, ID=-1. |
1142 |
|
* @param ip Tracker plane (0-5) |
1143 |
|
*/ |
1144 |
|
Int_t TrkTrack::GetClusterX_ID(int ip){ |
1145 |
|
return ((Int_t)fabs(xgood[ip]))%10000000-1; |
1146 |
|
}; |
1147 |
|
/** |
1148 |
|
* Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track. |
1149 |
|
* If no cluster is associated, ID=-1. |
1150 |
|
* @param ip Tracker plane (0-5) |
1151 |
|
*/ |
1152 |
|
Int_t TrkTrack::GetClusterY_ID(int ip){ |
1153 |
|
return ((Int_t)fabs(ygood[ip]))%10000000-1; |
1154 |
|
}; |
1155 |
|
|
1156 |
|
/** |
1157 |
|
* Method to retrieve the ladder (0-2, increasing x) traversed by the track on this plane. |
1158 |
|
* If no ladder is traversed (dead area) the metod retuns -1. |
1159 |
|
* @param ip Tracker plane (0-5) |
1160 |
|
*/ |
1161 |
|
Int_t TrkTrack::GetLadder(int ip){ |
1162 |
|
if(XGood(ip))return (Int_t)fabs(xgood[ip]/100000000)-1; |
1163 |
|
if(YGood(ip))return (Int_t)fabs(ygood[ip]/100000000)-1; |
1164 |
|
return -1; |
1165 |
|
}; |
1166 |
|
/** |
1167 |
|
* Method to retrieve the sensor (0-1, increasing y) traversed by the track on this plane. |
1168 |
|
* If no sensor is traversed (dead area) the metod retuns -1. |
1169 |
|
* @param ip Tracker plane (0-5) |
1170 |
|
*/ |
1171 |
|
Int_t TrkTrack::GetSensor(int ip){ |
1172 |
|
if(XGood(ip))return (Int_t)((Int_t)fabs(xgood[ip]/10000000)%10)-1; |
1173 |
|
if(YGood(ip))return (Int_t)((Int_t)fabs(ygood[ip]/10000000)%10)-1; |
1174 |
|
return -1; |
1175 |
|
}; |
1176 |
|
|
1177 |
|
/** |
1178 |
|
* \brief Method to include a x-cluster to the track. |
1179 |
|
* @param ip Tracker plane (0-5) |
1180 |
|
* @param clid Cluster ID (0 = no-cluster, 1,2,... otherwise ) |
1181 |
|
* @param il Ladder (0-2, increasing x, -1 if no sensitive area is hit) |
1182 |
|
* @param is Sensor (0-1, increasing y, -1 if no sensitive area is hit) |
1183 |
|
* @param bad True if the cluster contains bad strips |
1184 |
|
* @see Fit(double pfixed, int& fail, int iprint, int froml1) |
1185 |
|
*/ |
1186 |
|
void TrkTrack::SetXGood(int ip, int clid, int il, int is, bool bad){ |
1187 |
|
// int il=0; //ladder (temporary) |
1188 |
|
// bool bad=false; //ladder (temporary) |
1189 |
|
if(ip<0||ip>5||clid<1||il<-1||il>2||is<-1||is>1) |
1190 |
|
cout << " void TrkTrack::SetXGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl; |
1191 |
|
xgood[ip]=(il+1)*100000000+(is+1)*10000000+clid; |
1192 |
|
if(bad)xgood[ip]=-xgood[ip]; |
1193 |
|
}; |
1194 |
|
/** |
1195 |
|
* \brief Method to include a y-cluster to the track. |
1196 |
|
* @param ip Tracker plane (0-5) |
1197 |
|
* @param clid Cluster ID (0 = no-cluster, 1,2,... otherwise ) |
1198 |
|
* @param il Ladder (0-2, increasing x, -1 if no sensitive area is hit) |
1199 |
|
* @param is Sensor (0-1, increasing y, -1 if no sensitive area is hit) |
1200 |
|
* @param bad True if the cluster contains bad strips |
1201 |
|
* @see Fit(double pfixed, int& fail, int iprint, int froml1) |
1202 |
|
*/ |
1203 |
|
void TrkTrack::SetYGood(int ip, int clid, int il, int is, bool bad){ |
1204 |
|
// int il=0; //ladder (temporary) |
1205 |
|
// bool bad=false; //ladder (temporary) |
1206 |
|
if(ip<0||ip>5||clid<1||il<-1||il>2||is<-1||is>1) |
1207 |
|
cout << " void TrkTrack::SetYGood(int,int,int,int,bool) --> MA SEI DI COCCIO?!?!"<<endl; |
1208 |
|
ygood[ip]=(il+1)*100000000+(is+1)*10000000+clid; |
1209 |
|
if(bad)ygood[ip]=-ygood[ip]; |
1210 |
|
}; |
1211 |
|
|
1212 |
|
/** |
1213 |
|
* \brief Average X |
1214 |
|
* Average value of <xv>, evaluated from the first to the last hit x view. |
1215 |
|
*/ |
1216 |
|
Float_t TrkTrack::GetXav(){ |
1217 |
|
|
1218 |
|
int first_plane = -1; |
1219 |
|
int last_plane = -1; |
1220 |
|
for(Int_t ip=0; ip<6; ip++){ |
1221 |
|
if( XGood(ip) && first_plane == -1 )first_plane = ip; |
1222 |
|
if( XGood(ip) && first_plane != -1 )last_plane = ip; |
1223 |
|
} |
1224 |
|
if( first_plane == -1 || last_plane == -1){ |
1225 |
|
return -100; |
1226 |
|
} |
1227 |
|
if( last_plane-first_plane+1 ==0 )return -100; |
1228 |
|
|
1229 |
|
Float_t av = 0; |
1230 |
|
for(int ip=first_plane; ip<=last_plane; ip++)av+=xv[ip]; |
1231 |
|
|
1232 |
|
return (av/(last_plane-first_plane+1)); |
1233 |
|
} |
1234 |
|
/** |
1235 |
|
* \brief Average Y |
1236 |
|
* Average value of <yv>, evaluated from the first to the last hit x view. |
1237 |
|
*/ |
1238 |
|
Float_t TrkTrack::GetYav(){ |
1239 |
|
|
1240 |
|
int first_plane = -1; |
1241 |
|
int last_plane = -1; |
1242 |
|
for(Int_t ip=0; ip<6; ip++){ |
1243 |
|
if( XGood(ip) && first_plane == -1 )first_plane = ip; |
1244 |
|
if( XGood(ip) && first_plane != -1 )last_plane = ip; |
1245 |
|
} |
1246 |
|
if( first_plane == -1 || last_plane == -1){ |
1247 |
|
return -100; |
1248 |
|
} |
1249 |
|
if( last_plane-first_plane+1 ==0 )return -100; |
1250 |
|
|
1251 |
|
Float_t av = 0; |
1252 |
|
for(int ip=first_plane; ip<=last_plane; ip++)av+=yv[ip]; |
1253 |
|
|
1254 |
|
return (av/(last_plane-first_plane+1)); |
1255 |
|
} |
1256 |
|
/** |
1257 |
|
* \brief Average Z |
1258 |
|
* Average value of <zv>, evaluated from the first to the last hit x view. |
1259 |
|
*/ |
1260 |
|
Float_t TrkTrack::GetZav(){ |
1261 |
|
|
1262 |
|
int first_plane = -1; |
1263 |
|
int last_plane = -1; |
1264 |
|
for(Int_t ip=0; ip<6; ip++){ |
1265 |
|
if( XGood(ip) && first_plane == -1 )first_plane = ip; |
1266 |
|
if( XGood(ip) && first_plane != -1 )last_plane = ip; |
1267 |
|
} |
1268 |
|
if( first_plane == -1 || last_plane == -1){ |
1269 |
|
return -100; |
1270 |
|
} |
1271 |
|
if( last_plane-first_plane+1 ==0 )return -100; |
1272 |
|
|
1273 |
|
Float_t av = 0; |
1274 |
|
for(int ip=first_plane; ip<=last_plane; ip++)av+=zv[ip]; |
1275 |
|
|
1276 |
|
return (av/(last_plane-first_plane+1)); |
1277 |
|
} |
1278 |
|
|
1279 |
|
/** |
1280 |
|
* \brief Number of column traversed |
1281 |
|
*/ |
1282 |
|
Int_t TrkTrack::GetNColumns(){ |
1283 |
|
int sensors[] = {0,0,0,0,0,0}; |
1284 |
|
for(int ip=0; ip<6; ip++){ |
1285 |
|
int sensorid = GetLadder(ip)+3*GetSensor(ip); |
1286 |
|
if(XGood(ip)||YGood(ip)) |
1287 |
|
if(sensorid>=0 && sensorid<6)sensors[sensorid]=1; |
1288 |
|
} |
1289 |
|
int nsensors=0; |
1290 |
|
for(int is=0; is<6; is++)nsensors += sensors[is]; |
1291 |
|
return nsensors; |
1292 |
|
}; |
1293 |
|
/** |
1294 |
|
* \brief Give the maximum energy release |
1295 |
|
*/ |
1296 |
|
Float_t TrkTrack::GetDEDX_max(int ip, int iv){ |
1297 |
|
Float_t max=0; |
1298 |
|
int pfrom = 0; |
1299 |
|
int pto = 6; |
1300 |
|
int vfrom = 0; |
1301 |
|
int vto = 2; |
1302 |
|
if(ip>=0&&ip<6){ |
1303 |
|
pfrom = ip; |
1304 |
|
pto = ip+1; |
1305 |
|
} |
1306 |
|
if(iv>=0&&iv<2){ |
1307 |
|
vfrom = iv; |
1308 |
|
vto = iv+1; |
1309 |
|
} |
1310 |
|
for(int i=pfrom; i<pto; i++) |
1311 |
|
for(int j=vfrom; j<vto; j++){ |
1312 |
|
if(j==0 && XGood(i) && GetDEDX(i,j)>max)max=GetDEDX(i,j); |
1313 |
|
if(j==1 && YGood(i) && GetDEDX(i,j)>max)max=GetDEDX(i,j); |
1314 |
|
} |
1315 |
|
return max; |
1316 |
|
|
1317 |
|
}; |
1318 |
|
|
1319 |
|
/** |
1320 |
|
* \brief Give the minimum energy release |
1321 |
|
*/ |
1322 |
|
Float_t TrkTrack::GetDEDX_min(int ip, int iv){ |
1323 |
|
Float_t min=100000000; |
1324 |
|
int pfrom = 0; |
1325 |
|
int pto = 6; |
1326 |
|
int vfrom = 0; |
1327 |
|
int vto = 2; |
1328 |
|
if(ip>=0&&ip<6){ |
1329 |
|
pfrom = ip; |
1330 |
|
pto = ip+1; |
1331 |
|
} |
1332 |
|
if(iv>=0&&iv<2){ |
1333 |
|
vfrom = iv; |
1334 |
|
vto = iv+1; |
1335 |
|
} |
1336 |
|
for(int i=pfrom; i<pto; i++) |
1337 |
|
for(int j=vfrom; j<vto; j++){ |
1338 |
|
if(j==0 && XGood(i) && GetDEDX(i,j)<min)min=GetDEDX(i,j); |
1339 |
|
if(j==1 && YGood(i) && GetDEDX(i,j)<min)min=GetDEDX(i,j); |
1340 |
|
} |
1341 |
|
return min; |
1342 |
|
|
1343 |
|
}; |
1344 |
|
|
1345 |
|
/** |
1346 |
|
* \brief Give the maximum spatial residual |
1347 |
|
*/ |
1348 |
|
Float_t TrkTrack::GetResidual_max(int ip, int iv){ |
1349 |
|
Float_t max=0; |
1350 |
|
int pfrom = 0; |
1351 |
|
int pto = 6; |
1352 |
|
int vfrom = 0; |
1353 |
|
int vto = 2; |
1354 |
|
if(ip>=0&&ip<6){ |
1355 |
|
pfrom = ip; |
1356 |
|
pto = ip+1; |
1357 |
|
} |
1358 |
|
if(iv>=0&&iv<2){ |
1359 |
|
vfrom = iv; |
1360 |
|
vto = iv+1; |
1361 |
|
} |
1362 |
|
for(int i=pfrom; i<pto; i++){ |
1363 |
|
for(int j=vfrom; j<vto; j++){ |
1364 |
|
if(j==0 && XGood(i) && fabs(xm[i]-xv[i])>fabs(max))max=xm[i]-xv[i]; |
1365 |
|
if(j==1 && YGood(i) && fabs(ym[i]-yv[i])>fabs(max))max=ym[i]-yv[i]; |
1366 |
|
} |
1367 |
|
} |
1368 |
|
return max; |
1369 |
|
|
1370 |
|
}; |
1371 |
|
/** |
1372 |
|
* \brief Give the anerage spatial residual |
1373 |
|
*/ |
1374 |
|
Float_t TrkTrack::GetResidual_av(int ip, int iv){ |
1375 |
|
// |
1376 |
|
//Sum$((xm>-50)*(xm-xv)/resx)/sqrt(TrkTrack.GetNX()*TrkTrack.GetChi2X())<0.3 |
1377 |
|
|
1378 |
|
Float_t av = 0.; |
1379 |
|
int nav = 0; |
1380 |
|
// |
1381 |
|
int pfrom = 0; |
1382 |
|
int pto = 6; |
1383 |
|
int vfrom = 0; |
1384 |
|
int vto = 2; |
1385 |
|
if(ip>=0&&ip<6){ |
1386 |
|
pfrom = ip; |
1387 |
|
pto = ip+1; |
1388 |
|
} |
1389 |
|
if(iv>=0&&iv<2){ |
1390 |
|
vfrom = iv; |
1391 |
|
vto = iv+1; |
1392 |
|
} |
1393 |
|
for(int i=pfrom; i<pto; i++){ |
1394 |
|
for(int j=vfrom; j<vto; j++){ |
1395 |
|
nav++; |
1396 |
|
if(j==0 && XGood(i)) av += (xm[i]-xv[i])/resx[i]; |
1397 |
|
if(j==1 && YGood(i)) av += (ym[i]-yv[i])/resy[i]; |
1398 |
|
} |
1399 |
|
} |
1400 |
|
if(nav==0)return -100.; |
1401 |
|
return av/nav; |
1402 |
|
|
1403 |
|
}; |
1404 |
|
|
1405 |
|
|
1406 |
|
/** |
1407 |
|
* \brief Give the maximum multiplicity on the x view |
1408 |
|
*/ |
1409 |
|
Int_t TrkTrack::GetClusterX_Multiplicity_max(){ |
1410 |
|
int max=0; |
1411 |
|
for(int ip=0; ip<6; ip++) |
1412 |
|
if(GetClusterX_Multiplicity(ip)>max)max=GetClusterX_Multiplicity(ip); |
1413 |
|
return max; |
1414 |
|
}; |
1415 |
|
/** |
1416 |
|
* \brief Give the minimum multiplicity on the x view |
1417 |
|
*/ |
1418 |
|
Int_t TrkTrack::GetClusterX_Multiplicity_min(){ |
1419 |
|
int min=50; |
1420 |
|
for(int ip=0; ip<6; ip++) |
1421 |
|
if(GetClusterX_Multiplicity(ip)<min)min=GetClusterX_Multiplicity(ip); |
1422 |
|
return min; |
1423 |
|
}; |
1424 |
|
/** |
1425 |
|
* \brief Give the maximum multiplicity on the x view |
1426 |
|
*/ |
1427 |
|
Int_t TrkTrack::GetClusterY_Multiplicity_max(){ |
1428 |
|
int max=0; |
1429 |
|
for(int ip=0; ip<6; ip++) |
1430 |
|
if(GetClusterY_Multiplicity(ip)>max)max=GetClusterY_Multiplicity(ip); |
1431 |
|
return max; |
1432 |
|
}; |
1433 |
|
/** |
1434 |
|
* \brief Give the minimum multiplicity on the x view |
1435 |
|
*/ |
1436 |
|
Int_t TrkTrack::GetClusterY_Multiplicity_min(){ |
1437 |
|
int min=50; |
1438 |
|
for(int ip=0; ip<6; ip++) |
1439 |
|
if(GetClusterY_Multiplicity(ip)<min)min=GetClusterY_Multiplicity(ip); |
1440 |
|
return min; |
1441 |
|
}; |
1442 |
|
|
1443 |
return ifail; |
/** |
1444 |
|
* \brief Give the minimum seed on the x view |
1445 |
|
*/ |
1446 |
|
Float_t TrkTrack::GetClusterX_Seed_min(){ |
1447 |
|
Float_t min=100000; |
1448 |
|
for(int ip=0; ip<6; ip++) |
1449 |
|
if(XGood(ip) && GetClusterX_Seed(ip)<min)min=GetClusterX_Seed(ip); |
1450 |
|
return min; |
1451 |
|
}; |
1452 |
|
/** |
1453 |
|
* \brief Give the minimum seed on the x view |
1454 |
|
*/ |
1455 |
|
Float_t TrkTrack::GetClusterY_Seed_min(){ |
1456 |
|
Float_t min=100000; |
1457 |
|
for(int ip=0; ip<6; ip++) |
1458 |
|
if(YGood(ip) && GetClusterY_Seed(ip)<min)min=GetClusterY_Seed(ip); |
1459 |
|
return min; |
1460 |
}; |
}; |
1461 |
|
|
1462 |
|
|
1463 |
//-------------------------------------- |
//-------------------------------------- |
1464 |
// |
// |
1465 |
// |
// |
1466 |
//-------------------------------------- |
//-------------------------------------- |
1467 |
//float TrkTrack::BdL(){ |
void TrkTrack::Clear(){ |
1468 |
//}; |
// cout << "TrkTrack::Clear()"<<endl; |
1469 |
|
seqno = -1; |
1470 |
|
image = -1; |
1471 |
|
chi2 = 0; |
1472 |
|
nstep = 0; |
1473 |
|
for(int it1=0;it1<5;it1++){ |
1474 |
|
al[it1] = 0; |
1475 |
|
for(int it2=0;it2<5;it2++)coval[it1][it2] = 0; |
1476 |
|
}; |
1477 |
|
for(int ip=0;ip<6;ip++){ |
1478 |
|
xgood[ip] = 0; |
1479 |
|
ygood[ip] = 0; |
1480 |
|
xm[ip] = 0; |
1481 |
|
ym[ip] = 0; |
1482 |
|
zm[ip] = 0; |
1483 |
|
resx[ip] = 0; |
1484 |
|
resy[ip] = 0; |
1485 |
|
tailx[ip] = 0; |
1486 |
|
taily[ip] = 0; |
1487 |
|
xv[ip] = 0; |
1488 |
|
yv[ip] = 0; |
1489 |
|
zv[ip] = 0; |
1490 |
|
axv[ip] = 0; |
1491 |
|
ayv[ip] = 0; |
1492 |
|
dedx_x[ip] = 0; |
1493 |
|
dedx_y[ip] = 0; |
1494 |
|
|
1495 |
|
}; |
1496 |
|
int ngf = TrkParams::nGF; |
1497 |
|
for(int i=0; i<ngf; i++){ |
1498 |
|
xGF[i] = 0.; |
1499 |
|
yGF[i] = 0.; |
1500 |
|
} |
1501 |
|
// if(clx)clx->Clear(); |
1502 |
|
// if(cly)cly->Clear(); |
1503 |
|
// clx.Clear(); |
1504 |
|
// cly.Clear(); |
1505 |
|
}; |
1506 |
//-------------------------------------- |
//-------------------------------------- |
1507 |
// |
// |
1508 |
// |
// |
1509 |
//-------------------------------------- |
//-------------------------------------- |
1510 |
Float_t TrkTrack::GetRigidity(){ |
void TrkTrack::Delete(){ |
1511 |
Float_t rig=0; |
// cout << "TrkTrack::Delete()"<<endl; |
1512 |
if(chi2>0)rig=1./al[4]; |
Clear(); |
1513 |
if(rig<0) rig=-rig; |
// if(clx)delete clx; |
1514 |
return rig; |
// if(cly)delete cly; |
|
}; |
|
|
// |
|
|
Float_t TrkTrack::GetDeflection(){ |
|
|
Float_t def=0; |
|
|
if(chi2>0)def=al[4]; |
|
|
return def; |
|
|
}; |
|
|
// |
|
|
Float_t TrkTrack::GetDEDX(){ |
|
|
Float_t dedx=0; |
|
|
for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*xgood[i]+dedx_y[i]*ygood[i]; |
|
|
dedx = dedx/(this->GetNX()+this->GetNY()); |
|
|
return dedx; |
|
1515 |
}; |
}; |
|
|
|
1516 |
//-------------------------------------- |
//-------------------------------------- |
1517 |
// |
// |
1518 |
// |
// |
1519 |
//-------------------------------------- |
//-------------------------------------- |
1520 |
void TrkTrack::Dump(){ |
|
|
cout << endl << "========== Track " ; |
|
|
cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " "; |
|
|
cout << endl << "chi^2 : "<< chi2; |
|
|
cout << endl << "xgood : "; for(int i=0; i<6; i++)cout << xgood[i] ; |
|
|
cout << endl << "ygood : "; for(int i=0; i<6; i++)cout << ygood[i] ; |
|
|
cout << endl << "xm : "; for(int i=0; i<6; i++)cout << xm[i] << " "; |
|
|
cout << endl << "ym : "; for(int i=0; i<6; i++)cout << ym[i] << " "; |
|
|
cout << endl << "zm : "; for(int i=0; i<6; i++)cout << zm[i] << " "; |
|
|
cout << endl << "dedx_x : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " "; |
|
|
cout << endl << "dedx_y : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " "; |
|
|
} |
|
1521 |
//-------------------------------------- |
//-------------------------------------- |
1522 |
// |
// |
1523 |
// |
// |
1524 |
//-------------------------------------- |
//-------------------------------------- |
1525 |
TrkSinglet::TrkSinglet(){ |
TrkSinglet::TrkSinglet(){ |
1526 |
plane = 0; |
// cout << "TrkSinglet::TrkSinglet() " << GetUniqueID()<<endl; |
1527 |
coord[0] = 0; |
// plane = 0; |
1528 |
coord[1] = 0; |
// coord[0] = 0; |
1529 |
sgnl = 0; |
// coord[1] = 0; |
1530 |
|
// sgnl = 0; |
1531 |
|
// multmax = 0; |
1532 |
|
// cls = 0; |
1533 |
|
Clear(); |
1534 |
}; |
}; |
1535 |
//-------------------------------------- |
//-------------------------------------- |
1536 |
// |
// |
1537 |
// |
// |
1538 |
//-------------------------------------- |
//-------------------------------------- |
1539 |
TrkSinglet::TrkSinglet(const TrkSinglet& s){ |
TrkSinglet::TrkSinglet(const TrkSinglet& s){ |
1540 |
|
// cout << "TrkSinglet::TrkSinglet(const TrkSinglet& s) " << GetUniqueID()<<endl; |
1541 |
plane = s.plane; |
plane = s.plane; |
1542 |
coord[0] = s.coord[0]; |
coord[0] = s.coord[0]; |
1543 |
coord[1] = s.coord[1]; |
coord[1] = s.coord[1]; |
1544 |
sgnl = s.sgnl; |
sgnl = s.sgnl; |
1545 |
|
multmax = s.multmax; |
1546 |
|
// cls = 0;//<<<<pointer |
1547 |
|
// cls = TRef(s.cls); |
1548 |
}; |
}; |
1549 |
//-------------------------------------- |
//-------------------------------------- |
1550 |
// |
// |
1553 |
void TrkSinglet::Dump(){ |
void TrkSinglet::Dump(){ |
1554 |
int i=0; |
int i=0; |
1555 |
cout << endl << "========== Singlet " ; |
cout << endl << "========== Singlet " ; |
1556 |
cout << endl << "plane : " << plane; |
cout << endl << "plane : " << plane; |
1557 |
cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++; |
cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++; |
1558 |
cout << endl << "sgnl : " << sgnl; |
cout << endl << "sgnl : " << sgnl; |
1559 |
|
cout << endl << "max.strip : "; |
1560 |
|
cout << endl << "multiplicity : "; |
1561 |
|
} |
1562 |
|
//-------------------------------------- |
1563 |
|
// |
1564 |
|
// |
1565 |
|
//-------------------------------------- |
1566 |
|
void TrkSinglet::Clear(){ |
1567 |
|
// cout << "TrkSinglet::Clear() " << GetUniqueID()<<endl; |
1568 |
|
// cls=0; |
1569 |
|
plane=-1; |
1570 |
|
coord[0]=-999; |
1571 |
|
coord[1]=-999; |
1572 |
|
sgnl=0; |
1573 |
|
multmax = 0; |
1574 |
|
|
1575 |
} |
} |
1576 |
//-------------------------------------- |
//-------------------------------------- |
1577 |
// |
// |
1578 |
// |
// |
1579 |
//-------------------------------------- |
//-------------------------------------- |
1580 |
TrkLevel2::TrkLevel2(){ |
TrkLevel2::TrkLevel2(){ |
1581 |
good2 = -1; |
// cout <<"TrkLevel2::TrkLevel2()"<<endl; |
1582 |
for(Int_t i=0; i<12 ; i++){ |
for(Int_t i=0; i<12 ; i++){ |
1583 |
crc[i] = -1; |
good[i] = -1; |
1584 |
|
VKmask[i] = 0; |
1585 |
|
VKflag[i] = 0; |
1586 |
}; |
}; |
1587 |
Track = new TClonesArray("TrkTrack"); |
Track = 0; |
1588 |
SingletX = new TClonesArray("TrkSinglet"); |
SingletX = 0; |
1589 |
SingletY = new TClonesArray("TrkSinglet"); |
SingletY = 0; |
1590 |
// Track = 0; |
|
1591 |
// Singlet = 0; |
} |
1592 |
// SingletY = 0; |
//-------------------------------------- |
1593 |
|
// |
1594 |
|
// |
1595 |
|
//-------------------------------------- |
1596 |
|
void TrkLevel2::Set(){ |
1597 |
|
if(!Track)Track = new TClonesArray("TrkTrack"); |
1598 |
|
if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); |
1599 |
|
if(!SingletY)SingletY = new TClonesArray("TrkSinglet"); |
1600 |
|
} |
1601 |
|
//-------------------------------------- |
1602 |
|
// |
1603 |
|
// |
1604 |
|
//-------------------------------------- |
1605 |
|
void TrkLevel2::SetTrackArray(TClonesArray *track){ |
1606 |
|
if(track && strcmp(track->GetClass()->GetName(),"TrkTrack")==0){ |
1607 |
|
if(Track)Track->Clear("C"); |
1608 |
|
Track = track; |
1609 |
|
} |
1610 |
} |
} |
1611 |
//-------------------------------------- |
//-------------------------------------- |
1612 |
// |
// |
1613 |
// |
// |
1614 |
//-------------------------------------- |
//-------------------------------------- |
1615 |
void TrkLevel2::Dump(){ |
void TrkLevel2::Dump(){ |
1616 |
TClonesArray &t = *Track; |
|
1617 |
TClonesArray &sx = *SingletX; |
// |
|
TClonesArray &sy = *SingletY; |
|
|
|
|
1618 |
cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-"; |
cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-"; |
1619 |
cout << endl << "good2 : " << good2; |
cout << endl << "good : "; for(int i=0; i<12; i++) cout << hex <<" 0x"<< good[i]<<dec; |
1620 |
cout << endl << "crc : "; for(int i=0; i<12; i++) cout << crc[i]; |
cout << endl << "ntrk() : " << ntrk() ; |
1621 |
cout << endl << "ntrk() : " << this->ntrk() ; |
cout << endl << "nclsx() : " << nclsx(); |
1622 |
cout << endl << "nclsx() : " << this->nclsx(); |
cout << endl << "nclsy() : " << nclsy(); |
1623 |
cout << endl << "nclsy() : " << this->nclsy(); |
if(Track){ |
1624 |
for(int i=0; i<this->ntrk(); i++) ((TrkTrack *)t[i])->Dump(); |
TClonesArray &t = *Track; |
1625 |
for(int i=0; i<this->nclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); |
for(int i=0; i<ntrk(); i++) ((TrkTrack *)t[i])->Dump(); |
1626 |
for(int i=0; i<this->nclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); |
} |
1627 |
|
// if(SingletX){ |
1628 |
|
// TClonesArray &sx = *SingletX; |
1629 |
|
// for(int i=0; i<nclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); |
1630 |
|
// } |
1631 |
|
// if(SingletY){ |
1632 |
|
// TClonesArray &sy = *SingletY; |
1633 |
|
// for(int i=0; i<nclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); |
1634 |
|
// } |
1635 |
|
cout << endl; |
1636 |
} |
} |
1637 |
|
/** |
1638 |
|
* \brief Dump processing status |
1639 |
|
*/ |
1640 |
|
void TrkLevel2::StatusDump(int view){ |
1641 |
|
cout << "DSP n. "<<view+1<<" status: "<<hex<<good[view]<<endl; |
1642 |
|
}; |
1643 |
|
/** |
1644 |
|
* \brief Check event status |
1645 |
|
* |
1646 |
|
* Check the event status, according to a flag-mask given as input. |
1647 |
|
* Return true if the view passes the check. |
1648 |
|
* |
1649 |
|
* @param view View number (0-11) |
1650 |
|
* @param flagmask Mask of flags to check (eg. flagmask=0x111 no missing packet, |
1651 |
|
* no crc error, no software alarm) |
1652 |
|
* |
1653 |
|
* @see TrkLevel2 class definition to know how the status flag is defined |
1654 |
|
* |
1655 |
|
*/ |
1656 |
|
Bool_t TrkLevel2::StatusCheck(int view, int flagmask){ |
1657 |
|
|
1658 |
|
if( view<0 || view >= 12)return false; |
1659 |
|
return !(good[view]&flagmask); |
1660 |
|
|
1661 |
|
}; |
1662 |
|
|
1663 |
|
|
1664 |
//-------------------------------------- |
//-------------------------------------- |
1665 |
// |
// |
1666 |
// |
// |
1667 |
//-------------------------------------- |
//-------------------------------------- |
1668 |
/** |
/** |
1669 |
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). |
* The method returns false if the viking-chip was masked |
1670 |
|
* either apriori ,on the basis of the mask read from the DB, |
1671 |
|
* or run-by-run, on the basis of the calibration parameters) |
1672 |
|
* @param iv Tracker view (0-11) |
1673 |
|
* @param ivk Viking-chip number (0-23) |
1674 |
*/ |
*/ |
1675 |
void TrkLevel2::FillCommonVar(cTrkLevel2 *l2){ |
Bool_t TrkLevel2::GetVKMask(int iv, int ivk){ |
1676 |
// |
Int_t whichbit = (Int_t)pow(2,ivk); |
1677 |
// Track = new TClonesArray("TrkTrack"); |
return (whichbit&VKmask[iv])!=0; |
1678 |
// SingletX = new TClonesArray("TrkSinglet"); |
} |
1679 |
// SingletY = new TClonesArray("TrkSinglet"); |
/** |
1680 |
|
* The method returns false if the viking-chip was masked |
1681 |
|
* for this event due to common-noise computation failure. |
1682 |
|
* @param iv Tracker view (0-11) |
1683 |
|
* @param ivk Viking-chip number (0-23) |
1684 |
|
*/ |
1685 |
|
Bool_t TrkLevel2::GetVKFlag(int iv, int ivk){ |
1686 |
|
Int_t whichbit = (Int_t)pow(2,ivk); |
1687 |
|
return (whichbit&VKflag[iv])!=0; |
1688 |
|
} |
1689 |
|
/** |
1690 |
|
* The method returns true if the viking-chip was masked, either |
1691 |
|
* forced (see TrkLevel2::GetVKMask(int,int)) or |
1692 |
|
* for this event only (TrkLevel2::GetVKFlag(int,int)). |
1693 |
|
* @param iv Tracker view (0-11) |
1694 |
|
* @param ivk Viking-chip number (0-23) |
1695 |
|
*/ |
1696 |
|
Bool_t TrkLevel2::IsMaskedVK(int iv, int ivk){ |
1697 |
|
return !(GetVKMask(iv,ivk)&&GetVKFlag(iv,ivk) ); |
1698 |
|
}; |
1699 |
|
|
1700 |
|
//-------------------------------------- |
1701 |
|
// |
1702 |
|
// |
1703 |
|
//-------------------------------------- |
1704 |
|
/** |
1705 |
|
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). |
1706 |
|
* Ref to Level1 data (clusters) is also set. If l1==NULL no references are set. |
1707 |
|
* (NB It make sense to set references only if events are stored in a tree that contains also the Level1 branch) |
1708 |
|
*/ |
1709 |
|
void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){ |
1710 |
|
|
1711 |
|
// cout << "void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1)"<<endl; |
1712 |
|
Clear(); |
1713 |
|
|
1714 |
// temporary objects: |
// temporary objects: |
1715 |
TrkSinglet* t_singlet = new TrkSinglet(); |
TrkSinglet* t_singlet = new TrkSinglet(); |
1716 |
TrkTrack* t_track = new TrkTrack(); |
TrkTrack* t_track = new TrkTrack(); |
1717 |
|
|
1718 |
|
// ----------------- |
1719 |
// general variables |
// general variables |
1720 |
good2 = l2->good2; |
// ----------------- |
1721 |
for(Int_t i=0; i<12 ; i++){ |
for(Int_t i=0; i<12 ; i++){ |
1722 |
crc[i] = l2->crc[i]; |
good[i] = l2->good[i]; |
1723 |
|
VKmask[i]=0; |
1724 |
|
VKflag[i]=0; |
1725 |
|
for(Int_t ii=0; ii<24 ; ii++){ |
1726 |
|
Int_t setbit = (Int_t)pow(2,ii); |
1727 |
|
if( l2->vkflag[ii][i]!=-1 )VKmask[i]=VKmask[i]|setbit; |
1728 |
|
if( l2->vkflag[ii][i]!=0 )VKflag[i]=VKflag[i]|setbit; |
1729 |
|
}; |
1730 |
}; |
}; |
1731 |
|
// -------------- |
1732 |
// *** TRACKS *** |
// *** TRACKS *** |
1733 |
|
// -------------- |
1734 |
|
if(!Track) Track = new TClonesArray("TrkTrack"); |
1735 |
TClonesArray &t = *Track; |
TClonesArray &t = *Track; |
1736 |
|
|
1737 |
for(int i=0; i<l2->ntrk; i++){ |
for(int i=0; i<l2->ntrk; i++){ |
1738 |
t_track->seqno = i; |
t_track->seqno = i;// NBNBNBNB deve sempre essere = i |
1739 |
t_track->image = l2->image[i]-1; |
t_track->image = l2->image[i]-1; |
|
// cout << "track "<<i<<t_track->seqno << t_track->image<<endl; |
|
1740 |
t_track->chi2 = l2->chi2_nt[i]; |
t_track->chi2 = l2->chi2_nt[i]; |
1741 |
|
t_track->nstep = l2->nstep_nt[i]; |
1742 |
for(int it1=0;it1<5;it1++){ |
for(int it1=0;it1<5;it1++){ |
1743 |
t_track->al[it1] = l2->al_nt[i][it1]; |
t_track->al[it1] = l2->al_nt[i][it1]; |
1744 |
for(int it2=0;it2<5;it2++) |
for(int it2=0;it2<5;it2++) |
1745 |
t_track->coval[it1][it2] = l2->coval[i][it2][it1]; |
t_track->coval[it1][it2] = l2->coval[i][it2][it1]; |
1746 |
}; |
}; |
1747 |
for(int ip=0;ip<6;ip++){ |
for(int ip=0;ip<6;ip++){ |
1748 |
t_track->xgood[ip] = l2->xgood_nt[i][ip]; |
// --------------------------------- |
1749 |
t_track->ygood[ip] = l2->ygood_nt[i][ip]; |
// new implementation of xgood/ygood |
1750 |
|
// --------------------------------- |
1751 |
|
t_track->xgood[ip] = l2->cltrx[i][ip]; //cluster ID |
1752 |
|
t_track->ygood[ip] = l2->cltry[i][ip]; //cluster ID |
1753 |
|
t_track->xgood[ip] += 10000000*l2->ls[i][ip]; // ladder+sensor |
1754 |
|
t_track->ygood[ip] += 10000000*l2->ls[i][ip]; // ladder+sensor |
1755 |
|
if(l2->xbad[i][ip]>0)t_track->xgood[ip]=-t_track->xgood[ip]; |
1756 |
|
if(l2->ybad[i][ip]>0)t_track->ygood[ip]=-t_track->ygood[ip]; |
1757 |
|
// if(l2->xbad[i][ip]>0 || l2->ybad[i][ip]>0){ |
1758 |
|
// if(l2->dedx_x[i][ip]<0 || l2->dedx_y[i][ip]<0){ |
1759 |
|
// cout << ip << " - "<< l2->cltrx[i][ip] << " "<<l2->cltry[i][ip]<<" "<<l2->ls[i][ip]<<endl; |
1760 |
|
// cout << ip << " - "<<t_track->xgood[ip]<<" "<<t_track->ygood[ip]<<endl; |
1761 |
|
// cout << ip << " - "<<t_track->GetClusterX_ID(ip)<<" "<<t_track->GetClusterY_ID(ip)<<" "<<t_track->GetLadder(ip)<<" "<<t_track->GetSensor(ip)<<endl; |
1762 |
|
// cout << ip << " - "<<t_track->BadClusterX(ip)<<" "<<t_track->BadClusterY(ip)<<endl; |
1763 |
|
// cout << ip << " - "<<t_track->SaturatedClusterX(ip)<<" "<<t_track->SaturatedClusterY(ip)<<endl; |
1764 |
|
// } |
1765 |
t_track->xm[ip] = l2->xm_nt[i][ip]; |
t_track->xm[ip] = l2->xm_nt[i][ip]; |
1766 |
t_track->ym[ip] = l2->ym_nt[i][ip]; |
t_track->ym[ip] = l2->ym_nt[i][ip]; |
1767 |
t_track->zm[ip] = l2->zm_nt[i][ip]; |
t_track->zm[ip] = l2->zm_nt[i][ip]; |
1768 |
t_track->resx[ip] = l2->resx_nt[i][ip]; |
t_track->resx[ip] = l2->resx_nt[i][ip]; |
1769 |
t_track->resy[ip] = l2->resy_nt[i][ip]; |
t_track->resy[ip] = l2->resy_nt[i][ip]; |
1770 |
|
t_track->tailx[ip] = l2->tailx[i][ip]; |
1771 |
|
t_track->taily[ip] = l2->taily[i][ip]; |
1772 |
t_track->xv[ip] = l2->xv_nt[i][ip]; |
t_track->xv[ip] = l2->xv_nt[i][ip]; |
1773 |
t_track->yv[ip] = l2->yv_nt[i][ip]; |
t_track->yv[ip] = l2->yv_nt[i][ip]; |
1774 |
t_track->zv[ip] = l2->zv_nt[i][ip]; |
t_track->zv[ip] = l2->zv_nt[i][ip]; |
1776 |
t_track->ayv[ip] = l2->ayv_nt[i][ip]; |
t_track->ayv[ip] = l2->ayv_nt[i][ip]; |
1777 |
t_track->dedx_x[ip] = l2->dedx_x[i][ip]; |
t_track->dedx_x[ip] = l2->dedx_x[i][ip]; |
1778 |
t_track->dedx_y[ip] = l2->dedx_y[i][ip]; |
t_track->dedx_y[ip] = l2->dedx_y[i][ip]; |
1779 |
|
t_track->multmaxx[ip] = l2->multmaxx[i][ip]; |
1780 |
|
t_track->multmaxy[ip] = l2->multmaxy[i][ip]; |
1781 |
|
t_track->seedx[ip] = l2->seedx[i][ip]; |
1782 |
|
t_track->seedy[ip] = l2->seedy[i][ip]; |
1783 |
|
t_track->xpu[ip] = l2->xpu[i][ip]; |
1784 |
|
t_track->ypu[ip] = l2->ypu[i][ip]; |
1785 |
|
//----------------------------------------------------- |
1786 |
|
//----------------------------------------------------- |
1787 |
|
//----------------------------------------------------- |
1788 |
|
//----------------------------------------------------- |
1789 |
}; |
}; |
1790 |
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
1791 |
|
// evaluated coordinates (to define GF) |
1792 |
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
1793 |
|
int ngf = TrkParams::nGF; |
1794 |
|
float *zgf = TrkParams::zGF; |
1795 |
|
Trajectory tgf = Trajectory(ngf,zgf); |
1796 |
|
tgf.DoTrack(t_track->al);//<<<< integrate the trajectory |
1797 |
|
for(int ip=0; ip<ngf; ip++){ |
1798 |
|
t_track->xGF[ip] = tgf.x[ip]; |
1799 |
|
t_track->yGF[ip] = tgf.y[ip]; |
1800 |
|
} |
1801 |
|
|
1802 |
|
// if(t_track->IsSaturated())t_track->Dump(); |
1803 |
new(t[i]) TrkTrack(*t_track); |
new(t[i]) TrkTrack(*t_track); |
1804 |
t_track->Clear(); |
t_track->Clear(); |
1805 |
}; |
};//end loop over track |
1806 |
|
|
1807 |
|
// ---------------- |
1808 |
// *** SINGLETS *** |
// *** SINGLETS *** |
1809 |
|
// ---------------- |
1810 |
|
if(!SingletX)SingletX = new TClonesArray("TrkSinglet"); |
1811 |
TClonesArray &sx = *SingletX; |
TClonesArray &sx = *SingletX; |
1812 |
for(int i=0; i<l2->nclsx; i++){ |
for(int i=0; i<l2->nclsx; i++){ |
1813 |
t_singlet->plane = l2->planex[i]; |
t_singlet->plane = l2->planex[i]; |
1814 |
t_singlet->coord[0] = l2->xs[i][0]; |
t_singlet->coord[0] = l2->xs[i][0]; |
1815 |
t_singlet->coord[1] = l2->xs[i][1]; |
t_singlet->coord[1] = l2->xs[i][1]; |
1816 |
t_singlet->sgnl = l2->signlxs[i]; |
t_singlet->sgnl = l2->signlxs[i]; |
1817 |
|
t_singlet->multmax = l2->multmaxsx[i]; |
1818 |
|
if(l2->sxbad[i]>0) t_singlet->multmax = -1*t_singlet->multmax; |
1819 |
|
//----------------------------------------------------- |
1820 |
|
// if(l1) t_singlet->cls = l1->GetCluster(l2->clsx[i]-1); |
1821 |
|
//----------------------------------------------------- |
1822 |
new(sx[i]) TrkSinglet(*t_singlet); |
new(sx[i]) TrkSinglet(*t_singlet); |
1823 |
t_singlet->Clear(); |
t_singlet->Clear(); |
1824 |
} |
} |
1825 |
|
if(!SingletY)SingletY = new TClonesArray("TrkSinglet"); |
1826 |
TClonesArray &sy = *SingletY; |
TClonesArray &sy = *SingletY; |
1827 |
for(int i=0; i<l2->nclsy; i++){ |
for(int i=0; i<l2->nclsy; i++){ |
1828 |
t_singlet->plane = l2->planey[i]; |
t_singlet->plane = l2->planey[i]; |
1829 |
t_singlet->coord[0] = l2->ys[i][0]; |
t_singlet->coord[0] = l2->ys[i][0]; |
1830 |
t_singlet->coord[1] = l2->ys[i][1]; |
t_singlet->coord[1] = l2->ys[i][1]; |
1831 |
t_singlet->sgnl = l2->signlys[i]; |
t_singlet->sgnl = l2->signlys[i]; |
1832 |
|
t_singlet->multmax = l2->multmaxsy[i]; |
1833 |
|
if(l2->sybad[i]>0) t_singlet->multmax = -1*t_singlet->multmax; |
1834 |
|
//----------------------------------------------------- |
1835 |
|
// if(l1) t_singlet->cls = l1->GetCluster(l2->clsy[i]-1); |
1836 |
|
//----------------------------------------------------- |
1837 |
new(sy[i]) TrkSinglet(*t_singlet); |
new(sy[i]) TrkSinglet(*t_singlet); |
1838 |
t_singlet->Clear(); |
t_singlet->Clear(); |
1839 |
}; |
}; |
1840 |
|
|
1841 |
|
|
1842 |
|
|
1843 |
|
delete t_track; |
1844 |
|
delete t_singlet; |
1845 |
|
} |
1846 |
|
/** |
1847 |
|
* Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common). |
1848 |
|
*/ |
1849 |
|
|
1850 |
|
void TrkLevel2::GetLevel2Struct(cTrkLevel2 *l2) const { |
1851 |
|
|
1852 |
|
// general variables |
1853 |
|
// l2->good2 = good2 ; |
1854 |
|
for(Int_t i=0; i<12 ; i++){ |
1855 |
|
// l2->crc[i] = crc[i]; |
1856 |
|
l2->good[i] = good[i]; |
1857 |
|
}; |
1858 |
|
// *** TRACKS *** |
1859 |
|
|
1860 |
|
if(Track){ |
1861 |
|
l2->ntrk = Track->GetEntries(); |
1862 |
|
for(Int_t i=0;i<l2->ntrk;i++){ |
1863 |
|
l2->image[i] = 1 + ((TrkTrack *)Track->At(i))->image; |
1864 |
|
l2->chi2_nt[i] = ((TrkTrack *)Track->At(i))->chi2; |
1865 |
|
l2->nstep_nt[i] = ((TrkTrack *)Track->At(i))->nstep; |
1866 |
|
for(int it1=0;it1<5;it1++){ |
1867 |
|
l2->al_nt[i][it1] = ((TrkTrack *)Track->At(i))->al[it1]; |
1868 |
|
for(int it2=0;it2<5;it2++) |
1869 |
|
l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2]; |
1870 |
|
}; |
1871 |
|
for(int ip=0;ip<6;ip++){ |
1872 |
|
l2->xgood_nt[i][ip] = ((TrkTrack *)Track->At(i))->XGood(ip); |
1873 |
|
l2->ygood_nt[i][ip] = ((TrkTrack *)Track->At(i))->YGood(ip); |
1874 |
|
l2->xm_nt[i][ip] = ((TrkTrack *)Track->At(i))->xm[ip]; |
1875 |
|
l2->ym_nt[i][ip] = ((TrkTrack *)Track->At(i))->ym[ip]; |
1876 |
|
l2->zm_nt[i][ip] = ((TrkTrack *)Track->At(i))->zm[ip]; |
1877 |
|
l2->resx_nt[i][ip] = ((TrkTrack *)Track->At(i))->resx[ip]; |
1878 |
|
l2->resy_nt[i][ip] = ((TrkTrack *)Track->At(i))->resy[ip]; |
1879 |
|
l2->tailx[i][ip] = ((TrkTrack *)Track->At(i))->tailx[ip]; |
1880 |
|
l2->taily[i][ip] = ((TrkTrack *)Track->At(i))->taily[ip]; |
1881 |
|
l2->xv_nt[i][ip] = ((TrkTrack *)Track->At(i))->xv[ip]; |
1882 |
|
l2->yv_nt[i][ip] = ((TrkTrack *)Track->At(i))->yv[ip]; |
1883 |
|
l2->zv_nt[i][ip] = ((TrkTrack *)Track->At(i))->zv[ip]; |
1884 |
|
l2->axv_nt[i][ip] = ((TrkTrack *)Track->At(i))->axv[ip]; |
1885 |
|
l2->ayv_nt[i][ip] = ((TrkTrack *)Track->At(i))->ayv[ip]; |
1886 |
|
l2->dedx_x[i][ip] = ((TrkTrack *)Track->At(i))->dedx_x[ip]; |
1887 |
|
l2->dedx_y[i][ip] = ((TrkTrack *)Track->At(i))->dedx_y[ip]; |
1888 |
|
}; |
1889 |
|
} |
1890 |
|
} |
1891 |
|
// *** SINGLETS *** |
1892 |
|
if(SingletX){ |
1893 |
|
l2->nclsx = SingletX->GetEntries(); |
1894 |
|
for(Int_t i=0;i<l2->nclsx;i++){ |
1895 |
|
l2->planex[i] = ((TrkSinglet *)SingletX->At(i))->plane; |
1896 |
|
l2->xs[i][0] = ((TrkSinglet *)SingletX->At(i))->coord[0]; |
1897 |
|
l2->xs[i][1] = ((TrkSinglet *)SingletX->At(i))->coord[1]; |
1898 |
|
l2->signlxs[i] = ((TrkSinglet *)SingletX->At(i))->sgnl; |
1899 |
|
} |
1900 |
|
} |
1901 |
|
|
1902 |
|
if(SingletY){ |
1903 |
|
l2->nclsy = SingletY->GetEntries(); |
1904 |
|
for(Int_t i=0;i<l2->nclsy;i++){ |
1905 |
|
l2->planey[i] = ((TrkSinglet *)SingletY->At(i))->plane; |
1906 |
|
l2->ys[i][0] = ((TrkSinglet *)SingletY->At(i))->coord[0]; |
1907 |
|
l2->ys[i][1] = ((TrkSinglet *)SingletY->At(i))->coord[1]; |
1908 |
|
l2->signlys[i] = ((TrkSinglet *)SingletY->At(i))->sgnl; |
1909 |
|
} |
1910 |
|
} |
1911 |
} |
} |
1912 |
//-------------------------------------- |
//-------------------------------------- |
1913 |
// |
// |
1914 |
// |
// |
1915 |
//-------------------------------------- |
//-------------------------------------- |
1916 |
void TrkLevel2::Clear(){ |
void TrkLevel2::Clear(){ |
|
good2 = -1; |
|
1917 |
for(Int_t i=0; i<12 ; i++){ |
for(Int_t i=0; i<12 ; i++){ |
1918 |
crc[i] = -1; |
good[i] = -1; |
1919 |
|
VKflag[i] = 0; |
1920 |
|
VKmask[i] = 0; |
1921 |
}; |
}; |
1922 |
Track->RemoveAll(); |
// if(Track)Track->Clear("C"); |
1923 |
SingletX->RemoveAll(); |
// if(SingletX)SingletX->Clear("C"); |
1924 |
SingletY->RemoveAll(); |
// if(SingletY)SingletY->Clear("C"); |
1925 |
|
if(Track)Track->Delete(); |
1926 |
|
if(SingletX)SingletX->Delete(); |
1927 |
|
if(SingletY)SingletY->Delete(); |
1928 |
|
} |
1929 |
|
// //-------------------------------------- |
1930 |
|
// // |
1931 |
|
// // |
1932 |
|
// //-------------------------------------- |
1933 |
|
void TrkLevel2::Delete(){ |
1934 |
|
|
1935 |
|
// cout << "void TrkLevel2::Delete()"<<endl; |
1936 |
|
Clear(); |
1937 |
|
if(Track)delete Track; |
1938 |
|
if(SingletX)delete SingletX; |
1939 |
|
if(SingletY)delete SingletY; |
1940 |
|
|
1941 |
} |
} |
1942 |
//-------------------------------------- |
//-------------------------------------- |
1943 |
// |
// |
1947 |
* 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). |
* 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). |
1948 |
* This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used. |
* This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used. |
1949 |
*/ |
*/ |
1950 |
TClonesArray *TrkLevel2::GetTracks(){ |
TRefArray *TrkLevel2::GetTracks_NFitSorted(){ |
|
TClonesArray *sorted = GetTracks_NFitSorted(); |
|
|
return sorted; |
|
|
}; |
|
|
TClonesArray *TrkLevel2::GetTracks_Chi2Sorted(){ |
|
1951 |
|
|
1952 |
TClonesArray *sorted = new TClonesArray("TrkTrack"); |
if(!Track)return 0; |
|
TClonesArray &t = *Track; |
|
|
TClonesArray &ts = *sorted; |
|
|
int N=this->ntrk(); |
|
|
vector<int> m(N); for(int i=0; i<N; i++)m[i]=1; |
|
|
|
|
|
int indo=0; |
|
|
int indi=0; |
|
|
while(N != 0){ |
|
|
float chi2ref=1000000; |
|
|
for(int i=0; i<this->ntrk(); i++){ |
|
|
if(((TrkTrack *)t[i])->chi2 < chi2ref && m[i]==1){ |
|
|
chi2ref = ((TrkTrack *)t[i])->chi2; |
|
|
indi = i; |
|
|
} |
|
|
} |
|
|
if( ((TrkTrack *)t[indi])->image != -1 ){ |
|
|
m[((TrkTrack *)t[indi])->image] = 0; |
|
|
N--; |
|
|
} |
|
|
new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]); |
|
|
m[indi] = 0; |
|
|
N--; |
|
|
indo++; |
|
|
} |
|
|
return sorted; |
|
|
} |
|
|
TClonesArray *TrkLevel2::GetTracks_NFitSorted(){ |
|
1953 |
|
|
1954 |
TClonesArray *sorted = new TClonesArray("TrkTrack"); |
// TRefArray *sorted = new TRefArray(); |
1955 |
TClonesArray &t = *Track; |
TRefArray *sorted = NULL; |
1956 |
TClonesArray &ts = *sorted; |
|
1957 |
int N=this->ntrk(); |
TClonesArray &t = *Track; |
1958 |
|
// TClonesArray &ts = *PhysicalTrack; |
1959 |
|
int N = ntrk(); |
1960 |
vector<int> m(N); for(int i=0; i<N; i++)m[i]=1; |
vector<int> m(N); for(int i=0; i<N; i++)m[i]=1; |
1961 |
|
// int m[50]; for(int i=0; i<N; i++)m[i]=1; |
1962 |
|
|
1963 |
int indo=0; |
int indo=0; |
1964 |
int indi=0; |
int indi=0; |
1965 |
while(N != 0){ |
while(N > 0){ |
1966 |
|
// while(N != 0){ |
1967 |
int nfit =0; |
int nfit =0; |
1968 |
float chi2ref=1000000; |
float chi2ref = numeric_limits<float>::max(); |
1969 |
|
|
1970 |
// first loop to search maximum num. of fit points |
// first loop to search maximum num. of fit points |
1971 |
for(int i=0; i<this->ntrk(); i++){ |
for(int i=0; i < ntrk(); i++){ |
1972 |
if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){ |
if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){ |
1973 |
nfit = ((TrkTrack *)t[i])->GetNtot(); |
nfit = ((TrkTrack *)t[i])->GetNtot(); |
|
// cout << "1** "<<i<< " " << nfit<<endl; |
|
1974 |
} |
} |
1975 |
} |
} |
1976 |
//second loop to search minimum chi2 among selected |
//second loop to search minimum chi2 among selected |
1977 |
for(int i=0; i<this->ntrk(); i++){ |
for(int i=0; i<ntrk(); i++){ |
1978 |
if( ((TrkTrack *)t[i])->chi2 < chi2ref |
Float_t chi2 = ((TrkTrack *)t[i])->chi2; |
1979 |
&& ((TrkTrack *)t[i])->GetNtot()== nfit |
if(chi2 < 0) chi2 = -chi2*1000; |
1980 |
&& m[i]==1){ |
if( chi2 < chi2ref |
1981 |
|
&& ((TrkTrack *)t[i])->GetNtot() == nfit |
1982 |
|
&& m[i]==1){ |
1983 |
chi2ref = ((TrkTrack *)t[i])->chi2; |
chi2ref = ((TrkTrack *)t[i])->chi2; |
1984 |
indi = i; |
indi = i; |
1985 |
// cout << "2** "<<i<< " " << nfit <<" "<<chi2ref<<endl; |
}; |
1986 |
} |
}; |
|
} |
|
1987 |
if( ((TrkTrack *)t[indi])->HasImage() ){ |
if( ((TrkTrack *)t[indi])->HasImage() ){ |
1988 |
m[((TrkTrack *)t[indi])->image] = 0; |
m[((TrkTrack *)t[indi])->image] = 0; |
1989 |
N--; |
N--; |
1990 |
|
|
|
// Int_t nfiti=((TrkTrack *)t[((TrkTrack *)t[indi])->image ])->GetNtot(); |
|
|
// Float_t chi2i=((TrkTrack *)t[((TrkTrack *)t[indi])->image ])->chi2; |
|
|
|
|
1991 |
// cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl; |
// cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl; |
1992 |
} |
}; |
1993 |
new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]); |
if(!sorted)sorted = new TRefArray( TProcessID::GetProcessWithUID(t[indi])); |
1994 |
|
sorted->Add( (TrkTrack*)t[indi] ); |
1995 |
|
|
1996 |
m[indi] = 0; |
m[indi] = 0; |
1997 |
|
// cout << "SORTED "<< indo << " "<< indi << " "<< N << " "<<((TrkTrack *)t[indi])->image<<" "<<chi2ref<<endl; |
1998 |
N--; |
N--; |
1999 |
indo++; |
indo++; |
2000 |
} |
} |
2001 |
|
m.clear(); |
2002 |
|
// cout << "GetTracks_NFitSorted(it): Done"<< endl; |
2003 |
|
|
2004 |
return sorted; |
return sorted; |
2005 |
|
// return PhysicalTrack; |
2006 |
} |
} |
2007 |
//-------------------------------------- |
//-------------------------------------- |
2008 |
// |
// |
2016 |
TrkTrack *TrkLevel2::GetStoredTrack(int is){ |
TrkTrack *TrkLevel2::GetStoredTrack(int is){ |
2017 |
|
|
2018 |
if(is >= this->ntrk()){ |
if(is >= this->ntrk()){ |
2019 |
cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl; |
cout << "TrkTrack *TrkLevel2::GetStoredTrack(int) >> Track "<< is << "doen not exits! " << endl; |
2020 |
cout << " Stored tracks ntrk() = "<< this->ntrk() << endl; |
cout << "Stored tracks ntrk() = "<< this->ntrk() << endl; |
2021 |
return 0; |
return 0; |
2022 |
} |
} |
2023 |
|
if(!Track){ |
2024 |
|
cout << "TrkTrack *TrkLevel2::GetStoredTrack(int is) >> (TClonesArray*) Track ==0 "<<endl; |
2025 |
|
}; |
2026 |
TClonesArray &t = *(Track); |
TClonesArray &t = *(Track); |
2027 |
TrkTrack *track = (TrkTrack*)t[is]; |
TrkTrack *track = (TrkTrack*)t[is]; |
2028 |
return track; |
return track; |
2032 |
// |
// |
2033 |
//-------------------------------------- |
//-------------------------------------- |
2034 |
/** |
/** |
2035 |
|
* Retrieves the is-th stored X singlet. |
2036 |
|
* @param it Singlet number, ranging from 0 to nclsx(). |
2037 |
|
*/ |
2038 |
|
TrkSinglet *TrkLevel2::GetSingletX(int is){ |
2039 |
|
|
2040 |
|
if(is >= this->nclsx()){ |
2041 |
|
cout << "TrkSinglet *TrkLevel2::GetSingletX(int) >> Singlet "<< is << "doen not exits! " << endl; |
2042 |
|
cout << "Stored x-singlets nclsx() = "<< this->nclsx() << endl; |
2043 |
|
return 0; |
2044 |
|
} |
2045 |
|
if(!SingletX)return 0; |
2046 |
|
TClonesArray &t = *(SingletX); |
2047 |
|
TrkSinglet *singlet = (TrkSinglet*)t[is]; |
2048 |
|
return singlet; |
2049 |
|
} |
2050 |
|
//-------------------------------------- |
2051 |
|
// |
2052 |
|
// |
2053 |
|
//-------------------------------------- |
2054 |
|
/** |
2055 |
|
* Retrieves the is-th stored Y singlet. |
2056 |
|
* @param it Singlet number, ranging from 0 to nclsx(). |
2057 |
|
*/ |
2058 |
|
TrkSinglet *TrkLevel2::GetSingletY(int is){ |
2059 |
|
|
2060 |
|
if(is >= this->nclsy()){ |
2061 |
|
cout << "TrkSinglet *TrkLevel2::GetSingletY(int) >> Singlet "<< is << "doen not exits! " << endl; |
2062 |
|
cout << "Stored y-singlets nclsx() = "<< this->nclsx() << endl; |
2063 |
|
return 0; |
2064 |
|
} |
2065 |
|
if(!SingletY)return 0; |
2066 |
|
TClonesArray &t = *(SingletY); |
2067 |
|
TrkSinglet *singlet = (TrkSinglet*)t[is]; |
2068 |
|
return singlet; |
2069 |
|
} |
2070 |
|
//-------------------------------------- |
2071 |
|
// |
2072 |
|
// |
2073 |
|
//-------------------------------------- |
2074 |
|
/** |
2075 |
* Retrieves the it-th "physical" track, sorted by the method GetNTracks(). |
* Retrieves the it-th "physical" track, sorted by the method GetNTracks(). |
2076 |
* @param it Track number, ranging from 0 to GetNTracks(). |
* @param it Track number, ranging from 0 to GetNTracks(). |
2077 |
*/ |
*/ |
2078 |
|
|
2079 |
TrkTrack *TrkLevel2::GetTrack(int it){ |
TrkTrack *TrkLevel2::GetTrack(int it){ |
2080 |
|
|
2081 |
if(it >= this->GetNTracks()){ |
if(it >= this->GetNTracks()){ |
2082 |
cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; |
cout << "TrkTrack *TrkLevel2::GetTrack(int) >> Track "<< it << "does not exits! " << endl; |
2083 |
cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl; |
cout << "Physical tracks GetNTracks() = "<< this->ntrk() << endl; |
2084 |
return 0; |
return 0; |
2085 |
} |
} |
2086 |
TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it]; |
|
2087 |
return track; |
TRefArray *sorted = GetTracks(); //TEMPORANEO |
2088 |
|
if(!sorted)return 0; |
2089 |
|
TrkTrack *track = (TrkTrack*)sorted->At(it); |
2090 |
|
sorted->Clear(); |
2091 |
|
delete sorted; |
2092 |
|
return track; |
2093 |
} |
} |
2094 |
|
/** |
2095 |
|
* Give the number of "physical" tracks, sorted by the method GetTracks(). |
2096 |
|
*/ |
2097 |
|
Int_t TrkLevel2::GetNTracks(){ |
2098 |
|
|
2099 |
|
Float_t ntot=0; |
2100 |
|
if(!Track)return 0; |
2101 |
|
TClonesArray &t = *Track; |
2102 |
|
for(int i=0; i<ntrk(); i++) { |
2103 |
|
if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.; |
2104 |
|
else ntot+=0.5; |
2105 |
|
} |
2106 |
|
return (Int_t)ntot; |
2107 |
|
|
2108 |
|
}; |
2109 |
//-------------------------------------- |
//-------------------------------------- |
2110 |
// |
// |
2111 |
// |
// |
2117 |
TrkTrack *TrkLevel2::GetTrackImage(int it){ |
TrkTrack *TrkLevel2::GetTrackImage(int it){ |
2118 |
|
|
2119 |
if(it >= this->GetNTracks()){ |
if(it >= this->GetNTracks()){ |
2120 |
cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; |
cout << "TrkTrack *TrkLevel2::GetTrackImage(int) >> Track "<< it << "does not exits! " << endl; |
2121 |
cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl; |
cout << "Physical tracks GetNTracks() = "<< this->ntrk() << endl; |
2122 |
return 0; |
return 0; |
2123 |
} |
} |
2124 |
TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it]; |
|
2125 |
|
TRefArray* sorted = GetTracks(); //TEMPORANEO |
2126 |
|
if(!sorted)return 0; |
2127 |
|
TrkTrack *track = (TrkTrack*)sorted->At(it); |
2128 |
|
|
2129 |
if(!track->HasImage()){ |
if(!track->HasImage()){ |
2130 |
cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl; |
cout << "TrkTrack *TrkLevel2::GetTrackImage(int) >> Track "<< it << "does not have image! " << endl; |
2131 |
return 0; |
return 0; |
2132 |
} |
} |
2133 |
|
if(!Track)return 0; |
2134 |
TrkTrack *image = (TrkTrack*)(*Track)[track->image]; |
TrkTrack *image = (TrkTrack*)(*Track)[track->image]; |
2135 |
|
|
2136 |
|
sorted->Delete(); |
2137 |
|
delete sorted; |
2138 |
|
|
2139 |
return image; |
return image; |
2140 |
|
|
2141 |
} |
} |
2147 |
* Loads the magnetic field. |
* Loads the magnetic field. |
2148 |
* @param s Path of the magnetic-field files. |
* @param s Path of the magnetic-field files. |
2149 |
*/ |
*/ |
2150 |
void TrkLevel2::LoadField(TString s){ |
void TrkLevel2::LoadField(TString path){ |
2151 |
readb_(s.Data()); |
// |
2152 |
|
// strcpy(path_.path,path.Data()); |
2153 |
|
// path_.pathlen = path.Length(); |
2154 |
|
// path_.error = 0; |
2155 |
|
// readb_(); |
2156 |
|
|
2157 |
|
// TrkParams::SetTrackingMode(); |
2158 |
|
// TrkParams::SetPrecisionFactor(); |
2159 |
|
// TrkParams::SetStepMin(); |
2160 |
|
TrkParams::SetMiniDefault(); |
2161 |
|
|
2162 |
|
TrkParams::Set(path,1); |
2163 |
|
TrkParams::Load(1); |
2164 |
|
if( !TrkParams::IsLoaded(1) ){ |
2165 |
|
cout << "void TrkLevel2::LoadField(TString path) --- ERROR --- m.field not loaded"<<endl; |
2166 |
|
} |
2167 |
|
|
2168 |
|
// |
2169 |
|
}; |
2170 |
|
// /** |
2171 |
|
// * Get BY (kGauss) |
2172 |
|
// * @param v (x,y,z) coordinates in cm |
2173 |
|
// */ |
2174 |
|
// float TrkLevel2::GetBX(float* v){ |
2175 |
|
// float b[3]; |
2176 |
|
// gufld_(v,b); |
2177 |
|
// return b[0]/10.; |
2178 |
|
// } |
2179 |
|
// /** |
2180 |
|
// * Get BY (kGauss) |
2181 |
|
// * @param v (x,y,z) coordinates in cm |
2182 |
|
// */ |
2183 |
|
// float TrkLevel2::GetBY(float* v){ |
2184 |
|
// float b[3]; |
2185 |
|
// gufld_(v,b); |
2186 |
|
// return b[1]/10.; |
2187 |
|
// } |
2188 |
|
// /** |
2189 |
|
// * Get BY (kGauss) |
2190 |
|
// * @param v (x,y,z) coordinates in cm |
2191 |
|
// */ |
2192 |
|
// float TrkLevel2::GetBZ(float* v){ |
2193 |
|
// float b[3]; |
2194 |
|
// gufld_(v,b); |
2195 |
|
// return b[2]/10.; |
2196 |
|
// } |
2197 |
|
//-------------------------------------- |
2198 |
|
// |
2199 |
|
// |
2200 |
|
//-------------------------------------- |
2201 |
|
/** |
2202 |
|
* Get tracker-plane (mechanical) z-coordinate |
2203 |
|
* @param plane_id plane index (1=TOP,2,3,4,5,6=BOTTOM) |
2204 |
|
*/ |
2205 |
|
Float_t TrkLevel2::GetZTrk(Int_t plane_id){ |
2206 |
|
switch(plane_id){ |
2207 |
|
case 1: return ZTRK1; |
2208 |
|
case 2: return ZTRK2; |
2209 |
|
case 3: return ZTRK3; |
2210 |
|
case 4: return ZTRK4; |
2211 |
|
case 5: return ZTRK5; |
2212 |
|
case 6: return ZTRK6; |
2213 |
|
default: return 0.; |
2214 |
|
}; |
2215 |
}; |
}; |
2216 |
//-------------------------------------- |
//-------------------------------------- |
2217 |
// |
// |
2222 |
* (By default is created with z-coordinates inside the tracking volume) |
* (By default is created with z-coordinates inside the tracking volume) |
2223 |
*/ |
*/ |
2224 |
Trajectory::Trajectory(){ |
Trajectory::Trajectory(){ |
2225 |
npoint = 10; |
npoint = 6; |
2226 |
x = new float[npoint]; |
x = new float[npoint]; |
2227 |
y = new float[npoint]; |
y = new float[npoint]; |
2228 |
z = new float[npoint]; |
z = new float[npoint]; |
2229 |
thx = new float[npoint]; |
thx = new float[npoint]; |
2230 |
thy = new float[npoint]; |
thy = new float[npoint]; |
2231 |
tl = new float[npoint]; |
tl = new float[npoint]; |
2232 |
float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1); |
float dz = ((ZTRK1)-(ZTRK6))/(npoint-1); |
2233 |
for(int i=0; i<npoint; i++){ |
for(int i=0; i<npoint; i++){ |
2234 |
x[i] = 0; |
x[i] = 0; |
2235 |
y[i] = 0; |
y[i] = 0; |
2236 |
z[i] = (ZTRKUP) - i*dz; |
z[i] = (ZTRK1) - i*dz; |
2237 |
thx[i] = 0; |
thx[i] = 0; |
2238 |
thy[i] = 0; |
thy[i] = 0; |
2239 |
tl[i] = 0; |
tl[i] = 0; |
2260 |
thx = new float[npoint]; |
thx = new float[npoint]; |
2261 |
thy = new float[npoint]; |
thy = new float[npoint]; |
2262 |
tl = new float[npoint]; |
tl = new float[npoint]; |
2263 |
float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1); |
float dz = ((ZTRK1)-(ZTRK6))/(npoint-1); |
2264 |
for(int i=0; i<npoint; i++){ |
for(int i=0; i<npoint; i++){ |
2265 |
x[i] = 0; |
x[i] = 0; |
2266 |
y[i] = 0; |
y[i] = 0; |
2267 |
z[i] = (ZTRKUP) - i*dz; |
z[i] = (ZTRK1) - i*dz; |
2268 |
thx[i] = 0; |
thx[i] = 0; |
2269 |
thy[i] = 0; |
thy[i] = 0; |
2270 |
tl[i] = 0; |
tl[i] = 0; |
2289 |
thy = new float[npoint]; |
thy = new float[npoint]; |
2290 |
tl = new float[npoint]; |
tl = new float[npoint]; |
2291 |
int i=0; |
int i=0; |
2292 |
do{ |
do{ |
2293 |
x[i] = 0; |
x[i] = 0.; |
2294 |
y[i] = 0; |
y[i] = 0.; |
2295 |
z[i] = zin[i]; |
z[i] = zin[i]; |
2296 |
thx[i] = 0; |
thx[i] = 0.; |
2297 |
thy[i] = 0; |
thy[i] = 0.; |
2298 |
tl[i] = 0; |
tl[i] = 0.; |
2299 |
i++; |
i++; |
2300 |
}while(zin[i-1] > zin[i] && i < npoint); |
}while(zin[i-1] > zin[i] && i < npoint); |
2301 |
npoint=i; |
npoint=i; |
2302 |
if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl; |
if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points instean of "<<n<<endl; |
2303 |
|
// Dump(); |
2304 |
|
} |
2305 |
|
void Trajectory::Delete(){ |
2306 |
|
|
2307 |
|
if(x) delete [] x; |
2308 |
|
if(y) delete [] y; |
2309 |
|
if(z) delete [] z; |
2310 |
|
if(thx) delete [] thx; |
2311 |
|
if(thy) delete [] thy; |
2312 |
|
if(tl) delete [] tl; |
2313 |
|
|
2314 |
} |
} |
2315 |
//-------------------------------------- |
//-------------------------------------- |
2316 |
// |
// |
2349 |
return l; |
return l; |
2350 |
|
|
2351 |
} |
} |
2352 |
|
|
2353 |
|
/** |
2354 |
|
* Evaluates the trajectory in the apparatus associated to the track. |
2355 |
|
* It integrates the equations of motion in the magnetic field. |
2356 |
|
* @param al Track state-vector (X0,Y0,sin(theta),phi,deflection). |
2357 |
|
* @param zini z-coordinate of the reference plane (Z0). |
2358 |
|
* @return error flag. |
2359 |
|
* |
2360 |
|
* This method is needed when you want to integrate the particle trajectory |
2361 |
|
* starting from a track state-vector relative to an arbitrary reference plane. |
2362 |
|
* The default reference plane, used by the tracker routines, is at zini=23.5. |
2363 |
|
* If you give as input the track state-vector from a TrkTrack object, |
2364 |
|
* you can use Trajectory::DoTrack(float* al) instead. |
2365 |
|
*/ |
2366 |
|
int Trajectory::DoTrack(float* al, float zini){ |
2367 |
|
|
2368 |
|
// double *dxout = new double[npoint]; |
2369 |
|
// double *dyout = new double[npoint]; |
2370 |
|
// double *dthxout = new double[npoint]; |
2371 |
|
// double *dthyout = new double[npoint]; |
2372 |
|
// double *dtlout = new double[npoint]; |
2373 |
|
// double *dzin = new double[npoint]; |
2374 |
|
|
2375 |
|
double *dxout; |
2376 |
|
double *dyout; |
2377 |
|
double *dthxout; |
2378 |
|
double *dthyout; |
2379 |
|
double *dtlout; |
2380 |
|
double *dzin; |
2381 |
|
|
2382 |
|
dxout = (double*) malloc(npoint*sizeof(double)); |
2383 |
|
dyout = (double*) malloc(npoint*sizeof(double)); |
2384 |
|
dthxout = (double*) malloc(npoint*sizeof(double)); |
2385 |
|
dthyout = (double*) malloc(npoint*sizeof(double)); |
2386 |
|
dtlout = (double*) malloc(npoint*sizeof(double)); |
2387 |
|
dzin = (double*) malloc(npoint*sizeof(double)); |
2388 |
|
|
2389 |
|
double dal[5]; |
2390 |
|
|
2391 |
|
double dzini = (double)zini; |
2392 |
|
|
2393 |
|
int ifail = 0; |
2394 |
|
|
2395 |
|
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
2396 |
|
for (int i=0; i<npoint; i++) dzin[i] = (double)z[i]; |
2397 |
|
|
2398 |
|
TrkParams::Load(1); |
2399 |
|
if( !TrkParams::IsLoaded(1) ){ |
2400 |
|
cout << "int Trajectory::DoTrack(float* al) --- ERROR --- m.field not loaded"<<endl; |
2401 |
|
return 0; |
2402 |
|
} |
2403 |
|
// dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
2404 |
|
dotrack3_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&dzini,&ifail); |
2405 |
|
|
2406 |
|
for (int i=0; i<npoint; i++){ |
2407 |
|
x[i] = (float)*(dxout+i); |
2408 |
|
y[i] = (float)*(dyout+i); |
2409 |
|
thx[i] = (float)*(dthxout+i); |
2410 |
|
thy[i] = (float)*(dthyout+i); |
2411 |
|
tl[i] = (float)*(dtlout+i); |
2412 |
|
} |
2413 |
|
|
2414 |
|
if(dxout) free( dxout ); |
2415 |
|
if(dyout) free( dyout ); |
2416 |
|
if(dthxout)free( dthxout ); |
2417 |
|
if(dthyout)free( dthyout ); |
2418 |
|
if(dtlout) free( dtlout ); |
2419 |
|
if(dzin) free( dzin ); |
2420 |
|
|
2421 |
|
// delete [] dxout; |
2422 |
|
// delete [] dyout; |
2423 |
|
// delete [] dthxout; |
2424 |
|
// delete [] dthyout; |
2425 |
|
// delete [] dtlout; |
2426 |
|
// delete [] dzin; |
2427 |
|
|
2428 |
|
|
2429 |
|
return ifail; |
2430 |
|
}; |
2431 |
|
|
2432 |
|
/** |
2433 |
|
* |
2434 |
|
* >>> OBSOLETE !!! use Trajectory::DoTrack(float* al, float zini) instead |
2435 |
|
* |
2436 |
|
*/ |
2437 |
|
int Trajectory::DoTrack2(float* al, float zini){ |
2438 |
|
|
2439 |
|
cout << endl; |
2440 |
|
cout << " int Trajectory::DoTrack2(float* al, float zini) --->> NB NB !! this method is going to be eliminated !!! "<<endl; |
2441 |
|
cout << " >>>> replace it with TrkTrack::DoTrack(Trajectory* t) <<<<"<<endl; |
2442 |
|
cout << " (Sorry Wolfgang!! Don't be totally confused!! By Elena)"<<endl; |
2443 |
|
cout << endl; |
2444 |
|
|
2445 |
|
return DoTrack(al,zini); |
2446 |
|
|
2447 |
|
}; |
2448 |
|
|
2449 |
|
|
2450 |
|
|
2451 |
ClassImp(TrkLevel2); |
ClassImp(TrkLevel2); |
2452 |
ClassImp(TrkSinglet); |
ClassImp(TrkSinglet); |
2453 |
ClassImp(TrkTrack); |
ClassImp(TrkTrack); |