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