TRD/AliTRDpidUtil.cxx

   1 #include "TObject.h"
   2 #include "TMath.h"
   3 #include "TF1.h"
   4 #include "TH1F.h"
   5 #include "TGraph.h"
   6 #include "TGraphErrors.h"
   7
   8 #include "AliLog.h"
   9 #include "Cal/AliTRDCalPID.h"
  10 #include "AliTRDpidUtil.h"
  11
  12 ClassImp(AliTRDpidUtil)
  13
  14
  15 Float_t AliTRDpidUtil::fEleEffi = 0.9;
  16
  17 //________________________________________________________________________
  18 AliTRDpidUtil::AliTRDpidUtil()
  19   : TObject()
  20   ,fCalcEleEffi(0.)
  21   ,fPionEffi(-1.)
  22   ,fError(-1.)
  23   ,fThreshold(-1.)
  24 {
  25   //
  26   // Default constructor
  27   //
  28 }
  29
  30
  31
  32 //________________________________________________________________________
  33 void  AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
  34 // Double_t  AliTRDpidUtil::GetError()
  35 {
  36   //
  37   // Calculates the pion efficiency
  38   //
  39
  40   histo1 -> SetLineColor(kRed);
  41   histo2 -> SetLineColor(kBlue);
  42   AliInfo(Form("Histo1[%d] Histo2[%d]", (Int_t)histo1 -> GetEntries(), (Int_t)histo2 -> GetEntries()));
  43   if(!(histo1 -> GetEntries() > 0 && histo2 -> GetEntries() > 0)){
  44     AliWarning("Histo has no Entries !");
  45     return;
  46   }
  47   histo1 -> Scale(1./histo1->GetEntries());
  48   histo2 -> Scale(1./histo2->GetEntries());
  49
  50   Int_t abinE, bbinE, cbinE = -1;
  51   Double_t aBinSumE, bBinSumE;                  // content of a single bin
  52   Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
  53   Double_t SumElecsE[kBins+2], SumPionsE[kBins+2];  // array of the integrated sum in each bin
  54   memset(SumElecsE, 0, (kBins+2)*sizeof(Double_t));
  55   memset(SumPionsE, 0, (kBins+2)*sizeof(Double_t));
  56
  57
  58   // calculate electron efficiency of each bin
  59   for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){
  60    aBinSumE = 0;
  61     aBinSumE = histo1 -> GetBinContent(abinE);
  62
  63     SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
  64
  65     if((SumElecsE[abinE] >= fEleEffi) && (bFirst == 1)){
  66       bFirst = 0;
  67       cbinE = abinE;
  68       fCalcEleEffi = (SumElecsE[cbinE]);
  69     }
  70   }
  71
  72   fThreshold = histo1 -> GetBinCenter(cbinE);
  73
  74   // calculate pion efficiency of each bin
  75   for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
  76     bBinSumE = 0;
  77     bBinSumE = histo2 -> GetBinContent(bbinE);
  78
  79     SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
  80     if(bbinE == cbinE){
  81       fPionEffi = (SumPionsE[cbinE]);
  82     }
  83   }
  84
  85
  86   // pion efficiency vs electron efficiency
  87   TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE);
  88
  89   // use fit function to get derivate of the TGraph for error calculation
  90   TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
  91   gEffis -> Fit("f1","Q","",fEleEffi-.05, fEleEffi+.05);
  92
  93   // return the error of the pion efficiency
  94   if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
  95     AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
  96     return;
  97   }
  98   fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1 -> Derivative(fEleEffi))*(f1 -> Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
  99
 100   AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
 101   AliInfo(Form("Derivative at %4.2f : %f\n", fEleEffi, f1 -> Derivative(fEleEffi)));
 102
 103 }
 104
 105
 106 //__________________________________________________________________________
 107 Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
 108 {
 109   //
 110   // returns the momentum bin for a given momentum
 111   //
 112
 113   Int_t pBin1 = -1;                                    // check bin1
 114   Int_t pBin2 = -1;                                    // check bin2
 115
 116   if(p < 0) return -1;                                 // return -1 if momentum < 0
 117   if(p < AliTRDCalPID::GetMomentum(0)) return 0;                                      // smallest momentum bin
 118   if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin
 119
 120
 121   // calculate momentum bin for non extremal momenta
 122   for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 123     if(p < AliTRDCalPID::GetMomentum(iMomBin)){
 124       pBin1 = iMomBin - 1;
 125       pBin2 = iMomBin;
 126     }
 127     else
 128       continue;
 129
 130     if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
 131        return pBin2;
 132     }
 133     else{
 134       return pBin1;
 135     }
 136   }
 137
 138   return -1;
 139 }
 140