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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide annotations) (download)
Thu Jun 1 09:03:09 2006 UTC (18 years, 6 months ago) by pam-fi
Branch: MAIN
Changes since 1.1: +121 -6 lines
Tracking routine updated in order to get track projected angles and track length

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

  ViewVC Help
Powered by ViewVC 1.1.23