/[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.11 by pam-fi, Wed Oct 11 06:53:01 2006 UTC revision 1.15 by pam-fi, Tue Nov 7 15:55:11 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
# Line 12  extern "C" {     Line 13  extern "C" {    
13      void dotrack_(int*, double*, double*, double*, double*, int*);      void dotrack_(int*, double*, double*, double*, double*, int*);
14      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
15      int  readb_(const char*);      int  readb_(const char*);
16        void mini2_(int*,int*,int*);
17        void guess_();
18  }  }
19  //--------------------------------------  //--------------------------------------
20  //  //
# Line 21  TrkTrack::TrkTrack(){ Line 24  TrkTrack::TrkTrack(){
24      seqno = -1;      seqno = -1;
25      image = -1;      image = -1;
26      chi2  = 0;      chi2  = 0;
27          nstep = 0;      nstep = 0;
28          for(int it1=0;it1<5;it1++){      for(int it1=0;it1<5;it1++){
29                  al[it1] = 0;          al[it1] = 0;
30                  for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;          for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
31      };      };
32      for(int ip=0;ip<6;ip++){      for(int ip=0;ip<6;ip++){
33                  xgood[ip]  = 0;                  xgood[ip]  = 0;
# Line 55  TrkTrack::TrkTrack(const TrkTrack& t){ Line 58  TrkTrack::TrkTrack(const TrkTrack& t){
58      seqno = t.seqno;      seqno = t.seqno;
59      image = t.image;      image = t.image;
60      chi2  = t.chi2;      chi2  = t.chi2;
61          nstep = t.nstep;      nstep = t.nstep;
62          for(int it1=0;it1<5;it1++){      for(int it1=0;it1<5;it1++){
63                  al[it1] = t.al[it1];          al[it1] = t.al[it1];
64                  for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];          for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
65      };      };
66      for(int ip=0;ip<6;ip++){      for(int ip=0;ip<6;ip++){
67                  xgood[ip]  = t.xgood[ip];                  xgood[ip]  = t.xgood[ip];
# Line 196  Float_t TrkTrack::GetDEDX(){ Line 199  Float_t TrkTrack::GetDEDX(){
199  //--------------------------------------  //--------------------------------------
200  void TrkTrack::Dump(){  void TrkTrack::Dump(){
201      cout << endl << "========== Track " ;      cout << endl << "========== Track " ;
202          cout << endl << "seq.  n. : "<< seqno;      cout << endl << "seq.  n. : "<< seqno;
203          cout << endl << "image n. : "<< image;      cout << endl << "image n. : "<< image;
204          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] << " ";
205      cout << endl << "chi^2    : "<< chi2;      cout << endl << "chi^2    : "<< chi2;
206          cout << endl << "n.step   : "<< nstep;      cout << endl << "n.step   : "<< nstep;
207          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] ;
208      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] ;
209      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] << " ";
210      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] << " ";
211      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] << " ";
212        cout << endl << "xv       : "; for(int i=0; i<6; i++)cout << xv[i] << " ";
213        cout << endl << "yv       : "; for(int i=0; i<6; i++)cout << yv[i] << " ";
214        cout << endl << "zv       : "; for(int i=0; i<6; i++)cout << zv[i] << " ";
215        cout << endl << "resx     : "; for(int i=0; i<6; i++)cout << resx[i] << " ";
216        cout << endl << "resy     : "; for(int i=0; i<6; i++)cout << resy[i] << " ";
217        cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
218        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
219        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
220        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
221        cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
222      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] << " ";
223      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] << " ";
224        cout << endl;
225    }
226    /**
227     * Set the TrkTrack position measurements
228     */
229    void TrkTrack::SetMeasure(double *xmeas, double *ymeas, double *zmeas){
230        for(int i=0; i<6; i++) xm[i]=*xmeas++;
231        for(int i=0; i<6; i++) ym[i]=*ymeas++;
232        for(int i=0; i<6; i++) zm[i]=*zmeas++;
233    }
234    /**
235     * Set the TrkTrack position resolution
236     */
237    void TrkTrack::SetResolution(double *rx, double *ry){
238        for(int i=0; i<6; i++) resx[i]=*rx++;
239        for(int i=0; i<6; i++) resy[i]=*ry++;
240  }  }
241    /**
242     * Set the TrkTrack good measurement
243     */
244    void TrkTrack::SetGood(int *xg, int *yg){
245        for(int i=0; i<6; i++) xgood[i]=*xg++;
246        for(int i=0; i<6; i++) ygood[i]=*yg++;
247    }
248    
249    /**
250     * Load the magnetic field
251     */
252    void TrkTrack::LoadField(TString s){
253        readb_(s.Data());
254    };
255    /**
256     * Tracking method. It calls F77 mini routine.
257     */
258    void TrkTrack::Fit(double pfixed, int& fail, int iprint){
259    
260        float al_ini[] = {0.,0.,0.,0.,0.};
261    
262        extern cMini2track track_;
263        fail = 0;
264    //    extern cMini2fitinfo fit_info_;
265    //    extern void mini_2_(int*,int*);
266    //    track_.xm[0]=1.0;
267    
268        cout << "fitting routine "<<endl;
269    
270        for(int i=0; i<6; i++) track_.xm[i]=xm[i];
271        for(int i=0; i<6; i++) track_.ym[i]=ym[i];
272        for(int i=0; i<6; i++) track_.zm[i]=zm[i];
273        for(int i=0; i<6; i++) track_.resx[i]=resx[i];
274        for(int i=0; i<6; i++) track_.resy[i]=resy[i];
275        for(int i=0; i<6; i++) track_.xgood[i]=xgood[i];
276        for(int i=0; i<6; i++) track_.ygood[i]=ygood[i];
277    
278    // initial guess of "al" with linear fit
279    //     if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.){
280    
281    //      cout << "initial guess "<<endl;
282    
283    //      double szz=0., szx=0., szy=0., ssx=0., ssy=0., sz=0., s1=0.;
284    //      double det, ax, ay, bx, by;
285    //      for(int i=0; i<NPLANE; i++) {
286    //          szz=szz+zm[i]*zm[i];
287    //          szx=szx+zm[i]*xm[i];
288    //          szy=szy+zm[i]*ym[i];
289    //          ssx=ssx+xm[i];
290    //          ssy=ssy+ym[i];
291    //          sz=sz+zm[i];
292    //          s1=s1+1.;
293    //      }
294    //      det=szz*s1-sz*sz;
295    //      ax=(szx*s1-sz*ssx)/det;
296    //      bx=(szz*ssx-szx*sz)/det;
297    //      ay=(szy*s1-sz*ssy)/det;
298    //      by=(szz*ssy-szy*sz)/det;
299    //      al[0]=ax*23.5+bx; // ZINI = 23.5 !!! it should be the same parameter in all codes
300    //      al[1]=ay*23.5+by; //  "  
301    //      al[2]= sqrt(pow(ax,2)+pow(ay,2))/ sqrt(pow(ax,2)+pow(ay,2)+1.);
302    //      al[3]=0.;
303    //      if( (ax!=0.)||(ay!=0.) ) {
304    //          al[3]= asin(ay/ sqrt(pow(ax,2)+pow(ay,2)));
305    //          if(ax<0.) al[3]=acos(-1.)-al[3];
306    //      }
307    //      al[4]=0.;
308    
309    //     }
310    // end guess
311        for(int i=0; i<5; i++) track_.al[i]=al[i];
312        track_.zini = 23.5; // ZINI = 23.5 !!! it should be the same parameter in all codes
313    
314    // elena -----    
315        if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_();
316    // elena -----    
317    
318        if(pfixed==0.) {
319    //      al[4]=0.;         // free momentum
320            track_.pfixed=0.; //         "
321        }
322        if(pfixed!=0.) {
323            al[4]=1./pfixed;      // to fix the momentum
324            track_.pfixed=pfixed; //         "
325        }
326    
327    //  store temporarily the initial guess
328        for(int i=0; i<5; i++) al_ini[i]=track_.al[i];
329    
330    
331        int istep=0;
332        int ifail=0;
333        mini2_(&istep,&ifail, &iprint);
334        if(ifail!=0) {
335            if(iprint==1)cout << "ERROR: ifail= " << ifail << endl;
336            fail = 1;
337    //      return;
338        }
339        
340    
341    //    cout << endl << "eta ===> " << track_.al[4] << endl;
342    
343        for(int i=0; i<5; i++) al[i]=track_.al[i];
344        chi2=track_.chi2;
345        nstep=track_.nstep;
346        for(int i=0; i<6; i++) xv[i]=track_.xv[i];
347        for(int i=0; i<6; i++) yv[i]=track_.yv[i];
348        for(int i=0; i<6; i++) zv[i]=track_.zv[i];
349        for(int i=0; i<6; i++) axv[i]=track_.axv[i];
350        for(int i=0; i<6; i++) ayv[i]=track_.ayv[i];
351        for(int i=0; i<5; i++) {
352            for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j];
353        }
354    
355        if(fail){
356            cout << " >>>> fit failed >>>> drawing initial par"<<endl;
357            for(int i=0; i<5; i++) al[i]=al_ini[i];
358        }
359    
360    };
361    /*
362     * Reset the fit parameters
363     */
364    void TrkTrack::FitReset(){
365        for(int i=0; i<5; i++) al[i]=-9999.;
366        chi2=0.;
367        nstep=0;
368        for(int i=0; i<6; i++) xv[i]=0.;
369        for(int i=0; i<6; i++) yv[i]=0.;
370        for(int i=0; i<6; i++) zv[i]=0.;
371        for(int i=0; i<6; i++) axv[i]=0.;
372        for(int i=0; i<6; i++) ayv[i]=0.;
373        for(int i=0; i<5; i++) {
374            for(int j=0; j<5; j++) coval[i][j]=0.;
375        }
376    }
377    
378  //--------------------------------------  //--------------------------------------
379  //  //
380  //  //
# Line 732  Int_t TrkLevel2::GetNTracks(){ Line 898  Int_t TrkLevel2::GetNTracks(){
898                                    
899          Float_t ntot=0;          Float_t ntot=0;
900          TClonesArray &t = *Track;          TClonesArray &t = *Track;
901          for(int i=0; i<ntrk(); i++) {          for(int i=0; i<ntrk(); i++) {    
902                  if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.;                  if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.;
903                  else ntot+=0.5;                  else ntot+=0.5;
904          }          }
Legend:
Removed from v.1.11  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.23