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

Annotation of /PamCut/Collections/BinnedEffCollection/BinnedEffCollection.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide annotations) (download)
Tue Jul 12 17:35:27 2011 UTC (13 years, 5 months ago) by pam-fi
Branch: MAIN
Added to repository. Mother class for binned efficiencies (reworked from EffRigCollection).

1 pam-fi 1.1 /*
2     * BinnedEffCollection.cpp
3     *
4     * Created on: 10/ago/2009
5     * Author: Nicola Mori
6     */
7    
8     /*! @file BinnedEffCollection.cpp The BinnedEffCollection class implementation file. */
9    
10     #include "BinnedEffCollection.h"
11    
12     #include "TGraphAsymmErrors.h"
13    
14     extern "C" {
15     /*! @brief Fortran routine by Sergio Ricciarini which computes errors on the efficiency.
16     *
17     * This routine has been developed to carefully evaluate the error bars on the efficiency,
18     * which must have a value between 0 and 1. See the appendix of Sergio's PhD thesis for more
19     * details.
20     * The important thing is that the computed efficiency will have the unphysical value 1.1
21     * if there are less than 8 events in the efficiency sample (eg., the sample used to compute the
22     * efficiency). This is necessary to have a reliable estimate of the errors.
23     *
24     * @param sel Pointer to an Int_t variable containing the number of events in the efficiency sample.
25     * @param det Pointer to an Int_t variable containing the number of events in the efficiency sample
26     * which satisfied the set of cuts for which the efficiency is being measured.
27     * @param eff Pointer to a Double_t variable where the computed efficiency will be returned.
28     * @param errLow Pointer to a Double_t variable where the length of the lower error bar will
29     * be returned.
30     * @param errHigh Pointer to a Double_t variable where the length of the upper error bar will
31     * be returned.
32     * @return Not clear at this point... ignore it.
33     */
34     bool efficiency_(Int_t* sel, Int_t* det, Double_t* eff, Double_t* errLow, Double_t* errHigh);
35     }
36    
37     BinnedEffCollection::BinnedEffCollection(const char *collectionName, TString outFileBase, TString binsFile, int errMethod,
38     bool owns) :
39     EffCollection(collectionName, outFileBase, errMethod, owns), _bins(0), _selVector(0), _detVector(0), _outUp(0),
40     _outDown(0) {
41    
42     ifstream binListFile;
43     binListFile.open(binsFile);
44    
45     TString auxString;
46     while (!binListFile.eof()) {
47     binListFile >> auxString;
48     if (auxString != "") {
49     _bins.push_back(auxString.Atof());
50     }
51     }
52     binListFile.close();
53    
54     _selVector.resize(_bins.size() - 1, 0);
55     _detVector.resize(_bins.size() - 1, 0);
56     }
57    
58     int BinnedEffCollection::ApplyCut(PamLevel2 *event) {
59    
60     _nEv++;
61     if (_selCollection.ApplyCut(event) == CUTOK) {
62     // Check if the event is inside the rigidity range
63     // NOTE: at this point a TrkPhSinCut should be already performed,
64     // since we are going to retrieve rigidity.
65    
66     float binVal = GetBinValue(event);
67    
68     if (binVal >= _bins[0]) {
69     int i = 1;
70     while (binVal >= _bins[i] && i < (int) _bins.size()) {
71     i++;
72     }
73     i--;
74    
75     if (i < (int) (_selVector.size())) {
76     _selVector[i]++;
77     _sel++;
78     if (_detCollection.ApplyCut(event) == CUTOK) {
79     _detVector[i]++;
80     _det++;
81     _nGood++;
82     return CUTOK;
83     }
84     }
85     else
86     _outUp++;
87    
88     }
89     else {
90     _outDown++;
91     return 0;
92     }
93     }
94    
95     return 0;
96     }
97    
98     void BinnedEffCollection::Finalize() {
99    
100     // Print the report
101     EffCollection::Finalize();
102     cout << " Events below the minimum rigidity: " << _outDown << "\n";
103     cout << " Events above the maximum rigidity: " << _outUp << "\n";
104    
105     // Compute the error
106     Int_t sel[_selVector.size()];
107     Int_t det[_detVector.size()];
108     Double_t eff[_selVector.size()];
109     Double_t errLow[_selVector.size()];
110     Double_t errHigh[_selVector.size()];
111     for (unsigned int i = 0; i < _selVector.size(); i++) {
112     sel[i] = (Int_t) _selVector[i];
113     det[i] = (Int_t) _detVector[i];
114     }
115    
116     TGraphAsymmErrors errGraph;
117     errGraph.SetName(GetName());
118     errGraph.SetMarkerColor(kRed);
119     errGraph.SetMarkerStyle(7);
120     errGraph.GetYaxis()->SetRangeUser(0, 1.2);
121    
122     if (_errMethod == EFFERR_ROOT) {
123     double binning[_bins.size()];
124     for (unsigned int i = 0; i < _bins.size(); i++)
125     binning[i] = _bins[i];
126    
127     TH1I pass("pass", "pass", _bins.size() - 1, binning);
128     TH1I total("total", "total", _bins.size() - 1, binning);
129     for (unsigned int i = 0; i < _selVector.size(); i++) {
130     for (int j = 0; j < det[i]; j++)
131     pass.Fill((binning[i + 1] + binning[i]) / 2.);
132     for (int j = 0; j < sel[i]; j++)
133     total.Fill((binning[i + 1] + binning[i]) / 2.);
134     }
135    
136     errGraph.BayesDivide(&pass, &total);
137     Double_t graphX;
138     double currBin;
139     int currPoint = 0;
140    
141     for (unsigned int i = 0; i < _selVector.size(); i++) {
142     errGraph.GetPoint(currPoint, graphX, eff[i]);
143     currBin = (binning[i + 1] + binning[i]) / 2.;
144     if (currBin == graphX) {
145     if (_selVector[i] < 8) {
146     eff[i] = 1.1;
147     errLow[i] = errHigh[i] = 0.;
148     errGraph.SetPoint(currPoint, graphX, 1.1);
149     float halfBin = (binning[i + 1] - binning[i]) / 2.;
150     errGraph.SetPointError(currPoint, halfBin, halfBin, 0., 0.);
151     }
152     currPoint++;
153     }
154     }
155    
156     }
157     if (_errMethod == EFFERR_SERGIO) {
158     for (unsigned int i = 0; i < _selVector.size(); i++) {
159     efficiency_(&(sel[i]), &(det[i]), &(eff[i]), &(errLow[i]), &(errHigh[i]));
160    
161     float R = (_bins[i] + _bins[i + 1]) / 2.;
162     float halfBin = (_bins[i + 1] - _bins[i]) / 2.;
163     errGraph.SetPoint(i, R, eff[i]);
164     errGraph.SetPointError(i, halfBin, halfBin, errLow[i], errHigh[i]);
165    
166     }
167    
168     }
169    
170     // Write the output files
171     if (_outFileBase != "") {
172     ofstream outTextFile((_outFileBase + "-" + GetName() + "-rig.txt").Data(), ios_base::out);
173     streamsize newPrec = 4;
174     outTextFile.precision(newPrec);
175     outTextFile.setf(ios::fixed, ios::floatfield);
176     for (unsigned int i = 0; i < _selVector.size(); i++) {
177     outTextFile << setw(10) << _bins[i] << " " << setw(10) << _bins[i + 1] << " " << setw(10) << _detVector[i] << " "
178     << setw(10) << _selVector[i] << " ";
179     if (_selVector[i] != 0)
180     outTextFile << setw(10) << eff[i] << " " << setw(10) << errLow[i] << " " << setw(10) << errHigh[i];
181     else
182     outTextFile << setw(10) << 0. << " " << setw(10) << 0. << " " << setw(10) << 0.;
183     if (i < _selVector.size() - 1) //Avoids to print an empty line at the end of the file
184     outTextFile << endl;
185    
186     }
187    
188     outTextFile.close();
189    
190     TFile outRootFile((_outFileBase + "-" + GetName() + "-rig.root"), "RECREATE");
191     outRootFile.cd();
192     errGraph.Write();
193     outRootFile.Close();
194    
195     }
196    
197     }

  ViewVC Help
Powered by ViewVC 1.1.23