From b765fc7f5af75320161501c99830e6ade2645351 Mon Sep 17 00:00:00 2001 From: Andrea Trapote Date: Mon, 13 Mar 2023 17:55:49 +0100 Subject: [PATCH 1/3] removing non-used muon MVAs --- DataFormats/PatCandidates/interface/Muon.h | 10 -- DataFormats/PatCandidates/src/Muon.cc | 8 - .../PatCandidates/src/classes_def_objects.xml | 3 +- PhysicsTools/NanoAOD/python/muons_cff.py | 2 - PhysicsTools/NanoAOD/python/nanoDQM_cfi.py | 2 - .../PatAlgos/interface/CalculatePtRatioRel.h | 33 ++++ .../PatAlgos/interface/MuonMvaEstimator.h | 45 ----- .../PatAlgos/plugins/PATMuonProducer.cc | 97 +++-------- .../displacedMuonProducer_cff.py | 5 +- .../producersLayer1/muonProducer_cfi.py | 5 +- .../PatAlgos/python/slimming/miniAOD_tools.py | 4 +- .../PatAlgos/src/CalculatePtRatioRel.cc | 92 ++++++++++ PhysicsTools/PatAlgos/src/MuonMvaEstimator.cc | 162 ------------------ 13 files changed, 151 insertions(+), 317 deletions(-) create mode 100644 PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h delete mode 100644 PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h create mode 100644 PhysicsTools/PatAlgos/src/CalculatePtRatioRel.cc delete mode 100644 PhysicsTools/PatAlgos/src/MuonMvaEstimator.cc diff --git a/DataFormats/PatCandidates/interface/Muon.h b/DataFormats/PatCandidates/interface/Muon.h index 5406d7bc92df0..5f34670e14e69 100644 --- a/DataFormats/PatCandidates/interface/Muon.h +++ b/DataFormats/PatCandidates/interface/Muon.h @@ -282,14 +282,6 @@ namespace pat { void setJetPtRatio(float jetPtRatio) { jetPtRatio_ = jetPtRatio; } void setJetPtRel(float jetPtRel) { jetPtRel_ = jetPtRel; } - /// Muon MVA - float mvaValue() const { return mvaValue_; } - void setMvaValue(float mva) { mvaValue_ = mva; } - - // Low pt Muon MVA - float lowptMvaValue() const { return lowptMvaValue_; } - void setLowPtMvaValue(float lowptmva) { lowptMvaValue_ = lowptmva; } - /// Soft Muon MVA float softMvaValue() const { return softMvaValue_; } void setSoftMvaValue(float softmva) { softMvaValue_ = softmva; } @@ -420,8 +412,6 @@ namespace pat { float jetPtRel_; /// Muon MVA - float mvaValue_; - float lowptMvaValue_; float mvaIDValue_; float softMvaValue_; diff --git a/DataFormats/PatCandidates/src/Muon.cc b/DataFormats/PatCandidates/src/Muon.cc index 5a0b5fd4ca960..cfdaa2253d900 100644 --- a/DataFormats/PatCandidates/src/Muon.cc +++ b/DataFormats/PatCandidates/src/Muon.cc @@ -31,8 +31,6 @@ Muon::Muon() pfEcalEnergy_(0), jetPtRatio_(0), jetPtRel_(0), - mvaValue_(0), - lowptMvaValue_(0), mvaIDValue_(0), softMvaValue_(0), inverseBeta_(0), @@ -63,8 +61,6 @@ Muon::Muon(const reco::Muon& aMuon) pfEcalEnergy_(0), jetPtRatio_(0), jetPtRel_(0), - mvaValue_(0), - lowptMvaValue_(0), mvaIDValue_(0), softMvaValue_(0), inverseBeta_(0), @@ -95,8 +91,6 @@ Muon::Muon(const edm::RefToBase& aMuonRef) pfEcalEnergy_(0), jetPtRatio_(0), jetPtRel_(0), - mvaValue_(0), - lowptMvaValue_(0), mvaIDValue_(0), softMvaValue_(0), inverseBeta_(0), @@ -127,8 +121,6 @@ Muon::Muon(const edm::Ptr& aMuonRef) pfEcalEnergy_(0), jetPtRatio_(0), jetPtRel_(0), - mvaValue_(0), - lowptMvaValue_(0), mvaIDValue_(0), softMvaValue_(0), inverseBeta_(0), diff --git a/DataFormats/PatCandidates/src/classes_def_objects.xml b/DataFormats/PatCandidates/src/classes_def_objects.xml index 88fb44593ba8d..a984a6d278381 100644 --- a/DataFormats/PatCandidates/src/classes_def_objects.xml +++ b/DataFormats/PatCandidates/src/classes_def_objects.xml @@ -67,7 +67,8 @@ - + + diff --git a/PhysicsTools/NanoAOD/python/muons_cff.py b/PhysicsTools/NanoAOD/python/muons_cff.py index 23cfb633c79a4..8beaff505fb2c 100644 --- a/PhysicsTools/NanoAOD/python/muons_cff.py +++ b/PhysicsTools/NanoAOD/python/muons_cff.py @@ -162,8 +162,6 @@ highPtId = Var("?passed('CutBasedIdGlobalHighPt')?2:passed('CutBasedIdTrkHighPt')","uint8",doc="high-pT cut-based ID (1 = tracker high pT, 2 = global high pT, which includes tracker high pT)"), pfIsoId = Var("passed('PFIsoVeryLoose')+passed('PFIsoLoose')+passed('PFIsoMedium')+passed('PFIsoTight')+passed('PFIsoVeryTight')+passed('PFIsoVeryVeryTight')","uint8",doc="PFIso ID from miniAOD selector (1=PFIsoVeryLoose, 2=PFIsoLoose, 3=PFIsoMedium, 4=PFIsoTight, 5=PFIsoVeryTight, 6=PFIsoVeryVeryTight)"), tkIsoId = Var("?passed('TkIsoTight')?2:passed('TkIsoLoose')","uint8",doc="TkIso ID (1=TkIsoLoose, 2=TkIsoTight)"), - mvaId = Var("passed('MvaLoose')+passed('MvaMedium')+passed('MvaTight')+passed('MvaVTight')+passed('MvaVVTight')","uint8",doc="Mva for ID of prompt leptons from miniAOD selector (1=MvaLoose, 2=MvaMedium, 3=MvaTight, 4=MvaVTight, 5=MvaVVTight)"), - mvaLowPtId = Var("passed('LowPtMvaLoose')+passed('LowPtMvaMedium')","uint8", doc="Low Pt Mva ID from miniAOD selector (1=LowPtMvaLoose, 2=LowPtMvaMedium)"), miniIsoId = Var("passed('MiniIsoLoose')+passed('MiniIsoMedium')+passed('MiniIsoTight')+passed('MiniIsoVeryTight')","uint8",doc="MiniIso ID from miniAOD selector (1=MiniIsoLoose, 2=MiniIsoMedium, 3=MiniIsoTight, 4=MiniIsoVeryTight)"), mvaMuID = Var("mvaIDValue()",float,doc="MVA-based ID score ",precision=6), mvaMuID_WP = Var("userFloat('mvaIDMuon_wpMedium') + userFloat('mvaIDMuon_wpTight')","uint8",doc="MVA-based ID selector WPs (1=MVAIDwpMedium,2=MVAIDwpTight)"), diff --git a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py index f10a05ebe642c..746f3a640fb00 100644 --- a/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py +++ b/PhysicsTools/NanoAOD/python/nanoDQM_cfi.py @@ -523,9 +523,7 @@ Plot1D('miniPFRelIso_all', 'miniPFRelIso_all', 20, 0, 1, 'mini PF relative isolation, total (with scaled rho*EA PU corrections)'), Plot1D('miniPFRelIso_chg', 'miniPFRelIso_chg', 20, 0, 1, 'mini PF relative isolation, charged component'), Plot1D('multiIsoId', 'multiIsoId', 3, -0.5, 2.5, 'MultiIsoId from miniAOD selector (1=MultiIsoLoose, 2=MultiIsoMedium)'), - Plot1D('mvaId', 'mvaId', 6, -0.5, 5.5, 'Mva ID from miniAOD selector (1=MvaLoose, 2=MvaMedium, 3=MvaTight, 4=MvaVTight, 5=MvaVVTight)'), Plot1D('mvaLowPt', 'mvaLowPt', 20, -1, 1, 'Low pt muon ID score'), - Plot1D('mvaLowPtId', 'mvaLowPtId', 3, -0.5, 2.5, 'Low Pt Mva ID from miniAOD selector (1=LowPtMvaLoose, 2=LowPtMvaMedium)'), Plot1D('mvaTTH', 'mvaTTH', 20, -1, 1, 'TTH MVA lepton ID score'), Plot1D('mvaMuID', 'mvaMuID', 20, 0, 1, 'Score of MVA-based muon ID'), Plot1D('mvaMuID_WP', 'mvaMuID_WP', 3, -0.5, 2.5, 'MVA-based ID selector WPs (1=MVAIDwpMedium,2=MVAIDwpTight)'), diff --git a/PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h b/PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h new file mode 100644 index 0000000000000..a43d9e543c00b --- /dev/null +++ b/PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h @@ -0,0 +1,33 @@ +#ifndef __PhysicsTools_PatAlgos_CalculatePtRatioRel__ +#define __PhysicsTools_PatAlgos_CalculatePtRatioRel__ + +#include "DataFormats/BTauReco/interface/JetTag.h" + +#include +#include + +namespace pat { + class Muon; +} + +namespace reco { + class JetCorrector; +} // namespace reco + +namespace pat { + class CalculatePtRatioRel { + public: + CalculatePtRatioRel(float dRmax); + + ~CalculatePtRatioRel(); + + std::vector computePtRatioRel(const pat::Muon& imuon, + const reco::JetTagCollection& bTags, + const reco::JetCorrector* correctorL1 = nullptr, + const reco::JetCorrector* correctorL1L2L3Res = nullptr) const; + + private: + float dRmax_; + }; +} // namespace pat +#endif diff --git a/PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h b/PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h deleted file mode 100644 index fc7f52edc1c8b..0000000000000 --- a/PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef __PhysicsTools_PatAlgos_MuonMvaEstimator__ -#define __PhysicsTools_PatAlgos_MuonMvaEstimator__ - -#include "DataFormats/BTauReco/interface/JetTag.h" - -#include -#include - -class GBRForest; - -namespace pat { - class Muon; -} - -namespace reco { - class JetCorrector; - class Vertex; -} // namespace reco - -namespace edm { - class FileInPath; -} - -namespace pat { - class MuonMvaEstimator { - public: - MuonMvaEstimator(const edm::FileInPath& weightsfile, float dRmax); - - ~MuonMvaEstimator(); - - float computeMva(const pat::Muon& imuon, - const reco::Vertex& vertex, - const reco::JetTagCollection& bTags, - float& jetPtRatio, - float& jetPtRel, - float& miniIsoValue, - const reco::JetCorrector* correctorL1 = nullptr, - const reco::JetCorrector* correctorL1L2L3Res = nullptr) const; - - private: - std::unique_ptr gbrForest_; - float dRmax_; - }; -} // namespace pat -#endif diff --git a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc index 50605edbd950a..b9dd7ab88f721 100644 --- a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc +++ b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc @@ -42,7 +42,7 @@ #include "PhysicsTools/PatAlgos/interface/EfficiencyLoader.h" #include "PhysicsTools/PatAlgos/interface/KinResolutionsLoader.h" #include "PhysicsTools/PatAlgos/interface/MultiIsolator.h" -#include "PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h" +#include "PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h" #include "PhysicsTools/PatAlgos/interface/MuonMvaIDEstimator.h" #include "PhysicsTools/PatAlgos/interface/PATUserDataHelper.h" #include "PhysicsTools/PatAlgos/interface/SoftMuonMvaEstimator.h" @@ -58,14 +58,12 @@ namespace pat { public: PATMuonHeavyObjectCache(const edm::ParameterSet&); - pat::MuonMvaEstimator const& muonMvaEstimator() const { return *muonMvaEstimator_; } - pat::MuonMvaEstimator const& muonLowPtMvaEstimator() const { return *muonLowPtMvaEstimator_; } + pat::CalculatePtRatioRel const& calculatePtRatioRel() const { return *calculatePtRatioRel_; } pat::MuonMvaIDEstimator const& muonMvaIDEstimator() const { return *muonMvaIDEstimator_; } pat::SoftMuonMvaEstimator const& softMuonMvaEstimator() const { return *softMuonMvaEstimator_; } private: - std::unique_ptr muonLowPtMvaEstimator_; - std::unique_ptr muonMvaEstimator_; + std::unique_ptr calculatePtRatioRel_; std::unique_ptr muonMvaIDEstimator_; std::unique_ptr softMuonMvaEstimator_; }; @@ -101,7 +99,7 @@ namespace pat { typedef std::vector>> IsolationValueMaps; typedef std::pair IsolationLabel; typedef std::vector IsolationLabels; - + /// common muon filling, for both the standard and PF2PAT case void fillMuon(Muon& aMuon, const MuonBaseRef& muonRef, @@ -234,11 +232,10 @@ namespace pat { edm::EDGetTokenT> PUPPINoLeptonsIsolation_neutral_hadrons_; edm::EDGetTokenT> PUPPINoLeptonsIsolation_photons_; /// standard muon selectors - bool computeMuonMVA_; bool computeMuonIDMVA_; bool computeSoftMuonMVA_; bool recomputeBasicSelectors_; - bool mvaUseJec_; + bool useJec_; edm::EDGetTokenT mvaBTagCollectionTag_; edm::EDGetTokenT mvaL1Corrector_; edm::EDGetTokenT mvaL1L2L3ResCorrector_; @@ -325,12 +322,9 @@ using namespace pat; using namespace std; PATMuonHeavyObjectCache::PATMuonHeavyObjectCache(const edm::ParameterSet& iConfig) { - if (iConfig.getParameter("computeMuonMVA")) { - edm::FileInPath mvaTrainingFile = iConfig.getParameter("mvaTrainingFile"); - edm::FileInPath mvaLowPtTrainingFile = iConfig.getParameter("lowPtmvaTrainingFile"); + if (iConfig.getParameter("computeMiniIso")) { float mvaDrMax = iConfig.getParameter("mvaDrMax"); - muonMvaEstimator_ = std::make_unique(mvaTrainingFile, mvaDrMax); - muonLowPtMvaEstimator_ = std::make_unique(mvaLowPtTrainingFile, mvaDrMax); + calculatePtRatioRel_ = std::make_unique(mvaDrMax); } if (iConfig.getParameter("computeMuonIDMVA")) { @@ -347,11 +341,10 @@ PATMuonHeavyObjectCache::PATMuonHeavyObjectCache(const edm::ParameterSet& iConfi PATMuonProducer::PATMuonProducer(const edm::ParameterSet& iConfig, PATMuonHeavyObjectCache const*) : relMiniIsoPUCorrected_(0), useUserData_(iConfig.exists("userData")), - computeMuonMVA_(false), computeMuonIDMVA_(false), computeSoftMuonMVA_(false), recomputeBasicSelectors_(false), - mvaUseJec_(false), + useJec_(false), isolator_(iConfig.getParameter("userIsolation"), consumesCollector(), false), geometryToken_{esConsumes()}, transientTrackBuilderToken_{esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))}, @@ -460,18 +453,15 @@ PATMuonProducer::PATMuonProducer(const edm::ParameterSet& iConfig, PATMuonHeavyO // standard selectors recomputeBasicSelectors_ = iConfig.getParameter("recomputeBasicSelectors"); - computeMuonMVA_ = iConfig.getParameter("computeMuonMVA"); computeMuonIDMVA_ = iConfig.getParameter("computeMuonIDMVA"); - if (computeMuonMVA_ and not computeMiniIso_) - throw cms::Exception("ConfigurationError") << "MiniIso is needed for Muon MVA calculation.\n"; - if (computeMuonMVA_) { + if (computeMiniIso_) { // pfCombinedInclusiveSecondaryVertexV2BJetTags mvaBTagCollectionTag_ = consumes(iConfig.getParameter("mvaJetTag")); mvaL1Corrector_ = consumes(iConfig.getParameter("mvaL1Corrector")); mvaL1L2L3ResCorrector_ = consumes(iConfig.getParameter("mvaL1L2L3ResCorrector")); rho_ = consumes(iConfig.getParameter("rho")); - mvaUseJec_ = iConfig.getParameter("mvaUseJec"); + useJec_ = iConfig.getParameter("useJec"); } computeSoftMuonMVA_ = iConfig.getParameter("computeSoftMuonMVA"); @@ -642,7 +632,7 @@ void PATMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) edm::Handle mvaBTagCollectionTag; edm::Handle mvaL1Corrector; edm::Handle mvaL1L2L3ResCorrector; - if (computeMuonMVA_) { + if (computeMiniIso_) { iEvent.getByToken(mvaBTagCollectionTag_, mvaBTagCollectionTag); iEvent.getByToken(mvaL1Corrector_, mvaL1Corrector); iEvent.getByToken(mvaL1L2L3ResCorrector_, mvaL1L2L3ResCorrector); @@ -932,7 +922,7 @@ void PATMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) // Need a separate loop over muons to have all inputs properly // computed and stored in the object. edm::Handle rho; - if (computeMuonMVA_) + if (computeMiniIso_) iEvent.getByToken(rho_, rho); const reco::Vertex* pv(nullptr); if (primaryVertexIsValid) @@ -980,68 +970,23 @@ void PATMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) muon.setSelector(reco::Muon::PuppiIsoTight, puppiCombinedIsolationPAT < 0.12); } - float jetPtRatio = 0.0; - float jetPtRel = 0.0; - float mva = 0.0; - float mva_lowpt = 0.0; - if (computeMuonMVA_ && primaryVertexIsValid && computeMiniIso_) { - if (mvaUseJec_) { - mva = globalCache()->muonMvaEstimator().computeMva(muon, - primaryVertex, - *(mvaBTagCollectionTag.product()), - jetPtRatio, - jetPtRel, - miniIsoValue, - mvaL1Corrector.product(), - mvaL1L2L3ResCorrector.product()); - mva_lowpt = globalCache()->muonLowPtMvaEstimator().computeMva(muon, - primaryVertex, - *(mvaBTagCollectionTag.product()), - jetPtRatio, - jetPtRel, - miniIsoValue, - mvaL1Corrector.product(), - mvaL1L2L3ResCorrector.product()); - + std::vector jetPtRatioRel = {0.0, 0.0}; + if (primaryVertexIsValid && computeMiniIso_) { + if (useJec_) { + jetPtRatioRel = globalCache()->calculatePtRatioRel().computePtRatioRel( + muon, *(mvaBTagCollectionTag.product()), mvaL1Corrector.product(), mvaL1L2L3ResCorrector.product()); } else { - mva = globalCache()->muonMvaEstimator().computeMva( - muon, primaryVertex, *mvaBTagCollectionTag, jetPtRatio, jetPtRel, miniIsoValue); - mva_lowpt = globalCache()->muonLowPtMvaEstimator().computeMva( - muon, primaryVertex, *mvaBTagCollectionTag, jetPtRatio, jetPtRel, miniIsoValue); + jetPtRatioRel = globalCache()->calculatePtRatioRel().computePtRatioRel(muon, *mvaBTagCollectionTag); } - muon.setMvaValue(mva); - muon.setLowPtMvaValue(mva_lowpt); - muon.setJetPtRatio(jetPtRatio); - muon.setJetPtRel(jetPtRel); - + muon.setJetPtRatio(jetPtRatioRel[0]); + muon.setJetPtRel(jetPtRatioRel[1]); + // multi-isolation if (computeMiniIso_) { muon.setSelector(reco::Muon::MultiIsoMedium, miniIsoValue < 0.11 && (muon.jetPtRatio() > 0.74 || muon.jetPtRel() > 6.8)); } - - // MVA working points - // https://twiki.cern.ch/twiki/bin/viewauth/CMS/LeptonMVA - const double dB2D = std::abs(muon.dB(pat::Muon::PV2D)); - const double dB3D = std::abs(muon.dB(pat::Muon::PV3D)); - const double edB3D = std::abs(muon.edB(pat::Muon::PV3D)); - const double sip3D = edB3D > 0 ? dB3D / edB3D : 0.0; - const double dz = std::abs(muon.muonBestTrack()->dz(primaryVertex.position())); - - // muon preselection - if (muon.pt() > 5 and muon.isLooseMuon() and muon.passed(reco::Muon::MiniIsoLoose) and sip3D < 8.0 and - dB2D < 0.05 and dz < 0.1) { - muon.setSelector(reco::Muon::MvaLoose, muon.mvaValue() > -0.60); - muon.setSelector(reco::Muon::MvaMedium, muon.mvaValue() > -0.20); - muon.setSelector(reco::Muon::MvaTight, muon.mvaValue() > 0.15); - muon.setSelector(reco::Muon::MvaVTight, muon.mvaValue() > 0.45); - muon.setSelector(reco::Muon::MvaVVTight, muon.mvaValue() > 0.9); - } - if (muon.pt() > 5 and muon.isLooseMuon() and sip3D < 4 and dB2D < 0.5 and dz < 1) { - muon.setSelector(reco::Muon::LowPtMvaLoose, muon.lowptMvaValue() > -0.60); - muon.setSelector(reco::Muon::LowPtMvaMedium, muon.lowptMvaValue() > -0.20); - } } // MVA ID diff --git a/PhysicsTools/PatAlgos/python/producersLayer1/displacedMuonProducer_cff.py b/PhysicsTools/PatAlgos/python/producersLayer1/displacedMuonProducer_cff.py index 4c160f238a501..6ddc7ee47ceff 100644 --- a/PhysicsTools/PatAlgos/python/producersLayer1/displacedMuonProducer_cff.py +++ b/PhysicsTools/PatAlgos/python/producersLayer1/displacedMuonProducer_cff.py @@ -60,11 +60,8 @@ # Standard Muon Selectors and Jet-related observables # Depends on MiniIsolation, so only works in miniaod # Don't forget to set flags properly in miniAOD_tools.py - computeMuonMVA = False, - mvaTrainingFile = "RecoMuon/MuonIdentification/data/mu_2017_BDTG.weights.xml", - lowPtmvaTrainingFile = "RecoMuon/MuonIdentification/data/mu_lowpt_BDTG.weights.xml", recomputeBasicSelectors = False, - mvaUseJec = False, + useJec = False, mvaDrMax = 0.4, mvaJetTag = "pfCombinedInclusiveSecondaryVertexV2BJetTags", mvaL1Corrector = "ak4PFCHSL1FastjetCorrector", diff --git a/PhysicsTools/PatAlgos/python/producersLayer1/muonProducer_cfi.py b/PhysicsTools/PatAlgos/python/producersLayer1/muonProducer_cfi.py index f95b31dd5a4d4..eb562b4b51340 100644 --- a/PhysicsTools/PatAlgos/python/producersLayer1/muonProducer_cfi.py +++ b/PhysicsTools/PatAlgos/python/producersLayer1/muonProducer_cfi.py @@ -109,13 +109,10 @@ # Standard Muon Selectors and Jet-related observables # Depends on MiniIsolation, so only works in miniaod # Don't forget to set flags properly in miniAOD_tools.py - computeMuonMVA = cms.bool(False), computeMuonIDMVA = cms.bool(False), - mvaTrainingFile = cms.FileInPath("RecoMuon/MuonIdentification/data/mu_2017_BDTG.weights.xml"), mvaIDTrainingFile = cms.FileInPath("RecoMuon/MuonIdentification/data/mvaID.onnx"), - lowPtmvaTrainingFile = cms.FileInPath("RecoMuon/MuonIdentification/data/mu_lowpt_BDTG.weights.xml"), recomputeBasicSelectors = cms.bool(True), - mvaUseJec = cms.bool(True), + useJec = cms.bool(True), mvaDrMax = cms.double(0.4), mvaJetTag = cms.InputTag("pfCombinedInclusiveSecondaryVertexV2BJetTags"), mvaL1Corrector = cms.InputTag("ak4PFCHSL1FastjetCorrector"), diff --git a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py index bff1bfb7ed5d8..669db367cfa86 100644 --- a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py +++ b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py @@ -30,7 +30,6 @@ def miniAOD_customizeCommon(process): process.patMuons.puppiNoLeptonsIsolationPhotons = cms.InputTag("muonPUPPINoLeptonsIsolation","gamma-DR040-ThresholdVeto000-ConeVeto001") process.patMuons.computeMiniIso = True - process.patMuons.computeMuonMVA = True process.patMuons.computeMuonIDMVA = True process.patMuons.computeSoftMuonMVA = True @@ -41,8 +40,7 @@ def miniAOD_customizeCommon(process): run2_muon_2016.toModify( process.patMuons, effectiveAreaVec = [0.0735,0.0619,0.0465,0.0433,0.0577]) run2_muon_2017.toModify( process.patMuons, effectiveAreaVec = [0.0566, 0.0562, 0.0363, 0.0119, 0.0064]) run2_muon_2018.toModify( process.patMuons, effectiveAreaVec = [0.0566, 0.0562, 0.0363, 0.0119, 0.0064]) - run2_muon_2016.toModify( process.patMuons, mvaTrainingFile = "RecoMuon/MuonIdentification/data/mu_2016_BDTG.weights.xml") - + process.patMuons.computePuppiCombinedIso = True # # disable embedding of electron and photon associated objects already stored by the ReducedEGProducer diff --git a/PhysicsTools/PatAlgos/src/CalculatePtRatioRel.cc b/PhysicsTools/PatAlgos/src/CalculatePtRatioRel.cc new file mode 100644 index 0000000000000..79e6f18986b71 --- /dev/null +++ b/PhysicsTools/PatAlgos/src/CalculatePtRatioRel.cc @@ -0,0 +1,92 @@ +#include "PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h" + +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "CommonTools/MVAUtils/interface/GBRForestTools.h" +#include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/MuonReco/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonSelectors.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/JetReco/interface/PFJet.h" +#include "DataFormats/JetReco/interface/PFJetCollection.h" +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "JetMETCorrections/JetCorrector/interface/JetCorrector.h" + +using namespace pat; + +CalculatePtRatioRel::CalculatePtRatioRel(float dRmax) : dRmax_(dRmax) {} + +CalculatePtRatioRel::~CalculatePtRatioRel() {} + +namespace { + + float ptRel(const reco::Candidate::LorentzVector& muP4, + const reco::Candidate::LorentzVector& jetP4, + bool subtractMuon = true) { + reco::Candidate::LorentzVector jp4 = jetP4; + if (subtractMuon) + jp4 -= muP4; + float dot = muP4.Vect().Dot(jp4.Vect()); + float ptrel = muP4.P2() - dot * dot / jp4.P2(); + ptrel = ptrel > 0 ? sqrt(ptrel) : 0.0; + return ptrel; + } +} // namespace + +std::vector CalculatePtRatioRel::computePtRatioRel(const pat::Muon& muon, + const reco::JetTagCollection& bTags, + const reco::JetCorrector* correctorL1, + const reco::JetCorrector* correctorL1L2L3Res) const { + //Initialise loop variables + double minDr = 9999; + double jecL1L2L3Res = 1.; + double jecL1 = 1.; + float JetPtRatio = 0.0; + float JetPtRel = 0.0; + + // Compute corrected isolation variables + double chIso = muon.pfIsolationR04().sumChargedHadronPt; + double nIso = muon.pfIsolationR04().sumNeutralHadronEt; + double phoIso = muon.pfIsolationR04().sumPhotonEt; + double puIso = muon.pfIsolationR04().sumPUPt; + double dbCorrectedIsolation = chIso + std::max(nIso + phoIso - .5 * puIso, 0.); + double dbCorrectedRelIso = dbCorrectedIsolation / muon.pt(); + + JetPtRatio = 1. / (1 + dbCorrectedRelIso); + JetPtRel = 0; + + for (const auto& tagI : bTags) { + // for each muon with the lepton + double dr = deltaR(*(tagI.first), muon); + if (dr > minDr) + continue; + minDr = dr; + + const reco::Candidate::LorentzVector& muP4(muon.p4()); + reco::Candidate::LorentzVector jetP4(tagI.first->p4()); + + if (correctorL1 && correctorL1L2L3Res) { + jecL1L2L3Res = correctorL1L2L3Res->correction(*(tagI.first)); + jecL1 = correctorL1->correction(*(tagI.first)); + } + + if (minDr < dRmax_) { + if ((jetP4 - muP4).Rho() < 0.0001) { + JetPtRel = 0; + JetPtRatio = 1; + } else { + jetP4 -= muP4 / jecL1; + jetP4 *= jecL1L2L3Res; + jetP4 += muP4; + + JetPtRatio = muP4.pt() / jetP4.pt(); + JetPtRel = ptRel(muP4, jetP4); + } + } + } + + if (JetPtRatio > 1.5) + JetPtRatio = 1.5; + + std::vector outputs = {JetPtRatio, JetPtRel}; + return outputs; +}; diff --git a/PhysicsTools/PatAlgos/src/MuonMvaEstimator.cc b/PhysicsTools/PatAlgos/src/MuonMvaEstimator.cc deleted file mode 100644 index 5bd3c4463c344..0000000000000 --- a/PhysicsTools/PatAlgos/src/MuonMvaEstimator.cc +++ /dev/null @@ -1,162 +0,0 @@ -#include "PhysicsTools/PatAlgos/interface/MuonMvaEstimator.h" - -#include "FWCore/ParameterSet/interface/FileInPath.h" -#include "CommonTools/MVAUtils/interface/GBRForestTools.h" -#include "DataFormats/Candidate/interface/Candidate.h" -#include "DataFormats/MuonReco/interface/Muon.h" -#include "DataFormats/MuonReco/interface/MuonSelectors.h" -#include "DataFormats/VertexReco/interface/Vertex.h" -#include "DataFormats/JetReco/interface/PFJet.h" -#include "DataFormats/JetReco/interface/PFJetCollection.h" -#include "DataFormats/PatCandidates/interface/Muon.h" -#include "JetMETCorrections/JetCorrector/interface/JetCorrector.h" - -using namespace pat; - -MuonMvaEstimator::MuonMvaEstimator(const edm::FileInPath& weightsfile, float dRmax) : dRmax_(dRmax) { - gbrForest_ = createGBRForest(weightsfile); -} - -MuonMvaEstimator::~MuonMvaEstimator() {} - -namespace { - - enum inputIndexes { - kPt, - kEta, - kJetNDauCharged, - kMiniRelIsoCharged, - kMiniRelIsoNeutral, - kJetPtRel, - kJetBTagCSV, - kJetPtRatio, - kLog_abs_dxyBS, - kSip, - kLog_abs_dzPV, - kSegmentCompatibility, - kLast - }; - - float ptRel(const reco::Candidate::LorentzVector& muP4, - const reco::Candidate::LorentzVector& jetP4, - bool subtractMuon = true) { - reco::Candidate::LorentzVector jp4 = jetP4; - if (subtractMuon) - jp4 -= muP4; - float dot = muP4.Vect().Dot(jp4.Vect()); - float ptrel = muP4.P2() - dot * dot / jp4.P2(); - ptrel = ptrel > 0 ? sqrt(ptrel) : 0.0; - return ptrel; - } -} // namespace - -float MuonMvaEstimator::computeMva(const pat::Muon& muon, - const reco::Vertex& vertex, - const reco::JetTagCollection& bTags, - float& jetPtRatio, - float& jetPtRel, - float& miniIsoValue, - const reco::JetCorrector* correctorL1, - const reco::JetCorrector* correctorL1L2L3Res) const { - float var[kLast]{}; - - var[kPt] = muon.pt(); - var[kEta] = muon.eta(); - var[kSegmentCompatibility] = muon.segmentCompatibility(); - var[kMiniRelIsoCharged] = muon.miniPFIsolation().chargedHadronIso() / muon.pt(); - var[kMiniRelIsoNeutral] = miniIsoValue - var[kMiniRelIsoCharged]; - - double dB2D = fabs(muon.dB(pat::Muon::PV2D)); - double dB3D = muon.dB(pat::Muon::PV3D); - double edB3D = muon.edB(pat::Muon::PV3D); - double dz = fabs(muon.muonBestTrack()->dz(vertex.position())); - var[kSip] = edB3D > 0 ? fabs(dB3D / edB3D) : 0.0; - var[kLog_abs_dxyBS] = dB2D > 0 ? log(dB2D) : 0; - var[kLog_abs_dzPV] = dz > 0 ? log(dz) : 0; - - //Initialise loop variables - double minDr = 9999; - double jecL1L2L3Res = 1.; - double jecL1 = 1.; - - // Compute corrected isolation variables - double chIso = muon.pfIsolationR04().sumChargedHadronPt; - double nIso = muon.pfIsolationR04().sumNeutralHadronEt; - double phoIso = muon.pfIsolationR04().sumPhotonEt; - double puIso = muon.pfIsolationR04().sumPUPt; - double dbCorrectedIsolation = chIso + std::max(nIso + phoIso - .5 * puIso, 0.); - double dbCorrectedRelIso = dbCorrectedIsolation / muon.pt(); - - var[kJetPtRatio] = 1. / (1 + dbCorrectedRelIso); - var[kJetPtRel] = 0; - var[kJetBTagCSV] = -999; - var[kJetNDauCharged] = -1; - - for (const auto& tagI : bTags) { - // for each muon with the lepton - double dr = deltaR(*(tagI.first), muon); - if (dr > minDr) - continue; - minDr = dr; - - const reco::Candidate::LorentzVector& muP4(muon.p4()); - reco::Candidate::LorentzVector jetP4(tagI.first->p4()); - - if (correctorL1 && correctorL1L2L3Res) { - jecL1L2L3Res = correctorL1L2L3Res->correction(*(tagI.first)); - jecL1 = correctorL1->correction(*(tagI.first)); - } - - // Get b-jet info - var[kJetBTagCSV] = tagI.second; - var[kJetNDauCharged] = 0; - for (auto jet : tagI.first->getJetConstituentsQuick()) { - const reco::PFCandidate* pfcand = dynamic_cast(jet); - if (pfcand == nullptr) - throw cms::Exception("ConfigurationError") << "Cannot get jet constituents"; - if (pfcand->charge() == 0) - continue; - auto bestTrackPtr = pfcand->bestTrack(); - if (!bestTrackPtr) - continue; - if (!bestTrackPtr->quality(reco::Track::highPurity)) - continue; - if (bestTrackPtr->pt() < 1.) - continue; - if (bestTrackPtr->hitPattern().numberOfValidHits() < 8) - continue; - if (bestTrackPtr->hitPattern().numberOfValidPixelHits() < 2) - continue; - if (bestTrackPtr->normalizedChi2() >= 5) - continue; - - if (std::fabs(bestTrackPtr->dxy(vertex.position())) > 0.2) - continue; - if (std::fabs(bestTrackPtr->dz(vertex.position())) > 17) - continue; - var[kJetNDauCharged]++; - } - - if (minDr < dRmax_) { - if ((jetP4 - muP4).Rho() < 0.0001) { - var[kJetPtRel] = 0; - var[kJetPtRatio] = 1; - } else { - jetP4 -= muP4 / jecL1; - jetP4 *= jecL1L2L3Res; - jetP4 += muP4; - - var[kJetPtRatio] = muP4.pt() / jetP4.pt(); - var[kJetPtRel] = ptRel(muP4, jetP4); - } - } - } - - if (var[kJetPtRatio] > 1.5) - var[kJetPtRatio] = 1.5; - if (var[kJetBTagCSV] < 0) - var[kJetBTagCSV] = 0; - jetPtRatio = var[kJetPtRatio]; - jetPtRel = var[kJetPtRel]; - return gbrForest_->GetClassifier(var); -}; From 5f24cbe31b1b5ee8670494d3e9032b8c7abaa790 Mon Sep 17 00:00:00 2001 From: Andrea Trapote Date: Wed, 22 Mar 2023 12:08:42 +0100 Subject: [PATCH 2/3] implementing A. Perrota comments --- .../PatAlgos/interface/CalculatePtRatioRel.h | 12 ++--- .../PatAlgos/plugins/PATMuonProducer.cc | 6 +-- .../PatAlgos/src/CalculatePtRatioRel.cc | 52 ++++++++----------- 3 files changed, 31 insertions(+), 39 deletions(-) diff --git a/PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h b/PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h index a43d9e543c00b..518e764c87d39 100644 --- a/PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h +++ b/PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h @@ -17,17 +17,17 @@ namespace reco { namespace pat { class CalculatePtRatioRel { public: - CalculatePtRatioRel(float dRmax); + CalculatePtRatioRel(float dR2max); ~CalculatePtRatioRel(); - std::vector computePtRatioRel(const pat::Muon& imuon, - const reco::JetTagCollection& bTags, - const reco::JetCorrector* correctorL1 = nullptr, - const reco::JetCorrector* correctorL1L2L3Res = nullptr) const; + std::array computePtRatioRel(const pat::Muon& imuon, + const reco::JetTagCollection& bTags, + const reco::JetCorrector* correctorL1 = nullptr, + const reco::JetCorrector* correctorL1L2L3Res = nullptr) const; private: - float dRmax_; + float dR2max_; }; } // namespace pat #endif diff --git a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc index b9dd7ab88f721..4a18c5feec5c4 100644 --- a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc +++ b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc @@ -99,7 +99,7 @@ namespace pat { typedef std::vector>> IsolationValueMaps; typedef std::pair IsolationLabel; typedef std::vector IsolationLabels; - + /// common muon filling, for both the standard and PF2PAT case void fillMuon(Muon& aMuon, const MuonBaseRef& muonRef, @@ -970,7 +970,7 @@ void PATMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) muon.setSelector(reco::Muon::PuppiIsoTight, puppiCombinedIsolationPAT < 0.12); } - std::vector jetPtRatioRel = {0.0, 0.0}; + std::array jetPtRatioRel = {{0.0, 0.0}}; if (primaryVertexIsValid && computeMiniIso_) { if (useJec_) { jetPtRatioRel = globalCache()->calculatePtRatioRel().computePtRatioRel( @@ -981,7 +981,7 @@ void PATMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) muon.setJetPtRatio(jetPtRatioRel[0]); muon.setJetPtRel(jetPtRatioRel[1]); - + // multi-isolation if (computeMiniIso_) { muon.setSelector(reco::Muon::MultiIsoMedium, diff --git a/PhysicsTools/PatAlgos/src/CalculatePtRatioRel.cc b/PhysicsTools/PatAlgos/src/CalculatePtRatioRel.cc index 79e6f18986b71..fa2c85905947a 100644 --- a/PhysicsTools/PatAlgos/src/CalculatePtRatioRel.cc +++ b/PhysicsTools/PatAlgos/src/CalculatePtRatioRel.cc @@ -1,11 +1,5 @@ #include "PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h" -#include "FWCore/ParameterSet/interface/FileInPath.h" -#include "CommonTools/MVAUtils/interface/GBRForestTools.h" -#include "DataFormats/Candidate/interface/Candidate.h" -#include "DataFormats/MuonReco/interface/Muon.h" -#include "DataFormats/MuonReco/interface/MuonSelectors.h" -#include "DataFormats/VertexReco/interface/Vertex.h" #include "DataFormats/JetReco/interface/PFJet.h" #include "DataFormats/JetReco/interface/PFJetCollection.h" #include "DataFormats/PatCandidates/interface/Muon.h" @@ -13,7 +7,7 @@ using namespace pat; -CalculatePtRatioRel::CalculatePtRatioRel(float dRmax) : dRmax_(dRmax) {} +CalculatePtRatioRel::CalculatePtRatioRel(float dR2max) : dR2max_(dR2max) {} CalculatePtRatioRel::~CalculatePtRatioRel() {} @@ -32,36 +26,34 @@ namespace { } } // namespace -std::vector CalculatePtRatioRel::computePtRatioRel(const pat::Muon& muon, - const reco::JetTagCollection& bTags, - const reco::JetCorrector* correctorL1, - const reco::JetCorrector* correctorL1L2L3Res) const { +std::array CalculatePtRatioRel::computePtRatioRel(const pat::Muon& muon, + const reco::JetTagCollection& bTags, + const reco::JetCorrector* correctorL1, + const reco::JetCorrector* correctorL1L2L3Res) const { //Initialise loop variables - double minDr = 9999; + float minDr2 = 9999; double jecL1L2L3Res = 1.; double jecL1 = 1.; - float JetPtRatio = 0.0; - float JetPtRel = 0.0; // Compute corrected isolation variables - double chIso = muon.pfIsolationR04().sumChargedHadronPt; - double nIso = muon.pfIsolationR04().sumNeutralHadronEt; - double phoIso = muon.pfIsolationR04().sumPhotonEt; - double puIso = muon.pfIsolationR04().sumPUPt; - double dbCorrectedIsolation = chIso + std::max(nIso + phoIso - .5 * puIso, 0.); - double dbCorrectedRelIso = dbCorrectedIsolation / muon.pt(); + const float chIso = muon.pfIsolationR04().sumChargedHadronPt; + const float nIso = muon.pfIsolationR04().sumNeutralHadronEt; + const float phoIso = muon.pfIsolationR04().sumPhotonEt; + const float puIso = muon.pfIsolationR04().sumPUPt; + const float dbCorrectedIsolation = chIso + std::max(nIso + phoIso - .5f * puIso, 0.f); + const float dbCorrectedRelIso = dbCorrectedIsolation / muon.pt(); - JetPtRatio = 1. / (1 + dbCorrectedRelIso); - JetPtRel = 0; + float JetPtRatio = 1. / (1 + dbCorrectedRelIso); + float JetPtRel = 0.; + const reco::Candidate::LorentzVector& muP4(muon.p4()); for (const auto& tagI : bTags) { // for each muon with the lepton - double dr = deltaR(*(tagI.first), muon); - if (dr > minDr) + float dr2 = deltaR2(*(tagI.first), muon); + if (dr2 > minDr2) continue; - minDr = dr; + minDr2 = dr2; - const reco::Candidate::LorentzVector& muP4(muon.p4()); reco::Candidate::LorentzVector jetP4(tagI.first->p4()); if (correctorL1 && correctorL1L2L3Res) { @@ -69,10 +61,10 @@ std::vector CalculatePtRatioRel::computePtRatioRel(const pat::Muon& muon, jecL1 = correctorL1->correction(*(tagI.first)); } - if (minDr < dRmax_) { + if (minDr2 < dR2max_) { if ((jetP4 - muP4).Rho() < 0.0001) { - JetPtRel = 0; - JetPtRatio = 1; + JetPtRel = 0.; + JetPtRatio = 1.; } else { jetP4 -= muP4 / jecL1; jetP4 *= jecL1L2L3Res; @@ -87,6 +79,6 @@ std::vector CalculatePtRatioRel::computePtRatioRel(const pat::Muon& muon, if (JetPtRatio > 1.5) JetPtRatio = 1.5; - std::vector outputs = {JetPtRatio, JetPtRel}; + std::array outputs = {{JetPtRatio, JetPtRel}}; return outputs; }; From e4ef45ec305fcbb9b9a1a7de3a08749fc86c761d Mon Sep 17 00:00:00 2001 From: Andrea Trapote Date: Wed, 22 Mar 2023 12:45:58 +0100 Subject: [PATCH 3/3] updating mvaDRMax * mvaDRMax --- PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc index 4a18c5feec5c4..9e3244ecce841 100644 --- a/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc +++ b/PhysicsTools/PatAlgos/plugins/PATMuonProducer.cc @@ -324,7 +324,7 @@ using namespace std; PATMuonHeavyObjectCache::PATMuonHeavyObjectCache(const edm::ParameterSet& iConfig) { if (iConfig.getParameter("computeMiniIso")) { float mvaDrMax = iConfig.getParameter("mvaDrMax"); - calculatePtRatioRel_ = std::make_unique(mvaDrMax); + calculatePtRatioRel_ = std::make_unique(mvaDrMax * mvaDrMax); } if (iConfig.getParameter("computeMuonIDMVA")) {