ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dasu/UltraFastSim/bbHAnalysis.cc
Revision: 1.7
Committed: Fri Feb 25 15:37:18 2011 UTC (14 years, 2 months ago) by dasu
Content type: text/plain
Branch: MAIN
CVS Tags: v2
Changes since 1.6: +6 -6 lines
Log Message:
Moved to storing the entire UFS data

File Contents

# User Rev Content
1 mmhl 1.5
2 dasu 1.1 #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 mmhl 1.5 #include "TH1D.h"
17 dasu 1.1
18 dasu 1.3 bbHAnalysis::bbHAnalysis(TFile *o, Pythia *p, UltraFastSim *u, bool v) : outFile(o), pythia(p), ufs(u), verbose_(v), iEvent(0) {
19     outFile->cd();
20 dasu 1.1 tree = new TTree("bbHTree","Tree containing bbH -> tau_mu, tau_h analysis objects");
21 dasu 1.2 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 dasu 1.1 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 dasu 1.2 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 dasu 1.1 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 mmhl 1.6
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 dasu 1.2 tree->Branch("gVisMass", &gVisMass, "gVisMass/D");
66     tree->Branch("rVisMass", &rVisMass, "rVisMass/D");
67 dasu 1.1 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 mmhl 1.5
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 mmhl 1.6 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 mmhl 1.5 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 mmhl 1.6 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 mmhl 1.5
90 mmhl 1.6 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 mmhl 1.5
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 mmhl 1.6 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 mmhl 1.5
115    
116 dasu 1.1 }
117    
118     bbHAnalysis::~bbHAnalysis() {
119 dasu 1.3 }
120    
121     bool bbHAnalysis::end() {
122     outFile->cd();
123 dasu 1.1 tree->Write();
124 mmhl 1.5 outFile->Close();
125 dasu 1.3 return true;
126 dasu 1.1 }
127    
128     bool bbHAnalysis::run() {
129    
130 mmhl 1.6 if(verbose_) {
131 dasu 1.1 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 dasu 1.2 cout << "Number of GenVisTau = " << ufs->visTauList().size() << endl;
138     cout << "Number of Rec Taus = " << ufs->tauList().size() << endl;
139 dasu 1.1 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 mmhl 1.6 cout << "Number of bJetsStdGeom= " << ufs->bJetListStdGeom().size()<<endl;
144 dasu 1.1 cout << endl;
145     }
146    
147 mmhl 1.5
148    
149 dasu 1.1 // Store information for top muon, hadronic tau (not matching muon), top bJets (not matching tau)
150    
151 dasu 1.2 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 dasu 1.1 nJets = ufs->jetList().size();
158     nBJets = ufs->bJetList().size();
159 mmhl 1.6 nBJetsStdGeom=ufs->bJetListStdGeom().size();
160 dasu 1.1 muonPt = 0;
161     muonEta = -999.;
162     muonPhi = -999.;
163     tauPt = 0;
164     tauEta = -999.;
165     tauPhi = -999.;
166 dasu 1.2 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 dasu 1.1 bJet1Pt = 0;
182     bJet1Eta = -999.;
183     bJet1Phi = -999.;
184     bJet2Pt = 0;
185     bJet2Eta = -999.;
186     bJet2Phi = -999.;
187 mmhl 1.6 bJet1PtStdGeom = 0;
188     bJet1EtaStdGeom = -999.;
189     bJet1PhiStdGeom = -999.;
190     bJet2PtStdGeom = 0;
191     bJet2EtaStdGeom = -999.;
192     bJet2PhiStdGeom = -999.;
193 dasu 1.2 gVisMass = -999.;
194     rVisMass = -999.;
195 dasu 1.4 ET = ufs->getET();
196     HT = ufs->getHT();
197     MET = ufs->getMET().Pt();
198     MHT = ufs->getMHT().Pt();
199 dasu 1.1
200 mmhl 1.5
201     //ml added
202 mmhl 1.6 nbJetCut=0;//passed bjet cut- 1 is passed
203     nbJetCutStdGeom=0;
204     nbasicCut=0;
205     nbasicCutStdGeom=0;
206    
207     bJetPtCut=30.;
208 mmhl 1.5 muonE=0;
209     muonPx=0;
210     muonPy=0;
211     muonPz=0;
212 mmhl 1.6 taumE=0;
213     taumPx=0;
214     taumPy=0;
215     taumPz=0;
216 mmhl 1.5
217     hadE=0;
218     hadPx=0;
219     hadPy=0;
220     hadPz=0;
221 mmhl 1.6 tauhE=0;
222     tauhPx=0;
223     tauhPy=0;
224     tauhPz=0;
225 mmhl 1.5
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 mmhl 1.6 rTau1E=0;
236     rTau1Px=0;
237     rTau1Py=0;
238     rTau1Pz=0;
239     rTau2E=0;
240     rTau2Px=0;
241     rTau2Py=0;
242     rTau2Pz=0;
243 mmhl 1.5
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 dasu 1.2 int muonMother;
259     TLorentzVector muonP;
260 dasu 1.1 for(unsigned int i = 0; i < ufs->muonList().size(); i++) {
261 mmhl 1.5 //loops over all the muons and stores highest pt one
262 dasu 1.1 if(ufs->muonList()[i].Pt() > muonPt) {
263 dasu 1.2 muonMother = ufs->muonList()[i].GetFirstMother();
264 dasu 1.1 muonPt = ufs->muonList()[i].Pt();
265     muonEta = ufs->muonList()[i].Eta();
266     muonPhi = ufs->muonList()[i].Phi();
267 dasu 1.2 ufs->muonList()[i].Momentum(muonP);
268 mmhl 1.5
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 dasu 1.1 }
275     }
276    
277 mmhl 1.5 //this is the loop over all generated taus
278 dasu 1.2 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 mmhl 1.5
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 dasu 1.2 }
296     }
297    
298     TLorentzVector taumP;
299     TLorentzVector tauhP;
300     for(unsigned int i = 0; i < ufs->visTauList().size(); i++) {
301 mmhl 1.5 //loop over all visible taus that decay to lepton and hadron
302 dasu 1.2 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 mmhl 1.5 //ml added
309 mmhl 1.6 taumE=ufs->visTauList()[i].Energy();
310     taumPx=ufs->visTauList()[i].Px();
311     taumPy=ufs->visTauList()[i].Py();
312     taumPz=ufs->visTauList()[i].Pz();
313 dasu 1.2 }
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 mmhl 1.5
322 mmhl 1.6 tauhE=ufs->visTauList()[i].Energy();
323     tauhPx=ufs->visTauList()[i].Px();
324     tauhPy=ufs->visTauList()[i].Py();
325     tauhPz=ufs->visTauList()[i].Pz();
326 dasu 1.2 }
327     }
328     }
329 mmhl 1.5
330 dasu 1.2 TLorentzVector gTauTauP = taumP + tauhP;
331     gVisMass = gTauTauP.M();
332    
333 mmhl 1.5
334     //full tau list highest taus rtau final state (i.e. no hadrons)
335 dasu 1.2 TLorentzVector rTauP;
336 dasu 1.1 for(unsigned int i = 0; i < ufs->tauList().size(); i++) {
337 dasu 1.2 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 mmhl 1.5
343 mmhl 1.6 //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 mmhl 1.5
354 dasu 1.2 }
355     }
356     TLorentzVector rTauTauP = muonP + rTauP;
357     rVisMass = rTauTauP.M();
358    
359 mmhl 1.5
360     //save the highest two pt of bquarks etc
361 dasu 1.2 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 dasu 1.1 }
370     }
371    
372     if(nBJets > 0) {
373 dasu 1.7 bJet1Pt = ufs->bJetList()[0].Pt();
374     bJet1Eta = ufs->bJetList()[0].Eta();
375     bJet1Phi = ufs->bJetList()[0].Phi();
376 dasu 1.1 if(nBJets > 1) {
377 dasu 1.7 bJet2Pt = ufs->bJetList()[1].Pt();
378     bJet2Eta = ufs->bJetList()[1].Eta();
379     bJet2Phi = ufs->bJetList()[1].Phi();
380 dasu 1.1 }
381     }
382 mmhl 1.6
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 mmhl 1.5 //ml added - passes secondary pt cut - saves 1 for pass, 0 otherwise
394     if(nMuons >0 && nTaus >1 && nBJets > 1){
395 mmhl 1.6 nbasicCut = 1;
396 mmhl 1.5
397     }
398     if(nBJets>1){
399 mmhl 1.6 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 mmhl 1.5 }
411     }
412 mmhl 1.6
413    
414    
415 dasu 1.1 tree->Fill();
416     iEvent++;
417     return true;
418     }
419