/[PAMELA software]/DarthVader/TrackerLevel2/src/TrkLevel2.cpp
ViewVC logotype

Contents of /DarthVader/TrackerLevel2/src/TrkLevel2.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (show annotations) (download)
Fri Jun 30 09:22:04 2006 UTC (18 years, 5 months ago) by mocchiut
Branch: MAIN
Changes since 1.3: +1 -0 lines
Memory leak bug fixed

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);

  ViewVC Help
Powered by ViewVC 1.1.23