Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:23:50

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