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"

//
// 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

//=== 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 );
  //  TNtuple * passtuple = new TNtuple( "passtuple", "passtuple", "run:evt:lumi:channel:mZ1:mZ2:m4l:pt4l:w" );
  TTree * passtuple = new TTree("passtuple", "passtuple" );
  unsigned run, evt, lumi, channel;
  float mZ1, mZ2, m4l, pt4l;
  unsigned tZ1, tZ2;
  double w;
  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();

  //
  // 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();
  TClonesArray *electronArr   = new TClonesArray("mithep::TElectron");
  TClonesArray *muonArr       = new TClonesArray("mithep::TMuon");
  

  TBranch *infoBr;
  TBranch *electronBr;
  TBranch *muonBr;

//   chain->SetBranchStatus("*",0); //disable all branches
//   chain->SetBranchStatus("Info", 1);
//   chain->SetBranchStatus("Muon", 1);
//   chain->SetBranchStatus("Electron", 1);

  chain->SetBranchAddress("Info",       &info);
  chain->SetBranchAddress("Electron", &electronArr);
  chain->SetBranchAddress("Muon", &muonArr);


  int count=0, pass=0;
  cout << "nEntries: " << chain->GetEntries() << endl;
  cerr << "nEntries: " << chain->GetEntries() << endl;
  UInt_t imax = chain->GetEntries();
  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);

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

    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;

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


Revision:	1.10
Committed:	Sun Oct 23 11:52:41 2011 UTC (13 years, 6 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
Changes since 1.9:	+27 -12 lines
Log Message:	* empty log message *
#	Content
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	//
33	// ntuple format headers
34	//
35	#include "HiggsAnaDefs.hh"
36	#include "TEventInfo.hh"
37	#include "TElectron.hh"
38	#include "TMuon.hh"
39	#include "TJet.hh"
40	#include "SiMVAElectronSelection.h"
41	#include "RunLumiRangeMap.h"
42
43	//
44	// utility headers
45	//
46	#include "ParseArgs.h"
47	#include "SampleWeight.h"
48	#include "Selection.h"
49	#include "HZZCiCElectronSelection.h"
50	#include "HZZLikelihoodElectronSelection.h"
51	#include "HZZBDTElectronSelection.h"
52	#include "SampleWeight.h"
53
54	//#define BDTFAIL_THEIR_EVENTS
55	//#define THEIR_EVENTS
56
57	//=== MAIN =================================================================================================
58	int main(int argc, char** argv)
59	{
60
61	vector<vector<string> > inputFiles;
62	inputFiles.push_back(vector<string>());
63
64	//
65	// args
66	//--------------------------------------------------------------------------------------------------------------
67	ControlFlags ctrl;
68	parse_args( argc, argv, ctrl );
69	if( ctrl.inputfile.empty() \|\|
70	ctrl.inputfile.empty() ) {
71	cerr << "usage: applySelection.exe <flags> " << endl;
72	cerr << "\tmandoatory flags : " << endl;
73	cerr << "\t--inputfile \| file containing a list of ntuples to run over" << endl;
74	cerr << "\t--outputfile \| file to store selected evet" << endl;
75	return 1;
76	}
77	ctrl.dump();
78
79	map<string,double> entrymap; // number of unskimmed entries in each file
80
81	//
82	// File I/O
83	//--------------------------------------------------------------------------------------------------------------
84	TChain * chain = new TChain("Events");
85	ifstream f(ctrl.inputfile.c_str());
86	string fname;
87	while (f >> fname) {
88	if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
89	cout << "adding inputfile : " << fname.c_str() << endl;
90	entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
91	chain->AddFile(fname.c_str());
92	}
93
94	// table of cross section values
95	SimpleTable xstab("./data/xs.dat");
96
97	//
98	// Setup
99	//--------------------------------------------------------------------------------------------------------------
100	const char * ofname;
101	if( strcmp( ctrl.outputfile.c_str(), "") ) {
102	ofname = ctrl.outputfile.c_str();
103	} else {
104	ofname = "tmp.root";
105	}
106	TFile *file = new TFile(ofname, "RECREATE");
107	TH1F * h_evt = new TH1F("hevt", "hevt", 2, 0, 2 );
108	TH1F * h_evtfail = new TH1F("hevtfail", "hevtfail", 1024, 0, 1024 );
109	// TNtuple * passtuple = new TNtuple( "passtuple", "passtuple", "run:evt:lumi:channel:mZ1:mZ2:m4l:pt4l:w" );
110	TTree * passtuple = new TTree("passtuple", "passtuple" );
111	unsigned run, evt, lumi, channel;
112	float mZ1, mZ2, m4l, pt4l;
113	unsigned tZ1, tZ2;
114	double w;
115	passtuple->Branch("run", &run);
116	passtuple->Branch("evt", &evt);
117	passtuple->Branch("lumi", &lumi);
118	passtuple->Branch("mZ1", &mZ1);
119	passtuple->Branch("mZ2", &mZ2);
120	passtuple->Branch("tZ1", &tZ1);
121	passtuple->Branch("tZ2", &tZ2);
122	passtuple->Branch("m4l", &m4l);
123	passtuple->Branch("pt4l", &pt4l);
124	passtuple->Branch("w", &w);
125	initCiCSelection();
126	if(ctrl.eleSele=="lik") initLikSelection();
127	if(ctrl.eleSele=="bdt") initBDTSelection();
128	// if(ctrl.eleSele=="si" ) initSiMVAElectronSelection();
129	initRunLumiRangeMap();
130
131	//
132	// Loop
133	//--------------------------------------------------------------------------------------------------------------
134	if(ctrl.fakeScheme!="none") {cout << "wrong executable!" << endl; assert(0);}
135
136	//
137	// Access samples and fill histograms
138	TFile *inputFile=0;
139	TTree *eventTree=0;
140
141	// Data structures to store info from TTrees
142	mithep::TEventInfo *info = new mithep::TEventInfo();
143	TClonesArray *electronArr = new TClonesArray("mithep::TElectron");
144	TClonesArray *muonArr = new TClonesArray("mithep::TMuon");
145
146
147	TBranch *infoBr;
148	TBranch *electronBr;
149	TBranch *muonBr;
150
151	// chain->SetBranchStatus("*",0); //disable all branches
152	// chain->SetBranchStatus("Info", 1);
153	// chain->SetBranchStatus("Muon", 1);
154	// chain->SetBranchStatus("Electron", 1);
155
156	chain->SetBranchAddress("Info", &info);
157	chain->SetBranchAddress("Electron", &electronArr);
158	chain->SetBranchAddress("Muon", &muonArr);
159
160
161	int count=0, pass=0;
162	cout << "nEntries: " << chain->GetEntries() << endl;
163	cerr << "nEntries: " << chain->GetEntries() << endl;
164	UInt_t imax = chain->GetEntries();
165	for(UInt_t ientry=0; ientry<imax; ientry++) {
166	chain->GetEntry(ientry);
167	if(!(ientry%100000)) cerr << "entry: " << ientry << endl;
168
169	// get event weight for this file
170	double eventweight=1;
171	if(ctrl.mc) eventweight = getWeight(xstab,entrymap,chain);
172
173	#ifdef THEIR_EVENTS
174	if( !( info->evtNum == 504867308 \|\|
175	info->evtNum == 368148849 \|\|
176	// info->evtNum == 129514273 \|\|
177	info->evtNum == 286336207 \|\|
178	info->evtNum == 344708580 \|\|
179	info->evtNum == 30998576 \|\|
180	info->evtNum == 155679852 \|\|
181	info->evtNum == 394010457 \|\|
182	info->evtNum == 917379387 \|\|
183	info->evtNum == 78213037 \|\|
184	info->evtNum == 862270386 \|\|
185	info->evtNum == 337493970 \|\| // not baseline anymore?
186	info->evtNum == 1491724484 \|\|
187	info->evtNum == 480301165 \|\|
188	info->evtNum == 1038911933 \|\|
189	info->evtNum == 876658967 \|\|
190	info->evtNum == 966824024 \|\|
191	info->evtNum == 141954801 \|\|
192	info->evtNum == 160966858 \|\|
193	info->evtNum == 191231387 \|\|
194	info->evtNum == 66033190 \|\|
195	info->evtNum == 10347106 \|\|
196	info->evtNum == 107360878 \|\|
197	info->evtNum == 2554393033 \|\|
198	info->evtNum == 933807102 \|\|
199	info->evtNum == 1188043146 \|\|
200	info->evtNum == 559839432 \|\|
201	info->evtNum == 16706390 \|\|
202	info->evtNum == 65557571 \|\|
203	info->evtNum == 389185367 \|\|
204	info->evtNum == 2722114329 ) ) continue;
205	#endif
206
207
208	#ifdef BDTFAIL_THEIR_EVENTS
209	if( !( info->evtNum == 78213037 \|\|
210	info->evtNum == 876658967 ) ) continue;
211
212	#endif
213
214	unsigned evtfail = fails_HZZ4L_selection(ctrl, info, electronArr, muonArr, eventweight, passtuple);
215
216	#ifdef THEIR_EVENTS
217	if(evtfail!=0) cout << "Failed!" << endl;
218	#endif
219
220	h_evtfail->Fill( evtfail, eventweight );
221
222	count++;
223	if( !evtfail ) pass++;
224
225	} //end loop over data
226
227	if( ctrl.mc ) {
228	cout << "--------------" << endl;
229	cout << "denom: " << count << endl;
230	cout << "pass: " << pass << endl;
231	cout << "axe: " << (float)pass/count << endl;
232	cout << "--------------" << endl;
233	}
234
235	delete info;
236	delete electronArr;
237	delete muonArr;
238
239
240
241	//--------------------------------------------------------------------------------------------------------------
242	// Save Histograms;
243	//==============================================================================================================
244	// const char * ofname;
245	// if( strcmp( ctrl.outputfile.c_str(), "") ) {
246	// ofname = ctrl.outputfile.c_str();
247	// } else {
248	// ofname = "tmp.root";
249	// }
250
251	// TFile *file = new TFile(ofname, "RECREATE");
252	gROOT->cd(file->GetPath());
253	h_evt->Write();
254	h_evtfail->Write();
255	passtuple->Write();
256	file->Close();
257	delete file;
258
259	// map<TString,TMVA::Reader*>::iterator it;
260	// for(it=readers.begin(); it!=readers.end(); it++) delete (*it).second;
261
262	cerr << "done!" << endl;
263	}
264
265
266