From 875eb9d65c1770e572ecd5b5b70888d837d0c3e1 Mon Sep 17 00:00:00 2001 From: Nicola Amapane Date: Mon, 24 Jul 2023 18:05:39 +0200 Subject: [PATCH] Add muon beamspot constrained pT. Increase muon eta and phi precision to 16 bits --- PhysicsTools/NanoAOD/plugins/BuildFile.xml | 1 + .../MuonBeamspotConstraintValueMapProducer.cc | 141 ++++++++++++++++++ PhysicsTools/NanoAOD/python/muons_cff.py | 16 +- 3 files changed, 154 insertions(+), 4 deletions(-) create mode 100644 PhysicsTools/NanoAOD/plugins/MuonBeamspotConstraintValueMapProducer.cc diff --git a/PhysicsTools/NanoAOD/plugins/BuildFile.xml b/PhysicsTools/NanoAOD/plugins/BuildFile.xml index aaffe2b44b891..032ad4a778ec9 100644 --- a/PhysicsTools/NanoAOD/plugins/BuildFile.xml +++ b/PhysicsTools/NanoAOD/plugins/BuildFile.xml @@ -24,6 +24,7 @@ + diff --git a/PhysicsTools/NanoAOD/plugins/MuonBeamspotConstraintValueMapProducer.cc b/PhysicsTools/NanoAOD/plugins/MuonBeamspotConstraintValueMapProducer.cc new file mode 100644 index 0000000000000..b20216cdaeb42 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/MuonBeamspotConstraintValueMapProducer.cc @@ -0,0 +1,141 @@ +/** \class MuonBeamspotConstraintValueMapProducer + * Compute muon pt and ptErr after beamspot constraint. + * + */ +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/Event.h" +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/TransientTrack/interface/TransientTrack.h" +#include "RecoVertex/KalmanVertexFit/interface/SingleTrackVertexConstraint.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +class MuonBeamspotConstraintValueMapProducer : public edm::global::EDProducer<> { +public: + explicit MuonBeamspotConstraintValueMapProducer(const edm::ParameterSet& config) + : muonToken_(consumes(config.getParameter("src"))), + beamSpotToken_(consumes(config.getParameter("beamspot"))), + PrimaryVertexToken_(consumes(config.getParameter("vertices"))), + ttbToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) { + produces>("muonBSConstrainedPt"); + produces>("muonBSConstrainedPtErr"); + produces>("muonBSConstrainedChi2"); + } + + ~MuonBeamspotConstraintValueMapProducer() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("src", edm::InputTag("muons"))->setComment("Muon collection"); + desc.add("beamspot", edm::InputTag("offlineBeamSpot"))->setComment("Beam spot collection"); + desc.add("vertices", edm::InputTag("offlineSlimmedPrimaryVertices")) + ->setComment("Primary vertex collection"); + + descriptions.addWithDefaultLabel(desc); + } + +private: + void produce(edm::StreamID streamID, edm::Event& event, const edm::EventSetup& setup) const override { + edm::Handle muons; + event.getByToken(muonToken_, muons); + + edm::Handle beamSpotHandle; + event.getByToken(beamSpotToken_, beamSpotHandle); + + edm::ESHandle ttkb = setup.getHandle(ttbToken_); + + std::vector pts, ptErrs, chi2s; + pts.reserve(muons->size()); + ptErrs.reserve(muons->size()); + chi2s.reserve(muons->size()); + + for (const auto& muon : *muons) { + bool tbd = true; + if (beamSpotHandle.isValid()) { + double BeamWidthX = beamSpotHandle->BeamWidthX(); + double BeamWidthXError = beamSpotHandle->BeamWidthXError(); + double BeamWidthY = beamSpotHandle->BeamWidthY(); + double BeamWidthYError = beamSpotHandle->BeamWidthYError(); + // Protect for mis-reconstructed beamspots (note that + // SingleTrackVertexConstraint uses the width for the constraint, + // not the error) + if ((BeamWidthXError / BeamWidthX < 0.3) && (BeamWidthYError / BeamWidthY < 0.3)) { + SingleTrackVertexConstraint::BTFtuple btft = + stvc.constrain(ttkb->build(muon.muonBestTrack()), *beamSpotHandle); + if (std::get<0>(btft)) { + const reco::Track& trkBS = std::get<1>(btft).track(); + pts.push_back(trkBS.pt()); + ptErrs.push_back(trkBS.ptError()); + chi2s.push_back(std::get<2>(btft)); + tbd = false; + } + } + } + + if (tbd) { + // Invalid BS; use PV instead + edm::Handle pvHandle; + event.getByToken(PrimaryVertexToken_, pvHandle); + + if (pvHandle.isValid() && !pvHandle->empty()) { + auto pv = pvHandle->at(0); + VertexState pvs = VertexState(GlobalPoint(Basic3DVector(pv.position())), GlobalError(pv.covariance())); + + SingleTrackVertexConstraint::BTFtuple btft = stvc.constrain(ttkb->build(muon.muonBestTrack()), pvs); + if (std::get<0>(btft)) { + // chi2 = std::get<2>(btft)); // should apply a cut, or store this as well? + const reco::Track& trkBS = std::get<1>(btft).track(); + pts.push_back(trkBS.pt()); + ptErrs.push_back(trkBS.ptError()); + chi2s.push_back(std::get<2>(btft)); + tbd = false; + } + } + } + + if (tbd) { + // Fall-back case, keep the unconstrained values + pts.push_back(muon.pt()); + ptErrs.push_back(muon.bestTrack()->ptError()); + chi2s.push_back(-1.f); + } + } + + { + std::unique_ptr> valueMap(new edm::ValueMap()); + edm::ValueMap::Filler filler(*valueMap); + filler.insert(muons, pts.begin(), pts.end()); + filler.fill(); + event.put(std::move(valueMap), "muonBSConstrainedPt"); + } + + { + std::unique_ptr> valueMap(new edm::ValueMap()); + edm::ValueMap::Filler filler(*valueMap); + filler.insert(muons, ptErrs.begin(), ptErrs.end()); + filler.fill(); + event.put(std::move(valueMap), "muonBSConstrainedPtErr"); + } + + { + std::unique_ptr> valueMap(new edm::ValueMap()); + edm::ValueMap::Filler filler(*valueMap); + filler.insert(muons, chi2s.begin(), chi2s.end()); + filler.fill(); + event.put(std::move(valueMap), "muonBSConstrainedChi2"); + } + } + + edm::EDGetTokenT muonToken_; + edm::EDGetTokenT beamSpotToken_; + edm::EDGetTokenT PrimaryVertexToken_; + edm::ESGetToken ttbToken_; + SingleTrackVertexConstraint stvc; +}; + +DEFINE_FWK_MODULE(MuonBeamspotConstraintValueMapProducer); diff --git a/PhysicsTools/NanoAOD/python/muons_cff.py b/PhysicsTools/NanoAOD/python/muons_cff.py index 473fe17b88d3c..9ea05412584ce 100644 --- a/PhysicsTools/NanoAOD/python/muons_cff.py +++ b/PhysicsTools/NanoAOD/python/muons_cff.py @@ -119,6 +119,11 @@ weightFile = "PhysicsTools/NanoAOD/data/mu_BDTG_2016.weights.xml", ) +from TrackingTools.TransientTrack.TransientTrackBuilder_cfi import * +muonBSConstrain = cms.EDProducer("MuonBeamspotConstraintValueMapProducer", + src = cms.InputTag("linkedObjects","muons"), +) + muonTable = simpleCandidateFlatTableProducer.clone( src = cms.InputTag("linkedObjects","muons"), name = cms.string("Muon"), @@ -175,9 +180,15 @@ mvaTTH = ExtVar(cms.InputTag("muonMVATTH"),float, doc="TTH MVA lepton ID score",precision=14), mvaLowPt = ExtVar(cms.InputTag("muonMVALowPt"),float, doc="Low pt muon ID score",precision=14), fsrPhotonIdx = ExtVar(cms.InputTag("leptonFSRphotons:muFsrIndex"), "int16", doc="Index of the lowest-dR/ET2 among associated FSR photons"), + bsConstrainedPt = ExtVar(cms.InputTag("muonBSConstrain:muonBSConstrainedPt"),float, doc="pT with beamspot constraint",precision=-1), + bsConstrainedPtErr = ExtVar(cms.InputTag("muonBSConstrain:muonBSConstrainedPtErr"),float, doc="pT error with beamspot constraint ",precision=6), + bsConstrainedChi2 = ExtVar(cms.InputTag("muonBSConstrain:muonBSConstrainedPtErr"),float, doc="pT error with beamspot constraint ",precision=6), ), ) +# Increase precision of eta and phi +muonTable.variables.eta.precision = 16 +muonTable.variables.phi.precision = 16 (run2_nanoAOD_106Xv2 | run3_nanoAOD_122).toModify(muonTable.variables,mvaMuID=None).toModify( muonTable.variables, mvaMuID = Var("userFloat('mvaIDMuon')", float, doc="MVA-based ID score",precision=6)) @@ -212,8 +223,5 @@ muonTask = cms.Task(slimmedMuonsUpdated,isoForMu,ptRatioRelForMu,slimmedMuonsWithUserData,finalMuons,finalLooseMuons ) muonMCTask = cms.Task(muonsMCMatchForTable,muonMCTable) -muonTablesTask = cms.Task(muonMVATTH,muonMVALowPt,muonTable,muonMVAID) - - - +muonTablesTask = cms.Task(muonMVATTH,muonMVALowPt,muonBSConstrain,muonTable,muonMVAID)