Collections/EffRigCollection/EffRigCollection.cpp

/*
 * EffRigCollection.cpp
 *
 *  Created on: 10/ago/2009
 *      Author: Nicola Mori
 */

/*! @file EffRigCollection.cpp The EffRigCollection class implementation file. */

#include "EffRigCollection.h"

#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);

  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);
}

int EffRigCollection::ApplyCut(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 = 1./event->GetTrack(0)->GetTrkTrack()->GetDeflection();

    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();

  }

}
1	/*
2	* EffRigCollection.cpp
3	*
4	* Created on: 10/ago/2009
5	* Author: Nicola Mori
6	*/
7
8	/! @file EffRigCollection.cpp The EffRigCollection class implementation file. /
9
10	#include "EffRigCollection.h"
11
12	#include "TGraphAsymmErrors.h"
13
14	extern "C" {
15	/*! @brief Fortran routine by Sergio Ricciarini which computes errors on the efficiency.
16	*
17	* This routine has been developed to carefully evaluate the error bars on the efficiency,
18	* which must have a value between 0 and 1. See the appendix of Sergio's PhD thesis for more
19	* details.
20	* The important thing is that the computed efficiency will have the unphysical value 1.1
21	* if there are less than 8 events in the efficiency sample (eg., the sample used to compute the
22	* efficiency). This is necessary to have a reliable estimate of the errors.
23	*
24	* @param sel Pointer to an Int_t variable containing the number of events in the efficiency sample.
25	* @param det Pointer to an Int_t variable containing the number of events in the efficiency sample
26	* which satisfied the set of cuts for which the efficiency is being measured.
27	* @param eff Pointer to a Double_t variable where the computed efficiency will be returned.
28	* @param errLow Pointer to a Double_t variable where the length of the lower error bar will
29	* be returned.
30	* @param errHigh Pointer to a Double_t variable where the length of the upper error bar will
31	* be returned.
32	* @return Not clear at this point... ignore it.
33	*/
34	bool efficiency_(Int_t* sel, Int_t* det, Double_t* eff, Double_t* errLow, Double_t* errHigh);
35	}
36
37	EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, int errMethod,
38	bool owns) :
39	EffCollection(collectionName, outFileBase, errMethod, owns), _bins(0), _selVector(0), _detVector(0), _outUp(0),
40	_outDown(0) {
41
42	ifstream rigBinListFile;
43	rigBinListFile.open(rigBinsFile);
44
45	TString auxString;
46	while (!rigBinListFile.eof()) {
47	rigBinListFile >> auxString;
48	if (auxString != "") {
49	_bins.push_back(auxString.Atof());
50	}
51	}
52	rigBinListFile.close();
53
54	_selVector.resize(_bins.size() - 1, 0);
55	_detVector.resize(_bins.size() - 1, 0);
56	}
57
58	int EffRigCollection::ApplyCut(PamLevel2 *event) {
59
60	_nEv++;
61	if (_selCollection.ApplyCut(event) == CUTOK) {
62	// Check if the event is inside the rigidity range
63	// NOTE: at this point a TrkPhSinCut should be already performed,
64	// since we are going to retrieve rigidity.
65
66	float rig = 1./event->GetTrack(0)->GetTrkTrack()->GetDeflection();
67
68	if (rig >= _bins[0]) {
69	int i = 1;
70	while (rig >= _bins[i] && i < (int) _bins.size()) {
71	i++;
72	}
73	i--;
74
75	if (i < (int) (_selVector.size())) {
76	_selVector[i]++;
77	_sel++;
78	if (_detCollection.ApplyCut(event) == CUTOK) {
79	_detVector[i]++;
80	_det++;
81	_nGood++;
82	return CUTOK;
83	}
84	}
85	else
86	_outUp++;
87
88	}
89	else {
90	_outDown++;
91	return 0;
92	}
93	}
94
95	return 0;
96	}
97
98	void EffRigCollection::Finalize() {
99
100	// Print the report
101	EffCollection::Finalize();
102	cout << " Events below the minimum rigidity: " << _outDown << "\n";
103	cout << " Events above the maximum rigidity: " << _outUp << "\n";
104
105	// Compute the error
106	Int_t sel[_selVector.size()];
107	Int_t det[_detVector.size()];
108	Double_t eff[_selVector.size()];
109	Double_t errLow[_selVector.size()];
110	Double_t errHigh[_selVector.size()];
111	for (unsigned int i = 0; i < _selVector.size(); i++) {
112	sel[i] = (Int_t) _selVector[i];
113	det[i] = (Int_t) _detVector[i];
114	}
115
116	TGraphAsymmErrors errGraph;
117	errGraph.SetName(GetName());
118	errGraph.SetMarkerColor(kRed);
119	errGraph.SetMarkerStyle(7);
120	errGraph.GetYaxis()->SetRangeUser(0, 1.2);
121
122	if (_errMethod == EFFERR_ROOT) {
123	double binning[_bins.size()];
124	for (unsigned int i = 0; i < _bins.size(); i++)
125	binning[i] = _bins[i];
126
127	TH1I pass("pass", "pass", _bins.size() - 1, binning);
128	TH1I total("total", "total", _bins.size() - 1, binning);
129	for (unsigned int i = 0; i < _selVector.size(); i++) {
130	for (int j = 0; j < det[i]; j++)
131	pass.Fill((binning[i + 1] + binning[i]) / 2.);
132	for (int j = 0; j < sel[i]; j++)
133	total.Fill((binning[i + 1] + binning[i]) / 2.);
134	}
135
136	errGraph.BayesDivide(&pass, &total);
137	Double_t graphX;
138	double currBin;
139	int currPoint = 0;
140
141	for (unsigned int i = 0; i < _selVector.size(); i++) {
142	errGraph.GetPoint(currPoint, graphX, eff[i]);
143	currBin = (binning[i + 1] + binning[i]) / 2.;
144	if (currBin == graphX) {
145	if (_selVector[i] < 8) {
146	eff[i] = 1.1;
147	errLow[i] = errHigh[i] = 0.;
148	errGraph.SetPoint(currPoint, graphX, 1.1);
149	float halfBin = (binning[i + 1] - binning[i]) / 2.;
150	errGraph.SetPointError(currPoint, halfBin, halfBin, 0., 0.);
151	}
152	currPoint++;
153	}
154	}
155
156	}
157	if (_errMethod == EFFERR_SERGIO) {
158	for (unsigned int i = 0; i < _selVector.size(); i++) {
159	efficiency_(&(sel[i]), &(det[i]), &(eff[i]), &(errLow[i]), &(errHigh[i]));
160
161	float R = (_bins[i] + _bins[i + 1]) / 2.;
162	float halfBin = (_bins[i + 1] - _bins[i]) / 2.;
163	errGraph.SetPoint(i, R, eff[i]);
164	errGraph.SetPointError(i, halfBin, halfBin, errLow[i], errHigh[i]);
165
166	}
167
168	}
169
170	// Write the output files
171	if (_outFileBase != "") {
172	ofstream outTextFile((_outFileBase + "-" + GetName() + "-rig.txt").Data(), ios_base::out);
173	streamsize newPrec = 4;
174	outTextFile.precision(newPrec);
175	outTextFile.setf(ios::fixed, ios::floatfield);
176	for (unsigned int i = 0; i < _selVector.size(); i++) {
177	outTextFile << setw(10) << _bins[i] << " " << setw(10) << _bins[i + 1] << " " << setw(10) << _detVector[i] << " "
178	<< setw(10) << _selVector[i] << " ";
179	if (_selVector[i] != 0)
180	outTextFile << setw(10) << eff[i] << " " << setw(10) << errLow[i] << " " << setw(10) << errHigh[i];
181	else
182	outTextFile << setw(10) << 0. << " " << setw(10) << 0. << " " << setw(10) << 0.;
183	if (i < _selVector.size() - 1) //Avoids to print an empty line at the end of the file
184	outTextFile << endl;
185
186	}
187
188	outTextFile.close();
189
190	TFile outRootFile((_outFileBase + "-" + GetName() + "-rig.root"), "RECREATE");
191	outRootFile.cd();
192	errGraph.Write();
193	outRootFile.Close();
194
195	}
196
197	}