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