RecoTau/plugins/PFRecoTauChargedHadronFromGenericTrackPlugin.cc

0001 /*
0002  * PFRecoTauChargedHadronFromGenericTrackPlugin
0003  *
0004  * Build PFRecoTauChargedHadron objects
0005  * using tracks as input, from either collection of RECO/AOD reco::Tracks
0006  * (PFRecoTauChargedHadronFromTrackPlugin) or from a collection of MINIAOD
0007  * pat::PackedCandidates (PFRecoTauChargedHadronFromLostTrackPlugin), typically
0008  * using the 'lostTracks' collection
0009  *
0010  * Author: Christian Veelken, LLR
0011  *
0012  * inclusion of lost tracks based on original implementation
0013  * by Michal Bluj, NCBJ, Poland
0014  */
0015
0016 #include "RecoTauTag/RecoTau/interface/PFRecoTauChargedHadronPlugins.h"
0017
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019
0020 #include "MagneticField/Engine/interface/MagneticField.h"
0021 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0022
0023 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0024 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0025 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "DataFormats/VertexReco/interface/Vertex.h"
0029 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0030 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
0031 #include "DataFormats/JetReco/interface/PFJet.h"
0032 #include "DataFormats/Common/interface/AssociationMap.h"
0033 #include "DataFormats/Common/interface/Ptr.h"
0034 #include "DataFormats/Math/interface/deltaR.h"
0035
0036 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0037 #include "CommonTools/BaseParticlePropagator/interface/RawParticle.h"
0038
0039 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0040 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0041 #include "RecoTauTag/RecoTau/interface/RecoTauVertexAssociator.h"
0042 #include "RecoTauTag/RecoTau/interface/pfRecoTauChargedHadronAuxFunctions.h"
0043
0044 #include <TMath.h>
0045
0046 #include <memory>
0047 #include <cmath>
0048 #include <algorithm>
0049 #include <atomic>
0050
0051 namespace reco {
0052   namespace tau {
0053
0054     template <class TrackClass>
0055     class PFRecoTauChargedHadronFromGenericTrackPlugin : public PFRecoTauChargedHadronBuilderPlugin {
0056     public:
0057       explicit PFRecoTauChargedHadronFromGenericTrackPlugin(const edm::ParameterSet&, edm::ConsumesCollector&& iC);
0058       ~PFRecoTauChargedHadronFromGenericTrackPlugin() override;
0059       // Return type is unique_ptr<ChargedHadronVector>
0060       return_type operator()(const reco::Jet&) const override;
0061       // Hook to update PV information
0062       void beginEvent() override;
0063
0064     private:
0065       bool filterTrack(const edm::Handle<std::vector<TrackClass> >&, size_t iTrack) const;
0066       void setChargedHadronTrack(PFRecoTauChargedHadron& chargedHadron, const edm::Ptr<TrackClass>& track) const;
0067       double getTrackPtError(const TrackClass& track) const;
0068       XYZTLorentzVector getTrackPos(const TrackClass& track) const;
0069
0070       RecoTauVertexAssociator vertexAssociator_;
0071
0072       std::unique_ptr<RecoTauQualityCuts> qcuts_;
0073
0074       edm::InputTag srcTracks_;
0075       edm::EDGetTokenT<std::vector<TrackClass> > Tracks_token;
0076       const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0077       double dRcone_;
0078       bool dRconeLimitedToJetArea_;
0079
0080       double dRmergeNeutralHadron_;
0081       double dRmergePhoton_;
0082
0083       math::XYZVector magneticFieldStrength_;
0084
0085       static std::atomic<unsigned int> numWarnings_;
0086       static constexpr unsigned int maxWarnings_ = 3;
0087
0088       int verbosity_;
0089     };
0090
0091     template <class TrackClass>
0092     std::atomic<unsigned int> PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::numWarnings_{0};
0093
0094     template <class TrackClass>
0095     PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::PFRecoTauChargedHadronFromGenericTrackPlugin(
0096         const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0097         : PFRecoTauChargedHadronBuilderPlugin(pset, std::move(iC)),
0098           vertexAssociator_(pset.getParameter<edm::ParameterSet>("qualityCuts"), std::move(iC)),
0099           qcuts_(nullptr),
0100           magneticFieldToken_(iC.esConsumes()) {
0101       edm::ParameterSet qcuts_pset = pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts");
0102       qcuts_ = std::make_unique<RecoTauQualityCuts>(qcuts_pset);
0103
0104       srcTracks_ = pset.getParameter<edm::InputTag>("srcTracks");
0105       Tracks_token = iC.consumes<std::vector<TrackClass> >(srcTracks_);
0106       dRcone_ = pset.getParameter<double>("dRcone");
0107       dRconeLimitedToJetArea_ = pset.getParameter<bool>("dRconeLimitedToJetArea");
0108
0109       dRmergeNeutralHadron_ = pset.getParameter<double>("dRmergeNeutralHadron");
0110       dRmergePhoton_ = pset.getParameter<double>("dRmergePhoton");
0111
0112       verbosity_ = pset.getParameter<int>("verbosity");
0113     }
0114
0115     template <class TrackClass>
0116     PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::~PFRecoTauChargedHadronFromGenericTrackPlugin() {}
0117
0118     // Update the primary vertex
0119     template <class TrackClass>
0120     void PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::beginEvent() {
0121       vertexAssociator_.setEvent(*this->evt());
0122
0123       magneticFieldStrength_ = evtSetup()->getData(magneticFieldToken_).inTesla(GlobalPoint(0., 0., 0.));
0124     }
0125
0126     template <>
0127     bool PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track>::filterTrack(
0128         const edm::Handle<std::vector<reco::Track> >& tracks, size_t iTrack) const {
0129       // ignore tracks which fail quality cuts
0130       reco::TrackRef trackRef(tracks, iTrack);
0131       return qcuts_->filterTrack(trackRef);
0132     }
0133
0134     template <>
0135     bool PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate>::filterTrack(
0136         const edm::Handle<std::vector<pat::PackedCandidate> >& tracks, size_t iTrack) const {
0137       // ignore tracks which fail quality cuts
0138       const pat::PackedCandidate& cand = (*tracks)[iTrack];
0139       if (cand.charge() == 0)
0140         return false;
0141
0142       return qcuts_->filterChargedCand(cand);
0143     }
0144
0145     template <>
0146     void PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track>::setChargedHadronTrack(
0147         PFRecoTauChargedHadron& chargedHadron, const edm::Ptr<reco::Track>& track) const {
0148       chargedHadron.track_ = track;
0149     }
0150
0151     template <>
0152     void PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate>::setChargedHadronTrack(
0153         PFRecoTauChargedHadron& chargedHadron, const edm::Ptr<pat::PackedCandidate>& track) const {
0154       chargedHadron.lostTrackCandidate_ = track;
0155     }
0156
0157     template <>
0158     double PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track>::getTrackPtError(const reco::Track& track) const {
0159       return track.ptError();
0160     }
0161
0162     template <>
0163     double PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate>::getTrackPtError(
0164         const pat::PackedCandidate& cand) const {
0165       double trackPtError =
0166           0.06;  // MB: Approximate avarage track PtError by 2.5% (barrel), 4% (transition), 6% (endcaps) lostTracks w/o detailed track information available (after TRK-11-001)
0167       const reco::Track* track(cand.bestTrack());
0168       if (track != nullptr) {
0169         trackPtError = track->ptError();
0170       } else {
0171         if (std::abs(cand.eta()) < 0.9)
0172           trackPtError = 0.025;
0173         else if (std::abs(cand.eta()) < 1.4)
0174           trackPtError = 0.04;
0175       }
0176       return trackPtError;
0177     }
0178
0179     template <>
0180     XYZTLorentzVector PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track>::getTrackPos(
0181         const reco::Track& track) const {
0182       return XYZTLorentzVector(track.referencePoint().x(), track.referencePoint().y(), track.referencePoint().z(), 0.);
0183     }
0184
0185     template <>
0186     XYZTLorentzVector PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate>::getTrackPos(
0187         const pat::PackedCandidate& track) const {
0188       return XYZTLorentzVector(track.vertex().x(), track.vertex().y(), track.vertex().z(), 0.);
0189     }
0190
0191     namespace {
0192       struct Candidate_withDistance {
0193         reco::CandidatePtr pfCandidate_;
0194         double distance_;
0195       };
0196
0197       bool isSmallerDistance(const Candidate_withDistance& cand1, const Candidate_withDistance& cand2) {
0198         return (cand1.distance_ < cand2.distance_);
0199       }
0200     }  // namespace
0201
0202     template <class TrackClass>
0203     typename PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::return_type
0204     PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::operator()(const reco::Jet& jet) const {
0205       if (verbosity_) {
0206         edm::LogPrint("TauChHFromTrack") << "<PFRecoTauChargedHadronFromGenericTrackPlugin::operator()>:";
0207         edm::LogPrint("TauChHFromTrack") << " pluginName = " << name();
0208       }
0209
0210       ChargedHadronVector output;
0211
0212       const edm::Event& evt = (*this->evt());
0213
0214       edm::Handle<std::vector<TrackClass> > tracks;
0215       evt.getByToken(Tracks_token, tracks);
0216
0217       qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0218       float jEta = jet.eta();
0219       float jPhi = jet.phi();
0220       size_t numTracks = tracks->size();
0221       for (size_t iTrack = 0; iTrack < numTracks; ++iTrack) {
0222         const TrackClass& track = (*tracks)[iTrack];
0223
0224         // consider tracks in vicinity of tau-jet candidate only
0225         double dR = deltaR(track.eta(), track.phi(), jEta, jPhi);
0226         double dRmatch = dRcone_;
0227         if (dRconeLimitedToJetArea_) {
0228           double jetArea = jet.jetArea();
0229           if (jetArea > 0.) {
0230             dRmatch = std::min(dRmatch, sqrt(jetArea / M_PI));
0231           } else {
0232             if (numWarnings_ < maxWarnings_) {
0233               edm::LogInfo("PFRecoTauChargedHadronFromGenericTrackPlugin::operator()")
0234                   << "Jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi()
0235                   << " has area = " << jetArea << " !!" << std::endl;
0236               ++numWarnings_;
0237             }
0238             dRmatch = 0.1;
0239           }
0240         }
0241         if (dR > dRmatch)
0242           continue;
0243
0244         if (!this->filterTrack(tracks, iTrack))
0245           continue;
0246
0247         reco::Candidate::Charge trackCharge_int = 0;
0248         if (track.charge() > 0.)
0249           trackCharge_int = +1;
0250         else if (track.charge() < 0.)
0251           trackCharge_int = -1;
0252
0253         const double chargedPionMass = 0.13957;  // GeV
0254         double chargedPionP = track.p();
0255         double chargedPionEn = TMath::Sqrt(chargedPionP * chargedPionP + chargedPionMass * chargedPionMass);
0256         reco::Candidate::LorentzVector chargedPionP4(track.px(), track.py(), track.pz(), chargedPionEn);
0257
0258         reco::Vertex::Point vtx(0., 0., 0.);
0259         if (vertexAssociator_.associatedVertex(jet).isNonnull())
0260           vtx = vertexAssociator_.associatedVertex(jet)->position();
0261
0262         auto chargedHadron = std::make_unique<PFRecoTauChargedHadron>(
0263             trackCharge_int, chargedPionP4, vtx, 0, true, PFRecoTauChargedHadron::kTrack);
0264
0265         setChargedHadronTrack(*chargedHadron, edm::Ptr<TrackClass>(tracks, iTrack));
0266
0267         // CV: Take code for propagating track to ECAL entrance
0268         //     from RecoParticleFlow/PFTracking/src/PFTrackTransformer.cc
0269         //     to make sure propagation is done in the same way as for charged PFCandidates.
0270         //
0271         //     The following replacements need to be made
0272         //       outerMomentum -> momentum
0273         //       outerPosition -> referencePoint
0274         //     in order to run on AOD input
0275         //    (outerMomentum and outerPosition require access to reco::TrackExtra objects, which are available in RECO only)
0276         //
0277         XYZTLorentzVector chargedPionPos(getTrackPos(track));
0278         RawParticle p(chargedPionP4, chargedPionPos);
0279         p.setCharge(track.charge());
0280         BaseParticlePropagator trackPropagator(p, 0., 0., magneticFieldStrength_.z());
0281         trackPropagator.propagateToEcalEntrance(false);
0282         if (trackPropagator.getSuccess() != 0) {
0283           chargedHadron->positionAtECALEntrance_ = trackPropagator.particle().vertex();
0284         } else {
0285           if (chargedPionP4.pt() > 2. and std::abs(chargedPionP4.eta()) < 3.) {
0286             edm::LogWarning("PFRecoTauChargedHadronFromGenericTrackPlugin::operator()")
0287                 << "Failed to propagate track: Pt = " << track.pt() << ", eta = " << track.eta()
0288                 << ", phi = " << track.phi() << " to ECAL entrance !!" << std::endl;
0289           }
0290           chargedHadron->positionAtECALEntrance_ = math::XYZPointF(0., 0., 0.);
0291         }
0292
0293         std::vector<Candidate_withDistance> neutralJetConstituents_withDistance;
0294         for (const auto& jetConstituent : jet.daughterPtrVector()) {
0295           int pdgId = jetConstituent->pdgId();
0296           if (!(pdgId == 130 || pdgId == 22))
0297             continue;
0298           double dR = deltaR(atECALEntrance(&*jetConstituent, magneticFieldStrength_.z()),
0299                              chargedHadron->positionAtECALEntrance_);
0300           double dRmerge = -1.;
0301           if (pdgId == 130)
0302             dRmerge = dRmergeNeutralHadron_;
0303           else if (pdgId == 22)
0304             dRmerge = dRmergePhoton_;
0305           if (dR < dRmerge) {
0306             Candidate_withDistance jetConstituent_withDistance;
0307             jetConstituent_withDistance.pfCandidate_ = jetConstituent;
0308             jetConstituent_withDistance.distance_ = dR;
0309             neutralJetConstituents_withDistance.push_back(jetConstituent_withDistance);
0310             chargedHadron->addDaughter(jetConstituent);
0311           }
0312         }
0313         std::sort(
0314             neutralJetConstituents_withDistance.begin(), neutralJetConstituents_withDistance.end(), isSmallerDistance);
0315
0316         const double caloResolutionCoeff =
0317             1.0;  // CV: approximate ECAL + HCAL calorimeter resolution for hadrons by 100%*sqrt(E)
0318         double resolutionTrackP = track.p() * (getTrackPtError(track) / track.pt());
0319         double neutralEnSum = 0.;
0320         for (std::vector<Candidate_withDistance>::const_iterator nextNeutral =
0321                  neutralJetConstituents_withDistance.begin();
0322              nextNeutral != neutralJetConstituents_withDistance.end();
0323              ++nextNeutral) {
0324           double nextNeutralEn = nextNeutral->pfCandidate_->energy();
0325           double resolutionCaloEn = caloResolutionCoeff * sqrt(neutralEnSum + nextNeutralEn);
0326           double resolution = sqrt(resolutionTrackP * resolutionTrackP + resolutionCaloEn * resolutionCaloEn);
0327           if ((neutralEnSum + nextNeutralEn) < (track.p() + 2. * resolution)) {
0328             chargedHadron->neutralPFCandidates_.push_back(nextNeutral->pfCandidate_);
0329             neutralEnSum += nextNeutralEn;
0330           } else {
0331             break;
0332           }
0333         }
0334
0335         setChargedHadronP4(*chargedHadron);
0336
0337         if (verbosity_) {
0338           edm::LogPrint("TauChHFromTrack") << *chargedHadron;
0339         }
0340
0341         output.push_back(std::move(chargedHadron));
0342       }
0343
0344       return output;
0345     }
0346
0347     typedef PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track> PFRecoTauChargedHadronFromTrackPlugin;
0348     typedef PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate> PFRecoTauChargedHadronFromLostTrackPlugin;
0349
0350   }  // namespace tau
0351 }  // namespace reco
0352
0353 #include "FWCore/Framework/interface/MakerMacros.h"
0354
0355 DEFINE_EDM_PLUGIN(PFRecoTauChargedHadronBuilderPluginFactory,
0356                   reco::tau::PFRecoTauChargedHadronFromTrackPlugin,
0357                   "PFRecoTauChargedHadronFromTrackPlugin");
0358 DEFINE_EDM_PLUGIN(PFRecoTauChargedHadronBuilderPluginFactory,
0359                   reco::tau::PFRecoTauChargedHadronFromLostTrackPlugin,
0360                   "PFRecoTauChargedHadronFromLostTrackPlugin");