/[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.1.1 by mocchiut, Fri May 19 13:15:54 2006 UTC revision 1.18 by pam-fi, Tue Nov 14 16:28:42 2006 UTC
# Line 4  Line 4 
4   */   */
5  #include <TrkLevel2.h>  #include <TrkLevel2.h>
6  #include <iostream>  #include <iostream>
7    #include <math.h>
8  using namespace std;  using namespace std;
9  //......................................  //......................................
10  // F77 routines  // F77 routines
11  //......................................  //......................................
12  extern "C" {      extern "C" {    
13      void dotrack_(int*, double*, double*, double*, double*, int*);      void dotrack_(int*, double*, double*, double*, double*, int*);
14      int  readb_(const char*);      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
15    //    int  readb_(const char*);
16        int  readb_();
17        void mini2_(int*,int*,int*);
18        void guess_();
19  }  }
20  //--------------------------------------  //--------------------------------------
21  //  //
22  //  //
23  //--------------------------------------  //--------------------------------------
24    /**
25     * Evaluates the trajectory in the apparatus associated to the track.
26     * 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.  
27     * @param t pointer to an object of the class Trajectory,
28     * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
29     * @return error flag.
30     */
31    int Trajectory::DoTrack2(float* al){
32    
33        double *dxout   = new double[npoint];
34        double *dyout   = new double[npoint];
35        double *dthxout = new double[npoint];
36        double *dthyout = new double[npoint];
37        double *dtlout  = new double[npoint];
38        double *dzin    = new double[npoint];
39        double dal[5];
40    
41        int ifail = 0;
42    
43        for (int i=0; i<5; i++)      dal[i]  = (double)al[i];
44        for (int i=0; i<npoint; i++) dzin[i] = (double)z[i];
45    
46        dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
47        
48        for (int i=0; i<npoint; i++){
49            x[i]   = (float)*dxout++;
50            y[i]   = (float)*dyout++;
51            thx[i] = (float)*dthxout++;
52            thy[i] = (float)*dthyout++;
53            tl[i]  = (float)*dtlout++;
54        }
55    
56        return ifail;
57    };
58    //---------------------------------------------
59    //---------------------------------------------
60  TrkTrack::TrkTrack(){  TrkTrack::TrkTrack(){
61      image = 0;      seqno = -1;
62        image = -1;
63      chi2  = 0;      chi2  = 0;
64        nstep = 0;
65      for(int it1=0;it1<5;it1++){      for(int it1=0;it1<5;it1++){
66          al[it1] = 0;          al[it1] = 0;
67          for(int it2=0;it2<5;it2++)          for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
             coval[it1][it2] = 0;  
68      };      };
69      for(int ip=0;ip<6;ip++){      for(int ip=0;ip<6;ip++){
70          xgood[ip]  = 0;                  xgood[ip]  = 0;
71          ygood[ip]  = 0;                  ygood[ip]  = 0;
72          xm[ip]     = 0;                  xm[ip]     = 0;
73          ym[ip]     = 0;                  ym[ip]     = 0;
74          zm[ip]     = 0;                  zm[ip]     = 0;
75          resx[ip]   = 0;                  resx[ip]   = 0;
76          resy[ip]   = 0;                  resy[ip]   = 0;
77          xv[ip]     = 0;                  xv[ip]     = 0;
78          yv[ip]     = 0;                  yv[ip]     = 0;
79          zv[ip]     = 0;                  zv[ip]     = 0;
80          axv[ip]    = 0;                  axv[ip]    = 0;
81          ayv[ip]    = 0;                  ayv[ip]    = 0;
82          dedx_x[ip] = 0;                  dedx_x[ip] = 0;
83          dedx_y[ip] = 0;                  dedx_y[ip] = 0;
84      };      //              clx[ip]    = 0;
85    //              cly[ip]    = 0;
86            };
87            clx = new TRefArray(6,0);
88            cly = new TRefArray(6,0);
89  };  };
90  //--------------------------------------  //--------------------------------------
91  //  //
92  //  //
93  //--------------------------------------  //--------------------------------------
94  TrkTrack::TrkTrack(const TrkTrack& t){  TrkTrack::TrkTrack(const TrkTrack& t){
95        seqno = t.seqno;
96      image = t.image;      image = t.image;
97      chi2  = t.chi2;      chi2  = t.chi2;
98        nstep = t.nstep;
99      for(int it1=0;it1<5;it1++){      for(int it1=0;it1<5;it1++){
100          al[it1] = t.al[it1];          al[it1] = t.al[it1];
101          for(int it2=0;it2<5;it2++)          for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
             coval[it1][it2] = t.coval[it1][it2];  
102      };      };
103      for(int ip=0;ip<6;ip++){      for(int ip=0;ip<6;ip++){
104          xgood[ip]  = t.xgood[ip];                  xgood[ip]  = t.xgood[ip];
105          ygood[ip]  = t.ygood[ip];                  ygood[ip]  = t.ygood[ip];
106          xm[ip]     = t.xm[ip];                  xm[ip]     = t.xm[ip];
107          ym[ip]     = t.ym[ip];                  ym[ip]     = t.ym[ip];
108          zm[ip]     = t.zm[ip];                  zm[ip]     = t.zm[ip];
109          resx[ip]   = t.resx[ip];                  resx[ip]   = t.resx[ip];
110          resy[ip]   = t.resy[ip];                  resy[ip]   = t.resy[ip];
111          xv[ip]     = t.xv[ip];                  xv[ip]     = t.xv[ip];
112          yv[ip]     = t.yv[ip];                  yv[ip]     = t.yv[ip];
113          zv[ip]     = t.zv[ip];                  zv[ip]     = t.zv[ip];
114          axv[ip]    = t.axv[ip];                  axv[ip]    = t.axv[ip];
115          ayv[ip]    = t.ayv[ip];                  ayv[ip]    = t.ayv[ip];
116          dedx_x[ip] = t.dedx_x[ip];                  dedx_x[ip] = t.dedx_x[ip];
117          dedx_y[ip] = t.dedx_y[ip];                  dedx_y[ip] = t.dedx_y[ip];
118      };              //      clx[ip]    = 0;//<<<<pointer
119            //      cly[ip]    = 0;//<<<<pointer
120            };
121            clx = new TRefArray(*(t.clx));
122            cly = new TRefArray(*(t.cly));
123            
124  };  };
125  //--------------------------------------  //--------------------------------------
126  //  //
# Line 110  int TrkTrack::DoTrack(Trajectory* t){ Line 162  int TrkTrack::DoTrack(Trajectory* t){
162  //  //
163  //  //
164  //--------------------------------------  //--------------------------------------
165    /**
166     * Evaluates the trajectory in the apparatus associated to the track.
167     * 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.  
168     * @param t pointer to an object of the class Trajectory,
169     * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
170     * @return error flag.
171     */
172    int TrkTrack::DoTrack2(Trajectory* t){
173    
174        double *dxout   = new double[t->npoint];
175        double *dyout   = new double[t->npoint];
176        double *dthxout = new double[t->npoint];
177        double *dthyout = new double[t->npoint];
178        double *dtlout  = new double[t->npoint];
179        double *dzin    = new double[t->npoint];
180        double dal[5];
181    
182        int ifail = 0;
183    
184        for (int i=0; i<5; i++)         dal[i]  = (double)al[i];
185        for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
186    
187        dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
188        
189        for (int i=0; i<t->npoint; i++){
190            t->x[i]   = (float)*dxout++;
191            t->y[i]   = (float)*dyout++;
192            t->thx[i] = (float)*dthxout++;
193            t->thy[i] = (float)*dthyout++;
194            t->tl[i]  = (float)*dtlout++;
195        }
196    
197    //    delete [] dxout;
198    //    delete [] dyout;
199    //    delete [] dzin;
200    
201        return ifail;
202    };
203    //--------------------------------------
204    //
205    //
206    //--------------------------------------
207  //float TrkTrack::BdL(){  //float TrkTrack::BdL(){
208  //};  //};
209  //--------------------------------------  //--------------------------------------
# Line 142  Float_t TrkTrack::GetDEDX(){ Line 236  Float_t TrkTrack::GetDEDX(){
236  //--------------------------------------  //--------------------------------------
237  void TrkTrack::Dump(){  void TrkTrack::Dump(){
238      cout << endl << "========== Track " ;      cout << endl << "========== Track " ;
239        cout << endl << "seq.  n. : "<< seqno;
240        cout << endl << "image n. : "<< image;
241      cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " ";      cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " ";
242      cout << endl << "chi^2    : "<< chi2;      cout << endl << "chi^2    : "<< chi2;
243        cout << endl << "n.step   : "<< nstep;
244      cout << endl << "xgood    : "; for(int i=0; i<6; i++)cout << xgood[i] ;      cout << endl << "xgood    : "; for(int i=0; i<6; i++)cout << xgood[i] ;
245      cout << endl << "ygood    : "; for(int i=0; i<6; i++)cout << ygood[i] ;      cout << endl << "ygood    : "; for(int i=0; i<6; i++)cout << ygood[i] ;
246      cout << endl << "xm       : "; for(int i=0; i<6; i++)cout << xm[i] << " ";      cout << endl << "xm       : "; for(int i=0; i<6; i++)cout << xm[i] << " ";
247      cout << endl << "ym       : "; for(int i=0; i<6; i++)cout << ym[i] << " ";      cout << endl << "ym       : "; for(int i=0; i<6; i++)cout << ym[i] << " ";
248      cout << endl << "zm       : "; for(int i=0; i<6; i++)cout << zm[i] << " ";      cout << endl << "zm       : "; for(int i=0; i<6; i++)cout << zm[i] << " ";
249        cout << endl << "xv       : "; for(int i=0; i<6; i++)cout << xv[i] << " ";
250        cout << endl << "yv       : "; for(int i=0; i<6; i++)cout << yv[i] << " ";
251        cout << endl << "zv       : "; for(int i=0; i<6; i++)cout << zv[i] << " ";
252        cout << endl << "resx     : "; for(int i=0; i<6; i++)cout << resx[i] << " ";
253        cout << endl << "resy     : "; for(int i=0; i<6; i++)cout << resy[i] << " ";
254        cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
255        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
256        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
257        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
258        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
259      cout << endl << "dedx_x   : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";      cout << endl << "dedx_x   : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";
260      cout << endl << "dedx_y   : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";      cout << endl << "dedx_y   : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";
261        cout << endl;
262    }
263    /**
264     * Set the TrkTrack position measurements
265     */
266    void TrkTrack::SetMeasure(double *xmeas, double *ymeas, double *zmeas){
267        for(int i=0; i<6; i++) xm[i]=*xmeas++;
268        for(int i=0; i<6; i++) ym[i]=*ymeas++;
269        for(int i=0; i<6; i++) zm[i]=*zmeas++;
270    }
271    /**
272     * Set the TrkTrack position resolution
273     */
274    void TrkTrack::SetResolution(double *rx, double *ry){
275        for(int i=0; i<6; i++) resx[i]=*rx++;
276        for(int i=0; i<6; i++) resy[i]=*ry++;
277    }
278    /**
279     * Set the TrkTrack good measurement
280     */
281    void TrkTrack::SetGood(int *xg, int *yg){
282        for(int i=0; i<6; i++) xgood[i]=*xg++;
283        for(int i=0; i<6; i++) ygood[i]=*yg++;
284  }  }
285    
286    /**
287     * Load the magnetic field
288     */
289    void TrkTrack::LoadField(TString path){
290        
291        strcpy(path_.path,path.Data());
292        path_.pathlen = path.Length();
293        path_.error   = 0;
294        readb_();
295    
296    };
297    /**
298     * Tracking method. It calls F77 mini routine.
299     */
300    void TrkTrack::Fit(double pfixed, int& fail, int iprint){
301    
302        float al_ini[] = {0.,0.,0.,0.,0.};
303    
304        extern cMini2track track_;
305        fail = 0;
306    
307        for(int i=0; i<6; i++) track_.xm[i]=xm[i];
308        for(int i=0; i<6; i++) track_.ym[i]=ym[i];
309        for(int i=0; i<6; i++) track_.zm[i]=zm[i];
310        for(int i=0; i<6; i++) track_.resx[i]=resx[i];
311        for(int i=0; i<6; i++) track_.resy[i]=resy[i];
312        for(int i=0; i<6; i++) track_.xgood[i]=xgood[i];
313        for(int i=0; i<6; i++) track_.ygood[i]=ygood[i];
314    
315    // initial guess of "al" with linear fit
316    //     if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.){
317    //      cout << "initial guess "<<endl;
318    //      double szz=0., szx=0., szy=0., ssx=0., ssy=0., sz=0., s1=0.;
319    //      double det, ax, ay, bx, by;
320    //      for(int i=0; i<NPLANE; i++) {
321    //          szz=szz+zm[i]*zm[i];
322    //          szx=szx+zm[i]*xm[i];
323    //          szy=szy+zm[i]*ym[i];
324    //          ssx=ssx+xm[i];
325    //          ssy=ssy+ym[i];
326    //          sz=sz+zm[i];
327    //          s1=s1+1.;
328    //      }
329    //      det=szz*s1-sz*sz;
330    //      ax=(szx*s1-sz*ssx)/det;
331    //      bx=(szz*ssx-szx*sz)/det;
332    //      ay=(szy*s1-sz*ssy)/det;
333    //      by=(szz*ssy-szy*sz)/det;
334    //      al[0]=ax*23.5+bx; // ZINI = 23.5 !!! it should be the same parameter in all codes
335    //      al[1]=ay*23.5+by; //  "  
336    //      al[2]= sqrt(pow(ax,2)+pow(ay,2))/ sqrt(pow(ax,2)+pow(ay,2)+1.);
337    //      al[3]=0.;
338    //      if( (ax!=0.)||(ay!=0.) ) {
339    //          al[3]= asin(ay/ sqrt(pow(ax,2)+pow(ay,2)));
340    //          if(ax<0.) al[3]=acos(-1.)-al[3];
341    //      }
342    //      al[4]=0.;
343    //     }
344    // end guess
345    
346        for(int i=0; i<5; i++) track_.al[i]=al[i];
347        track_.zini = 23.5;
348    // ZINI = 23.5 !!! it should be the same parameter in all codes
349    
350        if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_();
351    
352        // --------------------- free momentum
353        if(pfixed==0.) {
354    //      al[4]=0.;         // free momentum
355            track_.pfixed=0.; //         "
356        }
357        // --------------------- fixed momentum
358        if(pfixed!=0.) {
359            al[4]=1./pfixed;      // to fix the momentum
360            track_.pfixed=pfixed; //         "
361        }
362    
363    //  store temporarily the initial guess
364        for(int i=0; i<5; i++) al_ini[i]=track_.al[i];
365    
366    
367        int istep=0;
368        int ifail=0;
369        mini2_(&istep,&ifail, &iprint);
370        if(ifail!=0) {
371            if(iprint==1)cout << "ERROR: ifail= " << ifail << endl;
372            fail = 1;
373    //      return;
374        }
375        
376    
377    //    cout << endl << "eta ===> " << track_.al[4] << endl;
378    
379        for(int i=0; i<5; i++) al[i]=track_.al[i];
380        chi2=track_.chi2;
381        nstep=track_.nstep;
382        for(int i=0; i<6; i++) xv[i]=track_.xv[i];
383        for(int i=0; i<6; i++) yv[i]=track_.yv[i];
384        for(int i=0; i<6; i++) zv[i]=track_.zv[i];
385        for(int i=0; i<6; i++) axv[i]=track_.axv[i];
386        for(int i=0; i<6; i++) ayv[i]=track_.ayv[i];
387        for(int i=0; i<5; i++) {
388            for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j];
389        }
390    
391        if(fail){
392            cout << " >>>> fit failed >>>> drawing initial par"<<endl;
393            for(int i=0; i<5; i++) al[i]=al_ini[i];
394        }
395    
396    };
397    /*
398     * Reset the fit parameters
399     */
400    void TrkTrack::FitReset(){
401        for(int i=0; i<5; i++) al[i]=-9999.;
402        chi2=0.;
403        nstep=0;
404        for(int i=0; i<6; i++) xv[i]=0.;
405        for(int i=0; i<6; i++) yv[i]=0.;
406        for(int i=0; i<6; i++) zv[i]=0.;
407        for(int i=0; i<6; i++) axv[i]=0.;
408        for(int i=0; i<6; i++) ayv[i]=0.;
409        for(int i=0; i<5; i++) {
410            for(int j=0; j<5; j++) coval[i][j]=0.;
411        }
412    }
413    
414    //--------------------------------------
415    //
416    //
417    //--------------------------------------
418    void TrkTrack::Clear(){
419            seqno = -1;
420            image = -1;
421            chi2  = 0;
422            nstep = 0;
423            for(int it1=0;it1<5;it1++){
424                    al[it1] = 0;
425                    for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
426            };
427            for(int ip=0;ip<6;ip++){
428                    xgood[ip]  = 0;
429                    ygood[ip]  = 0;
430                    xm[ip]     = 0;
431                    ym[ip]     = 0;
432                    zm[ip]     = 0;
433                    resx[ip]   = 0;
434                    resy[ip]   = 0;
435                    xv[ip]     = 0;
436                    yv[ip]     = 0;
437                    zv[ip]     = 0;
438                    axv[ip]    = 0;
439                    ayv[ip]    = 0;
440                    dedx_x[ip] = 0;
441                    dedx_y[ip] = 0;
442    //              clx[ip]    = 0;
443    //              cly[ip]    = 0;
444            };
445            clx->Clear();
446            cly->Clear();
447    };
448    //--------------------------------------
449    //
450    //
451    //--------------------------------------
452    void TrkTrack::Delete(){
453            Clear();
454            clx->Delete();
455            cly->Delete();
456    };
457            //--------------------------------------
458    //
459    //
460    //--------------------------------------
461    
462  //--------------------------------------  //--------------------------------------
463  //  //
464  //  //
# Line 161  TrkSinglet::TrkSinglet(){ Line 468  TrkSinglet::TrkSinglet(){
468      coord[0] = 0;      coord[0] = 0;
469      coord[1] = 0;      coord[1] = 0;
470      sgnl     = 0;      sgnl     = 0;
471            cls      = 0;
472  };  };
473  //--------------------------------------  //--------------------------------------
474  //  //
# Line 171  TrkSinglet::TrkSinglet(const TrkSinglet& Line 479  TrkSinglet::TrkSinglet(const TrkSinglet&
479      coord[0] = s.coord[0];      coord[0] = s.coord[0];
480      coord[1] = s.coord[1];      coord[1] = s.coord[1];
481      sgnl     = s.sgnl;      sgnl     = s.sgnl;
482    //      cls      = 0;//<<<<pointer
483            cls      = TRef(s.cls);
484  };  };
485  //--------------------------------------  //--------------------------------------
486  //  //
# Line 188  void TrkSinglet::Dump(){ Line 498  void TrkSinglet::Dump(){
498  //  //
499  //--------------------------------------  //--------------------------------------
500  TrkLevel2::TrkLevel2(){  TrkLevel2::TrkLevel2(){
501      good2    = -1;  //    good2    = -1;
502      for(Int_t i=0; i<12 ; i++){      for(Int_t i=0; i<12 ; i++){
503          crc[i] = -1;  //      crc[i] = -1;
504      };                  good[i] = -1;
505            };
506      Track    = new TClonesArray("TrkTrack");      Track    = new TClonesArray("TrkTrack");
507      SingletX = new TClonesArray("TrkSinglet");      SingletX = new TClonesArray("TrkSinglet");
508      SingletY = new TClonesArray("TrkSinglet");      SingletY = new TClonesArray("TrkSinglet");
509    
510    //      PhysicalTrack = new TClonesArray("TrkTrack");
511            //sostituire con TRefArray... appena ho capito come si usa
512  }  }
513  //--------------------------------------  //--------------------------------------
514  //  //
515  //  //
516  //--------------------------------------  //--------------------------------------
517  void TrkLevel2::Dump(){  void TrkLevel2::Dump(){
518            
519      TClonesArray &t  = *Track;      TClonesArray &t  = *Track;
520      TClonesArray &sx = *SingletX;      TClonesArray &sx = *SingletX;
521      TClonesArray &sy = *SingletY;      TClonesArray &sy = *SingletY;
522            //
523      cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";      cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";
524      cout << endl << "good2    : " << good2;      cout << endl << "good     : "; for(int i=0; i<12; i++) cout << good[i]<<" ";
     cout << endl << "crc      : "; for(int i=0; i<12; i++) cout << crc[i];  
525      cout << endl << "ntrk()   : " << this->ntrk() ;      cout << endl << "ntrk()   : " << this->ntrk() ;
526      cout << endl << "nclsx()  : " << this->nclsx();      cout << endl << "nclsx()  : " << this->nclsx();
527      cout << endl << "nclsy()  : " << this->nclsy();      cout << endl << "nclsy()  : " << this->nclsy();
# Line 222  void TrkLevel2::Dump(){ Line 536  void TrkLevel2::Dump(){
536  /**  /**
537   * 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).
538   */   */
539  void TrkLevel2::FillCommonVar(cTrkLevel2 *l2){  void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){
540  //  temporary objects:  
541            //  temporary objects:
542      TrkSinglet* t_singlet = new TrkSinglet();      TrkSinglet* t_singlet = new TrkSinglet();
543      TrkTrack*   t_track   = new TrkTrack();      TrkTrack*   t_track   = new TrkTrack();
544  //  general variables  
545      good2 = l2->good2;          //  **** general variables ****
546    //    good2 = l2->good2;
547      for(Int_t i=0; i<12 ; i++){      for(Int_t i=0; i<12 ; i++){
548          crc[i] = l2->crc[i];  //              crc[i] = l2->crc[i];
549      };                  good[i] = l2->good[i];
550  //  *** TRACKS ***          };
551            //  *** TRACKS ***
552      TClonesArray &t = *Track;      TClonesArray &t = *Track;
553      for(int i=0; i<l2->ntrk; i++){      for(int i=0; i<l2->ntrk; i++){
554          t_track->image = l2->image[i]-1;                  t_track->seqno = i;// NBNBNBNB deve sempre essere = i
555          t_track->chi2  = l2->chi2_nt[i];                  t_track->image = l2->image[i]-1;
556          for(int it1=0;it1<5;it1++){          //      cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
557              t_track->al[it1] = l2->al_nt[i][it1];                  t_track->chi2  = l2->chi2_nt[i];
558              for(int it2=0;it2<5;it2++)                  t_track->nstep = l2->nstep_nt[i];
559                  t_track->coval[it1][it2] = l2->coval[i][it2][it1];                  for(int it1=0;it1<5;it1++){
560          };                          t_track->al[it1] = l2->al_nt[i][it1];
561          for(int ip=0;ip<6;ip++){                          for(int it2=0;it2<5;it2++)
562              t_track->xgood[ip]  = l2->xgood_nt[i][ip];                          t_track->coval[it1][it2] = l2->coval[i][it2][it1];
563              t_track->ygood[ip]  = l2->ygood_nt[i][ip];                  };
564              t_track->xm[ip]     = l2->xm_nt[i][ip];                  for(int ip=0;ip<6;ip++){
565              t_track->ym[ip]     = l2->ym_nt[i][ip];                          t_track->xgood[ip]  = l2->xgood_nt[i][ip];
566              t_track->zm[ip]     = l2->zm_nt[i][ip];                          t_track->ygood[ip]  = l2->ygood_nt[i][ip];
567              t_track->resx[ip]   = l2->resx_nt[i][ip];                          t_track->xm[ip]     = l2->xm_nt[i][ip];
568              t_track->resy[ip]   = l2->resy_nt[i][ip];                          t_track->ym[ip]     = l2->ym_nt[i][ip];
569              t_track->xv[ip]     = l2->xv_nt[i][ip];                          t_track->zm[ip]     = l2->zm_nt[i][ip];
570              t_track->yv[ip]     = l2->yv_nt[i][ip];                          t_track->resx[ip]   = l2->resx_nt[i][ip];
571              t_track->zv[ip]     = l2->zv_nt[i][ip];                          t_track->resy[ip]   = l2->resy_nt[i][ip];
572              t_track->axv[ip]    = l2->axv_nt[i][ip];                          t_track->xv[ip]     = l2->xv_nt[i][ip];
573              t_track->ayv[ip]    = l2->ayv_nt[i][ip];                          t_track->yv[ip]     = l2->yv_nt[i][ip];
574              t_track->dedx_x[ip] = l2->dedx_x[i][ip];                          t_track->zv[ip]     = l2->zv_nt[i][ip];
575              t_track->dedx_y[ip] = l2->dedx_y[i][ip];                          t_track->axv[ip]    = l2->axv_nt[i][ip];
576          };                          t_track->ayv[ip]    = l2->ayv_nt[i][ip];
577          new(t[i]) TrkTrack(*t_track);                          t_track->dedx_x[ip] = l2->dedx_x[i][ip];
578          t_track->Clear();                          t_track->dedx_y[ip] = l2->dedx_y[i][ip];
579    //                      t_track->clx[ip] = 0;
580    //                      t_track->cly[ip] = 0;
581                    };
582                    new(t[i]) TrkTrack(*t_track);
583                    t_track->Clear();
584      };      };
585  //  *** SINGLETS ***  //  *** SINGLETS ***
586      TClonesArray &sx = *SingletX;      TClonesArray &sx = *SingletX;
587      for(int i=0; i<l2->nclsx; i++){      for(int i=0; i<l2->nclsx; i++){
588          t_singlet->plane    = l2->planex[i];                  t_singlet->plane    = l2->planex[i];
589          t_singlet->coord[0] = l2->xs[i][0];                  t_singlet->coord[0] = l2->xs[i][0];
590          t_singlet->coord[1] = l2->xs[i][1];                  t_singlet->coord[1] = l2->xs[i][1];
591          t_singlet->sgnl     = l2->signlxs[i];                  t_singlet->sgnl     = l2->signlxs[i];
592          new(sx[i]) TrkSinglet(*t_singlet);  //              t_singlet->cls      = 0;
593          t_singlet->Clear();                  new(sx[i]) TrkSinglet(*t_singlet);
594                    t_singlet->Clear();
595      }      }
596      TClonesArray &sy = *SingletY;      TClonesArray &sy = *SingletY;
597      for(int i=0; i<l2->nclsy; i++){      for(int i=0; i<l2->nclsy; i++){
598          t_singlet->plane    = l2->planey[i];                  t_singlet->plane    = l2->planey[i];
599          t_singlet->coord[0] = l2->ys[i][0];                  t_singlet->coord[0] = l2->ys[i][0];
600          t_singlet->coord[1] = l2->ys[i][1];                  t_singlet->coord[1] = l2->ys[i][1];
601          t_singlet->sgnl     = l2->signlys[i];                  t_singlet->sgnl     = l2->signlys[i];
602          new(sy[i]) TrkSinglet(*t_singlet);  //              t_singlet->cls      = 0;
603          t_singlet->Clear();                  new(sy[i]) TrkSinglet(*t_singlet);
604                    t_singlet->Clear();
605          };          };
606            
607            delete t_track;
608            delete t_singlet;
609  }  }
610  //--------------------------------------  //--------------------------------------
611  //  //
612  //  //
613  //--------------------------------------  //--------------------------------------
614  void TrkLevel2::Clear(){  /**
615      good2    = -1;   * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
616     * Ref to Level1 data (clusters) is also set.
617     */
618    void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){
619    
620    //  temporary objects:
621            TrkSinglet* t_singlet = new TrkSinglet();
622            TrkTrack*   t_track   = new TrkTrack();
623    // general variables
624    //      good2 = l2->good2;
625            for(Int_t i=0; i<12 ; i++){
626    //              crc[i] = l2->crc[i];
627                    good[i] = l2->good[i];
628            };
629    // *** TRACKS ***
630            TClonesArray &t = *Track;
631            for(int i=0; i<l2->ntrk; i++){
632                    t_track->seqno = i;// NBNBNBNB deve sempre essere = i
633                    t_track->image = l2->image[i]-1;
634    //              cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
635                    t_track->chi2  = l2->chi2_nt[i];
636                    t_track->nstep = l2->nstep_nt[i];
637                    for(int it1=0;it1<5;it1++){
638                            t_track->al[it1] = l2->al_nt[i][it1];
639                            for(int it2=0;it2<5;it2++)
640                                    t_track->coval[it1][it2] = l2->coval[i][it2][it1];
641                    };
642                    for(int ip=0;ip<6;ip++){
643                            t_track->xgood[ip]  = l2->xgood_nt[i][ip];
644                            t_track->ygood[ip]  = l2->ygood_nt[i][ip];
645                            t_track->xm[ip]     = l2->xm_nt[i][ip];
646                            t_track->ym[ip]     = l2->ym_nt[i][ip];
647                            t_track->zm[ip]     = l2->zm_nt[i][ip];
648                            t_track->resx[ip]   = l2->resx_nt[i][ip];
649                            t_track->resy[ip]   = l2->resy_nt[i][ip];
650                            t_track->xv[ip]     = l2->xv_nt[i][ip];
651                            t_track->yv[ip]     = l2->yv_nt[i][ip];
652                            t_track->zv[ip]     = l2->zv_nt[i][ip];
653                            t_track->axv[ip]    = l2->axv_nt[i][ip];
654                            t_track->ayv[ip]    = l2->ayv_nt[i][ip];
655                            t_track->dedx_x[ip] = l2->dedx_x[i][ip];
656                            t_track->dedx_y[ip] = l2->dedx_y[i][ip];
657    //                      cout << "traccia "<<i<<"  --  "<< ip << " "<< l2->cltrx[i][ip] <<" "<< l2->cltry[i][ip] <<" "<< t_track->xgood[ip] << t_track->ygood[ip]<<endl;
658                            //-----------------------------------------------------
659    //                      t_track->clx[ip] = l1->GetCluster(l2->cltrx[i][ip]-1);
660    //                      t_track->cly[ip] = l1->GetCluster(l2->cltry[i][ip]-1);
661                            if(t_track->xgood[ip])t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip);
662                            if(t_track->ygood[ip])t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip);
663    //                      if(t_track->ygood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltry[i][ip] "<<l2->cltry[i][ip]<< " l1->GetCluster(l2->cltry[i][ip]-1) "<<l1->GetCluster(l2->cltry[i][ip]-1)<<endl;
664    //                      if(t_track->xgood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltrx[i][ip] "<<l2->cltrx[i][ip]<< " l1->GetCluster(l2->cltrx[i][ip]-1) "<<l1->GetCluster(l2->cltrx[i][ip]-1)<<endl;
665                            //-----------------------------------------------------
666                    };
667                    new(t[i]) TrkTrack(*t_track);
668                    t_track->Clear();
669            };
670    // *** SINGLETS ***
671            TClonesArray &sx = *SingletX;
672            for(int i=0; i<l2->nclsx; i++){
673                    t_singlet->plane    = l2->planex[i];
674                    t_singlet->coord[0] = l2->xs[i][0];
675                    t_singlet->coord[1] = l2->xs[i][1];
676                    t_singlet->sgnl     = l2->signlxs[i];
677                    //-----------------------------------------------------
678    //              cout << "singolo x "<<i<<"  --  "<<  l2->clsx[i] <<endl;
679                    t_singlet->cls      = l1->GetCluster(l2->clsx[i]-1);
680    //              cout<<" i "<<i<<" l2->clsx[i] "<<l2->clsx[i]<< " l1->GetCluster(l2->clsx[i]-1) "<<l1->GetCluster(l2->clsx[i]-1)<<endl;          
681                    //-----------------------------------------------------
682                    new(sx[i]) TrkSinglet(*t_singlet);
683                    t_singlet->Clear();
684            }
685            TClonesArray &sy = *SingletY;
686            for(int i=0; i<l2->nclsy; i++){
687                    t_singlet->plane    = l2->planey[i];
688                    t_singlet->coord[0] = l2->ys[i][0];
689                    t_singlet->coord[1] = l2->ys[i][1];
690                    t_singlet->sgnl     = l2->signlys[i];
691                    //-----------------------------------------------------
692    //              cout << "singolo y "<<i<<"  --  "<<  l2->clsy[i] <<endl;
693                    t_singlet->cls      = l1->GetCluster(l2->clsy[i]-1);
694    //              cout<<" i "<<i<<" l2->clsy[i] "<<l2->clsy[i]<< " l1->GetCluster(l2->clsy[i]-1) "<<l1->GetCluster(l2->clsy[i]-1)<<endl;          
695                    //-----------------------------------------------------
696                    new(sy[i]) TrkSinglet(*t_singlet);
697                    t_singlet->Clear();
698            };
699            
700            delete t_track;
701            delete t_singlet;
702    }
703    /**
704     * Fills a struct cTrkLevel2 with values from a TrkLevel2 object (to put data into a F77 common).
705     */
706    
707    void TrkLevel2::GetLevel2Struct(cTrkLevel2 *l2) const {
708      
709    //  general variables
710    //    l2->good2 = good2 ;
711      for(Int_t i=0; i<12 ; i++){      for(Int_t i=0; i<12 ; i++){
712          crc[i] = -1;  //      l2->crc[i] = crc[i];
713                    l2->good[i] = good[i];
714      };      };
715      Track->RemoveAll();  //  *** TRACKS ***
716    
717        l2->ntrk              =  Track->GetEntries();    
718        for(Int_t i=0;i<l2->ntrk;i++){
719          l2->image[i] = 1 + ((TrkTrack *)Track->At(i))->image;
720          l2->chi2_nt[i] =  ((TrkTrack *)Track->At(i))->chi2;
721              l2->nstep_nt[i] =  ((TrkTrack *)Track->At(i))->nstep;
722              for(int it1=0;it1<5;it1++){
723            l2->al_nt[i][it1] = ((TrkTrack *)Track->At(i))->al[it1];
724            for(int it2=0;it2<5;it2++)
725              l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2];
726          };
727          for(int ip=0;ip<6;ip++){
728            l2->xgood_nt[i][ip] = ((TrkTrack *)Track->At(i))->xgood[ip];
729            l2->ygood_nt[i][ip] = ((TrkTrack *)Track->At(i))->ygood[ip];
730            l2->xm_nt[i][ip]    = ((TrkTrack *)Track->At(i))->xm[ip];
731            l2->ym_nt[i][ip]    = ((TrkTrack *)Track->At(i))->ym[ip];
732            l2->zm_nt[i][ip]    = ((TrkTrack *)Track->At(i))->zm[ip];
733            l2->resx_nt[i][ip]  = ((TrkTrack *)Track->At(i))->resx[ip];
734            l2->resy_nt[i][ip]  = ((TrkTrack *)Track->At(i))->resy[ip];
735            l2->xv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->xv[ip];
736            l2->yv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->yv[ip];
737            l2->zv_nt[i][ip]    = ((TrkTrack *)Track->At(i))->zv[ip];
738            l2->axv_nt[i][ip]   = ((TrkTrack *)Track->At(i))->axv[ip];
739            l2->ayv_nt[i][ip]   = ((TrkTrack *)Track->At(i))->ayv[ip];
740            l2->dedx_x[i][ip]   = ((TrkTrack *)Track->At(i))->dedx_x[ip];
741            l2->dedx_y[i][ip]   = ((TrkTrack *)Track->At(i))->dedx_y[ip];
742          };
743        }
744    
745    //  *** SINGLETS ***    
746        l2->nclsx              = SingletX->GetEntries();
747        for(Int_t i=0;i<l2->nclsx;i++){
748          l2->planex[i]  = ((TrkSinglet *)SingletX->At(i))->plane;
749          l2->xs[i][0]   = ((TrkSinglet *)SingletX->At(i))->coord[0];
750          l2->xs[i][1]   = ((TrkSinglet *)SingletX->At(i))->coord[1];
751          l2->signlxs[i] = ((TrkSinglet *)SingletX->At(i))->sgnl;
752        }
753        l2->nclsy              = SingletY->GetEntries();
754        for(Int_t i=0;i<l2->nclsy;i++){
755          l2->planey[i]  = ((TrkSinglet *)SingletY->At(i))->plane;
756          l2->ys[i][0]   = ((TrkSinglet *)SingletY->At(i))->coord[0];
757          l2->ys[i][1]   = ((TrkSinglet *)SingletY->At(i))->coord[1];
758          l2->signlys[i] = ((TrkSinglet *)SingletY->At(i))->sgnl;
759        }
760    }
761    //--------------------------------------
762    //
763    //
764    //--------------------------------------
765    void TrkLevel2::Clear(){
766    //    good2    = -1;
767        for(Int_t i=0; i<12 ; i++){
768    //      crc[i] = -1;
769                    good[i] = -1;
770            };
771    /*    Track->RemoveAll();
772      SingletX->RemoveAll();      SingletX->RemoveAll();
773      SingletY->RemoveAll();      SingletY->RemoveAll();*/
774            // modify to avoid memory leakage
775            Track->Clear();
776            SingletX->Clear();
777            SingletY->Clear();
778    }
779    //--------------------------------------
780    //
781    //
782    //--------------------------------------
783    void TrkLevel2::Delete(){
784            
785            Clear();
786            Track->Delete();
787            SingletX->Delete();
788            SingletY->Delete();
789  }  }
790  //--------------------------------------  //--------------------------------------
791  //  //
# Line 301  void TrkLevel2::Clear(){ Line 795  void TrkLevel2::Clear(){
795   * 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).   * 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).
796   * 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.
797   */   */
798  TClonesArray *TrkLevel2::GetTracks(){  TRefArray *TrkLevel2::GetTracks_NFitSorted(){
799    
800      TClonesArray *sorted = new TClonesArray("TrkTrack");          TRefArray *sorted = new TRefArray();
801      TClonesArray &t = *Track;          
802      TClonesArray &ts = *sorted;          TClonesArray &t  = *Track;
803      int N=this->ntrk();  //    TClonesArray &ts = *PhysicalTrack;
804      vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;          int N = ntrk();
805            vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;
806      int indo=0;  //      int m[50]; for(int i=0; i<N; i++)m[i]=1;
807      int indi=0;          
808      while(N != 0){          int indo=0;
809          float chi2ref=1000000;          int indi=0;
810          for(int i=0; i<this->ntrk(); i++){          while(N != 0){
811              if(((TrkTrack *)t[i])->chi2 < chi2ref && m[i]==1){                  int nfit =0;
812                  chi2ref = ((TrkTrack *)t[i])->chi2;                  float chi2ref = numeric_limits<float>::max();
813                  indi = i;                  
814              }                  // first loop to search maximum num. of fit points
815                    for(int i=0; i < ntrk(); i++){
816                            if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){
817                                    nfit =    ((TrkTrack *)t[i])->GetNtot();
818                            }
819                    }
820                    //second loop to search minimum chi2 among selected
821                    for(int i=0; i<this->ntrk(); i++){
822                            Float_t chi2 = ((TrkTrack *)t[i])->chi2;
823                            if(chi2 < 0) chi2 = chi2*1000;
824                            if(    chi2 < chi2ref
825                                    && ((TrkTrack *)t[i])->GetNtot() == nfit
826                                    && m[i]==1){
827                                    chi2ref = ((TrkTrack *)t[i])->chi2;
828                                    indi = i;
829                            };
830                    };
831                    if( ((TrkTrack *)t[indi])->HasImage() ){
832                            m[((TrkTrack *)t[indi])->image] = 0;
833                            N--;
834            
835            //          cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl;
836                    };
837                    sorted->Add( (TrkTrack*)t[indi] );      
838                    
839                    m[indi] = 0;
840    //              cout << "SORTED "<< indo << " "<< indi << " "<< N << endl;
841                    N--;    
842                    indo++;
843          }          }
844          if( ((TrkTrack *)t[indi])->image != -1 ){          m.clear();
845              m[((TrkTrack *)t[indi])->image] = 0;  //      cout << "GetTracks_NFitSorted(it): Done"<< endl;
846              N--;  
847          }          return sorted;
848          new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]);  //    return PhysicalTrack;
         m[indi] = 0;  
         N--;      
         indo++;  
     }  
     return sorted;  
849  }  }
850  //--------------------------------------  //--------------------------------------
851  //  //
# Line 355  TrkTrack *TrkLevel2::GetStoredTrack(int Line 872  TrkTrack *TrkLevel2::GetStoredTrack(int
872  //  //
873  //--------------------------------------  //--------------------------------------
874  /**  /**
875     * Retrieves the is-th stored X singlet.
876     * @param it Singlet number, ranging from 0 to nclsx().
877     */
878    TrkSinglet *TrkLevel2::GetSingletX(int is){
879    
880            if(is >= this->nclsx()){
881                    cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;
882                    cout << "                Stored x-singlets nclsx() = "<< this->nclsx() << endl;
883                    return 0;
884            }
885            TClonesArray &t = *(SingletX);
886            TrkSinglet *singlet = (TrkSinglet*)t[is];
887            return singlet;
888    }
889    //--------------------------------------
890    //
891    //
892    //--------------------------------------
893    /**
894     * Retrieves the is-th stored Y singlet.
895     * @param it Singlet number, ranging from 0 to nclsx().
896     */
897    TrkSinglet *TrkLevel2::GetSingletY(int is){
898    
899            if(is >= this->nclsy()){
900                    cout << "** TrkLevel2 ** Singlet "<< is << "doen not exits! " << endl;
901                    cout << "                Stored y-singlets nclsy() = "<< this->nclsx() << endl;
902                    return 0;
903            }
904            TClonesArray &t = *(SingletY);
905            TrkSinglet *singlet = (TrkSinglet*)t[is];
906            return singlet;
907    }
908    //--------------------------------------
909    //
910    //
911    //--------------------------------------
912    /**
913   * Retrieves the it-th "physical" track, sorted by the method GetNTracks().   * Retrieves the it-th "physical" track, sorted by the method GetNTracks().
914   * @param it Track number, ranging from 0 to GetNTracks().   * @param it Track number, ranging from 0 to GetNTracks().
915   */   */
916    
917  TrkTrack *TrkLevel2::GetTrack(int it){  TrkTrack *TrkLevel2::GetTrack(int it){
918            
919      if(it >= this->GetNTracks()){          if(it >= this->GetNTracks()){
920          cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;                  cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
921          cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl;                  cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl;
922          return 0;                  return 0;
923      }          }
924      TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it];          
925      return track;          TRefArray *sorted = GetTracks();  //TEMPORANEO  
926            TrkTrack *track = (TrkTrack*)sorted->At(it);
927            sorted->Delete();
928            return track;
929  }  }
930    /**
931     * Give the number of "physical" tracks, sorted by the method GetTracks().
932     */
933    Int_t TrkLevel2::GetNTracks(){
934                    
935            Float_t ntot=0;
936            TClonesArray &t = *Track;
937            for(int i=0; i<ntrk(); i++) {    
938                    if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.;
939                    else ntot+=0.5;
940            }
941            return (Int_t)ntot;
942    
943    };
944  //--------------------------------------  //--------------------------------------
945  //  //
946  //  //
# Line 378  TrkTrack *TrkLevel2::GetTrack(int it){ Line 951  TrkTrack *TrkLevel2::GetTrack(int it){
951   */   */
952  TrkTrack *TrkLevel2::GetTrackImage(int it){  TrkTrack *TrkLevel2::GetTrackImage(int it){
953    
954      if(it >= this->GetNTracks()){          if(it >= this->GetNTracks()){
955          cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;                  cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl;
956          cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl;                  cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl;
957          return 0;                  return 0;
958      }          }
959      TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it];          
960      if(!track->HasImage()){          TRefArray* sorted = GetTracks(); //TEMPORANEO
961          cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl;          TrkTrack *track = (TrkTrack*)sorted->At(it);
962          return 0;          
963      }          if(!track->HasImage()){
964      TrkTrack *image = (TrkTrack*)(*Track)[track->image];                  cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl;
965      return image;                  return 0;
966            }
967            TrkTrack *image = (TrkTrack*)(*Track)[track->image];
968    
969            sorted->Delete();
970            
971            return image;
972            
973  }  }
974  //--------------------------------------  //--------------------------------------
# Line 400  TrkTrack *TrkLevel2::GetTrackImage(int i Line 979  TrkTrack *TrkLevel2::GetTrackImage(int i
979   * Loads the magnetic field.   * Loads the magnetic field.
980   * @param s Path of the magnetic-field files.   * @param s Path of the magnetic-field files.
981   */   */
982  void TrkLevel2::LoadField(TString s){  void TrkLevel2::LoadField(TString path){
983      readb_(s.Data());  //
984        strcpy(path_.path,path.Data());
985        path_.pathlen = path.Length();
986        path_.error   = 0;
987        readb_();
988    //
989  };  };
990  //--------------------------------------  //--------------------------------------
991  //  //
992  //  //
993  //--------------------------------------  //--------------------------------------
994  /**  /**
995     * Get tracker-plane (mechanical) z-coordinate
996     * @param plane_id plane index (1=TOP,2,3,4,5,6=BOTTOM)
997     */
998    Float_t TrkLevel2::GetZTrk(Int_t plane_id){
999            switch(plane_id){
1000                    case 1: return ZTRK1;
1001                    case 2: return ZTRK2;
1002                    case 3: return ZTRK3;
1003                    case 4: return ZTRK4;
1004                    case 5: return ZTRK5;
1005                    case 6: return ZTRK6;
1006                    default: return 0.;
1007            };
1008    };
1009    //--------------------------------------
1010    //
1011    //
1012    //--------------------------------------
1013    /**
1014     * Trajectory default constructor.
1015     * (By default is created with z-coordinates inside the tracking volume)
1016      */
1017    Trajectory::Trajectory(){
1018        npoint = 10;
1019        x = new float[npoint];
1020        y = new float[npoint];
1021        z = new float[npoint];
1022        thx = new float[npoint];
1023        thy = new float[npoint];
1024        tl = new float[npoint];
1025        float dz = ((ZTRK1)-(ZTRK6))/(npoint-1);
1026        for(int i=0; i<npoint; i++){
1027            x[i] = 0;
1028            y[i] = 0;
1029            z[i] = (ZTRK1) - i*dz;
1030            thx[i] = 0;
1031            thy[i] = 0;
1032            tl[i] = 0;
1033        }
1034    }
1035    //--------------------------------------
1036    //
1037    //
1038    //--------------------------------------
1039    /**
1040   * Trajectory constructor.   * Trajectory constructor.
1041     * (By default is created with z-coordinates inside the tracking volume)
1042   * \param n Number of points   * \param n Number of points
1043   */   */
1044  Trajectory::Trajectory(int n){  Trajectory::Trajectory(int n){
1045        if(n<=0){
1046            cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl;
1047            n=10;
1048        }
1049      npoint = n;      npoint = n;
1050      x = new float[npoint];      x = new float[npoint];
1051      y = new float[npoint];      y = new float[npoint];
1052      z = new float[npoint];      z = new float[npoint];
1053        thx = new float[npoint];
1054        thy = new float[npoint];
1055        tl = new float[npoint];
1056        float dz = ((ZTRK1)-(ZTRK6))/(npoint-1);
1057      for(int i=0; i<npoint; i++){      for(int i=0; i<npoint; i++){
1058          x[i] = 0;          x[i] = 0;
1059          y[i] = 0;          y[i] = 0;
1060          z[i] = 0;          z[i] = (ZTRK1) - i*dz;
1061            thx[i] = 0;
1062            thy[i] = 0;
1063            tl[i] = 0;
1064      }      }
1065  }  }
1066  //--------------------------------------  //--------------------------------------
# Line 432  Trajectory::Trajectory(int n){ Line 1073  Trajectory::Trajectory(int n){
1073   * \param pz Pointer to float array, defining z coordinates   * \param pz Pointer to float array, defining z coordinates
1074   */   */
1075  Trajectory::Trajectory(int n, float* zin){  Trajectory::Trajectory(int n, float* zin){
1076      npoint = n;      npoint = 10;
1077        if(n>0)npoint = n;
1078      x = new float[npoint];      x = new float[npoint];
1079      y = new float[npoint];      y = new float[npoint];
1080      z = new float[npoint];      z = new float[npoint];
1081      for(int i=0; i<npoint; i++){      thx = new float[npoint];
1082          x[i] = 0;      thy = new float[npoint];
1083          y[i] = 0;      tl = new float[npoint];
1084          z[i] = zin[i];      int i=0;
1085      }      do{
1086                    x[i] = 0;
1087                    y[i] = 0;
1088                    z[i] = zin[i];
1089                    thx[i] = 0;
1090                    thy[i] = 0;
1091                    tl[i] = 0;
1092                    i++;            
1093        }while(zin[i-1] > zin[i] && i < npoint);
1094        npoint=i;
1095        if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;
1096  }  }
1097  //--------------------------------------  //--------------------------------------
1098  //  //
# Line 452  Trajectory::Trajectory(int n, float* zin Line 1104  Trajectory::Trajectory(int n, float* zin
1104  void Trajectory::Dump(){  void Trajectory::Dump(){
1105      cout <<endl<< "Trajectory ========== "<<endl;      cout <<endl<< "Trajectory ========== "<<endl;
1106      for (int i=0; i<npoint; i++){      for (int i=0; i<npoint; i++){
1107          cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] << endl;;          cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ;
1108            cout <<" -- " << thx[i] <<" "<< thy[i] ;
1109            cout <<" -- " << tl[i] << endl;
1110      };      };
1111  }  }
1112    //--------------------------------------
1113    //
1114    //
1115    //--------------------------------------
1116    /**
1117     * Get trajectory length between two points
1118     * @param ifirst first point (default 0)
1119     * @param ilast last point (default npoint)
1120     */
1121    float Trajectory::GetLength(int ifirst, int ilast){
1122        if( ifirst<0    ) ifirst = 0;
1123        if( ilast>=npoint) ilast  = npoint-1;
1124        float l=0;
1125        for(int i=ifirst;i<=ilast;i++){
1126            l=l+tl[i];
1127        };
1128        if(z[ilast] > ZINI)l=l-tl[ilast];
1129        if(z[ifirst] < ZINI)   l=l-tl[ifirst];
1130        
1131        return l;
1132    
1133    }
1134    
1135    
1136  ClassImp(TrkLevel2);  ClassImp(TrkLevel2);
1137  ClassImp(TrkSinglet);  ClassImp(TrkSinglet);
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.18
  ViewVC Help
Powered by ViewVC 1.1.23