| 3 |
|
#include <iostream> |
| 4 |
|
#include <vector> |
| 5 |
|
#include <string> |
| 6 |
– |
|
| 6 |
|
using namespace std; |
| 7 |
|
|
| 8 |
+ |
#include "Pythia.h" |
| 9 |
+ |
using namespace Pythia8; |
| 10 |
+ |
|
| 11 |
|
#include "UltraFastSim.h" |
| 12 |
|
|
| 13 |
|
#include "TROOT.h" |
| 17 |
|
|
| 18 |
|
TROOT root("Mike's ROOT Session","PYTHIA Histograms"); |
| 19 |
|
|
| 20 |
< |
Analysis::Analysis(const char *r, const UltraFastSim *u) : runType(r), ufs(u) { |
| 20 |
> |
Analysis::Analysis(char *r, Pythia *p, UltraFastSim *u, bool v) : runType(r), pythia(p), ufs(u), verbose_(v), iEvent(0) { |
| 21 |
|
string name(runType); |
| 22 |
|
name += ".root"; |
| 23 |
|
outFile = TFile::Open(name.c_str(), "recreate"); |
| 22 |
– |
bool verbose_=false; |
| 23 |
– |
double _px,_py,_pz,_e,_m; |
| 24 |
– |
int _id,_np,_mother1,_mother2; //np is the number of particles... |
| 25 |
– |
|
| 24 |
|
bjet_pt = new TH1D("bjet_pt","bJet P_{T} [GeV]",20000,-.5,9999.5); |
| 25 |
|
tree = new TTree("T","All"); |
| 26 |
|
tree->Branch("_px",&_px,"_px/D"); |
| 34 |
|
tree->Branch("_mother2",&_mother2,"_mother2/I"); |
| 35 |
|
higgs_inv_mass = new TH1D("higgs_invm","Higgs Invariant Mass [GeV]", 20000, -0.5, 9999.5); |
| 36 |
|
|
| 39 |
– |
double _pxc,_pyc,_pzc,_ec,_mc; |
| 40 |
– |
int _idc,_npc,_mother1c,_mother2c; |
| 37 |
|
cuttree = new TTree("T_btagged","after b-tagging"); |
| 38 |
|
cuttree->Branch("_px",&_pxc,"_pxc/D"); |
| 39 |
|
cuttree->Branch("_py",&_pyc,"_pyc/D"); |
| 58 |
|
} |
| 59 |
|
|
| 60 |
|
bool Analysis::run() { |
| 65 |
– |
/* |
| 61 |
|
//ml added ------------------------------------------------------------// |
| 62 |
|
double bjet_ptcut(25.); //bjet pt cut. |
| 63 |
|
int bjet_counter25(0); |
| 64 |
|
int muon_counter(0),hadron_counter(0), tau_counter(0); |
| 65 |
< |
int mh_size = pythia.event.size(); //this is not a good idea as it will crash sometimes... |
| 65 |
> |
int mh_size = pythia->event.size(); //this is not a good idea as it will crash sometimes... |
| 66 |
|
double hadron_store[mh_size][4]; //obviously not a good solution... to set this.. |
| 67 |
|
double muon_store[mh_size][4]; |
| 68 |
|
|
| 86 |
|
_np=0; |
| 87 |
|
_npc=0; |
| 88 |
|
|
| 89 |
< |
vector<fastjet::PseudoJet> sortedbJets; //create a pseudojet |
| 95 |
< |
sortedbJets = sorted_by_pt(ufs->bJetList()); //psuedojet is filled with bjets that are sorted by pt |
| 96 |
< |
|
| 89 |
> |
vector<TLorentzVector> sortedbJets = ufs->bJetList(); |
| 90 |
|
// below is to work out what taulist and gentaulist is |
| 91 |
|
int countertau_fin(0), countertau_nf(0), countermu_fin(0), countermu_nf(0); |
| 92 |
< |
for (int i = 0; i < pythia.event.size(); ++i){ |
| 93 |
< |
Particle& theParticle = pythia.event[i]; |
| 92 |
> |
for (int i = 0; i < pythia->event.size(); ++i){ |
| 93 |
> |
Particle& theParticle = pythia->event[i]; |
| 94 |
|
//if(abs(theParticle.id())==15){cout<<theParticle.daughter1()<<" tau daughters " << theParticle.daughter2()<<endl;} |
| 95 |
|
|
| 96 |
|
// if((theParticle.id())==15 && theParticle.isFinal()){ countertau_fin++;} |
| 124 |
|
if(ufs->bJetList().size()>1 && (ufs->genTauList().size()>1 && countermu_fin>0) )//basic cuts |
| 125 |
|
{ |
| 126 |
|
//the first tree before the b tagging requirements |
| 127 |
< |
for (int i = 0; i < pythia.event.size(); ++i){ |
| 128 |
< |
Particle& theParticle = pythia.event[i]; |
| 127 |
> |
for (int i = 0; i < pythia->event.size(); ++i){ |
| 128 |
> |
Particle& theParticle = pythia->event[i]; |
| 129 |
|
// if(theParticle.isFinal()){ //displays only final state particles |
| 130 |
|
_np++; |
| 131 |
|
_px=theParticle.px(); |
| 177 |
|
higgs_inv_mass->Fill(Einv*Einv - Pxinv*Pxinv - Pyinv*Pyinv - Pzinv*Pzinv); //fills the inv mass for the higgs |
| 178 |
|
} |
| 179 |
|
} |
| 187 |
– |
if(muon_counter>0 && hadron_counter>0){cout<<"no of ev. w/ muon and hadron found = "<<event_counter_passed++<<" out of "<<iEvent<<endl;} |
| 188 |
– |
|
| 180 |
|
|
| 181 |
|
//-----------------NEXT DO FOR ADDITIONAL CUTS and repeat ------------- |
| 182 |
|
for(int i=0; i<mh_size; i++){ |
| 190 |
|
|
| 191 |
|
|
| 192 |
|
//checks the number of bjets with pt over cut |
| 193 |
< |
for(int bi =0; bi<ufs->bJetList().size(); bi++){ |
| 194 |
< |
if(abs(sortedbJets[bi].perp())>bjet_ptcut) bjet_counter25++; //this counts the bjets that pass cutt |
| 195 |
< |
bjet_pt->Fill(abs(sortedbJets[bi].perp())); |
| 193 |
> |
for(unsigned int bi =0; bi<ufs->bJetList().size(); bi++){ |
| 194 |
> |
if(abs(sortedbJets[bi].Et())>bjet_ptcut) bjet_counter25++; //this counts the bjets that pass cutt |
| 195 |
> |
bjet_pt->Fill(abs(sortedbJets[bi].Et())); |
| 196 |
|
} |
| 197 |
|
|
| 198 |
|
//implement the cut on event of btagging |
| 199 |
|
if(bjet_counter25>1){ |
| 200 |
|
if(verbose_) cout<<" event passed btagging selection!!!"<<endl; |
| 201 |
< |
for(int p=0; p< pythia.event.size(); p++){ |
| 202 |
< |
Particle& theParticle=pythia.event[p]; |
| 201 |
> |
for(int p=0; p< pythia->event.size(); p++){ |
| 202 |
> |
Particle& theParticle=pythia->event[p]; |
| 203 |
|
//if(theParticle.status()>0){ |
| 204 |
|
_npc++; |
| 205 |
|
_pxc=theParticle.px(); |
| 248 |
|
} |
| 249 |
|
} |
| 250 |
|
|
| 251 |
< |
|
| 252 |
< |
//ml end |
| 253 |
< |
if(verbose_){ |
| 254 |
< |
cout << "Number of Particles = " << ufs->particleList().size() << endl; |
| 255 |
< |
cout << "Number of Gen Elecs = " << ufs->electronList().size() << endl; |
| 251 |
> |
if(verbose_) { |
| 252 |
> |
cout << "Number of Gen Chrgd = " << ufs->chargedHadronList().size() << endl; |
| 253 |
> |
cout << "Number of Gen Neutr = " << ufs->neutralHadronList().size() << endl; |
| 254 |
> |
cout << "Number of Gen Phtns = " << ufs->photonList().size() << endl; |
| 255 |
> |
cout << "Number of Gen Elcns = " << ufs->electronList().size() << endl; |
| 256 |
|
cout << "Number of Gen Muons = " << ufs->muonList().size() << endl; |
| 257 |
|
cout << "Number of Gen Taus = " << ufs->genTauList().size() << endl; |
| 258 |
+ |
cout << "Number of taus = " << ufs->tauList().size() << endl; |
| 259 |
+ |
cout << "Number of c Quarks = " << ufs->cQuarkList().size() << endl; |
| 260 |
|
cout << "Number of b Quarks = " << ufs->bQuarkList().size() << endl; |
| 261 |
|
cout << "Number of Jets = " << ufs->jetList().size() << endl; |
| 262 |
|
cout << "Number of bJets = " << ufs->bJetList().size() << endl; |
| 270 |
– |
cout << "Number of taus = " << ufs->tauList().size() << endl; |
| 263 |
|
cout << endl; |
| 264 |
|
} |
| 265 |
< |
*/ |
| 265 |
> |
|
| 266 |
> |
iEvent++; |
| 267 |
|
return true; |
| 268 |
|
} |
| 269 |
|
|
