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 Bool_t  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   if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){
  43     AliWarning("Histo with no entries !");
  44     return kFALSE;
  45   }
  46   histo1 -> Scale(1./histo1->GetEntries());
  47   histo2 -> Scale(1./histo2->GetEntries());
  48
  49   Int_t abinE, bbinE, cbinE = -1;
  50   Double_t aBinSumE, bBinSumE;                  // content of a single bin
  51   Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
  52   Double_t SumElecsE[kBins+2], SumPionsE[kBins+2];  // array of the integrated sum in each bin
  53   memset(SumElecsE, 0, (kBins+2)*sizeof(Double_t));
  54   memset(SumPionsE, 0, (kBins+2)*sizeof(Double_t));
  55
  56
  57   // calculate electron efficiency of each bin
  58   for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){
  59    aBinSumE = 0;
  60     aBinSumE = histo1 -> GetBinContent(abinE);
  61
  62     SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
  63
  64     if((SumElecsE[abinE] >= fEleEffi) && (bFirst == 1)){
  65       bFirst = 0;
  66       cbinE = abinE;
  67       fCalcEleEffi = (SumElecsE[cbinE]);
  68     }
  69   }
  70
  71   fThreshold = histo1 -> GetBinCenter(cbinE);
  72
  73   // calculate pion efficiency of each bin
  74   for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
  75     bBinSumE = 0;
  76     bBinSumE = histo2 -> GetBinContent(bbinE);
  77
  78     SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
  79     if(bbinE == cbinE){
  80       fPionEffi = (SumPionsE[cbinE]);
  81     }
  82   }
  83
  84
  85   // pion efficiency vs electron efficiency
  86   TGraph gEffis(kBins, SumElecsE, SumPionsE);
  87
  88   // use fit function to get derivate of the TGraph for error calculation
  89   TF1 f1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
  90   gEffis.Fit(&f1, "Q", "", fEleEffi-.05, fEleEffi+.05);
  91
  92   // return the error of the pion efficiency
  93   if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
  94     AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
  95     return kFALSE;
  96   }
  97   fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fEleEffi))*(f1.Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
  98
  99 //   AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
 100 //   AliInfo(Form("Derivative at %4.2f : %f\n", fEleEffi, f1.Derivative(fEleEffi)));
 101   return kTRUE;
 102 }
 103
 104
 105 //__________________________________________________________________________
 106 Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
 107 {
 108   //
 109   // returns the momentum bin for a given momentum
 110   //
 111
 112   Int_t pBin1 = -1;                                    // check bin1
 113   Int_t pBin2 = -1;                                    // check bin2
 114
 115   if(p < 0) return -1;                                 // return -1 if momentum < 0
 116   if(p < AliTRDCalPID::GetMomentum(0)) return 0;                                      // smallest momentum bin
 117   if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin
 118
 119
 120   // calculate momentum bin for non extremal momenta
 121   for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 122     if(p < AliTRDCalPID::GetMomentum(iMomBin)){
 123       pBin1 = iMomBin - 1;
 124       pBin2 = iMomBin;
 125     }
 126     else
 127       continue;
 128
 129     if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
 130        return pBin2;
 131     }
 132     else{
 133       return pBin1;
 134     }
 135   }
 136
 137   return -1;
 138 }
 139