diff --git a/Calorimetry/include/JetAnalyzer.h b/Calorimetry/include/JetAnalyzer.h index 28ac8b7..0167001 100644 --- a/Calorimetry/include/JetAnalyzer.h +++ b/Calorimetry/include/JetAnalyzer.h @@ -54,12 +54,13 @@ class JetAnalyzer : public Processor { protected: std::string m_inputMCParticleCollection=""; - std::string m_inputRECOParticleCollection=""; - std::string m_genjetColName=""; - std::string m_recojetColName=""; +// std::string m_inputRECOParticleCollection=""; +// std::string m_genjetColName=""; +// std::string m_recojetColName=""; std::string m_rootFileName=""; + std::string m_processName=""; // Run and event counters int m_eventNumber=0; @@ -78,6 +79,70 @@ class JetAnalyzer : public Processor { float m_d2_mcPz=0; //particle level information sums --> visible and invisible energy +//#######NEW##################### + int m_H1_mcPDGID = 0; + float m_H1_mcE = 0; + float m_H1_mcPx = 0; + float m_H1_mcPy = 0; + float m_H1_mcPz = 0; + + int m_H2_mcPDGID = 0; + float m_H2_mcE = 0; + float m_H2_mcPx = 0; + float m_H2_mcPy = 0; + float m_H2_mcPz = 0; + + int m_b1_mcPDGID = 0; + float m_b1_mcE = 0; + float m_b1_mcPx = 0; + float m_b1_mcPy = 0; + float m_b1_mcPz = 0; + + int m_b2_mcPDGID = 0; + float m_b2_mcE = 0; + float m_b2_mcPx = 0; + float m_b2_mcPy = 0; + float m_b2_mcPz = 0; + + int m_b3_mcPDGID = 0; + float m_b3_mcE = 0; + float m_b3_mcPx = 0; + float m_b3_mcPy = 0; + float m_b3_mcPz = 0; + + int m_b4_mcPDGID = 0; + float m_b4_mcE = 0; + float m_b4_mcPx = 0; + float m_b4_mcPy = 0; + float m_b4_mcPz = 0; + + int m_b5_mcPDGID = 0; + float m_b5_mcE = 0; + float m_b5_mcPx = 0; + float m_b5_mcPy = 0; + float m_b5_mcPz = 0; + + int m_b6_mcPDGID = 0; + float m_b6_mcE = 0; + float m_b6_mcPx = 0; + float m_b6_mcPy = 0; + float m_b6_mcPz = 0; + + int m_b7_mcPDGID = 0; + float m_b7_mcE = 0; + float m_b7_mcPx = 0; + float m_b7_mcPy = 0; + float m_b7_mcPz = 0; + + int m_b8_mcPDGID = 0; + float m_b8_mcE = 0; + float m_b8_mcPx = 0; + float m_b8_mcPy = 0; + float m_b8_mcPz = 0; + + +//#######ENDNEW##################### + //only neutrinos float m_E_trueInv=0; float m_px_trueInv=0; @@ -90,12 +155,12 @@ class JetAnalyzer : public Processor { float m_py_trueVis=0; float m_pz_trueVis=0; - //all reconstructed particles +/* //all reconstructed particles float m_E_totPFO=0; float m_px_totPFO=0; float m_py_totPFO=0; float m_pz_totPFO=0; - +*/ //if extended parton quantities are saved: dibosons and their decay products std::vector* m_trueME_E=NULL; std::vector* m_trueME_Px=NULL; @@ -103,7 +168,7 @@ class JetAnalyzer : public Processor { std::vector* m_trueME_Pz=NULL; std::vector* m_trueME_PDGID=NULL; - std::vector* m_genJet_E=NULL; +/* std::vector* m_genJet_E=NULL; std::vector* m_genJet_Px=NULL; std::vector* m_genJet_Py=NULL; std::vector* m_genJet_Pz=NULL; @@ -112,7 +177,7 @@ class JetAnalyzer : public Processor { std::vector* m_recoJet_Px=NULL; std::vector* m_recoJet_Py=NULL; std::vector* m_recoJet_Pz=NULL; - +*/ // Call to get collections void getCollection(LCCollection*&, std::string, LCEvent*); diff --git a/Calorimetry/src/JetAnalyzer.cc b/Calorimetry/src/JetAnalyzer.cc index 6a3835d..baa2065 100644 --- a/Calorimetry/src/JetAnalyzer.cc +++ b/Calorimetry/src/JetAnalyzer.cc @@ -38,24 +38,29 @@ JetAnalyzer::JetAnalyzer() : Processor("JetAnalyzer") { m_rootFileName, std::string("JetAnalyzer.root")); - registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE, + registerProcessorParameter( "ProcessName", + "Which process to run", + m_processName, + std::string("ProcessName")); + +/* registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE, "RECOParticleCollectionName", "Name of the RECOParticle input collection", m_inputRECOParticleCollection, std::string("PandoraPFOs")); - - registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE, +*/ +/* registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE, "GenJetCollection" , "Name of the GenJet collection", m_genjetColName, std::string("GenJet_VLC")); - registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE, + registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE, "RecoJetCollection" , "Name of the RecoJet collection", m_recojetColName, std::string("RecoJet_VLC")); - +*/ registerProcessorParameter("doDiBosonChecks", "MC truth diboson checks for on-shell bosons", m_performDiBosonChecks, @@ -92,7 +97,7 @@ void JetAnalyzer::init() { m_trueME_PDGID=new std::vector(); } - m_genJet_E = new std::vector(); +/* m_genJet_E = new std::vector(); m_genJet_Px = new std::vector(); m_genJet_Py = new std::vector(); m_genJet_Pz = new std::vector(); @@ -101,7 +106,7 @@ void JetAnalyzer::init() { m_recoJet_Px = new std::vector(); m_recoJet_Py = new std::vector(); m_recoJet_Pz = new std::vector(); - +*/ m_E_trueInv = 0; m_px_trueInv = 0; m_py_trueInv = 0; @@ -112,11 +117,11 @@ void JetAnalyzer::init() { m_py_trueVis = 0; m_pz_trueVis = 0; - m_E_totPFO = 0; + /* m_E_totPFO = 0; m_px_totPFO = 0; m_py_totPFO = 0; m_pz_totPFO = 0; - +*/ m_d1_mcPDGID = 0; m_d1_mcE = 0; m_d1_mcPx = 0; @@ -129,7 +134,68 @@ void JetAnalyzer::init() { m_d2_mcPy = 0; m_d2_mcPz = 0; - m_genJet_E ->clear(); + m_H1_mcPDGID = 0; + m_H1_mcE = 0; + m_H1_mcPx = 0; + m_H1_mcPy = 0; + m_H1_mcPz = 0; + + m_H2_mcPDGID = 0; + m_H2_mcE = 0; + m_H2_mcPx = 0; + m_H2_mcPy = 0; + m_H2_mcPz = 0; + + m_b1_mcPDGID = 0; + m_b1_mcE = 0; + m_b1_mcPx = 0; + m_b1_mcPy = 0; + m_b1_mcPz = 0; + + m_b2_mcPDGID = 0; + m_b2_mcE = 0; + m_b2_mcPx = 0; + m_b2_mcPy = 0; + m_b2_mcPz = 0; + + m_b3_mcPDGID = 0; + m_b3_mcE = 0; + m_b3_mcPx = 0; + m_b3_mcPy = 0; + m_b3_mcPz = 0; + + m_b4_mcPDGID = 0; + m_b4_mcE = 0; + m_b4_mcPx = 0; + m_b4_mcPy = 0; + m_b4_mcPz = 0; + + m_b5_mcPDGID = 0; + m_b5_mcE = 0; + m_b5_mcPx = 0; + m_b5_mcPy = 0; + m_b5_mcPz = 0; + + m_b6_mcPDGID = 0; + m_b6_mcE = 0; + m_b6_mcPx = 0; + m_b6_mcPy = 0; + m_b6_mcPz = 0; + + m_b7_mcPDGID = 0; + m_b7_mcE = 0; + m_b7_mcPx = 0; + m_b7_mcPy = 0; + m_b7_mcPz = 0; + + m_b8_mcPDGID = 0; + m_b8_mcE = 0; + m_b8_mcPx = 0; + m_b8_mcPy = 0; + m_b8_mcPz = 0; + +//####################################### +/* m_genJet_E ->clear(); m_genJet_Px ->clear(); m_genJet_Py ->clear(); m_genJet_Pz ->clear(); @@ -138,7 +204,7 @@ void JetAnalyzer::init() { m_recoJet_Px ->clear(); m_recoJet_Py ->clear(); m_recoJet_Pz ->clear(); - +*/ if(m_fillMEInfo){ m_trueME_E->clear(); m_trueME_Px->clear(); @@ -167,6 +233,69 @@ void JetAnalyzer::init() { m_outputTree->Branch("d2_mcPy",&m_d2_mcPy,"d2_mcPy/F"); m_outputTree->Branch("d2_mcPz",&m_d2_mcPz,"d2_mcPz/F"); + //########################NEW FOR DOUBLE HIGGS ##########à + m_outputTree->Branch("H1_mcPDGID",&m_H1_mcPDGID,"H1_mcPDGID/I"); + m_outputTree->Branch("H1_mcE",&m_H1_mcE,"H1_mcE/F"); + m_outputTree->Branch("H1_mcPx",&m_H1_mcPx,"H1_mcPx/F"); + m_outputTree->Branch("H1_mcPy",&m_H1_mcPy,"H1_mcPy/F"); + m_outputTree->Branch("H1_mcPz",&m_H1_mcPz,"H1_mcPz/F"); + + m_outputTree->Branch("H2_mcPDGID",&m_H2_mcPDGID,"H2_mcPDGID/I"); + m_outputTree->Branch("H2_mcE",&m_H2_mcE,"H2_mcE/F"); + m_outputTree->Branch("H2_mcPx",&m_H2_mcPx,"H2_mcPx/F"); + m_outputTree->Branch("H2_mcPy",&m_H2_mcPy,"H2_mcPy/F"); + m_outputTree->Branch("H2_mcPz",&m_H2_mcPz,"H2_mcPz/F"); + + m_outputTree->Branch("b1_mcPDGID",&m_b1_mcPDGID,"b1_mcPDGID/I"); + m_outputTree->Branch("b1_mcE",&m_b1_mcE,"b1_mcE/F"); + m_outputTree->Branch("b1_mcPx",&m_b1_mcPx,"b1_mcPx/F"); + m_outputTree->Branch("b1_mcPy",&m_b1_mcPy,"b1_mcPy/F"); + m_outputTree->Branch("b1_mcPz",&m_b1_mcPz,"b1_mcPz/F"); + + m_outputTree->Branch("b2_mcPDGID",&m_b2_mcPDGID,"b2_mcPDGID/I"); + m_outputTree->Branch("b2_mcE",&m_b2_mcE,"b2_mcE/F"); + m_outputTree->Branch("b2_mcPx",&m_b2_mcPx,"b2_mcPx/F"); + m_outputTree->Branch("b2_mcPy",&m_b2_mcPy,"b2_mcPy/F"); + m_outputTree->Branch("b2_mcPz",&m_b2_mcPz,"b2_mcPz/F"); + + m_outputTree->Branch("b3_mcPDGID",&m_b3_mcPDGID,"b3_mcPDGID/I"); + m_outputTree->Branch("b3_mcE",&m_b3_mcE,"b3_mcE/F"); + m_outputTree->Branch("b3_mcPx",&m_b3_mcPx,"b3_mcPx/F"); + m_outputTree->Branch("b3_mcPy",&m_b3_mcPy,"b3_mcPy/F"); + m_outputTree->Branch("b3_mcPz",&m_b3_mcPz,"b3_mcPz/F"); + + m_outputTree->Branch("b4_mcPDGID",&m_b4_mcPDGID,"b4_mcPDGID/I"); + m_outputTree->Branch("b4_mcE",&m_b4_mcE,"b4_mcE/F"); + m_outputTree->Branch("b4_mcPx",&m_b4_mcPx,"b4_mcPx/F"); + m_outputTree->Branch("b4_mcPy",&m_b4_mcPy,"b4_mcPy/F"); + m_outputTree->Branch("b4_mcPz",&m_b4_mcPz,"b4_mcPz/F"); + + m_outputTree->Branch("b5_mcPDGID",&m_b5_mcPDGID,"b5_mcPDGID/I"); + m_outputTree->Branch("b5_mcE",&m_b5_mcE,"b5_mcE/F"); + m_outputTree->Branch("b5_mcPx",&m_b5_mcPx,"b5_mcPx/F"); + m_outputTree->Branch("b5_mcPy",&m_b5_mcPy,"b5_mcPy/F"); + m_outputTree->Branch("b5_mcPz",&m_b5_mcPz,"b5_mcPz/F"); + + m_outputTree->Branch("b6_mcPDGID",&m_b6_mcPDGID,"b6_mcPDGID/I"); + m_outputTree->Branch("b6_mcE",&m_b6_mcE,"b6_mcE/F"); + m_outputTree->Branch("b6_mcPx",&m_b6_mcPx,"b6_mcPx/F"); + m_outputTree->Branch("b6_mcPy",&m_b6_mcPy,"b6_mcPy/F"); + m_outputTree->Branch("b6_mcPz",&m_b6_mcPz,"b6_mcPz/F"); + + m_outputTree->Branch("b7_mcPDGID",&m_b7_mcPDGID,"b7_mcPDGID/I"); + m_outputTree->Branch("b7_mcE",&m_b7_mcE,"b7_mcE/F"); + m_outputTree->Branch("b7_mcPx",&m_b7_mcPx,"b7_mcPx/F"); + m_outputTree->Branch("b7_mcPy",&m_b7_mcPy,"b7_mcPy/F"); + m_outputTree->Branch("b7_mcPz",&m_b7_mcPz,"b7_mcPz/F"); + + m_outputTree->Branch("b8_mcPDGID",&m_b8_mcPDGID,"b8_mcPDGID/I"); + m_outputTree->Branch("b8_mcE",&m_b8_mcE,"b8_mcE/F"); + m_outputTree->Branch("b8_mcPx",&m_b8_mcPx,"b8_mcPx/F"); + m_outputTree->Branch("b8_mcPy",&m_b8_mcPy,"b8_mcPy/F"); + m_outputTree->Branch("b8_mcPz",&m_b8_mcPz,"b8_mcPz/F"); + +//###############################ENDNEW + //true particle level, exclude neutrinos, true visible content m_outputTree->Branch("E_trueVis" ,&m_E_trueVis, "E_trueVis/F"); m_outputTree->Branch("Px_trueVis",&m_px_trueVis,"Px_trueVis/F"); @@ -180,7 +309,7 @@ void JetAnalyzer::init() { m_outputTree->Branch("Pz_trueInv",&m_pz_trueInv,"Pz_trueInv/F"); //reconstructed leve - m_outputTree->Branch("E_totPFO" ,&m_E_totPFO, "E_totPFO/F"); +/* m_outputTree->Branch("E_totPFO" ,&m_E_totPFO, "E_totPFO/F"); m_outputTree->Branch("Px_totPFO",&m_px_totPFO,"Px_totPFO/F"); m_outputTree->Branch("Py_totPFO",&m_py_totPFO,"Py_totPFO/F"); m_outputTree->Branch("Pz_totPFO",&m_pz_totPFO,"Pz_totPFO/F"); @@ -195,7 +324,7 @@ void JetAnalyzer::init() { m_outputTree->Branch("recoJetPx", "std::vector< float >", &m_recoJet_Px); m_outputTree->Branch("recoJetPy", "std::vector< float >", &m_recoJet_Py); m_outputTree->Branch("recoJetPz", "std::vector< float >", &m_recoJet_Pz); - +*/ } @@ -219,11 +348,11 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_py_trueVis = 0; m_pz_trueVis = 0; - m_E_totPFO = 0; +/* m_E_totPFO = 0; m_px_totPFO = 0; m_py_totPFO = 0; m_pz_totPFO = 0; - +*/ m_d1_mcPDGID=-10; m_d1_mcE=-10; m_d1_mcPx=-10; @@ -236,7 +365,70 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_d2_mcPy=-10; m_d2_mcPz=-10; - m_genJet_E ->clear(); +//#######NEW##################### + m_H1_mcPDGID = -10; + m_H1_mcE = -10; + m_H1_mcPx = -10; + m_H1_mcPy = -10; + m_H1_mcPz = -10; + + m_H2_mcPDGID = -10; + m_H2_mcE = -10; + m_H2_mcPx = -10; + m_H2_mcPy = -10; + m_H2_mcPz = -10; + + + m_b1_mcPDGID = -10; + m_b1_mcE = -10; + m_b1_mcPx = -10; + m_b1_mcPy = -10; + m_b1_mcPz = -10; + + m_b2_mcPDGID = -10; + m_b2_mcE = -10; + m_b2_mcPx = -10; + m_b2_mcPy = -10; + m_b2_mcPz = -10; + + + m_b3_mcPDGID = -10; + m_b3_mcE = -10; + m_b3_mcPx = -10; + m_b3_mcPy = -10; + m_b3_mcPz = -10; + + m_b4_mcPDGID = -10; + m_b4_mcE = -10; + m_b4_mcPx = -10; + m_b4_mcPy = -10; + m_b4_mcPz = -10; + + m_b5_mcPDGID = -10; + m_b5_mcE = -10; + m_b5_mcPx = -10; + m_b5_mcPy = -10; + m_b5_mcPz = -10; + + m_b6_mcPDGID = -10; + m_b6_mcE = -10; + m_b6_mcPx = -10; + m_b6_mcPy = -10; + m_b6_mcPz = -10; + + m_b7_mcPDGID = -10; + m_b7_mcE = -10; + m_b7_mcPx = -10; + m_b7_mcPy = -10; + m_b7_mcPz = -10; + + m_b8_mcPDGID = -10; + m_b8_mcE = -10; + m_b8_mcPx = -10; + m_b8_mcPy = -10; + m_b8_mcPz = -10; +//#######ENDNEW##################### +/* m_genJet_E ->clear(); m_genJet_Px ->clear(); m_genJet_Py ->clear(); m_genJet_Pz ->clear(); @@ -245,7 +437,7 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_recoJet_Px ->clear(); m_recoJet_Py ->clear(); m_recoJet_Pz ->clear(); - +*/ if(m_fillMEInfo){ m_trueME_E->clear(); @@ -261,9 +453,19 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { bool pass_W_boson_mass=true; bool pass_Z_boson_mass=true; + if(mcColl!=NULL){ int boson_counter=0; int ind_second_boson=-1; + MCParticle* index_H1=NULL; //NEW + MCParticle* index_H2=NULL; + MCParticle* index_b1=NULL; + MCParticle* index_b2=NULL; + MCParticle* index_b3=NULL; + MCParticle* index_b4=NULL; + MCParticle* index_b5=NULL; + MCParticle* index_b7=NULL; + bool first_filled=true; for(int m =0; m< mcColl->getNumberOfElements(); m++){ MCParticle* mcp= dynamic_cast( mcColl->getElementAt(m) ) ; //boson checks have to be done for WW,WZ and ZZ configurations @@ -301,21 +503,366 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_trueME_PDGID->push_back(mcp->getPDG()); } //fill first quark anti quark pair which appears in the history - //this IS the hadronic W and Z (in case of a mixed di-boson event) - if(abs(mcp->getPDG())<7 && m_d1_mcE<0) { - m_d1_mcPDGID=mcp->getPDG(); - m_d1_mcE=mcp->getEnergy(); - m_d1_mcPx=mcp->getMomentum()[0]; - m_d1_mcPy=mcp->getMomentum()[1]; - m_d1_mcPz=mcp->getMomentum()[2]; + //this IS the hadronic W and Z (in case of a mixed di-boson event) + if (m_processName == "dijet") { + if(abs(mcp->getPDG())<7 && m_d1_mcE<0) { + m_d1_mcPDGID=mcp->getPDG(); + m_d1_mcE=mcp->getEnergy(); + m_d1_mcPx=mcp->getMomentum()[0]; + m_d1_mcPy=mcp->getMomentum()[1]; + m_d1_mcPz=mcp->getMomentum()[2]; + } + if(m_d2_mcE<0 && (abs(mcp->getPDG())<7 && mcp->getPDG()==(-m_d1_mcPDGID))){ + m_d2_mcPDGID=mcp->getPDG(); + m_d2_mcE=mcp->getEnergy(); + m_d2_mcPx=mcp->getMomentum()[0]; + m_d2_mcPy=mcp->getMomentum()[1]; + m_d2_mcPz=mcp->getMomentum()[2]; + } + } + + if (m_processName == "HHbbbb") { + if(mcp->getPDG()==25 && m_H1_mcE<0) { + index_H1=mcp; + m_H1_mcPDGID=mcp->getPDG(); + m_H1_mcE=mcp->getEnergy(); + m_H1_mcPx=mcp->getMomentum()[0]; + m_H1_mcPy=mcp->getMomentum()[1]; + m_H1_mcPz=mcp->getMomentum()[2]; + } + if(m_H2_mcE<0 && (mcp->getPDG()==25 && index_H1!=mcp )){ + index_H2=mcp; + m_H2_mcPDGID=mcp->getPDG(); + m_H2_mcE=mcp->getEnergy(); + m_H2_mcPx=mcp->getMomentum()[0]; + m_H2_mcPy=mcp->getMomentum()[1]; + m_H2_mcPz=mcp->getMomentum()[2]; + } + if(abs(mcp->getPDG())==5 && m_b1_mcE<0 && mcp->getParents()[0]==index_H1) { + m_b1_mcPDGID=mcp->getPDG(); + m_b1_mcE=mcp->getEnergy(); + m_b1_mcPx=mcp->getMomentum()[0]; + m_b1_mcPy=mcp->getMomentum()[1]; + m_b1_mcPz=mcp->getMomentum()[2]; + } + if(m_b2_mcE<0 && (abs(mcp->getPDG())==5 && mcp->getPDG()==(-m_b1_mcPDGID)) && mcp->getParents()[0]==index_H1 ){ + m_b2_mcPDGID=mcp->getPDG(); + m_b2_mcE=mcp->getEnergy(); + m_b2_mcPx=mcp->getMomentum()[0]; + m_b2_mcPy=mcp->getMomentum()[1]; + m_b2_mcPz=mcp->getMomentum()[2]; + } + if(abs(mcp->getPDG())==5 && m_b3_mcE<0 && mcp->getParents()[0]==index_H2) { + m_b3_mcPDGID=mcp->getPDG(); + m_b3_mcE=mcp->getEnergy(); + m_b3_mcPx=mcp->getMomentum()[0]; + m_b3_mcPy=mcp->getMomentum()[1]; + m_b3_mcPz=mcp->getMomentum()[2]; + } + if(m_b4_mcE<0 && (abs(mcp->getPDG())==5 && mcp->getPDG()==(-m_b3_mcPDGID))&& mcp->getParents()[0]==index_H2){ + m_b4_mcPDGID=mcp->getPDG(); + m_b4_mcE=mcp->getEnergy(); + m_b4_mcPx=mcp->getMomentum()[0]; + m_b4_mcPy=mcp->getMomentum()[1]; + m_b4_mcPz=mcp->getMomentum()[2]; + } + } + + // if (m_processName == "bb") { + // if(mcp->getPDG()==5 ) { + // m_b1_mcPDGID=mcp->getPDG(); + // m_b1_mcE=mcp->getEnergy(); + // m_b1_mcPx=mcp->getMomentum()[0]; + // m_b1_mcPy=mcp->getMomentum()[1]; + // m_b1_mcPz=mcp->getMomentum()[2]; + // } + // if(mcp->getPDG()==-5 ){ + // m_b2_mcPDGID=mcp->getPDG(); + // m_b2_mcE=mcp->getEnergy(); + // m_b2_mcPx=mcp->getMomentum()[0]; + // m_b2_mcPy=mcp->getMomentum()[1]; + // m_b2_mcPz=mcp->getMomentum()[2]; + // } + // } + + // if (m_processName == "cc") { + // if(mcp->getPDG()==4 ) { + // m_b1_mcPDGID=mcp->getPDG(); + // m_b1_mcE=mcp->getEnergy(); + // m_b1_mcPx=mcp->getMomentum()[0]; + // m_b1_mcPy=mcp->getMomentum()[1]; + // m_b1_mcPz=mcp->getMomentum()[2]; + // } + // if(mcp->getPDG()==-4 ){ + // m_b2_mcPDGID=mcp->getPDG(); + // m_b2_mcE=mcp->getEnergy(); + // m_b2_mcPx=mcp->getMomentum()[0]; + // m_b2_mcPy=mcp->getMomentum()[1]; + // m_b2_mcPz=mcp->getMomentum()[2]; + // } + // } + + // if (m_processName == "qq") { + // if(abs(mcp->getPDG())==3 || abs(mcp->getPDG())==2 || abs(mcp->getPDG())==1) { + // m_b1_mcPDGID=mcp->getPDG(); + // m_b1_mcE=mcp->getEnergy(); + // m_b1_mcPx=mcp->getMomentum()[0]; + // m_b1_mcPy=mcp->getMomentum()[1]; + // m_b1_mcPz=mcp->getMomentum()[2]; + // } + // if(abs(mcp->getPDG())==3 || abs(mcp->getPDG())==2 || abs(mcp->getPDG())==1){ + // m_b2_mcPDGID=mcp->getPDG(); + // m_b2_mcE=mcp->getEnergy(); + // m_b2_mcPx=mcp->getMomentum()[0]; + // m_b2_mcPy=mcp->getMomentum()[1]; + // m_b2_mcPz=mcp->getMomentum()[2]; + // } + // } + + + if (m_processName == "bbbb") { + if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==5000 ) { + index_H1=mcp; + m_H1_mcE=mcp->getEnergy(); + m_H1_mcPx=mcp->getMomentum()[0]; + m_H1_mcPy=mcp->getMomentum()[1]; + m_H1_mcPz=mcp->getMomentum()[2]; + } + if(mcp->getPDG()==-13 && m_H2_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==5000 ) { + index_H2=mcp; + m_H2_mcPDGID=mcp->getPDG(); + m_H2_mcE=mcp->getEnergy(); + m_H2_mcPx=mcp->getMomentum()[0]; + m_H2_mcPy=mcp->getMomentum()[1]; + m_H2_mcPz=mcp->getMomentum()[2]; + } + if(mcp->getPDG()==5 && m_b1_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==5000 && (mcp->getParents()[0]==index_H1 || mcp->getParents()[0]==index_H2) ) { + m_b1_mcPDGID=mcp->getPDG(); + m_b1_mcE=mcp->getEnergy(); + m_b1_mcPx=mcp->getMomentum()[0]; + m_b1_mcPy=mcp->getMomentum()[1]; + m_b1_mcPz=mcp->getMomentum()[2]; + continue; + } + if(mcp->getPDG()==-5 && m_b2_mcE<0 && (mcp->getParents()[0]==index_H1 || mcp->getParents()[0]==index_H2) ) { + m_b2_mcPDGID=mcp->getPDG(); + m_b2_mcE=mcp->getEnergy(); + m_b2_mcPx=mcp->getMomentum()[0]; + m_b2_mcPy=mcp->getMomentum()[1]; + m_b2_mcPz=mcp->getMomentum()[2]; + continue; + } + if(mcp->getPDG()==5 && m_b3_mcE<0&& (mcp->getParents()[0]==index_H1 || mcp->getParents()[0]==index_H2)) { + m_b3_mcPDGID=mcp->getPDG(); + m_b3_mcE=mcp->getEnergy(); + m_b3_mcPx=mcp->getMomentum()[0]; + m_b3_mcPy=mcp->getMomentum()[1]; + m_b3_mcPz=mcp->getMomentum()[2]; + continue; + } + if(mcp->getPDG()==-5 && m_b4_mcE<0&& (mcp->getParents()[0]==index_H1 || mcp->getParents()[0]==index_H2) ) { + m_b4_mcPDGID=mcp->getPDG(); + m_b4_mcE=mcp->getEnergy(); + m_b4_mcPx=mcp->getMomentum()[0]; + m_b4_mcPy=mcp->getMomentum()[1]; + m_b4_mcPz=mcp->getMomentum()[2]; + continue; + } + } + + if (m_processName == "bbH") { + if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==5000 ) { + index_H1=mcp; + m_H1_mcE=mcp->getEnergy(); + m_H1_mcPx=mcp->getMomentum()[0]; + m_H1_mcPy=mcp->getMomentum()[1]; + m_H1_mcPz=mcp->getMomentum()[2]; + } + if(mcp->getPDG()==25 && m_H2_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==5000 ) { + index_H2=mcp; + m_H2_mcPDGID=mcp->getPDG(); + m_H2_mcE=mcp->getEnergy(); + m_H2_mcPx=mcp->getMomentum()[0]; + m_H2_mcPy=mcp->getMomentum()[1]; + m_H2_mcPz=mcp->getMomentum()[2]; + } + if(mcp->getPDG()==5 && m_b1_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==5000 && (mcp->getParents()[0]==index_H1 || mcp->getParents()[1]==index_H1) ) { + m_b1_mcPDGID=mcp->getPDG(); + m_b1_mcE=mcp->getEnergy(); + m_b1_mcPx=mcp->getMomentum()[0]; + m_b1_mcPy=mcp->getMomentum()[1]; + m_b1_mcPz=mcp->getMomentum()[2]; + continue; + } + if(mcp->getPDG()==-5 && m_b2_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==5000 && (mcp->getParents()[0]==index_H1 || mcp->getParents()[1]==index_H1) ) { + m_b2_mcPDGID=mcp->getPDG(); + m_b2_mcE=mcp->getEnergy(); + m_b2_mcPx=mcp->getMomentum()[0]; + m_b2_mcPy=mcp->getMomentum()[1]; + m_b2_mcPz=mcp->getMomentum()[2]; + continue; + } + if(mcp->getPDG()==5 && m_b3_mcE<0&& (mcp->getParents()[0]==index_H2)) { //&& mcp->getGeneratorStatus()==23 + m_b3_mcPDGID=mcp->getPDG(); + m_b3_mcE=mcp->getEnergy(); + m_b3_mcPx=mcp->getMomentum()[0]; + m_b3_mcPy=mcp->getMomentum()[1]; + m_b3_mcPz=mcp->getMomentum()[2]; + continue; + } + if(mcp->getPDG()==-5 && m_b4_mcE<0&& (mcp->getParents()[0]==index_H2) ){ //&& mcp->getGeneratorStatus()==23 + m_b4_mcPDGID=mcp->getPDG(); + m_b4_mcE=mcp->getEnergy(); + m_b4_mcPx=mcp->getMomentum()[0]; + m_b4_mcPy=mcp->getMomentum()[1]; + m_b4_mcPz=mcp->getMomentum()[2]; + continue; + } + } + if (m_processName == "bc") { + if((abs(mcp->getPDG())==4 || abs(mcp->getPDG())==5) && mcp->getGeneratorStatus()==23 && first_filled==true) { + m_b1_mcPDGID=mcp->getPDG(); + m_b1_mcE=mcp->getEnergy(); + m_b1_mcPx=mcp->getMomentum()[0]; + m_b1_mcPy=mcp->getMomentum()[1]; + m_b1_mcPz=mcp->getMomentum()[2]; + first_filled=false; + continue; + } + if((abs(mcp->getPDG())==4 || abs(mcp->getPDG())==5) && mcp->getGeneratorStatus()==23 && first_filled==false ){ + m_b2_mcPDGID=mcp->getPDG(); + m_b2_mcE=mcp->getEnergy(); + m_b2_mcPx=mcp->getMomentum()[0]; + m_b2_mcPy=mcp->getMomentum()[1]; + m_b2_mcPz=mcp->getMomentum()[2]; + } + } + + if (m_processName == "qqH") { + if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==1500 ){ + index_H1=mcp; + m_H1_mcE=mcp->getEnergy(); + m_H1_mcPx=mcp->getMomentum()[0]; + m_H1_mcPy=mcp->getMomentum()[1]; + m_H1_mcPz=mcp->getMomentum()[2]; + } + if(mcp->getPDG()==25 && m_H2_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==1500 ) { + index_H2=mcp; + m_H2_mcPDGID=mcp->getPDG(); + m_H2_mcE=mcp->getEnergy(); + m_H2_mcPx=mcp->getMomentum()[0]; + m_H2_mcPy=mcp->getMomentum()[1]; + m_H2_mcPz=mcp->getMomentum()[2]; + } + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b1_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==1500 && (mcp->getParents()[0]==index_H1 || mcp->getParents()[1]==index_H1) ) { + index_b1=mcp; + m_b1_mcPDGID=mcp->getPDG(); + m_b1_mcE=mcp->getEnergy(); + m_b1_mcPx=mcp->getMomentum()[0]; + m_b1_mcPy=mcp->getMomentum()[1]; + m_b1_mcPz=mcp->getMomentum()[2]; + continue; } - if(m_d2_mcE<0 && (abs(mcp->getPDG())<7 && mcp->getPDG()==(-m_d1_mcPDGID))){ - m_d2_mcPDGID=mcp->getPDG(); - m_d2_mcE=mcp->getEnergy(); - m_d2_mcPx=mcp->getMomentum()[0]; - m_d2_mcPy=mcp->getMomentum()[1]; - m_d2_mcPz=mcp->getMomentum()[2]; + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b2_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==1500 && (mcp->getParents()[0]==index_H1 || mcp->getParents()[1]==index_H1) && mcp!=index_b1 ) { + m_b2_mcPDGID=mcp->getPDG(); + m_b2_mcE=mcp->getEnergy(); + m_b2_mcPx=mcp->getMomentum()[0]; + m_b2_mcPy=mcp->getMomentum()[1]; + m_b2_mcPz=mcp->getMomentum()[2]; + continue; } + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b5_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==1500 && (mcp->getParents()[0]==index_H1 || mcp->getParents()[1]==index_H1) && mcp!=index_b1 && mcp!=index_b2 ) { + index_b5=mcp; + m_b5_mcPDGID=mcp->getPDG(); + m_b5_mcE=mcp->getEnergy(); + m_b5_mcPx=mcp->getMomentum()[0]; + m_b5_mcPy=mcp->getMomentum()[1]; + m_b5_mcPz=mcp->getMomentum()[2]; + continue; + } + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b6_mcE<0 && mcp->getParents()[0]->getMomentum()[0]==0 && mcp->getParents()[0]->getMomentum()[1]==0 && abs(mcp->getParents()[0]->getMomentum()[2])==1500 && (mcp->getParents()[0]==index_H1 || mcp->getParents()[1]==index_H1) && mcp!=index_b1 && mcp!=index_b2 && mcp!=index_b5) { + m_b6_mcPDGID=mcp->getPDG(); + m_b6_mcE=mcp->getEnergy(); + m_b6_mcPx=mcp->getMomentum()[0]; + m_b6_mcPy=mcp->getMomentum()[1]; + m_b6_mcPz=mcp->getMomentum()[2]; + continue; + } + if(mcp->getPDG()==5 && m_b3_mcE<0&& (mcp->getParents()[0]==index_H2)) { + index_b3=mcp; + m_b3_mcPDGID=mcp->getPDG(); + m_b3_mcE=mcp->getEnergy(); + m_b3_mcPx=mcp->getMomentum()[0]; + m_b3_mcPy=mcp->getMomentum()[1]; + m_b3_mcPz=mcp->getMomentum()[2]; + continue; + } + if(mcp->getPDG()==-5 && m_b4_mcE<0&& (mcp->getParents()[0]==index_H2) ){ + index_b4=mcp; + m_b4_mcPDGID=mcp->getPDG(); + m_b4_mcE=mcp->getEnergy(); + m_b4_mcPx=mcp->getMomentum()[0]; + m_b4_mcPy=mcp->getMomentum()[1]; + m_b4_mcPz=mcp->getMomentum()[2]; + continue; + } + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b7_mcE<0&& (mcp->getParents()[0]==index_H2) && mcp!= index_b3 && mcp!= index_b4 ) { + index_b7=mcp; + m_b7_mcPDGID=mcp->getPDG(); + m_b7_mcE=mcp->getEnergy(); + m_b7_mcPx=mcp->getMomentum()[0]; + m_b7_mcPy=mcp->getMomentum()[1]; + m_b7_mcPz=mcp->getMomentum()[2]; + continue; + } + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b8_mcE<0&& (mcp->getParents()[0]==index_H2)&& mcp!= index_b3 && mcp!= index_b4 && mcp!=index_b7 ){ + m_b8_mcPDGID=mcp->getPDG(); + m_b8_mcE=mcp->getEnergy(); + m_b8_mcPx=mcp->getMomentum()[0]; + m_b8_mcPy=mcp->getMomentum()[1]; + m_b8_mcPz=mcp->getMomentum()[2]; + continue; + } + } + + if (m_processName == "qqqqlnu") { + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4||abs(mcp->getPDG())==3||abs(mcp->getPDG())==2||abs(mcp->getPDG())==1) && m_b1_mcPDGID == -10) { + index_b1=mcp; + m_b1_mcPDGID=mcp->getPDG(); + m_b1_mcE=mcp->getEnergy(); + m_b1_mcPx=mcp->getMomentum()[0]; + m_b1_mcPy=mcp->getMomentum()[1]; + m_b1_mcPz=mcp->getMomentum()[2]; + } + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4||abs(mcp->getPDG())==3||abs(mcp->getPDG())==2||abs(mcp->getPDG())==1) && + mcp!=index_b1 && m_b2_mcPDGID == -10) { + index_b2=mcp; + m_b2_mcPDGID=mcp->getPDG(); + m_b2_mcE=mcp->getEnergy(); + m_b2_mcPx=mcp->getMomentum()[0]; + m_b2_mcPy=mcp->getMomentum()[1]; + m_b2_mcPz=mcp->getMomentum()[2]; + } + if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4||abs(mcp->getPDG())==3||abs(mcp->getPDG())==2||abs(mcp->getPDG())==1) && + mcp!=index_b1 && mcp!=index_b2 && m_b3_mcPDGID == -10) { + index_b3=mcp; + m_b3_mcPDGID=mcp->getPDG(); + m_b3_mcE=mcp->getEnergy(); + m_b3_mcPx=mcp->getMomentum()[0]; + m_b3_mcPy=mcp->getMomentum()[1]; + m_b3_mcPz=mcp->getMomentum()[2]; + } + if (m_b3_mcE>0 && (abs(mcp->getPDG())==5||abs(mcp->getPDG())==4||abs(mcp->getPDG())==3||abs(mcp->getPDG())==2||abs(mcp->getPDG())==1) && + mcp!=index_b1 && mcp!=index_b2 && mcp!=index_b3 && m_b4_mcPDGID == -10) { + index_b4=mcp; + m_b4_mcPDGID=mcp->getPDG(); + m_b4_mcE=mcp->getEnergy(); + m_b4_mcPx=mcp->getMomentum()[0]; + m_b4_mcPy=mcp->getMomentum()[1]; + m_b4_mcPz=mcp->getMomentum()[2]; + } + } + if(mcp->getGeneratorStatus()==1){//visible sum of stable particles --> take neutrinos out if(abs(mcp->getPDG())!=12 && abs(mcp->getPDG())!=14 && abs(mcp->getPDG())!=16){ m_E_trueVis+=mcp->getEnergy(); @@ -332,7 +879,7 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { } } - LCCollection* recoparticlecol = NULL; +/* LCCollection* recoparticlecol = NULL; recoparticlecol = evt->getCollection(m_inputRECOParticleCollection) ; if(recoparticlecol!=NULL){ //PandoraCandidate loop @@ -344,8 +891,8 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_pz_totPFO+=pandorapart->getMomentum()[2]; } } - - LCCollection* recojets = NULL; +*/ +/* LCCollection* recojets = NULL; recojets = evt->getCollection(m_recojetColName) ; if(recojets!=NULL){ //PandoraCandidate loop @@ -357,6 +904,8 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_recoJet_Pz->push_back(recojet->getMomentum()[2]); } } +*/ +/* LCCollection* genjets = NULL; genjets = evt->getCollection(m_genjetColName) ; if(genjets!=NULL){ @@ -369,6 +918,7 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_genJet_Pz->push_back(genjet->getMomentum()[2]); } } +*/ if(pass_W_boson_mass && pass_Z_boson_mass){ m_outputTree->Fill();