Fork me on GitHub

Ticket #911: SRACodeB.C

File SRACodeB.C, 5.5 KB (added by J Dutta, 9 years ago)
Line 
1#include "TH1.h"
2#include "TSystem.h"
3#include <fstream>
4#include <iostream>
5#include <cmath>
6#include<cstdlib>
7
8using namespace std;
9
10template <typename T>
11inline T sqr(T a){return(a*a);}
12
13template <typename T>
14inline T sum(T a,T b){return(a+b);}
15
16template <typename Q>
17inline Q sqrsum(Q a,Q b){return(sqr(a)+sqr(b));}
18
19template <typename T>
20T DPhi(T a,T b);
21
22template <typename T>
23T DEta(T a,T b);
24
25template <typename T>
26T DR(T a,T b,T c,T d);
27
28
29struct MyPlots
30{
31
32 TH1 *fNLepton;
33 TH1 *fNJets;
34 TH1 *fNJetsI;
35 TH1 *fNPhoton;
36
37};
38
39class ExRootResult;
40
41class ExRootTreeReader;
42
43void BookHistograms(ExRootResult *results,MyPlots *plots);
44
45void PrintHistograms(ExRootResult *results,MyPlots *plots);
46
47void AnalyzeEvents(ExRootTreeReader *treeReader,MyPlots *plots);
48
49void SRACodeB(const char *inputFile1)
50{
51 gSystem->Load("libDelphes");
52
53 TChain *chain = new TChain("Delphes");
54
55 chain->Add(inputFile1);
56
57 ExRootTreeReader *treeReader = new ExRootTreeReader(chain);
58 ExRootResult *results = new ExRootResult();
59
60 MyPlots *plots = new MyPlots;
61
62 BookHistograms(results,plots);
63
64 AnalyzeEvents(treeReader,plots);
65
66 PrintHistograms(results,plots);
67
68 results->Write("CodeResults.root");
69
70 cout<<"***Exiting....."<<endl;
71
72 delete plots;
73 delete results;
74 delete treeReader;
75 delete chain;
76}
77
78void AnalyzeEvents(ExRootTreeReader *treeReader,MyPlots *plots)
79{
80
81
82 TClonesArray *branchElectron = treeReader->UseBranch("Electron");
83 TClonesArray *branchMuon = treeReader->UseBranch("Muon");
84 TClonesArray *branchJet = treeReader->UseBranch("Jet");
85 TClonesArray *branchPhoton = treeReader->UseBranch("Photon");
86
87 Long64_t allEntries = treeReader->GetEntries();
88 ofstream f1("CodeRun.dat");
89 f1<<"No. of Entries in sample = "<<allEntries<<endl;
90
91 Electron *el[30];
92 Muon *mu[30];
93 Jet *jets[30];
94 Photon *photon[30];
95
96 Long64_t entry;
97
98 cout<<"allEntries = "<<allEntries<<endl;
99 f1<<"allEntries = "<<allEntries<<endl;
100
101 Int_t i(0),j(0);
102
103 for(entry=0;entry<allEntries;entry++)
104 {
105 treeReader->ReadEntry(entry);
106
107 Int_t sl(0),sj(0),sp(0);
108
109 for(i=0;i<(branchElectron->GetEntriesFast());i++)
110 {
111 el[i]=(Electron*) branchElectron->At(i);
112 if((el[i]->PT)>10 && fabs(el[i]->Eta)<2.47 && (fabs(el[i]->Eta)<1.37 || fabs(el[i]->Eta)>1.52))
113 {
114 ++sl;
115
116 }
117
118
119 }
120
121 for(i=0;i<(branchMuon->GetEntriesFast());i++)
122 {
123 mu[i]=(Muon*) branchMuon->At(i);
124 if((mu[i]->PT)>10 && fabs(mu[i]->Eta)<2.4 && (fabs(mu[i]->Eta)<1.37 || fabs(mu[i]->Eta)>1.52))
125 {
126 ++sl;
127
128 }
129
130
131 }
132
133
134 for(i=0;i<(branchPhoton->GetEntriesFast());i++)
135 {
136 photon[i]=(Photon*) branchPhoton->At(i);
137 if((photon[i]->PT)>50 && fabs(photon[i]->Eta)<2.47 && (fabs(photon[i]->Eta)<1.37 || fabs(photon[i]->Eta)>1.52))
138 {
139 ++sp;
140
141 }
142
143
144 }
145
146 ofstream Jf1("JetPT.dat");
147 ofstream Jf2("JetPhi.dat");
148 ofstream Jf3("JetEta.dat");
149
150
151 for(i=0;i<(branchJet->GetEntriesFast());i++)
152 {
153 jets[i]=(Jet*) branchJet->At(i);
154 Double_t eta=jets[i]->Eta;
155 Double_t pt =jets[i]->PT;
156
157 if(pt>40 && fabs(eta)<2.5)
158 {
159
160 ++sj;
161 Jf1<<(jets[i]->PT)<<endl;
162 Jf2<<(jets[i]->Phi)<<endl;
163 Jf3<<(jets[i]->Eta)<<endl;
164
165
166 }
167
168
169 }
170
171
172 ifstream JIf1("JetPT.dat");
173 ifstream JIf2("JetPhi.dat");
174 ifstream JIf3("JetEta.dat");
175
176
177 Double_t * JPT = new Double_t [sj];
178 Double_t * JPhi = new Double_t [sj];
179 Double_t * JEta = new Double_t [sj];
180 Double_t * JMass = new Double_t [sj];
181
182 for(i=0;i<sj;i++)
183 {
184 JIf1>>JPT[i];
185 JIf2>>JPhi[i];
186 JIf3>>JEta[i];
187
188 }
189
190 Int_t isj(0);
191
192 for(i=0;i<sj;i++)
193 {
194 Int_t x(0);
195 for(j=0;j<sj;j++)
196 {
197 if(j!=i)
198 {
199 Double_t DRjj = DR(JPhi[i],JEta[i],JPhi[j],JEta[j]);
200
201 if(DRjj>0.4) ++x;
202
203 }
204
205 }
206
207
208
209 if(x==(sj-1))
210 {
211 ++isj;
212
213 }
214
215
216 }
217
218
219 plots->fNJetsI->Fill(isj);
220 plots->fNLepton->Fill(sl);
221 plots->fNPhoton->Fill(sp);
222
223
224
225
226 }
227
228
229}
230void BookHistograms(ExRootResult *results,MyPlots *plots)
231{
232 THStack *stack;
233 TLegend *legend;
234 TPaveText *comment;
235
236
237// books a histogram for multiplicites : NLepton,NJet,NPhoton
238
239 plots->fNLepton = results->AddHist1D(
240 "LeptonMultiplicity","Lepton Multiplicity",
241 "N_{Lepton}","Number of Events",
242 10,-0.5,9.5);
243
244// book general comment
245
246 comment = results->AddComment(0.64,0.86,0.98,0.98);
247 comment->AddText("demonstration plot ");
248 comment->AddText("produced by ex5.C");
249
250//Show histogram statistics :
251
252 plots->fNLepton->SetStats();
253
254// books a histogram for multiplicites : NJetsI
255
256 plots->fNJetsI = results->AddHist1D(
257 "IJetMultiplicity","IJet Multiplicity",
258 "N_{I_Jets}","Number of Events",
259 20,-0.5,19.5);
260
261 plots->fNJetsI->SetStats();
262
263// books a histogram for multiplicites : NPhoton
264
265 plots->fNPhoton = results->AddHist1D(
266 "PhotonMultiplicity","Photon Multiplicity",
267 "N_{Photon}","Number of Events",
268 10,-0.5,9.5);
269
270 plots->fNPhoton->SetStats();
271
272}
273
274
275template <typename T>
276T DPhi(T a,T b)
277{
278 T c = fabs(a-b);
279
280 Double_t PI = 4.0*atan(1.0);
281
282 if(c>PI) c = 2.0*PI-c;
283 return(c);
284}
285
286template <typename T>
287T DEta(T a,T b)
288{
289 T c = fabs(a-b);
290
291 return(c);
292}
293
294template <typename T>
295T DR(T a,T b,T c,T d)
296{
297
298 T e = DPhi(a,c);
299
300 T f = DEta(b,d);
301
302 T dr = sqrt(sqrsum(e,f));
303
304 return(dr);
305
306}
307
308
309void PrintHistograms(ExRootResult *results, MyPlots *plots)
310{
311
312 results->Print("C");
313
314 results->Print("jpeg");
315}
316