ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dasu/UltraFastSim/bbHAnalysis.cc
Revision: 1.8
Committed: Thu Jun 13 15:46:31 2013 UTC (11 years, 11 months ago) by dasu
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.7: +0 -0 lines
State: FILE REMOVED
Error occurred while calculating annotation data.
Log Message:
Pruned to the minimum needed for future

File Contents

# Content
1
2 #include "bbHAnalysis.h"
3
4 #include <iostream>
5 #include <vector>
6 #include <string>
7 using namespace std;
8
9 #include "Pythia.h"
10 using namespace Pythia8;
11
12 #include "UltraFastSim.h"
13
14 #include "TFile.h"
15 #include "TTree.h"
16 #include "TH1D.h"
17
18 bbHAnalysis::bbHAnalysis(TFile *o, Pythia *p, UltraFastSim *u, bool v) : outFile(o), pythia(p), ufs(u), verbose_(v), iEvent(0) {
19 outFile->cd();
20 tree = new TTree("bbHTree","Tree containing bbH -> tau_mu, tau_h analysis objects");
21 tree->Branch("nElecs", &nElecs, "nElecs/I");
22 tree->Branch("nMuons", &nMuons, "nMuons/I");
23 tree->Branch("nTaus", &nTaus, "nTaus/I");
24 tree->Branch("nRTaus", &nRTaus, "nRTaus/I");
25 tree->Branch("nbQrks", &nbQrks, "nbQrks/I");
26 tree->Branch("ncQrks", &ncQrks, "ncQrks/I");
27 tree->Branch("nJets", &nJets, "nJets/I");
28 tree->Branch("nBJets", &nBJets, "nBJets/I");
29 tree->Branch("muonPt", &muonPt, "muonPt/D");
30 tree->Branch("muonEta", &muonEta, "muonEta/D");
31 tree->Branch("muonPhi", &muonPhi, "muonPhi/D");
32 tree->Branch("tauPt", &tauPt, "tauPt/D");
33 tree->Branch("tauEta", &tauEta, "tauEta/D");
34 tree->Branch("tauPhi", &tauPhi, "tauPhi/D");
35 tree->Branch("taumPt", &taumPt, "taumPt/D");
36 tree->Branch("taumEta", &taumEta, "taumEta/D");
37 tree->Branch("taumPhi", &taumPhi, "taumPhi/D");
38 tree->Branch("tauhPt", &tauhPt, "tauhPt/D");
39 tree->Branch("tauhEta", &tauhEta, "tauhEta/D");
40 tree->Branch("tauhPhi", &tauhPhi, "tauhPhi/D");
41 tree->Branch("rTauPt", &rTauPt, "rTauPt/D");
42 tree->Branch("rTauEta", &rTauEta, "rTauEta/D");
43 tree->Branch("rTauPhi", &rTauPhi, "rTauPhi/D");
44 tree->Branch("bQrk1Pt", &bQrk1Pt, "bQrk1Pt/D");
45 tree->Branch("bQrk1Eta", &bQrk1Eta, "bQrk1Eta/D");
46 tree->Branch("bQrk1Phi", &bQrk1Phi, "bQrk1Phi/D");
47 tree->Branch("bQrk2Pt", &bQrk2Pt, "bQrk2Pt/D");
48 tree->Branch("bQrk2Eta", &bQrk2Eta, "bQrk2Eta/D");
49 tree->Branch("bQrk2Phi", &bQrk2Phi, "bQrk2Phi/D");
50 tree->Branch("bJet1Pt", &bJet1Pt, "bJet1Pt/D");
51 tree->Branch("bJet1Eta", &bJet1Eta, "bJet1Eta/D");
52 tree->Branch("bJet1Phi", &bJet1Phi, "bJet1Phi/D");
53 tree->Branch("bJet2Pt", &bJet2Pt, "bJet2Pt/D");
54 tree->Branch("bJet2Eta", &bJet2Eta, "bJet2Eta/D");
55 tree->Branch("bJet2Phi", &bJet2Phi, "bJet2Phi/D");
56
57 tree->Branch("bJet1PtStdGeom", &bJet1PtStdGeom, "bJet1PtStdGeom/D");
58 tree->Branch("bJet1EtaStdGeom", &bJet1EtaStdGeom, "bJet1EtaStdGeom/D");
59 tree->Branch("bJet1PhiStdGeom", &bJet1PhiStdGeom, "bJet1PhiStdGeom/D");
60 tree->Branch("bJet2PtStdGeom", &bJet2PtStdGeom, "bJet2PtStdGeom/D");
61 tree->Branch("bJet2EtaStdGeom", &bJet2EtaStdGeom, "bJet2EtaStdGeom/D");
62 tree->Branch("bJet2PhiStdGeom", &bJet2PhiStdGeom, "bJet2PhiStdGeom/D");
63
64
65 tree->Branch("gVisMass", &gVisMass, "gVisMass/D");
66 tree->Branch("rVisMass", &rVisMass, "rVisMass/D");
67 tree->Branch("ET", &ET, "ET/D");
68 tree->Branch("HT", &HT, "HT/D");
69 tree->Branch("MET", &MET, "MET/D");
70 tree->Branch("MHT", &MHT, "MHT/D");
71
72 //ml added branches for invariant mass analysis
73 tree->Branch("muonE", &muonE,"muonE/D");
74 tree->Branch("muonPx", &muonPx,"muonPx/D");
75 tree->Branch("muonPy", &muonPy,"muonPy/D");
76 tree->Branch("muonPz", &muonPz,"muonPz/D");
77 tree->Branch("taumE", &taumE,"taumuonE/D");
78 tree->Branch("taumPx", &taumPx,"taumPx/D");
79 tree->Branch("taumPy", &taumPy,"taumPy/D");
80 tree->Branch("taumPz", &taumPz,"taumPz/D");
81 tree->Branch("hadE", &hadE,"hadE/D");
82 tree->Branch("hadPx", &hadPx,"hadPx/D");
83 tree->Branch("hadPy", &hadPy,"hadPy/D");
84 tree->Branch("hadPz", &hadPz,"hadPz/D");
85 tree->Branch("tauhE", &tauhE,"tauhE/D");
86 tree->Branch("tauhPx", &tauhPx,"tauhPx/D");
87 tree->Branch("tauhPy", &tauhPy,"tauhPy/D");
88 tree->Branch("tauhPz", &tauhPz,"tauPz/D");
89
90 tree->Branch("nBJetsStdGeom", &nBJetsStdGeom,"nBJetStdsGeom/I");
91
92 tree->Branch("nbJetCut", &nbJetCut,"nbJetCut/I");
93 tree->Branch("nbJetCutStdGeom", &nbJetCutStdGeom,"nbJetCutStdGeom/I");
94 tree->Branch("nbasicCut",&nbasicCut,"nbasicCut/I");
95 tree->Branch("nbasicCutStdGeom",&nbasicCutStdGeom,"nbasicCutStdGeom/I");
96
97 tree->Branch("tau1E", &tau1E,"tau1E/D");
98 tree->Branch("tau1Px", &tau1Px,"tau1Px/D");
99 tree->Branch("tau1Py", &tau1Py,"tau1Py/D");
100 tree->Branch("tau1Pz", &tau1Pz,"tau1Pz/D");
101 tree->Branch("tau2E", &tau2E,"tau2E/D");
102 tree->Branch("tau2Px", &tau2Px,"tau2Px/D");
103 tree->Branch("tau2Py", &tau2Py,"tau2Py/D");
104 tree->Branch("tau2Pz", &tau2Pz,"tau2Pz/D");
105
106 tree->Branch("rTau1E", &rTau1E,"rTau1E/D");
107 tree->Branch("rTau1Px", &rTau1Px,"rTau1Px/D");
108 tree->Branch("rTau1Py", &rTau1Py,"rTau1Py/D");
109 tree->Branch("rTau1Pz", &rTau1Pz,"rTau1Pz/D");
110 tree->Branch("rTau2E", &rTau2E,"rTau2E/D");
111 tree->Branch("rTau2Px", &rTau2Px,"rTau2Px/D");
112 tree->Branch("rTau2Py", &rTau2Py,"rTau2Py/D");
113 tree->Branch("rTau2Pz", &rTau2Pz,"rTau2Pz/D");
114
115
116 }
117
118 bbHAnalysis::~bbHAnalysis() {
119 }
120
121 bool bbHAnalysis::end() {
122 outFile->cd();
123 tree->Write();
124 outFile->Close();
125 return true;
126 }
127
128 bool bbHAnalysis::run() {
129
130 if(verbose_) {
131 cout << "Number of Gen Chrgd = " << ufs->chargedHadronList().size() << endl;
132 cout << "Number of Gen Neutr = " << ufs->neutralHadronList().size() << endl;
133 cout << "Number of Gen Phtns = " << ufs->photonList().size() << endl;
134 cout << "Number of Gen Elcns = " << ufs->electronList().size() << endl;
135 cout << "Number of Gen Muons = " << ufs->muonList().size() << endl;
136 cout << "Number of Gen Taus = " << ufs->genTauList().size() << endl;
137 cout << "Number of GenVisTau = " << ufs->visTauList().size() << endl;
138 cout << "Number of Rec Taus = " << ufs->tauList().size() << endl;
139 cout << "Number of c Quarks = " << ufs->cQuarkList().size() << endl;
140 cout << "Number of b Quarks = " << ufs->bQuarkList().size() << endl;
141 cout << "Number of Jets = " << ufs->jetList().size() << endl;
142 cout << "Number of bJets = " << ufs->bJetList().size() << endl;
143 cout << "Number of bJetsStdGeom= " << ufs->bJetListStdGeom().size()<<endl;
144 cout << endl;
145 }
146
147
148
149 // Store information for top muon, hadronic tau (not matching muon), top bJets (not matching tau)
150
151 nElecs = ufs->electronList().size();
152 nMuons = ufs->muonList().size();
153 nTaus = ufs->genTauList().size();
154 nRTaus = ufs->tauList().size();
155 nbQrks = ufs->bQuarkList().size();
156 ncQrks = ufs->cQuarkList().size();
157 nJets = ufs->jetList().size();
158 nBJets = ufs->bJetList().size();
159 nBJetsStdGeom=ufs->bJetListStdGeom().size();
160 muonPt = 0;
161 muonEta = -999.;
162 muonPhi = -999.;
163 tauPt = 0;
164 tauEta = -999.;
165 tauPhi = -999.;
166 taumPt = 0;
167 taumEta = -999.;
168 taumPhi = -999.;
169 tauhPt = 0;
170 tauhEta = -999.;
171 tauhPhi = -999.;
172 rTauPt = 0;
173 rTauEta = -999.;
174 rTauPhi = -999.;
175 bQrk1Pt = 0;
176 bQrk1Eta = -999.;
177 bQrk1Phi = -999.;
178 bQrk2Pt = 0;
179 bQrk2Eta = -999.;
180 bQrk2Phi = -999.;
181 bJet1Pt = 0;
182 bJet1Eta = -999.;
183 bJet1Phi = -999.;
184 bJet2Pt = 0;
185 bJet2Eta = -999.;
186 bJet2Phi = -999.;
187 bJet1PtStdGeom = 0;
188 bJet1EtaStdGeom = -999.;
189 bJet1PhiStdGeom = -999.;
190 bJet2PtStdGeom = 0;
191 bJet2EtaStdGeom = -999.;
192 bJet2PhiStdGeom = -999.;
193 gVisMass = -999.;
194 rVisMass = -999.;
195 ET = ufs->getET();
196 HT = ufs->getHT();
197 MET = ufs->getMET().Pt();
198 MHT = ufs->getMHT().Pt();
199
200
201 //ml added
202 nbJetCut=0;//passed bjet cut- 1 is passed
203 nbJetCutStdGeom=0;
204 nbasicCut=0;
205 nbasicCutStdGeom=0;
206
207 bJetPtCut=30.;
208 muonE=0;
209 muonPx=0;
210 muonPy=0;
211 muonPz=0;
212 taumE=0;
213 taumPx=0;
214 taumPy=0;
215 taumPz=0;
216
217 hadE=0;
218 hadPx=0;
219 hadPy=0;
220 hadPz=0;
221 tauhE=0;
222 tauhPx=0;
223 tauhPy=0;
224 tauhPz=0;
225
226 tau1E=0;
227 tau1Px=0;
228 tau1Py=0;
229 tau1Pz=0;
230 tau2E=0;
231 tau2Px=0;
232 tau2Py=0;
233 tau2Pz=0;
234
235 rTau1E=0;
236 rTau1Px=0;
237 rTau1Py=0;
238 rTau1Pz=0;
239 rTau2E=0;
240 rTau2Px=0;
241 rTau2Py=0;
242 rTau2Pz=0;
243
244 hadronPt = 0;
245
246
247 for(unsigned int i =0; i<ufs->chargedHadronList().size();i++){
248 if(ufs->chargedHadronList()[i].Pt() > hadronPt) {
249 hadE=ufs->chargedHadronList()[i].Energy();
250 hadPx=ufs->chargedHadronList()[i].Px();
251 hadPy=ufs->chargedHadronList()[i].Py();
252 hadPz=ufs->chargedHadronList()[i].Pz();
253
254 hadronPt=ufs->chargedHadronList()[i].Pt();
255 }
256 }
257
258 int muonMother;
259 TLorentzVector muonP;
260 for(unsigned int i = 0; i < ufs->muonList().size(); i++) {
261 //loops over all the muons and stores highest pt one
262 if(ufs->muonList()[i].Pt() > muonPt) {
263 muonMother = ufs->muonList()[i].GetFirstMother();
264 muonPt = ufs->muonList()[i].Pt();
265 muonEta = ufs->muonList()[i].Eta();
266 muonPhi = ufs->muonList()[i].Phi();
267 ufs->muonList()[i].Momentum(muonP);
268
269 //ml added
270 muonE = ufs->muonList()[i].Energy();
271 muonPx = ufs->muonList()[i].Px();
272 muonPy = ufs->muonList()[i].Py();
273 muonPz = ufs->muonList()[i].Pz();
274 }
275 }
276
277 //this is the loop over all generated taus
278 for(unsigned int i = 0; i < ufs->genTauList().size(); i++) {
279 if(ufs->genTauList()[i].Pt() > tauPt) {
280 tauPt = ufs->genTauList()[i].Pt();
281 tauEta = ufs->genTauList()[i].Eta();
282 tauPhi = ufs->genTauList()[i].Phi();
283
284 tau2E=tau1E;
285 tau2Px=tau1Px;
286 tau2Py=tau1Py;
287 tau2Pz=tau1Pz;
288
289 tau1E=ufs->genTauList()[i].Energy();
290 tau1Px=ufs->genTauList()[i].Px();
291 tau1Py=ufs->genTauList()[i].Py();
292 tau1Pz=ufs->genTauList()[i].Pz();
293
294 // cout<<tau1E<<" "<<tau2E<<endl;
295 }
296 }
297
298 TLorentzVector taumP;
299 TLorentzVector tauhP;
300 for(unsigned int i = 0; i < ufs->visTauList().size(); i++) {
301 //loop over all visible taus that decay to lepton and hadron
302 if(muonMother == ufs->visTauList()[i].GetSecondMother()) { // We stole the second mother spot to store the pythia particle index
303 if(ufs->visTauList()[i].Pt() > taumPt) {
304 taumPt = ufs->visTauList()[i].Pt();
305 taumEta = ufs->visTauList()[i].Eta();
306 taumPhi = ufs->visTauList()[i].Phi();
307 ufs->visTauList()[i].Momentum(taumP);
308 //ml added
309 taumE=ufs->visTauList()[i].Energy();
310 taumPx=ufs->visTauList()[i].Px();
311 taumPy=ufs->visTauList()[i].Py();
312 taumPz=ufs->visTauList()[i].Pz();
313 }
314 }
315 else {
316 if(ufs->visTauList()[i].Pt() > tauhPt) {
317 tauhPt = ufs->visTauList()[i].Pt();
318 tauhEta = ufs->visTauList()[i].Eta();
319 tauhPhi = ufs->visTauList()[i].Phi();
320 ufs->visTauList()[i].Momentum(tauhP);
321
322 tauhE=ufs->visTauList()[i].Energy();
323 tauhPx=ufs->visTauList()[i].Px();
324 tauhPy=ufs->visTauList()[i].Py();
325 tauhPz=ufs->visTauList()[i].Pz();
326 }
327 }
328 }
329
330 TLorentzVector gTauTauP = taumP + tauhP;
331 gVisMass = gTauTauP.M();
332
333
334 //full tau list highest taus rtau final state (i.e. no hadrons)
335 TLorentzVector rTauP;
336 for(unsigned int i = 0; i < ufs->tauList().size(); i++) {
337 if(ufs->tauList()[i].Pt() > rTauPt) {
338 rTauPt = ufs->tauList()[i].Pt();
339 rTauEta = ufs->tauList()[i].Eta();
340 rTauPhi = ufs->tauList()[i].Phi();
341 ufs->tauList()[i].Momentum(rTauP);
342
343 //save two highest momentum rtaus
344 rTau2E=rTau1E;
345 rTau2Px=rTau1Px;
346 rTau2Py=rTau1Py;
347 rTau2Pz=rTau1Pz;
348
349 rTau1E=ufs->tauList()[i].Energy();
350 rTau1Px=ufs->tauList()[i].Px();
351 rTau1Py=ufs->tauList()[i].Py();
352 rTau1Pz=ufs->tauList()[i].Pz();
353
354 }
355 }
356 TLorentzVector rTauTauP = muonP + rTauP;
357 rVisMass = rTauTauP.M();
358
359
360 //save the highest two pt of bquarks etc
361 if(nbQrks > 0) {
362 bQrk1Pt = ufs->bQuarkList()[0].Pt();
363 bQrk1Eta = ufs->bQuarkList()[0].Eta();
364 bQrk1Phi = ufs->bQuarkList()[0].Phi();
365 if(nbQrks > 1) {
366 bQrk2Pt = ufs->bQuarkList()[1].Pt();
367 bQrk2Eta = ufs->bQuarkList()[1].Eta();
368 bQrk2Phi = ufs->bQuarkList()[1].Phi();
369 }
370 }
371
372 if(nBJets > 0) {
373 bJet1Pt = ufs->bJetList()[0].Pt();
374 bJet1Eta = ufs->bJetList()[0].Eta();
375 bJet1Phi = ufs->bJetList()[0].Phi();
376 if(nBJets > 1) {
377 bJet2Pt = ufs->bJetList()[1].Pt();
378 bJet2Eta = ufs->bJetList()[1].Eta();
379 bJet2Phi = ufs->bJetList()[1].Phi();
380 }
381 }
382
383 if(nBJetsStdGeom > 0) {
384 bJet1PtStdGeom = ufs->bJetListStdGeom()[0].Et();
385 bJet1EtaStdGeom = ufs->bJetListStdGeom()[0].rap();
386 bJet1PhiStdGeom = ufs->bJetListStdGeom()[0].phi();
387 if(nBJetsStdGeom > 1) {
388 bJet2PtStdGeom = ufs->bJetListStdGeom()[1].Et();
389 bJet2EtaStdGeom = ufs->bJetListStdGeom()[1].rap();
390 bJet2PhiStdGeom = ufs->bJetListStdGeom()[1].phi();
391 }
392 }
393 //ml added - passes secondary pt cut - saves 1 for pass, 0 otherwise
394 if(nMuons >0 && nTaus >1 && nBJets > 1){
395 nbasicCut = 1;
396
397 }
398 if(nBJets>1){
399 if(bJet1Pt > bJetPtCut && bJet2Pt>bJetPtCut){
400 nbJetCut=1;
401 }
402 }
403
404 if(nMuons>0 && nTaus>1 && nBJetsStdGeom > 1){
405 nbasicCutStdGeom =1;
406 }
407 if(nBJetsStdGeom > 1){
408 if(bJet1PtStdGeom > bJetPtCut && bJet2PtStdGeom > bJetPtCut){
409 nbJetCutStdGeom=1;
410 }
411 }
412
413
414
415 tree->Fill();
416 iEvent++;
417 return true;
418 }
419