RecoTau/plugins/PFTauTransverseImpactParameters.cc

0001 /* class PFTauTransverseImpactParamters
0002  * EDProducer of the
0003  * authors: Ian M. Nugent
0004  * This work is based on the impact parameter work by Rosamaria Venditti and reconstructing the 3 prong taus.
0005  * The idea of the fully reconstructing the tau using a kinematic fit comes from
0006  * Lars Perchalla and Philip Sauerland Theses under Achim Stahl supervision. This
0007  * work was continued by Ian M. Nugent and Vladimir Cherepanov.
0008  * Thanks goes to Christian Veelken and Evan Klose Friis for their help and suggestions.
0009  */
0010
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/stream/EDProducer.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Utilities/interface/Exception.h"
0019
0020 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0021 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0022
0023 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0024 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0025 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0026 #include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
0027 #include "TrackingTools/IPTools/interface/IPTools.h"
0028 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0029 #include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h"
0030 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0031
0032 #include "DataFormats/TauReco/interface/PFTau.h"
0033 #include "DataFormats/TauReco/interface/PFTauFwd.h"
0034 #include "DataFormats/VertexReco/interface/Vertex.h"
0035 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0036 #include "DataFormats/TrackReco/interface/Track.h"
0037 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0038 #include "DataFormats/TauReco/interface/PFTauTransverseImpactParameter.h"
0039 #include "DataFormats/TauReco/interface/PFTauTransverseImpactParameterFwd.h"
0040 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0041
0042 #include "DataFormats/Common/interface/Association.h"
0043 #include "DataFormats/Common/interface/AssociationVector.h"
0044 #include "DataFormats/Common/interface/RefProd.h"
0045 #include "TMath.h"
0046
0047 #include <memory>
0048
0049 using namespace reco;
0050 using namespace edm;
0051 using namespace std;
0052
0053 class PFTauTransverseImpactParameters : public edm::stream::EDProducer<> {
0054 public:
0055   enum CMSSWPerigee { aCurv = 0, aTheta, aPhi, aTip, aLip };
0056   explicit PFTauTransverseImpactParameters(const edm::ParameterSet& iConfig);
0057   ~PFTauTransverseImpactParameters() override;
0058   void produce(edm::Event&, const edm::EventSetup&) override;
0059   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060
0061 private:
0062   edm::EDGetTokenT<std::vector<reco::PFTau>> PFTauToken_;
0063   edm::EDGetTokenT<edm::AssociationVector<PFTauRefProd, std::vector<reco::VertexRef>>> PFTauPVAToken_;
0064   edm::EDGetTokenT<edm::AssociationVector<PFTauRefProd, std::vector<std::vector<reco::VertexRef>>>> PFTauSVAToken_;
0065   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transTrackBuilderToken_;
0066   bool useFullCalculation_;
0067 };
0068
0069 PFTauTransverseImpactParameters::PFTauTransverseImpactParameters(const edm::ParameterSet& iConfig)
0070     : PFTauToken_(consumes<std::vector<reco::PFTau>>(iConfig.getParameter<edm::InputTag>("PFTauTag"))),
0071       PFTauPVAToken_(consumes<edm::AssociationVector<PFTauRefProd, std::vector<reco::VertexRef>>>(
0072           iConfig.getParameter<edm::InputTag>("PFTauPVATag"))),
0073       PFTauSVAToken_(consumes<edm::AssociationVector<PFTauRefProd, std::vector<std::vector<reco::VertexRef>>>>(
0074           iConfig.getParameter<edm::InputTag>("PFTauSVATag"))),
0075       transTrackBuilderToken_(esConsumes(edm::ESInputTag{"", "TransientTrackBuilder"})),
0076       useFullCalculation_(iConfig.getParameter<bool>("useFullCalculation")) {
0077   produces<edm::AssociationVector<PFTauRefProd, std::vector<reco::PFTauTransverseImpactParameterRef>>>();
0078   produces<PFTauTransverseImpactParameterCollection>("PFTauTIP");
0079 }
0080
0081 PFTauTransverseImpactParameters::~PFTauTransverseImpactParameters() {}
0082
0083 namespace {
0084   inline const reco::Track* getTrack(const Candidate& cand) {
0085     const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0086     if (pfCandPtr != nullptr) {
0087       if (pfCandPtr->trackRef().isNonnull())
0088         return pfCandPtr->trackRef().get();
0089       else if (pfCandPtr->gsfTrackRef().isNonnull())
0090         return pfCandPtr->gsfTrackRef().get();
0091       else
0092         return nullptr;
0093     }
0094     const pat::PackedCandidate* packedCand = dynamic_cast<const pat::PackedCandidate*>(&cand);
0095     if (packedCand != nullptr && packedCand->hasTrackDetails())
0096       return &packedCand->pseudoTrack();
0097
0098     return nullptr;
0099   }
0100 }  // namespace
0101
0102 void PFTauTransverseImpactParameters::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0103   // Obtain Collections
0104   auto const& transTrackBuilder = iSetup.getData(transTrackBuilderToken_);
0105
0106   edm::Handle<std::vector<reco::PFTau>> Tau;
0107   iEvent.getByToken(PFTauToken_, Tau);
0108
0109   edm::Handle<edm::AssociationVector<PFTauRefProd, std::vector<reco::VertexRef>>> PFTauPVA;
0110   iEvent.getByToken(PFTauPVAToken_, PFTauPVA);
0111
0112   edm::Handle<edm::AssociationVector<PFTauRefProd, std::vector<std::vector<reco::VertexRef>>>> PFTauSVA;
0113   iEvent.getByToken(PFTauSVAToken_, PFTauSVA);
0114
0115   // Set Association Map
0116   auto AVPFTauTIP =
0117       std::make_unique<edm::AssociationVector<PFTauRefProd, std::vector<reco::PFTauTransverseImpactParameterRef>>>(
0118           PFTauRefProd(Tau));
0119   auto TIPCollection_out = std::make_unique<PFTauTransverseImpactParameterCollection>();
0120   reco::PFTauTransverseImpactParameterRefProd TIPRefProd_out =
0121       iEvent.getRefBeforePut<reco::PFTauTransverseImpactParameterCollection>("PFTauTIP");
0122
0123   // For each Tau Run Algorithim
0124   if (Tau.isValid()) {
0125     for (reco::PFTauCollection::size_type iPFTau = 0; iPFTau < Tau->size(); iPFTau++) {
0126       reco::PFTauRef RefPFTau(Tau, iPFTau);
0127       const reco::VertexRef PV = PFTauPVA->value(RefPFTau.key());
0128       const std::vector<reco::VertexRef> SV = PFTauSVA->value(RefPFTau.key());
0129       double dxy(-999), dxy_err(-999);
0130       reco::Vertex::Point poca(0, 0, 0);
0131       double ip3d(-999), ip3d_err(-999);
0132       reco::Vertex::Point ip3d_poca(0, 0, 0);
0133       if (RefPFTau->leadChargedHadrCand().isNonnull()) {
0134         const reco::Track* track = getTrack(*RefPFTau->leadChargedHadrCand());
0135         if (track != nullptr) {
0136           if (useFullCalculation_) {
0137             reco::TransientTrack transTrk = transTrackBuilder.build(*track);
0138             GlobalVector direction(
0139                 RefPFTau->p4().px(), RefPFTau->p4().py(), RefPFTau->p4().pz());  //To compute sign of IP
0140             std::pair<bool, Measurement1D> signed_IP2D =
0141                 IPTools::signedTransverseImpactParameter(transTrk, direction, (*PV));
0142             dxy = signed_IP2D.second.value();
0143             dxy_err = signed_IP2D.second.error();
0144             std::pair<bool, Measurement1D> signed_IP3D = IPTools::signedImpactParameter3D(transTrk, direction, (*PV));
0145             ip3d = signed_IP3D.second.value();
0146             ip3d_err = signed_IP3D.second.error();
0147             TransverseImpactPointExtrapolator extrapolator(transTrk.field());
0148             GlobalPoint pos =
0149                 extrapolator.extrapolate(transTrk.impactPointState(), RecoVertex::convertPos(PV->position()))
0150                     .globalPosition();
0151             poca = reco::Vertex::Point(pos.x(), pos.y(), pos.z());
0152             AnalyticalImpactPointExtrapolator extrapolator3D(transTrk.field());
0153             GlobalPoint pos3d =
0154                 extrapolator3D.extrapolate(transTrk.impactPointState(), RecoVertex::convertPos(PV->position()))
0155                     .globalPosition();
0156             ip3d_poca = reco::Vertex::Point(pos3d.x(), pos3d.y(), pos3d.z());
0157           } else {
0158             dxy_err = track->d0Error();
0159             dxy = track->dxy(PV->position());
0160             ip3d_err = track->dzError();       //store dz, ip3d not available
0161             ip3d = track->dz(PV->position());  //store dz, ip3d not available
0162           }
0163         }
0164       }
0165       if (!SV.empty()) {
0166         reco::Vertex::CovarianceMatrix cov;
0167         reco::Vertex::Point v(SV.at(0)->x() - PV->x(), SV.at(0)->y() - PV->y(), SV.at(0)->z() - PV->z());
0168         for (int i = 0; i < reco::Vertex::dimension; i++) {
0169           for (int j = 0; j < reco::Vertex::dimension; j++) {
0170             cov(i, j) = SV.at(0)->covariance(i, j) + PV->covariance(i, j);
0171           }
0172         }
0173         GlobalVector direction(RefPFTau->px(), RefPFTau->py(), RefPFTau->pz());
0174         double vSig = SecondaryVertex::computeDist3d(*PV, *SV.at(0), direction, true).significance();
0175         reco::PFTauTransverseImpactParameter TIPV(poca, dxy, dxy_err, ip3d_poca, ip3d, ip3d_err, PV, v, vSig, SV.at(0));
0176         reco::PFTauTransverseImpactParameterRef TIPVRef =
0177             reco::PFTauTransverseImpactParameterRef(TIPRefProd_out, TIPCollection_out->size());
0178         TIPCollection_out->push_back(TIPV);
0179         AVPFTauTIP->setValue(iPFTau, TIPVRef);
0180       } else {
0181         reco::PFTauTransverseImpactParameter TIPV(poca, dxy, dxy_err, ip3d_poca, ip3d, ip3d_err, PV);
0182         reco::PFTauTransverseImpactParameterRef TIPVRef =
0183             reco::PFTauTransverseImpactParameterRef(TIPRefProd_out, TIPCollection_out->size());
0184         TIPCollection_out->push_back(TIPV);
0185         AVPFTauTIP->setValue(iPFTau, TIPVRef);
0186       }
0187     }
0188   }
0189   iEvent.put(std::move(TIPCollection_out), "PFTauTIP");
0190   iEvent.put(std::move(AVPFTauTIP));
0191 }
0192
0193 void PFTauTransverseImpactParameters::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0194   // PFTauTransverseImpactParameters
0195   edm::ParameterSetDescription desc;
0196   desc.add<edm::InputTag>("PFTauPVATag", edm::InputTag("PFTauPrimaryVertexProducer"));
0197   desc.add<bool>("useFullCalculation", false);
0198   desc.add<edm::InputTag>("PFTauTag", edm::InputTag("hpsPFTauProducer"));
0199   desc.add<edm::InputTag>("PFTauSVATag", edm::InputTag("PFTauSecondaryVertexProducer"));
0200   descriptions.add("PFTauTransverseImpactParameters", desc);
0201 }
0202
0203 DEFINE_FWK_MODULE(PFTauTransverseImpactParameters);