10 |
//...................................... |
//...................................... |
11 |
extern "C" { |
extern "C" { |
12 |
void dotrack_(int*, double*, double*, double*, double*, int*); |
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*); |
int readb_(const char*); |
15 |
} |
} |
16 |
//-------------------------------------- |
//-------------------------------------- |
18 |
// |
// |
19 |
//-------------------------------------- |
//-------------------------------------- |
20 |
TrkTrack::TrkTrack(){ |
TrkTrack::TrkTrack(){ |
21 |
image = 0; |
seqno = -1; |
22 |
|
image = -1; |
23 |
chi2 = 0; |
chi2 = 0; |
24 |
for(int it1=0;it1<5;it1++){ |
for(int it1=0;it1<5;it1++){ |
25 |
al[it1] = 0; |
al[it1] = 0; |
48 |
// |
// |
49 |
//-------------------------------------- |
//-------------------------------------- |
50 |
TrkTrack::TrkTrack(const TrkTrack& t){ |
TrkTrack::TrkTrack(const TrkTrack& t){ |
51 |
|
seqno = t.seqno; |
52 |
image = t.image; |
image = t.image; |
53 |
chi2 = t.chi2; |
chi2 = t.chi2; |
54 |
for(int it1=0;it1<5;it1++){ |
for(int it1=0;it1<5;it1++){ |
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(){ |
//float TrkTrack::BdL(){ |
159 |
//}; |
//}; |
160 |
//-------------------------------------- |
//-------------------------------------- |
240 |
Track = new TClonesArray("TrkTrack"); |
Track = new TClonesArray("TrkTrack"); |
241 |
SingletX = new TClonesArray("TrkSinglet"); |
SingletX = new TClonesArray("TrkSinglet"); |
242 |
SingletY = new TClonesArray("TrkSinglet"); |
SingletY = new TClonesArray("TrkSinglet"); |
243 |
|
// Track = 0; |
244 |
|
// Singlet = 0; |
245 |
|
// SingletY = 0; |
246 |
} |
} |
247 |
//-------------------------------------- |
//-------------------------------------- |
248 |
// |
// |
271 |
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). |
* Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). |
272 |
*/ |
*/ |
273 |
void TrkLevel2::FillCommonVar(cTrkLevel2 *l2){ |
void TrkLevel2::FillCommonVar(cTrkLevel2 *l2){ |
274 |
|
// |
275 |
|
// Track = new TClonesArray("TrkTrack"); |
276 |
|
// SingletX = new TClonesArray("TrkSinglet"); |
277 |
|
// SingletY = new TClonesArray("TrkSinglet"); |
278 |
// temporary objects: |
// temporary objects: |
279 |
TrkSinglet* t_singlet = new TrkSinglet(); |
TrkSinglet* t_singlet = new TrkSinglet(); |
280 |
TrkTrack* t_track = new TrkTrack(); |
TrkTrack* t_track = new TrkTrack(); |
286 |
// *** TRACKS *** |
// *** TRACKS *** |
287 |
TClonesArray &t = *Track; |
TClonesArray &t = *Track; |
288 |
for(int i=0; i<l2->ntrk; i++){ |
for(int i=0; i<l2->ntrk; i++){ |
289 |
|
t_track->seqno = i; |
290 |
t_track->image = l2->image[i]-1; |
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]; |
t_track->chi2 = l2->chi2_nt[i]; |
293 |
for(int it1=0;it1<5;it1++){ |
for(int it1=0;it1<5;it1++){ |
294 |
t_track->al[it1] = l2->al_nt[i][it1]; |
t_track->al[it1] = l2->al_nt[i][it1]; |
356 |
* 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. |
357 |
*/ |
*/ |
358 |
TClonesArray *TrkLevel2::GetTracks(){ |
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"); |
TClonesArray *sorted = new TClonesArray("TrkTrack"); |
366 |
TClonesArray &t = *Track; |
TClonesArray &t = *Track; |
389 |
} |
} |
390 |
return sorted; |
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 |
// |
// |
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. |
* Trajectory constructor. |
543 |
|
* (By default is created with z-coordinates inside the tracking volume) |
544 |
* \param n Number of points |
* \param n Number of points |
545 |
*/ |
*/ |
546 |
Trajectory::Trajectory(int n){ |
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; |
npoint = n; |
552 |
x = new float[npoint]; |
x = new float[npoint]; |
553 |
y = new float[npoint]; |
y = new float[npoint]; |
554 |
z = new float[npoint]; |
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++){ |
for(int i=0; i<npoint; i++){ |
560 |
x[i] = 0; |
x[i] = 0; |
561 |
y[i] = 0; |
y[i] = 0; |
562 |
z[i] = 0; |
z[i] = (ZTRKUP) - i*dz; |
563 |
|
thx[i] = 0; |
564 |
|
thy[i] = 0; |
565 |
|
tl[i] = 0; |
566 |
} |
} |
567 |
} |
} |
568 |
//-------------------------------------- |
//-------------------------------------- |
575 |
* \param pz Pointer to float array, defining z coordinates |
* \param pz Pointer to float array, defining z coordinates |
576 |
*/ |
*/ |
577 |
Trajectory::Trajectory(int n, float* zin){ |
Trajectory::Trajectory(int n, float* zin){ |
578 |
npoint = n; |
npoint = 10; |
579 |
|
if(n>0)npoint = n; |
580 |
x = new float[npoint]; |
x = new float[npoint]; |
581 |
y = new float[npoint]; |
y = new float[npoint]; |
582 |
z = new float[npoint]; |
z = new float[npoint]; |
583 |
for(int i=0; i<npoint; i++){ |
thx = new float[npoint]; |
584 |
|
thy = new float[npoint]; |
585 |
|
tl = new float[npoint]; |
586 |
|
int i=0; |
587 |
|
do{ |
588 |
x[i] = 0; |
x[i] = 0; |
589 |
y[i] = 0; |
y[i] = 0; |
590 |
z[i] = zin[i]; |
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 |
// |
// |
606 |
void Trajectory::Dump(){ |
void Trajectory::Dump(){ |
607 |
cout <<endl<< "Trajectory ========== "<<endl; |
cout <<endl<< "Trajectory ========== "<<endl; |
608 |
for (int i=0; i<npoint; i++){ |
for (int i=0; i<npoint; i++){ |
609 |
cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] << endl;; |
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); |
ClassImp(TrkLevel2); |
637 |
ClassImp(TrkSinglet); |
ClassImp(TrkSinglet); |
638 |
ClassImp(TrkTrack); |
ClassImp(TrkTrack); |