RecoTau/plugins/PFRecoTauChargedHadronFromPFCandidatePlugin.cc

0001 /*
0002  * PFRecoTauChargedHadronFromPFCandidatePlugin
0003  *
0004  * Build PFRecoTauChargedHadron objects
0005  * using charged PFCandidates and/or PFNeutralHadrons as input
0006  *
0007  * Author: Christian Veelken, LLR
0008  *
0009  */
0010
0011 #include "RecoTauTag/RecoTau/interface/PFRecoTauChargedHadronPlugins.h"
0012
0013 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0014 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0020 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0021 #include "DataFormats/MuonReco/interface/Muon.h"
0022 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0023 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
0024 #include "DataFormats/JetReco/interface/Jet.h"
0025 #include "DataFormats/Math/interface/deltaR.h"
0026 #include "DataFormats/Common/interface/RefToPtr.h"
0027
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "MagneticField/Engine/interface/MagneticField.h"
0030 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0031
0032 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0033 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0034 #include "RecoTauTag/RecoTau/interface/RecoTauVertexAssociator.h"
0035 #include "RecoTauTag/RecoTau/interface/pfRecoTauChargedHadronAuxFunctions.h"
0036
0037 #include <memory>
0038
0039 namespace reco {
0040
0041   namespace tau {
0042
0043     class PFRecoTauChargedHadronFromPFCandidatePlugin : public PFRecoTauChargedHadronBuilderPlugin {
0044     public:
0045       explicit PFRecoTauChargedHadronFromPFCandidatePlugin(const edm::ParameterSet&, edm::ConsumesCollector&& iC);
0046       ~PFRecoTauChargedHadronFromPFCandidatePlugin() override {}
0047       // Return type is unique_ptr<ChargedHadronVector>
0048       return_type operator()(const reco::Jet&) const override;
0049       // Hook to update PV information
0050       void beginEvent() override;
0051
0052     private:
0053       typedef std::vector<reco::CandidatePtr> CandPtrs;
0054
0055       RecoTauVertexAssociator vertexAssociator_;
0056
0057       std::unique_ptr<RecoTauQualityCuts> qcuts_;
0058
0059       std::vector<int> inputParticleIds_;  // type of candidates to clusterize
0060
0061       double dRmergeNeutralHadronWrtChargedHadron_;
0062       double dRmergeNeutralHadronWrtNeutralHadron_;
0063       double dRmergeNeutralHadronWrtElectron_;
0064       double dRmergeNeutralHadronWrtOther_;
0065       int minBlockElementMatchesNeutralHadron_;
0066       int maxUnmatchedBlockElementsNeutralHadron_;
0067       double dRmergePhotonWrtChargedHadron_;
0068       double dRmergePhotonWrtNeutralHadron_;
0069       double dRmergePhotonWrtElectron_;
0070       double dRmergePhotonWrtOther_;
0071       int minBlockElementMatchesPhoton_;
0072       int maxUnmatchedBlockElementsPhoton_;
0073       double minMergeNeutralHadronEt_;
0074       double minMergeGammaEt_;
0075       double minMergeChargedHadronPt_;
0076
0077       const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> bFieldToken_;
0078       double bField_;
0079
0080       int verbosity_;
0081     };
0082
0083     PFRecoTauChargedHadronFromPFCandidatePlugin::PFRecoTauChargedHadronFromPFCandidatePlugin(
0084         const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0085         : PFRecoTauChargedHadronBuilderPlugin(pset, std::move(iC)),
0086           vertexAssociator_(pset.getParameter<edm::ParameterSet>("qualityCuts"), std::move(iC)),
0087           qcuts_(nullptr),
0088           bFieldToken_(iC.esConsumes()) {
0089       edm::ParameterSet qcuts_pset = pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts");
0090       qcuts_ = std::make_unique<RecoTauQualityCuts>(qcuts_pset);
0091
0092       inputParticleIds_ = pset.getParameter<std::vector<int> >("chargedHadronCandidatesParticleIds");
0093
0094       dRmergeNeutralHadronWrtChargedHadron_ = pset.getParameter<double>("dRmergeNeutralHadronWrtChargedHadron");
0095       dRmergeNeutralHadronWrtNeutralHadron_ = pset.getParameter<double>("dRmergeNeutralHadronWrtNeutralHadron");
0096       dRmergeNeutralHadronWrtElectron_ = pset.getParameter<double>("dRmergeNeutralHadronWrtElectron");
0097       dRmergeNeutralHadronWrtOther_ = pset.getParameter<double>("dRmergeNeutralHadronWrtOther");
0098       minBlockElementMatchesNeutralHadron_ = pset.getParameter<int>("minBlockElementMatchesNeutralHadron");
0099       maxUnmatchedBlockElementsNeutralHadron_ = pset.getParameter<int>("maxUnmatchedBlockElementsNeutralHadron");
0100       dRmergePhotonWrtChargedHadron_ = pset.getParameter<double>("dRmergePhotonWrtChargedHadron");
0101       dRmergePhotonWrtNeutralHadron_ = pset.getParameter<double>("dRmergePhotonWrtNeutralHadron");
0102       dRmergePhotonWrtElectron_ = pset.getParameter<double>("dRmergePhotonWrtElectron");
0103       dRmergePhotonWrtOther_ = pset.getParameter<double>("dRmergePhotonWrtOther");
0104       minBlockElementMatchesPhoton_ = pset.getParameter<int>("minBlockElementMatchesPhoton");
0105       maxUnmatchedBlockElementsPhoton_ = pset.getParameter<int>("maxUnmatchedBlockElementsPhoton");
0106       minMergeNeutralHadronEt_ = pset.getParameter<double>("minMergeNeutralHadronEt");
0107       minMergeGammaEt_ = pset.getParameter<double>("minMergeGammaEt");
0108       minMergeChargedHadronPt_ = pset.getParameter<double>("minMergeChargedHadronPt");
0109
0110       verbosity_ = pset.getParameter<int>("verbosity");
0111     }
0112
0113     // Update the primary vertex
0114     void PFRecoTauChargedHadronFromPFCandidatePlugin::beginEvent() {
0115       vertexAssociator_.setEvent(*evt());
0116
0117       bField_ = evtSetup()->getData(bFieldToken_).inTesla(GlobalPoint(0, 0, 0)).z();
0118     }
0119
0120     namespace {
0121       bool isMatchedByBlockElement(const reco::PFCandidate& pfCandidate1,
0122                                    const reco::PFCandidate& pfCandidate2,
0123                                    int minMatches1,
0124                                    int minMatches2,
0125                                    int maxUnmatchedBlockElements1plus2) {
0126         reco::PFCandidate::ElementsInBlocks blockElements1 = pfCandidate1.elementsInBlocks();
0127         int numBlocks1 = blockElements1.size();
0128         reco::PFCandidate::ElementsInBlocks blockElements2 = pfCandidate2.elementsInBlocks();
0129         int numBlocks2 = blockElements2.size();
0130         int numBlocks_matched = 0;
0131         for (reco::PFCandidate::ElementsInBlocks::const_iterator blockElement1 = blockElements1.begin();
0132              blockElement1 != blockElements1.end();
0133              ++blockElement1) {
0134           bool isMatched = false;
0135           for (reco::PFCandidate::ElementsInBlocks::const_iterator blockElement2 = blockElements2.begin();
0136                blockElement2 != blockElements2.end();
0137                ++blockElement2) {
0138             if (blockElement1->first.id() == blockElement2->first.id() &&
0139                 blockElement1->first.key() == blockElement2->first.key() &&
0140                 blockElement1->second == blockElement2->second) {
0141               isMatched = true;
0142             }
0143           }
0144           if (isMatched)
0145             ++numBlocks_matched;
0146         }
0147         assert(numBlocks_matched <= numBlocks1);
0148         assert(numBlocks_matched <= numBlocks2);
0149         if (numBlocks_matched >= minMatches1 && numBlocks_matched >= minMatches2 &&
0150             ((numBlocks1 - numBlocks_matched) + (numBlocks2 - numBlocks_matched)) <= maxUnmatchedBlockElements1plus2) {
0151           return true;
0152         } else {
0153           return false;
0154         }
0155       }
0156     }  // namespace
0157
0158     PFRecoTauChargedHadronFromPFCandidatePlugin::return_type PFRecoTauChargedHadronFromPFCandidatePlugin::operator()(
0159         const reco::Jet& jet) const {
0160       if (verbosity_) {
0161         edm::LogPrint("TauChHadronFromPF") << "<PFRecoTauChargedHadronFromPFCandidatePlugin::operator()>:";
0162         edm::LogPrint("TauChHadronFromPF") << " pluginName = " << name();
0163       }
0164
0165       ChargedHadronVector output;
0166
0167       // Get the candidates passing our quality cuts
0168       qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0169       CandPtrs candsVector = qcuts_->filterCandRefs(pfCandidates(jet, inputParticleIds_));
0170
0171       for (CandPtrs::iterator cand = candsVector.begin(); cand != candsVector.end(); ++cand) {
0172         if (verbosity_) {
0173           edm::LogPrint("TauChHadronFromPF")
0174               << "processing PFCandidate: Pt = " << (*cand)->pt() << ", eta = " << (*cand)->eta()
0175               << ", phi = " << (*cand)->phi() << " (pdgId = " << (*cand)->pdgId() << ", charge = " << (*cand)->charge()
0176               << ")";
0177         }
0178
0179         PFRecoTauChargedHadron::PFRecoTauChargedHadronAlgorithm algo = PFRecoTauChargedHadron::kUndefined;
0180         if (std::abs((*cand)->charge()) > 0.5)
0181           algo = PFRecoTauChargedHadron::kChargedPFCandidate;
0182         else
0183           algo = PFRecoTauChargedHadron::kPFNeutralHadron;
0184         auto chargedHadron = std::make_unique<PFRecoTauChargedHadron>(**cand, algo);
0185
0186         const reco::PFCandidate* pfCand = dynamic_cast<const reco::PFCandidate*>(&**cand);
0187         if (pfCand) {
0188           if (pfCand->trackRef().isNonnull())
0189             chargedHadron->track_ = edm::refToPtr(pfCand->trackRef());
0190           else if (pfCand->muonRef().isNonnull() && pfCand->muonRef()->innerTrack().isNonnull())
0191             chargedHadron->track_ = edm::refToPtr(pfCand->muonRef()->innerTrack());
0192           else if (pfCand->muonRef().isNonnull() && pfCand->muonRef()->globalTrack().isNonnull())
0193             chargedHadron->track_ = edm::refToPtr(pfCand->muonRef()->globalTrack());
0194           else if (pfCand->muonRef().isNonnull() && pfCand->muonRef()->outerTrack().isNonnull())
0195             chargedHadron->track_ = edm::refToPtr(pfCand->muonRef()->outerTrack());
0196           else if (pfCand->gsfTrackRef().isNonnull())
0197             chargedHadron->track_ = edm::refToPtr(pfCand->gsfTrackRef());
0198         }  // TauReco@MiniAOD: Tracks only available dynamically, so no possiblity to save ref here; checked by code downstream
0199
0200         chargedHadron->positionAtECALEntrance_ = atECALEntrance(&**cand, bField_);
0201         chargedHadron->chargedPFCandidate_ = (*cand);
0202         chargedHadron->addDaughter(*cand);
0203
0204         int pdgId = std::abs(chargedHadron->chargedPFCandidate_->pdgId());
0205
0206         if (chargedHadron->pt() > minMergeChargedHadronPt_) {
0207           for (const auto& jetConstituent : jet.daughterPtrVector()) {
0208             // CV: take care of not double-counting energy in case "charged" PFCandidate is in fact a PFNeutralHadron
0209             if (jetConstituent == chargedHadron->chargedPFCandidate_)
0210               continue;
0211
0212             int jetConstituentPdgId = std::abs(jetConstituent->pdgId());
0213             if (!(jetConstituentPdgId == 130 || jetConstituentPdgId == 22))
0214               continue;
0215
0216             double dR = deltaR(atECALEntrance(jetConstituent.get(), bField_),
0217                                atECALEntrance(chargedHadron->chargedPFCandidate_.get(), bField_));
0218             double dRmerge = -1.;
0219             int minBlockElementMatches = 1000;
0220             int maxUnmatchedBlockElements = 0;
0221             double minMergeEt = 1.e+6;
0222             if (jetConstituentPdgId == 130) {
0223               if (pdgId == 211)
0224                 dRmerge = dRmergeNeutralHadronWrtChargedHadron_;
0225               else if (pdgId == 130)
0226                 dRmerge = dRmergeNeutralHadronWrtNeutralHadron_;
0227               else if (pdgId == 11)
0228                 dRmerge = dRmergeNeutralHadronWrtElectron_;
0229               else
0230                 dRmerge = dRmergeNeutralHadronWrtOther_;
0231               minBlockElementMatches = minBlockElementMatchesNeutralHadron_;
0232               maxUnmatchedBlockElements = maxUnmatchedBlockElementsNeutralHadron_;
0233               minMergeEt = minMergeNeutralHadronEt_;
0234             } else if (jetConstituentPdgId == 22) {
0235               if (pdgId == 211)
0236                 dRmerge = dRmergePhotonWrtChargedHadron_;
0237               else if (pdgId == 130)
0238                 dRmerge = dRmergePhotonWrtNeutralHadron_;
0239               else if (pdgId == 11)
0240                 dRmerge = dRmergePhotonWrtElectron_;
0241               else
0242                 dRmerge = dRmergePhotonWrtOther_;
0243               minBlockElementMatches = minBlockElementMatchesPhoton_;
0244               maxUnmatchedBlockElements = maxUnmatchedBlockElementsPhoton_;
0245               minMergeEt = minMergeGammaEt_;
0246             }
0247
0248             if (jetConstituent->et() > minMergeEt) {
0249               if (dR < dRmerge) {
0250                 chargedHadron->neutralPFCandidates_.push_back(jetConstituent);
0251                 chargedHadron->addDaughter(jetConstituent);
0252               } else {
0253                 // TauReco@MiniAOD: No access to PF blocks at MiniAOD level, but the code below seems to have very minor impact
0254                 const reco::PFCandidate* pfJetConstituent =
0255                     dynamic_cast<const reco::PFCandidate*>(jetConstituent.get());
0256                 if (pfCand != nullptr && pfJetConstituent != nullptr) {
0257                   if (isMatchedByBlockElement(*pfJetConstituent,
0258                                               *pfCand,
0259                                               minBlockElementMatches,
0260                                               minBlockElementMatches,
0261                                               maxUnmatchedBlockElements)) {
0262                     chargedHadron->neutralPFCandidates_.push_back(jetConstituent);
0263                     chargedHadron->addDaughter(jetConstituent);
0264                   }
0265                 }
0266               }
0267             }
0268           }
0269         }
0270
0271         setChargedHadronP4(*chargedHadron);
0272
0273         if (verbosity_) {
0274           edm::LogPrint("TauChHadronFromPF") << *chargedHadron;
0275         }
0276         // Update the vertex
0277         if (chargedHadron->daughterPtr(0).isNonnull())
0278           chargedHadron->setVertex(chargedHadron->daughterPtr(0)->vertex());
0279         output.push_back(std::move(chargedHadron));
0280       }
0281
0282       return output;
0283     }
0284
0285   }  // namespace tau
0286 }  // namespace reco
0287
0288 #include "FWCore/Framework/interface/MakerMacros.h"
0289
0290 DEFINE_EDM_PLUGIN(PFRecoTauChargedHadronBuilderPluginFactory,
0291                   reco::tau::PFRecoTauChargedHadronFromPFCandidatePlugin,
0292                   "PFRecoTauChargedHadronFromPFCandidatePlugin");