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

Diff of /PamCut/Collections/EffCollection/EffCollection.cpp

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

revision 1.2 by pam-fi, Tue Aug 11 13:30:16 2009 UTC revision 1.12 by pam-fi, Mon Sep 5 12:59:52 2011 UTC
# Line 8  Line 8 
8  /*! @file EffCollection.cpp The EffCollection class implementation file. */  /*! @file EffCollection.cpp The EffCollection class implementation file. */
9    
10  #include "EffCollection.h"  #include "EffCollection.h"
11    #include "TGraphAsymmErrors.h"
12    
13  EffCollection::EffCollection(const char *collectionName, TString rigBinsFile, TString &outFileBase, bool absRig) :  extern "C" {
14    VerboseCollection(collectionName), _selCollection("selCollection"), _detCollection("detCollection"), _outFileBase(  /*! @brief Fortran routine by Sergio Ricciarini which computes errors on the efficiency.
15        outFileBase), _absRig(absRig), _bins(0), _sel(0), _det(0), _outUp(0), _outDown(0) {   *
16     * This routine has been developed to carefully evaluate the error bars on the efficiency,
17    ifstream rigBinListFile;   * which must have a value between 0 and 1. See the appendix of Sergio's PhD thesis for more
18    rigBinListFile.open(rigBinsFile);   * details.
19     * The important thing is that the computed efficiency will have the unphysical value 1.1
20    TString auxString;   * if there are less than 8 events in the efficiency sample (eg., the sample used to compute the
21    while (!rigBinListFile.eof()) {   * efficiency). This is necessary to have a reliable estimate of the errors.
22      rigBinListFile >> auxString;   *
23      if (auxString != "") {   * @param sel Pointer to an Int_t variable containing the number of events in the efficiency sample.
24        _bins.push_back(auxString.Atof());   * @param det Pointer to an Int_t variable containing the number of events in the efficiency sample
25      }   *            which satisfied the set of cuts  for which the efficiency is being measured.
26    }   * @param eff Pointer to a Double_t variable where the computed efficiency will be returned.
27    rigBinListFile.close();   * @param errLow Pointer to a Double_t variable where the length of the lower error bar will
28     *               be returned.
29     * @param errHigh Pointer to a Double_t variable where the length of the upper error bar will
30     *               be returned.
31     * @return Not clear at this point... ignore it.
32     */
33    bool efficiency_(Int_t* sel, Int_t* det, Double_t* eff, Double_t* errLow, Double_t* errHigh);
34    }
35    
36    EffCollection::EffCollection(const char *collectionName, TString outFileBase, int errMethod, bool owns) :
37      VerboseCollection(collectionName, false), _selCollection("selCollection",owns), _detCollection("detCollection",
38          owns), _outFileBase(outFileBase), _errMethod(errMethod), _det(0), _sel(0) {
39    
   _sel.resize(_bins.size() - 1, 0);  
   _det.resize(_bins.size() - 1, 0);  
40  }  }
41    
42  void EffCollection::AddDetectorCut(PamCut &cut) {  void EffCollection::AddDetectorCut(PamCut *cut) {
43    _detCollection.AddCut(cut);    _detCollection.AddCut(cut);
44  }  }
45    
46  void EffCollection::AddSelectionCut(PamCut &cut) {  void EffCollection::AddSelectionCut(PamCut *cut) {
47    _selCollection.AddCut(cut);    _selCollection.AddCut(cut);
48  }  }
49    
50  void EffCollection::AddDetectorAction(CollectionAction &action) {  void EffCollection::AddDetectorAction(CollectionAction *action) {
51    _detCollection.AddAction(action);    _detCollection.AddAction(action);
52  }  }
53    
54  void EffCollection::AddSelectionAction(CollectionAction &action) {  void EffCollection::AddSelectionAction(CollectionAction *action) {
55    _selCollection.AddAction(action);    _selCollection.AddAction(action);
56  }  }
57    
58    void EffCollection::Setup(PamLevel2 *events){
59      // Base class have a single vector for cuts and another for actions. Here the cuts and actions
60      // are not contained inside these vectors but rather inside two SmartCollection object members.
61      // So we must call their Setup().
62      _selCollection.Setup(events);
63      _detCollection.Setup(events);
64    
65      // We call also base class' Setup(), which will likely do nothing since _sel and _det are empty.
66      VerboseCollection::Setup(events);
67    
68    }
69    
70  int EffCollection::ApplyCut(PamLevel2 *event) {  int EffCollection::ApplyCut(PamLevel2 *event) {
71    
72    // See if the rigidity of the event is between the limits    _nEv++;
73    float rig;    if (_selCollection.ApplyCut(event) == CUTOK) {
74    if (_absRig)      _sel++;
75      rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity();      if (_detCollection.ApplyCut(event) == CUTOK) {
76    else        _det++;
77      rig = 1. / event->GetTrack(0)->GetTrkTrack()->GetDeflection();        _nGood++;
78    if (rig >= _bins[0]) {        return CUTOK;
     int i = 1;  
     while (rig >= _bins[i] && i < (int) _bins.size()) {  
       i++;  
79      }      }
     i--;  
     if (i < (int) (_sel.size())) {  
       _nEv++;  
       // Rigidity is OK, let's apply the cuts  
       if (_selCollection.ApplyCut(event) == CUTOK) {  
         _sel[i]++;  
         if (_detCollection.ApplyCut(event) == CUTOK) {  
           _det[i]++;  
           _nGood++;  
           return CUTOK;  
         }  
       }  
     }  
     else  
       _outUp++;  
80    }    }
   else  
     _outDown++;  
81    
82    return 0;    return 0;
83  }  }
# Line 84  void EffCollection::Finalize() { Line 86  void EffCollection::Finalize() {
86    // Let's add all the cuts to the vector of the collection before calling VerboseCollection::Finalize    // Let's add all the cuts to the vector of the collection before calling VerboseCollection::Finalize
87    for (unsigned int i = 0; i < _selCollection.GetSize(); i++)    for (unsigned int i = 0; i < _selCollection.GetSize(); i++)
88      _cuts.push_back(_selCollection.GetCut(i));      _cuts.push_back(_selCollection.GetCut(i));
89      for (unsigned int i = 0; i < _selCollection.GetActionsSize(); i++){
90        _actions.push_back(_selCollection.GetAction(i));
91        _actionsPositions.push_back(_selCollection.GetActionPosition(i));
92      }
93    for (unsigned int i = 0; i < _detCollection.GetSize(); i++)    for (unsigned int i = 0; i < _detCollection.GetSize(); i++)
94      _cuts.push_back(_detCollection.GetCut(i));      _cuts.push_back(_detCollection.GetCut(i));
95    // Now all the cuts are in place, and VerboseCollection can print its report    for (unsigned int i = 0; i < _detCollection.GetActionsSize(); i++){
96        _actions.push_back(_detCollection.GetAction(i));
97        _actionsPositions.push_back(_detCollection.GetActionPosition(i) + _selCollection.GetSize());
98      }
99    
100      // Now all the cuts and actions are in place, and VerboseCollection can print its report and call Finalize() for
101      // every cut and action (calling SmartCollection::Finalize().
102    VerboseCollection::Finalize();    VerboseCollection::Finalize();
103    // Add some info about efficiency  
104      // Compute the error
105      Int_t sel = (Int_t) _sel;
106      Int_t det = (Int_t) _det;
107      Double_t eff, errLow, errHigh;
108      if (_errMethod == EFFERR_ROOT) {
109        if (sel < 8) {
110          eff = 1.1;
111          errLow = errHigh = 0.;
112        }
113        else {
114          TH1I pass("pass", "pass", 1, 0., 1.);
115          TH1I total("total", "total", 1, 0., 1.);
116          pass.Fill(0.5, det);
117          total.Fill(0.5, sel);
118          TGraphAsymmErrors errGraph;
119          errGraph.BayesDivide(&pass, &total);
120          Double_t dummy;
121          errGraph.GetPoint(0, dummy, eff);
122          errLow = errGraph.GetErrorYlow(0);
123          errHigh = errGraph.GetErrorYhigh(0);
124        }
125      }
126      if (_errMethod == EFFERR_SERGIO) {
127        efficiency_(&sel, &det, &eff, &errLow, &errHigh);
128      }
129    
130      // Add some info about efficiency to the stdout
131    cout << "    ****** Efficiency informations ******\n";    cout << "    ****** Efficiency informations ******\n";
132    cout << "    Detector cuts:\n";    cout << "    Detector cuts:\n";
133    for (unsigned int i = 0; i < _detCollection.GetSize(); i++)    for (unsigned int i = 0; i < _detCollection.GetSize(); i++)
134      cout << "      - " << _detCollection.GetCut(i)->GetName() << "\n";      cout << "      - " << _detCollection.GetCut(i)->GetName() << "\n";
135    cout << "    Events below the minimum rigidity: " << _outDown << "\n";    cout << "    Total detector efficiency: " << eff << " +" << errHigh << "-" << errLow << endl;
   cout << "    Events above the maximum rigidity: " << _outUp << "\n";  
   cout << "    Total detector efficiency: "  
        << (float) _detCollection.GetNGood() /  
           (float) (_selCollection.GetCut(_selCollection.GetSize() - 1)->GetNGood())  
        << "\n" << endl;  
136    
137    // Write the output files    // Write the output file
138    if (_outFileBase != "") {    if (_outFileBase != "") {
     ofstream outTextFile((_outFileBase + "-sel.txt").Data(), ios_base::out);  
     for (unsigned int i = 0; i < _sel.size(); i++) {  
       outTextFile << _sel[i] << "\n";  
     }  
     outTextFile.close();  
139    
140      outTextFile.open((_outFileBase + "-det.txt").Data(), ios_base::out);      ofstream outTextFile((_outFileBase + "-" + GetName() + ".txt").Data(), ios_base::out);
141      for (unsigned int i = 0; i < _det.size(); i++) {      streamsize newPrec = 4;
142        outTextFile << _det[i] << "\n";      outTextFile.precision(newPrec);
143      }      outTextFile.setf(ios::fixed, ios::floatfield);
144        outTextFile << setw(10) << det << setw(10) << sel << setw(10) << eff << setw(10) << errLow << setw(10) << errHigh
145            << endl;
146      outTextFile.close();      outTextFile.close();
   
     outTextFile.open((_outFileBase + "-eff.txt").Data(), ios_base::out);  
     for (unsigned int i = 0; i < _det.size(); i++) {  
       if (_sel[i] != 0)  
         outTextFile << (float) _det[i] / (float) _sel[i] << "\n";  
       else  
         outTextFile << "0.\n";  
     }  
     outTextFile.close();  
   
147    }    }
148    
149  }  }
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23