/[PAMELA software]/PamCut/Collections/EffRigCollection/EffRigCollection.cpp
ViewVC logotype

Contents of /PamCut/Collections/EffRigCollection/EffRigCollection.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (show annotations) (download)
Thu Oct 1 10:33:49 2009 UTC (15 years, 2 months ago) by pam-fi
Branch: MAIN
Changes since 1.3: +18 -12 lines
Output format changed: now text output comprehends also the limits of the bin and the graph is produced even if outRoot == false.

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 }

  ViewVC Help
Powered by ViewVC 1.1.23