/[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.7 by pam-fi, Thu Mar 18 14:40:51 2010 UTC revision 1.9 by pam-fi, Tue Jul 12 17:36:03 2011 UTC
# Line 9  Line 9 
9    
10  #include "EffRigCollection.h"  #include "EffRigCollection.h"
11    
 #include "TGraphAsymmErrors.h"  
   
 extern "C" {  
 /*! @brief Fortran routine by Sergio Ricciarini which computes errors on the efficiency.  
  *  
  * This routine has been developed to carefully evaluate the error bars on the efficiency,  
  * which must have a value between 0 and 1. See the appendix of Sergio's PhD thesis for more  
  * details.  
  * The important thing is that the computed efficiency will have the unphysical value 1.1  
  * if there are less than 8 events in the efficiency sample (eg., the sample used to compute the  
  * efficiency). This is necessary to have a reliable estimate of the errors.  
  *  
  * @param sel Pointer to an Int_t variable containing the number of events in the efficiency sample.  
  * @param det Pointer to an Int_t variable containing the number of events in the efficiency sample  
  *            which satisfied the set of cuts  for which the efficiency is being measured.  
  * @param eff Pointer to a Double_t variable where the computed efficiency will be returned.  
  * @param errLow Pointer to a Double_t variable where the length of the lower error bar will  
  *               be returned.  
  * @param errHigh Pointer to a Double_t variable where the length of the upper error bar will  
  *               be returned.  
  * @return Not clear at this point... ignore it.  
  */  
 bool efficiency_(Int_t* sel, Int_t* det, Double_t* eff, Double_t* errLow, Double_t* errHigh);  
 }  
   
12  EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, int errMethod,  EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, int errMethod,
13      bool owns) :      bool owns) :
14    EffCollection(collectionName, outFileBase, errMethod, owns), _bins(0), _selVector(0), _detVector(0), _outUp(0),      BinnedEffCollection(collectionName, outFileBase, rigBinsFile, errMethod, owns) {
       _outDown(0) {  
   
   ifstream rigBinListFile;  
   rigBinListFile.open(rigBinsFile);  
   
   TString auxString;  
   while (!rigBinListFile.eof()) {  
     rigBinListFile >> auxString;  
     if (auxString != "") {  
       _bins.push_back(auxString.Atof());  
     }  
   }  
   rigBinListFile.close();  
   
   _selVector.resize(_bins.size() - 1, 0);  
   _detVector.resize(_bins.size() - 1, 0);  
15  }  }
16    
17  int EffRigCollection::ApplyCut(PamLevel2 *event) {  float EffRigCollection::GetBinValue(PamLevel2 *event) {
   
   _nEv++;  
   if (_selCollection.ApplyCut(event) == CUTOK) {  
     // Check if the event is inside the rigidity range  
     // NOTE: at this point a TrkPhSinCut should be already performed,  
     //       since we are going to retrieve rigidity.  
   
     float rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity();  
   
     if (rig >= _bins[0]) {  
       int i = 1;  
       while (rig >= _bins[i] && i < (int) _bins.size()) {  
         i++;  
       }  
       i--;  
   
       if (i < (int) (_selVector.size())) {  
         _selVector[i]++;  
         _sel++;  
         if (_detCollection.ApplyCut(event) == CUTOK) {  
           _detVector[i]++;  
           _det++;  
           _nGood++;  
           return CUTOK;  
         }  
       }  
       else  
         _outUp++;  
   
     }  
     else {  
       _outDown++;  
       return 0;  
     }  
   }  
   
   return 0;  
 }  
   
 void EffRigCollection::Finalize() {  
   
   // Print the report  
   EffCollection::Finalize();  
   cout << "    Events below the minimum rigidity: " << _outDown << "\n";  
   cout << "    Events above the maximum rigidity: " << _outUp << "\n";  
   
   // Compute the error  
   Int_t sel[_selVector.size()];  
   Int_t det[_detVector.size()];  
   Double_t eff[_selVector.size()];  
   Double_t errLow[_selVector.size()];  
   Double_t errHigh[_selVector.size()];  
   for (unsigned int i = 0; i < _selVector.size(); i++) {  
     sel[i] = (Int_t) _selVector[i];  
     det[i] = (Int_t) _detVector[i];  
   }  
   
   TGraphAsymmErrors errGraph;  
   errGraph.SetName(GetName());  
   errGraph.SetMarkerColor(kRed);  
   errGraph.SetMarkerStyle(7);  
   errGraph.GetYaxis()->SetRangeUser(0, 1.2);  
   
   if (_errMethod == EFFERR_ROOT) {  
     double binning[_bins.size()];  
     for (unsigned int i = 0; i < _bins.size(); i++)  
       binning[i] = _bins[i];  
   
     TH1I pass("pass", "pass", _bins.size() - 1, binning);  
     TH1I total("total", "total", _bins.size() - 1, binning);  
     for (unsigned int i = 0; i < _selVector.size(); i++) {  
       for (int j = 0; j < det[i]; j++)  
         pass.Fill((binning[i + 1] + binning[i]) / 2.);  
       for (int j = 0; j < sel[i]; j++)  
         total.Fill((binning[i + 1] + binning[i]) / 2.);  
     }  
   
     errGraph.BayesDivide(&pass, &total);  
     Double_t graphX;  
     double currBin;  
     int currPoint = 0;  
   
     for (unsigned int i = 0; i < _selVector.size(); i++) {  
       errGraph.GetPoint(currPoint, graphX, eff[i]);  
       currBin = (binning[i + 1] + binning[i]) / 2.;  
       if (currBin == graphX) {  
         if (_selVector[i] < 8) {  
           eff[i] = 1.1;  
           errLow[i] = errHigh[i] = 0.;  
           errGraph.SetPoint(currPoint, graphX, 1.1);  
           float halfBin = (binning[i + 1] - binning[i]) / 2.;  
           errGraph.SetPointError(currPoint, halfBin, halfBin, 0., 0.);  
         }  
         currPoint++;  
       }  
     }  
   
   }  
   if (_errMethod == EFFERR_SERGIO) {  
     for (unsigned int i = 0; i < _selVector.size(); i++) {  
       efficiency_(&(sel[i]), &(det[i]), &(eff[i]), &(errLow[i]), &(errHigh[i]));  
   
       float R = (_bins[i] + _bins[i + 1]) / 2.;  
       float halfBin = (_bins[i + 1] - _bins[i]) / 2.;  
       errGraph.SetPoint(i, R, eff[i]);  
       errGraph.SetPointError(i, halfBin, halfBin, errLow[i], errHigh[i]);  
   
     }  
   
   }  
   
   // Write the output files  
   if (_outFileBase != "") {  
     ofstream outTextFile((_outFileBase + "-" + GetName() + "-rig.txt").Data(), ios_base::out);  
     streamsize newPrec = 4;  
     outTextFile.precision(newPrec);  
     outTextFile.setf(ios::fixed, ios::floatfield);  
     for (unsigned int i = 0; i < _selVector.size(); i++) {  
       outTextFile << setw(10) << _bins[i] << " " << setw(10) << _bins[i + 1] << " " << setw(10) << _detVector[i] << " "  
           << setw(10) << _selVector[i] << " ";  
       if (_selVector[i] != 0)  
         outTextFile << setw(10) << eff[i] << " " << setw(10) << errLow[i] << " " << setw(10) << errHigh[i];  
       else  
         outTextFile << setw(10) << 0. << " " << setw(10) << 0. << " " << setw(10) << 0.;  
       if (i < _selVector.size() - 1) //Avoids to print an empty line at the end of the file  
         outTextFile << endl;  
   
     }  
   
     outTextFile.close();  
   
     TFile outRootFile((_outFileBase + "-" + GetName() + "-rig.root"), "RECREATE");  
     outRootFile.cd();  
     errGraph.Write();  
     outRootFile.Close();  
   
   }  
18    
19      return 1. / event->GetTrack(0)->GetTrkTrack()->GetDeflection();
20  }  }
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.23