The present commit corresponds to an important change in the way the
[u/mrichter/AliRoot.git] / TRD / AliTRDpidUtil.cxx
CommitLineData
0e60774c 1#include "TObject.h"
0affc9e5 2#include "TObjArray.h"
0e60774c 3#include "TMath.h"
4#include "TF1.h"
5#include "TH1F.h"
6#include "TGraph.h"
7#include "TGraphErrors.h"
0d83b3a5 8#include "TPDGCode.h"
0e60774c 9
10#include "AliLog.h"
11#include "Cal/AliTRDCalPID.h"
0affc9e5 12#include "AliCDBManager.h"
13#include "AliCDBEntry.h"
14#include "AliESDtrack.h"
15#include "AliPID.h"
0e60774c 16#include "AliTRDpidUtil.h"
17
18ClassImp(AliTRDpidUtil)
19
20
21Float_t AliTRDpidUtil::fEleEffi = 0.9;
22
23//________________________________________________________________________
24AliTRDpidUtil::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//________________________________________________________________________
c46a7947 39Bool_t AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
0e60774c 40// Double_t AliTRDpidUtil::GetError()
41{
42 //
43 // Calculates the pion efficiency
44 //
45
46 histo1 -> SetLineColor(kRed);
47 histo2 -> SetLineColor(kBlue);
c46a7947 48 if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){
49 AliWarning("Histo with no entries !");
50 return kFALSE;
0e60774c 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
422a2dc0 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));
0e60774c 61
62
63 // calculate electron efficiency of each bin
422a2dc0 64 for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){
65 aBinSumE = 0;
0e60774c 66 aBinSumE = histo1 -> GetBinContent(abinE);
67
68 SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
422a2dc0 69
0e60774c 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
422a2dc0 80 for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
0e60774c 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
c46a7947 92 TGraph gEffis(kBins, SumElecsE, SumPionsE);
0e60774c 93
94 // use fit function to get derivate of the TGraph for error calculation
c46a7947 95 TF1 f1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
96 gEffis.Fit(&f1, "Q", "", fEleEffi-.05, fEleEffi+.05);
422a2dc0 97
0e60774c 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!");
c46a7947 101 return kFALSE;
0e60774c 102 }
c46a7947 103 fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fEleEffi))*(f1.Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
422a2dc0 104
c46a7947 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;
0e60774c 108}
109
110
111//__________________________________________________________________________
112Int_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
0affc9e5 146
147//__________________________________________________________________________
0d83b3a5 148Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){
0affc9e5 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//__________________________________________________________________________
0d83b3a5 170Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){
0affc9e5 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}
0d83b3a5 184
185//________________________________________________________________________
186Int_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