/[PAMELA software]/DarthVader/ToFLevel2/src/ToFLevel2.cpp
ViewVC logotype

Diff of /DarthVader/ToFLevel2/src/ToFLevel2.cpp

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

revision 1.17 by mocchiut, Wed Oct 10 16:01:33 2007 UTC revision 1.20 by mocchiut, Thu Mar 6 14:51:07 2008 UTC
# Line 25  ToFPMT::ToFPMT(const ToFPMT &t){ Line 25  ToFPMT::ToFPMT(const ToFPMT &t){
25    tdc = t.tdc;    tdc = t.tdc;
26  }  }
27    
28  void ToFPMT::Clear(){  void ToFPMT::Clear(Option_t *t){
29    pmt_id = 0;    pmt_id = 0;
30    adc = 0.;    adc = 0.;
31    tdc_tw = 0.;    tdc_tw = 0.;
# Line 53  ToFTrkVar::ToFTrkVar() { Line 53  ToFTrkVar::ToFTrkVar() {
53    //    //
54  };  };
55    
56  void ToFTrkVar::Clear() {  void ToFTrkVar::Clear(Option_t *t) {
57    trkseqno = 0;    trkseqno = 0;
58    npmttdc = 0;    npmttdc = 0;
59    npmtadc = 0;    npmtadc = 0;
# Line 107  void ToFLevel2::Set(){//ELENA Line 107  void ToFLevel2::Set(){//ELENA
107      if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA      if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
108  }//ELENA  }//ELENA
109    
110  void ToFLevel2::Clear(){  void ToFLevel2::Clear(Option_t *t){
111    //    //
112    if(ToFTrk)ToFTrk->Delete(); //ELENA    if(ToFTrk)ToFTrk->Delete(); //ELENA
113    if(PMT)PMT->Delete(); //ELENA    if(PMT)PMT->Delete(); //ELENA
# Line 116  void ToFLevel2::Clear(){ Line 116  void ToFLevel2::Clear(){
116    //    //
117  };  };
118    
119  void ToFLevel2::Delete(){ //ELENA  void ToFLevel2::Delete(Option_t *t){ //ELENA
120    //    //
121    if(ToFTrk){    if(ToFTrk){
122        ToFTrk->Delete(); //ELENA        ToFTrk->Delete(); //ELENA
# Line 696  void ToFLevel2::GetPMTPaddle(Int_t pmt_i Line 696  void ToFLevel2::GetPMTPaddle(Int_t pmt_i
696  // gf Apr 07  // gf Apr 07
697    
698  void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){  void ToFLevel2::GetPaddlePMT(Int_t paddle, Int_t &pmtleft, Int_t &pmtright){
699      pmtleft=paddle*2;
700    if(paddle==0){    pmtright= pmtleft+1;
     pmtleft=0;  
     pmtright=1;  
   }  
   
   if(paddle==1){  
     pmtleft=2;  
     pmtright=3;  
   }  
   
   if(paddle==2){  
     pmtleft=4;  
     pmtright=5;  
   }  
   
   if(paddle==3){  
     pmtleft=6;  
     pmtright=7;  
   }  
   
   if(paddle==4){  
     pmtleft=8;  
     pmtright=9;  
   }  
   
   if(paddle==5){  
     pmtleft=10;  
     pmtright=11;  
   }  
   
   if(paddle==6){  
     pmtleft=12;  
     pmtright=13;  
   }  
   
   if(paddle==7){  
     pmtleft=14;  
     pmtright=15;  
   }  
   
   if(paddle==8){  
     pmtleft=16;  
     pmtright=17;  
   }  
   
   if(paddle==9){  
     pmtleft=18;  
     pmtright=19;  
   }  
   
   if(paddle==10){  
     pmtleft=20;  
     pmtright=21;  
   }  
   
   if(paddle==11){  
     pmtleft=22;  
     pmtright=23;  
   }  
   
   if(paddle==12){  
     pmtleft=24;  
     pmtright=25;  
   }  
   
   if(paddle==13){  
     pmtleft=26;  
     pmtright=27;  
   }  
   
   if(paddle==14){  
     pmtleft=28;  
     pmtright=29;  
   }  
   
   if(paddle==15){  
     pmtleft=30;  
     pmtright=31;  
   }  
   
   if(paddle==16){  
     pmtleft=32;  
     pmtright=33;  
   }  
   
   if(paddle==17){  
     pmtleft=34;  
     pmtright=35;  
   }  
   
   if(paddle==18){  
     pmtleft=36;  
     pmtright=37;  
   }  
   
   if(paddle==19){  
     pmtleft=38;  
     pmtright=39;  
   }  
   
   if(paddle==20){  
     pmtleft=40;  
     pmtright=41;  
   }  
   
   if(paddle==21){  
     pmtleft=42;  
     pmtright=43;  
   }  
   
   if(paddle==22){  
     pmtleft=44;  
     pmtright=45;  
   }  
   
   if(paddle==23){  
     pmtleft=46;  
     pmtright=47;  
   }  
     
701    return;    return;
702  }  }
703    
# Line 930  void ToFLevel2::GetPaddleGeometry(Int_t Line 811  void ToFLevel2::GetPaddleGeometry(Int_t
811   */   */
812  Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)  Int_t ToFLevel2::GetPaddleid(Int_t plane, Int_t paddle)
813  {  {
   
814    Int_t padid=-1;    Int_t padid=-1;
815    Int_t pads11=8;    Int_t pads[6]={8,6,2,2,3,3};
   Int_t pads12=6;  
   Int_t pads21=2;  
   Int_t pads22=2;  
   Int_t pads31=3;  
   //  Int_t pads32=3;  
   
816    
817    if(plane == 0){    int somma=0;
818      padid=paddle;    int np=plane;
819    }    for(Int_t j=0; j<np; j++){
820        somma+=pads[j];
   if(plane == 1){  
     padid=pads11+paddle;  
   }  
   
   if(plane == 2){  
     padid=pads11+pads12+paddle;  
821    }    }
822      padid=paddle+somma;
   if(plane == 3){  
     padid=pads11+pads12+pads21+paddle;  
   }  
   
   if(plane == 4){  
     padid=pads11+pads12+pads21+pads22+paddle;  
   }  
   
   if(plane == 5){  
     padid=pads11+pads12+pads21+pads22+pads31+paddle;  
   }  
   
823    return padid;    return padid;
   
824  }  }
825    
826    
# Line 1049  Int_t ToFLevel2::GetNPaddle(Int_t plane) Line 904  Int_t ToFLevel2::GetNPaddle(Int_t plane)
904    
905  }  }
906    
907    /// wm feb 08
908    
909    /**
910     * Method to calculate Beta from the 12 single measurements
911     * we check the individual weights for artificial TDC values, then calculate
912     * am mean beta for the first time. In a second step we loop again through
913     * the single measurements, checking for the residual from the mean
914     * The cut on the residual reject measurements > "x"-sigma. A chi2 value is
915     * calculated, furthermore a "quality" value by adding the weights which
916     * are finally used. If all measurements are taken, "quality" will be = 22.47.
917     * A chi2 cut around 3-4 and a quality-cut > 20 is needed for clean beta
918     * measurements like antiprotons etc.
919     * The Level2 output is derived in the fortran routines using: 10.,10.,20.
920     * @param notrack Track Number
921     * @param cut on residual: difference between single measurement and mean
922     * @param cut on "quality"
923     * @param cut on chi2
924     */
925    
926    Float_t ToFLevel2::CalcBeta(Int_t notrack, Float_t resmax, Float_t qualitycut, Float_t chi2cut){
927    
928    //  cout<<" in CalcBeta "<<resmax<<" "<<chi2cut<<" "<<qualitycut<<endl;
929    
930      Float_t bxx = 100.;
931      //
932      ToFTrkVar *trk = GetToFTrkVar(notrack);
933      if(!trk) return 0; //ELENA
934    
935    
936      Float_t chi2,xhelp,beta_mean;
937      Float_t w_i[12],quality,sw,sxw,res,betachi,beta_mean_inv;
938      Float_t b[12],tdcfl;
939      Int_t  pmt_id,pmt_plane;
940    
941      for (Int_t i=0; i<12; i++){
942        b[i] = trk->beta[i];
943                                  }
944          
945    
946    //========================================================================
947    //---  Find out ToF layers with artificial TDC values & fill vector    ---
948    //========================================================================
949    
950    Float_t  w_il[6];
951    
952         for (Int_t jj=0; jj<6;jj++) {
953             w_il[jj] = 1000.;
954                                     }
955    
956    
957      for (Int_t i=0; i<trk->npmttdc; i++){
958        //
959        pmt_id = (trk->pmttdc).At(i);
960        pmt_plane = GetPlaneIndex(pmt_id);
961        tdcfl = (trk->tdcflag).At(i);
962        if (w_il[pmt_plane] != 1.) w_il[pmt_plane] = tdcfl; //tdcflag
963                                         };
964      
965    //========================================================================
966    //---  Set weights for the 12 measurements using information for top and bottom:
967    //---  if no measurements: weight = set to very high value=> not used
968    //---  top or bottom artificial: weight*sqrt(2)
969    //---  top and bottom artificial: weight*sqrt(2)*sqrt(2)
970    //========================================================================
971    
972    Int_t itop[12] = {0,0,1,1,2,2,3,3,0,0,1,1};
973    Int_t ibot[12] = {4,5,4,5,4,5,4,5,2,3,2,3};
974    
975         xhelp= 1E09;
976      
977         for (Int_t jj=0; jj<12;jj++) {
978         if (jj<4)           xhelp = 0.11;    // S1-S3
979         if ((jj>3)&&(jj<8)) xhelp = 0.18;    // S2-S3
980         if (jj>7)           xhelp = 0.28;    // S1-S2
981         if ((w_il[itop[jj]] == 1000.) && (w_il[ibot[jj]] == 1000.)) xhelp = 1E09;
982         if ((w_il[itop[jj]] == 1) || (w_il[ibot[jj]] == 1.)) xhelp = xhelp*1.414 ;
983         if ((w_il[itop[jj]] == 1) && (w_il[ibot[jj]] == 1.)) xhelp = xhelp*2. ;
984    
985         w_i[jj] = 1./xhelp;
986                                      }
987    
988    
989    //========================================================================
990    //--- Calculate mean beta for the first time -----------------------------
991    //--- We are using "1/beta" since its error is gaussian ------------------
992    //========================================================================
993    
994          Int_t icount=0;
995          sw=0.;
996          sxw=0.;
997          beta_mean=100.;
998    
999              for (Int_t jj=0; jj<12;jj++){
1000            if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.))
1001             {
1002                icount= icount+1;
1003                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1004                sw =sw + w_i[jj]*w_i[jj] ;
1005    
1006             }
1007             }
1008    
1009          if (icount>0) beta_mean=1./(sxw/sw);
1010          beta_mean_inv = 1./beta_mean;
1011    
1012    //========================================================================
1013    //--- Calculate beta for the second time, use residuals of the single
1014    //--- measurements to get a chi2 value
1015    //========================================================================
1016    
1017          icount=0;
1018          sw=0.;
1019          sxw=0.;
1020          betachi = 100.;
1021          chi2 = 0.;
1022          quality=0.;
1023    
1024    
1025              for (Int_t jj=0; jj<12;jj++){
1026           if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)&&(w_i[jj]>0.01)) {
1027                res = beta_mean_inv - (1./b[jj]) ;
1028                if (fabs(res*w_i[jj])<resmax)          {;
1029                chi2 = chi2 + pow((res*w_i[jj]),2) ;
1030                icount= icount+1;
1031                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1032                sw =sw + w_i[jj]*w_i[jj] ;
1033                                                   }
1034                                                                            }
1035                                          }
1036          quality = sqrt(sw) ;
1037    
1038          if (icount==0) chi2 = 1000.;
1039          if (icount>0) chi2 = chi2/(icount) ;
1040          if (icount>0) betachi=1./(sxw/sw);
1041    
1042       bxx = 100.;
1043       if ((chi2 < chi2cut)&&(quality>qualitycut)) bxx = betachi;
1044      //
1045      return(bxx);
1046    };
1047    
1048    
1049    ////////////////////////////////////////////////////
1050  ////////////////////////////////////////////////////  ////////////////////////////////////////////////////
1051    
1052    
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.23