Calibrate Digits before putting to CaloCells
[u/mrichter/AliRoot.git] / TRD / AliTRDpidUtil.cxx
CommitLineData
0e60774c 1#include "TObject.h"
0affc9e5 2#include "TObjArray.h"
0e60774c 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"
0affc9e5 11#include "AliCDBManager.h"
12#include "AliCDBEntry.h"
13#include "AliESDtrack.h"
14#include "AliPID.h"
0e60774c 15#include "AliTRDpidUtil.h"
16
17ClassImp(AliTRDpidUtil)
18
19
20Float_t AliTRDpidUtil::fEleEffi = 0.9;
21
22//________________________________________________________________________
23AliTRDpidUtil::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//________________________________________________________________________
c46a7947 38Bool_t AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
0e60774c 39// Double_t AliTRDpidUtil::GetError()
40{
41 //
42 // Calculates the pion efficiency
43 //
44
45 histo1 -> SetLineColor(kRed);
46 histo2 -> SetLineColor(kBlue);
c46a7947 47 if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){
48 AliWarning("Histo with no entries !");
49 return kFALSE;
0e60774c 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
422a2dc0 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));
0e60774c 60
61
62 // calculate electron efficiency of each bin
422a2dc0 63 for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){
64 aBinSumE = 0;
0e60774c 65 aBinSumE = histo1 -> GetBinContent(abinE);
66
67 SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
422a2dc0 68
0e60774c 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
422a2dc0 79 for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
0e60774c 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
c46a7947 91 TGraph gEffis(kBins, SumElecsE, SumPionsE);
0e60774c 92
93 // use fit function to get derivate of the TGraph for error calculation
c46a7947 94 TF1 f1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
95 gEffis.Fit(&f1, "Q", "", fEleEffi-.05, fEleEffi+.05);
422a2dc0 96
0e60774c 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!");
c46a7947 100 return kFALSE;
0e60774c 101 }
c46a7947 102 fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fEleEffi))*(f1.Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
422a2dc0 103
c46a7947 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;
0e60774c 107}
108
109
110//__________________________________________________________________________
111Int_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
0affc9e5 145
146//__________________________________________________________________________
147Bool_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//__________________________________________________________________________
169Double_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}