diff --git a/interface/Defines.h b/interface/Defines.h index ece4138..033d798 100644 --- a/interface/Defines.h +++ b/interface/Defines.h @@ -2,10 +2,10 @@ // Debugging macros -// #define _TT_DEBUG_ -#define TT_MTT_DEBUG (false) -#define TT_HLT_DEBUG (false) -#define TT_GEN_DEBUG (false) +//#define _TT_DEBUG_ +//#define TT_MTT_DEBUG +//#define TT_HLT_DEBUG +//#define TT_GEN_DEBUG #if TT_GEN_DEBUG diff --git a/interface/Indices.h b/interface/Indices.h index da47f99..4740ac3 100644 --- a/interface/Indices.h +++ b/interface/Indices.h @@ -46,35 +46,39 @@ namespace TTAnalysis { } // Combination of Jet IDs for two jets (NOTE: NOT USED FOR NOW) - uint16_t JetJetID(const JetID::JetID& id1, const JetID::JetID& id2); - std::string JetJetIDStr(const JetID::JetID& id1, const JetID::JetID& id2); + //uint16_t JetJetID(const JetID::JetID& id1, const JetID::JetID& id2); + //std::string JetJetIDStr(const JetID::JetID& id1, const JetID::JetID& id2); // B-tagging working points namespace BWP { - enum BWP{ L, M, T, Count }; - const std::array it = {{ L, M, T }}; - const std::map map = { {L, "L"}, {M, "M"}, {T, "T"} }; + enum BWP{ A, L, M, Count }; + const std::array it = {{ A, L, M }}; + const std::map map = { {A, "A"}, {L, "L"}, {M, "M"} }; } // Combination of Jet ID and B-tagging working point (NOTE: NOT USED FOR NOW) - uint16_t JetIDBWP(const JetID::JetID& id, const BWP::BWP& wp); - std::string JetIDBWPStr(const JetID::JetID& id, const BWP::BWP& wp); + //uint16_t JetIDBWP(const JetID::JetID& id, const BWP::BWP& wp); + //std::string JetIDBWPStr(const JetID::JetID& id, const BWP::BWP& wp); - // Combination of Lepton ID + Lepton Isolation (one lepton) and B-tagging working point for one jet - uint16_t LepIDIsoJetBWP(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp); - std::string LepIDIsoJetBWPStr(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp); + // Combination of Lepton ID + Lepton Isolation (on lepton) and B-tagging working point for one jet (NOTE: NOT USED FOR NOW) + //uint16_t LepIDIsoJetBWP(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp); + //std::string LepIDIsoJetBWPStr(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp); + + // Combination of Lepton ID + Lepton Isolation (two leptons) and B-tagging working point for one jet + uint16_t LepLepIDIsoJetBWP(const LepID::LepID& id1, const LepIso::LepIso& iso1, const LepID::LepID& id2, const LepIso::LepIso& iso2, const BWP::BWP& wp); + std::string LepLepIDIsoJetBWPStr(const LepID::LepID& id1, const LepIso::LepIso& iso1, const LepID::LepID& id2, const LepIso::LepIso& iso2, const BWP::BWP& wp); // Combination of B-tagging working points for two jets uint16_t JetJetBWP(const BWP::BWP& wp1, const BWP::BWP& wp2); std::string JetJetBWPStr(const BWP::BWP& wp1, const BWP::BWP& wp2); // Combination of Jet ID and B-tagging working points for two jets (NOTE: NOT USED FOR NOW) - uint16_t JetJetIDBWP(const JetID::JetID& id1, const BWP::BWP& wp1, const JetID::JetID& id2, const BWP::BWP& wp2); - std::string JetJetIDBWPStr(const JetID::JetID& id1, const BWP::BWP wp1, const JetID::JetID& id2, const BWP::BWP wp2); + //uint16_t JetJetIDBWP(const JetID::JetID& id1, const BWP::BWP& wp1, const JetID::JetID& id2, const BWP::BWP& wp2); + //std::string JetJetIDBWPStr(const JetID::JetID& id1, const BWP::BWP wp1, const JetID::JetID& id2, const BWP::BWP wp2); - // Combination of Lepton ID + Lepton Isolation (one lepton) and B-tagging working points for two jets - uint16_t LepIDIsoJetJetBWP(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp1, const BWP::BWP& wp2); - std::string LepIDIsoJetJetBWPStr(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp1, const BWP::BWP& wp2); + // Combination of Lepton ID + Lepton Isolation (one lepton) and B-tagging working points for two jets (NOTE: NOT USED FOR NOW) + //uint16_t LepIDIsoJetJetBWP(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp1, const BWP::BWP& wp2); + //std::string LepIDIsoJetJetBWPStr(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp1, const BWP::BWP& wp2); // Combination of Lepton ID, Lepton Isolation, and B-tagging working points for a two-lepton-two-b-jets object uint16_t LepLepIDIsoJetJetBWP(const LepID::LepID& id1, const LepIso::LepIso& iso1, const LepID::LepID& id2, const LepIso::LepIso& iso2, const BWP::BWP& wp1, const BWP::BWP& wp2); diff --git a/interface/NeutrinosSolver.h b/interface/NeutrinosSolver.h index b60f815..86f931b 100644 --- a/interface/NeutrinosSolver.h +++ b/interface/NeutrinosSolver.h @@ -24,8 +24,7 @@ class NeutrinosSolver { public: using LorentzVector = ROOT::Math::LorentzVector>; - NeutrinosSolver(float top_mass, float w_mass): - t_mass(top_mass), w_mass(w_mass) { + NeutrinosSolver(){ // Empty } @@ -33,9 +32,7 @@ class NeutrinosSolver { const LorentzVector& lepton2_p4, const LorentzVector& bjet1_p4, const LorentzVector& bjet2_p4, - const LorentzVector& met); - - private: - float t_mass = 172.5; - float w_mass = 80.4; + const LorentzVector& met, + const double wplus_mass=80.4, const double wminus_mass=80.4, + const double t_mass=172.5, const double tbar_mass=172.5); }; diff --git a/interface/TTAnalyzer.h b/interface/TTAnalyzer.h index 6a3cfe0..21a9798 100644 --- a/interface/TTAnalyzer.h +++ b/interface/TTAnalyzer.h @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -24,7 +25,8 @@ class TTAnalyzer: public Framework::Analyzer { m_met_producer(config.getParameter("metProducer")), m_electronPtCut( config.getUntrackedParameter("electronPtCut", 20) ), - m_electronEtaCut( config.getUntrackedParameter("electronEtaCut", 2.5) ), + m_electronEtaCut( config.getUntrackedParameter("electronEtaCut", 2.4) ), + m_electronRemoveGap( config.getUntrackedParameter("electronRemoveGap", true) ), m_electronVetoIDName( config.getUntrackedParameter("electronVetoIDName") ), m_electronLooseIDName( config.getUntrackedParameter("electronLooseIDName") ), m_electronMediumIDName( config.getUntrackedParameter("electronMediumIDName") ), @@ -36,9 +38,9 @@ class TTAnalyzer: public Framework::Analyzer { m_muonTightIsoCut( config.getUntrackedParameter("muonTightIsoCut", 0.12) ), m_jetPtCut( config.getUntrackedParameter("jetPtCut", 30) ), - m_jetEtaCut( config.getUntrackedParameter("jetEtaCut", 2.5) ), + m_jetEtaCut( config.getUntrackedParameter("jetEtaCut", 2.4) ), m_bJetEtaCut( config.getUntrackedParameter("bJetEtaCut", 2.4) ), - m_jetPUID( config.getUntrackedParameter("jetPUID", std::numeric_limits::min()) ), + m_jetPUID( config.getUntrackedParameter("jetPUID", std::numeric_limits::min()) ), m_jetDRleptonCut( config.getUntrackedParameter("jetDRleptonCut", 0.3) ), m_jetID( config.getUntrackedParameter("jetID", "loose") ), m_jetCSVv2Name( config.getUntrackedParameter("jetCSVv2Name", "pfCombinedInclusiveSecondaryVertexV2BJetTags") ), @@ -46,9 +48,17 @@ class TTAnalyzer: public Framework::Analyzer { m_jetCSVv2M( config.getUntrackedParameter("jetCSVv2M", 0.89) ), m_jetCSVv2T( config.getUntrackedParameter("jetCSVv2T", 0.97) ), - m_hltDRCut( config.getUntrackedParameter("hltDRCut", std::numeric_limits::max()) ), - m_hltDPtCut( config.getUntrackedParameter("hltDPtCut", std::numeric_limits::max()) ) + m_hltDRCut( config.getUntrackedParameter("hltDRCut", std::numeric_limits::max()) ), + m_hltDPtCut( config.getUntrackedParameter("hltDPtCut", std::numeric_limits::max()) ) { + if(config.exists("hltScaleFactors")){ + const edm::ParameterSet& hltSFPSet = config.getUntrackedParameter("hltScaleFactors"); + std::vector hltSFNames = hltSFPSet.getParameterNames(); + for(const std::string& sf: hltSFNames){ + BinnedValuesJSONParser parser(hltSFPSet.getUntrackedParameter(sf).fullPath()); + m_hltSF.emplace(sf, std::move(parser.get_values())); + } + } } virtual void analyze(const edm::Event&, const edm::EventSetup&, const ProducersManager&, const AnalyzersManager&, const CategoryManager&) override; @@ -64,33 +74,24 @@ class TTAnalyzer: public Framework::Analyzer { BRANCH(diLeptons_IDIso, std::vector>); BRANCH(selJets, std::vector); - BRANCH(selJets_selID, std::vector); - // ex.: selectedJets_..._DRCut[X][0] is the highest Pt selected jet with minDRjl>0.3 taking into account ID/Iso-X Leptons - BRANCH(selJets_selID_DRCut, std::vector>); - // ex.: selectedBJets_..._PtOrdered[X][0] is the highest Pt selected jet with minDRjl>0.3 taking into account ID/Iso/Btag-X combination - BRANCH(selBJets_DRCut_BWP_PtOrdered, std::vector>); - BRANCH(selBJets_DRCut_BWP_CSVv2Ordered, std::vector>); + // ex.: selJets_..._PtOrdered[X][0] is the highest Pt selected jet with minDRjl>0.3 taking into account leptonID/Iso pair & Btag-X combination + BRANCH(selJets_DRCut_BWP_PtOrdered, std::vector>); + BRANCH(selJets_DRCut_BWP_CSVv2Ordered, std::vector>); BRANCH(diJets, std::vector); - // ex.: diJets_DRCut[X][0] is first diJet with minDRjl>0.3 taking into account ID/Iso-X Leptons - BRANCH(diJets_DRCut, std::vector>); - // ex.: diBJets_..._CSVv2Ordered[X][0] is the b-jet pair with highest CSVv2 values and with minDRjl>0.3 taking into account the leptonID/Iso/Btag-X combination - BRANCH(diBJets_DRCut_BWP_PtOrdered, std::vector>); - BRANCH(diBJets_DRCut_BWP_CSVv2Ordered, std::vector>); + // ex.: diJets_..._CSVv2Ordered[X][0] is the jet pair with highest CSVv2 values and with minDRjl>0.3 taking into account the leptonID/Iso pair & Btag-X combination + BRANCH(diJets_DRCut_BWP_PtOrdered, std::vector>); + BRANCH(diJets_DRCut_BWP_CSVv2Ordered, std::vector>); // For all the following: indices are combinations of LeptonID/LeptonIso/(B-tagging working point) - + BRANCH(diLepDiJets, std::vector); + BRANCH(diLepDiJets_DRCut_BWP_PtOrdered, std::vector>); // di-leptons of combined ID/Iso with di-jets built out of jets having minDRjl>cut taking into account lepton ID/Iso corresponding to the chosen combination of the two leptons of the object + BRANCH(diLepDiJets_DRCut_BWP_CSVv2Ordered, std::vector>); - BRANCH(diLepDiJets_DRCut, std::vector>); // di-leptons of combined ID/Iso with di-jets built out of jets having minDRjl>cut taking into account lepton ID/Iso corresponding to the loosest combination of the two leptons of the object - BRANCH(diLepDiBJets_DRCut_BWP_PtOrdered, std::vector>); - BRANCH(diLepDiBJets_DRCut_BWP_CSVv2Ordered, std::vector>); - BRANCH(diLepDiJetsMet, std::vector); - - BRANCH(diLepDiJetsMet_DRCut, std::vector>); - BRANCH(diLepDiBJetsMet_DRCut_BWP_PtOrdered, std::vector>); - BRANCH(diLepDiBJetsMet_DRCut_BWP_CSVv2Ordered, std::vector>); + BRANCH(diLepDiJetsMet_DRCut_BWP_PtOrdered, std::vector>); // di-leptons of combined ID/Iso with di-jets built out of jets having minDRjl>cut taking into account lepton ID/Iso corresponding to the chosen combination of the two leptons of the object + BRANCH(diLepDiJetsMet_DRCut_BWP_CSVv2Ordered, std::vector>); BRANCH(ttbar, std::vector>>); @@ -140,28 +141,32 @@ class TTAnalyzer: public Framework::Analyzer { BRANCH(gen_b_lepton_t_deltaR, float); // DeltaR between the b quark and the lepton coming from the top decay chain BRANCH(gen_bbar_lepton_tbar_deltaR, float); // DeltaR between the b quark and the lepton coming from the top decay chain - // These two vectors are indexed wrt LepLepId enum - BRANCH(gen_b_deltaR, std::vector>); // DeltaR between the gen b coming from the top decay and each selected jets. Indexed as `selectedJets_tightID_DRcut` array - BRANCH(gen_bbar_deltaR, std::vector>); // DeltaR between the gen bbar coming from the anti-top decay chain and each selected jets. Indexed as `selectedJets_tightID_DRcut` array - - // These two vectors are indexed wrt LepLepId enum - BRANCH(gen_b_beforeFSR_deltaR, std::vector>); // DeltaR between the gen b coming from the top decay and each selected jets. Indexed as `selectedJets_tightID_DRcut` array - BRANCH(gen_bbar_beforeFSR_deltaR, std::vector>); // DeltaR between the gen bbar coming from the anti-top decay chain and each selected jets. Indexed as `selectedJets_tightID_DRcut` array - BRANCH(gen_lepton_t_deltaR, std::vector); // DeltaR between the gen lepton coming from the top decay chain and each selected lepton. Indexed as `leptons` array BRANCH(gen_lepton_tbar_deltaR, std::vector); // DeltaR between the gen lepton coming from the anti-top decay chain and each selected lepton. Indexed as `leptons` array - // These two vectors are indexed wrt LepLepId enum - BRANCH(gen_matched_b, std::vector); // Index inside the `selectedJets_tightID_DRcut` collection of the jet with the smallest deltaR with the gen b coming from the top decay - BRANCH(gen_matched_bbar, std::vector); // Index inside the `selectedJets_tightID_DRcut` collection of the jet with the smallest deltaR with the gen bbar coming from the anti-top decay - - // These two vectors are indexed wrt LepLepId enum - BRANCH(gen_matched_b_beforeFSR, std::vector); // Index inside the `selectedJets_tightID_DRcut` collection of the jet with the smallest deltaR with the gen b coming from the top decay - BRANCH(gen_matched_bbar_beforeFSR, std::vector); // Index inside the `selectedJets_tightID_DRcut` collection of the jet with the smallest deltaR with the gen bbar coming from the anti-top decay + BRANCH(gen_lepton_beforeFSR_t_deltaR, std::vector); // DeltaR between the gen lepton coming from the top decay chain and each selected lepton. Indexed as `leptons` array + BRANCH(gen_lepton_beforeFSR_tbar_deltaR, std::vector); // DeltaR between the gen lepton coming from the anti-top decay chain and each selected lepton. Indexed as `leptons` array BRANCH(gen_matched_lepton_t, int16_t); // Index inside the `leptons` collection of the lepton with the smallest deltaR with the gen lepton coming from the top decay chain BRANCH(gen_matched_lepton_tbar, int16_t); // Index inside the `leptons` collection of the lepton with the smallest deltaR with the gen lepton coming from the anti-top decay chain + BRANCH(gen_matched_lepton_beforeFSR_t, int16_t); // Index inside the `leptons` collection of the lepton with the smallest deltaR with the gen lepton coming from the top decay chain + BRANCH(gen_matched_lepton_beforeFSR_tbar, int16_t); // Index inside the `leptons` collection of the lepton with the smallest deltaR with the gen lepton coming from the anti-top decay chain + + // The following vectors are indexed wrt. LepLepIdJetJetBWP enum (as, for instance, selJets_DRCut_BWP_PtOrdered), giving indices pointing to the `selJets` collection + BRANCH(gen_b_deltaR, std::vector>); // DeltaR between the gen b coming from the top decay and each selected jets. + BRANCH(gen_bbar_deltaR, std::vector>); // DeltaR between the gen bbar coming from the anti-top decay chain and each selected jets. + + BRANCH(gen_b_beforeFSR_deltaR, std::vector>); // DeltaR between the gen b coming from the top decay and each selected jets. + BRANCH(gen_bbar_beforeFSR_deltaR, std::vector>); // DeltaR between the gen bbar coming from the anti-top decay chain and each selected jets. + + // The following vectors are indexed wrt. LepLepIdJetJetBWP enum (as, for instance, selJets_DRCut_BWP_PtOrdered) + BRANCH(gen_matched_b, std::vector); // Index inside the `selJets` collection of the jet with the smallest deltaR with the gen b coming from the top decay + BRANCH(gen_matched_bbar, std::vector); // Index inside the `selJets` collection of the jet with the smallest deltaR with the gen bbar coming from the anti-top decay + + BRANCH(gen_matched_b_beforeFSR, std::vector); // Index inside the `selJets` collection of the jet with the smallest deltaR with the gen b coming from the top decay + BRANCH(gen_matched_bbar_beforeFSR, std::vector); // Index inside the `selJets` collection of the jet with the smallest deltaR with the gen bbar coming from the anti-top decay + private: // Producers name @@ -171,6 +176,7 @@ class TTAnalyzer: public Framework::Analyzer { const std::string m_met_producer; const float m_electronPtCut, m_electronEtaCut; + const bool m_electronRemoveGap; const std::string m_electronVetoIDName; const std::string m_electronLooseIDName; const std::string m_electronMediumIDName; @@ -183,6 +189,7 @@ class TTAnalyzer: public Framework::Analyzer { const float m_jetCSVv2L, m_jetCSVv2M, m_jetCSVv2T; const float m_hltDRCut, m_hltDPtCut; + std::map m_hltSF; std::shared_ptr m_neutrinos_solver; diff --git a/interface/TTDileptonCategories.h b/interface/TTDileptonCategories.h index fa67522..4714e3b 100644 --- a/interface/TTDileptonCategories.h +++ b/interface/TTDileptonCategories.h @@ -14,8 +14,8 @@ class DileptonCategory: public Category { virtual void configure(const edm::ParameterSet& conf) override { m_MllCutSF = conf.getUntrackedParameter("MllCutSF", 20); m_MllCutDF = conf.getUntrackedParameter("MllCutDF", 20); - m_MllZVetoCutLow = conf.getUntrackedParameter("MllZVetoCutLow", 86); - m_MllZVetoCutHigh = conf.getUntrackedParameter("MllZVetoCutHigh", 116); + m_MllZVetoCutLow = conf.getUntrackedParameter("MllZVetoCutLow", 76); + m_MllZVetoCutHigh = conf.getUntrackedParameter("MllZVetoCutHigh", 106); m_HLTDoubleMuon = conf.getUntrackedParameter>("HLTDoubleMuon"); m_HLTDoubleEG = conf.getUntrackedParameter>("HLTDoubleEG"); m_HLTMuonEG = conf.getUntrackedParameter>("HLTMuonEG"); @@ -31,6 +31,7 @@ class DileptonCategory: public Category { DileptonCategory(): baseStrCategory("Category"), baseStrExtraDiLeptonVeto("ExtraDiLeptonVeto"), + baseStrDiLeptonTriggerBit("DiLeptonTriggerBit"), baseStrDiLeptonTriggerMatch("DiLeptonTriggerMatch"), baseStrMllCut("Mll"), baseStrMllZVetoCut("MllZVeto"), @@ -46,6 +47,7 @@ class DileptonCategory: public Category { std::string baseStrCategory; std::string baseStrExtraDiLeptonVeto; + std::string baseStrDiLeptonTriggerBit; std::string baseStrDiLeptonTriggerMatch; std::string baseStrMllCut; std::string baseStrMllZVetoCut; @@ -56,30 +58,39 @@ class DileptonCategory: public Category { std::vector m_HLTMuonEGRegex; enum class HLT { DoubleMuon, DoubleEG, MuonEG }; - - // Check that the hlt objects at indices hltIdx1, hltIdx2 have fired at least one and the same - // of the trigger paths in the group specified by pathGroup. - bool checkHLT(const HLTProducer& hlt, uint16_t hltIdx1, uint16_t hltIdx2, HLT pathGroup) const { - const std::vector* chosenPathGroup(nullptr); - + const std::vector* getPathGroup(HLT pathGroup) const { switch(pathGroup){ - case HLT::DoubleMuon: - chosenPathGroup = &m_HLTDoubleMuonRegex; - break; + return &m_HLTDoubleMuonRegex; case HLT::DoubleEG: - chosenPathGroup = &m_HLTDoubleEGRegex; - break; + return &m_HLTDoubleEGRegex; case HLT::MuonEG: - chosenPathGroup = &m_HLTMuonEGRegex; - break; - - default: - break; + return &m_HLTMuonEGRegex; + } + return nullptr; + }; + + // Simply check that a trigger specified in pathGroup has fired + bool checkHLT(const HLTProducer& hlt, HLT pathGroup) const { + const std::vector* chosenPathGroup = getPathGroup(pathGroup); + for(const auto& checkPath: *chosenPathGroup){ + for(const auto& firedPath: hlt.paths){ + if( boost::regex_match(firedPath, checkPath) ) + return true; + } } + + return false; + } + + // Check that the hlt objects at indices hltIdx1, hltIdx2 have fired at least one and the same + // of the trigger paths in the group specified by pathGroup. + bool checkHLT(const HLTProducer& hlt, uint16_t hltIdx1, uint16_t hltIdx2, HLT pathGroup) const { + + const std::vector* chosenPathGroup = getPathGroup(pathGroup); if( !chosenPathGroup || hltIdx1 >= hlt.object_paths.size() || hltIdx2 >= hlt.object_paths.size() ) return false; diff --git a/interface/Types.h b/interface/Types.h index 24a84ab..e5f8bf1 100644 --- a/interface/Types.h +++ b/interface/Types.h @@ -89,7 +89,8 @@ namespace TTAnalysis { struct DiLepton: BaseObject { DiLepton(): ID(LepID::Count*LepID::Count, false), - iso(LepIso::Count*LepIso::Count, false) + iso(LepIso::Count*LepIso::Count, false), + hlt_SF({1,1,1}) {} std::pair idxs; // stores indices to electron/muon arrays @@ -100,6 +101,7 @@ namespace TTAnalysis { bool isSF; // same flavour std::vector ID; // combination of two lepton IDs std::vector iso; // combination of two lepton isolations + std::vector hlt_SF; // scale factors for the corresponding trigger -> nominal/down/up float DR; float DEta; float DPhi; @@ -108,26 +110,26 @@ namespace TTAnalysis { struct Jet: BaseObject { Jet(): ID(JetID::Count, false), - minDRjl_lepIDIso(LepID::Count*LepIso::Count, std::numeric_limits::max()), + minDRjl_lepLepIDIso(LepID::Count*LepIso::Count*LepID::Count*LepIso::Count, std::numeric_limits::max()), BWP(BWP::Count, false) {} uint16_t idx; // index to jet array std::vector ID; - std::vector minDRjl_lepIDIso; // defined for each combination of a lepton ID and isolation + std::vector minDRjl_lepLepIDIso; // defined for each combination of two leptons' ID and isolation float CSVv2; std::vector BWP; }; struct DiJet: BaseObject { DiJet(): - minDRjl_lepIDIso(LepID::Count*LepIso::Count, std::numeric_limits::max()), + minDRjl_lepLepIDIso(LepID::Count*LepIso::Count*LepID::Count*LepIso::Count, std::numeric_limits::max()), BWP(BWP::Count*BWP::Count, false) {} std::pair idxs; // stores indices to jets array std::pair jidxs; // stores indices to TTAnalysis::Jet array - std::vector minDRjl_lepIDIso; // defined for each combination of a lepton ID and isolation + std::vector minDRjl_lepLepIDIso; // defined for each combination of two leptons' ID and isolation std::vector BWP; // combination of two b-tagging working points float DR; float DEta; @@ -141,10 +143,8 @@ namespace TTAnalysis { {} DiLepDiJet(const DiLepton& _diLepton, const int diLepIdx, const DiJet& _diJet, const int diJetIdx): BaseObject(_diLepton.p4 + _diJet.p4), - diLepton(&_diLepton), - diLepIdx(diLepIdx), - diJet(&_diJet), - diJetIdx(diJetIdx), + diLepton(&_diLepton), diLepIdx(diLepIdx), lidxs(_diLepton.lidxs), fwk_lidxs(_diLepton.idxs), + diJet(&_diJet), diJetIdx(diJetIdx), jidxs(_diJet.jidxs), fwk_jidxs(_diJet.idxs), DR_ll_jj( ROOT::Math::VectorUtil::DeltaR(_diLepton.p4, _diJet.p4) ), DEta_ll_jj( DeltaEta(_diLepton.p4, _diJet.p4) ), DPhi_ll_jj( ROOT::Math::VectorUtil::DeltaPhi(_diLepton.p4, _diJet.p4) ) @@ -152,8 +152,12 @@ namespace TTAnalysis { const DiLepton* diLepton; uint16_t diLepIdx; + std::pair lidxs; + std::pair fwk_lidxs; const DiJet* diJet; uint16_t diJetIdx; + std::pair jidxs; + std::pair fwk_jidxs; float DR_ll_jj, DEta_ll_jj, DPhi_ll_jj; @@ -164,10 +168,9 @@ namespace TTAnalysis { struct DiLepDiJetMet: DiLepDiJet { DiLepDiJetMet() {} - DiLepDiJetMet(const DiLepDiJet& diLepDiJet, uint16_t diLepDiJetIdx, const myLorentzVector& MetP4, bool hasNoHFMet = false): + DiLepDiJetMet(const DiLepDiJet& diLepDiJet, uint16_t diLepDiJetIdx, const myLorentzVector& MetP4): DiLepDiJet(*diLepDiJet.diLepton, diLepDiJet.diLepIdx, *diLepDiJet.diJet, diLepDiJet.diJetIdx), - diLepDiJetIdx(diLepDiJetIdx), - hasNoHFMet(hasNoHFMet) + diLepDiJetIdx(diLepDiJetIdx) { DiLepDiJet::minDRjl = diLepDiJet.minDRjl; DiLepDiJet::maxDRjl = diLepDiJet.maxDRjl; @@ -193,7 +196,6 @@ namespace TTAnalysis { } uint16_t diLepDiJetIdx; - bool hasNoHFMet; float DR_ll_Met, DR_jj_Met; float DEta_ll_Met, DEta_jj_Met; diff --git a/plugins/Indices.cc b/plugins/Indices.cc index 9adf1cd..1c48df9 100644 --- a/plugins/Indices.cc +++ b/plugins/Indices.cc @@ -37,29 +37,42 @@ namespace TTAnalysis { } // Combination of Jet IDs for two jets (NOTE: NOT USED FOR NOW) - uint16_t JetJetID(const JetID::JetID& id1, const JetID::JetID& id2){ + /*uint16_t JetJetID(const JetID::JetID& id1, const JetID::JetID& id2){ return JetID::Count * id1 + id2; } std::string JetJetIDStr(const JetID::JetID& id1, const JetID::JetID& id2){ return JetID::map.at(id1) + JetID::map.at(id2); - } + }*/ // Combination of Jet ID and B-tagging working point (NOTE: NOT USED FOR NOW) - uint16_t JetIDBWP(const JetID::JetID& id, const BWP::BWP& wp){ + /*uint16_t JetIDBWP(const JetID::JetID& id, const BWP::BWP& wp){ return BWP::Count * id + wp; } std::string JetIDBWPStr(const JetID::JetID& id, const BWP::BWP& wp){ return "ID" + JetID::map.at(id) + "_B" + BWP::map.at(wp); - } + }*/ - // Combination of Lepton ID + Lepton Isolation (one lepton) and B-tagging working point for one jet - uint16_t LepIDIsoJetBWP(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp){ + // Combination of Lepton ID + Lepton Isolation (one lepton) and B-tagging working point for one jet (NOTE: NOT USED FOR NOW) + /*uint16_t LepIDIsoJetBWP(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp){ return LepIso::Count*BWP::Count * id + BWP::Count * iso + wp; } std::string LepIDIsoJetBWPStr(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp){ return "ID" + LepID::map.at(id) + "_Iso" + LepIso::map.at(iso) + "_B" + BWP::map.at(wp); - } + }*/ + // Combination of Lepton ID + Lepton Isolation (two leptons) and B-tagging working point for one jet + uint16_t LepLepIDIsoJetBWP(const LepID::LepID& id1, const LepIso::LepIso& iso1, const LepID::LepID& id2, const LepIso::LepIso& iso2, const BWP::BWP& wp){ + return + LepIso::Count*LepID::Count*LepIso::Count*BWP::Count * id1 + + LepID::Count*LepIso::Count*BWP::Count * iso1 + + LepIso::Count*BWP::Count * id2 + + BWP::Count * iso2 + + wp; + } + std::string LepLepIDIsoJetBWPStr(const LepID::LepID& id1, const LepIso::LepIso& iso1, const LepID::LepID& id2, const LepIso::LepIso& iso2, const BWP::BWP& wp){ + return "ID" + LepID::map.at(id1) + LepID::map.at(id2) + "_Iso" + LepIso::map.at(iso1) + LepIso::map.at(iso2) + "_B" + BWP::map.at(wp); + } + // Combination of B-tagging working points for two jets uint16_t JetJetBWP(const BWP::BWP& wp1, const BWP::BWP& wp2){ return BWP::Count * wp1 + wp2; @@ -69,7 +82,7 @@ namespace TTAnalysis { } // Combination of Jet ID and B-tagging working points for two jets (NOTE: NOT USED FOR NOW) - uint16_t JetJetIDBWP(const JetID::JetID& id1, const BWP::BWP& wp1, const JetID::JetID& id2, const BWP::BWP& wp2){ + /*uint16_t JetJetIDBWP(const JetID::JetID& id1, const BWP::BWP& wp1, const JetID::JetID& id2, const BWP::BWP& wp2){ return BWP::Count*JetID::Count*BWP::Count * id1 + JetID::Count*BWP::Count * wp1 + @@ -78,15 +91,15 @@ namespace TTAnalysis { } std::string JetJetIDBWPStr(const JetID::JetID& id1, const BWP::BWP wp1, const JetID::JetID& id2, const BWP::BWP wp2){ return "ID" + JetID::map.at(id1) + JetID::map.at(id2) + "_B" + BWP::map.at(wp1) + BWP::map.at(wp2); - } + }*/ - // Combination of Lepton ID + Lepton Isolation (one lepton) and B-tagging working points for two jets - uint16_t LepIDIsoJetJetBWP(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp1, const BWP::BWP& wp2){ + // Combination of Lepton ID + Lepton Isolation (one lepton) and B-tagging working points for two jets (NOTE: NOT USED FOR NOW) + /*uint16_t LepIDIsoJetJetBWP(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp1, const BWP::BWP& wp2){ return LepIso::Count*BWP::Count*BWP::Count * id + BWP::Count*BWP::Count * iso + BWP::Count * wp1 + wp2; } std::string LepIDIsoJetJetBWPStr(const LepID::LepID& id, const LepIso::LepIso& iso, const BWP::BWP& wp1, const BWP::BWP& wp2){ return "ID" + LepID::map.at(id) + "_Iso" + LepIso::map.at(iso) + "_B" + BWP::map.at(wp1) + BWP::map.at(wp2); - } + }*/ // Combination of Lepton ID, Lepton Isolation, and B-tagging working points for a two-lepton-two-b-jets object uint16_t LepLepIDIsoJetJetBWP(const LepID::LepID& id1, const LepIso::LepIso& iso1, const LepID::LepID& id2, const LepIso::LepIso& iso2, const BWP::BWP& wp1, const BWP::BWP& wp2){ diff --git a/plugins/TTAnalyzer.cc b/plugins/TTAnalyzer.cc index 4757d0e..e66f8e9 100644 --- a/plugins/TTAnalyzer.cc +++ b/plugins/TTAnalyzer.cc @@ -16,6 +16,9 @@ #include #include +#define ECAL_GAP_LOW 1.4442 +#define ECAL_GAP_HIGH 1.566 + // To access VectorUtil::DeltaR() more easily using namespace ROOT::Math; @@ -40,41 +43,36 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, diLeptons_IDIso.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count ); - selJets_selID_DRCut.resize( LepID::Count * LepIso::Count ); - selBJets_DRCut_BWP_PtOrdered.resize( LepID::Count * LepIso::Count * BWP::Count ); - selBJets_DRCut_BWP_CSVv2Ordered.resize( LepID::Count * LepIso::Count * BWP::Count ); + selJets_DRCut_BWP_PtOrdered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count ); + selJets_DRCut_BWP_CSVv2Ordered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count ); - diJets_DRCut.resize( LepID::Count * LepIso::Count ); - diBJets_DRCut_BWP_PtOrdered.resize( LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); - diBJets_DRCut_BWP_CSVv2Ordered.resize( LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); + diJets_DRCut_BWP_PtOrdered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); + diJets_DRCut_BWP_CSVv2Ordered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); - diLepDiJets_DRCut.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count ); - diLepDiBJets_DRCut_BWP_PtOrdered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); - diLepDiBJets_DRCut_BWP_CSVv2Ordered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); + diLepDiJets_DRCut_BWP_PtOrdered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); + diLepDiJets_DRCut_BWP_CSVv2Ordered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); - diLepDiJetsMet_DRCut.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count ); - diLepDiBJetsMet_DRCut_BWP_PtOrdered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); - diLepDiBJetsMet_DRCut_BWP_CSVv2Ordered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); + diLepDiJetsMet_DRCut_BWP_PtOrdered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); + diLepDiJetsMet_DRCut_BWP_CSVv2Ordered.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); ttbar.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count * BWP::Count ); - gen_matched_b.resize( LepID::Count * LepIso::Count , -1); - gen_matched_b_beforeFSR.resize( LepID::Count * LepIso::Count , -1); - gen_matched_bbar.resize( LepID::Count * LepIso::Count , -1); - gen_matched_bbar_beforeFSR.resize( LepID::Count * LepIso::Count , -1); - gen_b_deltaR.resize( LepID::Count * LepIso::Count ); - gen_b_beforeFSR_deltaR.resize( LepID::Count * LepIso::Count ); - gen_bbar_deltaR.resize( LepID::Count * LepIso::Count ); - gen_bbar_beforeFSR_deltaR.resize( LepID::Count * LepIso::Count ); - - if (!m_neutrinos_solver.get()) { - const float topMass = event.isRealData() ? 173.34 : 172.5; - // const float topWidth = event.isRealData() ? 1.41 : 1.50833649; - - const float wMass = event.isRealData() ? 80.385 : 80.419002; + gen_matched_b.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count, -1); + gen_matched_b_beforeFSR.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count, -1); + gen_matched_bbar.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count, -1); + gen_matched_bbar_beforeFSR.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count, -1); + gen_b_deltaR.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count); + gen_b_beforeFSR_deltaR.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count); + gen_bbar_deltaR.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count); + gen_bbar_beforeFSR_deltaR.resize( LepID::Count * LepIso::Count * LepID::Count * LepIso::Count * BWP::Count); + + const float topMass = event.isRealData() ? 173.34 : 172.5; + // const float topWidth = event.isRealData() ? 1.41 : 1.50833649; + const float wMass = event.isRealData() ? 80.385 : 80.419002; // const float wWidth = event.isRealData() ? 2.085 : 2.04759951; - m_neutrinos_solver.reset(new NeutrinosSolver(topMass, wMass)); + if (!m_neutrinos_solver.get()) { + m_neutrinos_solver.reset(new NeutrinosSolver()); } /////////////////////////// @@ -88,6 +86,9 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, const ElectronsProducer& electrons = producers.get(m_electrons_producer); for(uint16_t ielectron = 0; ielectron < electrons.p4.size(); ielectron++){ + if( m_electronRemoveGap && std::abs(electrons.p4[ielectron].Eta()) > ECAL_GAP_LOW && std::abs(electrons.p4[ielectron].Eta()) < ECAL_GAP_HIGH ) + continue; + if( electrons.p4[ielectron].Pt() > m_electronPtCut && std::abs(electrons.p4[ielectron].Eta()) < m_electronEtaCut ){ Lepton m_lepton( @@ -210,7 +211,15 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, m_diLepton.iso[idx] = l1.iso[iso1] && l2.iso[iso2]; } } - + + // Retrieve HLT SF + if(m_diLepton.isElEl && m_hltSF.find("DoubleEG") != m_hltSF.end()) + m_diLepton.hlt_SF = m_hltSF["DoubleEG"].get( { std::abs(l1.p4.Eta()), std::abs(l2.p4.Eta()) } ); + if((m_diLepton.isElMu || m_diLepton.isMuEl) && m_hltSF.find("MuonEG") != m_hltSF.end()) + m_diLepton.hlt_SF = m_hltSF["MuonEG"].get( { std::abs(l1.p4.Eta()), std::abs(l2.p4.Eta()) } ); + if(m_diLepton.isMuMu && m_hltSF.find("DoubleMuon") != m_hltSF.end()) + m_diLepton.hlt_SF = m_hltSF["DoubleMuon"].get( { std::abs(l1.p4.Eta()), std::abs(l2.p4.Eta()) } ); + m_diLepton.DR = VectorUtil::DeltaR(l1.p4, l2.p4); m_diLepton.DEta = TTAnalysis::DeltaEta(l1.p4, l2.p4); m_diLepton.DPhi = VectorUtil::DeltaPhi(l1.p4, l2.p4); @@ -266,54 +275,66 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, m_jet.ID[JetID::T] = jets.passTightID[ijet]; m_jet.ID[JetID::TLV] = jets.passTightLeptonVetoID[ijet]; m_jet.CSVv2 = jets.getBTagDiscriminant(ijet, m_jetCSVv2Name); - m_jet.BWP[BWP::L] = m_jet.CSVv2 > m_jetCSVv2L; - m_jet.BWP[BWP::M] = m_jet.CSVv2 > m_jetCSVv2M; - m_jet.BWP[BWP::T] = m_jet.CSVv2 > m_jetCSVv2T; + m_jet.BWP[BWP::A] = true; + m_jet.BWP[BWP::L] = m_jet.CSVv2 > m_jetCSVv2L && std::abs(jets.p4[ijet].Eta()) < m_bJetEtaCut; + m_jet.BWP[BWP::M] = m_jet.CSVv2 > m_jetCSVv2M && std::abs(jets.p4[ijet].Eta()) < m_bJetEtaCut; + //m_jet.BWP[BWP::T] = m_jet.CSVv2 > m_jetCSVv2T && std::abs(jets.p4[ijet].Eta()) < m_bJetEtaCut; - // Save minimal DR(l,j) using selected leptons, for each Lepton ID/Iso - for(const LepID::LepID& id: LepID::it){ - for(const LepIso::LepIso& iso: LepIso::it){ + // Save minimal DR(l,j) using diLeptons, for each Lepton ID/Iso pair + for(const LepID::LepID& id1: LepID::it){ + for(const LepID::LepID& id2: LepID::it){ + for(const LepIso::LepIso& iso1: LepIso::it){ + for(const LepIso::LepIso& iso2: LepIso::it){ + + uint16_t idx_comb = LepLepIDIso(id1, iso1, id2, iso2); - uint16_t idx_comb = LepIDIso(id, iso); - - for(const uint16_t& lepIdx: leptons_IDIso[idx_comb]){ - const Lepton& m_lepton = leptons[lepIdx]; - float DR = (float) VectorUtil::DeltaR(jets.p4[ijet], m_lepton.p4); - if( DR < m_jet.minDRjl_lepIDIso[idx_comb] ) - m_jet.minDRjl_lepIDIso[idx_comb] = DR; - } - - // Save the indices to Jets passing the selected jetID and minDRjl > cut for this lepton ID/Iso - if( m_jet.minDRjl_lepIDIso[idx_comb] > m_jetDRleptonCut && jetIDAccessor(jets, ijet, m_jetID) ){ - selJets_selID_DRCut[idx_comb].push_back(jetCounter); - - // Out of these, save the indices for different b-tagging working points - for(const BWP::BWP& wp: BWP::it){ - uint16_t idx_comb_b = LepIDIsoJetBWP(id, iso, wp); - if ((m_jet.BWP[wp]) && (std::abs(m_jet.p4.Eta()) < m_bJetEtaCut)) - selBJets_DRCut_BWP_PtOrdered[idx_comb_b].push_back(jetCounter); + for(const uint16_t& diLepIdx: diLeptons_IDIso[idx_comb]){ + const DiLepton& m_diLepton = diLeptons[diLepIdx]; + float DR1 = (float) VectorUtil::DeltaR(jets.p4[ijet], leptons[m_diLepton.lidxs.first].p4); + float DR2 = (float) VectorUtil::DeltaR(jets.p4[ijet], leptons[m_diLepton.lidxs.second].p4); + float DR = std::min(DR1, DR2); + if( DR < m_jet.minDRjl_lepLepIDIso[idx_comb] ) + m_jet.minDRjl_lepLepIDIso[idx_comb] = DR; + } + + // Save the indices to Jets passing the selected jetID and minDRjl > cut for this lepton ID/Iso, including the indices for different b-tagging working points + if( m_jet.minDRjl_lepLepIDIso[idx_comb] > m_jetDRleptonCut && jetIDAccessor(jets, ijet, m_jetID) ){ + for(const BWP::BWP& wp: BWP::it){ + uint16_t idx_comb_b = LepLepIDIsoJetBWP(id1, iso1, id2, iso2, wp); + if (m_jet.BWP[wp]) + selJets_DRCut_BWP_PtOrdered[idx_comb_b].push_back(jetCounter); + } + } } } } } - if(jetIDAccessor(jets, ijet, m_jetID)) // Save the indices to Jets passing the selected jet ID - selJets_selID.push_back(jetCounter); - selJets.push_back(m_jet); jetCounter++; } } - // Sort the b-jets according to decreasing CSVv2 value - selBJets_DRCut_BWP_CSVv2Ordered = selBJets_DRCut_BWP_PtOrdered; - for(const LepID::LepID& id: LepID::it){ - for(const LepIso::LepIso& iso: LepIso::it){ - for(const BWP::BWP& wp: BWP::it){ - uint16_t idx_comb_b = LepIDIsoJetBWP(id, iso, wp); - std::sort(selBJets_DRCut_BWP_CSVv2Ordered[idx_comb_b].begin(), selBJets_DRCut_BWP_CSVv2Ordered[idx_comb_b].end(), jetBTagDiscriminantSorter(jets, m_jetCSVv2Name, selJets)); + // Sort the jets according to decreasing CSVv2 value + selJets_DRCut_BWP_CSVv2Ordered = selJets_DRCut_BWP_PtOrdered; + + for(const LepID::LepID& id1: LepID::it){ + for(const LepID::LepID& id2: LepID::it){ + + for(const LepIso::LepIso& iso1: LepIso::it){ + for(const LepIso::LepIso& iso2: LepIso::it){ + + for(const BWP::BWP& wp: BWP::it){ + + uint16_t idx_comb_b = LepLepIDIsoJetBWP(id1, iso1, id2, iso2, wp); + std::sort(selJets_DRCut_BWP_CSVv2Ordered[idx_comb_b].begin(), selJets_DRCut_BWP_CSVv2Ordered[idx_comb_b].end(), jetBTagDiscriminantSorter(jets, m_jetCSVv2Name, selJets)); + + } + + } } + } } @@ -325,21 +346,19 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, std::cout << "Dijets" << std::endl; #endif - // Next, construct DiJets out of selected jets with selected ID (not accounting for minDRjl here) + // Next, construct DiJets out of selected jets uint16_t diJetCounter(0); - for(uint16_t j1 = 0; j1 < selJets_selID.size(); j1++){ - for(uint16_t j2 = j1 + 1; j2 < selJets_selID.size(); j2++){ - const uint16_t jidx1 = selJets_selID[j1]; - const Jet& jet1 = selJets[jidx1]; - const uint16_t jidx2 = selJets_selID[j2]; - const Jet& jet2 = selJets[jidx2]; + for(uint16_t j1 = 0; j1 < selJets.size(); j1++){ + for(uint16_t j2 = j1 + 1; j2 < selJets.size(); j2++){ + const Jet& jet1 = selJets[j1]; + const Jet& jet2 = selJets[j2]; DiJet m_diJet; m_diJet.p4 = jet1.p4 + jet2.p4; m_diJet.idxs = std::make_pair(jet1.idx, jet2.idx); - m_diJet.jidxs = std::make_pair(jidx1, jidx2); + m_diJet.jidxs = std::make_pair(j1, j2); m_diJet.DR = VectorUtil::DeltaR(jet1.p4, jet2.p4); m_diJet.DEta = DeltaEta(jet1.p4, jet2.p4); @@ -352,30 +371,30 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, } } - for(const LepID::LepID& id: LepID::it){ - for(const LepIso::LepIso& iso: LepIso::it){ - uint16_t combIDIso = LepIDIso(id, iso); + for(const LepID::LepID& id1: LepID::it){ + for(const LepID::LepID& id2: LepID::it){ + for(const LepIso::LepIso& iso1: LepIso::it){ + for(const LepIso::LepIso& iso2: LepIso::it){ + uint16_t combIDIso = LepLepIDIso(id1, iso1, id2, iso2); - m_diJet.minDRjl_lepIDIso[combIDIso] = std::min(jet1.minDRjl_lepIDIso[combIDIso], jet2.minDRjl_lepIDIso[combIDIso]); - - // Save the DiJets which have minDRjl>cut, for each leptonIDIso - if(m_diJet.minDRjl_lepIDIso[combIDIso] > m_jetDRleptonCut){ - diJets_DRCut[combIDIso].push_back(diJetCounter); - - // Out of these, save di-b-jets for each combination of b-tagging working points - for(const BWP::BWP& wp1: BWP::it){ - for(const BWP::BWP& wp2: BWP::it){ - uint16_t combB = JetJetBWP(wp1, wp2); - uint16_t combAll = LepIDIsoJetJetBWP(id, iso, wp1, wp2); - if ((m_diJet.BWP[combB]) - && (std::abs(jet1.p4.Eta()) < m_bJetEtaCut) - && (std::abs(jet2.p4.Eta()) < m_bJetEtaCut)) - diBJets_DRCut_BWP_PtOrdered[combAll].push_back(diJetCounter); + m_diJet.minDRjl_lepLepIDIso[combIDIso] = std::min(jet1.minDRjl_lepLepIDIso[combIDIso], jet2.minDRjl_lepLepIDIso[combIDIso]); + + // Save the DiJets which have minDRjl>cut, for each leptonIDIso and combination of b-tagging working points + if(m_diJet.minDRjl_lepLepIDIso[combIDIso] > m_jetDRleptonCut){ + + for(const BWP::BWP& wp1: BWP::it){ + for(const BWP::BWP& wp2: BWP::it){ + uint16_t combB = JetJetBWP(wp1, wp2); + uint16_t combAll = LepLepIDIsoJetJetBWP(id1, iso1, id2, iso2, wp1, wp2); + if (m_diJet.BWP[combB]) + diJets_DRCut_BWP_PtOrdered[combAll].push_back(diJetCounter); + } + } + } + } - } - } } @@ -384,16 +403,27 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, } } - // Order selected di-b-jets according to decreasing CSVv2 discriminant - diBJets_DRCut_BWP_CSVv2Ordered = diBJets_DRCut_BWP_PtOrdered; - for(const LepID::LepID& id: LepID::it){ - for(const LepIso::LepIso& iso: LepIso::it){ - for(const BWP::BWP& wp1: BWP::it){ - for(const BWP::BWP& wp2: BWP::it){ - uint16_t idx_comb_b = LepIDIsoJetJetBWP(id, iso, wp1, wp2); - std::sort(diBJets_DRCut_BWP_CSVv2Ordered[idx_comb_b].begin(), diBJets_DRCut_BWP_CSVv2Ordered[idx_comb_b].end(), diJetBTagDiscriminantSorter(jets, m_jetCSVv2Name, diJets)); + // Order selected di-jets according to decreasing CSVv2 discriminant + diJets_DRCut_BWP_CSVv2Ordered = diJets_DRCut_BWP_PtOrdered; + + for(const LepID::LepID& id1: LepID::it){ + for(const LepID::LepID& id2: LepID::it){ + + for(const LepIso::LepIso& iso1: LepIso::it){ + for(const LepIso::LepIso& iso2: LepIso::it){ + + for(const BWP::BWP& wp1: BWP::it){ + for(const BWP::BWP& wp2: BWP::it){ + + uint16_t idx_comb_all = LepLepIDIsoJetJetBWP(id1, iso1, id2, iso2, wp1, wp2); + std::sort(diJets_DRCut_BWP_CSVv2Ordered[idx_comb_all].begin(), diJets_DRCut_BWP_CSVv2Ordered[idx_comb_all].end(), diJetBTagDiscriminantSorter(jets, m_jetCSVv2Name, diJets)); + + } + } + } } + } } @@ -414,7 +444,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, for(uint16_t dijet = 0; dijet < diJets.size(); dijet++){ const DiJet& m_diJet = diJets[dijet]; - + DiLepDiJet m_diLepDiJet(m_diLepton, dilep, m_diJet, dijet); m_diLepDiJet.minDRjl = std::min( { @@ -462,27 +492,19 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, for(const LepIso::LepIso& iso2: LepIso::it){ uint16_t combID = LepLepID(id1, id2); - LepID::LepID minID = std::min(id1, id2); - uint16_t combIso = LepLepIso(iso1, iso2); - LepIso::LepIso minIso = std::min(iso1, iso2); - - uint16_t minCombIDIso = LepIDIso(minID, minIso); - uint16_t diLepCombIDIso = LepLepIDIso(id1, iso1, id2, iso2); + uint16_t diLepIDIso = LepLepIDIso(id1, iso1, id2, iso2); // Store objects for each combined lepton ID/Iso, with jets having minDRjl>cut for leptons corresponding to the loosest combination of the aforementioned ID/Iso - if(m_diLepton.ID[combID] && m_diLepton.iso[combIso] && m_diJet.minDRjl_lepIDIso[minCombIDIso] > m_jetDRleptonCut){ - diLepDiJets_DRCut[diLepCombIDIso].push_back(diLepDiJetCounter); + if(m_diLepton.ID[combID] && m_diLepton.iso[combIso] && m_diJet.minDRjl_lepLepIDIso[diLepIDIso] > m_jetDRleptonCut){ // Out of these, store combinations of b-tagging working points for(const BWP::BWP& wp1: BWP::it){ for(const BWP::BWP& wp2: BWP::it){ uint16_t combB = JetJetBWP(wp1, wp2); uint16_t combAll = LepLepIDIsoJetJetBWP(id1, iso1, id2, iso2, wp1, wp2); - if ((m_diJet.BWP[combB]) - && (std::abs(jets.p4[m_diJet.idxs.first].Eta()) < m_bJetEtaCut) - && (std::abs(jets.p4[m_diJet.idxs.second].Eta()) < m_bJetEtaCut)) - diLepDiBJets_DRCut_BWP_PtOrdered[combAll].push_back(diLepDiJetCounter); + if (m_diJet.BWP[combB]) + diLepDiJets_DRCut_BWP_PtOrdered[combAll].push_back(diLepDiJetCounter); } } // end b-jet loops @@ -497,8 +519,8 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, } // end dijet loop } // end dilepton loop - // Order selected di-lepton-di-b-jets according to decreasing CSVv2 discriminant - diLepDiBJets_DRCut_BWP_CSVv2Ordered = diLepDiBJets_DRCut_BWP_PtOrdered; + // Order selected di-lepton-di-jets according to decreasing CSVv2 discriminant + diLepDiJets_DRCut_BWP_CSVv2Ordered = diLepDiJets_DRCut_BWP_PtOrdered; for(const LepID::LepID& id1: LepID::it){ for(const LepID::LepID& id2: LepID::it){ @@ -510,7 +532,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, for(const BWP::BWP& wp2: BWP::it){ uint16_t idx_comb_all = LepLepIDIsoJetJetBWP(id1, iso1, id2, iso2, wp1, wp2); - std::sort(diLepDiBJets_DRCut_BWP_CSVv2Ordered[idx_comb_all].begin(), diLepDiBJets_DRCut_BWP_CSVv2Ordered[idx_comb_all].end(), diJetBTagDiscriminantSorter(jets, m_jetCSVv2Name, diLepDiJets)); + std::sort(diLepDiJets_DRCut_BWP_CSVv2Ordered[idx_comb_all].begin(), diLepDiJets_DRCut_BWP_CSVv2Ordered[idx_comb_all].end(), diJetBTagDiscriminantSorter(jets, m_jetCSVv2Name, diLepDiJets)); } } @@ -591,29 +613,19 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, for(const LepIso::LepIso& iso2: LepIso::it){ uint16_t combID = LepLepID(id1, id2); - LepID::LepID minID = std::min(id1, id2); - uint16_t combIso = LepLepIso(iso1, iso2); - LepIso::LepIso minIso = std::min(iso1, iso2); - - uint16_t minCombIDIso = LepIDIso(minID, minIso); - uint16_t diLepCombIDIso = LepLepIDIso(id1, iso1, id2, iso2); + uint16_t diLepIDIso = LepLepIDIso(id1, iso1, id2, iso2); // Store objects for each combined lepton ID/Iso, with jets having minDRjl>cut for leptons corresponding to the loosest combination of the aforementioned ID/Iso - // First regular MET - if(m_diLepDiJetMet.diLepton->ID[combID] && m_diLepDiJetMet.diLepton->iso[combIso] && m_diLepDiJetMet.diJet->minDRjl_lepIDIso[minCombIDIso] > m_jetDRleptonCut){ - diLepDiJetsMet_DRCut[diLepCombIDIso].push_back(i); - + if(m_diLepDiJetMet.diLepton->ID[combID] && m_diLepDiJetMet.diLepton->iso[combIso] && m_diLepDiJetMet.diJet->minDRjl_lepLepIDIso[diLepIDIso] > m_jetDRleptonCut){ // Out of these, store combinations of b-tagging working points for(const BWP::BWP& wp1: BWP::it){ for(const BWP::BWP& wp2: BWP::it){ uint16_t combB = JetJetBWP(wp1, wp2); uint16_t combAll = LepLepIDIsoJetJetBWP(id1, iso1, id2, iso2, wp1, wp2); - if ((m_diLepDiJetMet.diJet->BWP[combB]) - && (std::abs(jets.p4[m_diLepDiJetMet.diJet->idxs.first].Eta()) < m_bJetEtaCut) - && (std::abs(jets.p4[m_diLepDiJetMet.diJet->idxs.second].Eta()) < m_bJetEtaCut)) - diLepDiBJetsMet_DRCut_BWP_PtOrdered[combAll].push_back(i); + if (m_diLepDiJetMet.diJet->BWP[combB]) + diLepDiJetsMet_DRCut_BWP_PtOrdered[combAll].push_back(i); } } // end b-jet loops @@ -627,8 +639,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, } // end diLepDiJet loop // Store objects according to CSVv2 - // First regular MET - diLepDiBJetsMet_DRCut_BWP_CSVv2Ordered = diLepDiBJetsMet_DRCut_BWP_PtOrdered; + diLepDiJetsMet_DRCut_BWP_CSVv2Ordered = diLepDiJetsMet_DRCut_BWP_PtOrdered; for(const LepID::LepID& id1: LepID::it){ for(const LepID::LepID& id2: LepID::it){ @@ -639,7 +650,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, for(const BWP::BWP& wp2: BWP::it){ uint16_t idx_comb_all = LepLepIDIsoJetJetBWP(id1, iso1, id2, iso2, wp1, wp2); - std::sort(diLepDiBJetsMet_DRCut_BWP_CSVv2Ordered[idx_comb_all].begin(), diLepDiBJetsMet_DRCut_BWP_CSVv2Ordered[idx_comb_all].end(), diJetBTagDiscriminantSorter(jets, m_jetCSVv2Name, diLepDiJetsMet)); + std::sort(diLepDiJetsMet_DRCut_BWP_CSVv2Ordered[idx_comb_all].begin(), diLepDiJetsMet_DRCut_BWP_CSVv2Ordered[idx_comb_all].end(), diJetBTagDiscriminantSorter(jets, m_jetCSVv2Name, diLepDiJetsMet)); } } @@ -654,11 +665,11 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, // MTT // /////////////////////////// - #ifdef _TT_DEBUG_ +#ifdef _TT_DEBUG_ std::cout << "Reconstructing mtt" << std::endl; - #endif +#endif -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "Reconstructing ttbar system" << std::endl; #endif @@ -675,7 +686,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, std::vector> ttbar_event_sols; - for (const auto& idx: diLepDiBJetsMet_DRCut_BWP_CSVv2Ordered[idx_comb_all]) { + for (const auto& idx: diLepDiJetsMet_DRCut_BWP_CSVv2Ordered[idx_comb_all]) { using namespace TTAnalysis; @@ -686,7 +697,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, NeutrinosSolver::LorentzVector met_p4(met.p4); -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "Objects:" << std::endl; std::cout << "\t Lepton 1: " << lepton1_p4 << std::endl; std::cout << "\t b-jet 1: " << bjet1_p4 << std::endl; @@ -694,43 +705,43 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, std::cout << "\t b-jet 2: " << bjet2_p4 << std::endl; #endif - auto sols = m_neutrinos_solver->getNeutrinos(lepton1_p4, lepton2_p4, bjet1_p4, bjet2_p4, met_p4); + auto sols = m_neutrinos_solver->getNeutrinos(lepton1_p4, lepton2_p4, bjet1_p4, bjet2_p4, met_p4, wMass, wMass, topMass, topMass); -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "Got " << sols.size() << " solutions for neutrinos" << std::endl; #endif std::vector ttbar_sols; for (auto& sol: sols) { -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "\t Neutrino 1: " << sol.first << std::endl; std::cout << "\t Neutrino 2: " << sol.second << std::endl; #endif ttbar_sols.push_back(TTBar(idx, myLorentzVector(lepton1_p4 + bjet1_p4 + sol.first), myLorentzVector(lepton2_p4 + bjet2_p4 + sol.second))); -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "mtt: " << ttbar_sols.back().p4.M() << std::endl; #endif } -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "Swapping b-jets and recomputing solutions" << std::endl; #endif // Swap b-jets std::swap(bjet1_p4, bjet2_p4); - sols = m_neutrinos_solver->getNeutrinos(lepton1_p4, lepton2_p4, bjet1_p4, bjet2_p4, met_p4); + sols = m_neutrinos_solver->getNeutrinos(lepton1_p4, lepton2_p4, bjet1_p4, bjet2_p4, met_p4, wMass, wMass, topMass, topMass); -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "Got " << sols.size() << " solutions for neutrinos" << std::endl; #endif for (auto& sol: sols) { -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "\t Neutrino 1: " << sol.first << std::endl; std::cout << "\t Neutrino 2: " << sol.second << std::endl; #endif ttbar_sols.push_back(TTBar(idx, myLorentzVector(lepton1_p4 + bjet1_p4 + sol.first), myLorentzVector(lepton2_p4 + bjet2_p4 + sol.second))); -#if TT_MTT_DEBUG +#ifdef TT_MTT_DEBUG std::cout << "mtt: " << ttbar_sols.back().p4.M() << std::endl; #endif } @@ -756,22 +767,22 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, // TRIGGER // /////////////////////////// - #ifdef _TT_DEBUG_ +#ifdef _TT_DEBUG_ std::cout << "Trigger" << std::endl; - #endif +#endif if (producers.exists("hlt")) { const HLTProducer& hlt = producers.get("hlt"); if (hlt.paths.empty()) { -#if TT_HLT_DEBUG +#ifdef TT_HLT_DEBUG std::cout << "No HLT path triggered for this event. Skipping HLT matching." << std::endl; #endif goto after_hlt_matching; } -#if TT_HLT_DEBUG +#ifdef TT_HLT_DEBUG std::cout << "HLT path triggered for this event:" << std::endl; for (const std::string& path: hlt.paths) { std::cout << "\t" << path << std::endl; @@ -787,7 +798,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, if (lepton.hlt_already_tried_matching) return lepton.hlt_idx; -#if TT_HLT_DEBUG +#ifdef TT_HLT_DEBUG std::cout << "Trying to match offline lepton: " << std::endl; std::cout << "\tMuon? " << lepton.isMu << " ; Pt: " << lepton.p4.Pt() << " ; Eta: " << lepton.p4.Eta() << " ; Phi: " << lepton.p4.Phi() << " ; E: " << lepton.p4.E() << std::endl; #endif @@ -812,7 +823,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, } } -#if TT_HLT_DEBUG +#ifdef TT_HLT_DEBUG if (index != -1) { std::cout << "\033[32mMatched with online object:\033[00m" << std::endl; std::cout << "\tPDG Id: " << hlt.object_pdg_id[index] << " ; Pt: " << hlt.object_p4[index].Pt() << " ; Eta: " << hlt.object_p4[index].Eta() << " ; Phi: " << hlt.object_p4[index].Phi() << " ; E: " << hlt.object_p4[index].E() << std::endl; @@ -847,9 +858,9 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, // GEN INFO // /////////////////////////// - #ifdef _TT_DEBUG_ - std::cout << "Generator" << std::endl; - #endif +#ifdef _TT_DEBUG_ + std::cout << "Generator" << std::endl; +#endif if (event.isRealData()) return; @@ -876,7 +887,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, return false; }; -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::function print_mother_chain = [&gen_particles, &print_mother_chain](size_t p) { if (gen_particles.pruned_mothers_index[p].empty()) { @@ -920,7 +931,9 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, gen_neutrino_tbar_beforeFSR = -1; gen_matched_lepton_t = -1; + gen_matched_lepton_beforeFSR_t = -1; gen_matched_lepton_tbar = -1; + gen_matched_lepton_beforeFSR_tbar = -1; size_t gen_index(0); for (size_t i = 0; i < gen_particles.pruned_pdg_id.size(); i++) { @@ -940,7 +953,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, if (! flags.fromHardProcess()) continue; -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "---" << std::endl; std::cout << "Gen particle #" << i << ": PDG id: " << gen_particles.pruned_pdg_id[i]; print_mother_chain(i); @@ -980,7 +993,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, (gen_jet1_tbar_beforeFSR != -1 && std::abs(genParticles[gen_jet1_tbar_beforeFSR].pdg_id) == 5) || (gen_jet2_tbar_beforeFSR != -1 && std::abs(genParticles[gen_jet2_tbar_beforeFSR].pdg_id) == 5)) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "A quark coming from W decay is a b" << std::endl; #endif @@ -988,19 +1001,19 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, ! (gen_jet2_tbar_beforeFSR != -1 && pruned_decays_from(i, genParticles[gen_jet2_tbar_beforeFSR].pruned_idx)) && ! (gen_jet1_t_beforeFSR != -1 && pruned_decays_from(i, genParticles[gen_jet1_t_beforeFSR].pruned_idx)) && ! (gen_jet2_t_beforeFSR != -1 && pruned_decays_from(i, genParticles[gen_jet2_t_beforeFSR].pruned_idx))) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "This after-FSR b quark is not coming from a W decay" << std::endl; #endif FILL_GEN_COLL(b); continue; } -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG else { std::cout << "This after-FSR b quark comes from a W decay" << std::endl; } #endif } else { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "Assigning gen_b" << std::endl; #endif FILL_GEN_COLL(b); @@ -1010,7 +1023,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, FILL_GEN_COLL(b); continue; } else { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "This should not happen!" << std::endl; #endif } @@ -1026,7 +1039,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, (gen_jet1_tbar_beforeFSR != -1 && std::abs(genParticles[gen_jet1_tbar_beforeFSR].pdg_id) == 5) || (gen_jet2_tbar_beforeFSR != -1 && std::abs(genParticles[gen_jet2_tbar_beforeFSR].pdg_id) == 5)) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "A quark coming from W decay is a bbar" << std::endl; #endif @@ -1034,19 +1047,19 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, ! (gen_jet2_tbar_beforeFSR != -1 && pruned_decays_from(i, genParticles[gen_jet2_tbar_beforeFSR].pruned_idx)) && ! (gen_jet1_t_beforeFSR != -1 && pruned_decays_from(i, genParticles[gen_jet1_t_beforeFSR].pruned_idx)) && ! (gen_jet2_t_beforeFSR != -1 && pruned_decays_from(i, genParticles[gen_jet2_t_beforeFSR].pruned_idx))) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "This after-fsr b anti-quark is not coming from a W decay" << std::endl; #endif FILL_GEN_COLL(bbar); continue; } -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG else { std::cout << "This after-fsr b anti-quark comes from a W decay" << std::endl; } #endif } else { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "Assigning gen_bbar" << std::endl; #endif FILL_GEN_COLL(bbar); @@ -1062,7 +1075,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, continue; if (gen_t != -1 && from_t_decay) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "Coming from the top chain decay" << std::endl; #endif if (a_pdg_id >= 1 && a_pdg_id <= 5) { @@ -1075,7 +1088,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, std::cout << "Error: unknown particle coming from top decay - #" << i << " ; PDG Id: " << pdg_id << std::endl; } } else if (gen_tbar != -1 && from_tbar_decay) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "Coming from the anti-top chain decay" << std::endl; #endif if (a_pdg_id >= 1 && a_pdg_id <= 5) { @@ -1091,7 +1104,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, } if (!gen_t || !gen_tbar) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "This is not a ttbar event" << std::endl; #endif gen_ttbar_decay_type = NotTT; @@ -1099,7 +1112,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, } if ((gen_jet1_t != -1) && (gen_jet2_t != -1) && (gen_jet1_tbar != -1) && (gen_jet2_tbar != -1)) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "Hadronic ttbar decay" << std::endl; #endif gen_ttbar_decay_type = Hadronic; @@ -1108,7 +1121,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, ((gen_lepton_t == -1) && (gen_lepton_tbar != -1)) ) { -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "Semileptonic ttbar decay" << std::endl; #endif @@ -1128,7 +1141,7 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, uint16_t lepton_t_pdg_id = std::abs(genParticles[gen_lepton_t].pdg_id); uint16_t lepton_tbar_pdg_id = std::abs(genParticles[gen_lepton_tbar].pdg_id); -#if TT_GEN_DEBUG +#ifdef TT_GEN_DEBUG std::cout << "Dileptonic ttbar decay" << std::endl; #endif @@ -1178,7 +1191,9 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, if (gen_ttbar_decay_type > Hadronic) { float min_dr_lepton_t = std::numeric_limits::max(); + float min_dr_lepton_beforeFSR_t = std::numeric_limits::max(); float min_dr_lepton_tbar = std::numeric_limits::max(); + float min_dr_lepton_beforeFSR_tbar = std::numeric_limits::max(); size_t lepton_index = 0; for (auto& lepton: leptons) { @@ -1203,63 +1218,86 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, gen_lepton_tbar_deltaR.push_back(dr); } + if (gen_lepton_t_beforeFSR != -1) { + float dr = VectorUtil::DeltaR(genParticles[gen_lepton_t_beforeFSR].p4, lepton.p4); + if (dr < min_dr_lepton_beforeFSR_t && + (std::abs(lepton.pdg_id()) == std::abs(genParticles[gen_lepton_t_beforeFSR].pdg_id))) { + min_dr_lepton_beforeFSR_t = dr; + gen_matched_lepton_beforeFSR_t = lepton_index; + } + gen_lepton_beforeFSR_t_deltaR.push_back(dr); + } + + if (gen_lepton_tbar_beforeFSR != -1) { + float dr = VectorUtil::DeltaR(genParticles[gen_lepton_tbar_beforeFSR].p4, lepton.p4); + if (dr < min_dr_lepton_beforeFSR_tbar && + (std::abs(lepton.pdg_id()) == std::abs(genParticles[gen_lepton_tbar_beforeFSR].pdg_id))) { + min_dr_lepton_beforeFSR_tbar = dr; + gen_matched_lepton_beforeFSR_tbar = lepton_index; + } + gen_lepton_beforeFSR_tbar_deltaR.push_back(dr); + } + lepton_index++; } } // Match b quarks to jets - const float MIN_DR_JETS = 0.8; - for (const auto& id: LepID::it) { - for (const auto& iso: LepIso::it) { - uint16_t IdWP = LepIDIso(id, iso); - - float min_dr_b = MIN_DR_JETS; - float min_dr_bbar = MIN_DR_JETS; - float min_dr_b_beforeFSR = MIN_DR_JETS; - float min_dr_bbar_beforeFSR = MIN_DR_JETS; - size_t jet_index = 0; - - int16_t local_gen_matched_b = -1; - int16_t local_gen_matched_bbar = -1; - int16_t local_gen_matched_b_beforeFSR = -1; - int16_t local_gen_matched_bbar_beforeFSR = -1; - for (auto& jet: selJets_selID_DRCut[IdWP]) { - float dr = VectorUtil::DeltaR(genParticles[gen_b].p4, jets.p4[jet]); - if (dr < min_dr_b) { - min_dr_b = dr; - local_gen_matched_b = jet_index; - } - gen_b_deltaR[IdWP].push_back(dr); + for (const auto& id1: LepID::it) { + for (const auto& id2: LepID::it) { + for (const auto& iso1: LepIso::it) { + for (const auto& iso2: LepIso::it) { + for (const auto& bwp: BWP::it) { + uint16_t combIdx = LepLepIDIsoJetBWP(id1, iso1, id2, iso2, bwp); + + float min_dr_b = MIN_DR_JETS; + float min_dr_bbar = MIN_DR_JETS; + float min_dr_b_beforeFSR = MIN_DR_JETS; + float min_dr_bbar_beforeFSR = MIN_DR_JETS; + + int16_t local_gen_matched_b = -1; + int16_t local_gen_matched_bbar = -1; + int16_t local_gen_matched_b_beforeFSR = -1; + int16_t local_gen_matched_bbar_beforeFSR = -1; + + for (auto& jet_idx: selJets_DRCut_BWP_PtOrdered[combIdx]) { + float dr = VectorUtil::DeltaR(genParticles[gen_b].p4, selJets[jet_idx].p4); + if (dr < min_dr_b) { + min_dr_b = dr; + local_gen_matched_b = jet_idx; + } + gen_b_deltaR[combIdx].push_back(dr); - dr = VectorUtil::DeltaR(genParticles[gen_b_beforeFSR].p4, jets.p4[jet]); - if (dr < min_dr_b_beforeFSR) { - min_dr_b_beforeFSR = dr; - local_gen_matched_b_beforeFSR = jet_index; - } - gen_b_beforeFSR_deltaR[IdWP].push_back(dr); + dr = VectorUtil::DeltaR(genParticles[gen_b_beforeFSR].p4, selJets[jet_idx].p4); + if (dr < min_dr_b_beforeFSR) { + min_dr_b_beforeFSR = dr; + local_gen_matched_b_beforeFSR = jet_idx; + } + gen_b_beforeFSR_deltaR[combIdx].push_back(dr); - dr = VectorUtil::DeltaR(genParticles[gen_bbar].p4, jets.p4[jet]); - if (dr < min_dr_bbar) { - min_dr_bbar = dr; - local_gen_matched_bbar = jet_index; - } - gen_bbar_deltaR[IdWP].push_back(dr); + dr = VectorUtil::DeltaR(genParticles[gen_bbar].p4, selJets[jet_idx].p4); + if (dr < min_dr_bbar) { + min_dr_bbar = dr; + local_gen_matched_bbar = jet_idx; + } + gen_bbar_deltaR[combIdx].push_back(dr); - dr = VectorUtil::DeltaR(genParticles[gen_bbar_beforeFSR].p4, jets.p4[jet]); - if (dr < min_dr_bbar_beforeFSR) { - min_dr_bbar_beforeFSR = dr; - local_gen_matched_bbar_beforeFSR = jet_index; + dr = VectorUtil::DeltaR(genParticles[gen_bbar_beforeFSR].p4, selJets[jet_idx].p4); + if (dr < min_dr_bbar_beforeFSR) { + min_dr_bbar_beforeFSR = dr; + local_gen_matched_bbar_beforeFSR = jet_idx; + } + gen_bbar_beforeFSR_deltaR[combIdx].push_back(dr); } - gen_bbar_beforeFSR_deltaR[IdWP].push_back(dr); - jet_index++; + gen_matched_b[combIdx] = local_gen_matched_b; + gen_matched_bbar[combIdx] = local_gen_matched_bbar; + gen_matched_b_beforeFSR[combIdx] = local_gen_matched_b_beforeFSR; + gen_matched_bbar_beforeFSR[combIdx] = local_gen_matched_bbar_beforeFSR; + } } - - gen_matched_b[IdWP] = local_gen_matched_b; - gen_matched_bbar[IdWP] = local_gen_matched_bbar; - gen_matched_b_beforeFSR[IdWP] = local_gen_matched_b_beforeFSR; - gen_matched_bbar_beforeFSR[IdWP] = local_gen_matched_bbar_beforeFSR; + } } } @@ -1271,9 +1309,9 @@ void TTAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup, gen_bbar_lepton_tbar_deltaR = VectorUtil::DeltaR(genParticles[gen_bbar].p4, genParticles[gen_lepton_tbar].p4); } - #ifdef _TT_DEBUG_ +#ifdef _TT_DEBUG_ std::cout << "End event." << std::endl; - #endif +#endif } diff --git a/plugins/TTDileptonCategories.cc b/plugins/TTDileptonCategories.cc index 2438054..26c2a20 100644 --- a/plugins/TTDileptonCategories.cc +++ b/plugins/TTDileptonCategories.cc @@ -55,6 +55,7 @@ void ElElCategory::register_cuts(CutManager& manager) { manager.new_cut(baseStrCategory + postFix, baseStrCategory + postFix); manager.new_cut(baseStrExtraDiLeptonVeto + postFix, baseStrExtraDiLeptonVeto + postFix); + manager.new_cut(baseStrDiLeptonTriggerBit + postFix, baseStrDiLeptonTriggerBit + postFix); manager.new_cut(baseStrDiLeptonTriggerMatch + postFix, baseStrDiLeptonTriggerMatch + postFix); manager.new_cut(baseStrMllCut + postFix, baseStrMllCut + postFix); manager.new_cut(baseStrMllZVetoCut + postFix, baseStrMllZVetoCut + postFix); @@ -87,8 +88,12 @@ void ElElCategory::evaluate_cuts_post_analyzers(CutManager& manager, const Produ if(m_diLepton.isElEl) { manager.pass_cut(baseStrCategory + postFix); + // Check that one of the corresponding triggers has fired + if( checkHLT(hlt, HLT::DoubleEG) ) + manager.pass_cut(baseStrDiLeptonTriggerBit + postFix); + if(m_diLepton.hlt_idxs.first >= 0 && m_diLepton.hlt_idxs.second >= 0){ - // We have fired a trigger. Now, check that it is actually a DoubleEG trigger + // We have fired a trigger. Now, check that the objects are correctly matched to the same DoubleEG trigger. if( checkHLT(hlt, m_diLepton.hlt_idxs.first, m_diLepton.hlt_idxs.second, HLT::DoubleEG) ) manager.pass_cut(baseStrDiLeptonTriggerMatch + postFix); } @@ -163,6 +168,7 @@ void ElMuCategory::register_cuts(CutManager& manager) { manager.new_cut(baseStrCategory + postFix, baseStrCategory + postFix); manager.new_cut(baseStrExtraDiLeptonVeto + postFix, baseStrExtraDiLeptonVeto + postFix); + manager.new_cut(baseStrDiLeptonTriggerBit + postFix, baseStrDiLeptonTriggerBit + postFix); manager.new_cut(baseStrDiLeptonTriggerMatch + postFix, baseStrDiLeptonTriggerMatch + postFix); manager.new_cut(baseStrMllCut + postFix, baseStrMllCut + postFix); manager.new_cut(baseStrMllZVetoCut + postFix, baseStrMllZVetoCut + postFix); @@ -195,6 +201,10 @@ void ElMuCategory::evaluate_cuts_post_analyzers(CutManager& manager, const Produ if(m_diLepton.isElMu) { manager.pass_cut(baseStrCategory + postFix); + // Check that one of the corresponding triggers has fired + if( checkHLT(hlt, HLT::MuonEG) ) + manager.pass_cut(baseStrDiLeptonTriggerBit + postFix); + if(m_diLepton.hlt_idxs.first >= 0 && m_diLepton.hlt_idxs.second >= 0){ // We have fired a trigger. Now, check that it is actually a MuonEG trigger if( checkHLT(hlt, m_diLepton.hlt_idxs.first, m_diLepton.hlt_idxs.second, HLT::MuonEG) ) @@ -271,6 +281,7 @@ void MuElCategory::register_cuts(CutManager& manager) { manager.new_cut(baseStrCategory + postFix, baseStrCategory + postFix); manager.new_cut(baseStrExtraDiLeptonVeto + postFix, baseStrExtraDiLeptonVeto + postFix); + manager.new_cut(baseStrDiLeptonTriggerBit + postFix, baseStrDiLeptonTriggerBit + postFix); manager.new_cut(baseStrDiLeptonTriggerMatch + postFix, baseStrDiLeptonTriggerMatch + postFix); manager.new_cut(baseStrMllCut + postFix, baseStrMllCut + postFix); manager.new_cut(baseStrMllZVetoCut + postFix, baseStrMllZVetoCut + postFix); @@ -303,6 +314,10 @@ void MuElCategory::evaluate_cuts_post_analyzers(CutManager& manager, const Produ if(m_diLepton.isMuEl) { manager.pass_cut(baseStrCategory + postFix); + // Check that one of the corresponding triggers has fired + if( checkHLT(hlt, HLT::MuonEG) ) + manager.pass_cut(baseStrDiLeptonTriggerBit + postFix); + if(m_diLepton.hlt_idxs.first >= 0 && m_diLepton.hlt_idxs.second >= 0){ // We have fired a trigger. Now, check that it is actually a MuonEG trigger if( checkHLT(hlt, m_diLepton.hlt_idxs.first, m_diLepton.hlt_idxs.second, HLT::MuonEG) ) @@ -378,6 +393,7 @@ void MuMuCategory::register_cuts(CutManager& manager) { manager.new_cut(baseStrCategory + postFix, baseStrCategory + postFix); manager.new_cut(baseStrExtraDiLeptonVeto + postFix, baseStrExtraDiLeptonVeto + postFix); + manager.new_cut(baseStrDiLeptonTriggerBit + postFix, baseStrDiLeptonTriggerBit + postFix); manager.new_cut(baseStrDiLeptonTriggerMatch + postFix, baseStrDiLeptonTriggerMatch + postFix); manager.new_cut(baseStrMllCut + postFix, baseStrMllCut + postFix); manager.new_cut(baseStrMllZVetoCut + postFix, baseStrMllZVetoCut + postFix); @@ -410,6 +426,10 @@ void MuMuCategory::evaluate_cuts_post_analyzers(CutManager& manager, const Produ if(m_diLepton.isMuMu) { manager.pass_cut(baseStrCategory + postFix); + // Check that one of the corresponding triggers has fired + if( checkHLT(hlt, HLT::DoubleMuon) ) + manager.pass_cut(baseStrDiLeptonTriggerBit + postFix); + if(m_diLepton.hlt_idxs.first >= 0 && m_diLepton.hlt_idxs.second >= 0){ // We have fired a trigger. Now, check that it is actually a DoubleMuon trigger if( checkHLT(hlt, m_diLepton.hlt_idxs.first, m_diLepton.hlt_idxs.second, HLT::DoubleMuon) ) diff --git a/src/NeutrinosSolver.cc b/src/NeutrinosSolver.cc index 7f765d2..259187b 100644 --- a/src/NeutrinosSolver.cc +++ b/src/NeutrinosSolver.cc @@ -6,25 +6,23 @@ std::vector