/[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.18 by mocchiut, Mon Nov 26 08:01:17 2007 UTC revision 1.19 by mocchiut, Mon Mar 3 09:51:02 2008 UTC
# Line 1049  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.18  
changed lines
  Added in v.1.19

  ViewVC Help
Powered by ViewVC 1.1.23