Fork me on GitHub

WorkBook/Tutorials/Mc4Bsm: macro.C

File macro.C, 4.3 KB (added by Michele Selvaggi, 8 years ago)
Line 
1/*
2plots mass and tau21 from two Delphes event files
3*/
4
5#ifdef __CLING__
6R__LOAD_LIBRARY(libDelphes)
7#include "classes/DelphesClasses.h"
8#include "external/ExRootAnalysis/ExRootTreeReader.h"
9#include "external/ExRootAnalysis/ExRootResult.h"
10#else
11class ExRootTreeReader;
12class ExRootResult;
13#endif
14
15
16template<typename T>
17void CollectionFilter(const TClonesArray& inColl ,vector<T*>& outColl, Double_t ptMin=30, Double_t etaMax=2.5)
18{
19
20 const TObject *object;
21
22 for (Int_t i = 0; i < inColl.GetEntriesFast(); i++)
23 {
24
25 object = inColl.At(i);
26 const T *t = static_cast<const T*>(object);
27
28 if(t->P4().Pt() < ptMin) continue;
29 if(TMath::Abs(t->P4().Eta()) > etaMax) continue;
30
31 outColl.push_back(t);
32
33 }
34}
35
36
37//------------------------------------------------------------------------------
38
39void macro(const char *inputFile1, const char *inputFile2)
40{
41 gSystem->Load("libDelphes");
42
43 TString inputFile[2];
44 inputFile[0] = TString(inputFile1);
45 inputFile[1] = TString(inputFile2);
46
47 TH1F *hMass[2];
48 TH1F *hTau21[2];
49
50
51 for (Int_t i = 0; i < 2; i++)
52 {
53
54 // Create chain of root trees
55 TChain chain("Delphes");
56 chain.Add(inputFile[i]);
57
58 // Create object of class ExRootTreeReader
59 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
60 Long64_t numberOfEntries = treeReader->GetEntries();
61
62 // Get pointers to branches used in this analysis
63 TClonesArray *branchJet = treeReader->UseBranch("Jet");
64 TClonesArray *branchElectron = treeReader->UseBranch("Electron");
65 TClonesArray *branchMuon = treeReader->UseBranch("Muon");
66
67 hMass[i] = new TH1F("mass", "", 25, 0.0, 500.0);
68 hTau21[i] = new TH1F("tau21", "", 25, 0.0, 1.0);
69
70 // Define variables
71 const Jet *jet;
72
73 vector<const Electron*> *electrons = new vector<const Electron*>();
74 vector<const Muon*> *muons = new vector<const Muon*>();
75 vector<const Jet*> *jets = new vector<const Jet*>();
76
77 Double_t mass, tau21;
78
79 // Loop over all events
80 for(Int_t entry = 0; entry < numberOfEntries; ++entry)
81 {
82 // Load selected branches with data from specified event
83 treeReader->ReadEntry(entry);
84
85 electrons -> clear();
86 muons -> clear();
87 jets -> clear();
88
89 // Select leptons with pT > 10 GeV and |eta| < 2.5
90 CollectionFilter(*branchElectron, *electrons , 20.0 , 2.5);
91 CollectionFilter(*branchMuon , *muons , 20.0 , 2.5);
92
93 // Select jet with pT > 400 GeV and |eta| < 2.5
94 CollectionFilter(*branchJet , *jets , 300.0 , 2.5);
95
96 // Select event with 0 or 1 leptons at most
97 if(electrons->size() + muons->size() != 0) continue;
98
99 // Select event with at least 1 hard jet
100 if(jets->size() == 0) continue;
101
102 jet = jets->at(0);
103
104 // Define electron pair invariant mass
105 mass = jet->Mass;
106 tau21 = jet->Tau[0] > 0 ? jet->Tau[1]/jet->Tau[0] : -999;
107
108 // Fill the above quantities in histograms
109 hMass[i] -> Fill(mass, 1/double(numberOfEntries));
110 if(tau21 > 0) hTau21[i] -> Fill(tau21, 1/double(numberOfEntries));
111
112 } // end of event loop*/
113 }
114
115 // Plot all the above observables
116
117 TCanvas *c1 = new TCanvas("c1","multipads",1200,600);
118 gStyle->SetOptStat(0);
119 c1->Divide(2,1);
120
121 c1->cd(1);
122
123 hMass[0] -> SetLineWidth(3);
124 hMass[0] -> SetLineColor(kBlue);
125
126 hMass[0] -> GetXaxis()->SetTitle("m [GeV]");
127
128 hMass[1] -> SetLineWidth(3);
129 hMass[1] -> SetLineColor(kRed);
130
131 hMass[0] -> Draw("hist");
132 hMass[1] -> Draw("hist same");
133
134 hMass[0] -> SetMaximum(0.50);
135
136 TLegend *leg1 = new TLegend(0.32,0.65,0.70,0.85);
137
138 leg1->AddEntry(hMass[0], "light jet", "l");
139 leg1->AddEntry(hMass[1], "w jet", "l");
140
141 leg1->SetBorderSize(0);
142 leg1->SetShadowColor(0);
143 leg1->SetFillColor(0);
144 leg1->Draw();
145
146 c1->cd(2);
147
148 hTau21[0] -> SetLineWidth(3);
149 hTau21[0] -> SetLineColor(kBlue);
150
151 hTau21[0] -> GetXaxis()->SetTitle("#tau_{2,1}");
152
153 hTau21[1] -> SetLineWidth(3);
154 hTau21[1] -> SetLineColor(kRed);
155
156 hTau21[0] -> Draw("hist");
157 hTau21[1] -> Draw("hist same");
158
159 hTau21[0] -> SetMaximum(0.2);
160
161 TLegend *leg2 = new TLegend(0.55,0.65,0.89,0.85);
162
163 leg2->AddEntry(hTau21[0], "light jet", "l");
164 leg2->AddEntry(hTau21[1], "w jet", "l");
165
166 leg2->SetBorderSize(0);
167 leg2->SetShadowColor(0);
168 leg2->SetFillColor(0);
169 leg2->Draw();
170
171 TString fname(inputFile2);
172 fname.ReplaceAll("root", "png");
173
174 c1->Print(fname,"png");
175
176}