RecoTau/plugins/RecoTauPiZeroProducer.cc

0001 /*
0002  * RecoTauPiZeroProducer
0003  *
0004  * Author: Evan K. Friis, UC Davis
0005  *
0006  * Associates reconstructed PiZeros to PFJets.  The PiZeros are built using one
0007  * or more RecoTauBuilder plugins.  Any overlaps (PiZeros sharing constituents)
0008  * are removed, with the best PiZero candidates taken.  The 'best' are defined
0009  * via the input list of RecoTauPiZeroQualityPlugins, which form a
0010  * lexicograpical ranking.
0011  *
0012  */
0013
0014 #include <algorithm>
0015 #include <functional>
0016 #include <memory>
0017
0018 #include "FWCore/Framework/interface/stream/EDProducer.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/Framework/interface/ConsumesCollector.h"
0024 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0025 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0026
0027 #include "RecoTauTag/RecoTau/interface/RecoTauPiZeroPlugins.h"
0028 #include "RecoTauTag/RecoTau/interface/RecoTauCleaningTools.h"
0029 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0030 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0031
0032 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0033 #include "DataFormats/TauReco/interface/JetPiZeroAssociation.h"
0034 #include "DataFormats/TauReco/interface/RecoTauPiZero.h"
0035 #include "DataFormats/Common/interface/Association.h"
0036
0037 #include "CommonTools/CandUtils/interface/AddFourMomenta.h"
0038 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0039
0040 class RecoTauPiZeroProducer : public edm::stream::EDProducer<> {
0041 public:
0042   typedef reco::tau::RecoTauPiZeroBuilderPlugin Builder;
0043   typedef reco::tau::RecoTauPiZeroQualityPlugin Ranker;
0044
0045   explicit RecoTauPiZeroProducer(const edm::ParameterSet& pset);
0046   ~RecoTauPiZeroProducer() override {}
0047   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0048   void print(const std::vector<reco::RecoTauPiZero>& piZeros, std::ostream& out);
0049
0050   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0051
0052 private:
0053   typedef std::vector<std::unique_ptr<Ranker>> RankerList;
0054   typedef Builder::return_type PiZeroVector;
0055   typedef std::list<std::unique_ptr<reco::RecoTauPiZero>> PiZeroList;
0056
0057   typedef reco::tau::RecoTauLexicographicalRanking<RankerList, reco::RecoTauPiZero> PiZeroPredicate;
0058
0059   std::vector<std::unique_ptr<Builder>> builders_;
0060   RankerList rankers_;
0061   std::unique_ptr<PiZeroPredicate> predicate_;
0062   double piZeroMass_;
0063
0064   // Output selector
0065   std::unique_ptr<StringCutObjectSelector<reco::RecoTauPiZero>> outputSelector_;
0066
0067   //consumes interface
0068   edm::EDGetTokenT<reco::JetView> cand_token;
0069
0070   double minJetPt_;
0071   double maxJetAbsEta_;
0072
0073   int verbosity_;
0074 };
0075
0076 RecoTauPiZeroProducer::RecoTauPiZeroProducer(const edm::ParameterSet& pset) {
0077   cand_token = consumes<reco::JetView>(pset.getParameter<edm::InputTag>("jetSrc"));
0078   minJetPt_ = pset.getParameter<double>("minJetPt");
0079   maxJetAbsEta_ = pset.getParameter<double>("maxJetAbsEta");
0080
0081   typedef std::vector<edm::ParameterSet> VPSet;
0082   // Get the mass hypothesis for the pizeros
0083   piZeroMass_ = pset.getParameter<double>("massHypothesis");
0084
0085   // Get each of our PiZero builders
0086   const VPSet& builders = pset.getParameter<VPSet>("builders");
0087
0088   for (VPSet::const_iterator builderPSet = builders.begin(); builderPSet != builders.end(); ++builderPSet) {
0089     // Get plugin name
0090     const std::string& pluginType = builderPSet->getParameter<std::string>("plugin");
0091     // Build the plugin
0092     builders_.emplace_back(
0093         RecoTauPiZeroBuilderPluginFactory::get()->create(pluginType, *builderPSet, consumesCollector()));
0094   }
0095
0096   // Get each of our quality rankers
0097   const VPSet& rankers = pset.getParameter<VPSet>("ranking");
0098   for (VPSet::const_iterator rankerPSet = rankers.begin(); rankerPSet != rankers.end(); ++rankerPSet) {
0099     const std::string& pluginType = rankerPSet->getParameter<std::string>("plugin");
0100     rankers_.emplace_back(RecoTauPiZeroQualityPluginFactory::get()->create(pluginType, *rankerPSet));
0101   }
0102
0103   // Build the sorting predicate
0104   predicate_ = std::make_unique<PiZeroPredicate>(rankers_);
0105
0106   // now all producers apply a final output selection
0107   std::string selection = pset.getParameter<std::string>("outputSelection");
0108   if (!selection.empty()) {
0109     outputSelector_ = std::make_unique<StringCutObjectSelector<reco::RecoTauPiZero>>(selection);
0110   }
0111
0112   verbosity_ = pset.getParameter<int>("verbosity");
0113
0114   produces<reco::JetPiZeroAssociation>();
0115 }
0116
0117 void RecoTauPiZeroProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0118   // Get a view of our jets via the base candidates
0119   edm::Handle<reco::JetView> jetView;
0120   evt.getByToken(cand_token, jetView);
0121
0122   // Give each of our plugins a chance at doing something with the edm::Event
0123   for (auto& builder : builders_) {
0124     builder->setup(evt, es);
0125   }
0126
0127   // Make our association
0128   auto association = std::make_unique<reco::JetPiZeroAssociation>(reco::JetRefBaseProd(jetView));
0129
0130   // Loop over our jets
0131   size_t nJets = jetView->size();
0132   for (size_t i = 0; i < nJets; ++i) {
0133     const reco::JetBaseRef jet(jetView->refAt(i));
0134
0135     if (jet->pt() - minJetPt_ < 1e-5)
0136       continue;
0137     if (std::abs(jet->eta()) - maxJetAbsEta_ > -1e-5)
0138       continue;
0139     // Build our global list of RecoTauPiZero
0140     PiZeroList dirtyPiZeros;
0141
0142     // Compute the pi zeros from this jet for all the desired algorithms
0143     for (auto const& builder : builders_) {
0144       try {
0145         PiZeroVector result((*builder)(*jet));
0146         std::move(result.begin(), result.end(), std::back_inserter(dirtyPiZeros));
0147       } catch (cms::Exception& exception) {
0148         edm::LogError("BuilderPluginException")
0149             << "Exception caught in builder plugin " << builder->name() << ", rethrowing" << std::endl;
0150         throw exception;
0151       }
0152     }
0153     // Rank the candidates according to our quality plugins
0154     dirtyPiZeros.sort([this](const auto& a, const auto& b) { return (*predicate_)(*a, *b); });
0155
0156     // Keep track of the photons in the clean collection
0157     std::vector<reco::RecoTauPiZero> cleanPiZeros;
0158     std::set<reco::CandidatePtr> photonsInCleanCollection;
0159     while (!dirtyPiZeros.empty()) {
0160       // Pull our candidate pi zero from the front of the list
0161       std::unique_ptr<reco::RecoTauPiZero> toAdd(dirtyPiZeros.front().release());
0162       dirtyPiZeros.pop_front();
0163       // If this doesn't pass our basic selection, discard it.
0164       if (!(*outputSelector_)(*toAdd)) {
0165         continue;
0166       }
0167       // Find the sub-gammas that are not already in the cleaned collection
0168       std::vector<reco::CandidatePtr> uniqueGammas;
0169       std::set_difference(toAdd->daughterPtrVector().begin(),
0170                           toAdd->daughterPtrVector().end(),
0171                           photonsInCleanCollection.begin(),
0172                           photonsInCleanCollection.end(),
0173                           std::back_inserter(uniqueGammas));
0174       // If the pi zero has no unique gammas, discard it.  Note toAdd is deleted
0175       // when it goes out of scope.
0176       if (uniqueGammas.empty()) {
0177         continue;
0178       } else if (uniqueGammas.size() == toAdd->daughterPtrVector().size()) {
0179         // Check if it is composed entirely of unique gammas.  In this case
0180         // immediately add it to the clean collection.
0181         photonsInCleanCollection.insert(toAdd->daughterPtrVector().begin(), toAdd->daughterPtrVector().end());
0182         cleanPiZeros.push_back(*toAdd);
0183       } else {
0184         // Otherwise update the pizero that contains only the unique gammas and
0185         // add it back into the sorted list of dirty PiZeros
0186         toAdd->clearDaughters();
0187         // Add each of the unique daughters back to the pizero
0188         for (auto const& gamma : uniqueGammas) {
0189           toAdd->addDaughter(gamma);
0190         }
0191         // Update the four vector
0192         AddFourMomenta p4Builder_;
0193         p4Builder_.set(*toAdd);
0194         // Put this pi zero back into the collection of sorted dirty pizeros
0195         auto insertionPoint =
0196             std::lower_bound(dirtyPiZeros.begin(), dirtyPiZeros.end(), *toAdd, [this](const auto& a, const auto& b) {
0197               return (*predicate_)(*a, b);
0198             });
0199         dirtyPiZeros.insert(insertionPoint, std::move(toAdd));
0200       }
0201     }
0202     // Apply the mass hypothesis if desired
0203     if (piZeroMass_ >= 0) {
0204       for (auto& cleanPiZero : cleanPiZeros) {
0205         cleanPiZero.setMass(this->piZeroMass_);
0206       };
0207     }
0208     // Add to association
0209     if (verbosity_ >= 2) {
0210       print(cleanPiZeros, std::cout);
0211     }
0212     association->setValue(jet.key(), cleanPiZeros);
0213   }
0214   evt.put(std::move(association));
0215 }
0216
0217 // Print some helpful information
0218 void RecoTauPiZeroProducer::print(const std::vector<reco::RecoTauPiZero>& piZeros, std::ostream& out) {
0219   const unsigned int width = 25;
0220   for (auto const& piZero : piZeros) {
0221     out << piZero;
0222     out << "* Rankers:" << std::endl;
0223     for (auto const& ranker : rankers_) {
0224       out << "* " << std::setiosflags(std::ios::left) << std::setw(width) << ranker->name() << " "
0225           << std::resetiosflags(std::ios::left) << std::setprecision(3) << (*ranker)(piZero);
0226       out << std::endl;
0227     }
0228   }
0229 }
0230
0231 void RecoTauPiZeroProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0232   // common parameter descriptions
0233   edm::ParameterSetDescription desc_ranking;
0234   desc_ranking.add<std::string>("selectionPassFunction", "Func");
0235   desc_ranking.add<double>("selectionFailValue", 1000);
0236   desc_ranking.add<std::string>("selection", "Sel");
0237   desc_ranking.add<std::string>("name", "name");
0238   desc_ranking.add<std::string>("plugin", "plugin");
0239   edm::ParameterSet pset_ranking;
0240   pset_ranking.addParameter<std::string>("selectionPassFunction", "");
0241   pset_ranking.addParameter<double>("selectionFailValue", 1000);
0242   pset_ranking.addParameter<std::string>("selection", "");
0243   pset_ranking.addParameter<std::string>("name", "");
0244   pset_ranking.addParameter<std::string>("plugin", "");
0245   std::vector<edm::ParameterSet> vpsd_ranking;
0246   vpsd_ranking.push_back(pset_ranking);
0247
0248   edm::ParameterSetDescription desc_qualityCuts;
0249   reco::tau::RecoTauQualityCuts::fillDescriptions(desc_qualityCuts);
0250
0251   edm::ParameterSet pset_builders;
0252   pset_builders.addParameter<std::string>("name", "");
0253   pset_builders.addParameter<std::string>("plugin", "");
0254   edm::ParameterSet qualityCuts;
0255   pset_builders.addParameter<edm::ParameterSet>("qualityCuts", qualityCuts);
0256   pset_builders.addParameter<int>("verbosity", 0);
0257
0258   {
0259     // Tailored on ak4PFJetsLegacyHPSPiZeros
0260     edm::ParameterSetDescription desc;
0261     desc.add<double>("massHypothesis", 0.136);
0262     desc.addVPSet("ranking", desc_ranking, vpsd_ranking);
0263     desc.add<int>("verbosity", 0);
0264     desc.add<double>("maxJetAbsEta", 2.5);
0265     desc.add<std::string>("outputSelection", "pt > 0");
0266     desc.add<double>("minJetPt", 14.0);
0267     desc.add<edm::InputTag>("jetSrc", edm::InputTag("ak4PFJets"));
0268
0269     edm::ParameterSetDescription desc_builders;
0270     {
0271       edm::ParameterSetDescription psd0;
0272       psd0.add<std::string>("function", "TMath::Min(0.3, TMath::Max(0.05, [0]*TMath::Power(pT, -[1])))");
0273       psd0.add<double>("par1", 0.707716);
0274       psd0.add<double>("par0", 0.352476);
0275       desc_builders.addOptional<edm::ParameterSetDescription>("stripPhiAssociationDistanceFunc", psd0);
0276     }
0277     {
0278       edm::ParameterSetDescription psd0;
0279       psd0.add<std::string>("function", "TMath::Min(0.15, TMath::Max(0.05, [0]*TMath::Power(pT, -[1])))");
0280       psd0.add<double>("par1", 0.658701);
0281       psd0.add<double>("par0", 0.197077);
0282       desc_builders.addOptional<edm::ParameterSetDescription>("stripEtaAssociationDistanceFunc", psd0);
0283     }
0284     desc_builders.addOptional<double>("stripEtaAssociationDistance", 0.05);
0285     desc_builders.addOptional<double>("stripPhiAssociationDistance", 0.2);
0286
0287     desc_builders.add<edm::ParameterSetDescription>("qualityCuts", desc_qualityCuts);
0288
0289     desc_builders.add<std::string>("name");
0290     desc_builders.add<std::string>("plugin");
0291     desc_builders.add<int>("verbosity", 0);
0292
0293     desc_builders.addOptional<bool>("makeCombinatoricStrips");
0294     desc_builders.addOptional<int>("maxStripBuildIterations");
0295     desc_builders.addOptional<double>("minGammaEtStripAdd");
0296     desc_builders.addOptional<double>("minGammaEtStripSeed");
0297     desc_builders.addOptional<double>("minStripEt");
0298     desc_builders.addOptional<std::vector<int>>("stripCandidatesParticleIds");
0299     desc_builders.addOptional<bool>("updateStripAfterEachDaughter");
0300     desc_builders.addOptional<bool>("applyElecTrackQcuts");
0301
0302     std::vector<edm::ParameterSet> vpsd_builders;
0303     vpsd_builders.push_back(pset_builders);
0304     desc.addVPSet("builders", desc_builders, vpsd_builders);
0305
0306     descriptions.add("recoTauPiZeroProducer", desc);
0307   }
0308 }
0309
0310 #include "FWCore/Framework/interface/MakerMacros.h"
0311 DEFINE_FWK_MODULE(RecoTauPiZeroProducer);