PatAlgos/plugins/PATLeptonTimeLifeInfoProducer.cc

0001 /**
0002   \class    PATLeptonTimeLifeInfoProducer
0003   \brief    Produces lepton life-time information
0004
0005   \author   Michal Bluj, NCBJ, Warsaw
0006 */
0007
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0014 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "DataFormats/Common/interface/ValueMap.h"
0018 #include "DataFormats/PatCandidates/interface/Electron.h"
0019 #include "DataFormats/PatCandidates/interface/Muon.h"
0020 #include "DataFormats/PatCandidates/interface/Tau.h"
0021 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0022 #include "DataFormats/TrackReco/interface/Track.h"
0023 #include "DataFormats/VertexReco/interface/Vertex.h"
0024 #include "DataFormats/VertexReco/interface/TrackTimeLifeInfo.h"
0025
0026 #include "MagneticField/Engine/interface/MagneticField.h"
0027 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0028 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0029 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0030 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0031 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0032 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0033 #include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h"
0034 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0035 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0036
0037 #include <cstring>
0038
0039 template <typename T>
0040 class PATLeptonTimeLifeInfoProducer : public edm::stream::EDProducer<> {
0041 public:
0042   explicit PATLeptonTimeLifeInfoProducer(const edm::ParameterSet&);
0043   ~PATLeptonTimeLifeInfoProducer() override {}
0044
0045   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0046   void produce(edm::Event&, const edm::EventSetup&) override;
0047
0048 private:
0049   //--- private utility methods
0050   const reco::Track* getTrack(const T&);
0051   void produceAndFillIPInfo(const T&, const TransientTrackBuilder&, const reco::Vertex&, TrackTimeLifeInfo&);
0052   void produceAndFillSVInfo(const T&, const TransientTrackBuilder&, const reco::Vertex&, TrackTimeLifeInfo&);
0053   static bool fitVertex(const std::vector<reco::TransientTrack>& transTrk, TransientVertex& transVtx) {
0054     if (transTrk.size() < 2)
0055       return false;
0056     try {
0057       KalmanVertexFitter kvf(true);
0058       transVtx = kvf.vertex(transTrk);
0059       return transVtx.hasRefittedTracks() && transVtx.refittedTracks().size() == transTrk.size();
0060     } catch (VertexException& e) {
0061       edm::LogWarning("PATLeptonTimeLifeInfoProducer") << " fitVertex failed: " << e.what();
0062       return false;
0063     }
0064   }
0065
0066   //--- configuration parameters
0067   edm::EDGetTokenT<std::vector<T>> leptonsToken_;
0068   edm::EDGetTokenT<reco::VertexCollection> pvToken_;
0069   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transTrackBuilderToken_;
0070   const StringCutObjectSelector<T> selector_;
0071   int pvChoice_;
0072
0073   enum PVChoice { useFront = 0, useClosestInDz };
0074   //--- value map for TrackTimeLifeInfo (to be stored into the event)
0075   using TrackTimeLifeInfoMap = edm::ValueMap<TrackTimeLifeInfo>;
0076 };
0077
0078 template <typename T>
0079 PATLeptonTimeLifeInfoProducer<T>::PATLeptonTimeLifeInfoProducer(const edm::ParameterSet& cfg)
0080     : leptonsToken_(consumes<std::vector<T>>(cfg.getParameter<edm::InputTag>("src"))),
0081       pvToken_(consumes<reco::VertexCollection>(cfg.getParameter<edm::InputTag>("pvSource"))),
0082       transTrackBuilderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0083       selector_(cfg.getParameter<std::string>("selection")),
0084       pvChoice_(cfg.getParameter<int>("pvChoice")) {
0085   produces<TrackTimeLifeInfoMap>();
0086 }
0087
0088 template <typename T>
0089 void PATLeptonTimeLifeInfoProducer<T>::produce(edm::Event& evt, const edm::EventSetup& es) {
0090   // Get leptons
0091   edm::Handle<std::vector<T>> leptons;
0092   evt.getByToken(leptonsToken_, leptons);
0093
0094   // Get the vertices
0095   edm::Handle<reco::VertexCollection> vertices;
0096   evt.getByToken(pvToken_, vertices);
0097
0098   // Get transient track builder
0099   const TransientTrackBuilder& transTrackBuilder = es.getData(transTrackBuilderToken_);
0100
0101   std::vector<TrackTimeLifeInfo> infos;
0102   infos.reserve(leptons->size());
0103
0104   for (const auto& lepton : *leptons) {
0105     TrackTimeLifeInfo info;
0106
0107     // Do nothing for lepton not passing selection
0108     if (!selector_(lepton)) {
0109       infos.push_back(info);
0110       continue;
0111     }
0112     size_t pv_idx = 0;
0113     if (pvChoice_ == useClosestInDz && getTrack(lepton) != nullptr) {
0114       float dz_min = 999;
0115       size_t vtx_idx = 0;
0116       for (const auto& vtx : *vertices) {
0117         float dz_tmp = std::abs(getTrack(lepton)->dz(vtx.position()));
0118         if (dz_tmp < dz_min) {
0119           dz_min = dz_tmp;
0120           pv_idx = vtx_idx;
0121         }
0122         vtx_idx++;
0123       }
0124     }
0125     const reco::Vertex& pv = !vertices->empty() ? (*vertices)[pv_idx] : reco::Vertex();
0126
0127     // Obtain IP vector and set related info into lepton
0128     produceAndFillIPInfo(lepton, transTrackBuilder, pv, info);
0129
0130     // Fit SV and set related info for taus or do nothing for other lepton types
0131     produceAndFillSVInfo(lepton, transTrackBuilder, pv, info);
0132     infos.push_back(info);
0133   }  // end of lepton loop
0134
0135   // Build the valuemap
0136   auto infoMap = std::make_unique<TrackTimeLifeInfoMap>();
0137   TrackTimeLifeInfoMap::Filler filler(*infoMap);
0138   filler.insert(leptons, infos.begin(), infos.end());
0139   filler.fill();
0140
0141   // Store output into the event
0142   evt.put(std::move(infoMap));
0143 }
0144
0145 template <>
0146 const reco::Track* PATLeptonTimeLifeInfoProducer<pat::Electron>::getTrack(const pat::Electron& electron) {
0147   return electron.gsfTrack().isNonnull() ? electron.gsfTrack().get() : nullptr;
0148 }
0149
0150 template <>
0151 const reco::Track* PATLeptonTimeLifeInfoProducer<pat::Muon>::getTrack(const pat::Muon& muon) {
0152   return muon.innerTrack().isNonnull() ? muon.innerTrack().get() : nullptr;
0153 }
0154
0155 template <>
0156 const reco::Track* PATLeptonTimeLifeInfoProducer<pat::Tau>::getTrack(const pat::Tau& tau) {
0157   const reco::Track* track = nullptr;
0158   if (tau.leadChargedHadrCand().isNonnull())
0159     track = tau.leadChargedHadrCand()->bestTrack();
0160   return track;
0161 }
0162
0163 template <typename T>
0164 void PATLeptonTimeLifeInfoProducer<T>::produceAndFillIPInfo(const T& lepton,
0165                                                             const TransientTrackBuilder& transTrackBuilder,
0166                                                             const reco::Vertex& pv,
0167                                                             TrackTimeLifeInfo& info) {
0168   const reco::Track* track = getTrack(lepton);
0169   if (track != nullptr) {
0170     // Extrapolate track to the point closest to PV
0171     reco::TransientTrack transTrack = transTrackBuilder.build(track);
0172     AnalyticalImpactPointExtrapolator extrapolator(transTrack.field());
0173     TrajectoryStateOnSurface closestState =
0174         extrapolator.extrapolate(transTrack.impactPointState(), RecoVertex::convertPos(pv.position()));
0175     if (!closestState.isValid()) {
0176       edm::LogError("PATLeptonTimeLifeInfoProducer")
0177           << "closestState not valid! From:\n"
0178           << "transTrack.impactPointState():\n"
0179           << transTrack.impactPointState() << "RecoVertex::convertPos(pv.position()):\n"
0180           << RecoVertex::convertPos(pv.position());
0181       return;
0182     }
0183     GlobalPoint pca = closestState.globalPosition();
0184     GlobalError pca_cov = closestState.cartesianError().position();
0185     GlobalVector ip_vec = GlobalVector(pca.x() - pv.x(), pca.y() - pv.y(), pca.z() - pv.z());
0186     GlobalError ip_cov = pca_cov + GlobalError(pv.covariance());
0187     VertexDistance3D pca_dist;
0188     Measurement1D ip_mes = pca_dist.distance(pv, VertexState(pca, pca_cov));
0189     if (ip_vec.dot(GlobalVector(lepton.px(), lepton.py(), lepton.pz())) < 0)
0190       ip_mes = Measurement1D(-1. * ip_mes.value(), ip_mes.error());
0191
0192     // Store Track and PCA info
0193     info.setTrack(track);
0194     info.setBField_z(transTrackBuilder.field()->inInverseGeV(GlobalPoint(track->vx(), track->vy(), track->vz())).z());
0195     info.setPCA(pca, pca_cov);
0196     info.setIP(ip_vec, ip_cov);
0197     info.setIPLength(ip_mes);
0198   }
0199 }
0200
0201 template <typename T>
0202 void PATLeptonTimeLifeInfoProducer<T>::produceAndFillSVInfo(const T& lepton,
0203                                                             const TransientTrackBuilder& transTrackBuilder,
0204                                                             const reco::Vertex& pv,
0205                                                             TrackTimeLifeInfo& info) {}
0206
0207 template <>
0208 void PATLeptonTimeLifeInfoProducer<pat::Tau>::produceAndFillSVInfo(const pat::Tau& tau,
0209                                                                    const TransientTrackBuilder& transTrackBuilder,
0210                                                                    const reco::Vertex& pv,
0211                                                                    TrackTimeLifeInfo& info) {
0212   // Fit SV with tracks of charged tau decay products
0213   int fitOK = 0;
0214   if (tau.signalChargedHadrCands().size() + tau.signalLostTracks().size() > 1) {
0215     // Get tracks from tau signal charged candidates
0216     std::vector<reco::TransientTrack> transTrks;
0217     TransientVertex transVtx;
0218     for (const auto& cand : tau.signalChargedHadrCands()) {
0219       if (cand.isNull())
0220         continue;
0221       const reco::Track* track = cand->bestTrack();
0222       if (track != nullptr)
0223         transTrks.push_back(transTrackBuilder.build(track));
0224     }
0225     for (const auto& cand : tau.signalLostTracks()) {
0226       if (cand.isNull())
0227         continue;
0228       const reco::Track* track = cand->bestTrack();
0229       if (track != nullptr)
0230         transTrks.push_back(transTrackBuilder.build(track));
0231     }
0232     // Fit SV with KalmanVertexFitter
0233     fitOK = fitVertex(transTrks, transVtx) ? 1 : -1;
0234     if (fitOK > 0) {
0235       reco::Vertex sv = transVtx;
0236       // Get flight-length
0237       // Full PV->SV flight vector with its covariance
0238       GlobalVector flight_vec = GlobalVector(sv.x() - pv.x(), sv.y() - pv.y(), sv.z() - pv.z());
0239       GlobalError flight_cov = transVtx.positionError() + GlobalError(pv.covariance());
0240       //MB: can be taken from tau itself (but with different fit of PV and SV) as follows:
0241       //tau.flightLength().mag2());
0242       //tau.flightLengthSig();
0243       VertexDistance3D sv_dist;
0244       Measurement1D flightLength_mes = sv_dist.signedDistance(pv, sv, GlobalVector(tau.px(), tau.py(), tau.pz()));
0245
0246       // Store SV info
0247       info.setSV(sv);
0248       info.setFlightVector(flight_vec, flight_cov);
0249       info.setFlightLength(flightLength_mes);
0250     }
0251   }
0252 }
0253
0254 template <typename T>
0255 void PATLeptonTimeLifeInfoProducer<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0256   // pat{Electron,Muon,Tau}TimeLifeInfoProducer
0257   edm::ParameterSetDescription desc;
0258
0259   std::string lepCollName;
0260   if (typeid(T) == typeid(pat::Electron))
0261     lepCollName = "slimmedElectrons";
0262   else if (typeid(T) == typeid(pat::Muon))
0263     lepCollName = "slimmedMuons";
0264   else if (typeid(T) == typeid(pat::Tau))
0265     lepCollName = "slimmedTaus";
0266   desc.add<edm::InputTag>("src", edm::InputTag(lepCollName));
0267   desc.add<edm::InputTag>("pvSource", edm::InputTag("offlineSlimmedPrimaryVertices"));
0268   desc.add<std::string>("selection", "")->setComment("Selection required to produce and store time-life information");
0269   desc.add<int>("pvChoice", useFront)
0270       ->setComment(
0271           "Define PV to compute IP: 0: first PV, 1: PV with the smallest dz of the tau leading track (default: " +
0272           std::to_string(useFront) + ")");
0273
0274   descriptions.addWithDefaultLabel(desc);
0275 }
0276
0277 #include "FWCore/Framework/interface/MakerMacros.h"
0278 typedef PATLeptonTimeLifeInfoProducer<pat::Electron> PATElectronTimeLifeInfoProducer;
0279 DEFINE_FWK_MODULE(PATElectronTimeLifeInfoProducer);
0280 typedef PATLeptonTimeLifeInfoProducer<pat::Muon> PATMuonTimeLifeInfoProducer;
0281 DEFINE_FWK_MODULE(PATMuonTimeLifeInfoProducer);
0282 typedef PATLeptonTimeLifeInfoProducer<pat::Tau> PATTauTimeLifeInfoProducer;
0283 DEFINE_FWK_MODULE(PATTauTimeLifeInfoProducer);