1 |
/** |
2 |
* \file TrkLevel2.cpp |
3 |
* \author Elena Vannuccini |
4 |
*/ |
5 |
#include <TrkLevel2.h> |
6 |
#include <iostream> |
7 |
using namespace std; |
8 |
//...................................... |
9 |
// F77 routines |
10 |
//...................................... |
11 |
extern "C" { |
12 |
void dotrack_(int*, double*, double*, double*, double*, int*); |
13 |
void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*); |
14 |
int readb_(const char*); |
15 |
} |
16 |
//-------------------------------------- |
17 |
// |
18 |
// |
19 |
//-------------------------------------- |
20 |
TrkTrack::TrkTrack(){ |
21 |
seqno = -1; |
22 |
image = -1; |
23 |
chi2 = 0; |
24 |
for(int it1=0;it1<5;it1++){ |
25 |
al[it1] = 0; |
26 |
for(int it2=0;it2<5;it2++) |
27 |
coval[it1][it2] = 0; |
28 |
}; |
29 |
for(int ip=0;ip<6;ip++){ |
30 |
xgood[ip] = 0; |
31 |
ygood[ip] = 0; |
32 |
xm[ip] = 0; |
33 |
ym[ip] = 0; |
34 |
zm[ip] = 0; |
35 |
resx[ip] = 0; |
36 |
resy[ip] = 0; |
37 |
xv[ip] = 0; |
38 |
yv[ip] = 0; |
39 |
zv[ip] = 0; |
40 |
axv[ip] = 0; |
41 |
ayv[ip] = 0; |
42 |
dedx_x[ip] = 0; |
43 |
dedx_y[ip] = 0; |
44 |
}; |
45 |
}; |
46 |
//-------------------------------------- |
47 |
// |
48 |
// |
49 |
//-------------------------------------- |
50 |
TrkTrack::TrkTrack(const TrkTrack& t){ |
51 |
seqno = t.seqno; |
52 |
image = t.image; |
53 |
chi2 = t.chi2; |
54 |
for(int it1=0;it1<5;it1++){ |
55 |
al[it1] = t.al[it1]; |
56 |
for(int it2=0;it2<5;it2++) |
57 |
coval[it1][it2] = t.coval[it1][it2]; |
58 |
}; |
59 |
for(int ip=0;ip<6;ip++){ |
60 |
xgood[ip] = t.xgood[ip]; |
61 |
ygood[ip] = t.ygood[ip]; |
62 |
xm[ip] = t.xm[ip]; |
63 |
ym[ip] = t.ym[ip]; |
64 |
zm[ip] = t.zm[ip]; |
65 |
resx[ip] = t.resx[ip]; |
66 |
resy[ip] = t.resy[ip]; |
67 |
xv[ip] = t.xv[ip]; |
68 |
yv[ip] = t.yv[ip]; |
69 |
zv[ip] = t.zv[ip]; |
70 |
axv[ip] = t.axv[ip]; |
71 |
ayv[ip] = t.ayv[ip]; |
72 |
dedx_x[ip] = t.dedx_x[ip]; |
73 |
dedx_y[ip] = t.dedx_y[ip]; |
74 |
}; |
75 |
}; |
76 |
//-------------------------------------- |
77 |
// |
78 |
// |
79 |
//-------------------------------------- |
80 |
/** |
81 |
* Evaluates the trajectory in the apparatus associated to the track. |
82 |
* 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. |
83 |
* @param t pointer to an object of the class Trajectory, |
84 |
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
85 |
* @return error flag. |
86 |
*/ |
87 |
int TrkTrack::DoTrack(Trajectory* t){ |
88 |
|
89 |
double *dxout = new double[t->npoint]; |
90 |
double *dyout = new double[t->npoint]; |
91 |
double *dzin = new double[t->npoint]; |
92 |
double dal[5]; |
93 |
|
94 |
int ifail = 0; |
95 |
|
96 |
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
97 |
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; |
98 |
|
99 |
dotrack_(&(t->npoint),dzin,dxout,dyout,dal,&ifail); |
100 |
|
101 |
for (int i=0; i<t->npoint; i++){ |
102 |
t->x[i] = (float)*dxout++; |
103 |
t->y[i] = (float)*dyout++; |
104 |
} |
105 |
|
106 |
// delete [] dxout; |
107 |
// delete [] dyout; |
108 |
// delete [] dzin; |
109 |
|
110 |
return ifail; |
111 |
}; |
112 |
//-------------------------------------- |
113 |
// |
114 |
// |
115 |
//-------------------------------------- |
116 |
/** |
117 |
* Evaluates the trajectory in the apparatus associated to the track. |
118 |
* 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. |
119 |
* @param t pointer to an object of the class Trajectory, |
120 |
* which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ). |
121 |
* @return error flag. |
122 |
*/ |
123 |
int TrkTrack::DoTrack2(Trajectory* t){ |
124 |
|
125 |
double *dxout = new double[t->npoint]; |
126 |
double *dyout = new double[t->npoint]; |
127 |
double *dthxout = new double[t->npoint]; |
128 |
double *dthyout = new double[t->npoint]; |
129 |
double *dtlout = new double[t->npoint]; |
130 |
double *dzin = new double[t->npoint]; |
131 |
double dal[5]; |
132 |
|
133 |
int ifail = 0; |
134 |
|
135 |
for (int i=0; i<5; i++) dal[i] = (double)al[i]; |
136 |
for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i]; |
137 |
|
138 |
dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail); |
139 |
|
140 |
for (int i=0; i<t->npoint; i++){ |
141 |
t->x[i] = (float)*dxout++; |
142 |
t->y[i] = (float)*dyout++; |
143 |
t->thx[i] = (float)*dthxout++; |
144 |
t->thy[i] = (float)*dthyout++; |
145 |
t->tl[i] = (float)*dtlout++; |
146 |
} |
147 |
|
148 |
// delete [] dxout; |
149 |
// delete [] dyout; |
150 |
// delete [] dzin; |
151 |
|
152 |
return ifail; |
153 |
}; |
154 |
//-------------------------------------- |
155 |
// |
156 |
// |
157 |
//-------------------------------------- |
158 |
//float TrkTrack::BdL(){ |
159 |
//}; |
160 |
//-------------------------------------- |
161 |
// |
162 |
// |
163 |
//-------------------------------------- |
164 |
Float_t TrkTrack::GetRigidity(){ |
165 |
Float_t rig=0; |
166 |
if(chi2>0)rig=1./al[4]; |
167 |
if(rig<0) rig=-rig; |
168 |
return rig; |
169 |
}; |
170 |
// |
171 |
Float_t TrkTrack::GetDeflection(){ |
172 |
Float_t def=0; |
173 |
if(chi2>0)def=al[4]; |
174 |
return def; |
175 |
}; |
176 |
// |
177 |
Float_t TrkTrack::GetDEDX(){ |
178 |
Float_t dedx=0; |
179 |
for(Int_t i=0; i<6; i++)dedx+=dedx_x[i]*xgood[i]+dedx_y[i]*ygood[i]; |
180 |
dedx = dedx/(this->GetNX()+this->GetNY()); |
181 |
return dedx; |
182 |
}; |
183 |
|
184 |
//-------------------------------------- |
185 |
// |
186 |
// |
187 |
//-------------------------------------- |
188 |
void TrkTrack::Dump(){ |
189 |
cout << endl << "========== Track " ; |
190 |
cout << endl << "al : "; for(int i=0; i<5; i++)cout << al[i] << " "; |
191 |
cout << endl << "chi^2 : "<< chi2; |
192 |
cout << endl << "xgood : "; for(int i=0; i<6; i++)cout << xgood[i] ; |
193 |
cout << endl << "ygood : "; for(int i=0; i<6; i++)cout << ygood[i] ; |
194 |
cout << endl << "xm : "; for(int i=0; i<6; i++)cout << xm[i] << " "; |
195 |
cout << endl << "ym : "; for(int i=0; i<6; i++)cout << ym[i] << " "; |
196 |
cout << endl << "zm : "; for(int i=0; i<6; i++)cout << zm[i] << " "; |
197 |
cout << endl << "dedx_x : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " "; |
198 |
cout << endl << "dedx_y : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " "; |
199 |
} |
200 |
//-------------------------------------- |
201 |
// |
202 |
// |
203 |
//-------------------------------------- |
204 |
TrkSinglet::TrkSinglet(){ |
205 |
plane = 0; |
206 |
coord[0] = 0; |
207 |
coord[1] = 0; |
208 |
sgnl = 0; |
209 |
}; |
210 |
//-------------------------------------- |
211 |
// |
212 |
// |
213 |
//-------------------------------------- |
214 |
TrkSinglet::TrkSinglet(const TrkSinglet& s){ |
215 |
plane = s.plane; |
216 |
coord[0] = s.coord[0]; |
217 |
coord[1] = s.coord[1]; |
218 |
sgnl = s.sgnl; |
219 |
}; |
220 |
//-------------------------------------- |
221 |
// |
222 |
// |
223 |
//-------------------------------------- |
224 |
void TrkSinglet::Dump(){ |
225 |
int i=0; |
226 |
cout << endl << "========== Singlet " ; |
227 |
cout << endl << "plane : " << plane; |
228 |
cout << endl << "coord[2] : "; while( i<2 && cout << coord[i] << " ") i++; |
229 |
cout << endl << "sgnl : " << sgnl; |
230 |
} |
231 |
//-------------------------------------- |
232 |
// |
233 |
// |
234 |
//-------------------------------------- |
235 |
TrkLevel2::TrkLevel2(){ |
236 |
good2 = -1; |
237 |
for(Int_t i=0; i<12 ; i++){ |
238 |
crc[i] = -1; |
239 |
}; |
240 |
Track = new TClonesArray("TrkTrack"); |
241 |
SingletX = new TClonesArray("TrkSinglet"); |
242 |
SingletY = new TClonesArray("TrkSinglet"); |
243 |
// Track = 0; |
244 |
// Singlet = 0; |
245 |
// SingletY = 0; |
246 |
} |
247 |
//-------------------------------------- |
248 |
// |
249 |
// |
250 |
//-------------------------------------- |
251 |
void TrkLevel2::Dump(){ |
252 |
TClonesArray &t = *Track; |
253 |
TClonesArray &sx = *SingletX; |
254 |
TClonesArray &sy = *SingletY; |
255 |
|
256 |
cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-"; |
257 |
cout << endl << "good2 : " << good2; |
258 |
cout << endl << "crc : "; for(int i=0; i<12; i++) cout << crc[i]; |
259 |
cout << endl << "ntrk() : " << this->ntrk() ; |
260 |
cout << endl << "nclsx() : " << this->nclsx(); |
261 |
cout << endl << "nclsy() : " << this->nclsy(); |
262 |
for(int i=0; i<this->ntrk(); i++) ((TrkTrack *)t[i])->Dump(); |
263 |
for(int i=0; i<this->nclsx(); i++) ((TrkSinglet *)sx[i])->Dump(); |
264 |
for(int i=0; i<this->nclsy(); i++) ((TrkSinglet *)sy[i])->Dump(); |
265 |
} |
266 |
//-------------------------------------- |
267 |
// |
268 |
// |
269 |
//-------------------------------------- |
270 |
/** |
271 |
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). |
272 |
*/ |
273 |
void TrkLevel2::FillCommonVar(cTrkLevel2 *l2){ |
274 |
// |
275 |
// Track = new TClonesArray("TrkTrack"); |
276 |
// SingletX = new TClonesArray("TrkSinglet"); |
277 |
// SingletY = new TClonesArray("TrkSinglet"); |
278 |
// temporary objects: |
279 |
TrkSinglet* t_singlet = new TrkSinglet(); |
280 |
TrkTrack* t_track = new TrkTrack(); |
281 |
// general variables |
282 |
good2 = l2->good2; |
283 |
for(Int_t i=0; i<12 ; i++){ |
284 |
crc[i] = l2->crc[i]; |
285 |
}; |
286 |
// *** TRACKS *** |
287 |
TClonesArray &t = *Track; |
288 |
for(int i=0; i<l2->ntrk; i++){ |
289 |
t_track->seqno = i; |
290 |
t_track->image = l2->image[i]-1; |
291 |
// cout << "track "<<i<<t_track->seqno << t_track->image<<endl; |
292 |
t_track->chi2 = l2->chi2_nt[i]; |
293 |
for(int it1=0;it1<5;it1++){ |
294 |
t_track->al[it1] = l2->al_nt[i][it1]; |
295 |
for(int it2=0;it2<5;it2++) |
296 |
t_track->coval[it1][it2] = l2->coval[i][it2][it1]; |
297 |
}; |
298 |
for(int ip=0;ip<6;ip++){ |
299 |
t_track->xgood[ip] = l2->xgood_nt[i][ip]; |
300 |
t_track->ygood[ip] = l2->ygood_nt[i][ip]; |
301 |
t_track->xm[ip] = l2->xm_nt[i][ip]; |
302 |
t_track->ym[ip] = l2->ym_nt[i][ip]; |
303 |
t_track->zm[ip] = l2->zm_nt[i][ip]; |
304 |
t_track->resx[ip] = l2->resx_nt[i][ip]; |
305 |
t_track->resy[ip] = l2->resy_nt[i][ip]; |
306 |
t_track->xv[ip] = l2->xv_nt[i][ip]; |
307 |
t_track->yv[ip] = l2->yv_nt[i][ip]; |
308 |
t_track->zv[ip] = l2->zv_nt[i][ip]; |
309 |
t_track->axv[ip] = l2->axv_nt[i][ip]; |
310 |
t_track->ayv[ip] = l2->ayv_nt[i][ip]; |
311 |
t_track->dedx_x[ip] = l2->dedx_x[i][ip]; |
312 |
t_track->dedx_y[ip] = l2->dedx_y[i][ip]; |
313 |
}; |
314 |
new(t[i]) TrkTrack(*t_track); |
315 |
t_track->Clear(); |
316 |
}; |
317 |
// *** SINGLETS *** |
318 |
TClonesArray &sx = *SingletX; |
319 |
for(int i=0; i<l2->nclsx; i++){ |
320 |
t_singlet->plane = l2->planex[i]; |
321 |
t_singlet->coord[0] = l2->xs[i][0]; |
322 |
t_singlet->coord[1] = l2->xs[i][1]; |
323 |
t_singlet->sgnl = l2->signlxs[i]; |
324 |
new(sx[i]) TrkSinglet(*t_singlet); |
325 |
t_singlet->Clear(); |
326 |
} |
327 |
TClonesArray &sy = *SingletY; |
328 |
for(int i=0; i<l2->nclsy; i++){ |
329 |
t_singlet->plane = l2->planey[i]; |
330 |
t_singlet->coord[0] = l2->ys[i][0]; |
331 |
t_singlet->coord[1] = l2->ys[i][1]; |
332 |
t_singlet->sgnl = l2->signlys[i]; |
333 |
new(sy[i]) TrkSinglet(*t_singlet); |
334 |
t_singlet->Clear(); |
335 |
}; |
336 |
} |
337 |
//-------------------------------------- |
338 |
// |
339 |
// |
340 |
//-------------------------------------- |
341 |
void TrkLevel2::Clear(){ |
342 |
good2 = -1; |
343 |
for(Int_t i=0; i<12 ; i++){ |
344 |
crc[i] = -1; |
345 |
}; |
346 |
Track->RemoveAll(); |
347 |
SingletX->RemoveAll(); |
348 |
SingletY->RemoveAll(); |
349 |
} |
350 |
//-------------------------------------- |
351 |
// |
352 |
// |
353 |
//-------------------------------------- |
354 |
/** |
355 |
* 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). |
356 |
* This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used. |
357 |
*/ |
358 |
TClonesArray *TrkLevel2::GetTracks(){ |
359 |
TClonesArray *sorted = GetTracks_NFitSorted(); |
360 |
return sorted; |
361 |
}; |
362 |
|
363 |
TClonesArray *TrkLevel2::GetTracks_Chi2Sorted(){ |
364 |
|
365 |
TClonesArray *sorted = new TClonesArray("TrkTrack"); |
366 |
TClonesArray &t = *Track; |
367 |
TClonesArray &ts = *sorted; |
368 |
int N=this->ntrk(); |
369 |
vector<int> m(N); for(int i=0; i<N; i++)m[i]=1; |
370 |
|
371 |
int indo=0; |
372 |
int indi=0; |
373 |
while(N != 0){ |
374 |
float chi2ref=1000000; |
375 |
for(int i=0; i<this->ntrk(); i++){ |
376 |
if(((TrkTrack *)t[i])->chi2 < chi2ref && m[i]==1){ |
377 |
chi2ref = ((TrkTrack *)t[i])->chi2; |
378 |
indi = i; |
379 |
} |
380 |
} |
381 |
if( ((TrkTrack *)t[indi])->image != -1 ){ |
382 |
m[((TrkTrack *)t[indi])->image] = 0; |
383 |
N--; |
384 |
} |
385 |
new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]); |
386 |
m[indi] = 0; |
387 |
N--; |
388 |
indo++; |
389 |
} |
390 |
return sorted; |
391 |
} |
392 |
TClonesArray *TrkLevel2::GetTracks_NFitSorted(){ |
393 |
|
394 |
TClonesArray *sorted = new TClonesArray("TrkTrack"); |
395 |
TClonesArray &t = *Track; |
396 |
TClonesArray &ts = *sorted; |
397 |
int N=this->ntrk(); |
398 |
vector<int> m(N); for(int i=0; i<N; i++)m[i]=1; |
399 |
|
400 |
int indo=0; |
401 |
int indi=0; |
402 |
while(N != 0){ |
403 |
int nfit =0; |
404 |
float chi2ref=1000000; |
405 |
// first loop to search maximum num. of fit points |
406 |
for(int i=0; i<this->ntrk(); i++){ |
407 |
if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){ |
408 |
nfit = ((TrkTrack *)t[i])->GetNtot(); |
409 |
// cout << "1** "<<i<< " " << nfit<<endl; |
410 |
} |
411 |
} |
412 |
//second loop to search minimum chi2 among selected |
413 |
for(int i=0; i<this->ntrk(); i++){ |
414 |
if( ((TrkTrack *)t[i])->chi2 < chi2ref |
415 |
&& ((TrkTrack *)t[i])->GetNtot()== nfit |
416 |
&& m[i]==1){ |
417 |
chi2ref = ((TrkTrack *)t[i])->chi2; |
418 |
indi = i; |
419 |
// cout << "2** "<<i<< " " << nfit <<" "<<chi2ref<<endl; |
420 |
} |
421 |
} |
422 |
if( ((TrkTrack *)t[indi])->HasImage() ){ |
423 |
m[((TrkTrack *)t[indi])->image] = 0; |
424 |
N--; |
425 |
|
426 |
// Int_t nfiti=((TrkTrack *)t[((TrkTrack *)t[indi])->image ])->GetNtot(); |
427 |
// Float_t chi2i=((TrkTrack *)t[((TrkTrack *)t[indi])->image ])->chi2; |
428 |
|
429 |
// cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl; |
430 |
} |
431 |
new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]); |
432 |
m[indi] = 0; |
433 |
N--; |
434 |
indo++; |
435 |
} |
436 |
return sorted; |
437 |
} |
438 |
//-------------------------------------- |
439 |
// |
440 |
// |
441 |
//-------------------------------------- |
442 |
/** |
443 |
* Retrieves the is-th stored track. |
444 |
* @param it Track number, ranging from 0 to ntrk(). |
445 |
* Fitted tracks ( images included ) are stored in a TObjectArray ( TrkLevel2::Track ) in the same order they are returned by the F77 fitting routine. |
446 |
*/ |
447 |
TrkTrack *TrkLevel2::GetStoredTrack(int is){ |
448 |
|
449 |
if(is >= this->ntrk()){ |
450 |
cout << "** TrkLevel2 ** Track "<< is << "doen not exits! " << endl; |
451 |
cout << " Stored tracks ntrk() = "<< this->ntrk() << endl; |
452 |
return 0; |
453 |
} |
454 |
TClonesArray &t = *(Track); |
455 |
TrkTrack *track = (TrkTrack*)t[is]; |
456 |
return track; |
457 |
} |
458 |
//-------------------------------------- |
459 |
// |
460 |
// |
461 |
//-------------------------------------- |
462 |
/** |
463 |
* Retrieves the it-th "physical" track, sorted by the method GetNTracks(). |
464 |
* @param it Track number, ranging from 0 to GetNTracks(). |
465 |
*/ |
466 |
TrkTrack *TrkLevel2::GetTrack(int it){ |
467 |
|
468 |
if(it >= this->GetNTracks()){ |
469 |
cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; |
470 |
cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl; |
471 |
return 0; |
472 |
} |
473 |
TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it]; |
474 |
return track; |
475 |
} |
476 |
//-------------------------------------- |
477 |
// |
478 |
// |
479 |
//-------------------------------------- |
480 |
/** |
481 |
* Retrieves (if present) the image of the it-th "physical" track, sorted by the method GetNTracks(). |
482 |
* @param it Track number, ranging from 0 to GetNTracks(). |
483 |
*/ |
484 |
TrkTrack *TrkLevel2::GetTrackImage(int it){ |
485 |
|
486 |
if(it >= this->GetNTracks()){ |
487 |
cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; |
488 |
cout << " Physical tracks GetNTracks() = "<< this->ntrk() << endl; |
489 |
return 0; |
490 |
} |
491 |
TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it]; |
492 |
if(!track->HasImage()){ |
493 |
cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl; |
494 |
return 0; |
495 |
} |
496 |
TrkTrack *image = (TrkTrack*)(*Track)[track->image]; |
497 |
return image; |
498 |
|
499 |
} |
500 |
//-------------------------------------- |
501 |
// |
502 |
// |
503 |
//-------------------------------------- |
504 |
/** |
505 |
* Loads the magnetic field. |
506 |
* @param s Path of the magnetic-field files. |
507 |
*/ |
508 |
void TrkLevel2::LoadField(TString s){ |
509 |
readb_(s.Data()); |
510 |
}; |
511 |
//-------------------------------------- |
512 |
// |
513 |
// |
514 |
//-------------------------------------- |
515 |
/** |
516 |
* Trajectory default constructor. |
517 |
* (By default is created with z-coordinates inside the tracking volume) |
518 |
*/ |
519 |
Trajectory::Trajectory(){ |
520 |
npoint = 10; |
521 |
x = new float[npoint]; |
522 |
y = new float[npoint]; |
523 |
z = new float[npoint]; |
524 |
thx = new float[npoint]; |
525 |
thy = new float[npoint]; |
526 |
tl = new float[npoint]; |
527 |
float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1); |
528 |
for(int i=0; i<npoint; i++){ |
529 |
x[i] = 0; |
530 |
y[i] = 0; |
531 |
z[i] = (ZTRKUP) - i*dz; |
532 |
thx[i] = 0; |
533 |
thy[i] = 0; |
534 |
tl[i] = 0; |
535 |
} |
536 |
} |
537 |
//-------------------------------------- |
538 |
// |
539 |
// |
540 |
//-------------------------------------- |
541 |
/** |
542 |
* Trajectory constructor. |
543 |
* (By default is created with z-coordinates inside the tracking volume) |
544 |
* \param n Number of points |
545 |
*/ |
546 |
Trajectory::Trajectory(int n){ |
547 |
if(n<=0){ |
548 |
cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl; |
549 |
n=10; |
550 |
} |
551 |
npoint = n; |
552 |
x = new float[npoint]; |
553 |
y = new float[npoint]; |
554 |
z = new float[npoint]; |
555 |
thx = new float[npoint]; |
556 |
thy = new float[npoint]; |
557 |
tl = new float[npoint]; |
558 |
float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1); |
559 |
for(int i=0; i<npoint; i++){ |
560 |
x[i] = 0; |
561 |
y[i] = 0; |
562 |
z[i] = (ZTRKUP) - i*dz; |
563 |
thx[i] = 0; |
564 |
thy[i] = 0; |
565 |
tl[i] = 0; |
566 |
} |
567 |
} |
568 |
//-------------------------------------- |
569 |
// |
570 |
// |
571 |
//-------------------------------------- |
572 |
/** |
573 |
* Trajectory constructor. |
574 |
* \param n Number of points |
575 |
* \param pz Pointer to float array, defining z coordinates |
576 |
*/ |
577 |
Trajectory::Trajectory(int n, float* zin){ |
578 |
npoint = 10; |
579 |
if(n>0)npoint = n; |
580 |
x = new float[npoint]; |
581 |
y = new float[npoint]; |
582 |
z = new float[npoint]; |
583 |
thx = new float[npoint]; |
584 |
thy = new float[npoint]; |
585 |
tl = new float[npoint]; |
586 |
int i=0; |
587 |
do{ |
588 |
x[i] = 0; |
589 |
y[i] = 0; |
590 |
z[i] = zin[i]; |
591 |
thx[i] = 0; |
592 |
thy[i] = 0; |
593 |
tl[i] = 0; |
594 |
i++; |
595 |
}while(zin[i-1] > zin[i] && i < npoint); |
596 |
npoint=i; |
597 |
if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl; |
598 |
} |
599 |
//-------------------------------------- |
600 |
// |
601 |
// |
602 |
//-------------------------------------- |
603 |
/** |
604 |
* Dump the trajectory coordinates. |
605 |
*/ |
606 |
void Trajectory::Dump(){ |
607 |
cout <<endl<< "Trajectory ========== "<<endl; |
608 |
for (int i=0; i<npoint; i++){ |
609 |
cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ; |
610 |
cout <<" -- " << thx[i] <<" "<< thy[i] ; |
611 |
cout <<" -- " << tl[i] << endl; |
612 |
}; |
613 |
} |
614 |
//-------------------------------------- |
615 |
// |
616 |
// |
617 |
//-------------------------------------- |
618 |
/** |
619 |
* Get trajectory length between two points |
620 |
* @param ifirst first point (default 0) |
621 |
* @param ilast last point (default npoint) |
622 |
*/ |
623 |
float Trajectory::GetLength(int ifirst, int ilast){ |
624 |
if( ifirst<0 ) ifirst = 0; |
625 |
if( ilast>=npoint) ilast = npoint-1; |
626 |
float l=0; |
627 |
for(int i=ifirst;i<=ilast;i++){ |
628 |
l=l+tl[i]; |
629 |
}; |
630 |
if(z[ilast] > ZINI)l=l-tl[ilast]; |
631 |
if(z[ifirst] < ZINI) l=l-tl[ifirst]; |
632 |
|
633 |
return l; |
634 |
|
635 |
} |
636 |
ClassImp(TrkLevel2); |
637 |
ClassImp(TrkSinglet); |
638 |
ClassImp(TrkTrack); |
639 |
ClassImp(Trajectory); |