Line data Source code
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 :
32 : #include "TObject.h"
33 : #include "TObjArray.h"
34 : #include "TMath.h"
35 : #include "TF1.h"
36 : #include "TH1F.h"
37 : #include "TGraph.h"
38 : #include "TGraphErrors.h"
39 : #include "TPDGCode.h"
40 :
41 : #include "AliLog.h"
42 : #include "AliTRDCalPID.h"
43 : #include "AliCDBManager.h"
44 : #include "AliCDBEntry.h"
45 : #include "AliESDtrack.h"
46 : #include "AliPID.h"
47 : #include "AliTRDpidUtil.h"
48 :
49 48 : ClassImp(AliTRDpidUtil)
50 :
51 : Float_t AliTRDpidUtil::fgEleEffi = 0.9;
52 :
53 : //________________________________________________________________________
54 : AliTRDpidUtil::AliTRDpidUtil()
55 0 : : TObject()
56 0 : ,fCalcEleEffi(0.)
57 0 : ,fPionEffi(-1.)
58 0 : ,fError(-1.)
59 0 : ,fThreshold(-1.)
60 0 : {
61 : //
62 : // Default constructor
63 : //
64 0 : }
65 :
66 : //________________________________________________________________________
67 : Bool_t AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
68 : // Double_t AliTRDpidUtil::GetError()
69 : {
70 : //
71 : // Calculates the pion efficiency
72 : //
73 :
74 0 : histo1 -> SetLineColor(kRed);
75 0 : histo2 -> SetLineColor(kBlue);
76 0 : if(!histo1->GetEntries() || !histo2 -> GetEntries()){
77 0 : AliError("Probability histos empty !");
78 0 : return kFALSE;
79 : }
80 0 : if(histo1->GetNbinsX() != histo2->GetNbinsX()){
81 0 : AliError(Form("Electron probability discretized differently from pions [%d %d] !", histo1->GetNbinsX(), histo2->GetNbinsX()));
82 0 : return kFALSE;
83 : }
84 0 : histo1 -> Scale(1./histo1->GetEntries());
85 0 : histo2 -> Scale(1./histo2->GetEntries());
86 :
87 : // array of the integrated sum in each bin
88 0 : Double_t sumElecE[kBins+2], sumPionsE[kBins+2];
89 0 : memset(sumElecE, 0, (kBins+2)*sizeof(Double_t));
90 0 : memset(sumPionsE, 0, (kBins+2)*sizeof(Double_t));
91 :
92 0 : 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 0 : for(Bool_t bFirst(kTRUE); abinE--;){
99 0 : sumElecE[abinE] = sumElecE[abinE+1] + histo1->GetBinContent(abinE+1);
100 0 : if((sumElecE[abinE] >= fgEleEffi) && bFirst){
101 : cbinE = abinE;
102 0 : fCalcEleEffi = sumElecE[cbinE];
103 : bFirst = kFALSE;
104 0 : }
105 : }
106 0 : fThreshold = histo1->GetBinCenter(cbinE);
107 :
108 : // calculate pion efficiency of each bin
109 : // and also integral distribution
110 0 : for (;bbinE--;){
111 0 : sumPionsE[bbinE] = sumPionsE[bbinE+1] + histo2->GetBinContent(bbinE+1);
112 0 : if(bbinE == cbinE) fPionEffi = sumPionsE[cbinE];
113 : }
114 :
115 :
116 : // pion efficiency vs electron efficiency
117 0 : TGraph gEffis(nbinE, sumElecE, sumPionsE);
118 :
119 : // use fit function to get derivate of the TGraph for error calculation
120 0 : TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05);
121 0 : gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.05);
122 :
123 : // return the error of the pion efficiency
124 0 : if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
125 0 : AliError(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
126 0 : return kFALSE;
127 : }
128 0 : fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
129 :
130 0 : AliDebug(2, Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
131 0 : AliDebug(2, Form("Derivative at %4.2f : %f\n", fgEleEffi, f1.Derivative(fgEleEffi)));
132 0 : return kTRUE;
133 0 : }
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 0 : if(p < 0) return -1; // return -1 if momentum < 0
147 0 : if(p < AliTRDCalPID::GetMomentum(0)) return 0; // smallest momentum bin
148 0 : 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 0 : for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
153 0 : if(p < AliTRDCalPID::GetMomentum(iMomBin)){
154 0 : pBin1 = iMomBin - 1;
155 : pBin2 = iMomBin;
156 : }
157 : else
158 : continue;
159 :
160 0 : if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
161 0 : return pBin2;
162 : }
163 : else{
164 0 : return pBin1;
165 : }
166 : }
167 :
168 0 : return -1;
169 0 : }
170 :
171 :
172 : //__________________________________________________________________________
173 : Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){
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 0 : if(method == kESD) method = kNN;
180 0 : TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"};
181 0 : AliCDBManager *cdb = AliCDBManager::Instance();
182 0 : AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
183 0 : if (!cdbThresholds) return kFALSE;
184 0 : TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
185 0 : if (!histos) return kFALSE;
186 0 : TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
187 0 : if (!thresholdHist) return kFALSE;
188 0 : Double_t threshold = thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
189 :
190 : // Do Decision
191 0 : Double_t pidProbs[5];
192 0 : track->GetTRDpid(pidProbs);
193 0 : if(pidProbs[AliPID::kElectron] >= threshold) return kTRUE;
194 0 : return kFALSE;
195 0 : }
196 :
197 : //__________________________________________________________________________
198 : Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){
199 : //
200 : // Returns the pion efficiency at 90% electron efficiency from the OCDB
201 : //
202 : // Author: Markus Fasel (M.Fasel@gsi.de)
203 : //
204 0 : if(method == kESD) method = kNN;
205 0 : TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"};
206 0 : AliCDBManager *cdb = AliCDBManager::Instance();
207 0 : AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
208 0 : if (!cdbThresholds) return kFALSE;
209 0 : TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
210 0 : if (!histos) return kFALSE;
211 0 : TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
212 0 : if (!thresholdHist) return kFALSE;
213 0 : return thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
214 0 : }
215 :
216 : //________________________________________________________________________
217 : Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
218 : //
219 : // Private Helper function to get the paticle species (ALICE notation)
220 : // from the Pdg code
221 : //
222 : Int_t species;
223 0 : switch(TMath::Abs(pdg)){
224 : case kElectron:
225 : species = AliPID::kElectron;
226 0 : break;
227 : case kMuonMinus:
228 : species = AliPID::kMuon;
229 0 : break;
230 : case kPiPlus:
231 : species = AliPID::kPion;
232 0 : break;
233 : case kKPlus:
234 : species = AliPID::kKaon;
235 0 : break;
236 : case kProton:
237 : species = AliPID::kProton;
238 0 : break;
239 : default:
240 : species = -1;
241 0 : }
242 0 : return species;
243 : }
244 :
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 0 : for(Int_t is(0); is<AliPID::kSPECIES; is++) if(TMath::Abs(m-AliPID::ParticleMass(is))<1.e-4) return is;
253 0 : return -1;
254 0 : }
255 :
|
