From 27b32e6c1d9deff13481b0fe32778c3e991d5184 Mon Sep 17 00:00:00 2001 From: Nazim Huseynov Date: Wed, 2 Oct 2019 09:29:39 +0200 Subject: [PATCH 1/2] tH additional truth variables --- .../TopPartons/Root/CalcThqPartonHistory.cxx | 298 ++++++++++++++++++ .../TopPartons/CalcThqPartonHistory.h | 80 +++++ 2 files changed, 378 insertions(+) create mode 100644 PhysicsAnalysis/TopPhys/xAOD/TopPartons/Root/CalcThqPartonHistory.cxx create mode 100644 PhysicsAnalysis/TopPhys/xAOD/TopPartons/TopPartons/CalcThqPartonHistory.h diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopPartons/Root/CalcThqPartonHistory.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopPartons/Root/CalcThqPartonHistory.cxx new file mode 100644 index 00000000000..9512248a356 --- /dev/null +++ b/PhysicsAnalysis/TopPhys/xAOD/TopPartons/Root/CalcThqPartonHistory.cxx @@ -0,0 +1,298 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#include "TopPartons/CalcThqPartonHistory.h" +#include "TopPartons/CalcTopPartonHistory.h" +#include "TopConfiguration/TopConfig.h" + + +namespace top { + + CalcThqPartonHistory::CalcThqPartonHistory(const std::string& name) : CalcTopPartonHistory(name) {} + const xAOD::TruthParticle* CalcThqPartonHistory::findAfterGamma(const xAOD::TruthParticle* particle) { + bool isAfterGamma(false); + const int particle_ID = particle->pdgId(); + int forLoop; + while (!isAfterGamma) { + forLoop = 0; + for (size_t j = 0; j < particle->nChildren(); j++) { + const xAOD::TruthParticle* tmp_children = particle->child(j); + if (tmp_children && tmp_children->pdgId() == particle_ID && tmp_children->pdgId() != 22) { + particle = particle->child(j); + forLoop++; + break; + }//if + }//for + + if (forLoop == 0) isAfterGamma = true; + }//while + return particle; + } + int CalcThqPartonHistory::sign(int a) { + if (a < 0) { return -1; } + else return 1; + } + bool CalcThqPartonHistory::bottom(const xAOD::TruthParticleContainer* truthParticles, int start) { + for (const xAOD::TruthParticle* particle : *truthParticles) { + if (particle->pdgId() != start) { continue; } + tH.b_p4 = particle->p4(); + tH.b_pdgId = particle->pdgId(); + return true; + } + return false; + } + bool CalcThqPartonHistory::Higgstautau(const xAOD::TruthParticleContainer* truthParticles, int start) { + bool has_Higgs = false; + bool has_tau1_neutrino = false; + bool has_tau2_neutrino = false; + bool hadr_tau1 = false; + bool hadr_tau2 = false; + for (const xAOD::TruthParticle* particle : *truthParticles) { + if (particle->pdgId() != start || particle->nChildren() != 2) { continue; } + tH.Higgs_p4 = particle->p4(); + has_Higgs = true; + for (size_t k = 0; k < particle->nChildren(); k++) { + const xAOD::TruthParticle* HiggsChildren = particle->child(k); + if (fabs(HiggsChildren->pdgId()) == 15) {// demanding tau-leptons as childen + const xAOD::TruthParticle* tau = CalcThqPartonHistory::findAfterGamma(HiggsChildren); + if (k == 0) { + tH.Tau1_from_Higgs_p4 = tau->p4(); + tH.Tau1_from_Higgs_pdgId = tau->pdgId(); + } + else { + tH.Tau2_from_Higgs_p4 = tau->p4(); + tH.Tau2_from_Higgs_pdgId = tau->pdgId(); + } + for (size_t q = 0; q < tau->nChildren(); ++q) { + const xAOD::TruthParticle* tauChildren = tau->child(q); + if (fabs(tauChildren->pdgId()) == 16) {// tau neutrino + if (k == 0) { + tH.nu_from_Tau1_p4 = tauChildren->p4(); + tH.nu_from_Tau1_pdgId = tauChildren->pdgId(); + has_tau1_neutrino = true; + } + else { + tH.nu_from_Tau2_p4 = tauChildren->p4(); + tH.nu_from_Tau2_pdgId = tauChildren->pdgId(); + has_tau2_neutrino = true; + } + } + else if (fabs(tauChildren->pdgId()) >= 11 && fabs(tauChildren->pdgId()) <= 14) { //light leptons + if (fabs(tauChildren->pdgId()) == 11 || fabs(tauChildren->pdgId()) == 13) { //electron or muon + if (k == 0) { + tH.W_decay1_from_Tau1_p4 = tauChildren->p4(); + tH.W_decay1_from_Tau1_pdgId = tauChildren->pdgId(); + } + else { + tH.W_decay1_from_Tau2_p4 = tauChildren->p4(); + tH.W_decay1_from_Tau2_pdgId = tauChildren->pdgId(); + } + } + else if (fabs(tauChildren->pdgId()) == 12 || fabs(tauChildren->pdgId()) == 14) { // electron or muon neutrino + if (k == 0) { + tH.W_decay2_from_Tau1_p4 = tauChildren->p4(); + tH.W_decay2_from_Tau1_pdgId = tauChildren->pdgId(); + } + else { + tH.W_decay2_from_Tau2_p4 = tauChildren->p4(); + tH.W_decay2_from_Tau2_pdgId = tauChildren->pdgId(); + } + } + } + else { // if a particle passes the criteria above, it has to be a hadron. + if (k == 0) { hadr_tau1 = true; } + else { hadr_tau2 = true; } + }// else + } // for + }//if + } //for + } + if (has_Higgs && has_tau1_neutrino && has_tau2_neutrino) { + if (hadr_tau1) { //convention: store hadr. decaying W-Boson as Wdecay1, set all parameters of Wdecay2 to 0. + tH.W_decay1_from_Tau1_p4 = tH.Tau1_from_Higgs_p4 - tH.nu_from_Tau1_p4; + tH.W_decay1_from_Tau1_pdgId = -24 * sign(tH.nu_from_Tau1_pdgId); + tH.W_decay2_from_Tau1_p4 = { 0, 0, 0, 0 }; + tH.W_decay2_from_Tau1_pdgId = 0; + tH.TauJets1 = 1; + } + else if (hadr_tau2) { + tH.W_decay1_from_Tau2_p4 = tH.Tau2_from_Higgs_p4 - tH.nu_from_Tau2_p4; + tH.W_decay1_from_Tau2_pdgId = -24 * sign(tH.nu_from_Tau2_pdgId); + tH.W_decay2_from_Tau2_p4 = { 0, 0, 0, 0 }; + tH.W_decay2_from_Tau2_pdgId = 0; + tH.TauJets2 = 1; + } + return true; + } + tH.TauJets1 = 0; + tH.TauJets2 = 0; + return false; + } + void CalcThqPartonHistory::THHistorySaver(const xAOD::TruthParticleContainer* truthParticles, xAOD::PartonHistory* ThqPartonHistory) { + + ThqPartonHistory->IniVarThq(); + TLorentzVector t_before, t_after, t_after_SC; + TLorentzVector Wp; + TLorentzVector b; + TLorentzVector WpDecay1; + TLorentzVector WpDecay2; + int WpDecay1_pdgId; + int WpDecay2_pdgId; + + bool event_top = CalcTopPartonHistory::topWb(truthParticles, 6, t_before, t_after, Wp, b, WpDecay1, WpDecay1_pdgId, WpDecay2, WpDecay2_pdgId); + bool event_top_SC = CalcTopPartonHistory::topAfterFSR_SC(truthParticles, 6, t_after_SC); + bool event_topbar = CalcTopPartonHistory::topWb(truthParticles, -6, t_before, t_after, Wp, b, WpDecay1, WpDecay1_pdgId, WpDecay2, WpDecay2_pdgId); + bool event_topbar_SC = CalcTopPartonHistory::topAfterFSR_SC(truthParticles, -6, t_after_SC); + bool event_Higgs = CalcThqPartonHistory::Higgstautau(truthParticles, 25); + bool event_bottom = CalcThqPartonHistory::bottom(truthParticles, 5); + bool event_bottombar = CalcThqPartonHistory::bottom(truthParticles, -5); + + + if (event_Higgs ) { + if ((event_top && event_bottombar) || (event_topbar && event_bottom)) + { + ThqPartonHistory->auxdecor< float >("MC_t_beforeFSR_m") = t_before.M(); + ThqPartonHistory->auxdecor< float >("MC_t_beforeFSR_pt") = t_before.Pt(); + ThqPartonHistory->auxdecor< float >("MC_t_beforeFSR_phi") = t_before.Phi(); + fillEtaBranch(ThqPartonHistory, "MC_t_beforeFSR_eta", t_before); + + ThqPartonHistory->auxdecor< float >("MC_t_afterFSR_m") = t_after.M(); + ThqPartonHistory->auxdecor< float >("MC_t_afterFSR_pt") = t_after.Pt(); + ThqPartonHistory->auxdecor< float >("MC_t_afterFSR_phi") = t_after.Phi(); + fillEtaBranch(ThqPartonHistory, "MC_t_afterFSR_eta", t_after); + + if (event_top_SC || event_topbar_SC) { + ThqPartonHistory->auxdecor< float >("MC_t_afterFSR_SC_m") = t_after_SC.M(); + ThqPartonHistory->auxdecor< float >("MC_t_afterFSR_SC_pt") = t_after_SC.Pt(); + ThqPartonHistory->auxdecor< float >("MC_t_afterFSR_SC_phi") = t_after_SC.Phi(); + fillEtaBranch(ThqPartonHistory, "MC_t_afterFSR_SC_eta", t_after_SC); + } + + ThqPartonHistory->auxdecor< float >("MC_W_from_t_m") = Wp.M(); + ThqPartonHistory->auxdecor< float >("MC_W_from_t_pt") = Wp.Pt(); + ThqPartonHistory->auxdecor< float >("MC_W_from_t_phi") = Wp.Phi(); + fillEtaBranch(ThqPartonHistory, "MC_W_from_t_eta", Wp); + + ThqPartonHistory->auxdecor< float >("MC_b_from_t_m") = b.M(); + ThqPartonHistory->auxdecor< float >("MC_b_from_t_pt") = b.Pt(); + ThqPartonHistory->auxdecor< float >("MC_b_from_t_phi") = b.Phi(); + fillEtaBranch(ThqPartonHistory, "MC_b_from_t_eta", b); + + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_t_m") = WpDecay1.M(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_t_pt") = WpDecay1.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_t_phi") = WpDecay1.Phi(); + ThqPartonHistory->auxdecor< int >("MC_Wdecay1_from_t_pdgId") = WpDecay1_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_Wdecay1_from_t_eta", WpDecay1); + + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_t_m") = WpDecay2.M(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_t_pt") = WpDecay2.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_t_phi") = WpDecay2.Phi(); + ThqPartonHistory->auxdecor< int >("MC_Wdecay2_from_t_pdgId") = WpDecay2_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_Wdecay2_from_t_eta", WpDecay2); + + //Higgs-Variables + ThqPartonHistory->auxdecor< float >("MC_Higgs_m") = tH.Higgs_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_Higgs_pt") = tH.Higgs_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Higgs_phi") = tH.Higgs_p4.Phi(); + fillEtaBranch(ThqPartonHistory, "MC_Higgs_eta", tH.Higgs_p4); + + //First Tau (Tau1) + ThqPartonHistory->auxdecor< float >("MC_Tau1_from_Higgs_m") = tH.Tau1_from_Higgs_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_Tau1_from_Higgs_pt") = tH.Tau1_from_Higgs_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Tau1_from_Higgs_phi") = tH.Tau1_from_Higgs_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_Tau1_from_Higgs_pdgId") = tH.Tau1_from_Higgs_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_Tau1_from_Higgs_eta", tH.Tau1_from_Higgs_p4); + + //Second Tau (Tau2) + ThqPartonHistory->auxdecor< float >("MC_Tau2_from_Higgs_m") = tH.Tau2_from_Higgs_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_Tau2_from_Higgs_pt") = tH.Tau2_from_Higgs_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Tau2_from_Higgs_phi") = tH.Tau2_from_Higgs_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_Tau2_from_Higgs_pdgId") = tH.Tau2_from_Higgs_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_Tau2_from_Higgs_eta", tH.Tau2_from_Higgs_p4); + + //neutrino from first Tau (Tau1) + ThqPartonHistory->auxdecor< float >("MC_nu_from_Tau1_m") = tH.nu_from_Tau1_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_nu_from_Tau1_pt") = tH.nu_from_Tau1_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_nu_from_Tau1_phi") = tH.nu_from_Tau1_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_nu_from_Tau1_pdgId") = tH.nu_from_Tau1_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_nu_from_Tau1_eta", tH.nu_from_Tau1_p4); + + //neutrino from second Tau (Tau2) + ThqPartonHistory->auxdecor< float >("MC_nu_from_Tau2_m") = tH.nu_from_Tau2_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_nu_from_Tau2_pt") = tH.nu_from_Tau2_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_nu_from_Tau2_phi") = tH.nu_from_Tau2_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_nu_from_Tau2_pdgId") = tH.nu_from_Tau2_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_nu_from_Tau2_eta", tH.nu_from_Tau2_p4); + + //First hadronic decay product from the W-Boson from Tau1 + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_Tau1_m") = tH.W_decay1_from_Tau1_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_Tau1_pt") = tH.W_decay1_from_Tau1_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_Tau1_phi") = tH.W_decay1_from_Tau1_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_Wdecay1_from_Tau1_pdgId") = tH.W_decay1_from_Tau1_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_Wdecay1_from_Tau1_eta", tH.W_decay1_from_Tau1_p4); + + //Second hadronic decay product from the W-Boson from Tau1 + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_Tau1_m") = tH.W_decay2_from_Tau1_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_Tau1_pt") = tH.W_decay2_from_Tau1_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_Tau1_phi") = tH.W_decay2_from_Tau1_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_Wdecay2_from_Tau1_pdgId") = tH.W_decay2_from_Tau1_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_Wdecay2_from_Tau1_eta", tH.W_decay2_from_Tau1_p4); + + //First leptonic decay product from the W-Boson from Tau2 + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_Tau2_m") = tH.W_decay1_from_Tau2_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_Tau2_pt") = tH.W_decay1_from_Tau2_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay1_from_Tau2_phi") = tH.W_decay1_from_Tau2_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_Wdecay1_from_Tau2_pdgId") = tH.W_decay1_from_Tau2_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_Wdecay1_from_Tau2_eta", tH.W_decay1_from_Tau2_p4); + + //Second leptonic decay product from the W-Boson from Tau2 + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_Tau2_m") = tH.W_decay2_from_Tau2_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_Tau2_pt") = tH.W_decay2_from_Tau2_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_Wdecay2_from_Tau2_phi") = tH.W_decay2_from_Tau2_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_Wdecay2_from_Tau2_pdgId") = tH.W_decay2_from_Tau2_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_Wdecay2_from_Tau2_eta", tH.W_decay2_from_Tau2_p4); + + //b quark + ThqPartonHistory->auxdecor< float >("MC_b_m") = tH.b_p4.M(); + ThqPartonHistory->auxdecor< float >("MC_b_pt") = tH.b_p4.Pt(); + ThqPartonHistory->auxdecor< float >("MC_b_phi") = tH.b_p4.Phi(); + ThqPartonHistory->auxdecor< int >("MC_b_pdgId") = tH.b_pdgId; + fillEtaBranch(ThqPartonHistory, "MC_b_eta", tH.b_p4); + + //Hadr. Variables + ThqPartonHistory->auxdecor< int >("MC_hadr_Tau_Jet1") = tH.TauJets1; + ThqPartonHistory->auxdecor< int >("MC_hadr_Tau_Jet2") = tH.TauJets2; + } + } + } + StatusCode CalcThqPartonHistory::execute() + { + //Get the Truth Particles + const xAOD::TruthParticleContainer* truthParticles(nullptr); + ATH_CHECK(evtStore()->retrieve(truthParticles, m_config->sgKeyMCParticle())); + + // Create the partonHistory xAOD object + xAOD::PartonHistoryAuxContainer* partonAuxCont = new xAOD::PartonHistoryAuxContainer{}; + xAOD::PartonHistoryContainer* partonCont = new xAOD::PartonHistoryContainer{}; + partonCont->setStore(partonAuxCont); + + xAOD::PartonHistory* ThqPartonHistory = new xAOD::PartonHistory{}; + partonCont->push_back(ThqPartonHistory); + + // Recover the parton history for TH events + THHistorySaver(truthParticles, ThqPartonHistory); + + // Save to StoreGate / TStore + std::string outputSGKey = m_config->sgKeyTopPartonHistory(); + std::string outputSGKeyAux = outputSGKey + "Aux."; + + xAOD::TReturnCode save = evtStore()->tds()->record(partonCont, outputSGKey); + xAOD::TReturnCode saveAux = evtStore()->tds()->record(partonAuxCont, outputSGKeyAux); + if (!save || !saveAux) { + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; + } +} diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopPartons/TopPartons/CalcThqPartonHistory.h b/PhysicsAnalysis/TopPhys/xAOD/TopPartons/TopPartons/CalcThqPartonHistory.h new file mode 100644 index 00000000000..3843ec9b079 --- /dev/null +++ b/PhysicsAnalysis/TopPhys/xAOD/TopPartons/TopPartons/CalcThqPartonHistory.h @@ -0,0 +1,80 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef ANALYSISTOP_TOPPARTONS_CALCThqPARTONHISTORY_H +#define ANALYSISTOP_TOPPARTONS_CALCThqARTONHISTORY_H + + +// Framework include(s): +#include "TopPartons/CalcTopPartonHistory.h" +#include "xAODTruth/TruthParticleContainer.h" +#include "TopPartons/PartonHistory.h" + +// forward declaration(s): +namespace top{ + class TopConfig; +} + +namespace top{ + + class CalcThqPartonHistory : public CalcTopPartonHistory{ + public: + explicit CalcThqPartonHistory( const std::string& name ); + virtual ~CalcThqPartonHistory() {} + + struct tH_values{ + //Higgs + TLorentzVector Higgs_p4; + TLorentzVector Tau1_from_Higgs_p4; + int Tau1_from_Higgs_pdgId; + TLorentzVector Tau2_from_Higgs_p4; + int Tau2_from_Higgs_pdgId; + TLorentzVector nu_from_Tau1_p4; + int nu_from_Tau1_pdgId; + TLorentzVector nu_from_Tau2_p4; + int nu_from_Tau2_pdgId; + TLorentzVector W_decay1_from_Tau1_p4; + int W_decay1_from_Tau1_pdgId; + TLorentzVector W_decay2_from_Tau1_p4; + int W_decay2_from_Tau1_pdgId; + TLorentzVector W_decay1_from_Tau2_p4; + int W_decay1_from_Tau2_pdgId; + TLorentzVector W_decay2_from_Tau2_p4; + int W_decay2_from_Tau2_pdgId; + + //Bools + int TauJets1; + int TauJets2; + + //b + TLorentzVector b_p4; + int b_pdgId; + + }tH; + //Storing parton history for ttbar resonance analysis + CalcThqPartonHistory(const CalcThqPartonHistory& rhs) = delete; + CalcThqPartonHistory(CalcThqPartonHistory&& rhs) = delete; + CalcThqPartonHistory& operator=(const CalcThqPartonHistory& rhs) = delete; + + void THHistorySaver(const xAOD::TruthParticleContainer * truthParticles, xAOD::PartonHistory * ThqPartonHistory); + + //handle gamma radiation of taus + const xAOD::TruthParticle* findAfterGamma(const xAOD::TruthParticle* particle); + + ///Store the four-momentum of several particles in the Higgs decay chain + bool Higgstautau( const xAOD::TruthParticleContainer* truthParticles, int start); + + //Store four-momentum of bottom quark + bool bottom(const xAOD::TruthParticleContainer* truthParticles, int start); + + int sign(int a); + + virtual StatusCode execute(); + + }; + + +} + +#endif -- GitLab From b99808f46692e1ddedb67370241a51c9f332ddd4 Mon Sep 17 00:00:00 2001 From: Nazim Huseynov Date: Wed, 2 Oct 2019 09:33:12 +0200 Subject: [PATCH 2/2] tH additional truth variables2 --- .../xAOD/TopAnalysis/util/top-xaod-v2.cxx | 5 - .../xAOD/TopAnalysis/util/top-xaod.cxx | 32 ++---- .../xAOD/TopPartons/Root/PartonHistory.cxx | 106 ++++++++++++++++++ .../TopPartons/TopPartons/PartonHistory.h | 1 + 4 files changed, 118 insertions(+), 26 deletions(-) diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod-v2.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod-v2.cxx index dd388433e98..b61f8934a43 100644 --- a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod-v2.cxx +++ b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod-v2.cxx @@ -727,11 +727,6 @@ void RunEventLoop(std::vector filenames, AnalysisTopTools analysisT // Get the event analysisTopTools.xaodEvent->getEntry(entry); - // Need to run mcEventWeight recalculation before processing non-reco events - // otherwise AnalysisTop_eventWeight is not available in EventInfo - if (topConfig->isMC()){ - analysisTopTools.topScaleFactorCalculator->mcEventWeight(); // Do not need the returned value - } // Need to run pileup calculation first as required inside truth event if (topConfig->doPileupReweighting() && !topConfig->isTruthDxAOD()) { top::check(analysisTopTools.topScaleFactorCalculator->executePileup(), "Failed to execute pileup reweighting"); diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod.cxx index d7c77d26910..0a11dc8a967 100644 --- a/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod.cxx +++ b/PhysicsAnalysis/TopPhys/xAOD/TopAnalysis/util/top-xaod.cxx @@ -59,6 +59,7 @@ #include "TopPartons/CalcTTZPartonHistory.h" #include "TopPartons/CalcTopPartonHistory.h" #include "TopPartons/CalcTtbarGammaPartonHistory.h" +#include "TopPartons/CalcThqPartonHistory.h" #include "TopParticleLevel/ParticleLevelLoader.h" @@ -332,6 +333,10 @@ int main(int argc, char** argv) { top::check(topPartonHistory->setProperty( "config" , topConfig ) , "Failed to setProperty of top::CalcTtbarGammaPartonHistory"); } + else if (settings->value("TopPartonHistory") == "tHq") { + topPartonHistory = std::unique_ptr(new top::CalcThqPartonHistory("top::CalcThqPartonHistory")); + top::check(topPartonHistory->setProperty("config", topConfig), "Failed to setProperty of top::CalcThqPartonHistory"); + } //LHAPDF SF calculation std::unique_ptr PDF_SF(nullptr); @@ -382,10 +387,10 @@ int main(int argc, char** argv) { } for (unsigned int mmi = 0; mmi("AsymptMatrixTool_"+mmi)); - top::check(topfakesMMWeightsIFF.back()->setProperty( "InputFiles", std::vector{FakesMMConfigIFF[mmi][0]} ) , "Failed To setProperty InputFiles of AsymptMatrixTool"); - top::check(topfakesMMWeightsIFF.back()->setProperty( "Selection", FakesMMConfigIFF[mmi][1] ),"Failed to set the selection FakesMMIFFConfigs for selection " + FakesMMConfigIFF[mmi][1]); - top::check(topfakesMMWeightsIFF.back()->setProperty( "Process", FakesMMConfigIFF[mmi][2] ),"Failed to set the selection FakesMMIFFConfigs for process " + FakesMMConfigIFF[mmi][2]); + top::check(topfakesMMWeightsIFF.back()->setProperty( "InputFiles", FakesMMConfigIFF[mmi][0] ) , "Failed To setProperty InputFiles of AsymptMatrixTool"); + top::check(topfakesMMWeightsIFF.back()->setProperty( "DefaultSelection", FakesMMConfigIFF[mmi][1] ),"Failed to set the selection FakesMMIFFConfigs for selection " + FakesMMConfigIFF[mmi][1]); top::check(topfakesMMWeightsIFF.back()->setProperty( "EnergyUnit", "GeV" ) , "Failed to setProperty EnergyUnit of AsymptMatrixTool"); + top::check(topfakesMMWeightsIFF.back()->setProperty( "SkipUncertainties", false ) , "Failed to setProperty SkipUncertainties of AsymptMatrixTool"); top::check(topfakesMMWeightsIFF.back()->setProperty( "ConvertWhenMissing", true ) , "Failed to setProperty ConvertWhenMissing of AsymptMatrixTool"); top::check(topfakesMMWeightsIFF.back()->setProperty( "TightDecoration", "passPreORSelection,as_char" ) , "Failed to setProperty TightDecoration of AsymptMatrixTool"); if(topConfig->FakesMMIFFDebug()) @@ -815,6 +820,7 @@ int main(int argc, char** argv) { } ///-- weights for matrix-method fakes estimate from IFF tools, only for nominal --/// if (!topConfig->isMC() && topConfig->doFakesMMWeightsIFF() && currentSystematic->hashValue() == topConfig->nominalHashValue()) { + //std::vector lepton; xAOD::IParticleContainer lepton(SG::VIEW_ELEMENTS); for (auto const & t : topEvent.m_electrons) lepton.push_back(static_cast(t)); @@ -834,30 +840,14 @@ int main(int argc, char** argv) { } std::vector mmweight; - std::vector > mmweight_var; - std::vector > mmweight_varname; for (unsigned int mmi = 0; mmiaddEvent(lepton, tightIFF) , "Failed to execute fakes mmweight IFF addEvent()"); top::check( topfakesMMWeightsIFF[mmi]->addEvent(lepton) , "Failed to execute fakes mmweight IFF addEvent()"); float asmWgt = 0.; - top::check( topfakesMMWeightsIFF[mmi]->applySystematicVariation({}) , "Failed to execute fakes mmweight IFF applySystematicVariation()"); top::check( topfakesMMWeightsIFF[mmi]->getEventWeight(asmWgt, FakesMMConfigIFF[mmi][1], FakesMMConfigIFF[mmi][2]) , "Failed to execute fakes mmweight IFF getEventWeight()"); - - std::vector asmWgt_var(topfakesMMWeightsIFF[mmi]->affectingSystematics().size()); - std::vector asmWgt_varname(topfakesMMWeightsIFF[mmi]->affectingSystematics().size()); - for(const auto& sysvar : topfakesMMWeightsIFF[mmi]->affectingSystematics()) { - float mmweight_syst; - top::check( topfakesMMWeightsIFF[mmi]->applySystematicVariation({sysvar}) , "Failed to execute fakes mmweight IFF applySystematicVariation()"); - top::check( topfakesMMWeightsIFF[mmi]->getEventWeight(mmweight_syst, FakesMMConfigIFF[mmi][1], FakesMMConfigIFF[mmi][2]) , "Failed to execute fakes mmweight IFF getEventWeight()"); - asmWgt_var.push_back(mmweight_syst); - asmWgt_varname.push_back(sysvar.name()); - } mmweight.push_back(asmWgt); - mmweight_var.push_back(asmWgt_var); - mmweight_varname.push_back(asmWgt_varname); } - topEvent.m_info->auxdecor >("ASMWeight") = mmweight; - topEvent.m_info->auxdecor > >("ASMWeight_Syst") = mmweight_var; - topEvent.m_info->auxdecor > >("ASMWeight_Systname") = mmweight_varname; + topEvent.m_info->auxdecor >("MMWeight_IFF") = mmweight; } ///-- Save event - we defer to eventSaver the decision to write or not --/// eventSaver->saveEvent(topEvent); diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopPartons/Root/PartonHistory.cxx b/PhysicsAnalysis/TopPhys/xAOD/TopPartons/Root/PartonHistory.cxx index fc45c614c6a..04344cd19a7 100644 --- a/PhysicsAnalysis/TopPhys/xAOD/TopPartons/Root/PartonHistory.cxx +++ b/PhysicsAnalysis/TopPhys/xAOD/TopPartons/Root/PartonHistory.cxx @@ -480,6 +480,112 @@ namespace xAOD{ } + // Initialize variables for thq events + void PartonHistory::IniVarThq() { + //t variables + this->auxdecor< float >("MC_t_beforeFSR_m") = -1000; + this->auxdecor< float >("MC_t_beforeFSR_pt") = -1000; + this->auxdecor< float >("MC_t_beforeFSR_eta") = -1000; + this->auxdecor< float >("MC_t_beforeFSR_phi") = -1000; + + this->auxdecor< float >("MC_t_afterFSR_m") = -1000; + this->auxdecor< float >("MC_t_afterFSR_pt") = -1000; + this->auxdecor< float >("MC_t_afterFSR_eta") = -1000; + this->auxdecor< float >("MC_t_afterFSR_phi") = -1000; + + this->auxdecor< float >("MC_t_afterFSR_SC_m") = -1000; + this->auxdecor< float >("MC_t_afterFSR_SC_pt") = -1000; + this->auxdecor< float >("MC_t_afterFSR_SC_eta") = -1000; + this->auxdecor< float >("MC_t_afterFSR_SC_phi") = -1000; + + this->auxdecor< float >("MC_W_from_t_m") = -1000; + this->auxdecor< float >("MC_W_from_t_pt") = -1000; + this->auxdecor< float >("MC_W_from_t_eta") = -1000; + this->auxdecor< float >("MC_W_from_t_phi") = -1000; + + this->auxdecor< float >("MC_b_from_t_m") = -1000; + this->auxdecor< float >("MC_b_from_t_pt") = -1000; + this->auxdecor< float >("MC_b_from_t_eta") = -1000; + this->auxdecor< float >("MC_b_from_t_phi") = -1000; + + this->auxdecor< float >("MC_Wdecay1_from_t_m") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_t_pt") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_t_eta") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_t_phi") = -1000; + this->auxdecor< int >("MC_Wdecay1_from_t_pdgId") = -9999; + + this->auxdecor< float >("MC_Wdecay2_from_t_m") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_t_pt") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_t_eta") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_t_phi") = -1000; + this->auxdecor< int >("MC_Wdecay2_from_t_pdgId") = -9999; + + // Higgs variables + this->auxdecor< float >("MC_Higgs_m") = -1000; + this->auxdecor< float >("MC_Higgs_pt") = -1000; + this->auxdecor< float >("MC_Higgs_eta") = -1000; + this->auxdecor< float >("MC_Higgs_phi") = -1000; + + this->auxdecor< float >("MC_Tau1_from_Higgs_m") = -1000; + this->auxdecor< float >("MC_Tau1_from_Higgs_pt") = -1000; + this->auxdecor< float >("MC_Tau1_from_Higgs_eta") = -1000; + this->auxdecor< float >("MC_Tau1_from_Higgs_phi") = -1000; + this->auxdecor< int >("MC_Tau1_from_Higgs_pdgId") = -9999; + + this->auxdecor< float >("MC_Tau2_from_Higgs_m") = -1000; + this->auxdecor< float >("MC_Tau2_from_Higgs_pt") = -1000; + this->auxdecor< float >("MC_Tau2_from_Higgs_eta") = -1000; + this->auxdecor< float >("MC_Tau2_from_Higgs_phi") = -1000; + this->auxdecor< int >("MC_Tau2_from_Higgs_pdgId") = -9999; + + + this->auxdecor< float >("MC_nu_from_Tau1_m") = -1000; + this->auxdecor< float >("MC_nu_from_Tau1_pt") = -1000; + this->auxdecor< float >("MC_nu_from_Tau1_eta") = -1000; + this->auxdecor< float >("MC_nu_from_Tau1_phi") = -1000; + this->auxdecor< int >("MC_nu_from_Tau1_pdgId") = -9999; + + this->auxdecor< float >("MC_nu_from_Tau2_m") = -1000; + this->auxdecor< float >("MC_nu_from_Tau2_pt") = -1000; + this->auxdecor< float >("MC_nu_from_Tau2_eta") = -1000; + this->auxdecor< float >("MC_nu_from_Tau2_phi") = -1000; + this->auxdecor< int >("MC_nu_from_Tau2_pdgId") = -9999; + + this->auxdecor< float >("MC_Wdecay1_from_Tau1_m") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_Tau1_pt") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_Tau1_eta") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_Tau1_phi") = -1000; + this->auxdecor< int >("MC_Wdecay1_from_Tau1_pdgId") = -9999; + + this->auxdecor< float >("MC_Wdecay2_from_Tau1_m") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_Tau1_pt") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_Tau1_eta") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_Tau1_phi") = -1000; + this->auxdecor< int >("MC_Wdecay2_from_Tau1_pdgId") = -9999; + + this->auxdecor< float >("MC_Wdecay1_from_Tau2_m") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_Tau2_pt") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_Tau2_eta") = -1000; + this->auxdecor< float >("MC_Wdecay1_from_Tau2_phi") = -1000; + this->auxdecor< int >("MC_Wdecay1_from_Tau2_pdgId") = -9999; + + this->auxdecor< float >("MC_Wdecay2_from_Tau2_m") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_Tau2_pt") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_Tau2_eta") = -1000; + this->auxdecor< float >("MC_Wdecay2_from_Tau2_phi") = -1000; + this->auxdecor< int >("MC_Wdecay2_from_Tau2_pdgId") = -9999; + + this->auxdecor< float >("MC_b_m") = -1000; + this->auxdecor< float >("MC_b_pt") = -1000; + this->auxdecor< float >("MC_b_eta") = -1000; + this->auxdecor< float >("MC_b_phi") = -1000; + this->auxdecor< int >("MC_b_pdgId") = -9999; + + this->auxdecor< int >("MC_hadr_Tau_Jet1") = 0; + this->auxdecor< int >("MC_hadr_Tau_Jet2") = 0; + + + } } diff --git a/PhysicsAnalysis/TopPhys/xAOD/TopPartons/TopPartons/PartonHistory.h b/PhysicsAnalysis/TopPhys/xAOD/TopPartons/TopPartons/PartonHistory.h index 014dbf038f3..9c3672f9e69 100644 --- a/PhysicsAnalysis/TopPhys/xAOD/TopPartons/TopPartons/PartonHistory.h +++ b/PhysicsAnalysis/TopPhys/xAOD/TopPartons/TopPartons/PartonHistory.h @@ -60,6 +60,7 @@ namespace xAOD{ void IniVarWtb(); void IniVarZ(); void IniVarTtGamma(); + void IniVarThq(); }; typedef DataVector < xAOD::PartonHistory > PartonHistoryContainer; -- GitLab