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(TH1F* histo1, TH1F* 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], SumPionsE[kBins];  // array of the integrated sum in each bin
  54   memset(SumElecsE, 0, kBins*sizeof(Double_t));
  55   memset(SumPionsE, 0, kBins*sizeof(Double_t));
  56
  57
  58   // calculate electron efficiency of each bin
  59   for (abinE = (histo1 -> GetNbinsX())-2; abinE >= 0; abinE--){
  60     aBinSumE = 0;
  61     aBinSumE = histo1 -> GetBinContent(abinE);
  62
  63     SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
  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())-2; 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 = new TGraph(kBins, SumElecsE, SumPionsE);
  87
  88   // use fit function to get derivate of the TGraph for error calculation
  89   TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
  90   gEffis -> Fit("f1","Q","",fEleEffi-.05, fEleEffi+.05);
  91   AliInfo(Form("Derivative at %4.2f : %f", fEleEffi, f1 -> Derivative(fEleEffi)));
  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
 101
 102 //__________________________________________________________________________
 103 Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
 104 {
 105   //
 106   // returns the momentum bin for a given momentum
 107   //
 108
 109   Int_t pBin1 = -1;                                    // check bin1
 110   Int_t pBin2 = -1;                                    // check bin2
 111
 112   if(p < 0) return -1;                                 // return -1 if momentum < 0
 113   if(p < AliTRDCalPID::GetMomentum(0)) return 0;                                      // smallest momentum bin
 114   if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin
 115
 116
 117   // calculate momentum bin for non extremal momenta
 118   for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 119     if(p < AliTRDCalPID::GetMomentum(iMomBin)){
 120       pBin1 = iMomBin - 1;
 121       pBin2 = iMomBin;
 122     }
 123     else
 124       continue;
 125
 126     if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
 127        return pBin2;
 128     }
 129     else{
 130       return pBin1;
 131     }
 132   }
 133
 134   return -1;
 135 }
 136