From 686014bfeb185295acf8741d667c0074dd644e90 Mon Sep 17 00:00:00 2001 From: gianelle Date: Tue, 14 Jan 2025 11:03:18 +0100 Subject: [PATCH 1/5] Add 4quarks data --- Calorimetry/include/JetAnalyzer.h | 39 +++++++ Calorimetry/src/JetAnalyzer.cc | 176 ++++++++++++++++++++++++++++++ 2 files changed, 215 insertions(+) diff --git a/Calorimetry/include/JetAnalyzer.h b/Calorimetry/include/JetAnalyzer.h index 28ac8b7..c6599d9 100644 --- a/Calorimetry/include/JetAnalyzer.h +++ b/Calorimetry/include/JetAnalyzer.h @@ -78,6 +78,45 @@ 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; + +//#######ENDNEW##################### + //only neutrinos float m_E_trueInv=0; float m_px_trueInv=0; diff --git a/Calorimetry/src/JetAnalyzer.cc b/Calorimetry/src/JetAnalyzer.cc index 6a3835d..e422d2f 100644 --- a/Calorimetry/src/JetAnalyzer.cc +++ b/Calorimetry/src/JetAnalyzer.cc @@ -129,6 +129,48 @@ void JetAnalyzer::init() { m_d2_mcPy = 0; m_d2_mcPz = 0; + +//#######NEW##################### + 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; + +//#######ENDNEW##################### + m_genJet_E ->clear(); m_genJet_Px ->clear(); m_genJet_Py ->clear(); @@ -167,6 +209,46 @@ 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 + 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"); + +//###############################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"); @@ -206,6 +288,9 @@ void JetAnalyzer::processRunHeader( LCRunHeader*) { void JetAnalyzer::processEvent( LCEvent* evt ) { + streamlog_out( MESSAGE )<<"!!! NEW VERION !!!" << std::endl; + + m_eventNumber=evt->getEventNumber(); m_runNumber=evt->getRunNumber(); @@ -236,6 +321,47 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_d2_mcPy=-10; m_d2_mcPz=-10; +//#######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; + +//#######ENDNEW##################### + m_genJet_E ->clear(); m_genJet_Px ->clear(); m_genJet_Py ->clear(); @@ -264,6 +390,8 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { if(mcColl!=NULL){ int boson_counter=0; int ind_second_boson=-1; + MCParticle* index_H1=NULL; //NEW + MCParticle* index_H2=NULL; 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 @@ -316,6 +444,54 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_d2_mcPy=mcp->getMomentum()[1]; m_d2_mcPz=mcp->getMomentum()[2]; } + + if(mcp->getPDG()==25 && m_H1_mcE<0) { + index_H1=mcp; //metto variabile H1 uguale all'indice corrente + m_H1_mcPDGID=mcp->getPDG(); + m_H1_mcE=mcp->getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m + 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 )){ //se è passato per l'if precedente, l'm è cambiato + index_H2=mcp; //metto variabile H2 uguale all'indice corrente + 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(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(); From b0408c14479717d20aea79c6254972d64d90cf9b Mon Sep 17 00:00:00 2001 From: Davide Zuliani INFN Padova - PostDoc Date: Thu, 10 Apr 2025 15:38:17 +0200 Subject: [PATCH 2/5] initial commit --- Calorimetry/include/JetAnalyzer.h | 40 ++- Calorimetry/src/JetAnalyzer.cc | 524 +++++++++++++++++++++++++++--- 2 files changed, 517 insertions(+), 47 deletions(-) diff --git a/Calorimetry/include/JetAnalyzer.h b/Calorimetry/include/JetAnalyzer.h index c6599d9..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; @@ -115,6 +116,31 @@ class JetAnalyzer : public Processor { 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 @@ -129,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; @@ -142,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; @@ -151,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 e422d2f..eefbf66 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; @@ -130,7 +135,7 @@ void JetAnalyzer::init() { m_d2_mcPz = 0; -//#######NEW##################### +//#######NEW for HH##################### m_H1_mcPDGID = 0; m_H1_mcE = 0; m_H1_mcPx = 0; @@ -169,9 +174,40 @@ void JetAnalyzer::init() { 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; + + + //#######ENDNEW##################### - m_genJet_E ->clear(); +//##########NEW FOR DIJET########### +/* 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_genJet_E ->clear(); m_genJet_Px ->clear(); m_genJet_Py ->clear(); m_genJet_Pz ->clear(); @@ -180,7 +216,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(); @@ -209,7 +245,7 @@ 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 + //########################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"); @@ -246,8 +282,31 @@ void JetAnalyzer::init() { m_outputTree->Branch("b4_mcPy",&m_b4_mcPy,"b4_mcPy/F"); m_outputTree->Branch("b4_mcPz",&m_b4_mcPz,"b4_mcPz/F"); -//###############################ENDNEW + 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"); @@ -262,7 +321,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"); @@ -277,7 +336,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); - +*/ } @@ -288,9 +347,6 @@ void JetAnalyzer::processRunHeader( LCRunHeader*) { void JetAnalyzer::processEvent( LCEvent* evt ) { - streamlog_out( MESSAGE )<<"!!! NEW VERION !!!" << std::endl; - - m_eventNumber=evt->getEventNumber(); m_runNumber=evt->getRunNumber(); @@ -304,11 +360,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; @@ -322,7 +378,7 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_d2_mcPz=-10; //#######NEW##################### - m_H1_mcPDGID = -10; + m_H1_mcPDGID = -10; m_H1_mcE = -10; m_H1_mcPx = -10; m_H1_mcPy = -10; @@ -360,9 +416,31 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { 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_E ->clear(); m_genJet_Px ->clear(); m_genJet_Py ->clear(); m_genJet_Pz ->clear(); @@ -371,7 +449,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(); @@ -387,11 +465,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 @@ -444,54 +530,409 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_d2_mcPy=mcp->getMomentum()[1]; m_d2_mcPz=mcp->getMomentum()[2]; } - - if(mcp->getPDG()==25 && m_H1_mcE<0) { - index_H1=mcp; //metto variabile H1 uguale all'indice corrente - m_H1_mcPDGID=mcp->getPDG(); +//#######################SAVE HH#################### + // if(mcp->getPDG()==25 && m_H1_mcE<0) { + // index_H1=mcp; //metto variabile H1 uguale all'indice corrente + // m_H1_mcPDGID=mcp->getPDG(); + // m_H1_mcE=mcp->getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m + // 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 )){ //se è passato per l'if precedente, l'm è cambiato + // index_H2=mcp; //metto variabile H2 uguale all'indice corrente + // 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]; + // } +//#######################SAVE bb FOR DIJET#################### +/* +if(mcp->getPDG()==4 ) { //&& mcp->getGeneratorStatus()==23 + 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 ){ //&& mcp->getGeneratorStatus()==23 + 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]; + } +*/ + +//#######################SAVE bbbb FOR HH bkg#################### + +// if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==5000 ) //&& mcp->getMomentum()[2]==1500 +// { std::cout <getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m +// 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 ) //&& mcp->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(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) ) { //&& mcp->getGeneratorStatus()==23 + +// 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) ) { //&& mcp->getGeneratorStatus()==23 +// 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)) { //&& 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_H1 || 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; +// } + + +//#####################################SAVE bbH##################################### +// if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==5000 ) //&& mcp->getMomentum()[2]==1500 +// { std::cout <getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m +// 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 ) //&& mcp->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(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) ) { //&& mcp->getGeneratorStatus()==23 + +// 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]; +// std::cout << m_b1_mcPx << std::endl; +// 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) ) { //&& mcp->getGeneratorStatus()==23 +// 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]; +// std::cout << m_b2_mcPx << std::endl; + +// 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]; +// std::cout << m_b3_mcPx << std::endl; + +// 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]; +// std::cout << m_b4_mcPx << std::endl; + +// continue; +// } + + +//###########################################SAVE bc for light jets################### + /* 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]; + } +*/ +//########################################### +//######################SAVE qqH #################################################### +/* + if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==1500 ) //&& mcp->getMomentum()[2]==1500 + { std::cout <getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m 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 )){ //se è passato per l'if precedente, l'm è cambiato - index_H2=mcp; //metto variabile H2 uguale all'indice corrente + } + + 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 ) //&& mcp->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) ) { //&& mcp->getGeneratorStatus()==23 + 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]; + std::cout << m_b1_mcPx << std::endl; + continue; + } +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 ) { //&& mcp->getGeneratorStatus()==23 + 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]; + std::cout << m_b2_mcPx << std::endl; + + 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 ) { //&& mcp->getGeneratorStatus()==23 + 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]; + std::cout << m_b5_mcPx << std::endl; + 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) { //&& mcp->getGeneratorStatus()==23 + 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]; + std::cout << m_b6_mcPx << std::endl; + + continue; } - if(abs(mcp->getPDG())==5 && m_b1_mcE<0 && mcp->getParents()[0]==index_H1) { + + + if(mcp->getPDG()==5 && m_b3_mcE<0&& (mcp->getParents()[0]==index_H2)) { //&& mcp->getGeneratorStatus()==23 + 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]; + std::cout << m_b3_mcPx << std::endl; + + continue; + } + if(mcp->getPDG()==-5 && m_b4_mcE<0&& (mcp->getParents()[0]==index_H2) ){ //&& mcp->getGeneratorStatus()==23 + 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]; + std::cout << m_b4_mcPx << std::endl; + + 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 ) { //&& mcp->getGeneratorStatus()==23 + 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]; + std::cout << m_b7_mcPx << std::endl; + + 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 ){ //&& mcp->getGeneratorStatus()==23 + 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]; + std::cout << m_b8_mcPx << std::endl; + + continue; + } +*/ + + + + + + + + +//#######################SAVE qqqqlnu FOR HH bkg#################### + + // if(mcp->getPDG()==23 && m_H1_mcE<0) { + // index_H1=mcp; //metto variabile H1 uguale all'indice corrente + // m_H1_mcPDGID=mcp->getPDG(); + // m_H1_mcE=mcp->getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m + // m_H1_mcPx=mcp->getMomentum()[0]; + // m_H1_mcPy=mcp->getMomentum()[1]; + // m_H1_mcPz=mcp->getMomentum()[2]; + // } + + + // if(abs(mcp->getPDG())==24 && index_H1!=mcp){ + // index_H2=mcp; //metto variabile H2 uguale all'indice corrente + // 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||mcp->getPDG()==4||mcp->getPDG()==3||mcp->getPDG()==2||mcp->getPDG()==1) && m_b1_mcE<0 && mcp->getParents()[0]==index_H1) { + 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]; + std::cout << "b1: " << m_b1_mcPDGID << " " << m_b1_mcE << " " << m_b1_mcPx << " " << m_b1_mcPy << " " << m_b1_mcPz << std::endl; } - if(m_b2_mcE<0 && (abs(mcp->getPDG())==5 && mcp->getPDG()==(-m_b1_mcPDGID)) && mcp->getParents()[0]==index_H1 ){ + // if(m_b2_mcE<0 && (mcp->getPDG()==-5||mcp->getPDG()==-4||mcp->getPDG()==-3||mcp->getPDG()==-2||mcp->getPDG()==-1) && mcp->getPDG()==(-m_b1_mcPDGID) && mcp->getParents()[0]==index_H1 ){ + 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]; + std::cout << "b2: " << m_b2_mcPDGID << " " << m_b2_mcE << " " << m_b2_mcPx << " " << m_b2_mcPy << " " << m_b2_mcPz << std::endl; } - if(abs(mcp->getPDG())==5 && m_b3_mcE<0 && mcp->getParents()[0]==index_H2) { + // if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4||abs(mcp->getPDG())==3||abs(mcp->getPDG())==2||abs(mcp->getPDG())==1) && m_b3_mcE<0 && mcp->getParents()[0]==index_H2) { + 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]; + std::cout << "b3: " << m_b3_mcPDGID << " " << m_b3_mcE << " " << m_b3_mcPx << " " << m_b3_mcPy << " " << m_b3_mcPz << std::endl; } - if(m_b4_mcE<0 && (abs(mcp->getPDG())==5 && mcp->getPDG()==(-m_b3_mcPDGID))&& mcp->getParents()[0]==index_H2){ - m_b4_mcPDGID=mcp->getPDG(); + // if(m_b4_mcE<0 && (abs(mcp->getPDG())==5||abs(mcp->getPDG())==4||abs(mcp->getPDG())==3||abs(mcp->getPDG())==2||abs(mcp->getPDG())==1) && (mcp->getPDG()!=m_b3_mcPDGID)&& mcp->getParents()[0]==index_H2){ + 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]; + std::cout << "b4: " << m_b4_mcPDGID << " " << m_b4_mcE << " " << m_b4_mcPx << " " << m_b4_mcPy << " " << m_b4_mcPz << std::endl; } + + +//############################################################################################### + + + + + 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(); @@ -508,7 +949,7 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { } } - LCCollection* recoparticlecol = NULL; +/* LCCollection* recoparticlecol = NULL; recoparticlecol = evt->getCollection(m_inputRECOParticleCollection) ; if(recoparticlecol!=NULL){ //PandoraCandidate loop @@ -520,8 +961,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 @@ -533,6 +974,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){ @@ -545,6 +988,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(); From fb13d4a451d93d2d54a1eafb795aac66d2fad7f3 Mon Sep 17 00:00:00 2001 From: Davide Zuliani INFN Padova - PostDoc Date: Mon, 14 Apr 2025 16:56:54 +0200 Subject: [PATCH 3/5] update JetAnalyzer tool --- Calorimetry/src/JetAnalyzer.cc | 498 +++++++++++++-------------------- 1 file changed, 196 insertions(+), 302 deletions(-) diff --git a/Calorimetry/src/JetAnalyzer.cc b/Calorimetry/src/JetAnalyzer.cc index eefbf66..436e0f6 100644 --- a/Calorimetry/src/JetAnalyzer.cc +++ b/Calorimetry/src/JetAnalyzer.cc @@ -134,8 +134,6 @@ void JetAnalyzer::init() { m_d2_mcPy = 0; m_d2_mcPz = 0; - -//#######NEW for HH##################### m_H1_mcPDGID = 0; m_H1_mcE = 0; m_H1_mcPx = 0; @@ -148,7 +146,6 @@ void JetAnalyzer::init() { m_H2_mcPy = 0; m_H2_mcPz = 0; - m_b1_mcPDGID = 0; m_b1_mcE = 0; m_b1_mcPx = 0; @@ -161,7 +158,6 @@ void JetAnalyzer::init() { m_b2_mcPy = 0; m_b2_mcPz = 0; - m_b3_mcPDGID = 0; m_b3_mcE = 0; m_b3_mcPx = 0; @@ -174,7 +170,6 @@ void JetAnalyzer::init() { m_b4_mcPy = 0; m_b4_mcPz = 0; - m_b5_mcPDGID = 0; m_b5_mcE = 0; m_b5_mcPx = 0; @@ -187,24 +182,17 @@ void JetAnalyzer::init() { 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; - -//#######ENDNEW##################### - -//##########NEW FOR DIJET########### -/* 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_b8_mcPDGID = 0; + m_b8_mcE = 0; + m_b8_mcPx = 0; + m_b8_mcPy = 0; + m_b8_mcPz = 0; //####################################### /* m_genJet_E ->clear(); @@ -530,204 +518,179 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_d2_mcPy=mcp->getMomentum()[1]; m_d2_mcPz=mcp->getMomentum()[2]; } -//#######################SAVE HH#################### - // if(mcp->getPDG()==25 && m_H1_mcE<0) { - // index_H1=mcp; //metto variabile H1 uguale all'indice corrente - // m_H1_mcPDGID=mcp->getPDG(); - // m_H1_mcE=mcp->getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m - // 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 )){ //se è passato per l'if precedente, l'm è cambiato - // index_H2=mcp; //metto variabile H2 uguale all'indice corrente - // 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]; - // } -//#######################SAVE bb FOR DIJET#################### -/* -if(mcp->getPDG()==4 ) { //&& mcp->getGeneratorStatus()==23 + + 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()==4 ) { //&& mcp->getGeneratorStatus()==23 + 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 ){ //&& mcp->getGeneratorStatus()==23 + 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]; + } + } + +//#######################SAVE bbbb FOR HH bkg#################### + 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()==-4 ){ //&& mcp->getGeneratorStatus()==23 + 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; } -*/ - -//#######################SAVE bbbb FOR HH bkg#################### - -// if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==5000 ) //&& mcp->getMomentum()[2]==1500 -// { std::cout <getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m -// 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 ) //&& mcp->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(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) ) { //&& mcp->getGeneratorStatus()==23 - -// 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) ) { //&& mcp->getGeneratorStatus()==23 -// 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)) { //&& 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_H1 || 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; -// } - - -//#####################################SAVE bbH##################################### -// if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==5000 ) //&& mcp->getMomentum()[2]==1500 -// { std::cout <getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m -// 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 ) //&& mcp->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(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) ) { //&& mcp->getGeneratorStatus()==23 - -// 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]; -// std::cout << m_b1_mcPx << std::endl; -// 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) ) { //&& mcp->getGeneratorStatus()==23 -// 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]; -// std::cout << m_b2_mcPx << std::endl; - -// 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]; -// std::cout << m_b3_mcPx << std::endl; - -// 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]; -// std::cout << m_b4_mcPx << std::endl; - -// continue; -// } - - -//###########################################SAVE bc for light jets################### - /* if((abs(mcp->getPDG())==4 || abs(mcp->getPDG())==5) && mcp->getGeneratorStatus()==23 && first_filled==true) { + 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; + continue; } if((abs(mcp->getPDG())==4 || abs(mcp->getPDG())==5) && mcp->getGeneratorStatus()==23 && first_filled==false ){ m_b2_mcPDGID=mcp->getPDG(); @@ -736,151 +699,97 @@ if(mcp->getPDG()==4 ) { //&& mcp->getGeneratorStatus()==23 m_b2_mcPy=mcp->getMomentum()[1]; m_b2_mcPz=mcp->getMomentum()[2]; } -*/ -//########################################### -//######################SAVE qqH #################################################### -/* - if(mcp->getPDG()==13 && m_H1_mcE<0 && mcp->getMomentum()[0]==0 && mcp->getMomentum()[1]==0 && abs(mcp->getMomentum()[2])==1500 ) //&& mcp->getMomentum()[2]==1500 - { std::cout <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_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m + 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 ) //&& mcp->getMomentum()[2]==1500 - { + } + 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) ) { //&& mcp->getGeneratorStatus()==23 + } + 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]; - std::cout << m_b1_mcPx << std::endl; - continue; + continue; } -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 ) { //&& mcp->getGeneratorStatus()==23 + 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]; - std::cout << m_b2_mcPx << std::endl; - - continue; + 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 ) { //&& mcp->getGeneratorStatus()==23 + 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]; - std::cout << m_b5_mcPx << std::endl; - continue; + 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) { //&& mcp->getGeneratorStatus()==23 + 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]; - std::cout << m_b6_mcPx << std::endl; - - continue; + continue; } - - - - if(mcp->getPDG()==5 && m_b3_mcE<0&& (mcp->getParents()[0]==index_H2)) { //&& mcp->getGeneratorStatus()==23 - index_b3=mcp; - m_b3_mcPDGID=mcp->getPDG(); + 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]; - std::cout << m_b3_mcPx << std::endl; - - continue; + continue; } - if(mcp->getPDG()==-5 && m_b4_mcE<0&& (mcp->getParents()[0]==index_H2) ){ //&& mcp->getGeneratorStatus()==23 - index_b4=mcp; + 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]; - std::cout << m_b4_mcPx << std::endl; - - continue; + 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 ) { //&& mcp->getGeneratorStatus()==23 - index_b7=mcp; + 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]; - std::cout << m_b7_mcPx << std::endl; - - continue; + 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 ){ //&& mcp->getGeneratorStatus()==23 + 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]; - std::cout << m_b8_mcPx << std::endl; - - continue; + continue; } -*/ - - - - - - - + } //#######################SAVE qqqqlnu FOR HH bkg#################### - - // if(mcp->getPDG()==23 && m_H1_mcE<0) { - // index_H1=mcp; //metto variabile H1 uguale all'indice corrente - // m_H1_mcPDGID=mcp->getPDG(); - // m_H1_mcE=mcp->getEnergy(); // m_H1_mcE non è più minore di 0, a m+1 anche se ho il secondo higgs, non entra in questo if e index_H1 resta pari a m - // m_H1_mcPx=mcp->getMomentum()[0]; - // m_H1_mcPy=mcp->getMomentum()[1]; - // m_H1_mcPz=mcp->getMomentum()[2]; - // } - - - // if(abs(mcp->getPDG())==24 && index_H1!=mcp){ - // index_H2=mcp; //metto variabile H2 uguale all'indice corrente - // 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||mcp->getPDG()==4||mcp->getPDG()==3||mcp->getPDG()==2||mcp->getPDG()==1) && m_b1_mcE<0 && mcp->getParents()[0]==index_H1) { + 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(); @@ -888,9 +797,7 @@ if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b6_mcE<0 && mcp->getParen m_b1_mcPx=mcp->getMomentum()[0]; m_b1_mcPy=mcp->getMomentum()[1]; m_b1_mcPz=mcp->getMomentum()[2]; - std::cout << "b1: " << m_b1_mcPDGID << " " << m_b1_mcE << " " << m_b1_mcPx << " " << m_b1_mcPy << " " << m_b1_mcPz << std::endl; } - // if(m_b2_mcE<0 && (mcp->getPDG()==-5||mcp->getPDG()==-4||mcp->getPDG()==-3||mcp->getPDG()==-2||mcp->getPDG()==-1) && mcp->getPDG()==(-m_b1_mcPDGID) && mcp->getParents()[0]==index_H1 ){ 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; @@ -899,10 +806,7 @@ if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b6_mcE<0 && mcp->getParen m_b2_mcPx=mcp->getMomentum()[0]; m_b2_mcPy=mcp->getMomentum()[1]; m_b2_mcPz=mcp->getMomentum()[2]; - std::cout << "b2: " << m_b2_mcPDGID << " " << m_b2_mcE << " " << m_b2_mcPx << " " << m_b2_mcPy << " " << m_b2_mcPz << std::endl; } - - // if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4||abs(mcp->getPDG())==3||abs(mcp->getPDG())==2||abs(mcp->getPDG())==1) && m_b3_mcE<0 && mcp->getParents()[0]==index_H2) { 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; @@ -911,9 +815,7 @@ if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b6_mcE<0 && mcp->getParen m_b3_mcPx=mcp->getMomentum()[0]; m_b3_mcPy=mcp->getMomentum()[1]; m_b3_mcPz=mcp->getMomentum()[2]; - std::cout << "b3: " << m_b3_mcPDGID << " " << m_b3_mcE << " " << m_b3_mcPx << " " << m_b3_mcPy << " " << m_b3_mcPz << std::endl; } - // if(m_b4_mcE<0 && (abs(mcp->getPDG())==5||abs(mcp->getPDG())==4||abs(mcp->getPDG())==3||abs(mcp->getPDG())==2||abs(mcp->getPDG())==1) && (mcp->getPDG()!=m_b3_mcPDGID)&& mcp->getParents()[0]==index_H2){ 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; @@ -922,16 +824,8 @@ if((abs(mcp->getPDG())==5||abs(mcp->getPDG())==4) && m_b6_mcE<0 && mcp->getParen m_b4_mcPx=mcp->getMomentum()[0]; m_b4_mcPy=mcp->getMomentum()[1]; m_b4_mcPz=mcp->getMomentum()[2]; - std::cout << "b4: " << m_b4_mcPDGID << " " << m_b4_mcE << " " << m_b4_mcPx << " " << m_b4_mcPy << " " << m_b4_mcPz << std::endl; } - - - -//############################################################################################### - - - - + } 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){ From c6971b6126576a9c6fdce846869a7a829bdac35e Mon Sep 17 00:00:00 2001 From: Davide Zuliani INFN Padova - PostDoc Date: Mon, 14 Apr 2025 17:46:00 +0200 Subject: [PATCH 4/5] add cc and qq --- Calorimetry/src/JetAnalyzer.cc | 42 ++++++++++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 4 deletions(-) diff --git a/Calorimetry/src/JetAnalyzer.cc b/Calorimetry/src/JetAnalyzer.cc index 436e0f6..a56206a 100644 --- a/Calorimetry/src/JetAnalyzer.cc +++ b/Calorimetry/src/JetAnalyzer.cc @@ -565,15 +565,33 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { m_b4_mcPz=mcp->getMomentum()[2]; } } + if (m_processName == "bb") { - if(mcp->getPDG()==4 ) { //&& mcp->getGeneratorStatus()==23 + 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 ){ //&& mcp->getGeneratorStatus()==23 + if(mcp->getPDG()==-4 ){ m_b2_mcPDGID=mcp->getPDG(); m_b2_mcE=mcp->getEnergy(); m_b2_mcPx=mcp->getMomentum()[0]; @@ -581,8 +599,25 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { 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]; + } + } + -//#######################SAVE bbbb FOR HH bkg#################### 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; @@ -788,7 +823,6 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { } } -//#######################SAVE qqqqlnu FOR HH bkg#################### 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; From b0a8130299902ae58a02581ce596cbac5861bd46 Mon Sep 17 00:00:00 2001 From: Davide Zuliani INFN Padova - PostDoc Date: Thu, 17 Apr 2025 09:52:03 +0200 Subject: [PATCH 5/5] =?UTF-8?q?=C3=82finished=20jetanalyzer,=20adding=20al?= =?UTF-8?q?so=20dijet?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Calorimetry/src/JetAnalyzer.cc | 128 +++++++++++++++++---------------- 1 file changed, 65 insertions(+), 63 deletions(-) diff --git a/Calorimetry/src/JetAnalyzer.cc b/Calorimetry/src/JetAnalyzer.cc index a56206a..baa2065 100644 --- a/Calorimetry/src/JetAnalyzer.cc +++ b/Calorimetry/src/JetAnalyzer.cc @@ -503,20 +503,22 @@ 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]; - } - 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]; + //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") { @@ -566,56 +568,56 @@ void JetAnalyzer::processEvent( LCEvent* evt ) { } } - 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 == "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 == "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") {