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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Fri May 19 13:15:54 2006 UTC revision 1.5 by pam-fi, Fri Jun 30 09:48:15 2006 UTC
# Line 10  using namespace std; Line 10  using namespace std;
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  //--------------------------------------  //--------------------------------------
# Line 17  extern "C" {     Line 18  extern "C" {    
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;
# Line 46  TrkTrack::TrkTrack(){ Line 48  TrkTrack::TrkTrack(){
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++){
# Line 110  int TrkTrack::DoTrack(Trajectory* t){ Line 113  int TrkTrack::DoTrack(Trajectory* t){
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  //--------------------------------------  //--------------------------------------
# Line 195  TrkLevel2::TrkLevel2(){ Line 240  TrkLevel2::TrkLevel2(){
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  //  //
# Line 223  void TrkLevel2::Dump(){ Line 271  void TrkLevel2::Dump(){
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();
# Line 234  void TrkLevel2::FillCommonVar(cTrkLevel2 Line 286  void TrkLevel2::FillCommonVar(cTrkLevel2
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];
# Line 279  void TrkLevel2::FillCommonVar(cTrkLevel2 Line 333  void TrkLevel2::FillCommonVar(cTrkLevel2
333          new(sy[i]) TrkSinglet(*t_singlet);          new(sy[i]) TrkSinglet(*t_singlet);
334          t_singlet->Clear();          t_singlet->Clear();
335          };          };
336            
337            delete t_track;
338            delete t_singlet;
339  }  }
340  //--------------------------------------  //--------------------------------------
341  //  //
# Line 289  void TrkLevel2::Clear(){ Line 346  void TrkLevel2::Clear(){
346      for(Int_t i=0; i<12 ; i++){      for(Int_t i=0; i<12 ; i++){
347          crc[i] = -1;          crc[i] = -1;
348      };      };
349      Track->RemoveAll();  /*    Track->RemoveAll();
350      SingletX->RemoveAll();      SingletX->RemoveAll();
351      SingletY->RemoveAll();      SingletY->RemoveAll();*/
352            // modify to avoid memory leakage
353            Track->Clear();
354            SingletX->Clear();
355            SingletY->Clear();
356  }  }
357  //--------------------------------------  //--------------------------------------
358  //  //
# Line 302  void TrkLevel2::Clear(){ Line 363  void TrkLevel2::Clear(){
363   * 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.
364   */   */
365  TClonesArray *TrkLevel2::GetTracks(){  TClonesArray *TrkLevel2::GetTracks(){
366        TClonesArray *sorted = GetTracks_NFitSorted();
367        return sorted;
368    };
369    
370    TClonesArray *TrkLevel2::GetTracks_Chi2Sorted(){
371    
372      TClonesArray *sorted = new TClonesArray("TrkTrack");      TClonesArray *sorted = new TClonesArray("TrkTrack");
373      TClonesArray &t = *Track;      TClonesArray &t = *Track;
# Line 330  TClonesArray *TrkLevel2::GetTracks(){ Line 396  TClonesArray *TrkLevel2::GetTracks(){
396      }      }
397      return sorted;      return sorted;
398  }  }
399    TClonesArray *TrkLevel2::GetTracks_NFitSorted(){
400    
401        TClonesArray *sorted = new TClonesArray("TrkTrack");
402        TClonesArray &t = *Track;
403        TClonesArray &ts = *sorted;
404        int N=this->ntrk();
405        vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;
406    
407        int indo=0;
408        int indi=0;
409        while(N != 0){
410            int nfit =0;
411            float chi2ref=1000000;
412            // first loop to search maximum num. of fit points
413            for(int i=0; i<this->ntrk(); i++){
414                if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){
415                    nfit =    ((TrkTrack *)t[i])->GetNtot();
416    //              cout << "1** "<<i<< " " << nfit<<endl;
417                }
418            }
419            //second loop to search minimum chi2 among selected
420            for(int i=0; i<this->ntrk(); i++){
421                if(    ((TrkTrack *)t[i])->chi2 < chi2ref
422                    && ((TrkTrack *)t[i])->GetNtot()== nfit
423                    && m[i]==1){
424                    chi2ref = ((TrkTrack *)t[i])->chi2;
425                    indi = i;
426    //              cout << "2** "<<i<< " " << nfit <<" "<<chi2ref<<endl;
427                }
428            }
429            if( ((TrkTrack *)t[indi])->HasImage() ){
430                m[((TrkTrack *)t[indi])->image] = 0;
431                N--;
432    
433    //          Int_t nfiti=((TrkTrack *)t[((TrkTrack *)t[indi])->image  ])->GetNtot();
434    //          Float_t chi2i=((TrkTrack *)t[((TrkTrack *)t[indi])->image  ])->chi2;
435                    
436    //          cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl;
437            }
438            new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]);
439            m[indi] = 0;
440            N--;    
441            indo++;
442        }
443        return sorted;
444    }
445  //--------------------------------------  //--------------------------------------
446  //  //
447  //  //
# Line 366  TrkTrack *TrkLevel2::GetTrack(int it){ Line 478  TrkTrack *TrkLevel2::GetTrack(int it){
478          return 0;          return 0;
479      }      }
480      TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it];      TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it];
481            GetTracks()->Delete();////TEMPORANEO
482      return track;      return track;
483  }  }
484    Int_t TrkLevel2::GetNTracks(){
485            Int_t ntot=0;
486            ntot = GetTracks()->GetEntries();
487            GetTracks()->Delete();////TEMPORANEO
488            return ntot;
489    };
490  //--------------------------------------  //--------------------------------------
491  //  //
492  //  //
# Line 408  void TrkLevel2::LoadField(TString s){ Line 527  void TrkLevel2::LoadField(TString s){
527  //  //
528  //--------------------------------------  //--------------------------------------
529  /**  /**
530     * Trajectory default constructor.
531     * (By default is created with z-coordinates inside the tracking volume)
532      */
533    Trajectory::Trajectory(){
534        npoint = 10;
535        x = new float[npoint];
536        y = new float[npoint];
537        z = new float[npoint];
538        thx = new float[npoint];
539        thy = new float[npoint];
540        tl = new float[npoint];
541        float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1);
542        for(int i=0; i<npoint; i++){
543            x[i] = 0;
544            y[i] = 0;
545            z[i] = (ZTRKUP) - i*dz;
546            thx[i] = 0;
547            thy[i] = 0;
548            tl[i] = 0;
549        }
550    }
551    //--------------------------------------
552    //
553    //
554    //--------------------------------------
555    /**
556   * Trajectory constructor.   * Trajectory constructor.
557     * (By default is created with z-coordinates inside the tracking volume)
558   * \param n Number of points   * \param n Number of points
559   */   */
560  Trajectory::Trajectory(int n){  Trajectory::Trajectory(int n){
561        if(n<=0){
562            cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl;
563            n=10;
564        }
565      npoint = n;      npoint = n;
566      x = new float[npoint];      x = new float[npoint];
567      y = new float[npoint];      y = new float[npoint];
568      z = new float[npoint];      z = new float[npoint];
569        thx = new float[npoint];
570        thy = new float[npoint];
571        tl = new float[npoint];
572        float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1);
573      for(int i=0; i<npoint; i++){      for(int i=0; i<npoint; i++){
574          x[i] = 0;          x[i] = 0;
575          y[i] = 0;          y[i] = 0;
576          z[i] = 0;          z[i] = (ZTRKUP) - i*dz;
577            thx[i] = 0;
578            thy[i] = 0;
579            tl[i] = 0;
580      }      }
581  }  }
582  //--------------------------------------  //--------------------------------------
# Line 432  Trajectory::Trajectory(int n){ Line 589  Trajectory::Trajectory(int n){
589   * \param pz Pointer to float array, defining z coordinates   * \param pz Pointer to float array, defining z coordinates
590   */   */
591  Trajectory::Trajectory(int n, float* zin){  Trajectory::Trajectory(int n, float* zin){
592      npoint = n;      npoint = 10;
593        if(n>0)npoint = n;
594      x = new float[npoint];      x = new float[npoint];
595      y = new float[npoint];      y = new float[npoint];
596      z = new float[npoint];      z = new float[npoint];
597      for(int i=0; i<npoint; i++){      thx = new float[npoint];
598        thy = new float[npoint];
599        tl = new float[npoint];
600        int i=0;
601        do{
602          x[i] = 0;          x[i] = 0;
603          y[i] = 0;          y[i] = 0;
604          z[i] = zin[i];          z[i] = zin[i];
605      }          thx[i] = 0;
606            thy[i] = 0;
607            tl[i] = 0;
608            i++;            
609        }while(zin[i-1] > zin[i] && i < npoint);
610        npoint=i;
611        if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;
612  }  }
613  //--------------------------------------  //--------------------------------------
614  //  //
# Line 452  Trajectory::Trajectory(int n, float* zin Line 620  Trajectory::Trajectory(int n, float* zin
620  void Trajectory::Dump(){  void Trajectory::Dump(){
621      cout <<endl<< "Trajectory ========== "<<endl;      cout <<endl<< "Trajectory ========== "<<endl;
622      for (int i=0; i<npoint; i++){      for (int i=0; i<npoint; i++){
623          cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] << endl;;          cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ;
624            cout <<" -- " << thx[i] <<" "<< thy[i] ;
625            cout <<" -- " << tl[i] << endl;
626      };      };
627  }  }
628    //--------------------------------------
629    //
630    //
631    //--------------------------------------
632    /**
633     * Get trajectory length between two points
634     * @param ifirst first point (default 0)
635     * @param ilast last point (default npoint)
636     */
637    float Trajectory::GetLength(int ifirst, int ilast){
638        if( ifirst<0    ) ifirst = 0;
639        if( ilast>=npoint) ilast  = npoint-1;
640        float l=0;
641        for(int i=ifirst;i<=ilast;i++){
642            l=l+tl[i];
643        };
644        if(z[ilast] > ZINI)l=l-tl[ilast];
645        if(z[ifirst] < ZINI)   l=l-tl[ifirst];
646        
647        return l;
648    
649    }
650  ClassImp(TrkLevel2);  ClassImp(TrkLevel2);
651  ClassImp(TrkSinglet);  ClassImp(TrkSinglet);
652  ClassImp(TrkTrack);  ClassImp(TrkTrack);
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.23