]>
Commit | Line | Data |
---|---|---|
4d6aee34 | 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 | ***********************************************************************************************************/ | |
0e60774c | 11 | #include "TObject.h" |
0affc9e5 | 12 | #include "TObjArray.h" |
0e60774c | 13 | #include "TMath.h" |
14 | #include "TF1.h" | |
15 | #include "TH1F.h" | |
16 | #include "TGraph.h" | |
17 | #include "TGraphErrors.h" | |
0d83b3a5 | 18 | #include "TPDGCode.h" |
0e60774c | 19 | |
20 | #include "AliLog.h" | |
21 | #include "Cal/AliTRDCalPID.h" | |
0affc9e5 | 22 | #include "AliCDBManager.h" |
23 | #include "AliCDBEntry.h" | |
24 | #include "AliESDtrack.h" | |
25 | #include "AliPID.h" | |
0e60774c | 26 | #include "AliTRDpidUtil.h" |
27 | ||
28 | ClassImp(AliTRDpidUtil) | |
29 | ||
30 | ||
4d6aee34 | 31 | Float_t AliTRDpidUtil::fgEleEffi = 0.9; |
0e60774c | 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 | ||
0e60774c | 46 | //________________________________________________________________________ |
c46a7947 | 47 | Bool_t AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2) |
0e60774c | 48 | // Double_t AliTRDpidUtil::GetError() |
49 | { | |
50 | // | |
51 | // Calculates the pion efficiency | |
52 | // | |
53 | ||
54 | histo1 -> SetLineColor(kRed); | |
55 | histo2 -> SetLineColor(kBlue); | |
c46a7947 | 56 | if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){ |
57 | AliWarning("Histo with no entries !"); | |
58 | return kFALSE; | |
0e60774c | 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 | |
4d6aee34 | 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)); | |
0e60774c | 69 | |
70 | ||
71 | // calculate electron efficiency of each bin | |
422a2dc0 | 72 | for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){ |
73 | aBinSumE = 0; | |
0e60774c | 74 | aBinSumE = histo1 -> GetBinContent(abinE); |
75 | ||
4d6aee34 | 76 | sumElecE[abinE] = sumElecE[abinE+1] + aBinSumE; |
422a2dc0 | 77 | |
4d6aee34 | 78 | if((sumElecE[abinE] >= fgEleEffi) && (bFirst == 1)){ |
0e60774c | 79 | bFirst = 0; |
80 | cbinE = abinE; | |
4d6aee34 | 81 | fCalcEleEffi = (sumElecE[cbinE]); |
0e60774c | 82 | } |
83 | } | |
84 | ||
85 | fThreshold = histo1 -> GetBinCenter(cbinE); | |
86 | ||
87 | // calculate pion efficiency of each bin | |
422a2dc0 | 88 | for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){ |
0e60774c | 89 | bBinSumE = 0; |
90 | bBinSumE = histo2 -> GetBinContent(bbinE); | |
91 | ||
4d6aee34 | 92 | sumPionsE[bbinE] = sumPionsE[bbinE+1] + bBinSumE; |
0e60774c | 93 | if(bbinE == cbinE){ |
4d6aee34 | 94 | fPionEffi = (sumPionsE[cbinE]); |
0e60774c | 95 | } |
96 | } | |
97 | ||
98 | ||
99 | // pion efficiency vs electron efficiency | |
4d6aee34 | 100 | TGraph gEffis(kBins, sumElecE, sumPionsE); |
0e60774c | 101 | |
102 | // use fit function to get derivate of the TGraph for error calculation | |
4d6aee34 | 103 | TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05); |
104 | gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.05); | |
422a2dc0 | 105 | |
0e60774c | 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!"); | |
c46a7947 | 109 | return kFALSE; |
0e60774c | 110 | } |
4d6aee34 | 111 | fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi))); |
422a2dc0 | 112 | |
c46a7947 | 113 | // AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold)); |
4d6aee34 | 114 | // AliInfo(Form("Derivative at %4.2f : %f\n", fgEleEffi, f1.Derivative(fgEleEffi))); |
c46a7947 | 115 | return kTRUE; |
0e60774c | 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 | ||
0affc9e5 | 154 | |
155 | //__________________________________________________________________________ | |
0d83b3a5 | 156 | Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){ |
0affc9e5 | 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(); | |
4d6aee34 | 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); | |
0affc9e5 | 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 | //__________________________________________________________________________ | |
0d83b3a5 | 178 | Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){ |
0affc9e5 | 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(); | |
4d6aee34 | 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); | |
0affc9e5 | 191 | } |
0d83b3a5 | 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 |