Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-25 04:55:28

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       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_ = new 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       delete qcuts_;
0118     }
0119 
0120     // Update the primary vertex
0121     template <class TrackClass>
0122     void PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::beginEvent() {
0123       vertexAssociator_.setEvent(*this->evt());
0124 
0125       magneticFieldStrength_ = evtSetup()->getData(magneticFieldToken_).inTesla(GlobalPoint(0., 0., 0.));
0126     }
0127 
0128     template <>
0129     bool PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track>::filterTrack(
0130         const edm::Handle<std::vector<reco::Track> >& tracks, size_t iTrack) const {
0131       // ignore tracks which fail quality cuts
0132       reco::TrackRef trackRef(tracks, iTrack);
0133       return qcuts_->filterTrack(trackRef);
0134     }
0135 
0136     template <>
0137     bool PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate>::filterTrack(
0138         const edm::Handle<std::vector<pat::PackedCandidate> >& tracks, size_t iTrack) const {
0139       // ignore tracks which fail quality cuts
0140       const pat::PackedCandidate& cand = (*tracks)[iTrack];
0141       if (cand.charge() == 0)
0142         return false;
0143 
0144       return qcuts_->filterChargedCand(cand);
0145     }
0146 
0147     template <>
0148     void PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track>::setChargedHadronTrack(
0149         PFRecoTauChargedHadron& chargedHadron, const edm::Ptr<reco::Track>& track) const {
0150       chargedHadron.track_ = track;
0151     }
0152 
0153     template <>
0154     void PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate>::setChargedHadronTrack(
0155         PFRecoTauChargedHadron& chargedHadron, const edm::Ptr<pat::PackedCandidate>& track) const {
0156       chargedHadron.lostTrackCandidate_ = track;
0157     }
0158 
0159     template <>
0160     double PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track>::getTrackPtError(const reco::Track& track) const {
0161       return track.ptError();
0162     }
0163 
0164     template <>
0165     double PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate>::getTrackPtError(
0166         const pat::PackedCandidate& cand) const {
0167       double trackPtError =
0168           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)
0169       const reco::Track* track(cand.bestTrack());
0170       if (track != nullptr) {
0171         trackPtError = track->ptError();
0172       } else {
0173         if (std::abs(cand.eta()) < 0.9)
0174           trackPtError = 0.025;
0175         else if (std::abs(cand.eta()) < 1.4)
0176           trackPtError = 0.04;
0177       }
0178       return trackPtError;
0179     }
0180 
0181     template <>
0182     XYZTLorentzVector PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track>::getTrackPos(
0183         const reco::Track& track) const {
0184       return XYZTLorentzVector(track.referencePoint().x(), track.referencePoint().y(), track.referencePoint().z(), 0.);
0185     }
0186 
0187     template <>
0188     XYZTLorentzVector PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate>::getTrackPos(
0189         const pat::PackedCandidate& track) const {
0190       return XYZTLorentzVector(track.vertex().x(), track.vertex().y(), track.vertex().z(), 0.);
0191     }
0192 
0193     namespace {
0194       struct Candidate_withDistance {
0195         reco::CandidatePtr pfCandidate_;
0196         double distance_;
0197       };
0198 
0199       bool isSmallerDistance(const Candidate_withDistance& cand1, const Candidate_withDistance& cand2) {
0200         return (cand1.distance_ < cand2.distance_);
0201       }
0202     }  // namespace
0203 
0204     template <class TrackClass>
0205     typename PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::return_type
0206     PFRecoTauChargedHadronFromGenericTrackPlugin<TrackClass>::operator()(const reco::Jet& jet) const {
0207       if (verbosity_) {
0208         edm::LogPrint("TauChHFromTrack") << "<PFRecoTauChargedHadronFromGenericTrackPlugin::operator()>:";
0209         edm::LogPrint("TauChHFromTrack") << " pluginName = " << name();
0210       }
0211 
0212       ChargedHadronVector output;
0213 
0214       const edm::Event& evt = (*this->evt());
0215 
0216       edm::Handle<std::vector<TrackClass> > tracks;
0217       evt.getByToken(Tracks_token, tracks);
0218 
0219       qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0220       float jEta = jet.eta();
0221       float jPhi = jet.phi();
0222       size_t numTracks = tracks->size();
0223       for (size_t iTrack = 0; iTrack < numTracks; ++iTrack) {
0224         const TrackClass& track = (*tracks)[iTrack];
0225 
0226         // consider tracks in vicinity of tau-jet candidate only
0227         double dR = deltaR(track.eta(), track.phi(), jEta, jPhi);
0228         double dRmatch = dRcone_;
0229         if (dRconeLimitedToJetArea_) {
0230           double jetArea = jet.jetArea();
0231           if (jetArea > 0.) {
0232             dRmatch = std::min(dRmatch, sqrt(jetArea / M_PI));
0233           } else {
0234             if (numWarnings_ < maxWarnings_) {
0235               edm::LogInfo("PFRecoTauChargedHadronFromGenericTrackPlugin::operator()")
0236                   << "Jet: Pt = " << jet.pt() << ", eta = " << jet.eta() << ", phi = " << jet.phi()
0237                   << " has area = " << jetArea << " !!" << std::endl;
0238               ++numWarnings_;
0239             }
0240             dRmatch = 0.1;
0241           }
0242         }
0243         if (dR > dRmatch)
0244           continue;
0245 
0246         if (!this->filterTrack(tracks, iTrack))
0247           continue;
0248 
0249         reco::Candidate::Charge trackCharge_int = 0;
0250         if (track.charge() > 0.)
0251           trackCharge_int = +1;
0252         else if (track.charge() < 0.)
0253           trackCharge_int = -1;
0254 
0255         const double chargedPionMass = 0.13957;  // GeV
0256         double chargedPionP = track.p();
0257         double chargedPionEn = TMath::Sqrt(chargedPionP * chargedPionP + chargedPionMass * chargedPionMass);
0258         reco::Candidate::LorentzVector chargedPionP4(track.px(), track.py(), track.pz(), chargedPionEn);
0259 
0260         reco::Vertex::Point vtx(0., 0., 0.);
0261         if (vertexAssociator_.associatedVertex(jet).isNonnull())
0262           vtx = vertexAssociator_.associatedVertex(jet)->position();
0263 
0264         std::unique_ptr<PFRecoTauChargedHadron> chargedHadron(
0265             new PFRecoTauChargedHadron(trackCharge_int, chargedPionP4, vtx, 0, true, PFRecoTauChargedHadron::kTrack));
0266 
0267         setChargedHadronTrack(*chargedHadron, edm::Ptr<TrackClass>(tracks, iTrack));
0268 
0269         // CV: Take code for propagating track to ECAL entrance
0270         //     from RecoParticleFlow/PFTracking/src/PFTrackTransformer.cc
0271         //     to make sure propagation is done in the same way as for charged PFCandidates.
0272         //
0273         //     The following replacements need to be made
0274         //       outerMomentum -> momentum
0275         //       outerPosition -> referencePoint
0276         //     in order to run on AOD input
0277         //    (outerMomentum and outerPosition require access to reco::TrackExtra objects, which are available in RECO only)
0278         //
0279         XYZTLorentzVector chargedPionPos(getTrackPos(track));
0280         RawParticle p(chargedPionP4, chargedPionPos);
0281         p.setCharge(track.charge());
0282         BaseParticlePropagator trackPropagator(p, 0., 0., magneticFieldStrength_.z());
0283         trackPropagator.propagateToEcalEntrance(false);
0284         if (trackPropagator.getSuccess() != 0) {
0285           chargedHadron->positionAtECALEntrance_ = trackPropagator.particle().vertex();
0286         } else {
0287           if (chargedPionP4.pt() > 2. and std::abs(chargedPionP4.eta()) < 3.) {
0288             edm::LogWarning("PFRecoTauChargedHadronFromGenericTrackPlugin::operator()")
0289                 << "Failed to propagate track: Pt = " << track.pt() << ", eta = " << track.eta()
0290                 << ", phi = " << track.phi() << " to ECAL entrance !!" << std::endl;
0291           }
0292           chargedHadron->positionAtECALEntrance_ = math::XYZPointF(0., 0., 0.);
0293         }
0294 
0295         std::vector<Candidate_withDistance> neutralJetConstituents_withDistance;
0296         for (const auto& jetConstituent : jet.daughterPtrVector()) {
0297           int pdgId = jetConstituent->pdgId();
0298           if (!(pdgId == 130 || pdgId == 22))
0299             continue;
0300           double dR = deltaR(atECALEntrance(&*jetConstituent, magneticFieldStrength_.z()),
0301                              chargedHadron->positionAtECALEntrance_);
0302           double dRmerge = -1.;
0303           if (pdgId == 130)
0304             dRmerge = dRmergeNeutralHadron_;
0305           else if (pdgId == 22)
0306             dRmerge = dRmergePhoton_;
0307           if (dR < dRmerge) {
0308             Candidate_withDistance jetConstituent_withDistance;
0309             jetConstituent_withDistance.pfCandidate_ = jetConstituent;
0310             jetConstituent_withDistance.distance_ = dR;
0311             neutralJetConstituents_withDistance.push_back(jetConstituent_withDistance);
0312             chargedHadron->addDaughter(jetConstituent);
0313           }
0314         }
0315         std::sort(
0316             neutralJetConstituents_withDistance.begin(), neutralJetConstituents_withDistance.end(), isSmallerDistance);
0317 
0318         const double caloResolutionCoeff =
0319             1.0;  // CV: approximate ECAL + HCAL calorimeter resolution for hadrons by 100%*sqrt(E)
0320         double resolutionTrackP = track.p() * (getTrackPtError(track) / track.pt());
0321         double neutralEnSum = 0.;
0322         for (std::vector<Candidate_withDistance>::const_iterator nextNeutral =
0323                  neutralJetConstituents_withDistance.begin();
0324              nextNeutral != neutralJetConstituents_withDistance.end();
0325              ++nextNeutral) {
0326           double nextNeutralEn = nextNeutral->pfCandidate_->energy();
0327           double resolutionCaloEn = caloResolutionCoeff * sqrt(neutralEnSum + nextNeutralEn);
0328           double resolution = sqrt(resolutionTrackP * resolutionTrackP + resolutionCaloEn * resolutionCaloEn);
0329           if ((neutralEnSum + nextNeutralEn) < (track.p() + 2. * resolution)) {
0330             chargedHadron->neutralPFCandidates_.push_back(nextNeutral->pfCandidate_);
0331             neutralEnSum += nextNeutralEn;
0332           } else {
0333             break;
0334           }
0335         }
0336 
0337         setChargedHadronP4(*chargedHadron);
0338 
0339         if (verbosity_) {
0340           edm::LogPrint("TauChHFromTrack") << *chargedHadron;
0341         }
0342 
0343         output.push_back(std::move(chargedHadron));
0344       }
0345 
0346       return output;
0347     }
0348 
0349     typedef PFRecoTauChargedHadronFromGenericTrackPlugin<reco::Track> PFRecoTauChargedHadronFromTrackPlugin;
0350     typedef PFRecoTauChargedHadronFromGenericTrackPlugin<pat::PackedCandidate> PFRecoTauChargedHadronFromLostTrackPlugin;
0351 
0352   }  // namespace tau
0353 }  // namespace reco
0354 
0355 #include "FWCore/Framework/interface/MakerMacros.h"
0356 
0357 DEFINE_EDM_PLUGIN(PFRecoTauChargedHadronBuilderPluginFactory,
0358                   reco::tau::PFRecoTauChargedHadronFromTrackPlugin,
0359                   "PFRecoTauChargedHadronFromTrackPlugin");
0360 DEFINE_EDM_PLUGIN(PFRecoTauChargedHadronBuilderPluginFactory,
0361                   reco::tau::PFRecoTauChargedHadronFromLostTrackPlugin,
0362                   "PFRecoTauChargedHadronFromLostTrackPlugin");