]>
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); | |
d7aee360 | 56 | if(!histo1->GetEntries() || !histo2 -> GetEntries()){ |
57 | AliError("Probability histos empty !"); | |
58 | return kFALSE; | |
59 | } | |
60 | if(histo1->GetNbinsX() != histo2->GetNbinsX()){ | |
61 | AliError(Form("Electron probability discretized differently from pions [%d %d] !", histo1->GetNbinsX(), histo2->GetNbinsX())); | |
c46a7947 | 62 | return kFALSE; |
0e60774c | 63 | } |
64 | histo1 -> Scale(1./histo1->GetEntries()); | |
65 | histo2 -> Scale(1./histo2->GetEntries()); | |
66 | ||
d7aee360 | 67 | // array of the integrated sum in each bin |
68 | Double_t sumElecE[kBins+2], sumPionsE[kBins+2]; | |
4d6aee34 | 69 | memset(sumElecE, 0, (kBins+2)*sizeof(Double_t)); |
70 | memset(sumPionsE, 0, (kBins+2)*sizeof(Double_t)); | |
0e60774c | 71 | |
d7aee360 | 72 | Int_t nbinE(histo1->GetNbinsX()), |
73 | abinE(nbinE), | |
74 | bbinE(nbinE), | |
75 | cbinE(-1); | |
76 | // calculate electron efficiency for each bin | |
77 | // and also integral distribution | |
78 | for(Bool_t bFirst(kTRUE); abinE--;){ | |
79 | sumElecE[abinE] = sumElecE[abinE+1] + histo1->GetBinContent(abinE+1); | |
80 | if((sumElecE[abinE] >= fgEleEffi) && bFirst){ | |
0e60774c | 81 | cbinE = abinE; |
d7aee360 | 82 | fCalcEleEffi = sumElecE[cbinE]; |
83 | bFirst = kFALSE; | |
0e60774c | 84 | } |
85 | } | |
d7aee360 | 86 | fThreshold = histo1->GetBinCenter(cbinE); |
0e60774c | 87 | |
88 | // calculate pion efficiency of each bin | |
d7aee360 | 89 | // and also integral distribution |
90 | for (;bbinE--;){ | |
91 | sumPionsE[bbinE] = sumPionsE[bbinE+1] + histo2->GetBinContent(bbinE+1); | |
92 | if(bbinE == cbinE) fPionEffi = sumPionsE[cbinE]; | |
0e60774c | 93 | } |
94 | ||
95 | ||
96 | // pion efficiency vs electron efficiency | |
d7aee360 | 97 | TGraph gEffis(nbinE, sumElecE, sumPionsE); |
0e60774c | 98 | |
99 | // use fit function to get derivate of the TGraph for error calculation | |
4d6aee34 | 100 | TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05); |
101 | gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.05); | |
422a2dc0 | 102 | |
0e60774c | 103 | // return the error of the pion efficiency |
104 | if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){ | |
d7aee360 | 105 | AliError(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!"); |
c46a7947 | 106 | return kFALSE; |
0e60774c | 107 | } |
4d6aee34 | 108 | fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi))); |
422a2dc0 | 109 | |
d7aee360 | 110 | AliDebug(2, Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold)); |
111 | AliDebug(2, Form("Derivative at %4.2f : %f\n", fgEleEffi, f1.Derivative(fgEleEffi))); | |
c46a7947 | 112 | return kTRUE; |
0e60774c | 113 | } |
114 | ||
115 | ||
116 | //__________________________________________________________________________ | |
117 | Int_t AliTRDpidUtil::GetMomentumBin(Double_t p) | |
118 | { | |
119 | // | |
120 | // returns the momentum bin for a given momentum | |
121 | // | |
122 | ||
123 | Int_t pBin1 = -1; // check bin1 | |
124 | Int_t pBin2 = -1; // check bin2 | |
125 | ||
126 | if(p < 0) return -1; // return -1 if momentum < 0 | |
127 | if(p < AliTRDCalPID::GetMomentum(0)) return 0; // smallest momentum bin | |
128 | if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin | |
129 | ||
130 | ||
131 | // calculate momentum bin for non extremal momenta | |
132 | for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){ | |
133 | if(p < AliTRDCalPID::GetMomentum(iMomBin)){ | |
134 | pBin1 = iMomBin - 1; | |
135 | pBin2 = iMomBin; | |
136 | } | |
137 | else | |
138 | continue; | |
139 | ||
140 | if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){ | |
141 | return pBin2; | |
142 | } | |
143 | else{ | |
144 | return pBin1; | |
145 | } | |
146 | } | |
147 | ||
148 | return -1; | |
149 | } | |
150 | ||
0affc9e5 | 151 | |
152 | //__________________________________________________________________________ | |
0d83b3a5 | 153 | Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){ |
0affc9e5 | 154 | // |
155 | // Do PID decision for the TRD based on 90% Electron efficiency threshold | |
156 | // | |
157 | // Author: Markus Fasel (M.Fasel@gsi.de) | |
158 | // | |
159 | if(method == kESD) method = kNN; | |
160 | TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"}; | |
161 | AliCDBManager *cdb = AliCDBManager::Instance(); | |
4d6aee34 | 162 | AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds"); |
163 | TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject()); | |
164 | TH1 * thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data())); | |
165 | Double_t threshold = thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1); | |
0affc9e5 | 166 | |
167 | // Do Decision | |
168 | Double_t pid_probs[5]; | |
169 | track->GetTRDpid(pid_probs); | |
170 | if(pid_probs[AliPID::kElectron] >= threshold) return kTRUE; | |
171 | return kFALSE; | |
172 | } | |
173 | ||
174 | //__________________________________________________________________________ | |
0d83b3a5 | 175 | Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){ |
0affc9e5 | 176 | // |
177 | // Returns the pion efficiency at 90% electron efficiency from the OCDB | |
178 | // | |
179 | // Author: Markus Fasel (M.Fasel@gsi.de) | |
180 | // | |
181 | if(method == kESD) method = kNN; | |
182 | TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"}; | |
183 | AliCDBManager *cdb = AliCDBManager::Instance(); | |
4d6aee34 | 184 | AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds"); |
185 | TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject()); | |
186 | TH1 * thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data())); | |
187 | return thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1); | |
0affc9e5 | 188 | } |
0d83b3a5 | 189 | |
190 | //________________________________________________________________________ | |
191 | Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){ | |
192 | // | |
193 | // Private Helper function to get the paticle species (ALICE notation) | |
194 | // from the Pdg code | |
195 | // | |
196 | Int_t species; | |
197 | switch(TMath::Abs(pdg)){ | |
198 | case kElectron: | |
199 | species = AliPID::kElectron; | |
200 | break; | |
201 | case kMuonMinus: | |
202 | species = AliPID::kMuon; | |
203 | break; | |
204 | case kPiPlus: | |
205 | species = AliPID::kPion; | |
206 | break; | |
207 | case kKPlus: | |
208 | species = AliPID::kKaon; | |
209 | break; | |
210 | case kProton: | |
211 | species = AliPID::kProton; | |
212 | break; | |
213 | default: | |
214 | species = -1; | |
215 | } | |
216 | return species; | |
217 | } | |
218 |