]>
Commit | Line | Data |
---|---|---|
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 | ||
12 | ClassImp(AliTRDpidUtil) | |
13 | ||
14 | ||
15 | Float_t AliTRDpidUtil::fEleEffi = 0.9; | |
16 | ||
17 | //________________________________________________________________________ | |
18 | AliTRDpidUtil::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 | //________________________________________________________________________ | |
c46a7947 | 33 | Bool_t 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); | |
c46a7947 | 42 | if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){ |
43 | AliWarning("Histo with no entries !"); | |
44 | return kFALSE; | |
0e60774c | 45 | } |
46 | histo1 -> Scale(1./histo1->GetEntries()); | |
47 | histo2 -> Scale(1./histo2->GetEntries()); | |
48 | ||
49 | Int_t abinE, bbinE, cbinE = -1; | |
50 | Double_t aBinSumE, bBinSumE; // content of a single bin | |
51 | Bool_t bFirst = 1; // checks if threshold is crossed for the first time | |
422a2dc0 | 52 | Double_t SumElecsE[kBins+2], SumPionsE[kBins+2]; // array of the integrated sum in each bin |
53 | memset(SumElecsE, 0, (kBins+2)*sizeof(Double_t)); | |
54 | memset(SumPionsE, 0, (kBins+2)*sizeof(Double_t)); | |
0e60774c | 55 | |
56 | ||
57 | // calculate electron efficiency of each bin | |
422a2dc0 | 58 | for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){ |
59 | aBinSumE = 0; | |
0e60774c | 60 | aBinSumE = histo1 -> GetBinContent(abinE); |
61 | ||
62 | SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE; | |
422a2dc0 | 63 | |
0e60774c | 64 | if((SumElecsE[abinE] >= fEleEffi) && (bFirst == 1)){ |
65 | bFirst = 0; | |
66 | cbinE = abinE; | |
67 | fCalcEleEffi = (SumElecsE[cbinE]); | |
68 | } | |
69 | } | |
70 | ||
71 | fThreshold = histo1 -> GetBinCenter(cbinE); | |
72 | ||
73 | // calculate pion efficiency of each bin | |
422a2dc0 | 74 | for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){ |
0e60774c | 75 | bBinSumE = 0; |
76 | bBinSumE = histo2 -> GetBinContent(bbinE); | |
77 | ||
78 | SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE; | |
79 | if(bbinE == cbinE){ | |
80 | fPionEffi = (SumPionsE[cbinE]); | |
81 | } | |
82 | } | |
83 | ||
84 | ||
85 | // pion efficiency vs electron efficiency | |
c46a7947 | 86 | TGraph gEffis(kBins, SumElecsE, SumPionsE); |
0e60774c | 87 | |
88 | // use fit function to get derivate of the TGraph for error calculation | |
c46a7947 | 89 | TF1 f1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05); |
90 | gEffis.Fit(&f1, "Q", "", fEleEffi-.05, fEleEffi+.05); | |
422a2dc0 | 91 | |
0e60774c | 92 | // return the error of the pion efficiency |
93 | if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){ | |
94 | AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!"); | |
c46a7947 | 95 | return kFALSE; |
0e60774c | 96 | } |
c46a7947 | 97 | fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fEleEffi))*(f1.Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi))); |
422a2dc0 | 98 | |
c46a7947 | 99 | // AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold)); |
100 | // AliInfo(Form("Derivative at %4.2f : %f\n", fEleEffi, f1.Derivative(fEleEffi))); | |
101 | return kTRUE; | |
0e60774c | 102 | } |
103 | ||
104 | ||
105 | //__________________________________________________________________________ | |
106 | Int_t AliTRDpidUtil::GetMomentumBin(Double_t p) | |
107 | { | |
108 | // | |
109 | // returns the momentum bin for a given momentum | |
110 | // | |
111 | ||
112 | Int_t pBin1 = -1; // check bin1 | |
113 | Int_t pBin2 = -1; // check bin2 | |
114 | ||
115 | if(p < 0) return -1; // return -1 if momentum < 0 | |
116 | if(p < AliTRDCalPID::GetMomentum(0)) return 0; // smallest momentum bin | |
117 | if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin | |
118 | ||
119 | ||
120 | // calculate momentum bin for non extremal momenta | |
121 | for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){ | |
122 | if(p < AliTRDCalPID::GetMomentum(iMomBin)){ | |
123 | pBin1 = iMomBin - 1; | |
124 | pBin2 = iMomBin; | |
125 | } | |
126 | else | |
127 | continue; | |
128 | ||
129 | if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){ | |
130 | return pBin2; | |
131 | } | |
132 | else{ | |
133 | return pBin1; | |
134 | } | |
135 | } | |
136 | ||
137 | return -1; | |
138 | } | |
139 |