1 |
/* |
2 |
* EffRigCollection.cpp |
3 |
* |
4 |
* Created on: 10/ago/2009 |
5 |
* Author: Nicola Mori |
6 |
*/ |
7 |
|
8 |
/*! @file EffRigCollection.cpp The EffRigCollection class implementation file. */ |
9 |
|
10 |
#include "EffRigCollection.h" |
11 |
|
12 |
#include "TGraphAsymmErrors.h" |
13 |
|
14 |
extern "C" { |
15 |
bool efficiency_(Int_t*, Int_t*, Double_t*, Double_t*, Double_t*); |
16 |
} |
17 |
|
18 |
EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, int errMethod) : |
19 |
EffCollection(collectionName, outFileBase, errMethod), _bins(0), _selVector(0), _detVector(0), _outUp(0), _outDown(0) { |
20 |
|
21 |
ifstream rigBinListFile; |
22 |
rigBinListFile.open(rigBinsFile); |
23 |
|
24 |
TString auxString; |
25 |
while (!rigBinListFile.eof()) { |
26 |
rigBinListFile >> auxString; |
27 |
if (auxString != "") { |
28 |
_bins.push_back(auxString.Atof()); |
29 |
} |
30 |
} |
31 |
rigBinListFile.close(); |
32 |
|
33 |
_selVector.resize(_bins.size() - 1, 0); |
34 |
_detVector.resize(_bins.size() - 1, 0); |
35 |
} |
36 |
|
37 |
int EffRigCollection::ApplyCut(PamLevel2 *event) { |
38 |
|
39 |
_nEv++; |
40 |
if (_selCollection.ApplyCut(event) == CUTOK) { |
41 |
// Check if the event is inside the rigidity range |
42 |
// NOTE: at this point a TrkPhSinCut should be already performed, |
43 |
// since we are going to retrieve rigidity. |
44 |
|
45 |
float rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity(); |
46 |
|
47 |
if (rig >= _bins[0]) { |
48 |
int i = 1; |
49 |
while (rig >= _bins[i] && i < (int) _bins.size()) { |
50 |
i++; |
51 |
} |
52 |
i--; |
53 |
|
54 |
if (i < (int) (_selVector.size())) { |
55 |
_selVector[i]++; |
56 |
_sel++; |
57 |
if (_detCollection.ApplyCut(event) == CUTOK) { |
58 |
_detVector[i]++; |
59 |
_det++; |
60 |
_nGood++; |
61 |
return CUTOK; |
62 |
} |
63 |
} |
64 |
else |
65 |
_outUp++; |
66 |
|
67 |
} |
68 |
else { |
69 |
_outDown++; |
70 |
return 0; |
71 |
} |
72 |
} |
73 |
|
74 |
return 0; |
75 |
} |
76 |
|
77 |
void EffRigCollection::Finalize() { |
78 |
|
79 |
// Print the report |
80 |
EffCollection::Finalize(); |
81 |
cout << " Events below the minimum rigidity: " << _outDown << "\n"; |
82 |
cout << " Events above the maximum rigidity: " << _outUp << "\n"; |
83 |
|
84 |
// Compute the error |
85 |
Int_t sel[_selVector.size()]; |
86 |
Int_t det[_detVector.size()]; |
87 |
Double_t eff[_selVector.size()]; |
88 |
Double_t errLow[_selVector.size()]; |
89 |
Double_t errHigh[_selVector.size()]; |
90 |
for (unsigned int i = 0; i < _selVector.size(); i++) { |
91 |
sel[i] = (Int_t) _selVector[i]; |
92 |
det[i] = (Int_t) _detVector[i]; |
93 |
} |
94 |
|
95 |
TGraphAsymmErrors errGraph; |
96 |
errGraph.SetName(GetName()); |
97 |
errGraph.SetMarkerColor(kRed); |
98 |
errGraph.SetMarkerStyle(7); |
99 |
errGraph.GetYaxis()->SetRangeUser(0, 1.2); |
100 |
|
101 |
if (_errMethod == EFFERR_ROOT) { |
102 |
double binning[_bins.size()]; |
103 |
for (unsigned int i = 0; i < _bins.size(); i++) |
104 |
binning[i] = _bins[i]; |
105 |
|
106 |
TH1I pass("pass", "pass", _bins.size() - 1, binning); |
107 |
TH1I total("total", "total", _bins.size() - 1, binning); |
108 |
for (unsigned int i = 0; i < _selVector.size(); i++) { |
109 |
for (int j = 0; j < det[i]; j++) |
110 |
pass.Fill((binning[i + 1] + binning[i]) / 2.); |
111 |
for (int j = 0; j < sel[i]; j++) |
112 |
total.Fill((binning[i + 1] + binning[i]) / 2.); |
113 |
} |
114 |
|
115 |
errGraph.BayesDivide(&pass, &total); |
116 |
Double_t graphX; |
117 |
double currBin; |
118 |
int currPoint = 0; |
119 |
|
120 |
for (unsigned int i = 0; i < _selVector.size(); i++) { |
121 |
errGraph.GetPoint(currPoint, graphX, eff[i]); |
122 |
currBin = (binning[i + 1] + binning[i]) / 2.; |
123 |
if (currBin == graphX) { |
124 |
if (_selVector[i] < 8) { |
125 |
eff[i] = 1.1; |
126 |
errLow[i] = errHigh[i] = 0.; |
127 |
errGraph.SetPoint(currPoint, graphX, 1.1); |
128 |
float halfBin = (binning[i + 1] - binning[i]) / 2.; |
129 |
errGraph.SetPointError(currPoint, halfBin, halfBin, 0., 0.); |
130 |
} |
131 |
currPoint++; |
132 |
} |
133 |
} |
134 |
|
135 |
} |
136 |
if (_errMethod == EFFERR_SERGIO) { |
137 |
for (unsigned int i = 0; i < _selVector.size(); i++) { |
138 |
efficiency_(&(sel[i]), &(det[i]), &(eff[i]), &(errLow[i]), &(errHigh[i])); |
139 |
|
140 |
float R = (_bins[i] + _bins[i + 1]) / 2.; |
141 |
float halfBin = (_bins[i + 1] - _bins[i]) / 2.; |
142 |
errGraph.SetPoint(i, R, eff[i]); |
143 |
errGraph.SetPointError(i, halfBin, halfBin, errLow[i], errHigh[i]); |
144 |
|
145 |
} |
146 |
|
147 |
} |
148 |
|
149 |
// Write the output files |
150 |
if (_outFileBase != "") { |
151 |
ofstream outTextFile((_outFileBase + "-" + GetName() + "-rig.txt").Data(), ios_base::out); |
152 |
streamsize newPrec = 4; |
153 |
outTextFile.precision(newPrec); |
154 |
outTextFile.setf(ios::fixed, ios::floatfield); |
155 |
for (unsigned int i = 0; i < _selVector.size(); i++) { |
156 |
outTextFile << setw(10) << _bins[i] << setw(10) << _bins[i + 1] << setw(10) << _detVector[i] << setw(10) |
157 |
<< _selVector[i]; |
158 |
if (_selVector[i] != 0) |
159 |
outTextFile << setw(10) << eff[i] << setw(10) << errLow[i] << setw(10) << errHigh[i] << "\n"; |
160 |
else |
161 |
outTextFile << setw(10) << 0. << setw(10) << 0. << setw(10) << 0. << endl; |
162 |
|
163 |
} |
164 |
|
165 |
outTextFile.close(); |
166 |
|
167 |
TFile outRootFile((_outFileBase + "-" + GetName() + "-rig.root"), "RECREATE"); |
168 |
outRootFile.cd(); |
169 |
errGraph.Write(); |
170 |
outRootFile.Close(); |
171 |
|
172 |
} |
173 |
|
174 |
} |