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

Diff of /PamCut/Collections/EffRigCollection/EffRigCollection.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by pam-fi, Wed Aug 12 14:18:37 2009 UTC revision 1.7 by pam-fi, Thu Mar 18 14:40:51 2010 UTC
# Line 9  Line 9 
9    
10  #include "EffRigCollection.h"  #include "EffRigCollection.h"
11    
12  EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, bool absRig) :  #include "TGraphAsymmErrors.h"
13    EffCollection(collectionName, outFileBase), _bins(0), _selVector(0), _detVector(0), _outUp(0), _outDown(0) {  
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    EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, 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 rigBinListFile;    ifstream rigBinListFile;
43    rigBinListFile.open(rigBinsFile);    rigBinListFile.open(rigBinsFile);
# Line 35  int EffRigCollection::ApplyCut(PamLevel2 Line 62  int EffRigCollection::ApplyCut(PamLevel2
62      // Check if the event is inside the rigidity range      // Check if the event is inside the rigidity range
63      // NOTE: at this point a TrkPhSinCut should be already performed,      // NOTE: at this point a TrkPhSinCut should be already performed,
64      //       since we are going to retrieve rigidity.      //       since we are going to retrieve rigidity.
65      float rig;  
66      if (_absRig) {      float rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity();
67        rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity();  
     }  
     else  
       rig = 1. / event->GetTrack(0)->GetTrkTrack()->GetDeflection();  
68      if (rig >= _bins[0]) {      if (rig >= _bins[0]) {
69        int i = 1;        int i = 1;
70        while (rig >= _bins[i] && i < (int) _bins.size()) {        while (rig >= _bins[i] && i < (int) _bins.size()) {
# Line 50  int EffRigCollection::ApplyCut(PamLevel2 Line 74  int EffRigCollection::ApplyCut(PamLevel2
74    
75        if (i < (int) (_selVector.size())) {        if (i < (int) (_selVector.size())) {
76          _selVector[i]++;          _selVector[i]++;
77            _sel++;
78          if (_detCollection.ApplyCut(event) == CUTOK) {          if (_detCollection.ApplyCut(event) == CUTOK) {
79            _detVector[i]++;            _detVector[i]++;
80              _det++;
81            _nGood++;            _nGood++;
82            return CUTOK;            return CUTOK;
83          }          }
# Line 71  int EffRigCollection::ApplyCut(PamLevel2 Line 97  int EffRigCollection::ApplyCut(PamLevel2
97    
98  void EffRigCollection::Finalize() {  void EffRigCollection::Finalize() {
99    
   // Correct the number of selected events by subtracting events outside the rigidity interval  
   _sel -= _outUp + _outDown;  
100    // Print the report    // Print the report
101    EffCollection::Finalize();    EffCollection::Finalize();
102    cout << "    Events below the minimum rigidity: " << _outDown << "\n";    cout << "    Events below the minimum rigidity: " << _outDown << "\n";
103    cout << "    Events above the maximum rigidity: " << _outUp << "\n";    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    // Write the output files
171    if (_outFileBase != "") {    if (_outFileBase != "") {
172      ofstream outTextFile((_outFileBase + "-eff-rig.txt").Data(), ios_base::out);      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++) {      for (unsigned int i = 0; i < _selVector.size(); i++) {
177        outTextFile << _detVector[i] << "    " << _selVector[i] << "    ";        outTextFile << setw(10) << _bins[i] << " " << setw(10) << _bins[i + 1] << " " << setw(10) << _detVector[i] << " "
178              << setw(10) << _selVector[i] << " ";
179        if (_selVector[i] != 0)        if (_selVector[i] != 0)
180          outTextFile << (float) _detVector[i] / (float) _selVector[i] << "\n";          outTextFile << setw(10) << eff[i] << " " << setw(10) << errLow[i] << " " << setw(10) << errHigh[i];
181        else        else
182          outTextFile << "0.\n";          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();      outTextFile.close();
189    
190        TFile outRootFile((_outFileBase + "-" + GetName() + "-rig.root"), "RECREATE");
191        outRootFile.cd();
192        errGraph.Write();
193        outRootFile.Close();
194    
195    }    }
196    
197  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.23