/[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.3 by pam-fi, Fri Sep 25 15:02:02 2009 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    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;    ifstream rigBinListFile;
22    rigBinListFile.open(rigBinsFile);    rigBinListFile.open(rigBinsFile);
# Line 35  int EffRigCollection::ApplyCut(PamLevel2 Line 41  int EffRigCollection::ApplyCut(PamLevel2
41      // Check if the event is inside the rigidity range      // Check if the event is inside the rigidity range
42      // NOTE: at this point a TrkPhSinCut should be already performed,      // NOTE: at this point a TrkPhSinCut should be already performed,
43      //       since we are going to retrieve rigidity.      //       since we are going to retrieve rigidity.
44      float rig;  
45      if (_absRig) {      float rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity();
46        rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity();  
     }  
     else  
       rig = 1. / event->GetTrack(0)->GetTrkTrack()->GetDeflection();  
47      if (rig >= _bins[0]) {      if (rig >= _bins[0]) {
48        int i = 1;        int i = 1;
49        while (rig >= _bins[i] && i < (int) _bins.size()) {        while (rig >= _bins[i] && i < (int) _bins.size()) {
# Line 50  int EffRigCollection::ApplyCut(PamLevel2 Line 53  int EffRigCollection::ApplyCut(PamLevel2
53    
54        if (i < (int) (_selVector.size())) {        if (i < (int) (_selVector.size())) {
55          _selVector[i]++;          _selVector[i]++;
56            _sel++;
57          if (_detCollection.ApplyCut(event) == CUTOK) {          if (_detCollection.ApplyCut(event) == CUTOK) {
58            _detVector[i]++;            _detVector[i]++;
59              _det++;
60            _nGood++;            _nGood++;
61            return CUTOK;            return CUTOK;
62          }          }
# Line 71  int EffRigCollection::ApplyCut(PamLevel2 Line 76  int EffRigCollection::ApplyCut(PamLevel2
76    
77  void EffRigCollection::Finalize() {  void EffRigCollection::Finalize() {
78    
   // Correct the number of selected events by subtracting events outside the rigidity interval  
   _sel -= _outUp + _outDown;  
79    // Print the report    // Print the report
80    EffCollection::Finalize();    EffCollection::Finalize();
81    cout << "    Events below the minimum rigidity: " << _outDown << "\n";    cout << "    Events below the minimum rigidity: " << _outDown << "\n";
82    cout << "    Events above the maximum rigidity: " << _outUp << "\n";    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    // Write the output files
145    if (_outFileBase != "") {    if (_outFileBase != "") {
146      ofstream outTextFile((_outFileBase + "-eff-rig.txt").Data(), ios_base::out);      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++) {      for (unsigned int i = 0; i < _selVector.size(); i++) {
151        outTextFile << _detVector[i] << "    " << _selVector[i] << "    ";        outTextFile << setw(10) << _detVector[i] << setw(10) << _selVector[i];
152        if (_selVector[i] != 0)        if (_selVector[i] != 0)
153          outTextFile << (float) _detVector[i] / (float) _selVector[i] << "\n";          outTextFile << setw(10) << eff[i] << setw(10) << errLow[i] << setw(10) << errHigh[i] << "\n";
154        else        else
155          outTextFile << "0.\n";          outTextFile << setw(10) << 0. << setw(10) << 0. << setw(10) << 0. << endl;
156    
157      }      }
158    
159      outTextFile.close();      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  }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.23