--- PamCut/Collections/EffRigCollection/EffRigCollection.cpp 2009/09/24 18:17:26 1.2 +++ PamCut/Collections/EffRigCollection/EffRigCollection.cpp 2010/03/26 16:29:46 1.8 @@ -9,12 +9,35 @@ #include "EffRigCollection.h" +#include "TGraphAsymmErrors.h" + extern "C" { -bool efficiency_(Int_t*, Int_t*, Double_t*, Double_t*, Double_t*); +/*! @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, bool absRig) : - EffCollection(collectionName, outFileBase), _bins(0), _selVector(0), _detVector(0), _outUp(0), _outDown(0) { +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); @@ -39,12 +62,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 = 1./event->GetTrack(0)->GetTrkTrack()->GetDeflection(); + if (rig >= _bins[0]) { int i = 1; while (rig >= _bins[i] && i < (int) _bins.size()) { @@ -91,23 +111,87 @@ for (unsigned int i = 0; i < _selVector.size(); i++) { sel[i] = (Int_t) _selVector[i]; det[i] = (Int_t) _detVector[i]; - efficiency_(&(sel[i]), &(det[i]), &(eff[i]), &(errLow[i]), &(errHigh[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 << setw(10) << _detVector[i] << setw(10) << _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 << setw(10) << eff[i] << setw(10) << errLow[i] << setw(10) << errHigh[i] << "\n"; + outTextFile << setw(10) << eff[i] << " " << setw(10) << errLow[i] << " " << setw(10) << errHigh[i]; else - outTextFile << setw(10) << 0. << setw(10) << 0. << setw(10) << 0. << endl; + 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(); + } }