Context Navigation

Back to Ticket #319

Ticket #319: analyzeLambdab.C

File analyzeLambdab.C, 33.7 KB (added by Michele Selvaggi, 10 years ago)

Line
1
2	#include "TH1.h"
3	#include "THStack.h"
4	#include "TSystem.h"
5	#include "TLegend.h"
6	#include "TPaveText.h"
7	#include "TLorentzVector.h"
8
9	#include <utility> // needed to use "pair"
10	#include <vector> // needed to use vector of TLorentzVectors
11
12	//------------------------------------------------------------------------------
13
14	struct MyPlots
15	{
16	TH1 *cutFlow;
17
18	TH1 *fMissingET;
19	TH1 *fMuonPT;
20	TH1 *fSoftMuonPT;
21	TH1 *fMuonInJetPT;
22	TH1 *fMuonDxy;
23	TH1 *fSoftMuonDxy;
24	TH1 *fElectronPT;
25
26	TH1 *NHardLeptons;
27	TH1 *NJets;
28	TH1 *NBtags;
29	TH1 *NSoftMuons;
30	TH1 *NMuonsInJets;
31	TH1 *NTaggedSoftMuons;
32
33	TH1 *fDRPlusMuonBQuark;
34	TH1 *fDRPlusMuonAntiBQuark;
35	TH1 *fDRMinusMuonBQuark;
36	TH1 *fDRMinusMuonAntiBQuark;
37
38	TH1 *fPtLambdabOverB;
39	TH1 *fPxLambdabOverB;
40	TH1 *fPyLambdabOverB;
41	TH1 *fPzLambdabOverB;
42	TH1 *fELambdabOverB;
43	TH1 *fPtLambdabOverB_selected;
44	TH1 *fPxLambdabOverB_selected;
45	TH1 *fPyLambdabOverB_selected;
46	TH1 *fPzLambdabOverB_selected;
47	TH1 *fELambdabOverB_selected;
48
49	TH1 *fPtJetOverB;
50	TH1 *fPxJetOverB;
51	TH1 *fPyJetOverB;
52	TH1 *fPzJetOverB;
53	TH1 *fEJetOverB;
54
55	TH1 *fPtJetOverLambdab;
56	TH1 *fPxJetOverLambdab;
57	TH1 *fPyJetOverLambdab;
58	TH1 *fPzJetOverLambdab;
59	TH1 *fEJetOverLambdab;
60
61	TH1 *fPtCorrJetOverLambdab;
62	TH1 *fPxCorrJetOverLambdab;
63	TH1 *fPyCorrJetOverLambdab;
64	TH1 *fPzCorrJetOverLambdab;
65	TH1 *fECorrJetOverLambdab;
66
67	TH1 *fDeltaEtaJetLambdab;
68	TH1 *fDeltaPhiJetLambdab;
69
70	TH1 *fDeltaEtaCorrJetLambdab;
71	TH1 *fDeltaPhiCorrJetLambdab;
72
73	TH1 *fVertexAllParticles;
74	TH1 *fVertexLambdabDaughters;
75	};
76
77	//------------------------------------------------------------------------------
78
79	class ExRootResult;
80	class ExRootTreeReader;
81
82	//------------------------------------------------------------------------------
83
84	pair<Int_t,Double_t> findClosestObject(vector<Double_t> drV) {
85	Int_t index=0;
86	Double_t dRmin = 9999.;
87
88	for (Int_t i=0; i<drV.size(); i++) {
89	if (drV[i]<dRmin) {
90	index=i;
91	dRmin=drV[i];
92	}
93	}
94
95	return make_pair(index,dRmin);
96	}
97
98
99	void BookHistograms(ExRootResult result, MyPlots plots)
100	{
101	THStack *stack;
102	TLegend *legend;
103	TPaveText *comment;
104
105	// book more histograms
106
107	plots->cutFlow = result->AddHist1D(
108	"cut_flow", "Cut number",
109	"Cut number", "number of events",
110	4, 0., 4.);
111
112	plots->fMuonPT = result->AddHist1D(
113	"muon_pt", "muon P_{T}",
114	"muon P_{T}, GeV/c", "number of muons",
115	50, 0.0, 100.0);
116
117	plots->fSoftMuonPT = result->AddHist1D(
118	"soft_muon_pt", "soft muon P_{T}",
119	"muon P_{T}, GeV/c", "number of soft muons",
120	50, 0.0, 100.0);
121
122	plots->fMuonInJetPT = result->AddHist1D(
123	"muon_in_jet_pt", "muon in jet P_{T}",
124	"muon P_{T}, GeV/c", "number of muons in jets",
125	50, 0.0, 100.0);
126
127	plots->fMuonDxy = result->AddHist1D(
128	"iso_muon_dxy", "muon D_{xy}",
129	"iso muon I.P., cm", "number of iso muons",
130	50, -2.0, 2.0);
131
132	plots->fSoftMuonDxy = result->AddHist1D(
133	"soft_muon_dxy", "soft muon D_{xy}",
134	"muon I.P., cm", "number of soft muons",
135	50, -2.0, 2.0);
136
137	plots->fElectronPT = result->AddHist1D(
138	"electron_pt", "electron P_{T}",
139	"electron P_{T}, GeV/c", "number of electrons",
140	50, 0.0, 100.0);
141
142	plots->NHardLeptons = result->AddHist1D("n_hard_leptons", "number of hard leptons",
143	"number of hard leptons", "number of events",
144	5,0.,5.);
145
146	plots->NJets = result->AddHist1D("n_jets", "number of jets",
147	"number of jets", "number of events",
148	10,0.,10.);
149
150	plots->NBtags = result->AddHist1D("n_btags", "number of b-tagged jets",
151	"number of b-tagged jets", "number of events",
152	5,0.,5.);
153
154	plots->NSoftMuons = result->AddHist1D("n_soft_muons", "number of soft muons",
155	"number of soft muons", "number of events",
156	5,0.,5.);
157	plots->NMuonsInJets = result->AddHist1D("n_muons_in_jets", "number of muons in jets",
158	"number of muons in jets", "number of events",
159	5,0.,5.);
160
161	plots->NTaggedSoftMuons = result->AddHist1D("n_tagged_soft_muons", "number of tagged soft muons",
162	"number of tagged soft muons", "number of events",
163	5,0.,5.);
164
165	plots->fMissingET = result->AddHist1D("missing_et", "Missing E_{T}",
166	"Missing E_{T}, GeV", "number of events",
167	60, 0.0, 150.0);
168
169	// DR(mu,b)
170	plots->fDRPlusMuonBQuark = result->AddHist1D("dr_plusmuon_bquark", "DR(soft mu, b)",
171	"DR(soft mu, b)", "number of events",
172	50,0.,6.3);
173	plots->fDRPlusMuonAntiBQuark = result->AddHist1D("dr_plusmuon_antibquark", "DR(soft mu, b)",
174	"DR(soft mu, b)", "number of events",
175	50,0.,6.3);
176	plots->fDRMinusMuonBQuark = result->AddHist1D("dr_minusmuon_bquark", "DR(soft mu, b)",
177	"DR(soft mu, b)", "number of events",
178	50,0.,6.3);
179	plots->fDRMinusMuonAntiBQuark = result->AddHist1D("dr_minusmuon_antibquark", "DR(soft mu, b)",
180	"DR(soft mu, b)", "number of events",
181	50,0.,6.3);
182
183	// lambda_b over b
184	plots->fPtLambdabOverB = result->AddHist1D("Pt_ratio_lambdab_bquark", "Pt(lambda_b)/Pt(b)",
185	"Pt(lambda_b)/Pt(b)", "number of lambda_b's",
186	50,0.,2.);
187	plots->fPxLambdabOverB = result->AddHist1D("Px_ratio_lambdab_bquark", "Px(lambda_b)/Px(b)",
188	"Px(lambda_b)/Px(b)", "number of lambda_b's",
189	50,0.,2.);
190	plots->fPyLambdabOverB = result->AddHist1D("Py_ratio_lambdab_bquark", "Py(lambda_b)/Py(b)",
191	"Py(lambda_b)/Py(b)", "number of lambda_b's",
192	50,0.,2.);
193	plots->fPzLambdabOverB = result->AddHist1D("Pz_ratio_lambdab_bquark", "Pz(lambda_b)/Pz(b)",
194	"Pz(lambda_b)/Pz(b)", "number of lambda_b's",
195	50,0.,2.);
196	plots->fELambdabOverB = result->AddHist1D("E_ratio_lambdab_bquark", "E(lambda_b)/E(b)",
197	"E(lambda_b)/E(b)", "number of lambda_b's",
198	50,0.,2.);
199
200	plots->fELambdabOverB->SetStats();
201
202
203	plots->fPtLambdabOverB_selected = result->AddHist1D("Pt_ratio_lambdab_bquark_selected", "Pt(lambda_b)/Pt(b)",
204	"Pt(lambda_b)/Pt(b)", "number of lambda_b's",
205	50,0.,2.);
206	plots->fPxLambdabOverB_selected = result->AddHist1D("Px_ratio_lambdab_bquark_selected", "Px(lambda_b)/Px(b)",
207	"Px(lambda_b)/Px(b)", "number of lambda_b's",
208	50,0.,2.);
209	plots->fPyLambdabOverB_selected = result->AddHist1D("Py_ratio_lambdab_bquark_selected", "Py(lambda_b)/Py(b)",
210	"Py(lambda_b)/Py(b)", "number of lambda_b's",
211	50,0.,2.);
212	plots->fPzLambdabOverB_selected = result->AddHist1D("Pz_ratio_lambdab_bquark_selected", "Pz(lambda_b)/Pz(b)",
213	"Pz(lambda_b)/Pz(b)", "number of lambda_b's",
214	50,0.,2.);
215	plots->fELambdabOverB_selected = result->AddHist1D("E_ratio_lambdab_bquark_selected", "E(lambda_b)/E(b)",
216	"E(lambda_b)/E(b)", "number of lambda_b's",
217	50,0.,2.);
218
219	plots->fELambdabOverB_selected->SetStats();
220
221	// jet over b
222	plots->fPtJetOverB = result->AddHist1D("Pt_ratio_jet_bquark", "Pt(jet)/Pt(b)",
223	"Pt(jet)/Pt(b)", "number of jets",
224	50,0.,2.);
225	plots->fPxJetOverB = result->AddHist1D("Px_ratio_jet_bquark", "Px(jet)/Px(b)",
226	"Px(jet)/Px(b)", "number of jets",
227	50,0.,2.);
228	plots->fPyJetOverB = result->AddHist1D("Py_ratio_jet_bquark", "Py(jet)/Py(b)",
229	"Py(jet)/Py(b)", "number of jets",
230	50,0.,2.);
231	plots->fPzJetOverB = result->AddHist1D("Pz_ratio_jet_bquark", "Pz(jet)/Pz(b)",
232	"Pz(jet)/Pz(b)", "number of jets",
233	50,0.,2.);
234	plots->fEJetOverB = result->AddHist1D("E_ratio_jet_bquark", "E(jet)/E(b)",
235	"E(jet)/E(b)", "number of jets",
236	50,0.,2.);
237
238	// jet over lambda_b
239	plots->fPtJetOverLambdab = result->AddHist1D("Pt_ratio_jet_lambdab", "Pt(jet)/Pt(lambda_b)",
240	"Pt(jet)/Pt(lambda_b)", "number of jets",
241	50,0.,2.);
242	plots->fPxJetOverLambdab = result->AddHist1D("Px_ratio_jet_lambdab", "Px(jet)/Px(lambda_b)",
243	"Px(jet)/Px(lambda_b)", "number of jets",
244	50,0.,2.);
245	plots->fPyJetOverLambdab = result->AddHist1D("Py_ratio_jet_lambdab", "Py(jet)/Py(lambda_b)",
246	"Py(jet)/Py(lambda_b)", "number of jets",
247	50,0.,2.);
248	plots->fPzJetOverLambdab = result->AddHist1D("Pz_ratio_jet_lambdab", "Pz(jet)/Pz(lambda_b)",
249	"Pz(jet)/Pz(lambda_b)", "number of jets",
250	50,0.,2.);
251	plots->fEJetOverLambdab = result->AddHist1D("E_ratio_jet_lambdab", "E(jet)/E(lambda_b)",
252	"E(jet)/E(lambda_b)", "number of jets",
253	50,0.,2.);
254
255	plots->fDeltaEtaJetLambdab = result->AddHist1D("DeltaEta_jet_lambdab", "#eta(jet)-#eta(lambda_b)",
256	"#eta(jet)-#eta(lambda_b)", "number of jets",
257	50,-3.,3.);
258	plots->fDeltaPhiJetLambdab = result->AddHist1D("DeltaPhi_jet_lambdab", "#phi(jet)-#phi(lambda_b)",
259	"#phi(jet)-#phi(lambda_b)", "number of jets",
260	50,-3.,3.);
261
262	// corrected jet over lambda_b
263	plots->fPtCorrJetOverLambdab = result->AddHist1D("Pt_ratio_corrjet_lambdab", "Pt(jet)/Pt(lambda_b)",
264	"Pt(jet)/Pt(lambda_b)", "number of jets",
265	50,0.,2.);
266	plots->fPxCorrJetOverLambdab = result->AddHist1D("Px_ratio_corrjet_lambdab", "Px(jet)/Px(lambda_b)",
267	"Px(jet)/Px(lambda_b)", "number of jets",
268	50,0.,2.);
269	plots->fPyCorrJetOverLambdab = result->AddHist1D("Py_ratio_corrjet_lambdab", "Py(jet)/Py(lambda_b)",
270	"Py(jet)/Py(lambda_b)", "number of jets",
271	50,0.,2.);
272	plots->fPzCorrJetOverLambdab = result->AddHist1D("Pz_ratio_corrjet_lambdab", "Pz(jet)/Pz(lambda_b)",
273	"Pz(jet)/Pz(lambda_b)", "number of jets",
274	50,0.,2.);
275	plots->fECorrJetOverLambdab = result->AddHist1D("E_ratio_corrjet_lambdab", "E(jet)/E(lambda_b)",
276	"E(jet)/E(lambda_b)", "number of jets",
277	50,0.,2.);
278
279	plots->fDeltaEtaCorrJetLambdab = result->AddHist1D("DeltaEta_corrjet_lambdab", "#eta(jet)-#eta(lambda_b)",
280	"#eta(jet)-#eta(lambda_b)", "number of jets",
281	50,-3.,3.);
282	plots->fDeltaPhiCorrJetLambdab = result->AddHist1D("DeltaPhi_corrjet_lambdab", "#phi(jet)-#phi(lambda_b)",
283	"#phi(jet)-#phi(lambda_b)", "number of jets",
284	50,-3.,3.);
285
286	// distance from vertex
287	plots->fVertexAllParticles = result->AddHist1D("Vertex_all_particles", "vertex distance",
288	"vertex distance", "number of particles",
289	50,0.,200.);
290	plots->fVertexLambdabDaughters = result->AddHist1D("Vertex_lambdab_daughters", "vertex distance",
291	"vertex distance", "number of particles",
292	50,0.,200.);
293
294
295	// book 1 stack of 2 histograms
296
297	plots->fMuonDxy->SetLineColor(kRed);
298	plots->fSoftMuonDxy->SetLineColor(kBlue);
299	stack = result->AddHistStack("muon_dxy_all", "iso and soft muons D_{xy}");
300	stack->Add(plots->fMuonDxy);
301	stack->Add(plots->fSoftMuonDxy);
302
303	// book legend for stack of 2 histograms
304
305	legend = result->AddLegend(0.72, 0.86, 0.98, 0.98);
306	legend->AddEntry(plots->fMuonDxy, "iso muon", "l");
307	legend->AddEntry(plots->fSoftMuonDxy, "soft muon", "l");
308
309	// attach legend to stack (legend will be printed over stack in .eps file)
310
311	result->Attach(stack, legend);
312
313
314	// book general comment
315
316	comment = result->AddComment(0.64, 0.86, 0.98, 0.98);
317	comment->AddText("demonstration plot");
318	comment->AddText("produced by Example2.C");
319
320	// attach comment to single histograms
321	/*
322	result->Attach(plots->fJetPT[0], comment);
323	result->Attach(plots->fJetPT[1], comment);
324	result->Attach(plots->fMuonPT, comment);
325	result->Attach(plots->fElectronPT, comment);
326	*/
327
328	// show histogram statistics for MissingET
329	//plots->fMissingET->SetStats();
330	}
331
332	//------------------------------------------------------------------------------
333
334	void AnalyseEvents(ExRootTreeReader treeReader, MyPlots plots, bool debug, Int_t verbose)
335	{
336	// counters initialized globally
337	Int_t counter[12] = {};
338	Int_t nMatchedSoftMuons=0;
339
340	TClonesArray *branchParticle = treeReader->UseBranch("Particle");
341	TClonesArray *branchJet = treeReader->UseBranch("Jet");
342	TClonesArray *branchMuon = treeReader->UseBranch("Muon");
343	TClonesArray *branchSoftMuon = treeReader->UseBranch("SoftMuon");
344	TClonesArray *branchElectron = treeReader->UseBranch("Electron");
345	TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");
346	TClonesArray *branchTrack = treeReader->UseBranch("Track");
347
348	Long64_t allEntries = treeReader->GetEntries();
349	if (debug) allEntries = 50;
350
351	cout << "** Chain contains " << allEntries << " events" << endl;
352
353	TObject *object;
354	Jet *jet;
355	MissingET *met;
356	Electron *electron;
357	Muon *muon;
358	Muon *muonInJet;
359	Muon *softmuon;
360	GenParticle *particle;
361	Track *track;
362
363	TLorentzVector bquarkLV;
364	TLorentzVector antibquarkLV;
365	TLorentzVector LambdabLV;
366	TLorentzVector antiLambdabLV;
367	Int_t index_bquark;
368	Int_t index_antibquark;
369	Int_t index_lambdab;
370	Int_t index_antilambdab;
371
372	Long64_t entry;
373
374	Int_t i;
375
376	bool selected[12];
377	Int_t nHardLeptons;
378	Int_t nHardMuons;
379	Int_t nHardElectrons;
380	Int_t nJets;
381	Int_t nBtags;
382	Int_t nSoftMuons;
383	Int_t nMuonsInJets;
384	Int_t nTaggedSoftMuons;
385
386	Int_t chargeHardLepton;
387	Int_t chargeSoftMuon;
388
389	TLorentzVector theHardLeptonLV;
390	TLorentzVector theSoftMuonLV;
391	vector< TLorentzVector > theTrueMuons;
392	vector< TLorentzVector > theHardJets;
393	vector< TLorentzVector > theCorrHardJets;
394	TLorentzVector theJetWithSoftMuon;
395	TLorentzVector theCorrJetWithSoftMuon;
396
397	// Loop over all events
398	for(entry = 0; entry < allEntries; ++entry)
399	{
400	if (verbose>=0) cout << "******* event n." << entry << endl;
401
402	// Load selected branches with data from specified event
403	treeReader->ReadEntry(entry);
404
405	// things to be initialized per event
406	nHardLeptons = 0;
407	nHardMuons = 0;
408	nHardElectrons = 0;
409	nJets = 0;
410	nBtags = 0;
411	nSoftMuons = 0;
412	nTaggedSoftMuons = 0;
413	nMuonsInJets = 0;
414
415	index_bquark=-1;
416	index_antibquark=-1;
417	index_lambdab=-1;
418	index_antilambdab=-1;
419
420	theTrueMuons.clear();
421	theHardJets.clear();
422	theCorrHardJets.clear();
423
424	chargeHardLepton = 0;
425	chargeSoftMuon = 0;
426
427	bquarkLV.SetPxPyPzE(0,0,0,0);
428	antibquarkLV.SetPxPyPzE(0,0,0,0);
429	LambdabLV.SetPxPyPzE(0,0,0,0);
430	antiLambdabLV.SetPxPyPzE(0,0,0,0);
431
432	theHardLeptonLV.SetPxPyPzE(0,0,0,0);
433	theSoftMuonLV.SetPxPyPzE(0,0,0,0);
434	theJetWithSoftMuon.SetPxPyPzE(0,0,0,0);
435	theCorrJetWithSoftMuon.SetPxPyPzE(0,0,0,0);
436
437	// Find real b, anti-b, lambda_b, anti-lambda_b
438	for(i = 0; i < branchParticle->GetEntriesFast(); ++i) // note: if more than one is found, the assumption (used also in Delphes internally) is that the latest is the one to be used
439	{
440	particle = (GenParticle*) branchParticle->At(i);
441	Int_t pid = particle->PID;
442	if (pid==5 && index_bquark < 0) index_bquark = i;
443	if (pid==-5 && index_antibquark < 0) index_antibquark = i;
444	if (pid==5122) index_lambdab = i;
445	if (pid==-5122) index_antilambdab = i;
446	Float_t distanceFromVertex = sqrt( (particle->X)2 + (particle->Y)2 + (particle->Z)**2 );
447	plots->fVertexAllParticles->Fill(distanceFromVertex);
448	if (abs(pid)==13 && particle->Status == 1) { // the true final-state muons
449	if (verbose>=2) cout << " Muon-MC-truth pt: " << particle->PT << ", eta: " << particle->Eta << ", phi: " << particle->Phi << endl;
450	TLorentzVector trueMuonLV;
451	trueMuonLV.SetPxPyPzE(particle->Px,particle->Py,particle->Pz,particle->E);
452	theTrueMuons.push_back(trueMuonLV);
453	}
454	}
455	if (index_bquark > -1) {
456	GenParticle bquark = (GenParticle) branchParticle->At(index_bquark);
457	bquarkLV.SetPxPyPzE(bquark->Px,bquark->Py,bquark->Pz,bquark->E);
458	} else cout << "There is no b quark in the event!" << endl;
459
460	if (index_antibquark > -1) {
461	GenParticle antibquark = (GenParticle) branchParticle->At(index_antibquark);
462	antibquarkLV.SetPxPyPzE(antibquark->Px,antibquark->Py,antibquark->Pz,antibquark->E);
463	} else cout << "There is no anti-b quark in the event!" << endl;
464
465	if (index_lambdab > -1 ) {
466	GenParticle Lambdab = (GenParticle) branchParticle->At(index_lambdab);
467	LambdabLV.SetPxPyPzE(Lambdab->Px,Lambdab->Py,Lambdab->Pz,Lambdab->E);
468	for (i = Lambdab->D1; i<Lambdab->D2; ++i) {
469	GenParticle daughter = (GenParticle) branchParticle->At(i);
470	Float_t distanceFromVertex = sqrt( (daughter->X)2 + (daughter->Y)2 + (daughter->Z)**2 );
471	if (distanceFromVertex==0) cout << "daughter; distance == 0" << endl;
472	plots->fVertexLambdabDaughters->Fill(distanceFromVertex);
473	}
474	}
475	if (index_antilambdab > -1) {
476	GenParticle antiLambdab = (GenParticle) branchParticle->At(index_antilambdab);
477	antiLambdabLV.SetPxPyPzE(antiLambdab->Px,antiLambdab->Py,antiLambdab->Pz,antiLambdab->E);
478	for (i = antiLambdab->D1; i<antiLambdab->D2; ++i) {
479	GenParticle daughter = (GenParticle) branchParticle->At(i);
480	Float_t distanceFromVertex = sqrt( (daughter->X)2 + (daughter->Y)2 + (daughter->Z)**2 );
481	if (distanceFromVertex==0) cout << "daughter; distance == 0" << endl;
482	plots->fVertexLambdabDaughters->Fill(distanceFromVertex);
483	}
484	}
485
486
487
488	// Loop over iso muons in event
489	for(i = 0; i < branchMuon->GetEntriesFast(); ++i)
490	{
491	muon = (Muon*) branchMuon->At(i);
492	plots->fMuonPT->Fill(muon->PT);
493	plots->fMuonDxy->Fill(muon->Dxy);
494	if (muon->PT > 26 && fabs(muon->Eta)<2.1) {
495	nHardMuons++;
496	theHardLeptonLV.SetPtEtaPhiM(muon->PT,muon->Eta,muon->Phi,0.);
497	chargeHardLepton = muon->Charge;
498	cout << " Hard Muon pt: " << muon->PT << ", eta: " << muon->Eta << ", phi: " << muon->Phi << endl;
499	}
500	}
501
502	// Loop over iso electrons in event
503	for(i = 0; i < branchElectron->GetEntriesFast(); ++i)
504	{
505	electron = (Electron*) branchElectron->At(i);
506	plots->fElectronPT->Fill(electron->PT);
507	if (electron->PT > 26 && fabs(electron->Eta)<2.4) {
508	nHardElectrons++;
509	theHardLeptonLV.SetPtEtaPhiM(electron->PT,electron->Eta,electron->Phi,0.);
510	chargeHardLepton = electron->Charge;
511	}
512	}
513
514	nHardLeptons = nHardMuons + nHardElectrons;
515
516	if (verbose>=1) cout << "nHardLeptons = " << nHardLeptons << " (" << nHardMuons << " muons, " << nHardElectrons << " electrons)" << endl;
517	plots->NHardLeptons->Fill((Float_t)nHardLeptons);
518	if (nHardLeptons==1){
519	counter[0]++;
520	selected[0]=true;
521	} else selected[0]=false;
522
523	// Loop over jets in event
524	for(i = 0; i < branchJet->GetEntriesFast(); ++i)
525	{
526	jet = (Jet*) branchJet->At(i);
527	if (jet->PT > 30 && fabs(jet->Eta)<2.4){
528	TLorentzVector jetLV;
529	jetLV.SetPtEtaPhiM(jet->PT,jet->Eta,jet->Phi,0.);
530
531	cout << " Jet pt: " << jet->PT << ", eta: " << jet->Eta << ", phi: " << jet->Phi << endl;
532
533	theHardJets.push_back(jetLV);
534	// Correction by looping over constituents and subtracting those compatible with primary vertex
535	Int_t nTracksInJet = 0;
536	Int_t nTracksInJetFromPV = 0;
537	TLorentzVector momentumFromPV;
538	momentumFromPV.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
539	for(Int_t j = 0; j < jet->Constituents.GetEntriesFast(); ++j) {
540	object = jet->Constituents.At(j);
541	// Check if the constituent is accessible
542	if(object == 0) continue;
543
544	if(object->IsA() == Track::Class()) {
545	track = (Track*) object;
546	nTracksInJet++;
547	//cout << " Track X: " << track->X << ", Y: " << track->Y << ", Z: " << track->Z << endl;
548	if (track->X == 0 && track->Y == 0 && track->Z == 0) { // primary vertex has exactly coordinates (0,0,0), when the PU module is not executed
549	nTracksInJetFromPV++;
550	momentumFromPV += track->P4();
551	}
552
553	if( fabs(track->PID)==13 )
554	{
555	if (verbose>=2) cout << " Muon-in-jet pt: " << track->PT << ", eta: " << track->Eta << ", phi: " << track->Phi << endl;
556	plots->fMuonInJetPT->Fill(track->PT);
557	nMuonsInJets++;
558	// matching with true muons:
559	TLorentzVector muonInJetLV;
560	muonInJetLV.SetPtEtaPhiM(track->PT,track->Eta,track->Phi,0.);
561	vector<Double_t> drV;
562	drV.clear();
563	for (Int_t k=0; k<theTrueMuons.size(); ++k) drV.push_back(muonInJetLV.DeltaR(theTrueMuons[k]));
564	if (drV.size() > 0) {
565	pair<Int_t,Double_t> pair_index_dRmin = findClosestObject(drV);
566	// if (pair_index_dRmin.second < 0.5) { // dR minimum < 0.5 (i.e. jet size)
567	}
568	}
569
570	} // end check if track
571
572	} // end of loop on constituents
573
574	if (verbose>=3 && selected[0]) cout << "number of tracks: " << nTracksInJet
575	<< " - number of tracks from PV: " << nTracksInJetFromPV
576	<< " - number of constituents: " << jet->Constituents.GetEntriesFast()
577	<< " - number of particles: " << jet->Particles.GetEntriesFast() << endl;
578	theCorrHardJets.push_back(jetLV - momentumFromPV);
579	nJets++;
580	if (jet->BTag) nBtags++;
581	}
582	} // end loop over jets ...
583
584
585	if (verbose>=1) cout << "nJets = " << nJets << endl;
586	plots->NJets->Fill((Float_t)nJets);
587	if (verbose>=1) cout << "nBtags = " << nBtags << endl;
588	plots->NBtags->Fill((Float_t)nBtags);
589	if (verbose>=1) cout << "nMuonsInJets = " << nMuonsInJets << endl;
590	plots->NMuonsInJets->Fill((Float_t)nMuonsInJets);
591	if (nJets>=4){
592	counter[1]++;
593	selected[1]=true;
594	} else selected[1]=false;
595
596	// Loop over soft muons in event
597	for(i = 0; i < branchSoftMuon->GetEntriesFast(); ++i)
598	{
599	softmuon = (Muon*) branchSoftMuon->At(i);
600	plots->fSoftMuonPT->Fill(softmuon->PT);
601	plots->fSoftMuonDxy->Fill(softmuon->Dxy);
602	nSoftMuons++;
603	if (verbose>=2) cout << " Soft Muon pt: " << softmuon->PT << ", eta: " << softmuon->Eta << ", phi: " << softmuon->Phi << endl;
604	if (fabs(softmuon->Dxy)>0.17) {///// Cut chosen such to give efficiency close to 8% in an inclusive ttbar sample
605	theSoftMuonLV.SetPtEtaPhiM(softmuon->PT,softmuon->Eta,softmuon->Phi,0.);
606	chargeSoftMuon = softmuon->Charge;
607	vector<Double_t> drV;
608	drV.clear();
609	for (Int_t j=0; j<theHardJets.size(); ++j) drV.push_back(theSoftMuonLV.DeltaR(theHardJets[j]));
610	pair<Int_t,Double_t> pair_index_dRmin = findClosestObject(drV);
611	if (verbose>=2) cout << "Index of jet closest to soft muon: " << pair_index_dRmin.first << " - DRmin = " << pair_index_dRmin.second << endl;
612	if (pair_index_dRmin.second < 0.5) {
613	nTaggedSoftMuons++;
614	theJetWithSoftMuon = theHardJets[pair_index_dRmin.first];
615	theCorrJetWithSoftMuon = theCorrHardJets[pair_index_dRmin.first];
616	}
617	}
618	}
619	if (verbose>=1) cout << "nSoftMuons = " << nSoftMuons << endl;
620	plots->NSoftMuons->Fill((Float_t)nSoftMuons);
621	if (nSoftMuons>=1){
622	counter[2]++;
623	selected[2]=true;
624	} else selected[2]=false;
625	if (verbose>=1) cout << "nTaggedSoftMuons = " << nTaggedSoftMuons << endl;
626	plots->NTaggedSoftMuons->Fill((Float_t)nTaggedSoftMuons);
627	if (nTaggedSoftMuons==1){
628	counter[3]++;
629	selected[3]=true;
630	} else selected[3]=false;
631
632	vector<Double_t> drV;
633	// Matching between lambda_b and b/anti-b
634	if (index_bquark > -1 && index_antibquark > -1) { // just for simplicity, request both (it throws away very few events)
635	drV.clear();
636	if (index_lambdab > -1) {
637	drV.push_back(LambdabLV.DeltaR(bquarkLV));
638	drV.push_back(LambdabLV.DeltaR(antibquarkLV));
639	} else if (index_antilambdab > -1) { // note: by this logic, we don't profit from events with both lambda_b and anti-lambda_b; but we don't need much statistics for this plot
640	drV.push_back(antiLambdabLV.DeltaR(bquarkLV));
641	drV.push_back(antiLambdabLV.DeltaR(antibquarkLV));
642	}
643	if (drV.size() > 0) {
644	pair<Int_t,Double_t> pair_index_dRmin = findClosestObject(drV);
645	if (index_lambdab > -1 && pair_index_dRmin.first == 0) { // lambda_b + b quark
646	plots->fPtLambdabOverB->Fill(LambdabLV.Pt()/bquarkLV.Pt());
647	plots->fPxLambdabOverB->Fill(LambdabLV.Px()/bquarkLV.Px());
648	plots->fPyLambdabOverB->Fill(LambdabLV.Py()/bquarkLV.Py());
649	plots->fPzLambdabOverB->Fill(LambdabLV.Pz()/bquarkLV.Pz());
650	plots->fELambdabOverB->Fill(LambdabLV.E()/bquarkLV.E());
651	} else if (index_lambdab > -1 && pair_index_dRmin.first == 1) { // lambda_b + anti-b quark
652	plots->fPtLambdabOverB->Fill(LambdabLV.Pt()/antibquarkLV.Pt());
653	plots->fPxLambdabOverB->Fill(LambdabLV.Px()/antibquarkLV.Px());
654	plots->fPyLambdabOverB->Fill(LambdabLV.Py()/antibquarkLV.Py());
655	plots->fPzLambdabOverB->Fill(LambdabLV.Pz()/antibquarkLV.Pz());
656	plots->fELambdabOverB->Fill(LambdabLV.E()/antibquarkLV.E());
657	} else if (index_antilambdab > -1 && pair_index_dRmin.first == 0) { // anti-lambda_b + b quark
658	plots->fPtLambdabOverB->Fill(antiLambdabLV.Pt()/bquarkLV.Pt());
659	plots->fPxLambdabOverB->Fill(antiLambdabLV.Px()/bquarkLV.Px());
660	plots->fPyLambdabOverB->Fill(antiLambdabLV.Py()/bquarkLV.Py());
661	plots->fPzLambdabOverB->Fill(antiLambdabLV.Pz()/bquarkLV.Pz());
662	plots->fELambdabOverB->Fill(antiLambdabLV.E()/bquarkLV.E());
663	} else if (index_antilambdab > -1 && pair_index_dRmin.first == 1) { // anti-lambda_b + anti-b quark
664	plots->fPtLambdabOverB->Fill(antiLambdabLV.Pt()/antibquarkLV.Pt());
665	plots->fPxLambdabOverB->Fill(antiLambdabLV.Px()/antibquarkLV.Px());
666	plots->fPyLambdabOverB->Fill(antiLambdabLV.Py()/antibquarkLV.Py());
667	plots->fPzLambdabOverB->Fill(antiLambdabLV.Pz()/antibquarkLV.Pz());
668	plots->fELambdabOverB->Fill(antiLambdabLV.E()/antibquarkLV.E());
669	}
670	}
671	}
672
673	if (selected[3]) {
674
675	// Matching between soft muon and b/anti-b (used to estimate soft-muon-tag efficiency)
676	Float_t dR_sm_b = theSoftMuonLV.DeltaR(bquarkLV);
677	Float_t dR_sm_antib = theSoftMuonLV.DeltaR(antibquarkLV);
678	if (chargeSoftMuon > 0) {
679	plots->fDRPlusMuonBQuark->Fill(dR_sm_b);
680	plots->fDRPlusMuonAntiBQuark->Fill(dR_sm_antib);
681	} else if (chargeSoftMuon < 0) {
682	plots->fDRMinusMuonBQuark->Fill(dR_sm_b);
683	plots->fDRMinusMuonAntiBQuark->Fill(dR_sm_antib);
684	}
685	drV.clear();
686	if (index_bquark > -1) drV.push_back(theSoftMuonLV.DeltaR(bquarkLV));
687	if (index_antibquark > -1) drV.push_back(theSoftMuonLV.DeltaR(antibquarkLV));
688	if (drV.size() > 0) {
689	pair<Int_t,Double_t> pair_index_dRmin = findClosestObject(drV);
690	if (pair_index_dRmin.second < 0.5) nMatchedSoftMuons++; // dR minimum < 0.5 (i.e. jet size)
691	}
692
693	// Matching between lambda_b and b/anti-b
694	if (index_bquark > -1 && index_antibquark > -1) { // just for simplicity, request both (it throws away very few events)
695	drV.clear();
696	if (index_lambdab > -1) {
697	drV.push_back(LambdabLV.DeltaR(bquarkLV));
698	drV.push_back(LambdabLV.DeltaR(antibquarkLV));
699	} else if (index_antilambdab > -1) { // note: by this logic, we don't profit from events with both lambda_b and anti-lambda_b; but we don't need much statistics for this plot
700	drV.push_back(antiLambdabLV.DeltaR(bquarkLV));
701	drV.push_back(antiLambdabLV.DeltaR(antibquarkLV));
702	}
703	if (drV.size() > 0) {
704	pair<Int_t,Double_t> pair_index_dRmin = findClosestObject(drV);
705	if (index_lambdab > -1 && pair_index_dRmin.first == 0) { // lambda_b + b quark
706	plots->fPtLambdabOverB_selected->Fill(LambdabLV.Pt()/bquarkLV.Pt());
707	plots->fPxLambdabOverB_selected->Fill(LambdabLV.Px()/bquarkLV.Px());
708	plots->fPyLambdabOverB_selected->Fill(LambdabLV.Py()/bquarkLV.Py());
709	plots->fPzLambdabOverB_selected->Fill(LambdabLV.Pz()/bquarkLV.Pz());
710	plots->fELambdabOverB_selected->Fill(LambdabLV.E()/bquarkLV.E());
711	} else if (index_lambdab > -1 && pair_index_dRmin.first == 1) { // lambda_b + anti-b quark
712	plots->fPtLambdabOverB_selected->Fill(LambdabLV.Pt()/antibquarkLV.Pt());
713	plots->fPxLambdabOverB_selected->Fill(LambdabLV.Px()/antibquarkLV.Px());
714	plots->fPyLambdabOverB_selected->Fill(LambdabLV.Py()/antibquarkLV.Py());
715	plots->fPzLambdabOverB_selected->Fill(LambdabLV.Pz()/antibquarkLV.Pz());
716	plots->fELambdabOverB_selected->Fill(LambdabLV.E()/antibquarkLV.E());
717	} else if (index_antilambdab > -1 && pair_index_dRmin.first == 0) { // anti-lambda_b + b quark
718	plots->fPtLambdabOverB_selected->Fill(antiLambdabLV.Pt()/bquarkLV.Pt());
719	plots->fPxLambdabOverB_selected->Fill(antiLambdabLV.Px()/bquarkLV.Px());
720	plots->fPyLambdabOverB_selected->Fill(antiLambdabLV.Py()/bquarkLV.Py());
721	plots->fPzLambdabOverB_selected->Fill(antiLambdabLV.Pz()/bquarkLV.Pz());
722	plots->fELambdabOverB_selected->Fill(antiLambdabLV.E()/bquarkLV.E());
723	} else if (index_antilambdab > -1 && pair_index_dRmin.first == 1) { // anti-lambda_b + anti-b quark
724	plots->fPtLambdabOverB_selected->Fill(antiLambdabLV.Pt()/antibquarkLV.Pt());
725	plots->fPxLambdabOverB_selected->Fill(antiLambdabLV.Px()/antibquarkLV.Px());
726	plots->fPyLambdabOverB_selected->Fill(antiLambdabLV.Py()/antibquarkLV.Py());
727	plots->fPzLambdabOverB_selected->Fill(antiLambdabLV.Pz()/antibquarkLV.Pz());
728	plots->fELambdabOverB_selected->Fill(antiLambdabLV.E()/antibquarkLV.E());
729	}
730	}
731	}
732
733	// Matching between the jet with soft mu and b/anti-b
734	drV.clear();
735	if (index_bquark > -1) drV.push_back(theJetWithSoftMuon.DeltaR(bquarkLV));
736	if (index_antibquark > -1) drV.push_back(theJetWithSoftMuon.DeltaR(antibquarkLV));
737	if (drV.size() > 0) {
738	pair<Int_t,Double_t> pair_index_dRmin = findClosestObject(drV);
739	if (pair_index_dRmin.second < 0.5) { // dR minimum < 0.5 (i.e. jet size)
740	if (pair_index_dRmin.first == 0) { // b
741	plots->fPtJetOverB->Fill(theJetWithSoftMuon.Pt()/bquarkLV.Pt());
742	plots->fPxJetOverB->Fill(theJetWithSoftMuon.Px()/bquarkLV.Px());
743	plots->fPyJetOverB->Fill(theJetWithSoftMuon.Py()/bquarkLV.Py());
744	plots->fPzJetOverB->Fill(theJetWithSoftMuon.Pz()/bquarkLV.Pz());
745	plots->fEJetOverB->Fill(theJetWithSoftMuon.E()/bquarkLV.E());
746	} else if (pair_index_dRmin.first == 0) { // anti-b
747	plots->fPtJetOverB->Fill(theJetWithSoftMuon.Pt()/antibquarkLV.Pt());
748	plots->fPxJetOverB->Fill(theJetWithSoftMuon.Px()/antibquarkLV.Px());
749	plots->fPyJetOverB->Fill(theJetWithSoftMuon.Py()/antibquarkLV.Py());
750	plots->fPzJetOverB->Fill(theJetWithSoftMuon.Pz()/antibquarkLV.Pz());
751	plots->fEJetOverB->Fill(theJetWithSoftMuon.E()/antibquarkLV.E());
752	}
753	}
754	}
755
756	// Matching between the jet with soft mu and lambda_b/anti-lambda_b
757	drV.clear();
758	if (index_lambdab > -1) drV.push_back(theJetWithSoftMuon.DeltaR(LambdabLV));
759	if (index_antilambdab > -1) drV.push_back(theJetWithSoftMuon.DeltaR(antiLambdabLV));
760	if (drV.size() > 0) {
761	pair<Int_t,Double_t> pair_index_dRmin = findClosestObject(drV);
762	if (pair_index_dRmin.second < 0.5) { // dR minimum < 0.5 (i.e. jet size)
763	if (pair_index_dRmin.first == 0) { // b
764	plots->fPtJetOverLambdab->Fill(theJetWithSoftMuon.Pt()/LambdabLV.Pt());
765	plots->fPxJetOverLambdab->Fill(theJetWithSoftMuon.Px()/LambdabLV.Px());
766	plots->fPyJetOverLambdab->Fill(theJetWithSoftMuon.Py()/LambdabLV.Py());
767	plots->fPzJetOverLambdab->Fill(theJetWithSoftMuon.Pz()/LambdabLV.Pz());
768	plots->fEJetOverLambdab->Fill(theJetWithSoftMuon.E()/LambdabLV.E());
769	plots->fDeltaEtaJetLambdab->Fill(theJetWithSoftMuon.Eta()-LambdabLV.Eta());
770	plots->fDeltaPhiJetLambdab->Fill(theJetWithSoftMuon.Phi()-LambdabLV.Phi());
771
772	plots->fPtCorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.Pt()/LambdabLV.Pt());
773	plots->fPxCorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.Px()/LambdabLV.Px());
774	plots->fPyCorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.Py()/LambdabLV.Py());
775	plots->fPzCorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.Pz()/LambdabLV.Pz());
776	plots->fECorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.E()/LambdabLV.E());
777	plots->fDeltaEtaCorrJetLambdab->Fill(theCorrJetWithSoftMuon.Eta()-LambdabLV.Eta());
778	plots->fDeltaPhiCorrJetLambdab->Fill(theCorrJetWithSoftMuon.Phi()-LambdabLV.Phi());
779	} else if (pair_index_dRmin.first == 1) { // anti-b
780	plots->fPtJetOverLambdab->Fill(theJetWithSoftMuon.Pt()/antiLambdabLV.Pt());
781	plots->fPxJetOverLambdab->Fill(theJetWithSoftMuon.Px()/antiLambdabLV.Px());
782	plots->fPyJetOverLambdab->Fill(theJetWithSoftMuon.Py()/antiLambdabLV.Py());
783	plots->fPzJetOverLambdab->Fill(theJetWithSoftMuon.Pz()/antiLambdabLV.Pz());
784	plots->fEJetOverLambdab->Fill(theJetWithSoftMuon.E()/antiLambdabLV.E());
785	plots->fDeltaEtaJetLambdab->Fill(theJetWithSoftMuon.Eta()-antiLambdabLV.Eta());
786	plots->fDeltaPhiJetLambdab->Fill(theJetWithSoftMuon.Phi()-antiLambdabLV.Phi());
787
788	plots->fPtCorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.Pt()/antiLambdabLV.Pt());
789	plots->fPxCorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.Px()/antiLambdabLV.Px());
790	plots->fPyCorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.Py()/antiLambdabLV.Py());
791	plots->fPzCorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.Pz()/antiLambdabLV.Pz());
792	plots->fECorrJetOverLambdab->Fill(theCorrJetWithSoftMuon.E()/antiLambdabLV.E());
793	plots->fDeltaEtaCorrJetLambdab->Fill(theCorrJetWithSoftMuon.Eta()-antiLambdabLV.Eta());
794	plots->fDeltaPhiCorrJetLambdab->Fill(theCorrJetWithSoftMuon.Phi()-antiLambdabLV.Phi());
795	}
796	}
797	}
798
799	} // end of IF selected[3]
800
801
802	// Analyse missing ET
803	if(branchMissingET->GetEntriesFast() > 0)
804	{
805	met = (MissingET*) branchMissingET->At(0);
806	plots->fMissingET->Fill(met->MET);
807	}
808
809	} // end of loop over events
810	for (Int_t j=0; j<4; ++j) plots->cutFlow->Fill(j,(Float_t)counter[j]/(Float_t)allEntries);
811
812
813	cout << "Number of soft muons matched to b quarks: " << nMatchedSoftMuons << ", hence efficiency of soft-muon-tagging is " << (Float_t)nMatchedSoftMuons/(Float_t)allEntries << endl;
814	}
815
816	//------------------------------------------------------------------------------
817
818	void PrintHistograms(ExRootResult result, MyPlots plots)
819	{
820	result->Print("png");
821	}
822
823	//------------------------------------------------------------------------------
824
825	void analyzeLambdab(const char inputFile, const char outputString, bool debug, Int_t verbose = 1)
826	{
827	gSystem->Load("libDelphes");
828
829	TChain *chain = new TChain("Delphes");
830	chain->Add(inputFile);
831
832	ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
833	ExRootResult *result = new ExRootResult();
834
835	MyPlots *plots = new MyPlots;
836
837	BookHistograms(result, plots);
838
839	AnalyseEvents(treeReader, plots, debug, verbose);
840
841	PrintHistograms(result, plots);
842
843	TString outputRootFile("results_");
844	outputRootFile+=outputString;
845	outputRootFile+=".root";
846	result->Write(outputRootFile);
847
848	cout << "** Exiting..." << endl;
849
850	delete plots;
851	delete result;
852	delete treeReader;
853	delete chain;
854	}
855
856	//------------------------------------------------------------------------------
857

Download in other formats:

Original Format