/[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.32 by pam-fi, Fri Apr 27 10:39:57 2007 UTC revision 1.33 by pam-fi, Mon May 14 11:03:05 2007 UTC
# Line 14  extern "C" {     Line 14  extern "C" {    
14      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);      void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
15       void mini2_(int*,int*,int*);       void mini2_(int*,int*,int*);
16       void guess_();       void guess_();
      void gufld_(float*, float*);  
17  }  }
18    
19  //--------------------------------------  //--------------------------------------
# Line 56  TrkTrack::TrkTrack(){ Line 55  TrkTrack::TrkTrack(){
55  //    clx = TRefArray(6,0);  //    clx = TRefArray(6,0);
56  //    cly = TRefArray(6,0);  //    cly = TRefArray(6,0);
57    
58        TrkParams::SetTrackingMode();
59        TrkParams::SetPrecisionFactor();
60        TrkParams::SetStepMin();
61        TrkParams::SetPFA();
62    
63  };  };
64  //--------------------------------------  //--------------------------------------
65  //  //
# Line 95  TrkTrack::TrkTrack(const TrkTrack& t){ Line 99  TrkTrack::TrkTrack(const TrkTrack& t){
99  //    clx = TRefArray(t.clx);  //    clx = TRefArray(t.clx);
100  //    cly = TRefArray(t.cly);  //    cly = TRefArray(t.cly);
101    
102        TrkParams::SetTrackingMode();
103        TrkParams::SetPrecisionFactor();
104        TrkParams::SetStepMin();  
105        TrkParams::SetPFA();
106    
107  };  };
108  //--------------------------------------  //--------------------------------------
109  //  //
# Line 397  void TrkTrack::SetResolution(double *rx, Line 406  void TrkTrack::SetResolution(double *rx,
406      for(int i=0; i<6; i++) resy[i]=*ry++;      for(int i=0; i<6; i++) resy[i]=*ry++;
407  }  }
408  /**  /**
409   * Set the TrkTrack good measurement   * Set the TrkTrack good measurement.
410   */   */
411  void TrkTrack::SetGood(int *xg, int *yg){  void TrkTrack::SetGood(int *xg, int *yg){
412      // NB! si perdera` l'informazione sul numero del cluster  
413      for(int i=0; i<6; i++) xgood[i]=*xg++;      for(int i=0; i<6; i++) xgood[i]=*xg++;
414      for(int i=0; i<6; i++) ygood[i]=*yg++;      for(int i=0; i<6; i++) ygood[i]=*yg++;
415  }  }
# Line 415  void TrkTrack::LoadField(TString path){ Line 424  void TrkTrack::LoadField(TString path){
424  //     path_.error   = 0;  //     path_.error   = 0;
425  //     readb_();  //     readb_();
426    
427        TrkParams::SetTrackingMode();
428        TrkParams::SetPrecisionFactor();
429        TrkParams::SetStepMin();
430    
431      TrkParams::Set(path,1);      TrkParams::Set(path,1);
432      TrkParams::Load(1);      TrkParams::Load(1);
433    
# Line 428  void TrkTrack::FillMiniStruct(cMini2trac Line 441  void TrkTrack::FillMiniStruct(cMini2trac
441    
442      for(int i=0; i<6; i++){      for(int i=0; i<6; i++){
443    
444  //      track.xgood[i]=xgood[i];  //      cout << i<<" - "<<xgood[i]<<" "<<XGood(i)<<endl;
445  //      track.ygood[i]=ygood[i];  //      cout << i<<" - "<<ygood[i]<<" "<<YGood(i)<<endl;
446          track.xgood[i]=XGood(i);          track.xgood[i]=XGood(i);
447          track.ygood[i]=YGood(i);          track.ygood[i]=YGood(i);
448                    
# Line 486  void TrkTrack::SetFromMiniStruct(cMini2t Line 499  void TrkTrack::SetFromMiniStruct(cMini2t
499            
500  }  }
501  /**  /**
502   * Tracking method. It calls F77 mini routine.   * \brief Method to re-evaluate coordinates of clusters associated with a track.
503     *
504     * The method can be applied only after recovering level1 information
505     * (either by reprocessing single events from level0 or from  
506     * the TrkLevel1 branch, if present); it calls F77 subroutines that
507     * read the level1 common and fill the minimization-routine common.
508     * Some clusters can be excluded or added by means of the methods:
509     *
510     * TrkTrack::ResetXGood(int ip)
511     * TrkTrack::ResetYGood(int ip)
512     * TrkTrack::SetXGood(int ip, int cid, int is)
513     * TrkTrack::SetYGood(int ip, int cid, int is)
514     *
515     * NB! The method TrkTrack::SetGood(int *xg, int *yg) set the plane-mask (0-1)
516     * for the minimization-routine common. It deletes the cluster information
517     * (at least for the moment...) thus cannot be applied before
518     * TrkTrack::EvaluateClusterPositions().  
519     *
520     * Different p.f.a. can be applied by calling (once) the method:
521     *
522     * TrkParams::SetPFA(0); //Set ETA p.f.a.
523     *
524     * @see TrkParams::SetPFA(int)
525     */
526    void TrkTrack::EvaluateClusterPositions(){
527        
528    //     cout << "void TrkTrack::GetClusterPositions() "<<endl;
529    
530        TrkParams::Load( );
531        if( !TrkParams::IsLoaded() )return;
532        
533        for(int ip=0; ip<6; ip++){
534    //      cout << ip<<" ** "<<xm[ip]<<" / "<<ym[ip]<<endl;;
535            int icx = GetClusterX_ID(ip)+1;
536            int icy = GetClusterY_ID(ip)+1;
537            int sensor = GetSensor(ip)+1;//<< convenzione "Paolo"
538            if(ip==5 && sensor!=0)sensor=3-sensor;//<< convenzione "Elena"
539            int ladder = GetLadder(ip)+1;
540            float ax = axv[ip];
541            float ay = ayv[ip];
542            float v[3];
543            v[0]=xv[ip];
544            v[1]=yv[ip];
545            v[2]=zv[ip];
546            float bfx = 10*TrkParams::GetBX(v);//Tesla
547            float bfy = 10*TrkParams::GetBY(v);//Tesla
548            int ipp=ip+1;
549            xyzpam_(&ipp,&icx,&icy,&ladder,&sensor,&ax,&ay,&bfx,&bfy);
550        }
551    }
552    /**
553     * \brief Tracking method. It calls F77 mini routine.
554     *
555     * @param pfixed Particle momentum. If pfixed=0 the momentum
556     * is left as a free parameter, otherwise it is fixed to the input value.
557     * @param fail Output flag (!=0 if the fit failed).
558     * @param iprint Flag to set debug mode ( 0 = no output; 1 = verbose; 2 = debug).
559     * @param froml1 Flag to re-evaluate positions (see TrkTrack::GetClusterPositions()).
560     *
561     * The option to re-evaluate positions can be used only after recovering
562     * level1 information, eg. by reprocessing the single event.
563     *
564     * Example:
565     *
566     * if( !event->GetTrkLevel0() )return false;
567     * event->GetTrkLevel0()->ProcessEvent(); // re-processing level0->level1
568     * int fail=0;
569     * event->GetTrkLevel2()->GetTrack(0)->Fit(0.,fail,0,1);
570     *
571     * @see EvaluateClusterPositions()
572     *
573     * The fitting procedure can be varied by changing the tracking mode,
574     * the fit-precision factor and the minimum number of step.
575     * @see SetTrackingMode(int)
576     * @see SetPrecisionFactor(double)
577     * @see SetStepMin(int)
578   */   */
579  void TrkTrack::Fit(double pfixed, int& fail, int iprint){  void TrkTrack::Fit(double pfixed, int& fail, int iprint, int froml1){
580    
581      float al_ini[] = {0.,0.,0.,0.,0.};      float al_ini[] = {0.,0.,0.,0.,0.};
582    
583        TrkParams::Load( );
584        if( !TrkParams::IsLoaded() )return;
585    
586      extern cMini2track track_;      extern cMini2track track_;
587      fail = 0;      fail = 0;
588      FillMiniStruct(track_);      FillMiniStruct(track_);
589        
590        if(froml1!=0)EvaluateClusterPositions();
591        
592      // if fit variables have been reset, evaluate the initial guess      // if fit variables have been reset, evaluate the initial guess
593      if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_();      if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.)guess_();
594    
# Line 514  void TrkTrack::Fit(double pfixed, int& f Line 607  void TrkTrack::Fit(double pfixed, int& f
607    
608      //  ------------------------------------------      //  ------------------------------------------
609      //  call mini routine      //  call mini routine
610      TrkParams::Load(1);  //     TrkParams::Load(1);
611      if( !TrkParams::IsLoaded(1) ){  //     if( !TrkParams::IsLoaded(1) ){
612          cout << "void TrkTrack::Fit(double pfixed, int& fail, int iprint) --- ERROR --- m.field not loaded"<<endl;  //      cout << "void TrkTrack::Fit(double pfixed, int& fail, int iprint) --- ERROR --- m.field not loaded"<<endl;
613          return;  //      return;
614      }  //     }
615      int istep=0;      int istep=0;
616      int ifail=0;      int ifail=0;
617      mini2_(&istep,&ifail, &iprint);      mini2_(&istep,&ifail, &iprint);
# Line 529  void TrkTrack::Fit(double pfixed, int& f Line 622  void TrkTrack::Fit(double pfixed, int& f
622      //  ------------------------------------------      //  ------------------------------------------
623            
624      SetFromMiniStruct(&track_);      SetFromMiniStruct(&track_);
 //    cout << endl << "eta ===> " << track_.al[4] << endl;  
   
 //     for(int i=0; i<5; i++) al[i]=track_.al[i];  
 //     chi2=track_.chi2;  
 //     nstep=track_.nstep;  
 //     for(int i=0; i<6; i++) xv[i]=track_.xv[i];  
 //     for(int i=0; i<6; i++) yv[i]=track_.yv[i];  
 //     for(int i=0; i<6; i++) zv[i]=track_.zv[i];  
 //     for(int i=0; i<6; i++) axv[i]=track_.axv[i];  
 //     for(int i=0; i<6; i++) ayv[i]=track_.ayv[i];  
 //     for(int i=0; i<5; i++) {  
 //      for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j];  
 //     }  
625    
626      if(fail){      if(fail){
627          if(iprint)cout << " >>>> fit failed >>>> drawing initial par"<<endl;          if(iprint)cout << " >>>> fit failed "<<endl;
628          for(int i=0; i<5; i++) al[i]=al_ini[i];          for(int i=0; i<5; i++) al[i]=al_ini[i];
629      }      }
630    
631  };  };
632  /*  /**
633   * Reset the fit parameters   * Reset the fit parameters
634   */   */
635  void TrkTrack::FitReset(){  void TrkTrack::FitReset(){
636      for(int i=0; i<5; i++) al[i]=-9999.;      for(int i=0; i<5; i++) al[i]=-9999.;
637      chi2=0.;      chi2=0.;
638      nstep=0;      nstep=0;
639      for(int i=0; i<6; i++) xv[i]=0.;  //     for(int i=0; i<6; i++) xv[i]=0.;
640      for(int i=0; i<6; i++) yv[i]=0.;  //     for(int i=0; i<6; i++) yv[i]=0.;
641      for(int i=0; i<6; i++) zv[i]=0.;  //     for(int i=0; i<6; i++) zv[i]=0.;
642      for(int i=0; i<6; i++) axv[i]=0.;  //     for(int i=0; i<6; i++) axv[i]=0.;
643      for(int i=0; i<6; i++) ayv[i]=0.;  //     for(int i=0; i<6; i++) ayv[i]=0.;
644      for(int i=0; i<5; i++) {      for(int i=0; i<5; i++) {
645          for(int j=0; j<5; j++) coval[i][j]=0.;          for(int j=0; j<5; j++) coval[i][j]=0.;
646      }      }
647  }  }
648  /*  /**
649   * Set the tracking mode   * Set the tracking mode
650   */   */
651  void TrkTrack::SetTrackingMode(int trackmode){  void TrkTrack::SetTrackingMode(int trackmode){
652      extern cMini2track track_;      extern cMini2track track_;
653      track_.trackmode = trackmode;      track_.trackmode = trackmode;
654  }  }
655  /*  /**
656   * Set the factor scale for tracking precision   * Set the factor scale for tracking precision
657   */   */
658  void TrkTrack::SetPrecisionFactor(double fact){  void TrkTrack::SetPrecisionFactor(double fact){
659      extern cMini2track track_;      extern cMini2track track_;
660      track_.fact = fact;      track_.fact = fact;
661  }  }
662  /*  /**
663   * Set the factor scale for tracking precision   * Set the factor scale for tracking precision
664   */   */
665  void TrkTrack::SetStepMin(int istepmin){  void TrkTrack::SetStepMin(int istepmin){
# Line 588  void TrkTrack::SetStepMin(int istepmin){ Line 668  void TrkTrack::SetStepMin(int istepmin){
668  }  }
669    
670    
671  /*  /**
672   * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.   * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.
673   * If no cluster is associated, ID=-1.   * If no cluster is associated, ID=-1.
674   * @param ip Tracker plane (0-5)   * @param ip Tracker plane (0-5)
# Line 596  void TrkTrack::SetStepMin(int istepmin){ Line 676  void TrkTrack::SetStepMin(int istepmin){
676  Int_t TrkTrack::GetClusterX_ID(int ip){  Int_t TrkTrack::GetClusterX_ID(int ip){
677      return ((Int_t)fabs(xgood[ip]))%10000000-1;      return ((Int_t)fabs(xgood[ip]))%10000000-1;
678  };  };
679  /*  /**
680   * Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track.   * Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track.
681   * If no cluster is associated, ID=-1.   * If no cluster is associated, ID=-1.
682   * @param ip Tracker plane (0-5)   * @param ip Tracker plane (0-5)
# Line 604  Int_t TrkTrack::GetClusterX_ID(int ip){ Line 684  Int_t TrkTrack::GetClusterX_ID(int ip){
684  Int_t TrkTrack::GetClusterY_ID(int ip){  Int_t TrkTrack::GetClusterY_ID(int ip){
685      return ((Int_t)fabs(ygood[ip]))%10000000-1;      return ((Int_t)fabs(ygood[ip]))%10000000-1;
686  };  };
687  /*  /**
688   * Method to retrieve the ladder (0-4, increasing x) traversed by the track on this plane.   * Method to retrieve the ladder (0-4, increasing x) traversed by the track on this plane.
689   * If no ladder is traversed (dead area) the metod retuns -1.   * If no ladder is traversed (dead area) the metod retuns -1.
690   * @param ip Tracker plane (0-5)   * @param ip Tracker plane (0-5)
# Line 614  Int_t TrkTrack::GetLadder(int ip){ Line 694  Int_t TrkTrack::GetLadder(int ip){
694      if(YGood(ip))return (Int_t)fabs(ygood[ip]/100000000)-1;      if(YGood(ip))return (Int_t)fabs(ygood[ip]/100000000)-1;
695      return -1;      return -1;
696  };  };
697  /*  /**
698   * Method to retrieve the sensor (0-1, increasing y) traversed by the track on this plane.   * Method to retrieve the sensor (0-1, increasing y) traversed by the track on this plane.
699   * If no sensor is traversed (dead area) the metod retuns -1.   * If no sensor is traversed (dead area) the metod retuns -1.
700   * @param ip Tracker plane (0-5)   * @param ip Tracker plane (0-5)
# Line 625  Int_t TrkTrack::GetSensor(int ip){ Line 705  Int_t TrkTrack::GetSensor(int ip){
705      return -1;      return -1;
706  };  };
707    
708    /**
709     * \brief Method to include a x-cluster to the track.
710     * @param ip Tracker plane (0-5)
711     * @param clid Cluster ID (0,1,...)
712     * @param is Sensor (0-1, increasing y)
713     * @see Fit(double pfixed, int& fail, int iprint, int froml1)
714     */
715    void TrkTrack::SetXGood(int ip, int clid, int is){
716        int il=0;       //ladder (temporary)
717        bool bad=false; //ladder (temporary)
718        xgood[ip]=il*100000000+is*10000000+clid;
719        if(bad)xgood[ip]=-xgood[ip];
720    };
721    /**
722     * \brief Method to include a y-cluster to the track.
723     * @param ip Tracker plane (0-5)
724     * @param clid Cluster ID (0,1,...)
725     * @param is Sensor (0-1)
726     * @see Fit(double pfixed, int& fail, int iprint, int froml1)
727     */
728    void TrkTrack::SetYGood(int ip, int clid, int is){
729        int il=0;       //ladder (temporary)
730        bool bad=false; //ladder (temporary)
731        ygood[ip]=il*100000000+is*10000000+clid;
732        if(bad)ygood[ip]=-ygood[ip];
733    };
734    
735  //--------------------------------------  //--------------------------------------
736  //  //
# Line 1245  void TrkLevel2::LoadField(TString path){ Line 1351  void TrkLevel2::LoadField(TString path){
1351  //     path_.error   = 0;  //     path_.error   = 0;
1352  //     readb_();  //     readb_();
1353    
1354        TrkParams::SetTrackingMode();
1355        TrkParams::SetPrecisionFactor();
1356        TrkParams::SetStepMin();
1357    
1358      TrkParams::Set(path,1);      TrkParams::Set(path,1);
1359      TrkParams::Load(1);      TrkParams::Load(1);
1360    
1361  //  //
1362  };  };
1363  /**  // /**
1364   * Get BY (kGauss)  //  * Get BY (kGauss)
1365   * @param v (x,y,z) coordinates in cm  //  * @param v (x,y,z) coordinates in cm
1366   */  //  */
1367  float TrkLevel2::GetBX(float* v){  // float TrkLevel2::GetBX(float* v){
1368      float b[3];  //     float b[3];
1369      gufld_(v,b);  //     gufld_(v,b);
1370      return b[0]/10.;  //     return b[0]/10.;
1371  }  // }
1372  /**  // /**
1373   * Get BY (kGauss)  //  * Get BY (kGauss)
1374   * @param v (x,y,z) coordinates in cm  //  * @param v (x,y,z) coordinates in cm
1375   */  //  */
1376  float TrkLevel2::GetBY(float* v){  // float TrkLevel2::GetBY(float* v){
1377      float b[3];  //     float b[3];
1378      gufld_(v,b);  //     gufld_(v,b);
1379      return b[1]/10.;  //     return b[1]/10.;
1380  }  // }
1381  /**  // /**
1382   * Get BY (kGauss)  //  * Get BY (kGauss)
1383   * @param v (x,y,z) coordinates in cm  //  * @param v (x,y,z) coordinates in cm
1384   */  //  */
1385  float TrkLevel2::GetBZ(float* v){  // float TrkLevel2::GetBZ(float* v){
1386      float b[3];  //     float b[3];
1387      gufld_(v,b);  //     gufld_(v,b);
1388      return b[2]/10.;  //     return b[2]/10.;
1389  }  // }
1390  //--------------------------------------  //--------------------------------------
1391  //  //
1392  //  //
Legend:
Removed from v.1.32  
changed lines
  Added in v.1.33
  ViewVC Help
Powered by ViewVC 1.1.23