Coding rules
[u/mrichter/AliRoot.git] / TRD / AliTRDpidUtil.cxx
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  ***********************************************************************************************************/
11 #include "TObject.h"
12 #include "TObjArray.h"
13 #include "TMath.h"
14 #include "TF1.h"
15 #include "TH1F.h"
16 #include "TGraph.h"
17 #include "TGraphErrors.h"
18 #include "TPDGCode.h"
19
20 #include "AliLog.h"
21 #include "Cal/AliTRDCalPID.h"
22 #include "AliCDBManager.h"
23 #include "AliCDBEntry.h"
24 #include "AliESDtrack.h"
25 #include "AliPID.h"
26 #include "AliTRDpidUtil.h"
27
28 ClassImp(AliTRDpidUtil)
29
30
31 Float_t AliTRDpidUtil::fgEleEffi = 0.9;
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
46 //________________________________________________________________________
47 Bool_t  AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
48 // Double_t  AliTRDpidUtil::GetError()
49 {
50   //
51   // Calculates the pion efficiency
52   //
53
54   histo1 -> SetLineColor(kRed);
55   histo2 -> SetLineColor(kBlue); 
56   if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){
57     AliWarning("Histo with no entries !");
58     return kFALSE;
59   }
60   histo1 -> Scale(1./histo1->GetEntries());
61   histo2 -> Scale(1./histo2->GetEntries());
62
63   Int_t abinE, bbinE, cbinE = -1;                    
64   Double_t aBinSumE, bBinSumE;                  // content of a single bin
65   Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
66   Double_t sumElecE[kBins+2], sumPionsE[kBins+2];  // array of the integrated sum in each bin
67   memset(sumElecE, 0, (kBins+2)*sizeof(Double_t));
68   memset(sumPionsE, 0, (kBins+2)*sizeof(Double_t));
69
70
71   // calculate electron efficiency of each bin
72   for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){  
73    aBinSumE = 0;
74     aBinSumE = histo1 -> GetBinContent(abinE);
75     
76     sumElecE[abinE] = sumElecE[abinE+1] + aBinSumE;
77
78     if((sumElecE[abinE] >= fgEleEffi) && (bFirst == 1)){
79       bFirst = 0;
80       cbinE = abinE;
81       fCalcEleEffi = (sumElecE[cbinE]); 
82     }
83   }
84   
85   fThreshold = histo1 -> GetBinCenter(cbinE);
86
87   // calculate pion efficiency of each bin
88   for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){       
89     bBinSumE = 0;
90     bBinSumE = histo2 -> GetBinContent(bbinE);
91
92     sumPionsE[bbinE] = sumPionsE[bbinE+1] + bBinSumE;
93     if(bbinE == cbinE){
94       fPionEffi = (sumPionsE[cbinE]);
95     }
96   }
97   
98
99   // pion efficiency vs electron efficiency
100   TGraph gEffis(kBins, sumElecE, sumPionsE);
101
102   // use fit function to get derivate of the TGraph for error calculation
103   TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05);
104   gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.05);
105   
106   // return the error of the pion efficiency
107   if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
108     AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
109     return kFALSE;
110   }
111   fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
112
113 //   AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
114 //   AliInfo(Form("Derivative at %4.2f : %f\n", fgEleEffi, f1.Derivative(fgEleEffi)));
115   return kTRUE;
116 }
117
118
119 //__________________________________________________________________________
120 Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
121 {
122   //
123   // returns the momentum bin for a given momentum
124   //
125
126   Int_t pBin1 = -1;                                    // check bin1
127   Int_t pBin2 = -1;                                    // check bin2
128
129   if(p < 0) return -1;                                 // return -1 if momentum < 0
130   if(p < AliTRDCalPID::GetMomentum(0)) return 0;                                      // smallest momentum bin
131   if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin
132
133
134   // calculate momentum bin for non extremal momenta
135   for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
136     if(p < AliTRDCalPID::GetMomentum(iMomBin)){
137       pBin1 = iMomBin - 1;
138       pBin2 = iMomBin;
139     }
140     else
141       continue;
142
143     if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
144        return pBin2;
145     }
146     else{
147       return pBin1;
148     }
149   }
150
151   return -1;
152 }
153
154
155 //__________________________________________________________________________
156 Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){
157   //
158   // Do PID decision for the TRD based on 90% Electron efficiency threshold
159   //
160   // Author: Markus Fasel (M.Fasel@gsi.de)
161   //
162   if(method == kESD) method = kNN;
163   TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"};
164   AliCDBManager *cdb = AliCDBManager::Instance(); 
165   AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
166   TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
167   TH1 * thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
168   Double_t threshold = thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
169   
170   // Do Decision
171   Double_t pid_probs[5];
172   track->GetTRDpid(pid_probs);
173   if(pid_probs[AliPID::kElectron] >= threshold) return kTRUE;
174   return kFALSE; 
175 }
176
177 //__________________________________________________________________________
178 Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){
179   //
180   // Returns the pion efficiency at 90% electron efficiency from the OCDB
181   //
182   // Author: Markus Fasel (M.Fasel@gsi.de)
183   //
184   if(method == kESD) method = kNN;
185   TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"};
186   AliCDBManager *cdb = AliCDBManager::Instance(); 
187   AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
188   TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
189   TH1 * thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
190   return thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
191 }
192
193 //________________________________________________________________________
194 Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
195   //
196   // Private Helper function to get the paticle species (ALICE notation) 
197   // from the Pdg code
198   //
199   Int_t species;
200   switch(TMath::Abs(pdg)){
201   case kElectron:
202     species = AliPID::kElectron;
203     break;
204   case kMuonMinus:
205     species = AliPID::kMuon;
206     break;
207   case kPiPlus:
208     species = AliPID::kPion;
209     break;
210   case kKPlus:
211     species = AliPID::kKaon;
212     break;
213   case kProton:
214     species = AliPID::kProton;
215     break;
216   default:
217     species = -1;
218   }
219   return species;
220 }
221