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" {
bool efficiency_(Int_t*, Int_t*, Double_t*, Double_t*, Double_t*);
}

EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, int errMethod) :
  EffCollection(collectionName, outFileBase, errMethod), _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 = 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 == EFFRIG_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++) {
      cout << (binning[i + 1] + binning[i]) / 2. << "  " << det[i] << "  " << sel[i] << endl;
      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;
    //cout << errGraph.GetN() << "  " << pass.GetNbinsX() << "  " << total.GetNbinsX() << endl;
    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 == EFFRIG_SERGIO) {
    for (unsigned int i = 0; i < _selVector.size(); i++) {
      efficiency_(&(sel[i]), &(det[i]), &(eff[i]), &(errLow[i]), &(errHigh[i]));
    }

  }

  // Write the output files
  if (_outFileBase != "") {
    ofstream outTextFile((_outFileBase + "-eff-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];
      if (_selVector[i] != 0)
        outTextFile << setw(10) << eff[i] << setw(10) << errLow[i] << setw(10) << errHigh[i] << "\n";
      else
        outTextFile << setw(10) << 0. << setw(10) << 0. << setw(10) << 0. << endl;

    }

    outTextFile.close();
    if (_errMethod == EFFRIG_ROOT) {
      TFile outRootFile((_outFileBase + "-eff-rig.root"), "RECREATE");
      outRootFile.cd();
      errGraph.Write();
      outRootFile.Close();
    }
  }

}
1	pam-fi	1.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	pam-fi	1.3	#include "TGraphAsymmErrors.h"
13
14	pam-fi	1.2	extern "C" {
15			bool efficiency_(Int_t, Int_t, Double_t, Double_t, Double_t*);
16			}
17
18	pam-fi	1.3	EffRigCollection::EffRigCollection(const char *collectionName, TString outFileBase, TString rigBinsFile, int errMethod) :
19			EffCollection(collectionName, outFileBase, errMethod), _bins(0), _selVector(0), _detVector(0), _outUp(0), _outDown(0) {
20	pam-fi	1.1
21			ifstream rigBinListFile;
22			rigBinListFile.open(rigBinsFile);
23
24			TString auxString;
25			while (!rigBinListFile.eof()) {
26			rigBinListFile >> auxString;
27			if (auxString != "") {
28			_bins.push_back(auxString.Atof());
29			}
30			}
31			rigBinListFile.close();
32
33			_selVector.resize(_bins.size() - 1, 0);
34			_detVector.resize(_bins.size() - 1, 0);
35			}
36
37			int EffRigCollection::ApplyCut(PamLevel2 *event) {
38
39			_nEv++;
40			if (_selCollection.ApplyCut(event) == CUTOK) {
41			// Check if the event is inside the rigidity range
42			// NOTE: at this point a TrkPhSinCut should be already performed,
43			// since we are going to retrieve rigidity.
44	pam-fi	1.3
45			float rig = event->GetTrack(0)->GetTrkTrack()->GetRigidity();
46
47	pam-fi	1.1	if (rig >= _bins[0]) {
48			int i = 1;
49			while (rig >= _bins[i] && i < (int) _bins.size()) {
50			i++;
51			}
52			i--;
53
54			if (i < (int) (_selVector.size())) {
55			_selVector[i]++;
56	pam-fi	1.2	_sel++;
57	pam-fi	1.1	if (_detCollection.ApplyCut(event) == CUTOK) {
58			_detVector[i]++;
59	pam-fi	1.2	_det++;
60	pam-fi	1.1	_nGood++;
61			return CUTOK;
62			}
63			}
64			else
65			_outUp++;
66
67			}
68			else {
69			_outDown++;
70			return 0;
71			}
72			}
73
74			return 0;
75			}
76
77			void EffRigCollection::Finalize() {
78
79			// Print the report
80			EffCollection::Finalize();
81			cout << " Events below the minimum rigidity: " << _outDown << "\n";
82			cout << " Events above the maximum rigidity: " << _outUp << "\n";
83
84	pam-fi	1.2	// Compute the error
85			Int_t sel[_selVector.size()];
86			Int_t det[_detVector.size()];
87			Double_t eff[_selVector.size()];
88			Double_t errLow[_selVector.size()];
89			Double_t errHigh[_selVector.size()];
90			for (unsigned int i = 0; i < _selVector.size(); i++) {
91			sel[i] = (Int_t) _selVector[i];
92			det[i] = (Int_t) _detVector[i];
93			}
94	pam-fi	1.3
95			TGraphAsymmErrors errGraph;
96			errGraph.SetName(GetName());
97			errGraph.SetMarkerColor(kRed);
98			errGraph.SetMarkerStyle(7);
99			errGraph.GetYaxis()->SetRangeUser(0, 1.2);
100
101			if (_errMethod == EFFRIG_ROOT) {
102			double binning[_bins.size()];
103			for (unsigned int i = 0; i < _bins.size(); i++)
104			binning[i] = _bins[i];
105
106			TH1I pass("pass", "pass", _bins.size() - 1, binning);
107			TH1I total("total", "total", _bins.size() - 1, binning);
108			for (unsigned int i = 0; i < _selVector.size(); i++) {
109			cout << (binning[i + 1] + binning[i]) / 2. << " " << det[i] << " " << sel[i] << endl;
110			for (int j = 0; j < det[i]; j++)
111			pass.Fill((binning[i + 1] + binning[i]) / 2.);
112			for (int j = 0; j < sel[i]; j++)
113			total.Fill((binning[i + 1] + binning[i]) / 2.);
114			}
115
116			errGraph.BayesDivide(&pass, &total);
117			Double_t graphX;
118			double currBin;
119			int currPoint = 0;
120			//cout << errGraph.GetN() << " " << pass.GetNbinsX() << " " << total.GetNbinsX() << endl;
121			for (unsigned int i = 0; i < _selVector.size(); i++) {
122			errGraph.GetPoint(currPoint, graphX, eff[i]);
123			currBin = (binning[i + 1] + binning[i]) / 2.;
124			if (currBin == graphX) {
125			if (_selVector[i] < 8) {
126			eff[i] = 1.1;
127			errLow[i] = errHigh[i] = 0.;
128			errGraph.SetPoint(currPoint, graphX, 1.1);
129			float halfBin = (binning[i + 1] - binning[i]) / 2.;
130			errGraph.SetPointError(currPoint, halfBin, halfBin, 0., 0.);
131			}
132			currPoint++;
133			}
134			}
135
136			}
137			if (_errMethod == EFFRIG_SERGIO) {
138			for (unsigned int i = 0; i < _selVector.size(); i++) {
139			efficiency_(&(sel[i]), &(det[i]), &(eff[i]), &(errLow[i]), &(errHigh[i]));
140			}
141
142			}
143
144	pam-fi	1.1	// Write the output files
145			if (_outFileBase != "") {
146			ofstream outTextFile((_outFileBase + "-eff-rig.txt").Data(), ios_base::out);
147	pam-fi	1.2	streamsize newPrec = 4;
148			outTextFile.precision(newPrec);
149			outTextFile.setf(ios::fixed, ios::floatfield);
150	pam-fi	1.1	for (unsigned int i = 0; i < _selVector.size(); i++) {
151	pam-fi	1.2	outTextFile << setw(10) << _detVector[i] << setw(10) << _selVector[i];
152	pam-fi	1.1	if (_selVector[i] != 0)
153	pam-fi	1.2	outTextFile << setw(10) << eff[i] << setw(10) << errLow[i] << setw(10) << errHigh[i] << "\n";
154	pam-fi	1.1	else
155	pam-fi	1.2	outTextFile << setw(10) << 0. << setw(10) << 0. << setw(10) << 0. << endl;
156	pam-fi	1.1
157			}
158	pam-fi	1.3
159	pam-fi	1.1	outTextFile.close();
160	pam-fi	1.3	if (_errMethod == EFFRIG_ROOT) {
161			TFile outRootFile((_outFileBase + "-eff-rig.root"), "RECREATE");
162			outRootFile.cd();
163			errGraph.Write();
164			outRootFile.Close();
165			}
166	pam-fi	1.1	}
167
168			}