ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitHzz4l/Selection/src/applySelection.cc

Comparing UserCode/MitHzz4l/Selection/src/applySelection.cc (file contents):
Revision 1.2 by dkralph, Wed Sep 14 12:11:57 2011 UTC vs.
Revision 1.12 by khahn, Sun Oct 30 15:35:25 2011 UTC

#	Line 29 | Line 29 | using namespace std;
29		#include <TLorentzVector.h>         // 4-vector class
30		#include <TVector3.h>               // 3D vector class
31
32	+
33		//
34		// ntuple format headers
35		//
#	Line 37 | Line 38 | using namespace std;
38		#include "TElectron.hh"
39		#include "TMuon.hh"
40		#include "TJet.hh"
41	+	#include "SiMVAElectronSelection.h"
42		#include "RunLumiRangeMap.h"
43	+	#include "EfficiencyWeightsInterface.h"
44
45		//
46		// utility headers
#	Line 46 | Line 49 | using namespace std;
49		#include "SampleWeight.h"
50		#include "Selection.h"
51		#include "HZZCiCElectronSelection.h"
52	+	#include "HZZLikelihoodElectronSelection.h"
53	+	#include "HZZBDTElectronSelection.h"
54		#include "SampleWeight.h"
55
56	+	//#define BDTFAIL_THEIR_EVENTS
57		//#define THEIR_EVENTS
58
59	+	Double_t etaFromY(Double_t pt, Double_t y, Double_t phi, Double_t m)
60	+	{
61	+	Double_t a  = (1+exp(2*y))/(exp(2*y)-1); // intermediate term
62	+	if(a*a<1) { cout << "a too small" << endl; assert(0); }
63	+	Double_t E  = sqrt( a*a*(pt*pt+m*m)/(a*a-1) );
64	+	Double_t pz = E*E - pt*pt - 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(pt*cos(phi),pt*sin(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		//=== MAIN =================================================================================================
90		int main(int argc, char** argv)
91		{
#	Line 72 | Line 108 | int main(int argc, char** argv)
108		}
109		ctrl.dump();
110
111	+	map<string,double> entrymap; // number of unskimmed entries in each file
112	+
113		//
114		// File I/O
115		//--------------------------------------------------------------------------------------------------------------
#	Line 79 | Line 117 | int main(int argc, char** argv)
117		ifstream f(ctrl.inputfile.c_str());
118		string fname;
119		while (f >> fname) {
120	<	if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
120	>	if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines
121		cout << "adding inputfile : " << fname.c_str() << endl;
122	+	entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str());
123		chain->AddFile(fname.c_str());
124		}
125
#	Line 90 | Line 129 | int main(int argc, char** argv)
129		//
130		// Setup
131		//--------------------------------------------------------------------------------------------------------------
132	+	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		TH1F * h_evt = new TH1F("hevt", "hevt", 2, 0, 2 );
140	<	TH1F * h_evtfail = new TH1F("hevtfail", "hevtfail", 100, 0, 100 );
141	<	TNtuple * passtuple = new TNtuple( "passtuple", "passtuple", "run:evt:lumi:channel:mZ1:mZ2:m4l:pt4l:w" );
140	>	TH1F * h_evtfail = new TH1F("hevtfail", "hevtfail", 1024, 0, 1024 );
141	>	TTree * passtuple = new TTree("passtuple", "passtuple" );
142	>	unsigned run, evt, lumi, channel, gchannel;
143	>	float mZ1, mZ2, m4l, pt4l;
144	>	unsigned tZ1, tZ2;
145	>	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	>	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	>	passtuple->Branch("tZ1", &tZ1);
160	>	passtuple->Branch("tZ2", &tZ2);
161	>	passtuple->Branch("m4l", &m4l);
162	>	passtuple->Branch("pt4l", &pt4l);
163	>	passtuple->Branch("w", &w);
164	>
165		initCiCSelection();
166	+	if(ctrl.eleSele=="lik") initLikSelection();
167	+	if(ctrl.eleSele=="bdt") initBDTSelection();
168	+	// if(ctrl.eleSele=="si" ) initSiMVAElectronSelection();
169		initRunLumiRangeMap();
170	+	initEfficiencyWeights();
171	+
172	+
173
174
175		//
176		// Loop
177		//--------------------------------------------------------------------------------------------------------------
178	+	if(ctrl.fakeScheme!="none") {cout << "wrong executable!" << endl; assert(0);}
179
180		//
181		// Access samples and fill histograms
#	Line 108 | Line 184 | int main(int argc, char** argv)
184
185		// Data structures to store info from TTrees
186		mithep::TEventInfo *info    = new mithep::TEventInfo();
187	+	mithep::TGenInfo    *ginfo  = new mithep::TGenInfo();
188		TClonesArray *electronArr   = new TClonesArray("mithep::TElectron");
189		TClonesArray *muonArr       = new TClonesArray("mithep::TMuon");
190
191
115	–	TBranch *infoBr;
116	–	TBranch *electronBr;
117	–	TBranch *muonBr;
118	–
119	–	//   chain->SetBranchStatus("*",0); //disable all branches
120	–	//   chain->SetBranchStatus("Info", 1);
121	–	//   chain->SetBranchStatus("Muon", 1);
122	–	//   chain->SetBranchStatus("Electron", 1);
123	–
192		chain->SetBranchAddress("Info",       &info);
193		chain->SetBranchAddress("Electron", &electronArr);
194		chain->SetBranchAddress("Muon", &muonArr);
195	+	if( ctrl.mc ) { chain->SetBranchAddress("Gen",        &ginfo);}
196
197
198	+	int count=0, pass=0;
199	+	float denom[3]={0.,0.,0.};
200	+	float numer[3]={0.,0.,0.};
201	+	float numercorr[3]={0.,0.,0.};
202
203	<	cout << "nEntries: " << chain->GetEntries() << endl;
204	<	for(UInt_t ientry=0; ientry<chain->GetEntries(); ientry++) {
205	<
203	>	UInt_t imax = chain->GetEntries();
204	>	cout << "nEntries: " << imax << endl;
205	>
206	>
207	>	for(UInt_t ientry=0; ientry<imax; ientry++) {
208		chain->GetEntry(ientry);
209	+	if(!(ientry%100000)) cerr << "entry: " << ientry << endl;
210
211		// get event weight for this file
212	<	double xs=1,eventweight=1;
213	<	if(ctrl.mc) eventweight = xstab.Get(getname(chain).c_str()) / chain->GetEntries();
212	>	double eventweight=1;
213	>	if(ctrl.mc) eventweight = getWeight(xstab,entrymap,chain);
214
215		#ifdef THEIR_EVENTS
216		if( !( info->evtNum == 504867308  ||
217		info->evtNum == 368148849  ||
218	<	info->evtNum == 129514273  ||
218	>	//      info->evtNum == 129514273  ||
219		info->evtNum == 286336207  ||
220		info->evtNum == 344708580  ||
221		info->evtNum == 30998576   ||
#	Line 147 | Line 223 | int main(int argc, char** argv)
223		info->evtNum == 394010457  ||
224		info->evtNum == 917379387  ||
225		info->evtNum == 78213037   ||
226	<	info->evtNum == 337493970  ||
226	>	info->evtNum == 862270386  ||
227	>	info->evtNum == 337493970  ||  // not baseline anymore?
228		info->evtNum == 1491724484 ||
229		info->evtNum == 480301165  ||
230		info->evtNum == 1038911933 ||
#	Line 158 | Line 235 | int main(int argc, char** argv)
235		info->evtNum == 191231387  ||
236		info->evtNum == 66033190   ||
237		info->evtNum == 10347106   ||
238	<	info->evtNum == 107360878  ) ) continue;
238	>	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		#endif
248	+
249
250	<	unsigned evtfail = fails_HZZ4L_selection(ctrl, info, electronArr, muonArr, eventweight, passtuple );
250	>	#ifdef BDTFAIL_THEIR_EVENTS
251	>	if( !( info->evtNum == 78213037   ||
252	>	info->evtNum == 876658967  ) ) continue;
253	>
254	>	#endif
255	>
256	>	unsigned evtfail = fails_HZZ4L_selection(ctrl, info, electronArr, muonArr, eventweight, passtuple, NULL);
257	>	//    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	>
265	>	#ifdef THEIR_EVENTS
266	>	if(evtfail!=0) cout << "Failed!" << endl;
267	>	#endif
268	>
269		h_evtfail->Fill( evtfail, eventweight );
270	<	cout << endl << endl;
270	>
271	>	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	>	count++;
293	>	if( !evtfail ) pass++;
294
295		} //end loop over data
296
297	<
297	>	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
305		delete info;
306		delete electronArr;
#	Line 178 | Line 311 | int main(int argc, char** argv)
311		//--------------------------------------------------------------------------------------------------------------
312		// Save Histograms;
313		//==============================================================================================================
314	<	const char * ofname;
315	<	if( strcmp( ctrl.outputfile.c_str(), "") ) {
316	<	ofname = ctrl.outputfile.c_str();
317	<	} else {
318	<	ofname = "tmp.root";
319	<	}
314	>	// const char * ofname;
315	>	// if( strcmp( ctrl.outputfile.c_str(), "") ) {
316	>	//   ofname = ctrl.outputfile.c_str();
317	>	// } else {
318	>	//   ofname = "tmp.root";
319	>	// }
320
321	<	TFile *file = new TFile(ofname, "RECREATE");
321	>	// TFile *file = new TFile(ofname, "RECREATE");
322	>	gROOT->cd(file->GetPath());
323		h_evt->Write();
324		h_evtfail->Write();
325		passtuple->Write();
326		file->Close();
327		delete file;
328
329	+	// map<TString,TMVA::Reader*>::iterator it;
330	+	// for(it=readers.begin(); it!=readers.end(); it++) delete (*it).second;
331	+
332	+
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		cerr << "done!" << endl;
363		}
364

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines