| 1 |
#include "Analysis.h"
|
| 2 |
|
| 3 |
#include <iostream>
|
| 4 |
#include <vector>
|
| 5 |
#include <string>
|
| 6 |
using namespace std;
|
| 7 |
|
| 8 |
#include "Pythia.h"
|
| 9 |
using namespace Pythia8;
|
| 10 |
|
| 11 |
#include "UltraFastSim.h"
|
| 12 |
|
| 13 |
#include "TROOT.h"
|
| 14 |
#include "TFile.h"
|
| 15 |
#include "TH1F.h"
|
| 16 |
#include "TTree.h"
|
| 17 |
|
| 18 |
TROOT root("Mike's ROOT Session","PYTHIA Histograms");
|
| 19 |
|
| 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");
|
| 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");
|
| 27 |
tree->Branch("_py",&_py,"_py/D");
|
| 28 |
tree->Branch("_pz",&_pz,"_pz/D");
|
| 29 |
tree->Branch("_e",&_e,"_e/D");
|
| 30 |
tree->Branch("_m",&_m,"_m/D");
|
| 31 |
tree->Branch("_id",&_id,"_id/I");
|
| 32 |
tree->Branch("_np",&_np,"_np/I");
|
| 33 |
tree->Branch("_mother1",&_mother1,"_mother1/I");
|
| 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 |
|
| 37 |
cuttree = new TTree("T_btagged","after b-tagging");
|
| 38 |
cuttree->Branch("_px",&_pxc,"_pxc/D");
|
| 39 |
cuttree->Branch("_py",&_pyc,"_pyc/D");
|
| 40 |
cuttree->Branch("_pz",&_pzc,"_pzc/D");
|
| 41 |
cuttree->Branch("_e",&_ec,"_ec/D");
|
| 42 |
cuttree->Branch("_m",&_mc,"_mc/D");
|
| 43 |
cuttree->Branch("_id",&_idc,"_idc/I");
|
| 44 |
cuttree->Branch("_np",&_npc,"_npc/I");
|
| 45 |
cuttree->Branch("_mother1",&_mother1c,"_mother2c/I");
|
| 46 |
cuttree->Branch("_mother2",&_mother1c,"_mother2c/I");
|
| 47 |
higgs_inv_mass_cut = new TH1D("higgs_inv_mcut","Higgs Invariant Mass After Cut [GeV]", 10000, -0.5, 999.5);
|
| 48 |
}
|
| 49 |
|
| 50 |
Analysis::~Analysis() {
|
| 51 |
// Save output
|
| 52 |
bjet_pt->Write();
|
| 53 |
higgs_inv_mass->Write();
|
| 54 |
higgs_inv_mass_cut->Write();
|
| 55 |
tree->Write();
|
| 56 |
cuttree->Write();
|
| 57 |
outFile->Close();
|
| 58 |
}
|
| 59 |
|
| 60 |
bool Analysis::run() {
|
| 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...
|
| 66 |
double hadron_store[mh_size][4]; //obviously not a good solution... to set this..
|
| 67 |
double muon_store[mh_size][4];
|
| 68 |
|
| 69 |
int tau_daughters[mh_size][2];
|
| 70 |
|
| 71 |
double Einv=0.;
|
| 72 |
double Pxinv=0.;
|
| 73 |
double Pyinv=0.;
|
| 74 |
double Pzinv=0.;
|
| 75 |
|
| 76 |
for(int i=0; i<mh_size; i++){
|
| 77 |
for(int j=0; j<4; j++){
|
| 78 |
muon_store[i][j]=0.;
|
| 79 |
hadron_store[i][j]=0.;
|
| 80 |
}
|
| 81 |
for(int j=0; j<2;j++){
|
| 82 |
tau_daughters[i][j]=0;
|
| 83 |
}
|
| 84 |
|
| 85 |
}
|
| 86 |
_np=0;
|
| 87 |
_npc=0;
|
| 88 |
|
| 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];
|
| 94 |
//if(abs(theParticle.id())==15){cout<<theParticle.daughter1()<<" tau daughters " << theParticle.daughter2()<<endl;}
|
| 95 |
|
| 96 |
// if((theParticle.id())==15 && theParticle.isFinal()){ countertau_fin++;}
|
| 97 |
// else if(abs(theParticle.id())==15 && !theParticle.isFinal() ){countertau_nf++; }
|
| 98 |
if(abs(theParticle.id())==15){
|
| 99 |
if(verbose_) cout<<theParticle.daughter1()<<" tau daughters " << theParticle.daughter2()<<endl;
|
| 100 |
tau_daughters[tau_counter][0]=theParticle.daughter1();
|
| 101 |
tau_daughters[tau_counter][1]=theParticle.daughter2();
|
| 102 |
tau_counter++;
|
| 103 |
}
|
| 104 |
|
| 105 |
if((theParticle.id())==13 && theParticle.isFinal()){countermu_fin++;}
|
| 106 |
else if (abs(theParticle.id())==13 && !theParticle.isFinal()){countermu_nf++;}
|
| 107 |
|
| 108 |
//if(abs(theParticle.id())==13){cout<<theParticle.mother1()<<" mu mothers "<<theParticle.mother2()<<endl;}
|
| 109 |
}
|
| 110 |
|
| 111 |
if(verbose_){
|
| 112 |
cout<<" tauList.size()= "<<ufs->tauList().size()
|
| 113 |
<<" genTauList.size()= "<<ufs->genTauList().size()
|
| 114 |
<<" muList.size()= "<<ufs->muonList().size()
|
| 115 |
<<" countertau_fin "<< countertau_fin
|
| 116 |
<<" countertau_nf "<<countertau_nf
|
| 117 |
<<" countermu_fin "<<countermu_fin
|
| 118 |
<<" countermu_nf "<<countermu_nf
|
| 119 |
<<endl;
|
| 120 |
}
|
| 121 |
|
| 122 |
|
| 123 |
// int indextau(0);
|
| 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];
|
| 129 |
// if(theParticle.isFinal()){ //displays only final state particles
|
| 130 |
_np++;
|
| 131 |
_px=theParticle.px();
|
| 132 |
_py=theParticle.py();
|
| 133 |
_pz=theParticle.pz();
|
| 134 |
_e=theParticle.e();
|
| 135 |
_m=theParticle.m();
|
| 136 |
_id=theParticle.id();
|
| 137 |
_mother1=theParticle.mother1();
|
| 138 |
_mother2=theParticle.mother2();
|
| 139 |
tree->Fill();
|
| 140 |
|
| 141 |
//this stores the muons and hadrons 4mom -- checks if mother of particle is tau+- daughter by looping over taus
|
| 142 |
for(int i=0; i<tau_counter; i++){
|
| 143 |
//checks whether the particle is a muon or not
|
| 144 |
if( theParticle.mother1() == tau_daughters[i][0] || theParticle.mother1()==tau_daughters[i][1] || (theParticle.mother1()>tau_daughters[i][0] && theParticle.mother1()< tau_daughters[i][1]))
|
| 145 |
{
|
| 146 |
if(abs(theParticle.id())==13){ //if the particle is muon+-
|
| 147 |
if(verbose_)cout<<theParticle.mother1()<<" muon mothers "<<theParticle.mother2()<<endl;
|
| 148 |
muon_store[muon_counter][1]=theParticle.px();
|
| 149 |
muon_store[muon_counter][2]=theParticle.py();
|
| 150 |
muon_store[muon_counter][3]=theParticle.pz();
|
| 151 |
muon_store[muon_counter][0]=theParticle.e();
|
| 152 |
//cout<<muon_store[muon_counter][1]<<endl; //some issue with storing the data
|
| 153 |
muon_counter++;
|
| 154 |
}
|
| 155 |
else if( abs(theParticle.id())>18 && (abs(theParticle.id())<21 || abs(theParticle.id())>42) ){ //if the particle is neither lepton or sm gauge boson or single quark
|
| 156 |
// cout<<theParticle.mother1()<<" hadron mothers "<<theParticle.mother2()<<endl;
|
| 157 |
hadron_store[hadron_counter][1]=theParticle.px();
|
| 158 |
hadron_store[hadron_counter][2]=theParticle.py();
|
| 159 |
hadron_store[hadron_counter][3]=theParticle.pz();
|
| 160 |
hadron_store[hadron_counter][0]=theParticle.e();
|
| 161 |
hadron_counter++;
|
| 162 |
}
|
| 163 |
}
|
| 164 |
}//end loop over tau
|
| 165 |
// }//end loop over final particles only
|
| 166 |
}//end loop over particles in event
|
| 167 |
|
| 168 |
if(verbose_) cout<<"Muon counter "<<muon_counter<<" had counter "<<hadron_counter<<endl;
|
| 169 |
//find invariant mass for all combos
|
| 170 |
for(int i=0; i<muon_counter; i++){
|
| 171 |
for(int j=0; j<hadron_counter; j++){
|
| 172 |
Einv = muon_store[i][0]+hadron_store[j][0];
|
| 173 |
Pxinv= muon_store[i][1]+hadron_store[j][1];
|
| 174 |
Pyinv= muon_store[i][2]+hadron_store[j][2];
|
| 175 |
Pzinv= muon_store[i][3]+hadron_store[j][3];
|
| 176 |
//cout<<" inv mass "<<Einv*Einv - Pxinv*Pxinv - Pyinv*Pyinv - Pzinv*Pzinv<<endl;
|
| 177 |
higgs_inv_mass->Fill(Einv*Einv - Pxinv*Pxinv - Pyinv*Pyinv - Pzinv*Pzinv); //fills the inv mass for the higgs
|
| 178 |
}
|
| 179 |
}
|
| 180 |
|
| 181 |
//-----------------NEXT DO FOR ADDITIONAL CUTS and repeat -------------
|
| 182 |
for(int i=0; i<mh_size; i++){
|
| 183 |
for(int j=0; j<4; j++){
|
| 184 |
muon_store[i][j]=0.;
|
| 185 |
hadron_store[i][j]=0.;
|
| 186 |
}
|
| 187 |
}
|
| 188 |
hadron_counter=0;
|
| 189 |
muon_counter=0;
|
| 190 |
|
| 191 |
|
| 192 |
//checks the number of bjets with pt over cut
|
| 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];
|
| 203 |
//if(theParticle.status()>0){
|
| 204 |
_npc++;
|
| 205 |
_pxc=theParticle.px();
|
| 206 |
_pyc=theParticle.py();
|
| 207 |
_pzc=theParticle.pz();
|
| 208 |
_ec =theParticle.e();
|
| 209 |
_mc = theParticle.m();
|
| 210 |
_idc=theParticle.id();
|
| 211 |
_mother1c=theParticle.mother1();
|
| 212 |
_mother2c=theParticle.mother2();
|
| 213 |
cuttree->Fill();
|
| 214 |
|
| 215 |
for(int i=0; i<tau_counter; i++){
|
| 216 |
if( theParticle.mother1() == tau_daughters[i][0] || theParticle.mother1()==tau_daughters[i][1] || (theParticle.mother1()>tau_daughters[i][0] && theParticle.mother1()< tau_daughters[i][1]))
|
| 217 |
{
|
| 218 |
|
| 219 |
if(abs(theParticle.id())==13){ //if the particle is muon+-
|
| 220 |
muon_store[muon_counter][1]=theParticle.px();
|
| 221 |
muon_store[muon_counter][2]=theParticle.py();
|
| 222 |
muon_store[muon_counter][3]=theParticle.pz();
|
| 223 |
muon_store[muon_counter][0]=theParticle.e();
|
| 224 |
muon_counter++;
|
| 225 |
}
|
| 226 |
else if( (abs(theParticle.id())<11 || abs(theParticle.id())>16) && (abs(theParticle.id())<21 || abs(theParticle.id())>25)){ //if the particle is neither lepton or sm gauge boson
|
| 227 |
hadron_store[hadron_counter][1]=theParticle.px();
|
| 228 |
hadron_store[hadron_counter][2]=theParticle.py();
|
| 229 |
hadron_store[hadron_counter][3]=theParticle.pz();
|
| 230 |
hadron_store[hadron_counter][0]=theParticle.e();
|
| 231 |
hadron_counter++;
|
| 232 |
}
|
| 233 |
}
|
| 234 |
}//loop over taus
|
| 235 |
//} //loop over only final state particles
|
| 236 |
}//loop over particles in one event
|
| 237 |
}//bjet cut loop
|
| 238 |
}//basic cut loop
|
| 239 |
|
| 240 |
//find invariant mass for all combos after cuts
|
| 241 |
for(int i=0; i<muon_counter; i++){
|
| 242 |
for(int j=0; j<hadron_counter; j++){
|
| 243 |
Einv = muon_store[i][0]+hadron_store[j][0];
|
| 244 |
Pxinv= muon_store[i][1]+hadron_store[j][1];
|
| 245 |
Pyinv= muon_store[i][2]+hadron_store[j][2];
|
| 246 |
Pzinv= muon_store[i][3]+hadron_store[j][3];
|
| 247 |
higgs_inv_mass_cut->Fill(Einv*Einv - Pxinv*Pxinv - Pyinv*Pyinv - Pzinv*Pzinv); //fills the inv mass for the higgs
|
| 248 |
}
|
| 249 |
}
|
| 250 |
|
| 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;
|
| 263 |
cout << endl;
|
| 264 |
}
|
| 265 |
|
| 266 |
iEvent++;
|
| 267 |
return true;
|
| 268 |
}
|
| 269 |
|
