Commented unnecessary include AliLog
[u/mrichter/AliRoot.git] / TRD / AliTRDpidUtil.cxx
CommitLineData
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
49ClassImp(AliTRDpidUtil)
50
4d6aee34 51Float_t AliTRDpidUtil::fgEleEffi = 0.9;
0e60774c 52
53//________________________________________________________________________
54AliTRDpidUtil::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 67Bool_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//__________________________________________________________________________
137Int_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 173Bool_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 198Double_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//________________________________________________________________________
217Int_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//________________________________________________________________________
246Int_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