/[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.16 by mocchiut, Mon Apr 30 15:46:30 2007 UTC revision 1.19 by mocchiut, Mon Mar 3 09:51:02 2008 UTC
# Line 15  ToFPMT::ToFPMT(){ Line 15  ToFPMT::ToFPMT(){
15    pmt_id = 0;    pmt_id = 0;
16    adc = 0.;    adc = 0.;
17    tdc_tw = 0.;    tdc_tw = 0.;
18      tdc = 0.;
19  }  }
20    
21  ToFPMT::ToFPMT(const ToFPMT &t){  ToFPMT::ToFPMT(const ToFPMT &t){
22    pmt_id = t.pmt_id;    pmt_id = t.pmt_id;
23    adc = t.adc;    adc = t.adc;
24    tdc_tw = t.tdc_tw;    tdc_tw = t.tdc_tw;
25      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.;
32      tdc = 0.;
33  }  }
34    
35    
# Line 50  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 104  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 113  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 1046  Int_t ToFLevel2::GetNPaddle(Int_t plane) Line 1049  Int_t ToFLevel2::GetNPaddle(Int_t plane)
1049    
1050  }  }
1051    
 ////////////////////////////////////////////////////  
1052    
1053    
1054    /// wm feb 08
1055    
1056    /**
1057     * Method to calculate Beta from the 12 single measurements
1058     * we check the individual weights for artificial TDC values, then calculate
1059     * am mean beta for the first time. In a second step we loop again through
1060     * the single measurements, checking for the residual from the mean
1061     * The cut on the residual reject measurements > "x"-sigma. A chi2 value is
1062     * calculated, furthermore a "quality" value by adding the weights which
1063     * are finally used. If all measurements are taken, "quality" will be = 22.47.
1064     * A chi2 cut around 3-4 and a quality-cut > 20 is needed for clean beta
1065     * measurements like antiprotons etc.
1066     * The Level2 output is derived in the fortran routines using: 10.,10.,20.
1067     * @param notrack Track Number
1068     * @param cut on residual: difference between single measurement and mean
1069     * @param cut on "quality"
1070     * @param cut on chi2
1071     */
1072    
1073    Float_t ToFLevel2::CalcBeta(Int_t notrack, Float_t resmax, Float_t qualitycut, Float_t chi2cut){
1074    
1075    //  cout<<" in CalcBeta "<<resmax<<" "<<chi2cut<<" "<<qualitycut<<endl;
1076    
1077      Float_t bxx = 100.;
1078      //
1079      ToFTrkVar *trk = GetToFTrkVar(notrack);
1080      if(!trk) return 0; //ELENA
1081    
1082    
1083      Float_t chi2,xhelp,beta_mean;
1084      Float_t w_i[12],quality,sw,sxw,res,betachi,beta_mean_inv;
1085      Float_t b[12],tdcfl;
1086      Int_t  pmt_id,pmt_plane;
1087    
1088      for (Int_t i=0; i<12; i++){
1089        b[i] = trk->beta[i];
1090                                  }
1091          
1092    
1093    //========================================================================
1094    //---  Find out ToF layers with artificial TDC values & fill vector    ---
1095    //========================================================================
1096    
1097    Float_t  w_il[6];
1098    
1099         for (Int_t jj=0; jj<6;jj++) {
1100             w_il[jj] = 1000.;
1101                                     }
1102    
1103    
1104      for (Int_t i=0; i<trk->npmttdc; i++){
1105        //
1106        pmt_id = (trk->pmttdc).At(i);
1107        pmt_plane = GetPlaneIndex(pmt_id);
1108        tdcfl = (trk->tdcflag).At(i);
1109        if (w_il[pmt_plane] != 1.) w_il[pmt_plane] = tdcfl; //tdcflag
1110                                         };
1111      
1112    //========================================================================
1113    //---  Set weights for the 12 measurements using information for top and bottom:
1114    //---  if no measurements: weight = set to very high value=> not used
1115    //---  top or bottom artificial: weight*sqrt(2)
1116    //---  top and bottom artificial: weight*sqrt(2)*sqrt(2)
1117    //========================================================================
1118    
1119    Int_t itop[12] = {0,0,1,1,2,2,3,3,0,0,1,1};
1120    Int_t ibot[12] = {4,5,4,5,4,5,4,5,2,3,2,3};
1121    
1122         xhelp= 1E09;
1123      
1124         for (Int_t jj=0; jj<12;jj++) {
1125         if (jj<4)           xhelp = 0.11;    // S1-S3
1126         if ((jj>3)&&(jj<8)) xhelp = 0.18;    // S2-S3
1127         if (jj>7)           xhelp = 0.28;    // S1-S2
1128         if ((w_il[itop[jj]] == 1000.) && (w_il[ibot[jj]] == 1000.)) xhelp = 1E09;
1129         if ((w_il[itop[jj]] == 1) || (w_il[ibot[jj]] == 1.)) xhelp = xhelp*1.414 ;
1130         if ((w_il[itop[jj]] == 1) && (w_il[ibot[jj]] == 1.)) xhelp = xhelp*2. ;
1131    
1132         w_i[jj] = 1./xhelp;
1133                                      }
1134    
1135    
1136    //========================================================================
1137    //--- Calculate mean beta for the first time -----------------------------
1138    //--- We are using "1/beta" since its error is gaussian ------------------
1139    //========================================================================
1140    
1141          Int_t icount=0;
1142          sw=0.;
1143          sxw=0.;
1144          beta_mean=100.;
1145    
1146              for (Int_t jj=0; jj<12;jj++){
1147            if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.))
1148             {
1149                icount= icount+1;
1150                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1151                sw =sw + w_i[jj]*w_i[jj] ;
1152    
1153             }
1154             }
1155    
1156          if (icount>0) beta_mean=1./(sxw/sw);
1157          beta_mean_inv = 1./beta_mean;
1158    
1159    //========================================================================
1160    //--- Calculate beta for the second time, use residuals of the single
1161    //--- measurements to get a chi2 value
1162    //========================================================================
1163    
1164          icount=0;
1165          sw=0.;
1166          sxw=0.;
1167          betachi = 100.;
1168          chi2 = 0.;
1169          quality=0.;
1170    
1171    
1172              for (Int_t jj=0; jj<12;jj++){
1173           if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)&&(w_i[jj]>0.01)) {
1174                res = beta_mean_inv - (1./b[jj]) ;
1175                if (fabs(res*w_i[jj])<resmax)          {;
1176                chi2 = chi2 + pow((res*w_i[jj]),2) ;
1177                icount= icount+1;
1178                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1179                sw =sw + w_i[jj]*w_i[jj] ;
1180                                                   }
1181                                                                            }
1182                                          }
1183          quality = sqrt(sw) ;
1184    
1185          if (icount==0) chi2 = 1000.;
1186          if (icount>0) chi2 = chi2/(icount) ;
1187          if (icount>0) betachi=1./(sxw/sw);
1188    
1189       bxx = 100.;
1190       if ((chi2 < chi2cut)&&(quality>qualitycut)) bxx = betachi;
1191      //
1192      return(bxx);
1193    };
1194    
1195    
1196    ////////////////////////////////////////////////////
1197    ////////////////////////////////////////////////////
1198    
1199    
1200  /**  /**
1201   * Fills a struct cToFLevel2 with values from a ToFLevel2 object (to put data into a F77 common).   * Fills a struct cToFLevel2 with values from a ToFLevel2 object (to put data into a F77 common).
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.19
  ViewVC Help
Powered by ViewVC 1.1.23