Fork me on GitHub

Ticket #1061: DelphesFastAnalysis.C

File DelphesFastAnalysis.C, 18.0 KB (added by Soeren Stamm, 8 years ago)

Script for creating mtW plot

Line 
1#include <TROOT.h>
2#include <TSystem.h>
3#include <TFile.h>
4#include <TH1F.h>
5#include <TChain.h>
6#include <TStopwatch.h>
7#include <TStyle.h>
8#include <TClonesArray.h>
9#include <TInterpreter.h>
10#include <classes/DelphesClasses.h>
11#include <iostream>
12#include <TBranchElement.h>
13#include <TClonesArray.h>
14
15
16using namespace std;
17
18//_________________________________________________________________________
19
20void DelphesFastAnalysis(const char* InputFiles, const char* OutFile,
21 Bool_t WriteTree = kFALSE,
22 Bool_t tchannel_selection = kFALSE,
23 Long64_t NEvt = TChain::kBigNumber) {
24
25 //
26 // Delphes Analysis
27 //
28 ProcInfo_t *ProcInfo = new ProcInfo_t;
29 MemInfo_t *MemInfo = new MemInfo_t;
30 TStopwatch StopWatch;
31 StopWatch.Start();
32 cout << "DelphesAnalysis :\tStart analysis" << endl;
33
34 // ===================
35 // Initialize analysis
36 // ===================
37
38 // Set style
39 gStyle->SetOptStat(11);
40
41 // Create chain and add input files
42 // (wildcards in pathanme are allowed)
43 TChain *ch = new TChain("Delphes");
44
45 FILE *flist = gSystem->OpenPipe(Form("ls -1 %s", InputFiles), "r");
46 TString line;
47 TObjArray *FileList = new TObjArray;
48 while( !feof(flist) ) {
49 line.Gets(flist);
50 if ( line.Length() > 0 ) {
51 FileList->Add(new TObjString(line.Data()));
52 }
53 }
54 gSystem->ClosePipe(flist);
55 for ( Int_t i = 0; i < FileList->GetEntries(); i++ ) {
56 ch->Add(((TObjString*)FileList->At(i))->GetString().Data());
57 }
58 cout << "DelphesAnalysis :\tAdded " << FileList->GetEntries()
59 << " input files." << endl;
60
61 // Setup Chain without ExRootTreeReader
62
63 // Setup Branch Status
64 ch->SetBranchStatus("*", 0);
65 ch->SetBranchStatus("Event*", 0);
66 ch->SetBranchStatus("EventLHEF*", 0);
67 ch->SetBranchStatus("WeightLHEF*", 0);
68 ch->SetBranchStatus("Particle*", 0);
69 ch->SetBranchStatus("Track*", 0);
70 ch->SetBranchStatus("Tower*", 0);
71 ch->SetBranchStatus("EFlowTrack*", 0);
72 ch->SetBranchStatus("EFlowPhoton*", 0);
73 ch->SetBranchStatus("EFlowNeutralHadron*", 0);
74 ch->SetBranchStatus("GenJet*", 0);
75 ch->SetBranchStatus("GenMissingET*", 0);
76 ch->SetBranchStatus("Jet*", 1);
77 ch->SetBranchStatus("Electron*", 1);
78 ch->SetBranchStatus("Photon*", 0);
79 ch->SetBranchStatus("Muon*", 1);
80 ch->SetBranchStatus("MissingET*", 1);
81 ch->SetBranchStatus("ScalarHT*", 0);
82
83 ch->SetBranchStatus("EFlowNeutralHadron_size", 0);
84 ch->SetBranchStatus("EFlowPhoton_size", 0);
85 ch->SetBranchStatus("EFlowTrack_size", 0);
86 ch->SetBranchStatus("Electron_size", 1);
87 ch->SetBranchStatus("EventLHEF_size", 0);
88 ch->SetBranchStatus("Event_size", 0);
89 ch->SetBranchStatus("GenJet_size", 0);
90 ch->SetBranchStatus("GenMissingET_size", 0);
91 ch->SetBranchStatus("Jet_size", 1);
92 ch->SetBranchStatus("MissingET_size", 1);
93 ch->SetBranchStatus("Muon_size", 1);
94 ch->SetBranchStatus("Particle_size", 0);
95 ch->SetBranchStatus("Photon_size", 0);
96 ch->SetBranchStatus("ScalarHT_size", 0);
97 ch->SetBranchStatus("Tower_size", 0);
98 ch->SetBranchStatus("Track_size", 0);
99 ch->SetBranchStatus("WeightLHEF_size", 0);
100
101
102 // Setup Branch Address
103 TBranchElement *branchElectron = (TBranchElement*) ch->GetBranch("Electron");
104 TClonesArray* arrayElectron = new TClonesArray("Electron", branchElectron->GetMaximum());
105 ch->SetBranchAddress("Electron", &arrayElectron);
106
107 TBranchElement *branchMuon = (TBranchElement*) ch->GetBranch("Muon");
108 TClonesArray* arrayMuon = new TClonesArray("Muon", branchMuon->GetMaximum());
109 ch->SetBranchAddress("Muon", &arrayMuon);
110
111 TBranchElement *branchJet = (TBranchElement*) ch->GetBranch("Jet");
112 TClonesArray* arrayJet = new TClonesArray("Jet", branchJet->GetMaximum());
113 ch->SetBranchAddress("Jet", &arrayJet);
114
115 TBranchElement *branchMissingET = (TBranchElement*) ch->GetBranch("MissingET");
116 TClonesArray* arrayMissingET = new TClonesArray("MissingET", branchMissingET->GetMaximum());
117 ch->SetBranchAddress("MissingET", &arrayMissingET);
118
119 TBranchElement *branchParticle = (TBranchElement*) ch->GetBranch("Particle");
120 TClonesArray* arrayParticle = new TClonesArray("GenParticle", branchParticle->GetMaximum());
121 ch->SetBranchAddress("Particle", &arrayParticle);
122
123 // TClonesArray *arrayElectron = new TClonesArray( gROOT->GetClass("Electron"), );
124 // const char *classJet = element->GetClonesName();
125 // Int_t size = element->GetMaximum();
126 // TClass *cl = gROOT->GetClass(className);
127 // if(cl)
128 // {
129 // array = new TClonesArray(cl, size);
130 // array->SetName(branchName);
131 // fBranchMap.insert(make_pair(branchName, make_pair(branch, array)));
132 // branch->SetAddress(&array);
133 // }
134
135
136 // Open output file
137 cout << "DelphesAnalysis :\tOpen output file \"" << OutFile << "\"." << endl;
138 TFile *f_out = new TFile(OutFile, "recreate");
139
140 // =================
141 // Book histograms
142 // =================
143 cout << "DelphesAnalysis :\tBook histograms." << endl;
144
145 // Cut-Flow-Hist
146 TH1F *h_cutflow = new TH1F("h_cutflow",
147 "Cut Flow Histogramm", 8, -0.5, 7.5);
148 h_cutflow->SetYTitle("Number of Entries");
149 h_cutflow->GetXaxis()->SetBinLabel(1, "N_{gen.}");
150 h_cutflow->GetXaxis()->SetBinLabel(2, "N_{lep.} = 1");
151 h_cutflow->GetXaxis()->SetBinLabel(3, "p_{T}(l) > 25 GeV");
152 h_cutflow->GetXaxis()->SetBinLabel(4, "OR");
153 h_cutflow->GetXaxis()->SetBinLabel(5, "N_{jets} = 2");
154 if ( tchannel_selection )
155 h_cutflow->GetXaxis()->SetBinLabel(6, "N_{b-jets} = 1");
156 else
157 h_cutflow->GetXaxis()->SetBinLabel(6, "N_{b-jets} = 2");
158 h_cutflow->GetXaxis()->SetBinLabel(7, "#slash{E}_{T} > 30 GeV");
159 h_cutflow->GetXaxis()->SetBinLabel(8, "m_{T}(W) > 30 GeV");
160
161
162 // Histograms
163 f_out->mkdir("Hists");
164 f_out->cd("Hists");
165
166 // Lepton Histograms
167 // -----------------
168 TH1F *fHist_Lepton_pt = new TH1F("h_Lepton_pt", "Lepton p_{T}",
169 30, 0., 150.);
170 fHist_Lepton_pt->SetXTitle("Lepton p_{T} [GeV]");
171 fHist_Lepton_pt->SetYTitle("Number of Entries");
172
173 TH1F *fHist_Lepton_eta = new TH1F("h_Lepton_eta", "Lepton #eta, tag",
174 20, -4., 4.);
175 fHist_Lepton_eta->SetXTitle("Lepton #eta");
176 fHist_Lepton_eta->SetYTitle("Number of Entries");
177
178 TH1F *fHist_Lepton_phi = new TH1F("h_Lepton_phi", "Lepton #phi, tag",
179 20, -TMath::Pi(), TMath::Pi());
180 fHist_Lepton_phi->SetXTitle("Lepton #phi");
181 fHist_Lepton_phi->SetYTitle("Number of Entries");
182
183 // Jet Histograms
184 // --------------
185 TH1F *fHist_LeadingJet_pt = new TH1F("h_LeadingJet_pt", "Leading Jet p_{T}",
186 44, 0., 220.);
187 fHist_LeadingJet_pt->SetXTitle("p_{T} [GeV]");
188 fHist_LeadingJet_pt->SetYTitle("Number of Entries");
189
190 TH1F *fHist_LeadingJet_eta = new TH1F("h_LeadingJet_eta", "Leading Jet #eta",
191 25, -5., 5.);
192 fHist_LeadingJet_eta->SetXTitle("#eta");
193 fHist_LeadingJet_eta->SetYTitle("Number of Entries");
194
195 TH1F *fHist_LeadingJet_phi = new TH1F("h_LeadingJet_phi", "Leading Jet #phi, tag",
196 20, -TMath::Pi(), TMath::Pi());
197 fHist_LeadingJet_phi->SetXTitle("#phi");
198 fHist_LeadingJet_phi->SetYTitle("Number of Entries");
199
200 // 2nd jet
201 TH1F *fHist_2ndLeadingJet_pt = new TH1F("h_2ndLeadingJet_pt",
202 "2nd Leading Jet p_{T}",
203 44, 0., 220.);
204 fHist_2ndLeadingJet_pt->SetXTitle("p_{T} [GeV]");
205 fHist_2ndLeadingJet_pt->SetYTitle("Number of Entries");
206
207 TH1F *fHist_2ndLeadingJet_eta = new TH1F("h_2ndLeadingJet_eta",
208 "2nd Leading Jet #eta",
209 25, -5., 5.);
210 fHist_2ndLeadingJet_eta->SetXTitle("#eta");
211 fHist_2ndLeadingJet_eta->SetYTitle("Number of Entries");
212
213 TH1F *fHist_2ndLeadingJet_phi = new TH1F("h_2ndLeadingJet_phi",
214 "2nd Leading Jet #phi",
215 20, -TMath::Pi(), TMath::Pi());
216 fHist_2ndLeadingJet_phi->SetXTitle("#phi");
217 fHist_2ndLeadingJet_phi->SetYTitle("Number of Entries");
218
219 // Missing Et
220 TH1F *fHist_met = new TH1F("h_met", "Missing E_{T}, tag",
221 40, 0., 200.);
222 fHist_met->SetXTitle("E_{T, miss} [GeV]");
223 fHist_met->SetYTitle("Number of Entries");
224
225 // Mt_W
226 TH1F *fHist_MtW = new TH1F("h_MtW", "Transverse W-Boson mass, tag",
227 42, 0., 210);
228 fHist_MtW->SetXTitle("M_{T, W} [GeV]");
229 fHist_MtW->SetYTitle("Number of Entries");
230
231 f_out->cd("");
232
233
234 // ==================
235 // Book Output Tree
236 // ==================
237
238 // Define Tree Variables
239 static const Int_t fgMaxElectrons = 10;
240 static const Int_t fgMaxMuons = 10;
241 static const Int_t fgMaxJets = 30;
242
243 // Event Info
244 Int_t fRunNr; // Run number
245 Int_t fEventNr; // Event number
246 Float_t fEventWeight; // Total event weight
247
248 // MET
249 Float_t fMET_Et; // MET Et (GeV)
250 Float_t fMET_Phi; // MET Phi (rad)
251
252 // Electrons
253 Int_t fElectronN; // No. of electrons
254 Float_t fElectronE[fgMaxElectrons]; // Electron E (GeV)
255 Float_t fElectronPt[fgMaxElectrons]; // Electron Pt (GeV)
256 Float_t fElectronEta[fgMaxElectrons]; // Electron Eta
257 Float_t fElectronPhi[fgMaxElectrons]; // Electron Phi (rad)
258 Int_t fElectronChg[fgMaxElectrons]; // Electron charge (+1e/-1e)
259
260 // Muons
261 Int_t fMuonN; // No. of muons
262 Float_t fMuonE[fgMaxMuons]; // Muon E (GeV)
263 Float_t fMuonPt[fgMaxMuons]; // Muon Pt (GeV)
264 Float_t fMuonEta[fgMaxMuons]; // Muon Eta
265 Float_t fMuonPhi[fgMaxMuons]; // Muon Phi (rad)
266 Int_t fMuonChg[fgMaxMuons]; // Muon charge (+1e/-1e)
267
268 // Jets
269 Int_t fJetN; // No. of jets
270 Float_t fJetE[fgMaxJets]; // Jet E (GeV)
271 Float_t fJetPt[fgMaxJets]; // Jet Pt (GeV)
272 Float_t fJetEta[fgMaxJets]; // Jet Eta
273 Float_t fJetPhi[fgMaxJets]; // Jet Phi (rad)
274 Float_t fJetBTagWeight[fgMaxJets]; // B-tag weight (defined by SetBTagger())
275 Bool_t fJetBTagged[fgMaxJets]; // Jet is b-tagged ?
276 Char_t fJetFlav[fgMaxJets]; // Jet truth flavour: l=light lfavour u,d,s; c=charm; b=beauty; t=tau; u=undefined (for DATA)
277
278 TTree *fTree = 0;
279 if ( WriteTree ) {
280 cout << "DelphesAnalysis :\tBook OutputTree." << endl;
281
282 // Book the ntuple.
283 fTree = new TTree("DelphesTree", "Delphes Tree");
284
285 // Evt header
286 fTree->Branch("run_nr", &fRunNr, "run_nr/I");
287 fTree->Branch("evt_nr", &fEventNr, "evt_nr/I");
288 fTree->Branch("evt_weight", &fEventWeight, "evt_weight/F");
289
290 // Missing Et
291 fTree->Branch("met_et", &fMET_Et, "met_et/F");
292 fTree->Branch("met_phi", &fMET_Phi, "met_phi/F");
293
294 // Electrons
295 fTree->Branch("el_n", &fElectronN, "el_n/I");
296 fTree->Branch("el_e", fElectronE, "el_e[el_n]/F");
297 fTree->Branch("el_pt", fElectronPt, "el_pt[el_n]/F");
298 fTree->Branch("el_eta", fElectronEta, "el_eta[el_n]/F");
299 fTree->Branch("el_phi", fElectronPhi, "el_phi[el_n]/F");
300 fTree->Branch("el_chg", fElectronChg, "el_chg[el_n]/I");
301
302 // Muons
303 fTree->Branch("mu_n", &fMuonN, "mu_n/I");
304 fTree->Branch("mu_e", fMuonE, "mu_e[mu_n]/F");
305 fTree->Branch("mu_pt", fMuonPt, "mu_pt[mu_n]/F");
306 fTree->Branch("mu_eta", fMuonEta, "mu_eta[mu_n]/F");
307 fTree->Branch("mu_phi", fMuonPhi, "mu_phi[mu_n]/F");
308 fTree->Branch("mu_chg", fMuonChg, "mu_chg[mu_n]/I");
309
310 // Jets
311 fTree->Branch("jet_n", &fJetN, "jet_n/I");
312 fTree->Branch("jet_e", fJetE, "jet_e[jet_n]/F");
313 fTree->Branch("jet_pt", fJetPt, "jet_pt[jet_n]/F");
314 fTree->Branch("jet_eta", fJetEta, "jet_eta[jet_n]/F");
315 fTree->Branch("jet_phi", fJetPhi, "jet_phi[jet_n]/F");
316 fTree->Branch("jet_btagged", fJetBTagged, "jet_btagged[jet_n]/O");
317 fTree->Branch("jet_btagw", fJetBTagWeight, "jet_btagw[jet_n]/F");
318 fTree->Branch("jet_flav", fJetFlav, "jet_flav[jet_n]/B");
319 }
320
321
322 // ==========
323 // Event loop
324 // ==========
325 Long64_t nevt = TMath::Min(NEvt, ch->GetEntries());
326 for ( Int_t ievt = 0; ievt < nevt; ievt++ ) {
327
328 // Print info
329 Int_t mod = 0;
330 if ( ievt < 10000 ) {
331 mod = 1000;
332 } else {
333 mod = 10000;
334 }
335 if ( (ievt+1) % mod == 0 ) {
336 printf("DelphesAnalysis :\tProcessing chain entry %-d / %-d\n",
337 ievt+1, (Int_t)nevt);
338 gSystem->GetProcInfo(ProcInfo);
339 gSystem->GetMemInfo(MemInfo);
340 printf("DelphesAnalysis :\tMemory usage proc_res=%d proc_virt=%d total=%d Mb\n",
341 Int_t(ProcInfo->fMemResident/1024),
342 Int_t(ProcInfo->fMemVirtual/1024),
343 MemInfo->fMemTotal);
344 cout.flush();
345 }
346
347 // Get event
348 ch->GetEntry(ievt);
349
350 // Fill Event Header
351 fRunNr = 1;
352 fEventNr = ievt;
353 fEventWeight = 1.;
354
355 // Count all events
356 Int_t icut = 0;
357 h_cutflow->Fill(icut);
358 icut++;
359
360 // Number of Electrons, Muons, Jets
361 Int_t Nelectrons = arrayElectron->GetEntries();
362 Int_t Nmuons = arrayMuon->GetEntries();
363 Int_t Njets = arrayJet->GetEntries();
364
365 // cout << "Nelectrons " << Nelectrons << endl;
366 // cout << "Nmuons " << Nmuons << endl;
367
368 // ----------------
369 // Cut: N leptons
370 // ----------------
371 Int_t Nelmu = Nelectrons + Nmuons;
372 if ( Nelmu != 1 )
373 continue;
374 h_cutflow->Fill(icut);
375 icut++;
376
377 // ------------------------
378 // Cut: lepton pt > 25GeV
379 // ------------------------
380 Float_t LeptonPt;
381 Float_t LeptonEta;
382 Float_t LeptonPhi;
383 Int_t LeptonCharge;
384
385 // Reset
386 fElectronN = 0;
387 fMuonN = 0;
388
389 if ( Nelectrons > 0 ) {
390 Electron *el = (Electron*) arrayElectron->At(0);
391 // cout << "el pointer " << el << endl;
392 LeptonEta = el->Eta;
393 LeptonPt = el->PT;
394 LeptonPhi = el->Phi;
395 LeptonCharge = el->Charge;
396
397 // cout << el->PT << endl;
398 // return;
399
400 // Fill electron variables
401 fElectronE[fElectronN] = el->P4().E();
402 fElectronPt[fElectronN] = LeptonPt;
403 fElectronEta[fElectronN] = LeptonEta;
404 fElectronPhi[fElectronN] = LeptonPhi;
405 fElectronChg[fElectronN] = (LeptonCharge > 0.) ? 1 : -1;
406 fElectronN++;
407 } else {
408 Muon *mu = (Muon*) arrayMuon->At(0);
409 LeptonEta = mu->Eta;
410 LeptonPt = mu->PT;
411 LeptonPhi = mu->Phi;
412 LeptonCharge = mu->Charge;
413
414 // Fill muon variables
415 fMuonE[fMuonN] = mu->P4().E();
416 fMuonPt[fMuonN] = LeptonPt;
417 fMuonEta[fMuonN] = LeptonEta;
418 fMuonPhi[fMuonN] = LeptonPhi;
419 fMuonChg[fMuonN] = (LeptonCharge > 0.) ? 1 : -1;
420 fMuonN++;
421 }
422
423 if ( LeptonPt < 25. || TMath::Abs(LeptonEta) > 2.5 )
424 continue;
425 h_cutflow->Fill(icut);
426 icut++;
427
428 // Lepton vector
429 TLorentzVector Plep;
430 Plep.SetPtEtaPhiM(LeptonPt, LeptonEta, LeptonPhi, 0.);
431
432 // --------------
433 // Cut: Overlap
434 // --------------
435
436 Bool_t overlap = false;
437
438 // Select jets (pt > 25, |eta| < 2.5)
439 Int_t Njets_sel = 0;
440 Jet *jet1 = 0;
441 Jet *jet2 = 0;
442
443 for ( Int_t ijet = 0; ijet < Njets; ++ijet ) {
444 Jet *jet = (Jet*) arrayJet->At(ijet);
445
446 Float_t JetPt = jet->PT;
447 Float_t JetEta = jet->Eta;
448 Float_t JetPhi = jet->Phi;
449
450 if ( JetPt > 25. && TMath::Abs(JetEta) < 2.5 ) {
451 Njets_sel++;
452 if ( jet1 == 0 )
453 jet1 = jet;
454 else if ( jet2 == 0 )
455 jet2 = jet;
456
457 TLorentzVector Pjet;
458 Pjet.SetPtEtaPhiM(JetPt, JetEta, JetPhi, 0.);
459
460 if ( Pjet.DeltaR(Plep) < 0.4 )
461 overlap = true;
462 }
463 }
464 if ( overlap )
465 continue;
466 h_cutflow->Fill(icut);
467 icut++;
468
469 // -------------
470 // Cut: N jets
471 // -------------
472 if ( Njets_sel != 2 )
473 continue;
474 h_cutflow->Fill(icut);
475 icut++;
476
477 // Check if jets are pt ordered
478 if ( jet1->PT < jet2->PT ) {
479 // Apply Pt order
480 Jet *tmp1 = jet1;
481 jet1 = jet2;
482 jet2 = tmp1;
483 }
484
485
486 // ---------------
487 // Cut: N b-jets
488 // ---------------
489 Int_t Nbtags = 0;
490 if ( jet1->BTag == 1 )
491 Nbtags++;
492 if ( jet2->BTag == 1 )
493 Nbtags++;
494
495 // if ( tchannel_selection ) {
496 // if ( Nbtags != 1 )
497 // continue;
498 // } else {
499 // if ( Nbtags != 2 )
500 // continue;
501 // }
502 h_cutflow->Fill(icut);
503 icut++;
504
505
506 // Fill jet variables
507 fJetN = 0;
508
509 fJetE[fJetN] = jet1->P4().E();
510 fJetPt[fJetN] = jet1->PT;
511 fJetEta[fJetN] = jet1->Eta;
512 fJetPhi[fJetN] = jet1->Phi;
513 fJetBTagWeight[fJetN] = (jet1->BTag == 1) ? 1. : 0;
514 fJetBTagged[fJetN] = (jet1->BTag == 1) ? 1 : 0; // UInt > Int conversion?
515 fJetFlav[fJetN] = (jet1->BTag == 1) ? 'b' : 'l';
516 fJetN++;
517
518 fJetE[fJetN] = jet2->P4().E();
519 fJetPt[fJetN] = jet2->PT;
520 fJetEta[fJetN] = jet2->Eta;
521 fJetPhi[fJetN] = jet2->Phi;
522 fJetBTagWeight[fJetN] = (jet2->BTag == 1) ? 1. : 0;
523 fJetBTagged[fJetN] = (jet2->BTag == 1) ? 1 : 0; // UInt > Int conversion?
524 fJetFlav[fJetN] = (jet2->BTag == 1) ? 'b' : 'l';
525 fJetN++;
526
527 // ----------
528 // Cut: MET
529 // ----------
530
531 // Fill missing et variables
532 MissingET *MET_vect = (MissingET*)arrayMissingET->At(0);
533 fMET_Et = MET_vect->MET;
534 fMET_Phi = MET_vect->Phi;
535
536 TVector2 Pmet;
537 Pmet.SetMagPhi(fMET_Et, fMET_Phi);
538
539 // if ( fMET_Et < 30. )
540 // continue;
541 h_cutflow->Fill(icut);
542 icut++;
543
544 // ----------
545 // Cut: MtW
546 // ----------
547
548 Float_t DeltaPhi = Plep.Vect().XYvector().DeltaPhi(Pmet);
549 Float_t MtW = TMath::Sqrt(2.*Plep.Et()*Pmet.Mod()*(1. - TMath::Cos(DeltaPhi)));
550
551 // don't apply MtW cut
552 if ( MtW > 30. ) {
553 h_cutflow->Fill(icut);
554 icut++;
555 }
556
557
558 // -----------------
559 // Fill histograms
560 // -----------------
561
562 fHist_Lepton_pt->Fill(LeptonPt);
563 fHist_Lepton_eta->Fill(LeptonEta);
564 fHist_Lepton_phi->Fill(LeptonPhi);
565
566 fHist_LeadingJet_pt->Fill(jet1->PT);
567 fHist_LeadingJet_eta->Fill(jet1->Eta);
568 fHist_LeadingJet_phi->Fill(jet1->Phi);
569
570 fHist_2ndLeadingJet_pt->Fill(jet2->PT);
571 fHist_2ndLeadingJet_eta->Fill(jet2->Eta);
572 fHist_2ndLeadingJet_phi->Fill(jet2->Phi);
573
574 fHist_met->Fill(fMET_Et);
575 fHist_MtW->Fill(MtW);
576
577 if ( WriteTree ) {
578 fTree->Fill();
579 }
580 }
581
582 // ====================
583 // Analysis termination
584 // ====================
585
586 // Write and close output file
587 cout << "DelphesAnalysis :\tClose output file." << endl;
588 f_out->Write();
589 delete f_out;
590
591 // Print summary
592 cout << "DelphesAnalysis :\tAnalysis complete." << endl;
593 StopWatch.Stop();
594 cout << "DelphesAnalysis :\t";
595 StopWatch.Print();
596 cout << endl;
597}
598