/[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.3 - (show annotations) (download)
Fri Sep 25 15:02:02 2009 UTC (15 years, 2 months ago) by pam-fi
Branch: MAIN
Changes since 1.2: +64 -9 lines
Support for error computation using ROOT's TGraphAsymmErrors::BayesDivide() added.

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 == EFFRIG_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 cout << (binning[i + 1] + binning[i]) / 2. << " " << det[i] << " " << sel[i] << endl;
110 for (int j = 0; j < det[i]; j++)
111 pass.Fill((binning[i + 1] + binning[i]) / 2.);
112 for (int j = 0; j < sel[i]; j++)
113 total.Fill((binning[i + 1] + binning[i]) / 2.);
114 }
115
116 errGraph.BayesDivide(&pass, &total);
117 Double_t graphX;
118 double currBin;
119 int currPoint = 0;
120 //cout << errGraph.GetN() << " " << pass.GetNbinsX() << " " << total.GetNbinsX() << endl;
121 for (unsigned int i = 0; i < _selVector.size(); i++) {
122 errGraph.GetPoint(currPoint, graphX, eff[i]);
123 currBin = (binning[i + 1] + binning[i]) / 2.;
124 if (currBin == graphX) {
125 if (_selVector[i] < 8) {
126 eff[i] = 1.1;
127 errLow[i] = errHigh[i] = 0.;
128 errGraph.SetPoint(currPoint, graphX, 1.1);
129 float halfBin = (binning[i + 1] - binning[i]) / 2.;
130 errGraph.SetPointError(currPoint, halfBin, halfBin, 0., 0.);
131 }
132 currPoint++;
133 }
134 }
135
136 }
137 if (_errMethod == EFFRIG_SERGIO) {
138 for (unsigned int i = 0; i < _selVector.size(); i++) {
139 efficiency_(&(sel[i]), &(det[i]), &(eff[i]), &(errLow[i]), &(errHigh[i]));
140 }
141
142 }
143
144 // Write the output files
145 if (_outFileBase != "") {
146 ofstream outTextFile((_outFileBase + "-eff-rig.txt").Data(), ios_base::out);
147 streamsize newPrec = 4;
148 outTextFile.precision(newPrec);
149 outTextFile.setf(ios::fixed, ios::floatfield);
150 for (unsigned int i = 0; i < _selVector.size(); i++) {
151 outTextFile << setw(10) << _detVector[i] << setw(10) << _selVector[i];
152 if (_selVector[i] != 0)
153 outTextFile << setw(10) << eff[i] << setw(10) << errLow[i] << setw(10) << errHigh[i] << "\n";
154 else
155 outTextFile << setw(10) << 0. << setw(10) << 0. << setw(10) << 0. << endl;
156
157 }
158
159 outTextFile.close();
160 if (_errMethod == EFFRIG_ROOT) {
161 TFile outRootFile((_outFileBase + "-eff-rig.root"), "RECREATE");
162 outRootFile.cd();
163 errGraph.Write();
164 outRootFile.Close();
165 }
166 }
167
168 }

  ViewVC Help
Powered by ViewVC 1.1.23