--- PamCut/Collections/EffRigCollection/EffRigCollection.cpp 2009/08/12 14:18:37 1.1 +++ PamCut/Collections/EffRigCollection/EffRigCollection.cpp 2009/11/08 16:04:01 1.6 @@ -9,8 +9,33 @@ #include "EffRigCollection.h" -EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, bool absRig) : - EffCollection(collectionName, outFileBase), _bins(0), _selVector(0), _detVector(0), _outUp(0), _outDown(0) { +#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); +} + +EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, int errMethod, bool owns) : + EffCollection(collectionName, outFileBase, errMethod, owns), _bins(0), _selVector(0), _detVector(0), _outUp(0), _outDown(0) { ifstream rigBinListFile; rigBinListFile.open(rigBinsFile); @@ -35,12 +60,9 @@ // 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; - if (_absRig) { - rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity(); - } - else - rig = 1. / event->GetTrack(0)->GetTrkTrack()->GetDeflection(); + + float rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity(); + if (rig >= _bins[0]) { int i = 1; while (rig >= _bins[i] && i < (int) _bins.size()) { @@ -50,8 +72,10 @@ if (i < (int) (_selVector.size())) { _selVector[i]++; + _sel++; if (_detCollection.ApplyCut(event) == CUTOK) { _detVector[i]++; + _det++; _nGood++; return CUTOK; } @@ -71,25 +95,102 @@ void EffRigCollection::Finalize() { - // Correct the number of selected events by subtracting events outside the rigidity interval - _sel -= _outUp + _outDown; // 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 + "-eff-rig.txt").Data(), ios_base::out); + 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 << _detVector[i] << " " << _selVector[i] << " "; + outTextFile << setw(10) << _bins[i] << setw(10) << _bins[i + 1] << setw(10) << _detVector[i] << setw(10) + << _selVector[i]; if (_selVector[i] != 0) - outTextFile << (float) _detVector[i] / (float) _selVector[i] << "\n"; + outTextFile << setw(10) << eff[i] << setw(10) << errLow[i] << setw(10) << errHigh[i]; else - outTextFile << "0.\n"; + 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(); + } }