Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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       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_ = new 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     PFRecoTauChargedHadronFromPFCandidatePlugin::~PFRecoTauChargedHadronFromPFCandidatePlugin() { delete qcuts_; }
0114 
0115     // Update the primary vertex
0116     void PFRecoTauChargedHadronFromPFCandidatePlugin::beginEvent() {
0117       vertexAssociator_.setEvent(*evt());
0118 
0119       bField_ = evtSetup()->getData(bFieldToken_).inTesla(GlobalPoint(0, 0, 0)).z();
0120     }
0121 
0122     namespace {
0123       bool isMatchedByBlockElement(const reco::PFCandidate& pfCandidate1,
0124                                    const reco::PFCandidate& pfCandidate2,
0125                                    int minMatches1,
0126                                    int minMatches2,
0127                                    int maxUnmatchedBlockElements1plus2) {
0128         reco::PFCandidate::ElementsInBlocks blockElements1 = pfCandidate1.elementsInBlocks();
0129         int numBlocks1 = blockElements1.size();
0130         reco::PFCandidate::ElementsInBlocks blockElements2 = pfCandidate2.elementsInBlocks();
0131         int numBlocks2 = blockElements2.size();
0132         int numBlocks_matched = 0;
0133         for (reco::PFCandidate::ElementsInBlocks::const_iterator blockElement1 = blockElements1.begin();
0134              blockElement1 != blockElements1.end();
0135              ++blockElement1) {
0136           bool isMatched = false;
0137           for (reco::PFCandidate::ElementsInBlocks::const_iterator blockElement2 = blockElements2.begin();
0138                blockElement2 != blockElements2.end();
0139                ++blockElement2) {
0140             if (blockElement1->first.id() == blockElement2->first.id() &&
0141                 blockElement1->first.key() == blockElement2->first.key() &&
0142                 blockElement1->second == blockElement2->second) {
0143               isMatched = true;
0144             }
0145           }
0146           if (isMatched)
0147             ++numBlocks_matched;
0148         }
0149         assert(numBlocks_matched <= numBlocks1);
0150         assert(numBlocks_matched <= numBlocks2);
0151         if (numBlocks_matched >= minMatches1 && numBlocks_matched >= minMatches2 &&
0152             ((numBlocks1 - numBlocks_matched) + (numBlocks2 - numBlocks_matched)) <= maxUnmatchedBlockElements1plus2) {
0153           return true;
0154         } else {
0155           return false;
0156         }
0157       }
0158     }  // namespace
0159 
0160     PFRecoTauChargedHadronFromPFCandidatePlugin::return_type PFRecoTauChargedHadronFromPFCandidatePlugin::operator()(
0161         const reco::Jet& jet) const {
0162       if (verbosity_) {
0163         edm::LogPrint("TauChHadronFromPF") << "<PFRecoTauChargedHadronFromPFCandidatePlugin::operator()>:";
0164         edm::LogPrint("TauChHadronFromPF") << " pluginName = " << name();
0165       }
0166 
0167       ChargedHadronVector output;
0168 
0169       // Get the candidates passing our quality cuts
0170       qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0171       CandPtrs candsVector = qcuts_->filterCandRefs(pfCandidates(jet, inputParticleIds_));
0172 
0173       for (CandPtrs::iterator cand = candsVector.begin(); cand != candsVector.end(); ++cand) {
0174         if (verbosity_) {
0175           edm::LogPrint("TauChHadronFromPF")
0176               << "processing PFCandidate: Pt = " << (*cand)->pt() << ", eta = " << (*cand)->eta()
0177               << ", phi = " << (*cand)->phi() << " (pdgId = " << (*cand)->pdgId() << ", charge = " << (*cand)->charge()
0178               << ")";
0179         }
0180 
0181         PFRecoTauChargedHadron::PFRecoTauChargedHadronAlgorithm algo = PFRecoTauChargedHadron::kUndefined;
0182         if (std::abs((*cand)->charge()) > 0.5)
0183           algo = PFRecoTauChargedHadron::kChargedPFCandidate;
0184         else
0185           algo = PFRecoTauChargedHadron::kPFNeutralHadron;
0186         auto chargedHadron = std::make_unique<PFRecoTauChargedHadron>(**cand, algo);
0187 
0188         const reco::PFCandidate* pfCand = dynamic_cast<const reco::PFCandidate*>(&**cand);
0189         if (pfCand) {
0190           if (pfCand->trackRef().isNonnull())
0191             chargedHadron->track_ = edm::refToPtr(pfCand->trackRef());
0192           else if (pfCand->muonRef().isNonnull() && pfCand->muonRef()->innerTrack().isNonnull())
0193             chargedHadron->track_ = edm::refToPtr(pfCand->muonRef()->innerTrack());
0194           else if (pfCand->muonRef().isNonnull() && pfCand->muonRef()->globalTrack().isNonnull())
0195             chargedHadron->track_ = edm::refToPtr(pfCand->muonRef()->globalTrack());
0196           else if (pfCand->muonRef().isNonnull() && pfCand->muonRef()->outerTrack().isNonnull())
0197             chargedHadron->track_ = edm::refToPtr(pfCand->muonRef()->outerTrack());
0198           else if (pfCand->gsfTrackRef().isNonnull())
0199             chargedHadron->track_ = edm::refToPtr(pfCand->gsfTrackRef());
0200         }  // TauReco@MiniAOD: Tracks only available dynamically, so no possiblity to save ref here; checked by code downstream
0201 
0202         chargedHadron->positionAtECALEntrance_ = atECALEntrance(&**cand, bField_);
0203         chargedHadron->chargedPFCandidate_ = (*cand);
0204         chargedHadron->addDaughter(*cand);
0205 
0206         int pdgId = std::abs(chargedHadron->chargedPFCandidate_->pdgId());
0207 
0208         if (chargedHadron->pt() > minMergeChargedHadronPt_) {
0209           for (const auto& jetConstituent : jet.daughterPtrVector()) {
0210             // CV: take care of not double-counting energy in case "charged" PFCandidate is in fact a PFNeutralHadron
0211             if (jetConstituent == chargedHadron->chargedPFCandidate_)
0212               continue;
0213 
0214             int jetConstituentPdgId = std::abs(jetConstituent->pdgId());
0215             if (!(jetConstituentPdgId == 130 || jetConstituentPdgId == 22))
0216               continue;
0217 
0218             double dR = deltaR(atECALEntrance(jetConstituent.get(), bField_),
0219                                atECALEntrance(chargedHadron->chargedPFCandidate_.get(), bField_));
0220             double dRmerge = -1.;
0221             int minBlockElementMatches = 1000;
0222             int maxUnmatchedBlockElements = 0;
0223             double minMergeEt = 1.e+6;
0224             if (jetConstituentPdgId == 130) {
0225               if (pdgId == 211)
0226                 dRmerge = dRmergeNeutralHadronWrtChargedHadron_;
0227               else if (pdgId == 130)
0228                 dRmerge = dRmergeNeutralHadronWrtNeutralHadron_;
0229               else if (pdgId == 11)
0230                 dRmerge = dRmergeNeutralHadronWrtElectron_;
0231               else
0232                 dRmerge = dRmergeNeutralHadronWrtOther_;
0233               minBlockElementMatches = minBlockElementMatchesNeutralHadron_;
0234               maxUnmatchedBlockElements = maxUnmatchedBlockElementsNeutralHadron_;
0235               minMergeEt = minMergeNeutralHadronEt_;
0236             } else if (jetConstituentPdgId == 22) {
0237               if (pdgId == 211)
0238                 dRmerge = dRmergePhotonWrtChargedHadron_;
0239               else if (pdgId == 130)
0240                 dRmerge = dRmergePhotonWrtNeutralHadron_;
0241               else if (pdgId == 11)
0242                 dRmerge = dRmergePhotonWrtElectron_;
0243               else
0244                 dRmerge = dRmergePhotonWrtOther_;
0245               minBlockElementMatches = minBlockElementMatchesPhoton_;
0246               maxUnmatchedBlockElements = maxUnmatchedBlockElementsPhoton_;
0247               minMergeEt = minMergeGammaEt_;
0248             }
0249 
0250             if (jetConstituent->et() > minMergeEt) {
0251               if (dR < dRmerge) {
0252                 chargedHadron->neutralPFCandidates_.push_back(jetConstituent);
0253                 chargedHadron->addDaughter(jetConstituent);
0254               } else {
0255                 // TauReco@MiniAOD: No access to PF blocks at MiniAOD level, but the code below seems to have very minor impact
0256                 const reco::PFCandidate* pfJetConstituent =
0257                     dynamic_cast<const reco::PFCandidate*>(jetConstituent.get());
0258                 if (pfCand != nullptr && pfJetConstituent != nullptr) {
0259                   if (isMatchedByBlockElement(*pfJetConstituent,
0260                                               *pfCand,
0261                                               minBlockElementMatches,
0262                                               minBlockElementMatches,
0263                                               maxUnmatchedBlockElements)) {
0264                     chargedHadron->neutralPFCandidates_.push_back(jetConstituent);
0265                     chargedHadron->addDaughter(jetConstituent);
0266                   }
0267                 }
0268               }
0269             }
0270           }
0271         }
0272 
0273         setChargedHadronP4(*chargedHadron);
0274 
0275         if (verbosity_) {
0276           edm::LogPrint("TauChHadronFromPF") << *chargedHadron;
0277         }
0278         // Update the vertex
0279         if (chargedHadron->daughterPtr(0).isNonnull())
0280           chargedHadron->setVertex(chargedHadron->daughterPtr(0)->vertex());
0281         output.push_back(std::move(chargedHadron));
0282       }
0283 
0284       return output;
0285     }
0286 
0287   }  // namespace tau
0288 }  // namespace reco
0289 
0290 #include "FWCore/Framework/interface/MakerMacros.h"
0291 
0292 DEFINE_EDM_PLUGIN(PFRecoTauChargedHadronBuilderPluginFactory,
0293                   reco::tau::PFRecoTauChargedHadronFromPFCandidatePlugin,
0294                   "PFRecoTauChargedHadronFromPFCandidatePlugin");