Selection/src/applySelection.cc

//
// System headers
//
#include <vector>                   // STL vector class
#include <iostream>                 // standard I/O
#include <iomanip>                  // functions to format standard I/O
#include <fstream>                  // functions for file I/O
#include <string>                   // C++ string class
#include <sstream>                  // class for parsing strings
#include <assert.h>
#include <stdlib.h>   
#include <getopt.h>
using namespace std;

//
// ROOT headers
//
#include <TROOT.h>                  // access to gROOT, entry point to ROOT system
#include <TSystem.h>                // interface to OS
#include <TFile.h>                  // file handle class
#include <TNtuple.h> 
#include <TTree.h>                  // class to access ntuples
#include <TChain.h>                 // 
#include <TBranch.h>                // class to access branches in TTree
#include <TClonesArray.h>           // ROOT array class
#include <TCanvas.h>                // class for drawing
#include <TH1F.h>                   // 1D histograms
#include <TBenchmark.h>             // class to track macro running statistics
#include <TLorentzVector.h>         // 4-vector class
#include <TVector3.h>               // 3D vector class


//
// ntuple format headers
//
#include "HiggsAnaDefs.hh"
#include "TEventInfo.hh"
#include "TElectron.hh"
#include "TMuon.hh"
#include "TJet.hh"
#include "SiMVAElectronSelection.h"
#include "RunLumiRangeMap.h"
#include "EfficiencyWeightsInterface.h"

//
// utility headers
//
#include "ParseArgs.h"
#include "SampleWeight.h"
#include "Selection.h"
#include "HZZCiCElectronSelection.h"
#include "HZZLikelihoodElectronSelection.h"
#include "HZZBDTElectronSelection.h"
#include "SampleWeight.h"

//#define BDTFAIL_THEIR_EVENTS
//#define THEIR_EVENTS

Double_t etaFromY(Double_t pt, Double_t y, Double_t phi, Double_t m)
{
  Double_t a  = (1+exp(2*y))/(exp(2*y)-1); // intermediate term
  if(a*a<1) { cout << "a too small" << endl; assert(0); }
  Double_t E  = sqrt( a*a*(pt*pt+m*m)/(a*a-1) );
  Double_t pz = E*E - pt*pt - m*m;
  if(pz<0) { cout << "imag. pz" << endl; assert(0); }
  pz = sqrt(pz);
  if(y<0) pz *= -1;
  TLorentzVector v;
  v.SetPxPyPzE(pt*cos(phi),pt*sin(phi),pz,E);
  Double_t th = v.Theta();
  return -log(tan(th/2));
};


int getChannel(mithep::TGenInfo * ginfo) {

    int gchannel=-1;

    if( abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kElectron )  gchannel=0;
    else if( abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kMuon ) gchannel=1;
    else if( (abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kMuon) || 
             (abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kElectron) ) gchannel=2; 
    
    return gchannel;
}


//=== MAIN =================================================================================================
int main(int argc, char** argv) 
{
  
  vector<vector<string> > inputFiles;
  inputFiles.push_back(vector<string>());

  //
  // args
  //--------------------------------------------------------------------------------------------------------------
  ControlFlags ctrl;
  parse_args( argc, argv, ctrl );
  if( ctrl.inputfile.empty() || 
      ctrl.inputfile.empty() ) {
    cerr << "usage: applySelection.exe <flags> " << endl;
    cerr << "\tmandoatory flags : " << endl;
    cerr << "\t--inputfile |  file containing a list of ntuples to run over" << endl;
    cerr << "\t--outputfile | file to store  selected evet" << endl;
    return 1;
  }
  ctrl.dump();

  map<string,double> entrymap; // number of unskimmed entries in each file

  //
  // File I/O
  //--------------------------------------------------------------------------------------------------------------
  TChain * chain = new TChain("Events");
  ifstream f(ctrl.inputfile.c_str());
  string fname;
  while (f >> fname) { 
    if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
    cout << "adding inputfile : " << fname.c_str() << endl;
    entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
    chain->AddFile(fname.c_str());
  }

  // table of cross section values   
  SimpleTable xstab("./data/xs.dat");

  //
  // Setup
  //--------------------------------------------------------------------------------------------------------------
  const char * ofname;
  if( strcmp( ctrl.outputfile.c_str(), "") ) {
    ofname = ctrl.outputfile.c_str();
  } else {
    ofname = "tmp.root";
  }
  TFile *file = new TFile(ofname, "RECREATE");
  TH1F * h_evt = new TH1F("hevt", "hevt", 2, 0, 2 ); 
  TH1F * h_evtfail = new TH1F("hevtfail", "hevtfail", 1024, 0, 1024 );
  TTree * passtuple = new TTree("passtuple", "passtuple" );
  unsigned run, evt, lumi, channel, gchannel;
  float mZ1, mZ2, m4l, pt4l;
  unsigned tZ1, tZ2;
  double w, woff, werroff, won, werron;
  if( ctrl.mc )   { 
    passtuple->Branch("gchannel", &gchannel);
    passtuple->Branch("woff", &woff);
    passtuple->Branch("werroff", &werroff);
    passtuple->Branch("won", &won);
    passtuple->Branch("werron", &werron);
  }
  passtuple->Branch("channel", &channel);
  passtuple->Branch("run", &run);
  passtuple->Branch("evt", &evt);
  passtuple->Branch("lumi", &lumi);
  passtuple->Branch("mZ1", &mZ1);
  passtuple->Branch("mZ2", &mZ2);
  passtuple->Branch("tZ1", &tZ1);
  passtuple->Branch("tZ2", &tZ2);
  passtuple->Branch("m4l", &m4l);
  passtuple->Branch("pt4l", &pt4l);
  passtuple->Branch("w", &w);

  initCiCSelection();
  if(ctrl.eleSele=="lik") initLikSelection();
  if(ctrl.eleSele=="bdt") initBDTSelection();
  // if(ctrl.eleSele=="si" ) initSiMVAElectronSelection(); 
  initRunLumiRangeMap();
  initEfficiencyWeights();


  //
  // Loop
  //--------------------------------------------------------------------------------------------------------------
  if(ctrl.fakeScheme!="none") {cout << "wrong executable!" << endl; assert(0);}
  
  //
  // Access samples and fill histograms
  TFile *inputFile=0;
  TTree *eventTree=0;  
   
  // Data structures to store info from TTrees
  mithep::TEventInfo *info    = new mithep::TEventInfo();
  mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
  TClonesArray *electronArr   = new TClonesArray("mithep::TElectron");
  TClonesArray *muonArr       = new TClonesArray("mithep::TMuon");
  

  chain->SetBranchAddress("Info",       &info);
  chain->SetBranchAddress("Electron", &electronArr);
  chain->SetBranchAddress("Muon", &muonArr);
  if( ctrl.mc ) { chain->SetBranchAddress("Gen",        &ginfo);}


  int count=0, pass=0;
  float denom[3]={0.,0.,0.};
  float numer[3]={0.,0.,0.};
  float numercorr[3]={0.,0.,0.};

  UInt_t imax = chain->GetEntries();
  cout << "nEntries: " << imax << endl;


  for(UInt_t ientry=0; ientry<imax; ientry++) {
    chain->GetEntry(ientry);
    if(!(ientry%100000)) cerr << "entry: " << ientry << endl;

    // get event weight for this file                                                     
    double eventweight=1;
    if(ctrl.mc) eventweight = getWeight(xstab,entrymap,chain);

#ifdef THEIR_EVENTS
    if( !( info->evtNum == 504867308  ||
           info->evtNum == 368148849  ||
           //      info->evtNum == 129514273  ||
           info->evtNum == 286336207  ||
           info->evtNum == 344708580  ||
           info->evtNum == 30998576   ||
           info->evtNum == 155679852  ||
           info->evtNum == 394010457  ||
           info->evtNum == 917379387  ||
           info->evtNum == 78213037   ||
           info->evtNum == 862270386  ||
           info->evtNum == 337493970  ||  // not baseline anymore?
           info->evtNum == 1491724484 ||
           info->evtNum == 480301165  ||
           info->evtNum == 1038911933 ||
           info->evtNum == 876658967  ||
           info->evtNum == 966824024  ||
           info->evtNum == 141954801  ||
           info->evtNum == 160966858  ||
           info->evtNum == 191231387  ||
           info->evtNum == 66033190   ||
           info->evtNum == 10347106   ||
           info->evtNum == 107360878  ||
           info->evtNum == 2554393033 ||
           info->evtNum == 933807102  ||
           info->evtNum == 1188043146 ||
           info->evtNum == 559839432  ||
           info->evtNum == 16706390   ||
           info->evtNum == 65557571   ||
           info->evtNum == 389185367  ||
           info->evtNum == 2722114329 ) ) continue;
#endif
 

#ifdef BDTFAIL_THEIR_EVENTS
    if( !( info->evtNum == 78213037   ||
           info->evtNum == 876658967  ) ) continue;
           
#endif

    unsigned evtfail = fails_HZZ4L_selection(ctrl, info, electronArr, muonArr, eventweight, passtuple, NULL);
    //    if( !evtfail ) { 
//       cout << "run: " << info->runNum
//         << "\tevt: " << info->evtNum
//         << "\tlumi: " << info->lumiSec
//         << "\tfile: " << chain->GetFile()->GetEndpointUrl()->GetFile()
//         << endl;
//    }

#ifdef THEIR_EVENTS
    if(evtfail!=0) cout << "Failed!" << endl;
#endif
    
    h_evtfail->Fill( evtfail, eventweight );

    if( ctrl.mc ) {
      // BR is implicit here ...
      int gchannel = getChannel(ginfo); 
        if( gchannel >= 0 ) { 
//        denom[gchannel] += eventweight;
//        if( !evtfail ) numer[gchannel] += eventweight;
          denom[gchannel] += 1;
          if( !evtfail ) numer[gchannel] += 1;

          // get eff weight for current evt
          if( !evtfail ) { 
            double thewoff, thewon;
            passtuple->SetBranchAddress("woff", &thewoff); // must put this here ...
            passtuple->SetBranchAddress("won",  &thewon); // must put this here ...
            //      cout << "thewoff: " << thewoff << "\tthewon: " << thewon << endl; 
            int ret = passtuple->GetEntry(passtuple->GetEntries()-1); 
            numercorr[gchannel] += thewoff*thewon;        
          }
        }
    }

    count++;
    if( !evtfail ) pass++;
          
  } //end loop over data     

  if( ctrl.mc ) { 
    cout << "--------------" << endl;
    cout << "denom: " << count << endl; 
    cout << "pass: " << pass << endl; 
    cout << "axe: " << (float)pass/count << endl;
    cout << "--------------" << endl;
  }

  delete info;
  delete electronArr;
  delete muonArr;


  //--------------------------------------------------------------------------------------------------------------
  // Save Histograms;
  //============================================================================================================== 
  // const char * ofname;
  // if( strcmp( ctrl.outputfile.c_str(), "") ) {
  //   ofname = ctrl.outputfile.c_str();
  // } else {
  //   ofname = "tmp.root";
  // }

  // TFile *file = new TFile(ofname, "RECREATE");
  gROOT->cd(file->GetPath());
  h_evt->Write();
  h_evtfail->Write();
  passtuple->Write();
  file->Close();
  delete file;

  // map<TString,TMVA::Reader*>::iterator it;
  // for(it=readers.begin(); it!=readers.end(); it++) delete (*it).second;


  float acc_4e        = numer[0]/denom[0];
  float acc_4mu       = numer[1]/denom[1];
  float acc_2e2mu     = numer[2]/denom[2];
  float acc_err_4e    = (1./denom[0])*sqrt(numer[0]*(1-acc_4e));
  float acc_err_4mu   = (1./denom[1])*sqrt(numer[1]*(1-acc_4mu));
  float acc_err_2e2mu = (1./denom[2])*sqrt(numer[2]*(1-acc_2e2mu));

  float acc_c_4e        = numercorr[0]/denom[0];
  float acc_c_4mu       = numercorr[1]/denom[1];
  float acc_c_2e2mu     = numercorr[2]/denom[2];
  float acc_c_err_4e    = (1./denom[0])*sqrt(numer[0]*(1-acc_4e));
  float acc_c_err_4mu   = (1./denom[1])*sqrt(numer[1]*(1-acc_4mu));
  float acc_c_err_2e2mu = (1./denom[2])*sqrt(numer[2]*(1-acc_2e2mu));

  if( ctrl.mc ) {
  cerr << "---------------------" << endl;
//   cerr << "acc_4e: " << fixed << setw(8) << setprecision(7) << numer[0] << endl; // numer already weighted ...
//   cerr << "acc_4mu: " << fixed << setw(8) << setprecision(7) << numer[1] <<  endl;
//   cerr << "acc_2e2mu: " << fixed << setw(8) << setprecision(7) << numer[2] << endl;
  cerr << "acc_4e: " << fixed << setw(8) << setprecision(7) << acc_4e << " +/- " << acc_err_4e << endl; 
  cerr << "acc_4mu: " << fixed << setw(8) << setprecision(7) << acc_4mu << " +/- " << acc_err_4mu << endl; 
  cerr << "acc_2e2mu: " << fixed << setw(8) << setprecision(7) << acc_2e2mu << " +/- " << acc_err_2e2mu << endl; 
  cerr << "-------corrected------" << endl;
  cerr << "acc_c_4e: " << fixed << setw(8) << setprecision(7) << acc_c_4e << " +/- " << acc_c_err_4e << endl; 
  cerr << "acc_c_4mu: " << fixed << setw(8) << setprecision(7) << acc_c_4mu << " +/- " << acc_c_err_4mu << endl; 
  cerr << "acc_c_2e2mu: " << fixed << setw(8) << setprecision(7) << acc_c_2e2mu << " +/- " << acc_c_err_2e2mu << endl; 
  cerr << "---------------------" << endl;
}

  cerr << "done!" << endl;
} 


Revision:	1.12
Committed:	Sun Oct 30 15:35:25 2011 UTC (13 years, 6 months ago) by khahn
Content type:	text/plain
Branch:	MAIN
CVS Tags:	synched, AN490
Changes since 1.11:	+1 -1 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	khahn	1.1	//
2			// System headers
3			//
4			#include <vector> // STL vector class
5			#include <iostream> // standard I/O
6			#include <iomanip> // functions to format standard I/O
7			#include <fstream> // functions for file I/O
8			#include <string> // C++ string class
9			#include <sstream> // class for parsing strings
10			#include <assert.h>
11			#include <stdlib.h>
12			#include <getopt.h>
13			using namespace std;
14
15			//
16			// ROOT headers
17			//
18			#include <TROOT.h> // access to gROOT, entry point to ROOT system
19			#include <TSystem.h> // interface to OS
20			#include <TFile.h> // file handle class
21			#include <TNtuple.h>
22			#include <TTree.h> // class to access ntuples
23			#include <TChain.h> //
24			#include <TBranch.h> // class to access branches in TTree
25			#include <TClonesArray.h> // ROOT array class
26			#include <TCanvas.h> // class for drawing
27			#include <TH1F.h> // 1D histograms
28			#include <TBenchmark.h> // class to track macro running statistics
29			#include <TLorentzVector.h> // 4-vector class
30			#include <TVector3.h> // 3D vector class
31
32	khahn	1.11
33	khahn	1.1	//
34			// ntuple format headers
35			//
36	dkralph	1.2	#include "HiggsAnaDefs.hh"
37	khahn	1.1	#include "TEventInfo.hh"
38			#include "TElectron.hh"
39			#include "TMuon.hh"
40			#include "TJet.hh"
41	khahn	1.8	#include "SiMVAElectronSelection.h"
42	khahn	1.1	#include "RunLumiRangeMap.h"
43	khahn	1.11	#include "EfficiencyWeightsInterface.h"
44	khahn	1.1
45			//
46			// utility headers
47			//
48			#include "ParseArgs.h"
49	dkralph	1.2	#include "SampleWeight.h"
50	khahn	1.1	#include "Selection.h"
51			#include "HZZCiCElectronSelection.h"
52	dkralph	1.4	#include "HZZLikelihoodElectronSelection.h"
53	dkralph	1.5	#include "HZZBDTElectronSelection.h"
54	dkralph	1.2	#include "SampleWeight.h"
55	khahn	1.1
56	khahn	1.8	//#define BDTFAIL_THEIR_EVENTS
57	khahn	1.1	//#define THEIR_EVENTS
58
59	khahn	1.11	Double_t etaFromY(Double_t pt, Double_t y, Double_t phi, Double_t m)
60			{
61			Double_t a = (1+exp(2y))/(exp(2y)-1); // intermediate term
62			if(a*a<1) { cout << "a too small" << endl; assert(0); }
63			Double_t E = sqrt( aa(ptpt+mm)/(a*a-1) );
64			Double_t pz = EE - ptpt - m*m;
65			if(pz<0) { cout << "imag. pz" << endl; assert(0); }
66			pz = sqrt(pz);
67			if(y<0) pz *= -1;
68			TLorentzVector v;
69			v.SetPxPyPzE(ptcos(phi),ptsin(phi),pz,E);
70			Double_t th = v.Theta();
71			return -log(tan(th/2));
72			};
73
74
75			int getChannel(mithep::TGenInfo * ginfo) {
76
77			int gchannel=-1;
78
79			if( abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kElectron ) gchannel=0;
80			else if( abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kMuon ) gchannel=1;
81			else if( (abs(ginfo->id_1_a) == EGenType::kElectron && abs(ginfo->id_1_b) == EGenType::kMuon) \|\|
82			(abs(ginfo->id_1_a) == EGenType::kMuon && abs(ginfo->id_1_b) == EGenType::kElectron) ) gchannel=2;
83
84			return gchannel;
85			}
86
87
88
89	khahn	1.1	//=== MAIN =================================================================================================
90			int main(int argc, char** argv)
91			{
92
93			vector<vector<string> > inputFiles;
94			inputFiles.push_back(vector<string>());
95
96			//
97			// args
98			//--------------------------------------------------------------------------------------------------------------
99			ControlFlags ctrl;
100			parse_args( argc, argv, ctrl );
101			if( ctrl.inputfile.empty() \|\|
102			ctrl.inputfile.empty() ) {
103			cerr << "usage: applySelection.exe <flags> " << endl;
104			cerr << "\tmandoatory flags : " << endl;
105			cerr << "\t--inputfile \| file containing a list of ntuples to run over" << endl;
106			cerr << "\t--outputfile \| file to store selected evet" << endl;
107			return 1;
108			}
109			ctrl.dump();
110
111	dkralph	1.4	map<string,double> entrymap; // number of unskimmed entries in each file
112
113	khahn	1.1	//
114			// File I/O
115			//--------------------------------------------------------------------------------------------------------------
116			TChain * chain = new TChain("Events");
117			ifstream f(ctrl.inputfile.c_str());
118			string fname;
119			while (f >> fname) {
120	dkralph	1.4	if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
121	khahn	1.1	cout << "adding inputfile : " << fname.c_str() << endl;
122	dkralph	1.4	entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
123	khahn	1.1	chain->AddFile(fname.c_str());
124			}
125
126	dkralph	1.2	// table of cross section values
127			SimpleTable xstab("./data/xs.dat");
128
129	khahn	1.1	//
130			// Setup
131			//--------------------------------------------------------------------------------------------------------------
132	dkralph	1.10	const char * ofname;
133			if( strcmp( ctrl.outputfile.c_str(), "") ) {
134			ofname = ctrl.outputfile.c_str();
135			} else {
136			ofname = "tmp.root";
137			}
138			TFile *file = new TFile(ofname, "RECREATE");
139	khahn	1.1	TH1F * h_evt = new TH1F("hevt", "hevt", 2, 0, 2 );
140	dkralph	1.4	TH1F * h_evtfail = new TH1F("hevtfail", "hevtfail", 1024, 0, 1024 );
141	khahn	1.6	TTree * passtuple = new TTree("passtuple", "passtuple" );
142	khahn	1.11	unsigned run, evt, lumi, channel, gchannel;
143	khahn	1.7	float mZ1, mZ2, m4l, pt4l;
144	khahn	1.9	unsigned tZ1, tZ2;
145	khahn	1.11	double w, woff, werroff, won, werron;
146			if( ctrl.mc ) {
147			passtuple->Branch("gchannel", &gchannel);
148			passtuple->Branch("woff", &woff);
149			passtuple->Branch("werroff", &werroff);
150			passtuple->Branch("won", &won);
151			passtuple->Branch("werron", &werron);
152			}
153			passtuple->Branch("channel", &channel);
154	khahn	1.6	passtuple->Branch("run", &run);
155			passtuple->Branch("evt", &evt);
156			passtuple->Branch("lumi", &lumi);
157			passtuple->Branch("mZ1", &mZ1);
158			passtuple->Branch("mZ2", &mZ2);
159	khahn	1.9	passtuple->Branch("tZ1", &tZ1);
160			passtuple->Branch("tZ2", &tZ2);
161	khahn	1.6	passtuple->Branch("m4l", &m4l);
162			passtuple->Branch("pt4l", &pt4l);
163			passtuple->Branch("w", &w);
164	khahn	1.11
165	khahn	1.1	initCiCSelection();
166	dkralph	1.5	if(ctrl.eleSele=="lik") initLikSelection();
167			if(ctrl.eleSele=="bdt") initBDTSelection();
168	dkralph	1.10	// if(ctrl.eleSele=="si" ) initSiMVAElectronSelection();
169	khahn	1.1	initRunLumiRangeMap();
170	khahn	1.11	initEfficiencyWeights();
171
172
173
174	khahn	1.1
175			//
176			// Loop
177			//--------------------------------------------------------------------------------------------------------------
178	dkralph	1.10	if(ctrl.fakeScheme!="none") {cout << "wrong executable!" << endl; assert(0);}
179	khahn	1.1
180			//
181			// Access samples and fill histograms
182			TFile *inputFile=0;
183			TTree *eventTree=0;
184
185			// Data structures to store info from TTrees
186			mithep::TEventInfo *info = new mithep::TEventInfo();
187	khahn	1.11	mithep::TGenInfo *ginfo = new mithep::TGenInfo();
188	khahn	1.1	TClonesArray *electronArr = new TClonesArray("mithep::TElectron");
189			TClonesArray *muonArr = new TClonesArray("mithep::TMuon");
190
191
192			chain->SetBranchAddress("Info", &info);
193			chain->SetBranchAddress("Electron", &electronArr);
194			chain->SetBranchAddress("Muon", &muonArr);
195	khahn	1.11	if( ctrl.mc ) { chain->SetBranchAddress("Gen", &ginfo);}
196	khahn	1.1
197
198	khahn	1.9	int count=0, pass=0;
199	khahn	1.11	float denom[3]={0.,0.,0.};
200			float numer[3]={0.,0.,0.};
201			float numercorr[3]={0.,0.,0.};
202
203	dkralph	1.10	UInt_t imax = chain->GetEntries();
204	khahn	1.11	cout << "nEntries: " << imax << endl;
205
206
207	dkralph	1.10	for(UInt_t ientry=0; ientry<imax; ientry++) {
208	khahn	1.1	chain->GetEntry(ientry);
209	dkralph	1.10	if(!(ientry%100000)) cerr << "entry: " << ientry << endl;
210	khahn	1.1
211	dkralph	1.2	// get event weight for this file
212	dkralph	1.4	double eventweight=1;
213			if(ctrl.mc) eventweight = getWeight(xstab,entrymap,chain);
214	dkralph	1.2
215	khahn	1.1	#ifdef THEIR_EVENTS
216			if( !( info->evtNum == 504867308 \|\|
217			info->evtNum == 368148849 \|\|
218	khahn	1.9	// info->evtNum == 129514273 \|\|
219	khahn	1.1	info->evtNum == 286336207 \|\|
220			info->evtNum == 344708580 \|\|
221			info->evtNum == 30998576 \|\|
222			info->evtNum == 155679852 \|\|
223			info->evtNum == 394010457 \|\|
224			info->evtNum == 917379387 \|\|
225			info->evtNum == 78213037 \|\|
226	khahn	1.9	info->evtNum == 862270386 \|\|
227			info->evtNum == 337493970 \|\| // not baseline anymore?
228	khahn	1.1	info->evtNum == 1491724484 \|\|
229			info->evtNum == 480301165 \|\|
230			info->evtNum == 1038911933 \|\|
231			info->evtNum == 876658967 \|\|
232			info->evtNum == 966824024 \|\|
233			info->evtNum == 141954801 \|\|
234			info->evtNum == 160966858 \|\|
235			info->evtNum == 191231387 \|\|
236			info->evtNum == 66033190 \|\|
237			info->evtNum == 10347106 \|\|
238	khahn	1.9	info->evtNum == 107360878 \|\|
239			info->evtNum == 2554393033 \|\|
240			info->evtNum == 933807102 \|\|
241			info->evtNum == 1188043146 \|\|
242			info->evtNum == 559839432 \|\|
243			info->evtNum == 16706390 \|\|
244			info->evtNum == 65557571 \|\|
245			info->evtNum == 389185367 \|\|
246			info->evtNum == 2722114329 ) ) continue;
247	khahn	1.1	#endif
248	khahn	1.8
249
250			#ifdef BDTFAIL_THEIR_EVENTS
251			if( !( info->evtNum == 78213037 \|\|
252			info->evtNum == 876658967 ) ) continue;
253
254			#endif
255	khahn	1.1
256	khahn	1.12	unsigned evtfail = fails_HZZ4L_selection(ctrl, info, electronArr, muonArr, eventweight, passtuple, NULL);
257	khahn	1.11	// if( !evtfail ) {
258			// cout << "run: " << info->runNum
259			// << "\tevt: " << info->evtNum
260			// << "\tlumi: " << info->lumiSec
261			// << "\tfile: " << chain->GetFile()->GetEndpointUrl()->GetFile()
262			// << endl;
263			// }
264	dkralph	1.10
265			#ifdef THEIR_EVENTS
266			if(evtfail!=0) cout << "Failed!" << endl;
267			#endif
268
269	khahn	1.1	h_evtfail->Fill( evtfail, eventweight );
270	khahn	1.8
271	khahn	1.11	if( ctrl.mc ) {
272			// BR is implicit here ...
273			int gchannel = getChannel(ginfo);
274			if( gchannel >= 0 ) {
275			// denom[gchannel] += eventweight;
276			// if( !evtfail ) numer[gchannel] += eventweight;
277			denom[gchannel] += 1;
278			if( !evtfail ) numer[gchannel] += 1;
279
280			// get eff weight for current evt
281			if( !evtfail ) {
282			double thewoff, thewon;
283			passtuple->SetBranchAddress("woff", &thewoff); // must put this here ...
284			passtuple->SetBranchAddress("won", &thewon); // must put this here ...
285			// cout << "thewoff: " << thewoff << "\tthewon: " << thewon << endl;
286			int ret = passtuple->GetEntry(passtuple->GetEntries()-1);
287			numercorr[gchannel] += thewoff*thewon;
288			}
289			}
290			}
291
292	khahn	1.9	count++;
293			if( !evtfail ) pass++;
294	khahn	1.1
295			} //end loop over data
296
297	khahn	1.9	if( ctrl.mc ) {
298			cout << "--------------" << endl;
299			cout << "denom: " << count << endl;
300			cout << "pass: " << pass << endl;
301			cout << "axe: " << (float)pass/count << endl;
302			cout << "--------------" << endl;
303			}
304	khahn	1.1
305			delete info;
306			delete electronArr;
307			delete muonArr;
308
309
310
311			//--------------------------------------------------------------------------------------------------------------
312			// Save Histograms;
313			//==============================================================================================================
314	dkralph	1.10	// const char * ofname;
315			// if( strcmp( ctrl.outputfile.c_str(), "") ) {
316			// ofname = ctrl.outputfile.c_str();
317			// } else {
318			// ofname = "tmp.root";
319			// }
320	khahn	1.1
321	dkralph	1.10	// TFile *file = new TFile(ofname, "RECREATE");
322			gROOT->cd(file->GetPath());
323	khahn	1.1	h_evt->Write();
324			h_evtfail->Write();
325			passtuple->Write();
326			file->Close();
327			delete file;
328
329	dkralph	1.5	// map<TString,TMVA::Reader*>::iterator it;
330			// for(it=readers.begin(); it!=readers.end(); it++) delete (*it).second;
331
332	khahn	1.11
333			float acc_4e = numer[0]/denom[0];
334			float acc_4mu = numer[1]/denom[1];
335			float acc_2e2mu = numer[2]/denom[2];
336			float acc_err_4e = (1./denom[0])sqrt(numer[0](1-acc_4e));
337			float acc_err_4mu = (1./denom[1])sqrt(numer[1](1-acc_4mu));
338			float acc_err_2e2mu = (1./denom[2])sqrt(numer[2](1-acc_2e2mu));
339
340			float acc_c_4e = numercorr[0]/denom[0];
341			float acc_c_4mu = numercorr[1]/denom[1];
342			float acc_c_2e2mu = numercorr[2]/denom[2];
343			float acc_c_err_4e = (1./denom[0])sqrt(numer[0](1-acc_4e));
344			float acc_c_err_4mu = (1./denom[1])sqrt(numer[1](1-acc_4mu));
345			float acc_c_err_2e2mu = (1./denom[2])sqrt(numer[2](1-acc_2e2mu));
346
347			if( ctrl.mc ) {
348			cerr << "---------------------" << endl;
349			// cerr << "acc_4e: " << fixed << setw(8) << setprecision(7) << numer[0] << endl; // numer already weighted ...
350			// cerr << "acc_4mu: " << fixed << setw(8) << setprecision(7) << numer[1] << endl;
351			// cerr << "acc_2e2mu: " << fixed << setw(8) << setprecision(7) << numer[2] << endl;
352			cerr << "acc_4e: " << fixed << setw(8) << setprecision(7) << acc_4e << " +/- " << acc_err_4e << endl;
353			cerr << "acc_4mu: " << fixed << setw(8) << setprecision(7) << acc_4mu << " +/- " << acc_err_4mu << endl;
354			cerr << "acc_2e2mu: " << fixed << setw(8) << setprecision(7) << acc_2e2mu << " +/- " << acc_err_2e2mu << endl;
355			cerr << "-------corrected------" << endl;
356			cerr << "acc_c_4e: " << fixed << setw(8) << setprecision(7) << acc_c_4e << " +/- " << acc_c_err_4e << endl;
357			cerr << "acc_c_4mu: " << fixed << setw(8) << setprecision(7) << acc_c_4mu << " +/- " << acc_c_err_4mu << endl;
358			cerr << "acc_c_2e2mu: " << fixed << setw(8) << setprecision(7) << acc_c_2e2mu << " +/- " << acc_c_err_2e2mu << endl;
359			cerr << "---------------------" << endl;
360			}
361
362	khahn	1.1	cerr << "done!" << endl;
363			}
364
365
366