TRD/AliTRDpidUtil.cxx

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