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