ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dasu/UltraFastSim/bbHAnalysis.cc
(Generate patch)

Comparing UserCode/dasu/UltraFastSim/bbHAnalysis.cc (file contents):
Revision 1.4 by dasu, Sun Feb 20 02:01:28 2011 UTC vs.
Revision 1.7 by dasu, Fri Feb 25 15:37:18 2011 UTC

# Line 1 | Line 1
1 +
2   #include "bbHAnalysis.h"
3  
4   #include <iostream>
# Line 12 | Line 13 | using namespace Pythia8;
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();
# Line 51 | Line 53 | bbHAnalysis::bbHAnalysis(TFile *o, Pythi
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() {
# Line 65 | Line 121 | bbHAnalysis::~bbHAnalysis() {
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_) {
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;
# Line 83 | Line 140 | bool bbHAnalysis::run() {
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();
# Line 96 | Line 156 | bool bbHAnalysis::run() {
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.;
# Line 123 | Line 184 | bool bbHAnalysis::run() {
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();
# Line 130 | Line 197 | bool bbHAnalysis::run() {
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 {
# Line 167 | Line 318 | bool bbHAnalysis::run() {
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) {
# Line 180 | Line 339 | bool bbHAnalysis::run() {
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();
# Line 197 | Line 370 | bool bbHAnalysis::run() {
370    }
371  
372    if(nBJets > 0) {
373 <    bJet1Pt = ufs->bJetList()[0].Et();
374 <    bJet1Eta = ufs->bJetList()[0].rap();
375 <    bJet1Phi = ufs->bJetList()[0].phi();
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].Et();
378 <      bJet2Eta = ufs->bJetList()[1].rap();
379 <      bJet2Phi = ufs->bJetList()[1].phi();
377 >      bJet2Pt = ufs->bJetList()[1].Pt();
378 >      bJet2Eta = ufs->bJetList()[1].Eta();
379 >      bJet2Phi = ufs->bJetList()[1].Phi();
380      }
381    }
382  
383 <  tree->Fill();
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   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines