From 84637116006c50068010a54dfb208e128cb65285 Mon Sep 17 00:00:00 2001 From: Michal Bluj Date: Tue, 14 Nov 2023 11:26:56 +0100 Subject: [PATCH] Life-time information updater for pat lepton collections --- .../plugins/PATLeptonTimeLifeInfoUpdater.cc | 331 ++++++++++++++++++ 1 file changed, 331 insertions(+) create mode 100644 PhysicsTools/PatAlgos/plugins/PATLeptonTimeLifeInfoUpdater.cc diff --git a/PhysicsTools/PatAlgos/plugins/PATLeptonTimeLifeInfoUpdater.cc b/PhysicsTools/PatAlgos/plugins/PATLeptonTimeLifeInfoUpdater.cc new file mode 100644 index 0000000000000..1bb901c2cdb36 --- /dev/null +++ b/PhysicsTools/PatAlgos/plugins/PATLeptonTimeLifeInfoUpdater.cc @@ -0,0 +1,331 @@ +/** + \class PATLeptonTimeLifeInfoUpdater + \brief Produces updated collection of pat::Electron's or pat::Muon's + + The PATLeptonInfoUpdater produces analysis-level pat::Electron's or + pat::Muon's with updated time-life information. + + \author Michal Bluj, NCBJ, Warsaw +*/ + +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" + +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/PatCandidates/interface/Tau.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/TrackReco/interface/Track.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "RecoVertex/VertexTools/interface/VertexDistance3D.h" +#include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h" +#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h" +#include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h" + +#include + +template +class PATLeptonTimeLifeInfoUpdater : public edm::stream::EDProducer<> { +public: + explicit PATLeptonTimeLifeInfoUpdater(const edm::ParameterSet&); + ~PATLeptonTimeLifeInfoUpdater() override{}; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + void produce(edm::Event&, const edm::EventSetup&) override; + +private: + //--- private utility methods + const reco::Track* getTrack(const T&); + void produceAndFillIPInfo(T&, const TransientTrackBuilder&, const reco::Vertex&); + void produceAndFillSVInfo(T&, const TransientTrackBuilder&, const reco::Vertex&); + + //--- configuration parameters + edm::EDGetTokenT> leptonsToken_; + edm::EDGetTokenT pvToken_; + edm::ESGetToken transTrackBuilderToken_; + int pvChoice_; + + enum PVChoice { useFront = 0, useClosestInDz }; +}; + +template +PATLeptonTimeLifeInfoUpdater::PATLeptonTimeLifeInfoUpdater(const edm::ParameterSet& cfg) + : leptonsToken_(consumes>(cfg.getParameter("src"))), + pvToken_(consumes(cfg.getParameter("pvSource"))), + transTrackBuilderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))), + pvChoice_(cfg.getParameter("pvChoice")) { + produces>(); +} + +template +void PATLeptonTimeLifeInfoUpdater::produce(edm::Event& evt, const edm::EventSetup& es) { + // Get leptons + edm::Handle> inputLeptons; + evt.getByToken(leptonsToken_, inputLeptons); + + auto outputLeptons = std::make_unique>(); + outputLeptons->reserve(inputLeptons->size()); + + // Get the vertices + edm::Handle vertices; + evt.getByToken(pvToken_, vertices); + + // Get transient track builder + const TransientTrackBuilder& transTrackBuilder = es.getData(transTrackBuilderToken_); + + for (const auto& inputLepton : *inputLeptons) { + T outputLepton(inputLepton); + size_t pv_idx = 0; + if (pvChoice_ == useClosestInDz && getTrack(inputLepton) != nullptr) { + float dz_min = 999; + size_t vtx_idx = 0; + for (const auto& vtx : *vertices) { + float dz_tmp = std::abs(getTrack(inputLepton)->dz(vtx.position())); + if (dz_tmp < dz_min) { + dz_min = dz_tmp; + pv_idx = vtx_idx; + } + vtx_idx++; + } + } + const reco::Vertex& pv = !vertices->empty() ? (*vertices)[pv_idx] : reco::Vertex(); + // Obtain IP vector and set related info into lepton + produceAndFillIPInfo(outputLepton, transTrackBuilder, pv); + + // Fit SV and set related info in taus or do nothing for other lepton types + produceAndFillSVInfo(outputLepton, transTrackBuilder, pv); + + // Set updated lepton to output collection + outputLeptons->push_back(outputLepton); + } // end of lepton loop + + evt.put(std::move(outputLeptons)); +} + +template <> +const reco::Track* PATLeptonTimeLifeInfoUpdater::getTrack(const pat::Electron& electron) { + return electron.gsfTrack().isNonnull() ? electron.gsfTrack().get() : nullptr; +} + +template <> +const reco::Track* PATLeptonTimeLifeInfoUpdater::getTrack(const pat::Muon& muon) { + return muon.innerTrack().isNonnull() ? muon.innerTrack().get() : nullptr; +} + +template <> +const reco::Track* PATLeptonTimeLifeInfoUpdater::getTrack(const pat::Tau& tau) { + const reco::Track* track = nullptr; + if (tau.leadChargedHadrCand().isNonnull()) + track = tau.leadChargedHadrCand()->bestTrack(); + return track; +} + +template +void PATLeptonTimeLifeInfoUpdater::produceAndFillIPInfo(T& lepton, + const TransientTrackBuilder& transTrackBuilder, + const reco::Vertex& pv) { + GlobalPoint pca; + GlobalError pca_cov; + GlobalVector ip_vec; + GlobalError ip_cov; + float ipLength = -1; + float ipLength_sig = -1; + float bField_z = 0; + const reco::Track* track = getTrack(lepton); + if (track != nullptr) { + bField_z = transTrackBuilder.field()->inInverseGeV(GlobalPoint(track->vx(), track->vy(), track->vz())).z(); + // Extrapolate track to the point closest to PV + reco::TransientTrack transTrack = transTrackBuilder.build(track); + AnalyticalImpactPointExtrapolator extrapolator(transTrack.field()); + TrajectoryStateOnSurface closestState = + extrapolator.extrapolate(transTrack.impactPointState(), RecoVertex::convertPos(pv.position())); + pca = closestState.globalPosition(); + pca_cov = closestState.cartesianError().position(); + ip_vec = GlobalVector(pca.x() - pv.x(), pca.y() - pv.y(), pca.z() - pv.z()); + ip_cov = pca_cov + GlobalError(pv.covariance()); + VertexDistance3D pca_dist; + Measurement1D ip_mes = pca_dist.distance(pv, VertexState(pca, pca_cov)); + if (ip_vec.dot(GlobalVector(lepton.px(), lepton.py(), lepton.pz())) < 0) + ip_mes = Measurement1D(-1. * ip_mes.value(), ip_mes.error()); + + ipLength = ip_mes.value(); + ipLength_sig = ip_mes.significance(); + } + // Store PCA info + lepton.addUserFloat("refitPCA_x", pca.x()); + lepton.addUserFloat("refitPCA_y", pca.y()); + lepton.addUserFloat("refitPCA_z", pca.z()); + lepton.addUserFloat("refitPCACov_xx", pca_cov.cxx()); + lepton.addUserFloat("refitPCACov_yx", pca_cov.cyx()); + lepton.addUserFloat("refitPCACov_zx", pca_cov.czx()); + lepton.addUserFloat("refitPCACov_yy", pca_cov.cyy()); + lepton.addUserFloat("refitPCACov_zy", pca_cov.czy()); + lepton.addUserFloat("refitPCACov_zz", pca_cov.czz()); + // Store IP info + lepton.addUserFloat("IP", ipLength); + lepton.addUserFloat("IP_sig", ipLength_sig); + lepton.addUserFloat("IP_x", ip_vec.x()); + lepton.addUserFloat("IP_y", ip_vec.y()); + lepton.addUserFloat("IP_z", ip_vec.z()); + lepton.addUserFloat("IPCov_xx", ip_cov.cxx()); + lepton.addUserFloat("IPCov_yx", ip_cov.cyx()); + lepton.addUserFloat("IPCov_zx", ip_cov.czx()); + lepton.addUserFloat("IPCov_yy", ip_cov.cyy()); + lepton.addUserFloat("IPCov_zy", ip_cov.czy()); + lepton.addUserFloat("IPCov_zz", ip_cov.czz()); + // Store leading track parameters and covariance at reference point + if (track != nullptr) { + lepton.addUserFloat("track_qoverp", track->parameter(reco::TrackBase::i_qoverp)); + lepton.addUserFloat("track_lambda", track->parameter(reco::TrackBase::i_lambda)); + lepton.addUserFloat("track_phi", track->parameter(reco::TrackBase::i_phi)); + lepton.addUserFloat("track_dxy", track->parameter(reco::TrackBase::i_dxy)); + lepton.addUserFloat("track_dsz", track->parameter(reco::TrackBase::i_dsz)); + for (size_t i = 0; i < reco::TrackBase::dimension; ++i) + for (size_t j = i; j < reco::TrackBase::dimension; ++j) + lepton.addUserFloat("trackCov_" + std::to_string(j + i), track->covariance(j, i)); + } else { + lepton.addUserFloat("track_qoverp", 0); + lepton.addUserFloat("track_lambda", 0); + lepton.addUserFloat("track_phi", 0); + lepton.addUserFloat("track_dxy", 0); + lepton.addUserFloat("track_dsz", 0); + for (size_t i = 0; i < reco::TrackBase::dimension; ++i) + for (size_t j = i; j < reco::TrackBase::dimension; ++j) + lepton.addUserFloat("trackCov_" + std::to_string(j + i), 0); + } + // ... and z coordinate of magnetic field at track reference point + lepton.addUserFloat("bField_z", bField_z); +} + +template +void PATLeptonTimeLifeInfoUpdater::produceAndFillSVInfo(T& lepton, + const TransientTrackBuilder& transTrackBuilder, + const reco::Vertex& pv) {} + +template <> +void PATLeptonTimeLifeInfoUpdater::produceAndFillSVInfo(pat::Tau& tau, + const TransientTrackBuilder& transTrackBuilder, + const reco::Vertex& pv) { + //MB: Set initially SV to the tau vertex (tau.vertex()) or the PV choosen? + reco::Vertex sv; + GlobalVector flight_vec; + GlobalError flight_cov; + float flightLength = -1; + float flightLength_sig = -1; + + // Fit SV with tracks of charged tau decay products + int fitOK = 0; + if (tau.signalChargedHadrCands().size() + tau.signalLostTracks().size() > 1) { + // Get tracks form tau signal charged candidates + std::vector transTrks; + TransientVertex transVtx; + for (const auto& cand : tau.signalChargedHadrCands()) { + if (cand.isNull()) + continue; + const reco::Track* track = cand->bestTrack(); + if (track != nullptr) + transTrks.push_back(transTrackBuilder.build(track)); + } + for (const auto& cand : tau.signalLostTracks()) { + if (cand.isNull()) + continue; + const reco::Track* track = cand->bestTrack(); + if (track != nullptr) + transTrks.push_back(transTrackBuilder.build(track)); + } + // Fit SV with KalmanVertexFitter + fitOK = 1; + KalmanVertexFitter kvf(true); + if (transTrks.size() > 1) { + try { + transVtx = kvf.vertex(transTrks); + } catch (...) { + fitOK = -1; + } + } else { + fitOK = -1; + } + if (!transVtx.hasRefittedTracks() || transVtx.refittedTracks().size() != transTrks.size()) + fitOK = -1; + if (fitOK > 0) { + sv = transVtx; + // Get flight-length + // Full PV->SV flight vector with its covariance + flight_vec = GlobalVector(sv.x() - pv.x(), sv.y() - pv.y(), sv.z() - pv.z()); + flight_cov = transVtx.positionError() + GlobalError(pv.covariance()); + //MB: can be taken from tau itself (but with different fit of PV and SV) as follows: + //flightLength = sqrt(tau.flightLength().mag2()); + //flightLength_sig = tau.flightLengthSig(); + VertexDistance3D sv_dist; + Measurement1D flightLength_mes = sv_dist.signedDistance(pv, sv, GlobalVector(tau.px(), tau.py(), tau.pz())); + flightLength = flightLength_mes.value(); + flightLength_sig = flightLength_mes.significance(); + } + } + // Embed SV info in tau + // Store SV info + tau.addUserInt("refitSV_status", fitOK); + tau.addUserFloat("refitSV_x", sv.x()); + tau.addUserFloat("refitSV_y", sv.y()); + tau.addUserFloat("refitSV_z", sv.z()); + tau.addUserFloat("refitSV_chi2", sv.chi2()); + tau.addUserFloat("refitSV_ndof", sv.ndof()); + tau.addUserFloat("refitSVCov_xx", sv.covariance(0, 0)); + tau.addUserFloat("refitSVCov_yx", sv.covariance(1, 0)); + tau.addUserFloat("refitSVCov_zx", sv.covariance(2, 0)); + tau.addUserFloat("refitSVCov_yy", sv.covariance(1, 1)); + tau.addUserFloat("refitSVCov_zy", sv.covariance(2, 1)); + tau.addUserFloat("refitSVCov_zz", sv.covariance(2, 2)); + // Store flight-length info + tau.addUserFloat("flightLength", flightLength); + tau.addUserFloat("flightLength_sig", flightLength_sig); + tau.addUserFloat("flight_x", flight_vec.x()); + tau.addUserFloat("flight_y", flight_vec.y()); + tau.addUserFloat("flight_z", flight_vec.z()); + tau.addUserFloat("flightCov_xx", flight_cov.cxx()); + tau.addUserFloat("flightCov_yx", flight_cov.cyx()); + tau.addUserFloat("flightCov_zx", flight_cov.czx()); + tau.addUserFloat("flightCov_yy", flight_cov.cyy()); + tau.addUserFloat("flightCov_zy", flight_cov.czy()); + tau.addUserFloat("flightCov_zz", flight_cov.czz()); +} + +template +void PATLeptonTimeLifeInfoUpdater::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + // pat{Electron,Muon,Tau}TimeLifeInfoUpdater + edm::ParameterSetDescription desc; + + std::string lepCollName; + if (typeid(T) == typeid(pat::Electron)) + lepCollName = "slimmedElectrons"; + else if (typeid(T) == typeid(pat::Muon)) + lepCollName = "slimmedMuons"; + else if (typeid(T) == typeid(pat::Tau)) + lepCollName = "slimmedTaus"; + desc.add("src", edm::InputTag(lepCollName)); + desc.add("pvSource", edm::InputTag("offlineSlimmedPrimaryVertices")); + desc.add("pvChoice", useFront) + ->setComment( + "Define PV to compute IP: 0: first PV, 1: PV with the smallest dz of the tau leading track (default: 0)"); + + descriptions.addWithDefaultLabel(desc); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +typedef PATLeptonTimeLifeInfoUpdater PATElectronTimeLifeInfoUpdater; +DEFINE_FWK_MODULE(PATElectronTimeLifeInfoUpdater); +typedef PATLeptonTimeLifeInfoUpdater PATMuonTimeLifeInfoUpdater; +DEFINE_FWK_MODULE(PATMuonTimeLifeInfoUpdater); +typedef PATLeptonTimeLifeInfoUpdater PATTauTimeLifeInfoUpdater; +DEFINE_FWK_MODULE(PATTauTimeLifeInfoUpdater);