RecoTau/plugins/RecoTauCleaner.cc

0001 /*
0002  * RecoTauCleaner
0003  *
0004  * Author: Evan K. Friis, UC Davis
0005  *
0006  * Given a input collection of PFTaus, produces an output collection of PFTaus
0007  * such that no two PFTaus come from the same PFJet.  If multiple taus in the
0008  * collection come from the same PFJet, (dirty) they are ordered according to a
0009  * list of cleaners.  Each cleaner is a RecoTauCleanerPlugin, and returns a
0010  * double corresponding to the 'quality' for a given tau - an example would be
0011  * the level of isolation.  The set of dirty taus is then ranked
0012  * lexicographically by these cleaners, and the best one is placed in the
0013  * output collection.
0014  *
0015  */
0016
0017 #include <algorithm>
0018 #include <memory>
0019
0020 #include "FWCore/Framework/interface/stream/EDProducer.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024
0025 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0026 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0027
0028 #include "RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h"
0029 #include "RecoTauTag/RecoTau/interface/RecoTauCleaningTools.h"
0030 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0031
0032 #include "DataFormats/TauReco/interface/PFTau.h"
0033 #include "DataFormats/TauReco/interface/PFTauFwd.h"
0034 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
0035 #include "DataFormats/TauReco/interface/RecoTauPiZero.h"
0036 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0037 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0038
0039 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0040
0041 template <typename Prod>
0042 class RecoTauCleanerImpl : public edm::stream::EDProducer<> {
0043   typedef reco::tau::RecoTauCleanerPlugin Cleaner;
0044   struct CleanerEntryType {
0045     std::shared_ptr<Cleaner> plugin_;
0046     float tolerance_;
0047   };
0048   typedef std::vector<std::unique_ptr<CleanerEntryType>> CleanerList;
0049   // Define our output type - i.e. reco::PFTau OR reco::PFTauRef
0050   typedef typename Prod::value_type output_type;
0051
0052   // Predicate that determines if two taus 'overlap' i.e. share a base PFJet
0053   class RemoveDuplicateJets {
0054   public:
0055     bool operator()(const reco::PFTauRef& a, const reco::PFTauRef& b) const { return (a->jetRef() == b->jetRef()); }
0056   };
0057
0058 public:
0059   explicit RecoTauCleanerImpl(const edm::ParameterSet& pset);
0060   ~RecoTauCleanerImpl() override;
0061   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0062
0063   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0064
0065 private:
0066   edm::InputTag tauSrc_;
0067   CleanerList cleaners_;
0068   // Optional selection on the output of the taus
0069   std::unique_ptr<const StringCutObjectSelector<reco::PFTau>> outputSelector_;
0070   edm::EDGetTokenT<reco::PFTauCollection> tau_token;
0071   int verbosity_;
0072 };
0073
0074 template <typename Prod>
0075 RecoTauCleanerImpl<Prod>::RecoTauCleanerImpl(const edm::ParameterSet& pset) {
0076   tauSrc_ = pset.getParameter<edm::InputTag>("src");
0077   tau_token = consumes<reco::PFTauCollection>(tauSrc_);
0078   // Build our list of quality plugins
0079   typedef std::vector<edm::ParameterSet> VPSet;
0080   // Get each of our tau builders
0081   const VPSet& cleaners = pset.getParameter<VPSet>("cleaners");
0082   for (VPSet::const_iterator cleanerPSet = cleaners.begin(); cleanerPSet != cleaners.end(); ++cleanerPSet) {
0083     auto cleanerEntry = std::make_unique<CleanerEntryType>();
0084     // Get plugin name
0085     const std::string& pluginType = cleanerPSet->getParameter<std::string>("plugin");
0086     // Build the plugin
0087     cleanerEntry->plugin_ = RecoTauCleanerPluginFactory::get()->create(pluginType, *cleanerPSet, consumesCollector());
0088     cleanerEntry->tolerance_ = cleanerPSet->getParameter<double>("tolerance");
0089     cleaners_.emplace_back(std::move(cleanerEntry));
0090   }
0091
0092   // Check if we want to apply a final output selection
0093   std::string selection = pset.getParameter<std::string>("outputSelection");
0094   if (!selection.empty()) {
0095     outputSelector_ = std::make_unique<StringCutObjectSelector<reco::PFTau>>(selection);
0096   }
0097
0098   // Enable/disable debug output
0099   verbosity_ = pset.getParameter<int>("verbosity");
0100
0101   // Build the predicate that ranks our taus.
0102   produces<Prod>();
0103 }
0104
0105 template <typename Prod>
0106 RecoTauCleanerImpl<Prod>::~RecoTauCleanerImpl() {}
0107
0108 namespace {
0109   // Template to convert a ref to desired output type
0110   template <typename T>
0111   const T convert(const reco::PFTauRef& tau);
0112
0113   //template<> const edm::RefToBase<reco::PFTau>
0114   //convert<edm::RefToBase<reco::PFTau> >(const reco::PFTauRef &tau) {
0115   // return edm::RefToBase<reco::PFTau>(tau);
0116   //}
0117
0118   template <>
0119   const reco::PFTauRef convert<reco::PFTauRef>(const reco::PFTauRef& tau) {
0120     return tau;
0121   }
0122
0123   template <>
0124   const reco::PFTau convert<reco::PFTau>(const reco::PFTauRef& tau) {
0125     return *tau;
0126   }
0127 }  // namespace
0128
0129 namespace {
0130   template <typename T>
0131   std::string format_vT(const std::vector<T>& vT) {
0132     std::ostringstream os;
0133     os << "{ ";
0134     unsigned numEntries = vT.size();
0135     for (unsigned iEntry = 0; iEntry < numEntries; ++iEntry) {
0136       os << vT[iEntry];
0137       if (iEntry < (numEntries - 1))
0138         os << ", ";
0139     }
0140     os << " }";
0141     return os.str();
0142   }
0143
0144   struct PFTauRankType {
0145     PFTauRankType(const reco::PFTauRef& tauRef) : idx_(tauRef.key()), tauRef_(tauRef) {}
0146     ~PFTauRankType() {}
0147     void print(const std::string& label) const {
0148       std::cout << label << " (" << tauRef_.id() << ":" << tauRef_.key() << ", idx = " << idx_ << "):";
0149       assert(tauRef_.key() == idx_);
0150       std::cout << " Pt = " << tauRef_->pt() << ", eta = " << tauRef_->eta() << ", phi = " << tauRef_->phi()
0151                 << ", mass = " << tauRef_->mass() << " (decayMode = " << tauRef_->decayMode() << ")";
0152       std::cout << std::endl;
0153       std::cout << "associated jet:";
0154       if (tauRef_->jetRef().isNonnull()) {
0155         std::cout << " Pt = " << tauRef_->jetRef()->pt() << ", eta = " << tauRef_->jetRef()->eta()
0156                   << ", phi = " << tauRef_->jetRef()->phi() << ", mass = " << tauRef_->jetRef()->mass()
0157                   << ", area = " << tauRef_->jetRef()->jetArea();
0158       } else
0159         std::cout << " N/A";
0160       std::cout << std::endl;
0161       const std::vector<reco::PFRecoTauChargedHadron>& signalTauChargedHadronCandidates =
0162           tauRef_->signalTauChargedHadronCandidates();
0163       size_t numChargedHadrons = signalTauChargedHadronCandidates.size();
0164       for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
0165         const reco::PFRecoTauChargedHadron& chargedHadron = signalTauChargedHadronCandidates.at(iChargedHadron);
0166         std::cout << " chargedHadron #" << iChargedHadron << ":" << std::endl;
0167         chargedHadron.print(std::cout);
0168       }
0169       const std::vector<reco::RecoTauPiZero>& signalPiZeroCandidates = tauRef_->signalPiZeroCandidates();
0170       size_t numPiZeros = signalPiZeroCandidates.size();
0171       std::cout << "signalPiZeroCandidates = " << numPiZeros << std::endl;
0172       for (size_t iPiZero = 0; iPiZero < numPiZeros; ++iPiZero) {
0173         const reco::RecoTauPiZero& piZero = signalPiZeroCandidates.at(iPiZero);
0174         std::cout << " piZero #" << iPiZero << ": Pt = " << piZero.pt() << ", eta = " << piZero.eta()
0175                   << ", phi = " << piZero.phi() << ", mass = " << piZero.mass() << std::endl;
0176       }
0177       const auto& isolationCands = tauRef_->isolationCands();
0178       size_t numCands = isolationCands.size();
0179       std::cout << "isolationCands = " << numCands << std::endl;
0180       for (size_t iCand = 0; iCand < numCands; ++iCand) {
0181         const auto& cand = isolationCands.at(iCand);
0182         std::cout << " pfCand #" << iCand << " (" << cand.id() << ":" << cand.key() << "):"
0183                   << " Pt = " << cand->pt() << ", eta = " << cand->eta() << ", phi = " << cand->phi() << std::endl;
0184       }
0185       std::cout << " ranks = " << format_vT(ranks_) << std::endl;
0186       std::cout << " tolerances = " << format_vT(tolerances_) << std::endl;
0187     }
0188     size_t idx_;
0189     reco::PFTauRef tauRef_;
0190     size_t N_;
0191     std::vector<float> ranks_;
0192     std::vector<float> tolerances_;
0193   };
0194
0195   bool isHigherRank(const PFTauRankType* tau1, const PFTauRankType* tau2) {
0196     //std::cout << "<isHigherRank>:" << std::endl;
0197     //std::cout << "tau1 @ " << tau1;
0198     //tau1->print("");
0199     //std::cout << "tau2 @ " << tau2;
0200     //tau2->print("");
0201     assert(tau1->N_ == tau1->ranks_.size());
0202     assert(tau1->N_ == tau2->ranks_.size());
0203     assert(tau1->N_ == tau1->tolerances_.size());
0204     for (size_t i = 0; i < tau1->N_; ++i) {
0205       const float& val1 = tau1->ranks_[i];
0206       const float& val2 = tau2->ranks_[i];
0207       double av = 0.5 * (val1 + val2);
0208       double thresh = av * tau1->tolerances_[i];
0209       if (val1 < (val2 - thresh))
0210         return true;
0211       else if (val2 < (val1 - thresh))
0212         return false;
0213     }
0214     return true;
0215   }
0216 }  // namespace
0217
0218 template <typename Prod>
0219 void RecoTauCleanerImpl<Prod>::produce(edm::Event& evt, const edm::EventSetup& es) {
0220   if (verbosity_) {
0221     std::cout << "<RecoTauCleanerImpl::produce>:" << std::endl;
0222   }
0223
0224   // Update all our cleaners with the event info if they need it
0225   for (typename CleanerList::iterator cleaner = cleaners_.begin(); cleaner != cleaners_.end(); ++cleaner) {
0226     (*cleaner)->plugin_->setup(evt, es);
0227   }
0228
0229   // Get the input collection of all taus. Some are from the same PFJet. We must clean them.
0230   edm::Handle<reco::PFTauCollection> inputTaus;
0231   evt.getByToken(tau_token, inputTaus);
0232
0233   // Sort the input tau refs according to our predicate
0234   std::list<PFTauRankType*> rankedTaus;
0235   size_t N = inputTaus->size();
0236   for (size_t idx = 0; idx < N; ++idx) {
0237     reco::PFTauRef inputRef(inputTaus, idx);
0238     PFTauRankType* rankedTau = new PFTauRankType(inputRef);
0239     rankedTau->N_ = cleaners_.size();
0240     rankedTau->ranks_.reserve(rankedTau->N_);
0241     rankedTau->tolerances_.reserve(rankedTau->N_);
0242     for (typename CleanerList::const_iterator cleaner = cleaners_.begin(); cleaner != cleaners_.end(); ++cleaner) {
0243       rankedTau->ranks_.push_back((*(*cleaner)->plugin_)(inputRef));
0244       rankedTau->tolerances_.push_back((*cleaner)->tolerance_);
0245     }
0246     if (verbosity_) {
0247       std::ostringstream os;
0248       os << "rankedTau #" << idx;
0249       rankedTau->print(os.str());
0250     }
0251     rankedTaus.push_back(rankedTau);
0252   }
0253   rankedTaus.sort(isHigherRank);
0254
0255   // Make an STL algorithm friendly vector of refs
0256   typedef std::vector<reco::PFTauRef> PFTauRefs;
0257   PFTauRefs dirty(inputTaus->size());
0258   size_t idx_sorted = 0;
0259   for (std::list<PFTauRankType*>::const_iterator rankedTau = rankedTaus.begin(); rankedTau != rankedTaus.end();
0260        ++rankedTau) {
0261     dirty[idx_sorted] = (*rankedTau)->tauRef_;
0262     if (verbosity_) {
0263       std::cout << "dirty[" << idx_sorted << "] = " << dirty[idx_sorted].id() << ":" << dirty[idx_sorted].key()
0264                 << std::endl;
0265     }
0266     delete (*rankedTau);
0267     ++idx_sorted;
0268   }
0269
0270   // Clean the taus, ensuring that only one tau per jet is produced
0271   PFTauRefs cleanTaus = reco::tau::cleanOverlaps<PFTauRefs, RemoveDuplicateJets>(dirty);
0272
0273   // create output collection
0274   auto output = std::make_unique<Prod>();
0275   //output->reserve(cleanTaus.size());
0276
0277   // Copy clean refs into output
0278   for (PFTauRefs::const_iterator tau = cleanTaus.begin(); tau != cleanTaus.end(); ++tau) {
0279     // If we are applying an output selection, check if it passes
0280     bool selected = true;
0281     if (outputSelector_.get() && !(*outputSelector_)(**tau)) {
0282       selected = false;
0283     }
0284     if (selected) {
0285       output->push_back(convert<output_type>(*tau));
0286     }
0287   }
0288   evt.put(std::move(output));
0289 }
0290
0291 typedef RecoTauCleanerImpl<reco::PFTauCollection> RecoTauCleaner;
0292 typedef RecoTauCleanerImpl<reco::PFTauRefVector> RecoTauRefCleaner;
0293
0294 template <>
0295 void RecoTauCleanerImpl<reco::PFTauCollection>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0296   // RecoTauCleaner
0297   edm::ParameterSetDescription desc;
0298   desc.add<std::string>("outputSelection", "");
0299   {
0300     // this description is the validator for all PSets in the cleaners VPSet passed to the plugin from the python configuration
0301     edm::ParameterSetDescription vps_description_for_cleaners;
0302     // the common parameters for all cleaners
0303     // no default value is provided -- the user has to provide the values for these parameters
0304     vps_description_for_cleaners.add<std::string>("plugin");
0305     vps_description_for_cleaners.add<double>("tolerance", 0);
0306     vps_description_for_cleaners.add<std::string>("name");
0307
0308     // the following parameters are not common for all cleaners, they are needed for few PSets, therefore they are added as optional
0309     //     these optional parameters are used in the default cleaners vector
0310     vps_description_for_cleaners.addOptional<int>("passForCharge");
0311     vps_description_for_cleaners.addOptional<double>("selectionFailValue");
0312     vps_description_for_cleaners.addOptional<std::vector<unsigned int>>("nprongs");
0313     vps_description_for_cleaners.addOptional<edm::InputTag>("src");
0314     vps_description_for_cleaners.addOptional<double>("minTrackPt");
0315     vps_description_for_cleaners.addOptional<std::string>("selection");
0316     vps_description_for_cleaners.addOptional<std::string>("selectionPassFunction");
0317     //     more PSets for cleaners can be found in
0318     //     RecoTauTag/RecoTau/python/RecoTauCleanerPlugins.py
0319     //     however, at this moment (2018-11-09) they do not have any new optional parameters
0320
0321     // the cleaner defaults, as in RecoTauTag/RecoTau/python/RecoTauCleaner_cfi.py
0322     std::vector<edm::ParameterSet> default_cleaners;
0323     default_cleaners.reserve(7);
0324     {
0325       edm::ParameterSet cleaner_Charge;
0326       cleaner_Charge.addParameter<std::string>("name", "Charge");
0327       cleaner_Charge.addParameter<std::string>("plugin", "RecoTauChargeCleanerPlugin");
0328       cleaner_Charge.addParameter<int>("passForCharge", 1);
0329       cleaner_Charge.addParameter<double>("selectionFailValue", 0);
0330       cleaner_Charge.addParameter<std::vector<unsigned int>>("nprongs",
0331                                                              {
0332                                                                  1,
0333                                                                  3,
0334                                                              });
0335       cleaner_Charge.addParameter<double>("tolerance", 0);
0336       default_cleaners.push_back(cleaner_Charge);
0337     }
0338     {
0339       edm::ParameterSet temp2;
0340       temp2.addParameter<std::string>("name", "HPS_Select");
0341       temp2.addParameter<std::string>("plugin", "RecoTauDiscriminantCleanerPlugin");
0342       temp2.addParameter<edm::InputTag>("src", edm::InputTag("hpsSelectionDiscriminator"));
0343       temp2.addParameter<double>("tolerance", 0);
0344       default_cleaners.push_back(temp2);
0345     }
0346     {
0347       edm::ParameterSet temp2;
0348       temp2.addParameter<std::string>("name", "killSoftTwoProngTaus");
0349       temp2.addParameter<std::string>("plugin", "RecoTauSoftTwoProngTausCleanerPlugin");
0350       temp2.addParameter<double>("minTrackPt", 5.0);
0351       temp2.addParameter<double>("tolerance", 0);
0352       default_cleaners.push_back(temp2);
0353     }
0354     {
0355       edm::ParameterSet temp2;
0356       temp2.addParameter<std::string>("name", "ChargedHadronMultiplicity");
0357       temp2.addParameter<std::string>("plugin", "RecoTauChargedHadronMultiplicityCleanerPlugin");
0358       temp2.addParameter<double>("tolerance", 0);
0359       default_cleaners.push_back(temp2);
0360     }
0361     {
0362       edm::ParameterSet temp2;
0363       temp2.addParameter<std::string>("name", "Pt");
0364       temp2.addParameter<std::string>("plugin", "RecoTauStringCleanerPlugin");
0365       temp2.addParameter<std::string>("selectionPassFunction", "-pt()");
0366       temp2.addParameter<std::string>("selection", "leadPFCand().isNonnull()");
0367       temp2.addParameter<double>("selectionFailValue", 1000.0);
0368       temp2.addParameter<double>("tolerance", 0.01);
0369       default_cleaners.push_back(temp2);
0370     }
0371     {
0372       edm::ParameterSet temp2;
0373       temp2.addParameter<std::string>("name", "StripMultiplicity");
0374       temp2.addParameter<std::string>("plugin", "RecoTauStringCleanerPlugin");
0375       temp2.addParameter<std::string>("selectionPassFunction", "-signalPiZeroCandidates().size()");
0376       temp2.addParameter<std::string>("selection", "leadPFCand().isNonnull()");
0377       temp2.addParameter<double>("selectionFailValue", 1000.0);
0378       temp2.addParameter<double>("tolerance", 0);
0379       default_cleaners.push_back(temp2);
0380     }
0381     {
0382       edm::ParameterSet temp2;
0383       temp2.addParameter<std::string>("name", "CombinedIsolation");
0384       temp2.addParameter<std::string>("plugin", "RecoTauStringCleanerPlugin");
0385       temp2.addParameter<std::string>("selectionPassFunction",
0386                                       "isolationPFChargedHadrCandsPtSum() + isolationPFGammaCandsEtSum()");
0387       temp2.addParameter<std::string>("selection", "leadPFCand().isNonnull()");
0388       temp2.addParameter<double>("selectionFailValue", 1000.0);
0389       temp2.addParameter<double>("tolerance", 0);
0390       default_cleaners.push_back(temp2);
0391     }
0392
0393     desc.addVPSet("cleaners", vps_description_for_cleaners, default_cleaners);
0394   }
0395
0396   desc.add<int>("verbosity", 0);
0397   desc.add<edm::InputTag>("src", edm::InputTag("combinatoricRecoTaus"));
0398   descriptions.add("RecoTauCleaner", desc);
0399 }
0400
0401 template <>
0402 void RecoTauCleanerImpl<reco::PFTauRefVector>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0403   // there was no cfi file for this plugin
0404 }
0405
0406 #include "FWCore/Framework/interface/MakerMacros.h"
0407
0408 DEFINE_FWK_MODULE(RecoTauCleaner);
0409 DEFINE_FWK_MODULE(RecoTauRefCleaner);