/[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.23 by pamelats, Thu Dec 4 11:41:08 2008 UTC
# Line 1  Line 1 
1  /**  /**
2   * \file ToFLevel2.cpp   * \file ToFLevel2.cpp
3   * \author Gianfranca DeRosa, Wolfgang Menn   * \author Gianfranca DeRosa, Wolfgang Menn
4     *
5     * WM dec 2008: Description of "GetdEdx" changed
6     * WM dec 2008: "GetdEdxPaddle" modified: Now includes saturation limit
7     *              PMTs higher than the saturation limit are not used for dEdx
8   */   */
9    
10  #include <TObject.h>  #include <TObject.h>
# Line 15  ToFPMT::ToFPMT(){ Line 19  ToFPMT::ToFPMT(){
19    pmt_id = 0;    pmt_id = 0;
20    adc = 0.;    adc = 0.;
21    tdc_tw = 0.;    tdc_tw = 0.;
22      tdc = 0.;
23  }  }
24    
25  ToFPMT::ToFPMT(const ToFPMT &t){  ToFPMT::ToFPMT(const ToFPMT &t){
26    pmt_id = t.pmt_id;    pmt_id = t.pmt_id;
27    adc = t.adc;    adc = t.adc;
28    tdc_tw = t.tdc_tw;    tdc_tw = t.tdc_tw;
29      tdc = t.tdc;
30  }  }
31    
32  void ToFPMT::Clear(){  void ToFPMT::Clear(Option_t *t){
33    pmt_id = 0;    pmt_id = 0;
34    adc = 0.;    adc = 0.;
35    tdc_tw = 0.;    tdc_tw = 0.;
36      tdc = 0.;
37  }  }
38    
39    
# Line 50  ToFTrkVar::ToFTrkVar() { Line 57  ToFTrkVar::ToFTrkVar() {
57    //    //
58  };  };
59    
60  void ToFTrkVar::Clear() {  void ToFTrkVar::Clear(Option_t *t) {
61    trkseqno = 0;    trkseqno = 0;
62    npmttdc = 0;    npmttdc = 0;
63    npmtadc = 0;    npmtadc = 0;
# Line 104  void ToFLevel2::Set(){//ELENA Line 111  void ToFLevel2::Set(){//ELENA
111      if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA      if(!ToFTrk)ToFTrk = new TClonesArray("ToFTrkVar",2); //ELENA
112  }//ELENA  }//ELENA
113    
114  void ToFLevel2::Clear(){  void ToFLevel2::Clear(Option_t *t){
115    //    //
116    if(ToFTrk)ToFTrk->Delete(); //ELENA    if(ToFTrk)ToFTrk->Delete(); //ELENA
117    if(PMT)PMT->Delete(); //ELENA    if(PMT)PMT->Delete(); //ELENA
# Line 113  void ToFLevel2::Clear(){ Line 120  void ToFLevel2::Clear(){
120    //    //
121  };  };
122    
123  void ToFLevel2::Delete(){ //ELENA  void ToFLevel2::Delete(Option_t *t){ //ELENA
124    //    //
125    if(ToFTrk){    if(ToFTrk){
126        ToFTrk->Delete(); //ELENA        ToFTrk->Delete(); //ELENA
# Line 213  Int_t ToFLevel2::GetNHitPaddles(Int_t pl Line 220  Int_t ToFLevel2::GetNHitPaddles(Int_t pl
220      return npad;      return npad;
221  };  };
222    
223    //wm Nov 08
224  //gf Apr 07  //gf Apr 07
225  /**  /**
226   * Method to get the mean dEdx from a given ToF plane. This current version   * Method to get the mean dEdx from a ToF layer - ATTENTION:
227   * is just summing up all PMT signals, which will not give proper results,   * It will sum up the dEdx of all the paddles, but since by definition
228   *  and needs a revision.   * only the paddle hitted by the track gets a dEdx value and the other
229     * paddles are set to zero, the output is just the dEdx of the hitted
230     * paddle in each layer!
231     * The "adcfl" option is not very useful (an artificial dEdx is per
232     * definition= 1 mip and not a real measurement), anyway left in the code
233   * @param notrack Track Number   * @param notrack Track Number
234   * @param plane Plane index (0,1,2,3,4,5)   * @param plane Plane index (0,1,2,3,4,5)
235   * @param adcflag in the plane (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )   * @param adcflag in the plane (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
# Line 378  void ToFLevel2::GetPMTIndex(Int_t ind, I Line 389  void ToFLevel2::GetPMTIndex(Int_t ind, I
389    
390    
391    
392    //  wm Nov 08 revision - saturation values included
393  /// gf Apr 07  /// gf Apr 07
   
394  /**  /**
395   * Method to get the dEdx from a given ToF paddle.   * Method to get the dEdx from a given ToF paddle.
396     * If two PMTs are good, the mean dEdx of both PMTs is taken, otherwise
397     * just the dEdx of the "good" PMT. If both PMTs are above saturation => dEdx=1000
398   * @param notrack Track Number   * @param notrack Track Number
399   * @param Paddle index (0,1,...,23).   * @param Paddle index (0,1,...,23).
400   * @param adcflag in the paddle (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )   * @param adcflag in the paddle (100<-> independent of the adcflag; !=0&&!=100 <-> at least one PMT with adcflag!=0; )
# Line 391  void ToFLevel2::GetPMTIndex(Int_t ind, I Line 403  void ToFLevel2::GetPMTIndex(Int_t ind, I
403   */   */
404  void ToFLevel2::GetdEdxPaddle(Int_t notrack, Int_t paddleid, Int_t adcfl, Float_t &PadEdx, Int_t &SatWarning){  void ToFLevel2::GetdEdxPaddle(Int_t notrack, Int_t paddleid, Int_t adcfl, Float_t &PadEdx, Int_t &SatWarning){
405    
406    /*
407    Float_t  PMTsat[48] = {
408    3162.14, 3165.48, 3153.85, 3085.73, 3089.65, 3107.64, 3097.52, 3078.37,
409    3130.05, 3087.07, 3112.22, 3102.92, 3080.58, 3092.55, 3087.94, 3125.03,
410    3094.09, 3143.16, 3125.51, 3181.27, 3092.09, 3124.98, 3069.3, 3095.53,
411    3097.11, 3133.53, 3114.73, 3113.01, 3091.19, 3097.99, 3033.84, 3134.98,
412    3081.37, 3111.04, 3066.77, 3108.17, 3133, 3111.06, 3052.52, 3140.66,
413    3106.33, 3094.85, 3150.85, 3118.8, 3096.24, 3118.47,3111.36, 3117.11 } ;
414    */
415    
416    // new values from Napoli dec 2008
417    Float_t  PMTsat[48] = {
418    3176.35,3178.19,3167.38,3099.73,3117.00,3126.29,3111.44,3092.27,
419    3146.48,3094.41,3132.13,3115.37,3099.32,3110.97,3111.80,3143.14,
420    3106.72,3153.44,3136.00,3188.96,3104.73,3140.45,3073.18,3106.62,
421    3112.48,3146.92,3127.24,3136.52,3109.59,3112.89,3045.15,3147.26,
422    3095.92,3121.05,3083.25,3123.62,3150.92,3125.30,3067.60,3160.18,
423    3119.36,3108.92,3164.77,3133.64,3111.47,3131.98,3128.87,3135.56 };
424    
425    for (Int_t i=0; i<48;i++) PMTsat[i] = PMTsat[i] - 5.;  // safety margin
426    
427    
428    PadEdx = 0.;    PadEdx = 0.;
429    SatWarning = 1000;  //  SatWarning = 1000;
430      SatWarning = 0;   // 0=good, increase for each bad PMT
431    
432    Float_t dEdx[48] = {0};    Float_t dEdx[48] = {0};
433    Int_t pmt_id = -1;    Int_t pmt_id = -1;
# Line 424  void ToFLevel2::GetdEdxPaddle(Int_t notr Line 459  void ToFLevel2::GetdEdxPaddle(Int_t notr
459        adcraw[pmtright] = pmt->adc;        adcraw[pmtright] = pmt->adc;
460      }      }
461    }    }
462    
463        
464    for (Int_t i=0; i<trk->npmtadc; i++){    for (Int_t i=0; i<trk->npmtadc; i++){
465    
# Line 436  void ToFLevel2::GetdEdxPaddle(Int_t notr Line 472  void ToFLevel2::GetdEdxPaddle(Int_t notr
472      }      }
473    }    }
474    
475    if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1;  
476        //  if( adcraw[pmtleft] >3000 || adcraw[pmtright] >3000)SatWarning=1;  //old version
477    if(dEdx[pmtleft]!=0 && dEdx[pmtright]!=0){  
478      PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5;  // Increase SatWarning Counter for each PMT>Sat
479    }    if( adcraw[pmtleft] > PMTsat[pmtleft])SatWarning++;  
480    if(dEdx[pmtleft]==0 && dEdx[pmtright]!=0){    if( adcraw[pmtright] > PMTsat[pmtright])SatWarning++;
481      PadEdx = dEdx[pmtright];  
482    }  // if ADC  > sat set dEdx=1000
483    if(dEdx[pmtleft]!=0 && dEdx[pmtright]==0){    if( adcraw[pmtleft] > PMTsat[pmtleft]) dEdx[pmtleft] = 1000.;
484      PadEdx = dEdx[pmtleft];    if( adcraw[pmtright] > PMTsat[pmtright]) dEdx[pmtright] = 1000. ;
485    }  
486    // if two PMT are good, take mean dEdx, otherwise only the good dEdx
487      if(dEdx[pmtleft]<1000 && dEdx[pmtright]<1000) PadEdx = (dEdx[pmtleft]+dEdx[pmtright])*0.5;
488      if(dEdx[pmtleft]==1000 && dEdx[pmtright]<1000) PadEdx = dEdx[pmtright];  
489      if(dEdx[pmtleft]<1000 && dEdx[pmtright]==1000) PadEdx = dEdx[pmtleft];
490        
   return;  
491  };  };
492  //  //
493    
# Line 505  TString ToFLevel2::GetPMTName(Int_t ind) Line 544  TString ToFLevel2::GetPMTName(Int_t ind)
544        
545  };  };
546    
547    // wm jun 08
 // gf Apr 07  
548  Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane){  Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane){
549    return GetPaddleIdOfTrack(xtr ,ytr ,plane, 0.4);
550    }
551    
552    // gf Apr 07
553    Int_t ToFLevel2::GetPaddleIdOfTrack(Float_t xtr, Float_t ytr, Int_t plane, Float_t margin){
554      
555    Double_t xt,yt,xl,xh,yl,yh;    Double_t xt,yt,xl,xh,yl,yh;
556        
557    Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};    Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85};
# Line 538  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 581  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
581      yh =  33.0/2. ;      yh =  33.0/2. ;
582      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
583        for (Int_t i1=0; i1<8;i1++){        for (Int_t i1=0; i1<8;i1++){
584          xl = tof11_x[i1] - (5.1-0.4)/2. ;          xl = tof11_x[i1] - (5.1-margin)/2. ;
585          xh = tof11_x[i1] + (5.1-0.4)/2. ;          xh = tof11_x[i1] + (5.1-margin)/2. ;
586          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
587        }        }
588      }      }
# Line 556  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 599  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
599            
600      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
601        for (Int_t i1=0; i1<6;i1++){        for (Int_t i1=0; i1<6;i1++){
602          yl = tof12_y[i1] - (5.5-0.4)/2. ;          yl = tof12_y[i1] - (5.5-margin)/2. ;
603          yh = tof12_y[i1] + (5.5-0.4)/2. ;          yh = tof12_y[i1] + (5.5-margin)/2. ;
604          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
605        }        }
606      }      }
# Line 574  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 617  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
617            
618      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
619        for (Int_t i1=0; i1<2;i1++){        for (Int_t i1=0; i1<2;i1++){
620          yl = tof21_y[i1] - (7.5-0.4)/2. ;          yl = tof21_y[i1] - (7.5-margin)/2. ;
621          yh = tof21_y[i1] + (7.5-0.4)/2. ;          yh = tof21_y[i1] + (7.5-margin)/2. ;
622          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh))  paddleidoftrack=i1;
623        }        }
624      }      }
# Line 591  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 634  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
634            
635      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
636        for (Int_t i1=0; i1<2;i1++){        for (Int_t i1=0; i1<2;i1++){
637          xl = tof22_x[i1] - (9.0-0.4)/2. ;          xl = tof22_x[i1] - (9.0-margin)/2. ;
638          xh = tof22_x[i1] + (9.0-0.4)/2. ;          xh = tof22_x[i1] + (9.0-margin)/2. ;
639          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
640        }        }
641      }      }
# Line 608  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 651  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
651            
652      if ((yt>yl)&&(yt<yh)) {      if ((yt>yl)&&(yt<yh)) {
653        for (Int_t i1=0; i1<3;i1++){        for (Int_t i1=0; i1<3;i1++){
654          xl = tof31_x[i1] - (6.0-0.4)/2. ;          xl = tof31_x[i1] - (6.0-margin)/2. ;
655          xh = tof31_x[i1] + (6.0-0.4)/2. ;          xh = tof31_x[i1] + (6.0-margin)/2. ;
656          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;          if ((xt>xl)&&(xt<xh))  paddleidoftrack=i1;
657        }        }
658      }      }
# Line 625  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa Line 668  Int_t ToFLevel2::GetPaddleIdOfTrack(Floa
668            
669      if ((xt>xl)&&(xt<xh)) {      if ((xt>xl)&&(xt<xh)) {
670        for (Int_t i1=0; i1<3;i1++){        for (Int_t i1=0; i1<3;i1++){
671          yl = tof32_y[i1] - (5.0-0.4)/2. ;          yl = tof32_y[i1] - (5.0-margin)/2. ;
672          yh = tof32_y[i1] + (5.0-0.4)/2. ;          yh = tof32_y[i1] + (5.0-margin)/2. ;
673          if ((yt>yl)&&(yt<yh)) paddleidoftrack=i1;          if ((yt>yl)&&(yt<yh)) paddleidoftrack=i1;
674        }        }
675      }      }
# Line 1046  Int_t ToFLevel2::GetNPaddle(Int_t plane) Line 1089  Int_t ToFLevel2::GetNPaddle(Int_t plane)
1089    
1090  }  }
1091    
 ////////////////////////////////////////////////////  
1092    
1093    
1094    /// wm feb 08
1095    
1096    /**
1097     * Method to calculate Beta from the 12 single measurements
1098     * we check the individual weights for artificial TDC values, then calculate
1099     * am mean beta for the first time. In a second step we loop again through
1100     * the single measurements, checking for the residual from the mean
1101     * The cut on the residual reject measurements > "x"-sigma. A chi2 value is
1102     * calculated, furthermore a "quality" value by adding the weights which
1103     * are finally used. If all measurements are taken, "quality" will be = 22.47.
1104     * A chi2 cut around 3-4 and a quality-cut > 20 is needed for clean beta
1105     * measurements like antiprotons etc.
1106     * The Level2 output is derived in the fortran routines using: 10.,10.,20.
1107     * @param notrack Track Number
1108     * @param cut on residual: difference between single measurement and mean
1109     * @param cut on "quality"
1110     * @param cut on chi2
1111     */
1112    
1113    Float_t ToFLevel2::CalcBeta(Int_t notrack, Float_t resmax, Float_t qualitycut, Float_t chi2cut){
1114    
1115    //  cout<<" in CalcBeta "<<resmax<<" "<<chi2cut<<" "<<qualitycut<<endl;
1116    
1117      Float_t bxx = 100.;
1118      //
1119      ToFTrkVar *trk = GetToFTrkVar(notrack);
1120      if(!trk) return 0; //ELENA
1121    
1122    
1123      Float_t chi2,xhelp,beta_mean;
1124      Float_t w_i[12],quality,sw,sxw,res,betachi,beta_mean_inv;
1125      Float_t b[12],tdcfl;
1126      Int_t  pmt_id,pmt_plane;
1127    
1128      for (Int_t i=0; i<12; i++){
1129        b[i] = trk->beta[i];
1130                                  }
1131          
1132    
1133    //========================================================================
1134    //---  Find out ToF layers with artificial TDC values & fill vector    ---
1135    //========================================================================
1136    
1137    Float_t  w_il[6];
1138    
1139         for (Int_t jj=0; jj<6;jj++) {
1140             w_il[jj] = 1000.;
1141                                     }
1142    
1143    
1144      for (Int_t i=0; i<trk->npmttdc; i++){
1145        //
1146        pmt_id = (trk->pmttdc).At(i);
1147        pmt_plane = GetPlaneIndex(pmt_id);
1148        tdcfl = (trk->tdcflag).At(i);
1149        if (w_il[pmt_plane] != 1.) w_il[pmt_plane] = tdcfl; //tdcflag
1150                                         };
1151      
1152    //========================================================================
1153    //---  Set weights for the 12 measurements using information for top and bottom:
1154    //---  if no measurements: weight = set to very high value=> not used
1155    //---  top or bottom artificial: weight*sqrt(2)
1156    //---  top and bottom artificial: weight*sqrt(2)*sqrt(2)
1157    //========================================================================
1158    
1159    Int_t itop[12] = {0,0,1,1,2,2,3,3,0,0,1,1};
1160    Int_t ibot[12] = {4,5,4,5,4,5,4,5,2,3,2,3};
1161    
1162         xhelp= 1E09;
1163      
1164         for (Int_t jj=0; jj<12;jj++) {
1165         if (jj<4)           xhelp = 0.11;    // S1-S3
1166         if ((jj>3)&&(jj<8)) xhelp = 0.18;    // S2-S3
1167         if (jj>7)           xhelp = 0.28;    // S1-S2
1168         if ((w_il[itop[jj]] == 1000.) && (w_il[ibot[jj]] == 1000.)) xhelp = 1E09;
1169         if ((w_il[itop[jj]] == 1) || (w_il[ibot[jj]] == 1.)) xhelp = xhelp*1.414 ;
1170         if ((w_il[itop[jj]] == 1) && (w_il[ibot[jj]] == 1.)) xhelp = xhelp*2. ;
1171    
1172         w_i[jj] = 1./xhelp;
1173                                      }
1174    
1175    
1176    //========================================================================
1177    //--- Calculate mean beta for the first time -----------------------------
1178    //--- We are using "1/beta" since its error is gaussian ------------------
1179    //========================================================================
1180    
1181          Int_t icount=0;
1182          sw=0.;
1183          sxw=0.;
1184          beta_mean=100.;
1185    
1186              for (Int_t jj=0; jj<12;jj++){
1187            if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.))
1188             {
1189                icount= icount+1;
1190                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1191                sw =sw + w_i[jj]*w_i[jj] ;
1192    
1193             }
1194             }
1195    
1196          if (icount>0) beta_mean=1./(sxw/sw);
1197          beta_mean_inv = 1./beta_mean;
1198    
1199    //========================================================================
1200    //--- Calculate beta for the second time, use residuals of the single
1201    //--- measurements to get a chi2 value
1202    //========================================================================
1203    
1204          icount=0;
1205          sw=0.;
1206          sxw=0.;
1207          betachi = 100.;
1208          chi2 = 0.;
1209          quality=0.;
1210    
1211    
1212              for (Int_t jj=0; jj<12;jj++){
1213           if ((fabs(1./b[jj])>0.1)&&(fabs(1./b[jj])<15.)&&(w_i[jj]>0.01)) {
1214                res = beta_mean_inv - (1./b[jj]) ;
1215                if (fabs(res*w_i[jj])<resmax)          {;
1216                chi2 = chi2 + pow((res*w_i[jj]),2) ;
1217                icount= icount+1;
1218                sxw=sxw + (1./b[jj])*w_i[jj]*w_i[jj] ;
1219                sw =sw + w_i[jj]*w_i[jj] ;
1220                                                   }
1221                                                                            }
1222                                          }
1223          quality = sqrt(sw) ;
1224    
1225          if (icount==0) chi2 = 1000.;
1226          if (icount>0) chi2 = chi2/(icount) ;
1227          if (icount>0) betachi=1./(sxw/sw);
1228    
1229       bxx = 100.;
1230       if ((chi2 < chi2cut)&&(quality>qualitycut)) bxx = betachi;
1231      //
1232      return(bxx);
1233    };
1234    
1235    
1236    ////////////////////////////////////////////////////
1237    ////////////////////////////////////////////////////
1238    
1239    
1240  /**  /**
1241   * 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.23
  ViewVC Help
Powered by ViewVC 1.1.23