Collections/EffCollection/EffCollection.cpp

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

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

#include "EffCollection.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);
}

EffCollection::EffCollection(const char *collectionName, TString outFileBase, int errMethod, bool owns) :
  VerboseCollection(collectionName, false), _selCollection("selCollection",owns), _detCollection("detCollection",
      owns), _outFileBase(outFileBase), _errMethod(errMethod), _det(0), _sel(0) {

}

void EffCollection::AddDetectorCut(PamCut *cut) {
  _detCollection.AddCut(cut);
}

void EffCollection::AddSelectionCut(PamCut *cut) {
  _selCollection.AddCut(cut);
}

void EffCollection::AddDetectorAction(CollectionAction *action) {
  _detCollection.AddAction(action);
}

void EffCollection::AddSelectionAction(CollectionAction *action) {
  _selCollection.AddAction(action);
}

void EffCollection::Setup(PamLevel2 *events){
  // Base class have a single vector for cuts and another for actions. Here the cuts and actions
  // are not contained inside these vectors but rather inside two SmartCollection object members.
  // So we must call their Setup().
  _selCollection.Setup(events);
  _detCollection.Setup(events);

  // We call also base class' Setup(), which will likely do nothing since _sel and _det are empty.
  VerboseCollection::Setup(events);

}

int EffCollection::ApplyCut(PamLevel2 *event) {

  _nEv++;
  if (_selCollection.ApplyCut(event) == CUTOK) {
    _sel++;
    if (_detCollection.ApplyCut(event) == CUTOK) {
      _det++;
      _nGood++;
      return CUTOK;
    }
  }

  return 0;
}

void EffCollection::Finalize() {
  // Let's add all the cuts to the vector of the collection before calling VerboseCollection::Finalize
  for (unsigned int i = 0; i < _selCollection.GetSize(); i++)
    _cuts.push_back(_selCollection.GetCut(i));
  for (unsigned int i = 0; i < _selCollection.GetActionsSize(); i++){
    _actions.push_back(_selCollection.GetAction(i));
    _actionsPositions.push_back(_selCollection.GetActionPosition(i));
  }
  for (unsigned int i = 0; i < _detCollection.GetSize(); i++)
    _cuts.push_back(_detCollection.GetCut(i));
  for (unsigned int i = 0; i < _detCollection.GetActionsSize(); i++){
    _actions.push_back(_detCollection.GetAction(i));
    _actionsPositions.push_back(_detCollection.GetActionPosition(i) + _selCollection.GetSize());
  }

  // Now all the cuts and actions are in place, and VerboseCollection can print its report and call Finalize() for
  // every cut and action (calling SmartCollection::Finalize().
  VerboseCollection::Finalize();

  // Compute the error
  Int_t sel = (Int_t) _sel;
  Int_t det = (Int_t) _det;
  Double_t eff, errLow, errHigh;
  if (_errMethod == EFFERR_ROOT) {
    if (sel < 8) {
      eff = 1.1;
      errLow = errHigh = 0.;
    }
    else {
      TH1I pass("pass", "pass", 1, 0., 1.);
      TH1I total("total", "total", 1, 0., 1.);
      pass.Fill(0.5, det);
      total.Fill(0.5, sel);
      TGraphAsymmErrors errGraph;
      errGraph.BayesDivide(&pass, &total);
      Double_t dummy;
      errGraph.GetPoint(0, dummy, eff);
      errLow = errGraph.GetErrorYlow(0);
      errHigh = errGraph.GetErrorYhigh(0);
    }
  }
  if (_errMethod == EFFERR_SERGIO) {
    efficiency_(&sel, &det, &eff, &errLow, &errHigh);
  }

  // Add some info about efficiency to the stdout
  cout << "    ****** Efficiency informations ******\n";
  cout << "    Detector cuts:\n";
  for (unsigned int i = 0; i < _detCollection.GetSize(); i++)
    cout << "      - " << _detCollection.GetCut(i)->GetName() << "\n";
  cout << "    Total detector efficiency: " << eff << " +" << errHigh << "-" << errLow << endl;

  // Write the output file
  if (_outFileBase != "") {

    ofstream outTextFile((_outFileBase + "-" + GetName() + ".txt").Data(), ios_base::out);
    streamsize newPrec = 4;
    outTextFile.precision(newPrec);
    outTextFile.setf(ios::fixed, ios::floatfield);
    outTextFile << setw(10) << det << setw(10) << sel << setw(10) << eff << setw(10) << errLow << setw(10) << errHigh
        << endl;
    outTextFile.close();
  }

}
1	pam-fi	1.1	/*
2			* EffCollection.cpp
3			*
4			* Created on: 10/ago/2009
5			* Author: Nicola Mori
6			*/
7
8			/! @file EffCollection.cpp The EffCollection class implementation file. /
9
10			#include "EffCollection.h"
11	pam-fi	1.6	#include "TGraphAsymmErrors.h"
12	pam-fi	1.1
13	pam-fi	1.5	extern "C" {
14	pam-fi	1.8	/*! @brief Fortran routine by Sergio Ricciarini which computes errors on the efficiency.
15			*
16			* This routine has been developed to carefully evaluate the error bars on the efficiency,
17			* which must have a value between 0 and 1. See the appendix of Sergio's PhD thesis for more
18			* details.
19			* The important thing is that the computed efficiency will have the unphysical value 1.1
20			* if there are less than 8 events in the efficiency sample (eg., the sample used to compute the
21			* efficiency). This is necessary to have a reliable estimate of the errors.
22			*
23			* @param sel Pointer to an Int_t variable containing the number of events in the efficiency sample.
24			* @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			* @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	pam-fi	1.5	}
35
36	pam-fi	1.8	EffCollection::EffCollection(const char *collectionName, TString outFileBase, int errMethod, bool owns) :
37	pam-fi	1.9	VerboseCollection(collectionName, false), _selCollection("selCollection",owns), _detCollection("detCollection",
38			owns), _outFileBase(outFileBase), _errMethod(errMethod), _det(0), _sel(0) {
39	pam-fi	1.1
40			}
41
42	pam-fi	1.8	void EffCollection::AddDetectorCut(PamCut *cut) {
43	pam-fi	1.1	_detCollection.AddCut(cut);
44			}
45
46	pam-fi	1.8	void EffCollection::AddSelectionCut(PamCut *cut) {
47	pam-fi	1.1	_selCollection.AddCut(cut);
48			}
49
50	pam-fi	1.8	void EffCollection::AddDetectorAction(CollectionAction *action) {
51	pam-fi	1.1	_detCollection.AddAction(action);
52			}
53
54	pam-fi	1.8	void EffCollection::AddSelectionAction(CollectionAction *action) {
55	pam-fi	1.1	_selCollection.AddAction(action);
56			}
57
58	pam-fi	1.10	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	pam-fi	1.11	// We call also base class' Setup(), which will likely do nothing since _sel and _det are empty.
66	pam-fi	1.10	VerboseCollection::Setup(events);
67
68			}
69
70	pam-fi	1.1	int EffCollection::ApplyCut(PamLevel2 *event) {
71
72	pam-fi	1.3	_nEv++;
73			if (_selCollection.ApplyCut(event) == CUTOK) {
74			_sel++;
75			if (_detCollection.ApplyCut(event) == CUTOK) {
76			_det++;
77			_nGood++;
78			return CUTOK;
79	pam-fi	1.1	}
80			}
81
82			return 0;
83			}
84
85			void EffCollection::Finalize() {
86			// 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++)
88			_cuts.push_back(_selCollection.GetCut(i));
89	pam-fi	1.12	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	pam-fi	1.1	for (unsigned int i = 0; i < _detCollection.GetSize(); i++)
94			_cuts.push_back(_detCollection.GetCut(i));
95	pam-fi	1.12	for (unsigned int i = 0; i < _detCollection.GetActionsSize(); i++){
96	pam-fi	1.10	_actions.push_back(_detCollection.GetAction(i));
97	pam-fi	1.12	_actionsPositions.push_back(_detCollection.GetActionPosition(i) + _selCollection.GetSize());
98			}
99	pam-fi	1.10
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	pam-fi	1.1	VerboseCollection::Finalize();
103	pam-fi	1.3
104	pam-fi	1.5	// Compute the error
105			Int_t sel = (Int_t) _sel;
106			Int_t det = (Int_t) _det;
107			Double_t eff, errLow, errHigh;
108	pam-fi	1.7	if (_errMethod == EFFERR_ROOT) {
109	pam-fi	1.6	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	pam-fi	1.7	if (_errMethod == EFFERR_SERGIO) {
127	pam-fi	1.6	efficiency_(&sel, &det, &eff, &errLow, &errHigh);
128			}
129	pam-fi	1.5
130			// Add some info about efficiency to the stdout
131	pam-fi	1.2	cout << " **** Efficiency informations ****\n";
132			cout << " Detector cuts:\n";
133			for (unsigned int i = 0; i < _detCollection.GetSize(); i++)
134			cout << " - " << _detCollection.GetCut(i)->GetName() << "\n";
135	pam-fi	1.5	cout << " Total detector efficiency: " << eff << " +" << errHigh << "-" << errLow << endl;
136	pam-fi	1.1
137	pam-fi	1.3	// Write the output file
138	pam-fi	1.1	if (_outFileBase != "") {
139	pam-fi	1.5
140	pam-fi	1.7	ofstream outTextFile((_outFileBase + "-" + GetName() + ".txt").Data(), ios_base::out);
141	pam-fi	1.5	streamsize newPrec = 4;
142			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	pam-fi	1.1	outTextFile.close();
147			}
148
149			}