TRD/AliTRDpidUtil.cxx

   1 /***********************************************************************************************************
   2  *                                                                                                         *
   3  * Helper class for TRD PID efficiency calculation Calculation of the hadron efficiency depenedent on      *
   4  * momentum and of the errors implemented in function CalculatePionEff. The pion efficiency is based on a  *
   5  * predefined electron efficiency. The default is 90%. To change the, one has to call the function         *
   6  * SetElectronEfficiency.                                                                                  *
   7  * Other Helper functions decide based on 90% electron efficiency whether a certain track is accepted      *
   8  * as Electron Track. The reference data is stored in the TRD OCDB.                                        *
   9  *                                                                                                         *
  10  ***********************************************************************************************************/
  11 #include "TObject.h"
  12 #include "TObjArray.h"
  13 #include "TMath.h"
  14 #include "TF1.h"
  15 #include "TH1F.h"
  16 #include "TGraph.h"
  17 #include "TGraphErrors.h"
  18 #include "TPDGCode.h"
  19
  20 #include "AliLog.h"
  21 #include "Cal/AliTRDCalPID.h"
  22 #include "AliCDBManager.h"
  23 #include "AliCDBEntry.h"
  24 #include "AliESDtrack.h"
  25 #include "AliPID.h"
  26 #include "AliTRDpidUtil.h"
  27
  28 ClassImp(AliTRDpidUtil)
  29
  30
  31 Float_t AliTRDpidUtil::fgEleEffi = 0.9;
  32
  33 //________________________________________________________________________
  34 AliTRDpidUtil::AliTRDpidUtil()
  35   : TObject()
  36   ,fCalcEleEffi(0.)
  37   ,fPionEffi(-1.)
  38   ,fError(-1.)
  39   ,fThreshold(-1.)
  40 {
  41   //
  42   // Default constructor
  43   //
  44 }
  45
  46 //________________________________________________________________________
  47 Bool_t  AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
  48 // Double_t  AliTRDpidUtil::GetError()
  49 {
  50   //
  51   // Calculates the pion efficiency
  52   //
  53
  54   histo1 -> SetLineColor(kRed);
  55   histo2 -> SetLineColor(kBlue);
  56   if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){
  57     AliWarning("Histo with no entries !");
  58     return kFALSE;
  59   }
  60   histo1 -> Scale(1./histo1->GetEntries());
  61   histo2 -> Scale(1./histo2->GetEntries());
  62
  63   Int_t abinE, bbinE, cbinE = -1;
  64   Double_t aBinSumE, bBinSumE;                  // content of a single bin
  65   Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
  66   Double_t sumElecE[kBins+2], sumPionsE[kBins+2];  // array of the integrated sum in each bin
  67   memset(sumElecE, 0, (kBins+2)*sizeof(Double_t));
  68   memset(sumPionsE, 0, (kBins+2)*sizeof(Double_t));
  69
  70
  71   // calculate electron efficiency of each bin
  72   for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){
  73    aBinSumE = 0;
  74     aBinSumE = histo1 -> GetBinContent(abinE);
  75
  76     sumElecE[abinE] = sumElecE[abinE+1] + aBinSumE;
  77
  78     if((sumElecE[abinE] >= fgEleEffi) && (bFirst == 1)){
  79       bFirst = 0;
  80       cbinE = abinE;
  81       fCalcEleEffi = (sumElecE[cbinE]);
  82     }
  83   }
  84
  85   fThreshold = histo1 -> GetBinCenter(cbinE);
  86
  87   // calculate pion efficiency of each bin
  88   for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
  89     bBinSumE = 0;
  90     bBinSumE = histo2 -> GetBinContent(bbinE);
  91
  92     sumPionsE[bbinE] = sumPionsE[bbinE+1] + bBinSumE;
  93     if(bbinE == cbinE){
  94       fPionEffi = (sumPionsE[cbinE]);
  95     }
  96   }
  97
  98
  99   // pion efficiency vs electron efficiency
 100   TGraph gEffis(kBins, sumElecE, sumPionsE);
 101
 102   // use fit function to get derivate of the TGraph for error calculation
 103   TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05);
 104   gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.05);
 105
 106   // return the error of the pion efficiency
 107   if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
 108     AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
 109     return kFALSE;
 110   }
 111   fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
 112
 113 //   AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
 114 //   AliInfo(Form("Derivative at %4.2f : %f\n", fgEleEffi, f1.Derivative(fgEleEffi)));
 115   return kTRUE;
 116 }
 117
 118
 119 //__________________________________________________________________________
 120 Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
 121 {
 122   //
 123   // returns the momentum bin for a given momentum
 124   //
 125
 126   Int_t pBin1 = -1;                                    // check bin1
 127   Int_t pBin2 = -1;                                    // check bin2
 128
 129   if(p < 0) return -1;                                 // return -1 if momentum < 0
 130   if(p < AliTRDCalPID::GetMomentum(0)) return 0;                                      // smallest momentum bin
 131   if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin
 132
 133
 134   // calculate momentum bin for non extremal momenta
 135   for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 136     if(p < AliTRDCalPID::GetMomentum(iMomBin)){
 137       pBin1 = iMomBin - 1;
 138       pBin2 = iMomBin;
 139     }
 140     else
 141       continue;
 142
 143     if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
 144        return pBin2;
 145     }
 146     else{
 147       return pBin1;
 148     }
 149   }
 150
 151   return -1;
 152 }
 153
 154
 155 //__________________________________________________________________________
 156 Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){
 157   //
 158   // Do PID decision for the TRD based on 90% Electron efficiency threshold
 159   //
 160   // Author: Markus Fasel (M.Fasel@gsi.de)
 161   //
 162   if(method == kESD) method = kNN;
 163   TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"};
 164   AliCDBManager *cdb = AliCDBManager::Instance();
 165   AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
 166   TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
 167   TH1 * thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
 168   Double_t threshold = thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
 169
 170   // Do Decision
 171   Double_t pid_probs[5];
 172   track->GetTRDpid(pid_probs);
 173   if(pid_probs[AliPID::kElectron] >= threshold) return kTRUE;
 174   return kFALSE;
 175 }
 176
 177 //__________________________________________________________________________
 178 Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){
 179   //
 180   // Returns the pion efficiency at 90% electron efficiency from the OCDB
 181   //
 182   // Author: Markus Fasel (M.Fasel@gsi.de)
 183   //
 184   if(method == kESD) method = kNN;
 185   TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"};
 186   AliCDBManager *cdb = AliCDBManager::Instance();
 187   AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
 188   TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
 189   TH1 * thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
 190   return thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
 191 }
 192
 193 //________________________________________________________________________
 194 Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
 195   //
 196   // Private Helper function to get the paticle species (ALICE notation)
 197   // from the Pdg code
 198   //
 199   Int_t species;
 200   switch(TMath::Abs(pdg)){
 201   case kElectron:
 202     species = AliPID::kElectron;
 203     break;
 204   case kMuonMinus:
 205     species = AliPID::kMuon;
 206     break;
 207   case kPiPlus:
 208     species = AliPID::kPion;
 209     break;
 210   case kKPlus:
 211     species = AliPID::kKaon;
 212     break;
 213   case kProton:
 214     species = AliPID::kProton;
 215     break;
 216   default:
 217     species = -1;
 218   }
 219   return species;
 220 }
 221