]>
Commit | Line | Data |
---|---|---|
142229c3 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* $Id: AliTRDdigitizer.cxx 44182 2010-10-10 16:23:39Z cblume $ */ | |
17 | ||
18 | //////////////////////////////////////////////////////////////////////////// | |
19 | // // | |
20 | // Helper class for TRD PID efficiency calculation. // | |
21 | // Calculation of the hadron efficiency dependent on momentum and of // | |
22 | // the errors implemented in function CalculatePionEff. The pion // | |
23 | // efficiency is based on a predefined electron efficiency. // | |
24 | // The default is 90%. To change the, one has to call the function // | |
25 | // SetElectronEfficiency. // | |
26 | // Other Helper functions decide based on 90% electron efficiency // | |
27 | // whether a certain track is accepted as Electron Track. // | |
28 | // The reference data is stored in the TRD OCDB. // | |
29 | // // | |
30 | //////////////////////////////////////////////////////////////////////////// | |
31 | ||
0e60774c | 32 | #include "TObject.h" |
0affc9e5 | 33 | #include "TObjArray.h" |
0e60774c | 34 | #include "TMath.h" |
35 | #include "TF1.h" | |
36 | #include "TH1F.h" | |
37 | #include "TGraph.h" | |
38 | #include "TGraphErrors.h" | |
0d83b3a5 | 39 | #include "TPDGCode.h" |
0e60774c | 40 | |
41 | #include "AliLog.h" | |
42 | #include "Cal/AliTRDCalPID.h" | |
0affc9e5 | 43 | #include "AliCDBManager.h" |
44 | #include "AliCDBEntry.h" | |
45 | #include "AliESDtrack.h" | |
46 | #include "AliPID.h" | |
0e60774c | 47 | #include "AliTRDpidUtil.h" |
48 | ||
49 | ClassImp(AliTRDpidUtil) | |
50 | ||
4d6aee34 | 51 | Float_t AliTRDpidUtil::fgEleEffi = 0.9; |
0e60774c | 52 | |
53 | //________________________________________________________________________ | |
54 | AliTRDpidUtil::AliTRDpidUtil() | |
55 | : TObject() | |
56 | ,fCalcEleEffi(0.) | |
57 | ,fPionEffi(-1.) | |
58 | ,fError(-1.) | |
59 | ,fThreshold(-1.) | |
60 | { | |
61 | // | |
62 | // Default constructor | |
63 | // | |
64 | } | |
65 | ||
0e60774c | 66 | //________________________________________________________________________ |
c46a7947 | 67 | Bool_t AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2) |
0e60774c | 68 | // Double_t AliTRDpidUtil::GetError() |
69 | { | |
70 | // | |
71 | // Calculates the pion efficiency | |
72 | // | |
73 | ||
74 | histo1 -> SetLineColor(kRed); | |
75 | histo2 -> SetLineColor(kBlue); | |
d7aee360 | 76 | if(!histo1->GetEntries() || !histo2 -> GetEntries()){ |
77 | AliError("Probability histos empty !"); | |
78 | return kFALSE; | |
79 | } | |
80 | if(histo1->GetNbinsX() != histo2->GetNbinsX()){ | |
81 | AliError(Form("Electron probability discretized differently from pions [%d %d] !", histo1->GetNbinsX(), histo2->GetNbinsX())); | |
c46a7947 | 82 | return kFALSE; |
0e60774c | 83 | } |
84 | histo1 -> Scale(1./histo1->GetEntries()); | |
85 | histo2 -> Scale(1./histo2->GetEntries()); | |
86 | ||
d7aee360 | 87 | // array of the integrated sum in each bin |
88 | Double_t sumElecE[kBins+2], sumPionsE[kBins+2]; | |
4d6aee34 | 89 | memset(sumElecE, 0, (kBins+2)*sizeof(Double_t)); |
90 | memset(sumPionsE, 0, (kBins+2)*sizeof(Double_t)); | |
0e60774c | 91 | |
d7aee360 | 92 | Int_t nbinE(histo1->GetNbinsX()), |
93 | abinE(nbinE), | |
94 | bbinE(nbinE), | |
95 | cbinE(-1); | |
96 | // calculate electron efficiency for each bin | |
97 | // and also integral distribution | |
98 | for(Bool_t bFirst(kTRUE); abinE--;){ | |
99 | sumElecE[abinE] = sumElecE[abinE+1] + histo1->GetBinContent(abinE+1); | |
100 | if((sumElecE[abinE] >= fgEleEffi) && bFirst){ | |
0e60774c | 101 | cbinE = abinE; |
d7aee360 | 102 | fCalcEleEffi = sumElecE[cbinE]; |
103 | bFirst = kFALSE; | |
0e60774c | 104 | } |
105 | } | |
d7aee360 | 106 | fThreshold = histo1->GetBinCenter(cbinE); |
0e60774c | 107 | |
108 | // calculate pion efficiency of each bin | |
d7aee360 | 109 | // and also integral distribution |
110 | for (;bbinE--;){ | |
111 | sumPionsE[bbinE] = sumPionsE[bbinE+1] + histo2->GetBinContent(bbinE+1); | |
112 | if(bbinE == cbinE) fPionEffi = sumPionsE[cbinE]; | |
0e60774c | 113 | } |
114 | ||
115 | ||
116 | // pion efficiency vs electron efficiency | |
d7aee360 | 117 | TGraph gEffis(nbinE, sumElecE, sumPionsE); |
0e60774c | 118 | |
119 | // use fit function to get derivate of the TGraph for error calculation | |
4d6aee34 | 120 | TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05); |
121 | gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.05); | |
422a2dc0 | 122 | |
0e60774c | 123 | // return the error of the pion efficiency |
124 | if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){ | |
d7aee360 | 125 | AliError(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!"); |
c46a7947 | 126 | return kFALSE; |
0e60774c | 127 | } |
4d6aee34 | 128 | fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi))); |
422a2dc0 | 129 | |
d7aee360 | 130 | AliDebug(2, Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold)); |
131 | AliDebug(2, Form("Derivative at %4.2f : %f\n", fgEleEffi, f1.Derivative(fgEleEffi))); | |
c46a7947 | 132 | return kTRUE; |
0e60774c | 133 | } |
134 | ||
135 | ||
136 | //__________________________________________________________________________ | |
137 | Int_t AliTRDpidUtil::GetMomentumBin(Double_t p) | |
138 | { | |
139 | // | |
140 | // returns the momentum bin for a given momentum | |
141 | // | |
142 | ||
143 | Int_t pBin1 = -1; // check bin1 | |
144 | Int_t pBin2 = -1; // check bin2 | |
145 | ||
146 | if(p < 0) return -1; // return -1 if momentum < 0 | |
147 | if(p < AliTRDCalPID::GetMomentum(0)) return 0; // smallest momentum bin | |
148 | if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin | |
149 | ||
150 | ||
151 | // calculate momentum bin for non extremal momenta | |
152 | for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){ | |
153 | if(p < AliTRDCalPID::GetMomentum(iMomBin)){ | |
154 | pBin1 = iMomBin - 1; | |
155 | pBin2 = iMomBin; | |
156 | } | |
157 | else | |
158 | continue; | |
159 | ||
160 | if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){ | |
161 | return pBin2; | |
162 | } | |
163 | else{ | |
164 | return pBin1; | |
165 | } | |
166 | } | |
167 | ||
168 | return -1; | |
169 | } | |
170 | ||
0affc9e5 | 171 | |
172 | //__________________________________________________________________________ | |
0d83b3a5 | 173 | Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){ |
0affc9e5 | 174 | // |
175 | // Do PID decision for the TRD based on 90% Electron efficiency threshold | |
176 | // | |
177 | // Author: Markus Fasel (M.Fasel@gsi.de) | |
178 | // | |
179 | if(method == kESD) method = kNN; | |
180 | TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"}; | |
181 | AliCDBManager *cdb = AliCDBManager::Instance(); | |
4d6aee34 | 182 | AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds"); |
142229c3 | 183 | if (!cdbThresholds) return kFALSE; |
4d6aee34 | 184 | TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject()); |
142229c3 | 185 | if (!histos) return kFALSE; |
186 | TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data())); | |
187 | if (!thresholdHist) return kFALSE; | |
4d6aee34 | 188 | Double_t threshold = thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1); |
0affc9e5 | 189 | |
190 | // Do Decision | |
d8069611 | 191 | Double_t pidProbs[5]; |
192 | track->GetTRDpid(pidProbs); | |
193 | if(pidProbs[AliPID::kElectron] >= threshold) return kTRUE; | |
0affc9e5 | 194 | return kFALSE; |
195 | } | |
196 | ||
197 | //__________________________________________________________________________ | |
0d83b3a5 | 198 | Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){ |
0affc9e5 | 199 | // |
200 | // Returns the pion efficiency at 90% electron efficiency from the OCDB | |
201 | // | |
202 | // Author: Markus Fasel (M.Fasel@gsi.de) | |
203 | // | |
204 | if(method == kESD) method = kNN; | |
205 | TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"}; | |
206 | AliCDBManager *cdb = AliCDBManager::Instance(); | |
4d6aee34 | 207 | AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds"); |
142229c3 | 208 | if (!cdbThresholds) return kFALSE; |
4d6aee34 | 209 | TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject()); |
142229c3 | 210 | if (!histos) return kFALSE; |
211 | TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data())); | |
212 | if (!thresholdHist) return kFALSE; | |
4d6aee34 | 213 | return thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1); |
0affc9e5 | 214 | } |
0d83b3a5 | 215 | |
216 | //________________________________________________________________________ | |
217 | Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){ | |
218 | // | |
2f4384e6 | 219 | // Private Helper function to get the paticle species (ALICE notation) |
0d83b3a5 | 220 | // from the Pdg code |
221 | // | |
222 | Int_t species; | |
223 | switch(TMath::Abs(pdg)){ | |
224 | case kElectron: | |
225 | species = AliPID::kElectron; | |
226 | break; | |
227 | case kMuonMinus: | |
228 | species = AliPID::kMuon; | |
229 | break; | |
230 | case kPiPlus: | |
231 | species = AliPID::kPion; | |
232 | break; | |
233 | case kKPlus: | |
234 | species = AliPID::kKaon; | |
235 | break; | |
236 | case kProton: | |
237 | species = AliPID::kProton; | |
238 | break; | |
239 | default: | |
240 | species = -1; | |
241 | } | |
242 | return species; | |
243 | } | |
244 | ||
2f4384e6 | 245 | //________________________________________________________________________ |
246 | Int_t AliTRDpidUtil::Mass2Pid(Float_t m){ | |
247 | // | |
248 | // Private Helper function to get the paticle species (ALICE notation) | |
249 | // from the Pdg mass | |
250 | // | |
251 | ||
252 | for(Int_t is(0); is<AliPID::kSPECIES; is++) if(TMath::Abs(m-AliPID::ParticleMass(is))<1.e-4) return is; | |
253 | return -1; | |
254 | } | |
255 |