RecoTau/plugins/PFRecoTauChargedHadronProducer.cc

0001 /*
0002  * PFRecoTauChargedHadronProducer
0003  *
0004  * Author: Christian Veelken, LLR
0005  *
0006  * Associates reconstructed ChargedHadrons to PFJets.  The ChargedHadrons are built using one
0007  * or more RecoTauBuilder plugins.  Any overlaps (ChargedHadrons sharing tracks)
0008  * are removed, with the best ChargedHadron candidates taken.  The 'best' are defined
0009  * via the input list of PFRecoTauChargedHadronQualityPlugins, which form a
0010  * lexicograpical ranking.
0011  *
0012  */
0013
0014 #include "FWCore/Framework/interface/stream/EDProducer.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021
0022 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0023 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0024
0025 #include "RecoTauTag/RecoTau/interface/PFRecoTauChargedHadronPlugins.h"
0026 #include "RecoTauTag/RecoTau/interface/RecoTauCleaningTools.h"
0027 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0028 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0029 #include "RecoTauTag/RecoTau/interface/pfRecoTauChargedHadronAuxFunctions.h"
0030
0031 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0032 #include "DataFormats/TauReco/interface/PFJetChargedHadronAssociation.h"
0033 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
0034 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0035 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0036 #include "DataFormats/TrackReco/interface/Track.h"
0037 #include "DataFormats/MuonReco/interface/Muon.h"
0038 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0039 #include "DataFormats/Common/interface/Association.h"
0040 #include "DataFormats/Math/interface/deltaR.h"
0041
0042 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0043
0044 #include <algorithm>
0045 #include <cmath>
0046 #include <functional>
0047 #include <list>
0048 #include <memory>
0049 #include <set>
0050 #include <string>
0051 #include <vector>
0052
0053 class PFRecoTauChargedHadronProducer : public edm::stream::EDProducer<> {
0054 public:
0055   typedef reco::tau::PFRecoTauChargedHadronBuilderPlugin Builder;
0056   typedef reco::tau::PFRecoTauChargedHadronQualityPlugin Ranker;
0057
0058   explicit PFRecoTauChargedHadronProducer(const edm::ParameterSet& cfg);
0059   ~PFRecoTauChargedHadronProducer() override {}
0060   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0061   template <typename T>
0062   void print(const T& chargedHadron);
0063
0064   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0065
0066 private:
0067   typedef std::vector<std::unique_ptr<Ranker>> RankerList;
0068   typedef Builder::return_type ChargedHadronVector;
0069   typedef std::list<std::unique_ptr<reco::PFRecoTauChargedHadron>> ChargedHadronList;
0070
0071   typedef reco::tau::RecoTauLexicographicalRanking<RankerList, reco::PFRecoTauChargedHadron> ChargedHadronPredicate;
0072
0073   std::string moduleLabel_;
0074
0075   // input jet collection
0076   edm::InputTag srcJets_;
0077   edm::EDGetTokenT<reco::JetView> Jets_token;
0078   double minJetPt_;
0079   double maxJetAbsEta_;
0080
0081   // plugins for building and ranking ChargedHadron candidates
0082   std::vector<std::unique_ptr<Builder>> builders_;
0083   RankerList rankers_;
0084
0085   std::unique_ptr<ChargedHadronPredicate> predicate_;
0086
0087   // output selector
0088   std::unique_ptr<StringCutObjectSelector<reco::PFRecoTauChargedHadron>> outputSelector_;
0089
0090   // flag to enable/disable debug print-out
0091   int verbosity_;
0092 };
0093
0094 PFRecoTauChargedHadronProducer::PFRecoTauChargedHadronProducer(const edm::ParameterSet& cfg)
0095     : moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
0096   srcJets_ = cfg.getParameter<edm::InputTag>("jetSrc");
0097   Jets_token = consumes<reco::JetView>(srcJets_);
0098   minJetPt_ = cfg.getParameter<double>("minJetPt");
0099   maxJetAbsEta_ = cfg.getParameter<double>("maxJetAbsEta");
0100   verbosity_ = cfg.getParameter<int>("verbosity");
0101
0102   // get set of ChargedHadron builder plugins
0103   edm::VParameterSet psets_builders = cfg.getParameter<edm::VParameterSet>("builders");
0104   for (edm::VParameterSet::const_iterator pset = psets_builders.begin(); pset != psets_builders.end(); ++pset) {
0105     std::string pluginType = pset->getParameter<std::string>("plugin");
0106     edm::ParameterSet pset_modified = (*pset);
0107     pset_modified.addParameter<int>("verbosity", verbosity_);
0108     builders_.emplace_back(
0109         PFRecoTauChargedHadronBuilderPluginFactory::get()->create(pluginType, pset_modified, consumesCollector()));
0110   }
0111
0112   // get set of plugins for ranking ChargedHadrons in quality
0113   edm::VParameterSet psets_rankers = cfg.getParameter<edm::VParameterSet>("ranking");
0114   for (edm::VParameterSet::const_iterator pset = psets_rankers.begin(); pset != psets_rankers.end(); ++pset) {
0115     std::string pluginType = pset->getParameter<std::string>("plugin");
0116     edm::ParameterSet pset_modified = (*pset);
0117     pset_modified.addParameter<int>("verbosity", verbosity_);
0118     rankers_.emplace_back(PFRecoTauChargedHadronQualityPluginFactory::get()->create(pluginType, pset_modified));
0119   }
0120
0121   // build the sorting predicate
0122   predicate_ = std::make_unique<ChargedHadronPredicate>(rankers_);
0123
0124   // check if we want to apply a final output selection
0125   std::string selection = cfg.getParameter<std::string>("outputSelection");
0126   if (!selection.empty()) {
0127     outputSelector_ = std::make_unique<StringCutObjectSelector<reco::PFRecoTauChargedHadron>>(selection);
0128   }
0129
0130   produces<reco::PFJetChargedHadronAssociation>();
0131 }
0132
0133 void PFRecoTauChargedHadronProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0134   if (verbosity_) {
0135     edm::LogPrint("PFRecoTauChHProducer") << "<PFRecoTauChargedHadronProducer::produce>:";
0136     edm::LogPrint("PFRecoTauChHProducer") << " moduleLabel = " << moduleLabel_;
0137   }
0138
0139   // give each of our plugins a chance at doing something with the edm::Event
0140   for (auto& builder : builders_) {
0141     builder->setup(evt, es);
0142   }
0143
0144   // get a view of our jets via the base candidates
0145   edm::Handle<reco::JetView> jets;
0146   evt.getByToken(Jets_token, jets);
0147
0148   // convert the view to a RefVector of actual PFJets
0149   edm::RefToBaseVector<reco::Jet> pfJets;
0150   size_t nElements = jets->size();
0151   for (size_t i = 0; i < nElements; ++i) {
0152     pfJets.push_back(jets->refAt(i));
0153   }
0154
0155   // make our association
0156   std::unique_ptr<reco::PFJetChargedHadronAssociation> pfJetChargedHadronAssociations;
0157
0158   if (!pfJets.empty()) {
0159     pfJetChargedHadronAssociations = std::make_unique<reco::PFJetChargedHadronAssociation>(reco::JetRefBaseProd(jets));
0160   } else {
0161     pfJetChargedHadronAssociations = std::make_unique<reco::PFJetChargedHadronAssociation>();
0162   }
0163
0164   // loop over our jets
0165   for (const auto& pfJet : pfJets) {
0166     if (pfJet->pt() - minJetPt_ < 1e-5)
0167       continue;
0168     if (std::abs(pfJet->eta()) - maxJetAbsEta_ > -1e-5)
0169       continue;
0170
0171     // build global list of ChargedHadron candidates for each jet
0172     ChargedHadronList uncleanedChargedHadrons;
0173
0174     // merge candidates reconstructed by all desired algorithm plugins
0175     for (auto const& builder : builders_) {
0176       try {
0177         ChargedHadronVector result((*builder)(*pfJet));
0178         if (verbosity_) {
0179           edm::LogPrint("PFRecoTauChHProducer") << "result of builder = " << builder->name() << ":";
0180           for (auto const& had : result) {
0181             print(*had);
0182           }
0183         }
0184         std::move(result.begin(), result.end(), std::back_inserter(uncleanedChargedHadrons));
0185       } catch (cms::Exception& exception) {
0186         edm::LogError("BuilderPluginException")
0187             << "Exception caught in builder plugin " << builder->name() << ", rethrowing" << std::endl;
0188         throw exception;
0189       }
0190     }
0191
0192     // rank the candidates according to our quality plugins
0193     uncleanedChargedHadrons.sort([this](const auto& a, const auto& b) { return (*predicate_)(*a, *b); });
0194
0195     // define collection of cleaned ChargedHadrons;
0196     std::vector<reco::PFRecoTauChargedHadron> cleanedChargedHadrons;
0197
0198     // keep track of neutral PFCandidates, charged PFCandidates and tracks "used" by ChargedHadron candidates in the clean collection
0199     typedef std::pair<double, double> etaPhiPair;
0200     std::list<etaPhiPair> tracksInCleanCollection;
0201     std::set<reco::CandidatePtr> neutralPFCandsInCleanCollection;
0202
0203     while (!uncleanedChargedHadrons.empty()) {
0204       // get next best ChargedHadron candidate
0205       std::unique_ptr<reco::PFRecoTauChargedHadron> nextChargedHadron = std::move(uncleanedChargedHadrons.front());
0206       uncleanedChargedHadrons.pop_front();
0207       if (verbosity_) {
0208         edm::LogPrint("PFRecoTauChHProducer") << "processing nextChargedHadron:";
0209         edm::LogPrint("PFRecoTauChHProducer") << (*nextChargedHadron);
0210       }
0211
0212       // discard candidates which fail final output selection
0213       if (!(*outputSelector_)(*nextChargedHadron))
0214         continue;
0215
0216       const reco::Track* track = reco::tau::getTrackFromChargedHadron(*nextChargedHadron);
0217
0218       // discard candidate in case its track is "used" by any ChargedHadron in the clean collection
0219       bool isTrack_overlap = false;
0220       if (track) {
0221         double track_eta = track->eta();
0222         double track_phi = track->phi();
0223         for (std::list<etaPhiPair>::const_iterator trackInCleanCollection = tracksInCleanCollection.begin();
0224              trackInCleanCollection != tracksInCleanCollection.end();
0225              ++trackInCleanCollection) {
0226           double dR = deltaR(track_eta, track_phi, trackInCleanCollection->first, trackInCleanCollection->second);
0227           if (dR < 1.e-4)
0228             isTrack_overlap = true;
0229         }
0230       }
0231       if (verbosity_) {
0232         edm::LogPrint("PFRecoTauChHProducer") << "isTrack_overlap = " << isTrack_overlap;
0233       }
0234       if (isTrack_overlap)
0235         continue;
0236
0237       // discard ChargedHadron candidates without track in case they are close to neutral PFCandidates "used" by ChargedHadron candidates in the clean collection
0238       bool isNeutralPFCand_overlap = false;
0239       if (nextChargedHadron->algoIs(reco::PFRecoTauChargedHadron::kPFNeutralHadron)) {
0240         for (std::set<reco::CandidatePtr>::const_iterator neutralPFCandInCleanCollection =
0241                  neutralPFCandsInCleanCollection.begin();
0242              neutralPFCandInCleanCollection != neutralPFCandsInCleanCollection.end();
0243              ++neutralPFCandInCleanCollection) {
0244           if ((*neutralPFCandInCleanCollection) == nextChargedHadron->getChargedPFCandidate())
0245             isNeutralPFCand_overlap = true;
0246         }
0247       }
0248       if (verbosity_) {
0249         edm::LogPrint("PFRecoTauChHProducer") << "isNeutralPFCand_overlap = " << isNeutralPFCand_overlap;
0250       }
0251       if (isNeutralPFCand_overlap)
0252         continue;
0253
0254       // find neutral PFCandidates that are not "used" by any ChargedHadron in the clean collection
0255       std::vector<reco::CandidatePtr> uniqueNeutralPFCands;
0256       std::set_difference(nextChargedHadron->getNeutralPFCandidates().begin(),
0257                           nextChargedHadron->getNeutralPFCandidates().end(),
0258                           neutralPFCandsInCleanCollection.begin(),
0259                           neutralPFCandsInCleanCollection.end(),
0260                           std::back_inserter(uniqueNeutralPFCands));
0261
0262       if (uniqueNeutralPFCands.size() ==
0263           nextChargedHadron->getNeutralPFCandidates()
0264               .size()) {  // all neutral PFCandidates are unique, add ChargedHadron candidate to clean collection
0265         if (track)
0266           tracksInCleanCollection.push_back(std::make_pair(track->eta(), track->phi()));
0267         neutralPFCandsInCleanCollection.insert(nextChargedHadron->getNeutralPFCandidates().begin(),
0268                                                nextChargedHadron->getNeutralPFCandidates().end());
0269         if (verbosity_) {
0270           edm::LogPrint("PFRecoTauChHProducer") << "--> adding nextChargedHadron to output collection.";
0271         }
0272         cleanedChargedHadrons.push_back(*nextChargedHadron);
0273       } else {  // remove overlapping neutral PFCandidates, reevaluate ranking criterion and process ChargedHadron candidate again
0274         nextChargedHadron->neutralPFCandidates_.clear();
0275         for (auto const& neutralPFCand : uniqueNeutralPFCands) {
0276           nextChargedHadron->neutralPFCandidates_.push_back(neutralPFCand);
0277         }
0278         // update ChargedHadron four-momentum
0279         reco::tau::setChargedHadronP4(*nextChargedHadron);
0280         // reinsert ChargedHadron candidate into list of uncleaned candidates,
0281         // at position according to new rank
0282         auto insertionPoint = std::lower_bound(uncleanedChargedHadrons.begin(),
0283                                                uncleanedChargedHadrons.end(),
0284                                                *nextChargedHadron,
0285                                                [this](const auto& a, const auto& b) { return (*predicate_)(*a, b); });
0286         if (verbosity_) {
0287           edm::LogPrint("PFRecoTauChHProducer") << "--> removing non-unique neutral PFCandidates and reinserting "
0288                                                    "nextChargedHadron in uncleaned collection.";
0289         }
0290         uncleanedChargedHadrons.insert(insertionPoint, std::move(nextChargedHadron));
0291       }
0292     }
0293
0294     if (verbosity_) {
0295       for (auto const& had : cleanedChargedHadrons) {
0296         print(had);
0297       }
0298     }
0299
0300     // add ChargedHadron-to-jet association
0301     pfJetChargedHadronAssociations->setValue(pfJet.key(), cleanedChargedHadrons);
0302   }
0303
0304   evt.put(std::move(pfJetChargedHadronAssociations));
0305 }
0306
0307 template <typename T>
0308 void PFRecoTauChargedHadronProducer::print(const T& chargedHadron) {
0309   edm::LogPrint("PFRecoTauChHProducer") << chargedHadron;
0310   edm::LogPrint("PFRecoTauChHProducer") << "Rankers:";
0311   for (auto const& ranker : rankers_) {
0312     constexpr unsigned width = 25;
0313     edm::LogPrint("PFRecoTauChHProducer")
0314         << " " << std::setiosflags(std::ios::left) << std::setw(width) << ranker->name() << " "
0315         << std::resetiosflags(std::ios::left) << std::setprecision(3) << (*ranker)(chargedHadron) << std::endl;
0316   }
0317 }
0318
0319 void PFRecoTauChargedHadronProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0320   // ak4PFJetsRecoTauChargedHadrons
0321   edm::ParameterSetDescription desc;
0322   {
0323     edm::ParameterSetDescription desc_ranking;
0324     desc_ranking.add<std::string>("selectionPassFunction", "-pt");
0325     desc_ranking.add<double>("selectionFailValue", 1000.0);
0326     desc_ranking.add<std::string>("selection", "algoIs(\"kChargedPFCandidate\")");
0327     desc_ranking.add<std::string>("name", "ChargedPFCandidate");
0328     desc_ranking.add<std::string>("plugin", "PFRecoTauChargedHadronStringQuality");
0329
0330     edm::ParameterSet pset_ranking;
0331     pset_ranking.addParameter<std::string>("selectionPassFunction", "-pt");
0332     pset_ranking.addParameter<double>("selectionFailValue", 1000.0);
0333     pset_ranking.addParameter<std::string>("selection", "algoIs(\"kChargedPFCandidate\")");
0334     pset_ranking.addParameter<std::string>("name", "ChargedPFCandidate");
0335     pset_ranking.addParameter<std::string>("plugin", "PFRecoTauChargedHadronStringQuality");
0336     std::vector<edm::ParameterSet> vpsd_ranking;
0337     vpsd_ranking.push_back(pset_ranking);
0338
0339     desc.addVPSet("ranking", desc_ranking, vpsd_ranking);
0340   }
0341
0342   desc.add<int>("verbosity", 0);
0343   desc.add<double>("maxJetAbsEta", 2.5);
0344   desc.add<std::string>("outputSelection", "pt > 0.5");
0345   desc.add<double>("minJetPt", 14.0);
0346   desc.add<edm::InputTag>("jetSrc", edm::InputTag("ak4PFJets"));
0347
0348   {
0349     edm::ParameterSetDescription desc_builders;
0350     desc_builders.add<double>("minMergeChargedHadronPt");
0351     desc_builders.add<std::string>("name");
0352     desc_builders.add<std::string>("plugin");
0353     desc_builders.addOptional<double>("dRcone");
0354     desc_builders.addOptional<bool>("dRconeLimitedToJetArea");
0355     desc_builders.addOptional<double>("dRmergeNeutralHadron");
0356     desc_builders.addOptional<double>("dRmergePhoton");
0357     desc_builders.addOptional<edm::InputTag>("srcTracks");
0358
0359     edm::ParameterSetDescription desc_qualityCuts;
0360     reco::tau::RecoTauQualityCuts::fillDescriptions(desc_qualityCuts);
0361     desc_builders.add<edm::ParameterSetDescription>("qualityCuts", desc_qualityCuts);
0362
0363     desc_builders.add<double>("minMergeGammaEt");
0364     desc_builders.add<int>("verbosity", 0);
0365     desc_builders.add<double>("minMergeNeutralHadronEt");
0366
0367     desc_builders.addOptional<double>("dRmergePhotonWrtChargedHadron");
0368     desc_builders.addOptional<double>("dRmergePhotonWrtNeutralHadron");
0369     desc_builders.addOptional<int>("maxUnmatchedBlockElementsNeutralHadron");
0370     desc_builders.addOptional<double>("dRmergePhotonWrtElectron");
0371     desc_builders.addOptional<std::vector<int>>("chargedHadronCandidatesParticleIds");
0372     desc_builders.addOptional<int>("minBlockElementMatchesPhoton");
0373     desc_builders.addOptional<double>("dRmergeNeutralHadronWrtNeutralHadron");
0374     desc_builders.addOptional<int>("maxUnmatchedBlockElementsPhoton");
0375     desc_builders.addOptional<double>("dRmergeNeutralHadronWrtOther");
0376     desc_builders.addOptional<double>("dRmergeNeutralHadronWrtElectron");
0377     desc_builders.addOptional<int>("minBlockElementMatchesNeutralHadron");
0378     desc_builders.addOptional<double>("dRmergePhotonWrtOther");
0379     desc_builders.addOptional<double>("dRmergeNeutralHadronWrtChargedHadron");
0380
0381     edm::ParameterSet pset_builders;
0382     pset_builders.addParameter<std::string>("name", "");
0383     pset_builders.addParameter<std::string>("plugin", "");
0384     edm::ParameterSet qualityCuts;
0385     pset_builders.addParameter<edm::ParameterSet>("qualityCuts", qualityCuts);
0386     pset_builders.addParameter<int>("verbosity", 0);
0387     std::vector<edm::ParameterSet> vpsd_builders;
0388     vpsd_builders.push_back(pset_builders);
0389
0390     desc.addVPSet("builders", desc_builders, vpsd_builders);
0391   }
0392
0393   descriptions.add("pfRecoTauChargedHadronProducer", desc);
0394 }
0395
0396 #include "FWCore/Framework/interface/MakerMacros.h"
0397
0398 DEFINE_FWK_MODULE(PFRecoTauChargedHadronProducer);