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 |
|
---|
16 | using namespace std;
|
---|
17 |
|
---|
18 | //_________________________________________________________________________
|
---|
19 |
|
---|
20 | void 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 |
|
---|