Be sure to load mapping when needed
[u/mrichter/AliRoot.git] / TRD / AliTRDpidUtil.cxx
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 "Cal/AliTRDCalPID.h"
43 #include "AliCDBManager.h"
44 #include "AliCDBEntry.h"
45 #include "AliESDtrack.h"
46 #include "AliPID.h"
47 #include "AliTRDpidUtil.h"
48
49 ClassImp(AliTRDpidUtil)
50
51 Float_t AliTRDpidUtil::fgEleEffi = 0.9;
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
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   histo1 -> SetLineColor(kRed);
75   histo2 -> SetLineColor(kBlue); 
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()));
82     return kFALSE;
83   }
84   histo1 -> Scale(1./histo1->GetEntries());
85   histo2 -> Scale(1./histo2->GetEntries());
86
87   // array of the integrated sum in each bin
88   Double_t sumElecE[kBins+2], sumPionsE[kBins+2];  
89   memset(sumElecE, 0, (kBins+2)*sizeof(Double_t));
90   memset(sumPionsE, 0, (kBins+2)*sizeof(Double_t));
91
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){
101       cbinE = abinE;
102       fCalcEleEffi = sumElecE[cbinE];
103       bFirst = kFALSE;
104     }
105   }
106   fThreshold = histo1->GetBinCenter(cbinE);
107
108   // calculate pion efficiency of each bin
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];
113   }
114   
115
116   // pion efficiency vs electron efficiency
117   TGraph gEffis(nbinE, sumElecE, sumPionsE);
118
119   // use fit function to get derivate of the TGraph for error calculation
120   TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05);
121   gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.05);
122   
123   // return the error of the pion efficiency
124   if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
125     AliError(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
126     return kFALSE;
127   }
128   fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
129
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)));
132   return kTRUE;
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
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   if(method == kESD) method = kNN;
180   TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"};
181   AliCDBManager *cdb = AliCDBManager::Instance(); 
182   AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
183   if (!cdbThresholds) return kFALSE;
184   TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
185   if (!histos) return kFALSE;
186   TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
187   if (!thresholdHist) return kFALSE;
188   Double_t threshold = thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
189   
190   // Do Decision
191   Double_t pidProbs[5];
192   track->GetTRDpid(pidProbs);
193   if(pidProbs[AliPID::kElectron] >= threshold) return kTRUE;
194   return kFALSE; 
195 }
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   if(method == kESD) method = kNN;
205   TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"};
206   AliCDBManager *cdb = AliCDBManager::Instance(); 
207   AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
208   if (!cdbThresholds) return kFALSE;
209   TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
210   if (!histos) return kFALSE;
211   TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
212   if (!thresholdHist) return kFALSE;
213   return thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
214 }
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   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
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