1 /***********************************************************************************************************
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. *
10 ***********************************************************************************************************/
12 #include "TObjArray.h"
17 #include "TGraphErrors.h"
21 #include "Cal/AliTRDCalPID.h"
22 #include "AliCDBManager.h"
23 #include "AliCDBEntry.h"
24 #include "AliESDtrack.h"
26 #include "AliTRDpidUtil.h"
28 ClassImp(AliTRDpidUtil)
31 Float_t AliTRDpidUtil::fgEleEffi = 0.9;
33 //________________________________________________________________________
34 AliTRDpidUtil::AliTRDpidUtil()
42 // Default constructor
46 //________________________________________________________________________
47 Bool_t AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
48 // Double_t AliTRDpidUtil::GetError()
51 // Calculates the pion efficiency
54 histo1 -> SetLineColor(kRed);
55 histo2 -> SetLineColor(kBlue);
56 if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){
57 AliWarning("Histo with no entries !");
60 histo1 -> Scale(1./histo1->GetEntries());
61 histo2 -> Scale(1./histo2->GetEntries());
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));
71 // calculate electron efficiency of each bin
72 for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){
74 aBinSumE = histo1 -> GetBinContent(abinE);
76 sumElecE[abinE] = sumElecE[abinE+1] + aBinSumE;
78 if((sumElecE[abinE] >= fgEleEffi) && (bFirst == 1)){
81 fCalcEleEffi = (sumElecE[cbinE]);
85 fThreshold = histo1 -> GetBinCenter(cbinE);
87 // calculate pion efficiency of each bin
88 for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
90 bBinSumE = histo2 -> GetBinContent(bbinE);
92 sumPionsE[bbinE] = sumPionsE[bbinE+1] + bBinSumE;
94 fPionEffi = (sumPionsE[cbinE]);
99 // pion efficiency vs electron efficiency
100 TGraph gEffis(kBins, sumElecE, sumPionsE);
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);
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!");
111 fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
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)));
119 //__________________________________________________________________________
120 Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
123 // returns the momentum bin for a given momentum
126 Int_t pBin1 = -1; // check bin1
127 Int_t pBin2 = -1; // check bin2
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
134 // calculate momentum bin for non extremal momenta
135 for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
136 if(p < AliTRDCalPID::GetMomentum(iMomBin)){
143 if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
155 //__________________________________________________________________________
156 Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){
158 // Do PID decision for the TRD based on 90% Electron efficiency threshold
160 // Author: Markus Fasel (M.Fasel@gsi.de)
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);
171 Double_t pid_probs[5];
172 track->GetTRDpid(pid_probs);
173 if(pid_probs[AliPID::kElectron] >= threshold) return kTRUE;
177 //__________________________________________________________________________
178 Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){
180 // Returns the pion efficiency at 90% electron efficiency from the OCDB
182 // Author: Markus Fasel (M.Fasel@gsi.de)
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);
193 //________________________________________________________________________
194 Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
196 // Private Helper function to get the paticle species (ALICE notation)
200 switch(TMath::Abs(pdg)){
202 species = AliPID::kElectron;
205 species = AliPID::kMuon;
208 species = AliPID::kPion;
211 species = AliPID::kKaon;
214 species = AliPID::kProton;