ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/UNL/macros/MyHbbAnalyzer.C

Comparing UserCode/UNL/macros/MyHbbAnalyzer.C (file contents):
Revision 1.3 by malbouis, Mon Sep 5 18:18:12 2011 UTC vs.
Revision 1.5 by malbouis, Thu Oct 13 15:10:45 2011 UTC

#	Line 9 | Line 9
9		void MyHbbAnalyzer::Loop()
10		{
11		//   In a ROOT session, you can do:
12	<	//      Root > .L MyHbbAnalyzer.C
12	>	//      Root > .L MyHbbAnalyzer.C+
13		//      Root > MyHbbAnalyzer t
14		//      Root > t.GetEntry(12); // Fill t data members with entry number 12
15		//      Root > t.Show();       // Show values of entry 12
#	Line 37 | Line 37 | void MyHbbAnalyzer::Loop()
37
38		Long64_t nbytes = 0, nb = 0;
39
40	–	////////
41	–	// open output text file
42	–	////////
43	–	ofstream myfile;
44	–	myfile.open ("sync10K.txt");
40
41		////////
42		// get samples weight:
#	Line 55 | Line 50 | void MyHbbAnalyzer::Loop()
50		// initialize variables for significance
51		////////
52
58	–	// CSV vars:
59	–	// variables for significance:
60	–	// variables for significance:
61	–	double cut_CSV = 0, cut2_CSV = 0;
62	–	double cut_interval_CSV = 0.025;
63	–	int const max_cut_iter_CSV = 200;
64	–
65	–	int evt_count_CSV_zh[max_cut_iter_CSV][2];
66	–	//initializeArr(evt_count_CSV_zh, max_cut_iter_CSV);
67	–	int evt_count_CSV_zz[max_cut_iter_CSV][2];
68	–	//initializeArr(evt_count_CSV_zz, max_cut_iter_CSV);
69	–	int evt_count_CSV_z[max_cut_iter_CSV][2];
70	–	//initializeArr(evt_count_CSV_z, max_cut_iter_CSV);
71	–	// 2D significance:
72	–	int evt_count_CSV2D_zh[max_cut_iter_CSV][max_cut_iter_CSV];
73	–	//initializeArr2D(evt_count_CSV2D_zh, max_cut_iter_CSV);
74	–	int evt_count_CSV2D_zz[max_cut_iter_CSV][max_cut_iter_CSV];
75	–	//initializeArr2D(evt_count_CSV2D_zz, max_cut_iter_CSV);
76	–	int evt_count_CSV2D_z[max_cut_iter_CSV][max_cut_iter_CSV];
77	–	//initializeArr2D(evt_count_CSV2D_z, max_cut_iter_CSV);
78	–
79	–	for (int i = 0; i != max_cut_iter_CSV; i++){
80	–	for (int j = 0; j != max_cut_iter_CSV; j++){
81	–	evt_count_CSV2D_zh[i][j] = 0;
82	–	evt_count_CSV2D_zz[i][j] = 0;
83	–	evt_count_CSV2D_z[i][j] = 0;
84	–	if (j < 2){
85	–	evt_count_CSV_zh[i][j] = 0;
86	–	evt_count_CSV_zz[i][j] = 0;
87	–	evt_count_CSV_z[i][j] = 0;
88	–	}// j < 2
89	–	}// cuts
90	–	}// cuts
91	–
53		////////
54		// cut and count variables
55		////////
#	Line 101 | Line 62 | void MyHbbAnalyzer::Loop()
62		int allEvts_z = 0;
63		int preSelection_z = 0;
64		int pT_jj_z = 0, pT_Z_z = 0, CSV1_z = 0, CSV2_z = 0, DPhi_ZH_z = 0, N_aj_z = 0, M_jj_z = 0;
65	<
65	>
66		////////
67		// loop on entries
68		////////
69		pct = 0;
70		total_entries = fChain->GetEntries();
71
72	<	//nentries = 1000;
72	>	//nentries = 10000;
73		for (Long64_t jentry=0; jentry<nentries;jentry++) {
74		Long64_t ientry = LoadTree(jentry);
75		if (ientry < 0) break;
#	Line 123 | Line 84 | void MyHbbAnalyzer::Loop()
84		if (isZZ()) weight = ZZ_weight;
85		if (isDY()) weight = DY_weight;
86
126	–	//cout << "sample ZH: " << isZH() << " DY: " << isDY() << " ZZ: " << isZZ() << endl;
127	–	//cout << "weight = " << weight << endl;
128	–
87		// count events:
88		if (isZH()) allEvts++;
89		if (isZZ()) allEvts_zz++;
90		if (isDY()) allEvts_z++;
91		basicMuonDistributions("allEvts", weight);
92		basicJetDistributions("allEvts", weight);
93	<
93	>
94	>
95		// pre-selection cuts:
96		if (!preSelect()) continue;
97		if (isZH()) preSelection++;
#	Line 140 | Line 99 | void MyHbbAnalyzer::Loop()
99		if (isDY()) preSelection_z++;
100		basicMuonDistributions("preSel", weight);
101		basicJetDistributions("preSel", weight);
143	–	myfile << "passed pre-selection:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << "\n";
102
103		// pt_jj cut:
104		if (!Select_pTjj()) continue;
#	Line 149 | Line 107 | void MyHbbAnalyzer::Loop()
107		if (isDY()) pT_jj_z++;
108		basicMuonDistributions("pT_jj", weight);
109		basicJetDistributions("pT_jj", weight);
152	–	myfile << "passed pT(jj):" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
153	–
110
111		// pT_Z cut:
112		if (!Select_pTZ()) continue;
#	Line 159 | Line 115 | void MyHbbAnalyzer::Loop()
115		if (isDY()) pT_Z_z++;
116		basicMuonDistributions("pT_Z", weight);
117		basicJetDistributions("pT_Z", weight);
162	–	myfile << "passed pTZ:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
163	–
164	–	////////
165	–	// look at CSV1 optimization:
166	–	// fill arrays with discriminators after cuts
167	–	////////
168	–	int Jet1Idx, Jet2Idx;
169	–	jetIndices(Jet1Idx, Jet2Idx, bDisc_CSV);
170	–
171	–	for(int icut = 0; icut != max_cut_iter_CSV; icut++){
172	–	cut_CSV = icut*cut_interval_CSV;
173	–	if (bDisc_CSV[Jet1Idx] > cut_CSV){
174	–	if (isZH()) evt_count_CSV_zh[icut][0] = evt_count_CSV_zh[icut][0] + 1;
175	–	if (isZZ()) evt_count_CSV_zz[icut][0] = evt_count_CSV_zz[icut][0] + 1;
176	–	if (isDY()) evt_count_CSV_z[icut][0] = evt_count_CSV_z[icut][0] + 1;
177	–	}
178	–	}// cuts
118
119		// CSV1 cut:
120		if (!Select_CSV1()) continue;
#	Line 184 | Line 123 | void MyHbbAnalyzer::Loop()
123		if (isDY()) CSV1_z++;
124		basicMuonDistributions("CSV1", weight);
125		basicJetDistributions("CSV1", weight);
187	–	myfile << "passed CSV1:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
188	–
189	–	////////
190	–	// look at CSV2 optimization, AFTER CSV1 cut!
191	–	// fill arrays with discriminators after cuts
192	–	////////
193	–
194	–	for(int icut = 0; icut != max_cut_iter_CSV; icut++){
195	–	cut_CSV = icut*cut_interval_CSV;
196	–	if (bDisc_CSV[Jet2Idx] > cut_CSV){
197	–	if (isZH()) evt_count_CSV_zh[icut][1] = evt_count_CSV_zh[icut][1] + 1;
198	–	if (isZZ()) evt_count_CSV_zz[icut][1] = evt_count_CSV_zz[icut][1] + 1;
199	–	if (isDY()) evt_count_CSV_z[icut][1] = evt_count_CSV_z[icut][1] + 1;
200	–	}
201	–	}// cuts
126
203	–
127		// CSV2 cut:
128		if (!Select_CSV2()) continue;
129		if (isZH()) CSV2++;
#	Line 208 | Line 131 | void MyHbbAnalyzer::Loop()
131		if (isDY()) CSV2_z++;
132		basicMuonDistributions("CSV2", weight);
133		basicJetDistributions("CSV2", weight);
211	–	myfile << "passed CSV2:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
134
135		// DPhi cut:
136		if (!Select_DPhi()) continue;
#	Line 217 | Line 139 | void MyHbbAnalyzer::Loop()
139		if (isDY()) DPhi_ZH_z++;
140		basicMuonDistributions("DPhi", weight);
141		basicJetDistributions("DPhi", weight);
142	<	myfile << "passed DPhi:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
221	<
222	<
142	>
143		// Naj cut:
144		if (!Select_Naj()) continue;
145		if (isZH()) N_aj++;
#	Line 227 | Line 147 | void MyHbbAnalyzer::Loop()
147		if (isDY()) N_aj_z++;
148		basicMuonDistributions("Naj", weight);
149		basicJetDistributions("Naj", weight);
150	<	myfile << "passed Naj:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
231	<
232	<
150	>
151		// Mjj cut:
152		if (!Select_Mjj()) continue;
153		if (isZH()) M_jj++;
#	Line 237 | Line 155 | void MyHbbAnalyzer::Loop()
155		if (isDY()) M_jj_z++;
156		basicMuonDistributions("Mjj", weight);
157		basicJetDistributions("Mjj", weight);
240	–	myfile << "passed Mjj:" << runNumber << ":" << evtNumber << ":" << lumiBlock << ":" << jetPt[0] << ":" << jetEta[0] << ":" << muonPt[0] << ":" << muonEta[0] << ":" << computePtjj() << ":" << computePtZ() << ":" << computeDPhi() << endl;
158
242	–
159		}// loop on entries
160
245	–	//////
246	–	// write significance histos
247	–	////////
248	–	int nbins = max_cut_iter_CSV;
249	–	double binMax = (max_cut_iter_CSV*cut_interval_CSV)+cut_interval_CSV;
250	–	double binMin = cut_interval_CSV;
251	–
252	–	// fill cut histos:
253	–	for (int i = 0; i!=max_cut_iter_CSV; i++){
254	–	for (int ijet = 0; ijet != 2; ijet++){
255	–	string zh_csv_histo = Form("ZH_CSV_jet%i", ijet);
256	–	fillhisto(zh_csv_histo, i*cut_interval_CSV, ZH_weight*evt_count_CSV_zh[i][ijet], "", nbins, binMin, binMax);
257	–	string zz_csv_histo = Form("ZZ_CSV_jet%i", ijet);
258	–	fillhisto(zz_csv_histo, i*cut_interval_CSV, ZZ_weight*evt_count_CSV_zz[i][ijet], "", nbins, binMin, binMax);
259	–	string z_csv_histo = Form("DY_CSV_jet%i", ijet);
260	–	fillhisto(z_csv_histo, i*cut_interval_CSV, DY_weight*evt_count_CSV_z[i][ijet], "", nbins, binMin, binMax);
261	–	string sqrtBackground_CSV = Form("sqrtBackground_CSV_jet%i", ijet);
262	–	fillhisto(sqrtBackground_CSV, i*cut_interval_CSV, sqrt((DY_weight*evt_count_CSV_z[i][ijet])+(ZZ_weight*evt_count_CSV_zz[i][ijet])), "", nbins, binMin, binMax);
263	–	}// jets
264	–	}// cuts
265	–
161
162		////////
163		// write histos to file:
#	Line 275 | Line 170 | void MyHbbAnalyzer::Loop()
170		(*it2).second->Write();
171		delete (*it2).second;
172		}
173	<
173	>
174		f->Write();
175		f->Close();
176
282	–	// close output text file
283	–	myfile.close();
177
178		////////
179		// print stats
#	Line 288 | Line 181 | void MyHbbAnalyzer::Loop()
181		cout << endl;
182		cout << endl;
183		cout << "**** stats (no normalization) ****" << endl;
184	<	cout << "allEvts: " << allEvts <<  endl;
185	<	cout << "preSelection: " << preSelection << endl;
186	<	cout << "pT_jj: " << pT_jj << endl;
187	<	cout << "pT_Z: " << pT_Z << endl;
188	<	cout << "CSV1: " << CSV1 << endl;
189	<	cout << "CSV2: " << CSV2 << endl;
190	<	cout << "DPhi: " << DPhi_ZH << endl;
191	<	cout << "Naj: " << N_aj << endl;
192	<	cout << "Mjj: " << M_jj << endl;
184	>	cout << "allEvts, ZH: " << allEvts << " , ZZ: " << allEvts_zz << ", DY: " << allEvts_z << endl;
185	>	cout << "preSelection, ZH: " << preSelection << " , ZZ: " << preSelection_zz << ", DY: " << preSelection_z <<  endl;
186	>	cout << "pT_jj, ZH: " << pT_jj << " , ZZ: " << pT_jj_zz << ", DY: " << pT_jj_z <<  endl;
187	>	cout << "pT_Z, ZH: " << pT_Z << " , ZZ: " << pT_Z_zz << ", DY: " << pT_Z_z <<  endl;
188	>	cout << "CSV1, ZH: " << CSV1 << " , ZZ: " << CSV1_zz << ", DY: " << CSV1_z <<  endl;
189	>	cout << "CSV2, ZH: " << CSV2 << " , ZZ: " << CSV2_zz << ", DY: " << CSV2_z <<  endl;
190	>	cout << "DPhi, ZH: " << DPhi_ZH << " , ZZ: " << DPhi_ZH_zz << ", DY: " << DPhi_ZH_z <<  endl;
191	>	cout << "Naj, ZH: " << N_aj << " , ZZ: " << N_aj_zz << ", DY: " << N_aj_z <<  endl;
192	>	cout << "Mjj, ZH: " << M_jj << " , ZZ: " << M_jj_zz << ", DY: " << M_jj_z <<  endl;
193		cout << endl;
194		cout << endl;
195		cout << "**** stats normalized to 10 fb-1 ****" << endl;
#	Line 335 | Line 228 | void MyHbbAnalyzer::Loop()
228		cout << endl;
229		cout << "**** TOTAL SIGNIFICANCE before Mjj cut ****" << endl;
230		cout << "S = S/sqrt{B} = " << ((double)N_aj*ZH_weight)/sqrt((N_aj_zz*ZZ_weight) + (N_aj_z*DY_weight)) << endl;
231	<
231	>	cout << endl;
232	>	cout << endl;
233
234		}// Loop
235		bool MyHbbAnalyzer::preSelect(){
#	Line 344 | Line 238 | bool MyHbbAnalyzer::preSelect(){
238		if (!HLT_IsoMu17_v5_trig) return false;
239		if (nMuons < 2) passedPreSel = false;
240		if (nJets < 2) passedPreSel = false;
241	+	int firstJetIdx, secondJetIdx;
242	+	jetIndices(firstJetIdx, secondJetIdx, bDisc_CSV);
243	+	// reject events with less than 2 CSV tags
244	+	if (firstJetIdx == -1 || secondJetIdx == -1) return false;
245	+	double mz = computeMmumu();
246	+	if (mz <= 75. || mz >= 105.) passedPreSel = false;
247		return passedPreSel;
248		}// preSelect
249
#	Line 393 | Line 293 | bool MyHbbAnalyzer::Select_DPhi(){
293
294		bool MyHbbAnalyzer::Select_Naj(){
295		bool passedNaj = true;
296	<	if (nJets > 3) passedNaj = false;
296	>	////if (nJets > 3) passedNaj = false;
297	>	int bjets = 0;
298	>	for(int j = 0; j != nJets; ++j){
299	>	if (bDisc_CSV[j] < 0.244) continue;
300	>	bjets++;
301	>	}// jet loop
302	>	if (bjets > 3) passedNaj = false;
303		return passedNaj;
304		}// Naj
305
#	Line 403 | Line 309 | bool  MyHbbAnalyzer::Select_Mjj(){
309		if (Mjj < 95 || Mjj > 125) passedMjj = false;
310		return passedMjj;
311		}
312	<	void MyHbbAnalyzer::jetIndices(int &index1, int &index2, double bDisc[5]){
312	>	void MyHbbAnalyzer::jetIndices(int &index1, int &index2, double bDisc[9]){
313		int fJetIdx = -1, sJetIdx = -1;
314		double highestBdiscSum = 0;
315		for (int i = 0; i != nJets; ++i){
#	Line 447 | Line 353 | void MyHbbAnalyzer::basicMuonDistributio
353		if (isDY()) ph_process = "DY";
354
355		for (int i=0; i != nMuons; ++i) {
356	<	string muonpT = "muonPt_" + ph_process + cut;
357	<	string muoneta = "muonEta_" + ph_process + cut;
358	<	string muonphi = "muonPhi_" + ph_process + cut;
356	>	string muonpT = Form("muonPt_%i", i) + ph_process + cut;
357	>	string muoneta = Form("muonEta_%i", i) + ph_process + cut;
358	>	string muonphi = Form("muonPhi_%i", i) + ph_process + cut;
359
360		//cout << muonpT << endl;
361
#	Line 457 | Line 363 | void MyHbbAnalyzer::basicMuonDistributio
363		fillhisto(muoneta, muonEta[i], ph_weight, "#mu #eta", 100, -5, 5);
364		fillhisto(muonphi, muonPhi[i], ph_weight, "#mu #phi", 80, -4, 4);
365		}
366	+
367	+	double mz = computeMmumu();
368	+	string Mz = "MZ+" + ph_process + cut;
369	+	fillhisto(Mz, mz, ph_weight, "", 400, 0, 400);
370	+
371		}// basicMuonDistributions
372
373		void MyHbbAnalyzer::basicJetDistributions(string cut, double ph_weight){
#	Line 468 | Line 379 | void MyHbbAnalyzer::basicJetDistribution
379		if (isDY()){ ph_process = "DY";}
380
381		for (int i=0; i != nJets; ++i) {
382	<	string jetpT = "jetPt_" + ph_process + cut;
383	<	string jeteta = "jetEta_" + ph_process + cut;
384	<	string jetphi = "jetPhi_" + ph_process + cut;
474	<
475	<	//cout << muonpT << endl;
382	>	string jetpT = Form("jetPt_%i", i) + ph_process + cut;
383	>	string jeteta = Form("jetEta_%i", i) + ph_process + cut;
384	>	string jetphi = Form("jetPhi_%i", i) + ph_process + cut;
385
386		fillhisto(jetpT, jetPt[i], ph_weight, "jet pT", 400, 0, 400);
387		fillhisto(jeteta, jetEta[i], ph_weight, "jet #eta", 100, -5, 5);
#	Line 483 | Line 392 | void MyHbbAnalyzer::basicJetDistribution
392		double Mjj = computeMjj();
393		string mjj = "Mjj_" + ph_process + cut;
394		fillhisto(mjj, Mjj, ph_weight, "", 1000, 0, 1000);
395	+
396		}// basicJetDistributions
397
398		double MyHbbAnalyzer::computeMjj(){
399		int firstJetIdx, secondJetIdx;
400		jetIndices(firstJetIdx, secondJetIdx, bDisc_CSV);
401	<
401	>	// reject event with less than 2 CSV tags
402	>	if (firstJetIdx == -1 || secondJetIdx == -1) return false;
403		TLorentzVector firstJet_p4 = TLorentzVector(jetPx[firstJetIdx], jetPy[firstJetIdx], jetPz[firstJetIdx], jetP[firstJetIdx]);
404		TLorentzVector secondJet_p4 = TLorentzVector(jetPx[secondJetIdx], jetPy[secondJetIdx], jetPz[secondJetIdx], jetP[secondJetIdx]);
405
#	Line 497 | Line 408 | double MyHbbAnalyzer::computeMjj(){
408		return jj.M();
409		}// Mjj
410
411	+	double MyHbbAnalyzer::computeMmumu(){
412	+
413	+	double Mz = 0.;
414	+	TLorentzVector firstMu_p4 = TLorentzVector(muonPx[0], muonPy[0], muonPz[0], muonP[0]);
415	+	TLorentzVector secondMu_p4 = TLorentzVector(muonPx[1], muonPy[1], muonPz[1], muonP[1]);
416	+	TLorentzVector Z = firstMu_p4 + secondMu_p4;
417	+	Mz = Z.M();
418	+
419	+	return Mz;
420	+
421	+	}// Mmumu
422	+
423		double MyHbbAnalyzer::computePtjj(){
424		double ptjj = 0.;
425		// get jet indices that maximize the b-tag discriminant
#	Line 531 | Line 454 | double MyHbbAnalyzer::computeDPhi(){
454		TLorentzVector Z = firstMu_p4 + secondMu_p4;
455		int firstJetIdx, secondJetIdx;
456		jetIndices(firstJetIdx, secondJetIdx, bDisc_CSV);
457	+	// reject event with less than 2 CSV tags
458	+	if (firstJetIdx == -1 || secondJetIdx == -1) return false;
459		TLorentzVector firstJet_p4 = TLorentzVector(jetPx[firstJetIdx], jetPy[firstJetIdx], jetPz[firstJetIdx], jetP[firstJetIdx]);
460		TLorentzVector secondJet_p4 = TLorentzVector(jetPx[secondJetIdx], jetPy[secondJetIdx], jetPz[secondJetIdx], jetP[secondJetIdx]);
461		TLorentzVector jj = firstJet_p4 + secondJet_p4;
#	Line 539 | Line 464 | double MyHbbAnalyzer::computeDPhi(){
464		return dphi;
465
466		}// dPhi
467	+	double MyHbbAnalyzer::DeltaPhi(const double phi1, const double phi2){
468	+	double delta = phi1-phi2;
469	+	if (delta > TMath::Pi()) delta -= 2.0*TMath::Pi();
470	+	if (delta < (-TMath::Pi())) delta += 2.0*TMath::Pi();
471	+	return fabs(delta);
472	+	}
473
474
475

Diff Legend

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