GlobalTrackingTools/plugins/MuonBeamspotConstraintValueMapProducer.cc

0001 /** \class MuonBeamspotConstraintValueMapProducer
0002  * Compute muon pt and ptErr after beamspot constraint.
0003  *
0004  */
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "DataFormats/PatCandidates/interface/Muon.h"
0009 #include "DataFormats/Common/interface/ValueMap.h"
0010 #include "FWCore/Framework/interface/global/EDProducer.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0013 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0014 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0015 #include "RecoVertex/KalmanVertexFit/interface/SingleTrackVertexConstraint.h"
0016 #include "DataFormats/VertexReco/interface/Vertex.h"
0017
0018 class MuonBeamspotConstraintValueMapProducer : public edm::global::EDProducer<> {
0019 public:
0020   explicit MuonBeamspotConstraintValueMapProducer(const edm::ParameterSet& config)
0021       : muonToken_(consumes<pat::MuonCollection>(config.getParameter<edm::InputTag>("src"))),
0022         beamSpotToken_(consumes<reco::BeamSpot>(config.getParameter<edm::InputTag>("beamspot"))),
0023         PrimaryVertexToken_(consumes<reco::VertexCollection>(config.getParameter<edm::InputTag>("vertices"))),
0024         ttbToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) {
0025     produces<edm::ValueMap<float>>("muonBSConstrainedPt");
0026     produces<edm::ValueMap<float>>("muonBSConstrainedPtErr");
0027     produces<edm::ValueMap<float>>("muonBSConstrainedChi2");
0028   }
0029
0030   ~MuonBeamspotConstraintValueMapProducer() override = default;
0031
0032   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0033     edm::ParameterSetDescription desc;
0034     desc.add<edm::InputTag>("src", edm::InputTag("muons"))->setComment("Muon collection");
0035     desc.add<edm::InputTag>("beamspot", edm::InputTag("offlineBeamSpot"))->setComment("Beam spot collection");
0036     desc.add<edm::InputTag>("vertices", edm::InputTag("offlineSlimmedPrimaryVertices"))
0037         ->setComment("Primary vertex collection");
0038
0039     descriptions.addWithDefaultLabel(desc);
0040   }
0041
0042 private:
0043   void produce(edm::StreamID streamID, edm::Event& event, const edm::EventSetup& setup) const override {
0044     edm::Handle<pat::MuonCollection> muons;
0045     event.getByToken(muonToken_, muons);
0046
0047     edm::Handle<reco::BeamSpot> beamSpotHandle;
0048     event.getByToken(beamSpotToken_, beamSpotHandle);
0049
0050     edm::ESHandle<TransientTrackBuilder> ttkb = setup.getHandle(ttbToken_);
0051
0052     std::vector<float> pts, ptErrs, chi2s;
0053     pts.reserve(muons->size());
0054     ptErrs.reserve(muons->size());
0055     chi2s.reserve(muons->size());
0056
0057     for (const auto& muon : *muons) {
0058       bool tbd = true;
0059       if (beamSpotHandle.isValid()) {
0060         double BeamWidthX = beamSpotHandle->BeamWidthX();
0061         double BeamWidthXError = beamSpotHandle->BeamWidthXError();
0062         double BeamWidthY = beamSpotHandle->BeamWidthY();
0063         double BeamWidthYError = beamSpotHandle->BeamWidthYError();
0064         // Protect for mis-reconstructed beamspots (note that
0065         // SingleTrackVertexConstraint uses the width for the constraint,
0066         // not the error)
0067         if ((BeamWidthXError / BeamWidthX < 0.3) && (BeamWidthYError / BeamWidthY < 0.3)) {
0068           SingleTrackVertexConstraint::BTFtuple btft =
0069               stvc.constrain(ttkb->build(muon.muonBestTrack()), *beamSpotHandle);
0070           if (std::get<0>(btft)) {
0071             const reco::Track& trkBS = std::get<1>(btft).track();
0072             pts.push_back(trkBS.pt());
0073             ptErrs.push_back(trkBS.ptError());
0074             chi2s.push_back(std::get<2>(btft));
0075             tbd = false;
0076           }
0077         }
0078       }
0079
0080       if (tbd) {
0081         // Invalid BS; use PV instead
0082         edm::Handle<reco::VertexCollection> pvHandle;
0083         event.getByToken(PrimaryVertexToken_, pvHandle);
0084
0085         if (pvHandle.isValid() && !pvHandle->empty()) {
0086           auto pv = pvHandle->at(0);
0087           VertexState pvs = VertexState(GlobalPoint(Basic3DVector<float>(pv.position())), GlobalError(pv.covariance()));
0088
0089           SingleTrackVertexConstraint::BTFtuple btft = stvc.constrain(ttkb->build(muon.muonBestTrack()), pvs);
0090           if (std::get<0>(btft)) {
0091             const reco::Track& trkBS = std::get<1>(btft).track();
0092             pts.push_back(trkBS.pt());
0093             ptErrs.push_back(trkBS.ptError());
0094             chi2s.push_back(std::get<2>(btft));
0095             tbd = false;
0096           }
0097         }
0098       }
0099
0100       if (tbd) {
0101         // Fall-back case, keep the unconstrained values
0102         pts.push_back(muon.pt());
0103         ptErrs.push_back(muon.bestTrack()->ptError());
0104         chi2s.push_back(-1.f);
0105       }
0106     }
0107
0108     {
0109       std::unique_ptr<edm::ValueMap<float>> valueMap(new edm::ValueMap<float>());
0110       edm::ValueMap<float>::Filler filler(*valueMap);
0111       filler.insert(muons, pts.begin(), pts.end());
0112       filler.fill();
0113       event.put(std::move(valueMap), "muonBSConstrainedPt");
0114     }
0115
0116     {
0117       std::unique_ptr<edm::ValueMap<float>> valueMap(new edm::ValueMap<float>());
0118       edm::ValueMap<float>::Filler filler(*valueMap);
0119       filler.insert(muons, ptErrs.begin(), ptErrs.end());
0120       filler.fill();
0121       event.put(std::move(valueMap), "muonBSConstrainedPtErr");
0122     }
0123
0124     {
0125       std::unique_ptr<edm::ValueMap<float>> valueMap(new edm::ValueMap<float>());
0126       edm::ValueMap<float>::Filler filler(*valueMap);
0127       filler.insert(muons, chi2s.begin(), chi2s.end());
0128       filler.fill();
0129       event.put(std::move(valueMap), "muonBSConstrainedChi2");
0130     }
0131   }
0132
0133   edm::EDGetTokenT<pat::MuonCollection> muonToken_;
0134   edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0135   edm::EDGetTokenT<reco::VertexCollection> PrimaryVertexToken_;
0136   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> ttbToken_;
0137   SingleTrackVertexConstraint stvc;
0138 };
0139
0140 DEFINE_FWK_MODULE(MuonBeamspotConstraintValueMapProducer);