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

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

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

revision 1.23 by pam-fi, Fri Aug 31 14:56:51 2007 UTC revision 1.29 by mocchiut, Thu Jul 3 10:38:00 2014 UTC
# Line 11  using namespace std; Line 11  using namespace std;
11  extern "C" {  extern "C" {
12                    
13  //      int readetaparam_();  //      int readetaparam_();
14      float cog_(int*,int*);    float cog_(int*,int*);
15      float pfaeta_(int*,float*);    float pfaeta_(int*,float*);
16      float pfaeta2_(int*,float*);    float pfaeta2_(int*,float*);
17      float pfaeta3_(int*,float*);    float pfaeta3_(int*,float*);
18      float pfaeta4_(int*,float*);    float pfaeta4_(int*,float*);
19      float pfaetal_(int*,float*);    float pfaetal_(int*,float*);
20      int   npfastrips_(int*,float*);    float digsat_(int*);
21              int   npfastrips_(int*,float*);
22      
23      float fbad_cog_(int*,int*);
24      float risx_cog_(float*);
25      float risy_cog_(float*);
26  }  }
27  //--------------------------------------  //--------------------------------------
28  //  //
# Line 407  void TrkCluster::GetLevel1Struct(cTrkLev Line 411  void TrkCluster::GetLevel1Struct(cTrkLev
411   *      @param ncog Number of strips to evaluate COG.     *      @param ncog Number of strips to evaluate COG.  
412   * If ncog=0, the COG of the cluster is evaluated according to the cluster multiplicity (defined by the inclusion cut).   * If ncog=0, the COG of the cluster is evaluated according to the cluster multiplicity (defined by the inclusion cut).
413   * If ncog>0, the COG is evaluated using ncog strips, even if they have a negative signal (according to G.Landi)   * If ncog>0, the COG is evaluated using ncog strips, even if they have a negative signal (according to G.Landi)
414     *
415     * (NB TrkCluster::GetLevel1Struct() showld be called first, in order to fill the F77 level1 common with this single cluster)
416   */   */
417  Float_t TrkCluster::GetCOG(Int_t ncog){  Float_t TrkCluster::GetCOG(Int_t ncog){
418                    
419      int ic = 1;      int ic = 1;
420      GetLevel1Struct();      //    GetLevel1Struct(); //Elena: dangerous...
421      return cog_(&ncog,&ic);      return cog_(&ncog,&ic);
422                    
423  };  };
# Line 455  Float_t TrkCluster::GetCOG(Float_t angle Line 461  Float_t TrkCluster::GetCOG(Float_t angle
461   *  @param angle Projected (effective) angle between particle track and detector plane.   *  @param angle Projected (effective) angle between particle track and detector plane.
462   *  @landi flag to apply Landi correction   *  @landi flag to apply Landi correction
463   * Implemented values of neta are 2,3,4. If neta=0, ETA2, ETA3 and ETA4 are applied according to the angle.   * Implemented values of neta are 2,3,4. If neta=0, ETA2, ETA3 and ETA4 are applied according to the angle.
464     * (NB TrkCluster::GetLevel1Struct() showld be called first, in order to fill the F77 level1 common with this single cluster)
465   */   */
466  Float_t TrkCluster::GetETA(Int_t neta, float angle, bool landi){  Float_t TrkCluster::GetETA(Int_t neta, float angle, bool landi){
467                    
# Line 469  Float_t TrkCluster::GetETA(Int_t neta, f Line 476  Float_t TrkCluster::GetETA(Int_t neta, f
476    
477      float ax = angle;      float ax = angle;
478      int ic = 1;      int ic = 1;
479      GetLevel1Struct();      //GetLevel1Struct(); //Elena: dangerous...
480      if(     neta == 0 && !landi) return pfaeta_(&ic,&ax);      if(     neta == 0 && !landi) return pfaeta_(&ic,&ax);
481      else if(neta == 0 && landi ) return pfaetal_(&ic,&ax);      else if(neta == 0 && landi ) return pfaetal_(&ic,&ax);
482      else if(neta == 2          ) return pfaeta2_(&ic,&ax);      else if(neta == 2          ) return pfaeta2_(&ic,&ax);
# Line 481  Float_t TrkCluster::GetETA(Int_t neta, f Line 488  Float_t TrkCluster::GetETA(Int_t neta, f
488  };  };
489    
490  /**  /**
491     * Evaluates the cluster position, in pitch units, relative to the strip
492     *  with the maximum signal (TrkCluster::maxs), by applying the digital
493     *  algorithm for saturated clusters.
494     *
495     *  @return The cluster position (0 also if if no saturated strip is found).
496     *
497     * (NB TrkCluster::GetLevel1Struct() showld be called first, in order to fill the F77 level1 common with this single cluster)
498     */
499    Float_t TrkCluster::GetDigSat() {
500    
501      //  GetLevel1Struct(); //Elena: dangerous...
502      int ic = 1;
503      return digsat_(&ic);
504    
505    }
506    
507    /**
508   * Evaluates the cluster position, in pitch unit, relative to the strip with   * Evaluates the cluster position, in pitch unit, relative to the strip with
509   * the maximum signal (TrkCluster::maxs), by applying the PFA set as default (see TrkParams).   * the maximum signal (TrkCluster::maxs), by applying the PFA set as default (see TrkParams).
510   *  @param angle Projected (effective) angle between particle track and detector plane.   *  @param angle Projected (effective) angle between particle track and detector plane.
# Line 508  Float_t TrkCluster::GetPositionPU(float Line 532  Float_t TrkCluster::GetPositionPU(float
532   * according to the p.f.a.   * according to the p.f.a.
533   * It returns 0 when the COG is used (in this case the number of strip used   * It returns 0 when the COG is used (in this case the number of strip used
534   * equals the multiplicity).   * equals the multiplicity).
535     * (NB TrkCluster::GetLevel1Struct() showld be called first, in order to fill the F77 level1 common with this single cluster)
536   */   */
537  Int_t TrkCluster::GetPFAstrips(float angle){  Int_t TrkCluster::GetPFAstrips(float angle){
538    
539      float ax = angle;      float ax = angle;
540      int ic = 1;      int ic = 1;
541      GetLevel1Struct();      //    GetLevel1Struct(); //Elena: dangerous...
542      return npfastrips_(&ic,&ax);      return npfastrips_(&ic,&ax);
543    
544  }  }
# Line 534  TrkLevel1::TrkLevel1(){ Line 559  TrkLevel1::TrkLevel1(){
559              cnn[j][i]=0;              cnn[j][i]=0;
560          };          };
561      };      };
562      TrkParams::SetTrackingMode();  //     TrkParams::SetTrackingMode();
563      TrkParams::SetPrecisionFactor();  //     TrkParams::SetPrecisionFactor();
564      TrkParams::SetStepMin();  //     TrkParams::SetStepMin();
565        TrkParams::SetMiniDefault();
566      TrkParams::SetPFA();      TrkParams::SetPFA();
567  }  }
568  //--------------------------------------  //--------------------------------------
# Line 761  TrkCluster *TrkLevel1::GetCluster(int is Line 787  TrkCluster *TrkLevel1::GetCluster(int is
787      TrkCluster *cluster = (TrkCluster*)t[is];      TrkCluster *cluster = (TrkCluster*)t[is];
788      return cluster;      return cluster;
789  }  }
 //--------------------------------------  
 //  
 //  
 //--------------------------------------  
 // /**  
 //  * Load Position-Finding-Algorythm parameters (call the F77 routine).  
 //  *  
 //  */  
 // int TrkLevel1::LoadPfaParam(TString path){  
           
 //     if( path.IsNull() ){  
 //      path = gSystem->Getenv("PAM_CALIB");  
 //      if(path.IsNull()){  
 //          cout << " TrkLevel1::LoadPfaParam() ==> No PAMELA environment variables defined "<<endl;  
 //          return 0;  
 //      }  
 //      path.Append("/trk-param/eta_param-0/");  
 //     }  
790    
 //     strcpy(path_.path,path.Data());  
 //     path_.pathlen = path.Length();  
 //     path_.error   = 0;  
 //     cout <<"Loading p.f.a. parameters: "<<path<<endl;  
 //     return readetaparam_();  
 // }  
791    
792  // /**  // int TrkLevel1::GetPfaNbinsAngle(){
793  //  * Load magnetic field parameters (call the F77 routine).  //     TrkParams::Load(4);
794  //  *  //     if( !TrkParams::IsLoaded(4) ){
795  //  */  //      cout << "int TrkLevel1::GetPfaNbinsAngle() --- ERROR --- p.f.a. parameters  not loaded"<<endl;
796  // int TrkLevel1::LoadFieldParam(TString path){  //      return 0;
           
 // //    if( strcmp(path_.path,path.Data()) ){  
 //     if( path.IsNull() ){  
 //      path = gSystem->Getenv("PAM_CALIB");  
 //      if(path.IsNull()){  
 //          cout << " TrkLevel1::LoadFieldParam() ==> No PAMELA environment variables defined "<<endl;  
 //          return 0;  
 //      }  
 //      path.Append("/trk-param/field_param-0/");  
797  //     }  //     }
798  //     cout <<"Loading magnetic field "<<path<<endl;  //     return pfa_.nangbin;
799  //     strcpy(path_.path,path.Data());  // };
800  //     path_.pathlen = path.Length();  
801  //     path_.error   = 0;  // int TrkLevel1::GetPfaNbinsETA(){
802  //     return readb_();  //     TrkParams::Load(4);
803  // //    }        //     if( !TrkParams::IsLoaded(4) ){
804  // //    return 0;  //      cout << "int TrkLevel1::GetPfaNbinsETA() --- ERROR --- p.f.a. parameters  not loaded"<<endl;
805  // }  //      return 0;
 // /**  
 //  * Load magnetic field parameters (call the F77 routine).  
 //  *  
 //  */  
 // int TrkLevel1::LoadChargeParam(TString path){  
           
 // //    if( strcmp(path_.path,path.Data()) ){  
 //     if( path.IsNull() ){  
 //      path = gSystem->Getenv("PAM_CALIB");  
 //      if(path.IsNull()){  
 //          cout << " TrkLevel1::LoadChargeParam() ==> No PAMELA environment variables defined "<<endl;  
 //          return 0;  
 //      }  
 //      path.Append("/trk-param/charge_param-1/");  
 //     }  
 //     cout <<"Loading charge-correlation parameters: "<<path<<endl;  
 //     strcpy(path_.path,path.Data());  
 //     path_.pathlen = path.Length();  
 //     path_.error   = 0;  
 //     return readchargeparam_();  
 // //    }        
 // //    return 0;  
 // }  
 // /**  
 //  * Load magnetic field parameters (call the F77 routine).  
 //  *  
 //  */  
 // int TrkLevel1::LoadAlignmentParam(TString path){  
           
 // //    if( strcmp(path_.path,path.Data()) ){  
 //     if( path.IsNull() ){  
 //      path = gSystem->Getenv("PAM_CALIB");  
 //      if(path.IsNull()){  
 //          cout << " TrkLevel1::LoadAlignmentParam() ==> No PAMELA environment variables defined "<<endl;  
 //          return 0;  
 //      }  
 //      path.Append("/trk-param/align_param-0/");  
 //     }  
 //     cout <<"Loading alignment parameters: "<<path<<endl;  
 //     strcpy(path_.path,path.Data());  
 //     path_.pathlen = path.Length();  
 //     path_.error   = 0;  
 //     return readalignparam_();  
 // //    }        
 // //    return 0;  
 // }  
 // /**  
 //  * Load magnetic field parameters (call the F77 routine).  
 //  *  
 //  */  
 // int TrkLevel1::LoadMipParam(TString path){  
           
 // //    if( strcmp(path_.path,path.Data()) ){  
 //     if( path.IsNull() ){  
 //      path = gSystem->Getenv("PAM_CALIB");  
 //      if(path.IsNull()){  
 //          cout << " TrkLevel1::LoadMipParam() ==> No PAMELA environment variables defined "<<endl;  
 //          return 0;  
 //      }  
 //      path.Append("/trk-param/mip_param-0/");  
 //     }  
 //     cout <<"Loading ADC-to-MIP conversion parameters: "<<path<<endl;  
 //     strcpy(path_.path,path.Data());  
 //     path_.pathlen = path.Length();  
 //     path_.error   = 0;  
 //     return readmipparam_();  
 // //    }        
 // //    return 0;  
 // }  
 // /**  
 //  * Load magnetic field parameters (call the F77 routine).  
 //  *  
 //  */  
 // int TrkLevel1::LoadVKMaskParam(TString path){  
           
 // //    if( strcmp(path_.path,path.Data()) ){  
 //     if( path.IsNull() ){  
 //      path = gSystem->Getenv("PAM_CALIB");  
 //      if(path.IsNull()){  
 //          cout << " TrkLevel1::LoadVKMaskParam() ==> No PAMELA environment variables defined "<<endl;  
 //          return 0;  
 //      }  
 //      path.Append("/trk-param/mask_param-1/");  
806  //     }  //     }
807  //     cout <<"Loading VK-mask parameters: "<<path<<endl;  //     return pfa_.netaval;
808  //     strcpy(path_.path,path.Data());  // };
 //     path_.pathlen = path.Length();  
 //     path_.error   = 0;  
 //     return readvkmask_();  
 // //    }        
 // //    return 0;  
 // }  
809    
810  // /**  // /**
811  //  * Load all (default) parameters. Environment variable must be defined.  //  *
812  //  *  //  *
813  //  */  //  */
814  // int TrkLevel1::LoadParams(){  // float* TrkLevel1::GetPfaCoord(TString pfa, int nview, int nladder, int nang){
815    
816  //     int result=0;  //     TrkParams::Load(4);
817        //     if( !TrkParams::IsLoaded(4) ){
818  //     result = result * LoadFieldParam();  //      cout << "float* TrkLevel1::GetPfaCoord(TString pfa, int nview, int nladder, int nang) --- ERROR --- p.f.a. parameters  not loaded"<<endl;
819  //     result = result * LoadPfaParam();  //      return 0;
820  //     result = result * LoadChargeParam();  //     }
 //     result = result * LoadAlignmentParam();  
 //     result = result * LoadMipParam();  
 //     result = result * LoadVKMaskParam();  
   
 //     return result;  
 // }  
   
   
   
 int TrkLevel1::GetPfaNbinsAngle(){  
     TrkParams::Load(4);  
     if( !TrkParams::IsLoaded(4) ){  
         cout << "int TrkLevel1::GetPfaNbinsAngle() --- ERROR --- p.f.a. parameters  not loaded"<<endl;  
         return 0;  
     }  
     return pfa_.nangbin;  
 };  
   
 int TrkLevel1::GetPfaNbinsETA(){  
     TrkParams::Load(4);  
     if( !TrkParams::IsLoaded(4) ){  
         cout << "int TrkLevel1::GetPfaNbinsETA() --- ERROR --- p.f.a. parameters  not loaded"<<endl;  
         return 0;  
     }  
     return pfa_.netaval;  
 };  
   
 /**  
  *  
  *  
  */  
 float* TrkLevel1::GetPfaCoord(TString pfa, int nview, int nladder, int nang){  
   
     TrkParams::Load(4);  
     if( !TrkParams::IsLoaded(4) ){  
         cout << "float* TrkLevel1::GetPfaCoord(TString pfa, int nview, int nladder, int nang) --- ERROR --- p.f.a. parameters  not loaded"<<endl;  
         return 0;  
     }  
821        
822      int nbins = GetPfaNbinsETA();  //     int nbins = GetPfaNbinsETA();
823      if(!nbins)return 0;  //     if(!nbins)return 0;
824    
825      float *fcorr = new float [nbins];  //     float *fcorr = new float [nbins];
826    
827      if(!pfa.CompareTo("ETA2",TString::kIgnoreCase)){  //     if(!pfa.CompareTo("ETA2",TString::kIgnoreCase)){
828          for(int ib=0; ib<nbins; ib++){  //      for(int ib=0; ib<nbins; ib++){
829              fcorr[ib] = pfa_.feta2[nang][nladder][nview][ib];  //          fcorr[ib] = pfa_.feta2[nang][nladder][nview][ib];
830              cout << pfa_.eta2[nang][ib] << " - " <<  pfa_.feta2[nang][nladder][nview][ib]<<endl;;  //          cout << pfa_.eta2[nang][ib] << " - " <<  pfa_.feta2[nang][nladder][nview][ib]<<endl;;
831          }  //      }
832      }else if (!pfa.CompareTo("ETA3",TString::kIgnoreCase)){  //     }else if (!pfa.CompareTo("ETA3",TString::kIgnoreCase)){
833          for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.feta3[nang][nladder][nview][ib];  //      for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.feta3[nang][nladder][nview][ib];
834      }else if (!pfa.CompareTo("ETA4",TString::kIgnoreCase)){  //     }else if (!pfa.CompareTo("ETA4",TString::kIgnoreCase)){
835          for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.feta4[nang][nladder][nview][ib];  //      for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.feta4[nang][nladder][nview][ib];
836      }else{  //     }else{
837          cout << pfa<<" pfa parameters not implemented "<<endl;  //      cout << pfa<<" pfa parameters not implemented "<<endl;
838          return 0;  //      return 0;
839      }      //     }    
840    
841      return fcorr;  //     return fcorr;
842    
843  };  // };
844    
845  float* TrkLevel1::GetPfaAbs(TString pfa, int nang){  // float* TrkLevel1::GetPfaAbs(TString pfa, int nang){
846        
847      TrkParams::Load(4);  //     TrkParams::Load(4);
848      if( !TrkParams::IsLoaded(4) ){  //     if( !TrkParams::IsLoaded(4) ){
849          cout << "float* TrkLevel1::GetPfaAbs(TString pfa, int nang) --- ERROR --- p.f.a. parameters  not loaded"<<endl;  //      cout << "float* TrkLevel1::GetPfaAbs(TString pfa, int nang) --- ERROR --- p.f.a. parameters  not loaded"<<endl;
850          return 0;  //      return 0;
851      }  //     }
852    
853      int nbins = GetPfaNbinsETA();  //     int nbins = GetPfaNbinsETA();
854      if(!nbins)return 0;  //     if(!nbins)return 0;
855    
856      float *fcorr = new float [nbins];  //     float *fcorr = new float [nbins];
857    
858      if(!pfa.CompareTo("ETA2",TString::kIgnoreCase)){  //     if(!pfa.CompareTo("ETA2",TString::kIgnoreCase)){
859          for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta2[nang][ib];  //      for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta2[nang][ib];
860      }else if (!pfa.CompareTo("ETA3",TString::kIgnoreCase)){  //     }else if (!pfa.CompareTo("ETA3",TString::kIgnoreCase)){
861          for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta3[nang][ib];  //      for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta3[nang][ib];
862      }else if (!pfa.CompareTo("ETA4",TString::kIgnoreCase)){  //     }else if (!pfa.CompareTo("ETA4",TString::kIgnoreCase)){
863          for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta4[nang][ib];  //      for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta4[nang][ib];
864      }else{  //     }else{
865          cout << pfa<<" pfa parameters not implemented "<<endl;  //      cout << pfa<<" pfa parameters not implemented "<<endl;
866          return 0;  //      return 0;
867      }      //     }    
868    
869      return fcorr;  //     return fcorr;
870    
871  };  // };
872    
873  /**  /**
874   * Method to call the F77 routine that performs level1->level2 processing.   * Method to call the F77 routine that performs level1->level2 processing.
# Line 1031  int TrkLevel1::ProcessEvent(){ Line 897  int TrkLevel1::ProcessEvent(){
897    
898  }  }
899    
900    //--------------------------------------
901    //
902    //
903    //--------------------------------------
904    /**
905     * Method to fill a TrkLevel1 object from an existing one, by cleaning low-signal clusters.
906     *
907     */
908    void TrkLevel1::Set(TrkLevel1 *trkl1, float mipCut, float fCut){
909    
910    
911        
912        if(!trkl1)return;
913    
914        //  -------------------------
915        //  ****general variables****
916        //  -------------------------    
917        for(Int_t i=0; i<12 ; i++){
918            good[i] = trkl1->good[i];
919            for(Int_t j=0; j<24 ; j++){
920                cn[j][i]     = trkl1->cn[j][i];
921                cnn[j][i]    = trkl1->cnn[j][i];
922            };
923        };
924        //  -------------------------
925        //  ****cluster array****
926        //  -------------------------    
927    
928        if(Cluster)Cluster->Clear("C");
929        Cluster = new TClonesArray("TrkCluster");
930        TClonesArray &t = *Cluster;
931    
932        if ( trkl1->GetClusters() ){ // EM: crash nuclei alg when  DSP status: -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
933          int isel=0;
934          for(int icl=0 ; icl< trkl1->GetClusters()->GetEntries(); icl++){
935            TrkCluster *cl = trkl1->GetCluster(icl);
936            
937            float mip = TrkParams::GetMIP(cl->GetLadder()-1,cl->view-1);
938            float smip = cl->GetSignal()/(mip>0.?mip:1.);
939            float smax =  cl->clsignal[cl->indmax]/(mip>0.?mip:1.);
940            if(smax/smip<fCut)continue;
941            if(smip<mipCut)continue;
942            if(smax<0.5*mipCut)continue;
943            
944            
945            
946            new(t[isel]) TrkCluster(*cl); // <<< store cluster
947            isel++;
948          }
949        }
950    
951    
952    
953    }
954    
955  ClassImp(TrkLevel1);  ClassImp(TrkLevel1);
956  ClassImp(TrkCluster);  ClassImp(TrkCluster);
Legend:
Removed from v.1.23  
changed lines
  Added in v.1.29
  ViewVC Help
Powered by ViewVC 1.1.23