More on char constness
[u/mrichter/AliRoot.git] / TRD / AliTRDpidUtil.cxx
CommitLineData
0e60774c 1#include "TObject.h"
2#include "TMath.h"
3#include "TF1.h"
4#include "TH1F.h"
5#include "TGraph.h"
6#include "TGraphErrors.h"
7
8#include "AliLog.h"
9#include "Cal/AliTRDCalPID.h"
10#include "AliTRDpidUtil.h"
11
12ClassImp(AliTRDpidUtil)
13
14
15Float_t AliTRDpidUtil::fEleEffi = 0.9;
16
17//________________________________________________________________________
18AliTRDpidUtil::AliTRDpidUtil()
19 : TObject()
20 ,fCalcEleEffi(0.)
21 ,fPionEffi(-1.)
22 ,fError(-1.)
23 ,fThreshold(-1.)
24{
25 //
26 // Default constructor
27 //
28}
29
30
31
32//________________________________________________________________________
422a2dc0 33void AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
0e60774c 34// Double_t AliTRDpidUtil::GetError()
35{
36 //
37 // Calculates the pion efficiency
38 //
39
40 histo1 -> SetLineColor(kRed);
41 histo2 -> SetLineColor(kBlue);
42 AliInfo(Form("Histo1[%d] Histo2[%d]", (Int_t)histo1 -> GetEntries(), (Int_t)histo2 -> GetEntries()));
43 if(!(histo1 -> GetEntries() > 0 && histo2 -> GetEntries() > 0)){
44 AliWarning("Histo has no Entries !");
45 return;
46 }
47 histo1 -> Scale(1./histo1->GetEntries());
48 histo2 -> Scale(1./histo2->GetEntries());
49
50 Int_t abinE, bbinE, cbinE = -1;
51 Double_t aBinSumE, bBinSumE; // content of a single bin
52 Bool_t bFirst = 1; // checks if threshold is crossed for the first time
422a2dc0 53 Double_t SumElecsE[kBins+2], SumPionsE[kBins+2]; // array of the integrated sum in each bin
54 memset(SumElecsE, 0, (kBins+2)*sizeof(Double_t));
55 memset(SumPionsE, 0, (kBins+2)*sizeof(Double_t));
0e60774c 56
57
58 // calculate electron efficiency of each bin
422a2dc0 59 for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){
60 aBinSumE = 0;
0e60774c 61 aBinSumE = histo1 -> GetBinContent(abinE);
62
63 SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
422a2dc0 64
0e60774c 65 if((SumElecsE[abinE] >= fEleEffi) && (bFirst == 1)){
66 bFirst = 0;
67 cbinE = abinE;
68 fCalcEleEffi = (SumElecsE[cbinE]);
69 }
70 }
71
72 fThreshold = histo1 -> GetBinCenter(cbinE);
73
74 // calculate pion efficiency of each bin
422a2dc0 75 for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
0e60774c 76 bBinSumE = 0;
77 bBinSumE = histo2 -> GetBinContent(bbinE);
78
79 SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
80 if(bbinE == cbinE){
81 fPionEffi = (SumPionsE[cbinE]);
82 }
83 }
84
85
86 // pion efficiency vs electron efficiency
87 TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE);
88
89 // use fit function to get derivate of the TGraph for error calculation
90 TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
91 gEffis -> Fit("f1","Q","",fEleEffi-.05, fEleEffi+.05);
422a2dc0 92
0e60774c 93 // return the error of the pion efficiency
94 if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
95 AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
96 return;
97 }
98 fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1 -> Derivative(fEleEffi))*(f1 -> Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
422a2dc0 99
100 AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
101 AliInfo(Form("Derivative at %4.2f : %f\n", fEleEffi, f1 -> Derivative(fEleEffi)));
102
0e60774c 103}
104
105
106//__________________________________________________________________________
107Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
108{
109 //
110 // returns the momentum bin for a given momentum
111 //
112
113 Int_t pBin1 = -1; // check bin1
114 Int_t pBin2 = -1; // check bin2
115
116 if(p < 0) return -1; // return -1 if momentum < 0
117 if(p < AliTRDCalPID::GetMomentum(0)) return 0; // smallest momentum bin
118 if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin
119
120
121 // calculate momentum bin for non extremal momenta
122 for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
123 if(p < AliTRDCalPID::GetMomentum(iMomBin)){
124 pBin1 = iMomBin - 1;
125 pBin2 = iMomBin;
126 }
127 else
128 continue;
129
130 if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
131 return pBin2;
132 }
133 else{
134 return pBin1;
135 }
136 }
137
138 return -1;
139}
140