RecoTau/plugins/RecoTauBuilderConePlugin.cc

0001 /*
0002  * RecoTauBuilderConePlugin
0003  *
0004  * Build a PFTau using cones defined in DeltaR.
0005  *
0006  * Original Authors: Ludovic Houchu, Simone Gennai
0007  * Modifications: Evan K. Friis
0008  *
0009  */
0010
0011 #include <vector>
0012 #include <algorithm>
0013
0014 #include "RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h"
0015 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0016 #include "RecoTauTag/RecoTau/interface/RecoTauConstructor.h"
0017 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0018
0019 #include "RecoTauTag/RecoTau/interface/ConeTools.h"
0020 #include "RecoTauTag/RecoTau/interface/RecoTauCrossCleaning.h"
0021
0022 #include "DataFormats/TauReco/interface/PFTau.h"
0023 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
0024 #include "DataFormats/TauReco/interface/RecoTauPiZero.h"
0025 #include "DataFormats/VertexReco/interface/Vertex.h"
0026 #include "DataFormats/Math/interface/deltaR.h"
0027
0028 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0029
0030 namespace reco {
0031   namespace tau {
0032
0033     typedef std::vector<RecoTauPiZero> PiZeroList;
0034
0035     class RecoTauBuilderConePlugin : public RecoTauBuilderPlugin {
0036     public:
0037       explicit RecoTauBuilderConePlugin(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC);
0038       ~RecoTauBuilderConePlugin() override {}
0039       // Build a tau from a jet
0040       return_type operator()(const reco::JetBaseRef& jet,
0041                              const std::vector<reco::PFRecoTauChargedHadron>& chargedHadrons,
0042                              const std::vector<RecoTauPiZero>& piZeros,
0043                              const std::vector<CandidatePtr>& regionalExtras) const override;
0044
0045     private:
0046       std::unique_ptr<RecoTauQualityCuts> qcuts_;
0047
0048       bool usePFLeptonsAsChargedHadrons_;
0049
0050       double leadObjecPtThreshold_;
0051
0052       /* String function to extract values from PFJets */
0053       typedef StringObjectFunction<reco::Jet> JetFunc;
0054
0055       // Cone defintions
0056       JetFunc matchingCone_;
0057       JetFunc signalConeChargedHadrons_;
0058       JetFunc isoConeChargedHadrons_;
0059       JetFunc signalConePiZeros_;
0060       StringObjectFunction<reco::PFTau> signalConeSizeToStore_;
0061       JetFunc isoConePiZeros_;
0062       JetFunc signalConeNeutralHadrons_;
0063       JetFunc isoConeNeutralHadrons_;
0064
0065       int maxSignalConeChargedHadrons_;
0066       double minAbsPhotonSumPt_insideSignalCone_;
0067       double minRelPhotonSumPt_insideSignalCone_;
0068
0069       void setTauQuantities(reco::PFTau& aTau,
0070                             double minAbsPhotonSumPt_insideSignalCone = 2.5,
0071                             double minRelPhotonSumPt_insideSignalCone = 0.,
0072                             double minAbsPhotonSumPt_outsideSignalCone = 1.e+9,
0073                             double minRelPhotonSumPt_outsideSignalCone = 1.e+9) const;
0074     };
0075
0076     // ctor - initialize all of our variables
0077     RecoTauBuilderConePlugin::RecoTauBuilderConePlugin(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0078         : RecoTauBuilderPlugin(pset, std::move(iC)),
0079           qcuts_(new RecoTauQualityCuts(pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts"))),
0080           usePFLeptonsAsChargedHadrons_(pset.getParameter<bool>("usePFLeptons")),
0081           leadObjecPtThreshold_(pset.getParameter<double>("leadObjectPt")),
0082           matchingCone_(pset.getParameter<std::string>("matchingCone")),
0083           signalConeChargedHadrons_(pset.getParameter<std::string>("signalConeChargedHadrons")),
0084           isoConeChargedHadrons_(pset.getParameter<std::string>("isoConeChargedHadrons")),
0085           signalConePiZeros_(pset.getParameter<std::string>("signalConePiZeros")),
0086           signalConeSizeToStore_(  //MB: stored inside PFTau and used to compte photon-pt-out-of-signal-cone and DM hypothsis
0087               pset.getParameter<std::string>("signalConePiZeros")),
0088           isoConePiZeros_(pset.getParameter<std::string>("isoConePiZeros")),
0089           signalConeNeutralHadrons_(pset.getParameter<std::string>("signalConeNeutralHadrons")),
0090           isoConeNeutralHadrons_(pset.getParameter<std::string>("isoConeNeutralHadrons")),
0091           maxSignalConeChargedHadrons_(pset.getParameter<int>("maxSignalConeChargedHadrons")),
0092           minAbsPhotonSumPt_insideSignalCone_(pset.getParameter<double>("minAbsPhotonSumPt_insideSignalCone")),
0093           minRelPhotonSumPt_insideSignalCone_(pset.getParameter<double>("minRelPhotonSumPt_insideSignalCone"))
0094
0095     {}
0096
0097     namespace xclean {
0098       // define template specialization for cross-cleaning
0099       template <>
0100       inline void CrossCleanPiZeros<cone::CandPtrDRFilterIter>::initialize(
0101           const cone::CandPtrDRFilterIter& signalTracksBegin, const cone::CandPtrDRFilterIter& signalTracksEnd) {
0102         // Get the list of objects we need to clean
0103         for (cone::CandPtrDRFilterIter signalTrack = signalTracksBegin; signalTrack != signalTracksEnd; ++signalTrack) {
0104           toRemove_.insert(reco::CandidatePtr(*signalTrack));
0105         }
0106       }
0107
0108       template <>
0109       inline void CrossCleanPtrs<PiZeroList::const_iterator>::initialize(const PiZeroList::const_iterator& piZerosBegin,
0110                                                                          const PiZeroList::const_iterator& piZerosEnd) {
0111         for (auto const& ptr : flattenPiZeros(piZerosBegin, piZerosEnd)) {
0112           toRemove_.insert(CandidatePtr(ptr));
0113         }
0114       }
0115     }  // namespace xclean
0116
0117     void RecoTauBuilderConePlugin::setTauQuantities(reco::PFTau& aTau,
0118                                                     double minAbsPhotonSumPt_insideSignalCone,
0119                                                     double minRelPhotonSumPt_insideSignalCone,
0120                                                     double minAbsPhotonSumPt_outsideSignalCone,
0121                                                     double minRelPhotonSumPt_outsideSignalCone) const {
0122       //MB: Set tau quantities which are computed by RecoTauConstructor::get()
0123       //    method with PFRecoTauChargedHadrons not used here
0124
0125       // Set charge
0126       int charge = 0;
0127       int leadCharge = aTau.leadChargedHadrCand().isNonnull() ? aTau.leadChargedHadrCand()->charge() : 0;
0128       const std::vector<reco::CandidatePtr>& pfChs = aTau.signalChargedHadrCands();
0129       for (const reco::CandidatePtr& pfCh : pfChs) {
0130         charge += pfCh->charge();
0131       }
0132       charge = charge == 0 ? leadCharge : charge;
0133       aTau.setCharge(charge);
0134
0135       // Set PDG id
0136       aTau.setPdgId(aTau.charge() < 0 ? 15 : -15);
0137
0138       // Set Decay Mode
0139       PFTau::hadronicDecayMode dm = PFTau::kNull;
0140       unsigned int nPiZeros = 0;
0141       const std::vector<RecoTauPiZero>& piZeros = aTau.signalPiZeroCandidates();
0142       for (const RecoTauPiZero& piZero : piZeros) {
0143         double photonSumPt_insideSignalCone = 0.;
0144         double photonSumPt_outsideSignalCone = 0.;
0145         int numPhotons = piZero.numberOfDaughters();
0146         for (int idxPhoton = 0; idxPhoton < numPhotons; ++idxPhoton) {
0147           const reco::Candidate* photon = piZero.daughter(idxPhoton);
0148           double dR = deltaR(photon->p4(), aTau.p4());
0149           if (dR < aTau.signalConeSize()) {
0150             photonSumPt_insideSignalCone += photon->pt();
0151           } else {
0152             photonSumPt_outsideSignalCone += photon->pt();
0153           }
0154         }
0155         if (photonSumPt_insideSignalCone > minAbsPhotonSumPt_insideSignalCone ||
0156             photonSumPt_insideSignalCone > (minRelPhotonSumPt_insideSignalCone * aTau.pt()) ||
0157             photonSumPt_outsideSignalCone > minAbsPhotonSumPt_outsideSignalCone ||
0158             photonSumPt_outsideSignalCone > (minRelPhotonSumPt_outsideSignalCone * aTau.pt()))
0159           ++nPiZeros;
0160       }
0161       // Find the maximum number of PiZeros our parameterization can hold
0162       const unsigned int maxPiZeros = PFTau::kOneProngNPiZero;
0163       // Determine our track index
0164       unsigned int nCharged = pfChs.size();
0165       if (nCharged > 0) {
0166         unsigned int trackIndex = (nCharged - 1) * (maxPiZeros + 1);
0167         // Check if we handle the given number of tracks
0168         if (trackIndex >= PFTau::kRareDecayMode)
0169           dm = PFTau::kRareDecayMode;
0170         else
0171           dm = static_cast<PFTau::hadronicDecayMode>(trackIndex + std::min(nPiZeros, maxPiZeros));
0172       } else {
0173         dm = PFTau::kNull;
0174       }
0175       aTau.setDecayMode(dm);
0176       return;
0177     }
0178
0179     RecoTauBuilderConePlugin::return_type RecoTauBuilderConePlugin::operator()(
0180         const reco::JetBaseRef& jet,
0181         const std::vector<reco::PFRecoTauChargedHadron>& chargedHadrons,
0182         const std::vector<RecoTauPiZero>& piZeros,
0183         const std::vector<CandidatePtr>& regionalExtras) const {
0184       //std::cout << "<RecoTauBuilderConePlugin::operator()>:" << std::endl;
0185       //std::cout << " jet: Pt = " << jet->pt() << ", eta = " << jet->eta() << ", phi = " << jet->phi() << std::endl;
0186
0187       // Get access to our cone tools
0188       using namespace cone;
0189       // Define output.  We only produce one tau per jet for the cone algo.
0190       output_type output;
0191
0192       // Our tau builder - the true indicates to automatically copy gamma candidates
0193       // from the pizeros.
0194       RecoTauConstructor tau(jet,
0195                              getPFCands(),
0196                              true,
0197                              &signalConeSizeToStore_,
0198                              minAbsPhotonSumPt_insideSignalCone_,
0199                              minRelPhotonSumPt_insideSignalCone_);
0200       // Setup our quality cuts to use the current vertex (supplied by base class)
0201       qcuts_->setPV(primaryVertex(jet));
0202
0203       typedef std::vector<CandidatePtr> CandPtrs;
0204
0205       // Get the PF Charged hadrons + quality cuts
0206       CandPtrs pfchs;
0207       if (!usePFLeptonsAsChargedHadrons_) {
0208         pfchs = qcuts_->filterCandRefs(pfCandidatesByPdgId(*jet, 211));
0209       } else {
0210         // Check if we want to include electrons in muons in "charged hadron"
0211         // collection.  This is the preferred behavior, as the PF lepton selections
0212         // are very loose.
0213         pfchs = qcuts_->filterCandRefs(pfChargedCands(*jet));
0214       }
0215
0216       // CV: sort collection of PF Charged hadrons by descending Pt
0217       std::sort(pfchs.begin(), pfchs.end(), SortPFCandsDescendingPt());
0218
0219       // Get the PF gammas
0220       CandPtrs pfGammas = qcuts_->filterCandRefs(pfCandidatesByPdgId(*jet, 22));
0221       // Neutral hadrons
0222       CandPtrs pfnhs = qcuts_->filterCandRefs(pfCandidatesByPdgId(*jet, 130));
0223
0224       // All the extra junk
0225       CandPtrs regionalJunk = qcuts_->filterCandRefs(regionalExtras);
0226
0227       /***********************************************
0228    ******     Lead Candidate Finding    **********
0229    ***********************************************/
0230
0231       // Define our matching cone and filters
0232       double matchingCone = matchingCone_(*jet);
0233       CandPtrDRFilter matchingConeFilter(jet->p4(), 0, matchingCone);
0234
0235       // Find the maximum PFCharged hadron in the matching cone.  The call to
0236       // PFCandidates always a sorted list, so we can just take the first if it
0237       // if it exists.
0238       CandidatePtr leadPFCH;
0239       CandPtrs::iterator leadPFCH_iter = std::find_if(pfchs.begin(), pfchs.end(), matchingConeFilter);
0240
0241       if (leadPFCH_iter != pfchs.end()) {
0242         leadPFCH = *leadPFCH_iter;
0243         // Set leading candidate
0244         tau.setleadChargedHadrCand(leadPFCH);
0245       } else {
0246         // If there is no leading charged candidate at all, return nothing - the
0247         // producer class that owns the plugin will build a null tau if desired.
0248         return output;
0249       }
0250
0251       // Find the leading neutral candidate
0252       CandidatePtr leadPFGamma;
0253       CandPtrs::iterator leadPFGamma_iter = std::find_if(pfGammas.begin(), pfGammas.end(), matchingConeFilter);
0254
0255       if (leadPFGamma_iter != pfGammas.end()) {
0256         leadPFGamma = *leadPFGamma_iter;
0257         // Set leading neutral candidate
0258         tau.setleadNeutralCand(leadPFGamma);
0259       }
0260
0261       CandidatePtr leadPFCand;
0262       // Always use the leadPFCH if it is above our threshold
0263       if (leadPFCH.isNonnull() && leadPFCH->pt() > leadObjecPtThreshold_) {
0264         leadPFCand = leadPFCH;
0265       } else if (leadPFGamma.isNonnull() && leadPFGamma->pt() > leadObjecPtThreshold_) {
0266         // Otherwise use the lead Gamma if it is above threshold
0267         leadPFCand = leadPFGamma;
0268       } else {
0269         // If both are too low PT, just take the charged one
0270         leadPFCand = leadPFCH;
0271       }
0272
0273       tau.setleadCand(leadPFCand);
0274
0275       // Our cone axis is defined about the lead charged hadron
0276       reco::Candidate::LorentzVector coneAxis = leadPFCH->p4();
0277
0278       /***********************************************
0279    ******     Cone Construction         **********
0280    ***********************************************/
0281
0282       // Define the signal and isolation cone sizes for this jet and build filters
0283       // to select elements in the given DeltaR regions
0284
0285       CandPtrDRFilter signalConePFCHFilter(coneAxis, -0.1, signalConeChargedHadrons_(*jet));
0286       CandPtrDRFilter signalConePFNHFilter(coneAxis, -0.1, signalConeNeutralHadrons_(*jet));
0287       PiZeroDRFilter signalConePiZeroFilter(coneAxis, -0.1, signalConePiZeros_(*jet));
0288
0289       CandPtrDRFilter isoConePFCHFilter(coneAxis, signalConeChargedHadrons_(*jet), isoConeChargedHadrons_(*jet));
0290       CandPtrDRFilter isoConePFGammaFilter(coneAxis, signalConePiZeros_(*jet), isoConePiZeros_(*jet));
0291       CandPtrDRFilter isoConePFNHFilter(coneAxis, signalConeNeutralHadrons_(*jet), isoConeNeutralHadrons_(*jet));
0292       PiZeroDRFilter isoConePiZeroFilter(coneAxis, signalConePiZeros_(*jet), isoConePiZeros_(*jet));
0293
0294       // Additionally make predicates to select the different PF object types
0295       // of the regional junk objects to add to the iso cone.
0296       typedef xclean::PredicateAND<xclean::FilterCandByAbsPdgId, CandPtrDRFilter> RegionalJunkConeAndIdFilter;
0297
0298       xclean::FilterCandByAbsPdgId pfchCandSelector(211);
0299       xclean::FilterCandByAbsPdgId pfgammaCandSelector(22);
0300       xclean::FilterCandByAbsPdgId pfnhCandSelector(130);
0301
0302       // Predicate to select the regional junk in the iso cone by PF id
0303       RegionalJunkConeAndIdFilter pfChargedJunk(pfchCandSelector,  // select charged stuff from junk
0304                                                 isoConePFCHFilter  // only take those in iso cone
0305       );
0306
0307       RegionalJunkConeAndIdFilter pfGammaJunk(pfgammaCandSelector,  // select gammas from junk
0308                                               isoConePFGammaFilter  // only take those in iso cone
0309       );
0310
0311       RegionalJunkConeAndIdFilter pfNeutralJunk(pfnhCandSelector,  // select neutral stuff from junk
0312                                                 isoConePFNHFilter  // select stuff in iso cone
0313       );
0314
0315       // Build filter iterators select the signal charged stuff.
0316       CandPtrDRFilterIter signalPFCHCands_begin(signalConePFCHFilter, pfchs.begin(), pfchs.end());
0317       CandPtrDRFilterIter signalPFCHCands_end(signalConePFCHFilter, pfchs.end(), pfchs.end());
0318       CandPtrs signalPFCHs;
0319       int numSignalPFCHs = 0;
0320       CandPtrs isolationPFCHs;
0321       for (CandPtrDRFilterIter iter = signalPFCHCands_begin; iter != signalPFCHCands_end; ++iter) {
0322         if (numSignalPFCHs < maxSignalConeChargedHadrons_ || maxSignalConeChargedHadrons_ == -1) {
0323           //std::cout << "adding signalPFCH #" << numSignalPFCHs << ": Pt = " << (*iter)->pt() << ", eta = " << (*iter)->eta() << ", phi = " << (*iter)->phi() << std::endl;
0324           signalPFCHs.push_back(*iter);
0325           ++numSignalPFCHs;
0326         } else {
0327           //std::cout << "maxSignalConeChargedHadrons reached"
0328           isolationPFCHs.push_back(*iter);
0329         }
0330       }
0331       CandPtrs::const_iterator signalPFCHs_begin = signalPFCHs.begin();
0332       CandPtrs::const_iterator signalPFCHs_end = signalPFCHs.end();
0333
0334       // Cross clean pi zero content using signal cone charged hadron constituents.
0335       xclean::CrossCleanPiZeros<CandPtrDRFilterIter> piZeroXCleaner(signalPFCHCands_begin, signalPFCHCands_end);
0336       std::vector<reco::RecoTauPiZero> cleanPiZeros = piZeroXCleaner(piZeros);
0337
0338       // For the rest of the constituents, we need to filter any constituents that
0339       // are already contained in the pizeros (i.e. electrons)
0340       xclean::CrossCleanPtrs<PiZeroList::const_iterator> pfCandXCleaner(cleanPiZeros.begin(), cleanPiZeros.end());
0341
0342       auto isolationPFCHCands_begin(boost::make_filter_iterator(
0343           xclean::makePredicateAND(isoConePFCHFilter, pfCandXCleaner), pfchs.begin(), pfchs.end()));
0344       auto isolationPFCHCands_end(boost::make_filter_iterator(
0345           xclean::makePredicateAND(isoConePFCHFilter, pfCandXCleaner), pfchs.end(), pfchs.end()));
0346       for (auto iter = isolationPFCHCands_begin; iter != isolationPFCHCands_end; ++iter) {
0347         isolationPFCHs.push_back(*iter);
0348       }
0349       CandPtrs::const_iterator isolationPFCHs_begin = isolationPFCHs.begin();
0350       CandPtrs::const_iterator isolationPFCHs_end = isolationPFCHs.end();
0351
0352       // Build signal charged hadrons
0353       tau.addPFCands(
0354           RecoTauConstructor::kSignal, RecoTauConstructor::kChargedHadron, signalPFCHs_begin, signalPFCHs_end);
0355
0356       tau.addPFCands(RecoTauConstructor::kSignal,
0357                      RecoTauConstructor::kNeutralHadron,
0358                      boost::make_filter_iterator(
0359                          xclean::makePredicateAND(signalConePFNHFilter, pfCandXCleaner), pfnhs.begin(), pfnhs.end()),
0360                      boost::make_filter_iterator(
0361                          xclean::makePredicateAND(signalConePFNHFilter, pfCandXCleaner), pfnhs.end(), pfnhs.end()));
0362
0363       // Build signal PiZeros
0364       tau.addPiZeros(RecoTauConstructor::kSignal,
0365                      PiZeroDRFilterIter(signalConePiZeroFilter, cleanPiZeros.begin(), cleanPiZeros.end()),
0366                      PiZeroDRFilterIter(signalConePiZeroFilter, cleanPiZeros.end(), cleanPiZeros.end()));
0367
0368       // Build isolation charged hadrons
0369       tau.addPFCands(
0370           RecoTauConstructor::kIsolation, RecoTauConstructor::kChargedHadron, isolationPFCHs_begin, isolationPFCHs_end);
0371
0372       // Add all the stuff in the isolation cone that wasn't in the jet constituents
0373       tau.addPFCands(RecoTauConstructor::kIsolation,
0374                      RecoTauConstructor::kChargedHadron,
0375                      boost::make_filter_iterator(pfChargedJunk, regionalJunk.begin(), regionalJunk.end()),
0376                      boost::make_filter_iterator(pfChargedJunk, regionalJunk.end(), regionalJunk.end()));
0377
0378       // Build isolation neutral hadrons
0379       tau.addPFCands(RecoTauConstructor::kIsolation,
0380                      RecoTauConstructor::kNeutralHadron,
0381                      boost::make_filter_iterator(
0382                          xclean::makePredicateAND(isoConePFNHFilter, pfCandXCleaner), pfnhs.begin(), pfnhs.end()),
0383                      boost::make_filter_iterator(
0384                          xclean::makePredicateAND(isoConePFNHFilter, pfCandXCleaner), pfnhs.end(), pfnhs.end()));
0385
0386       // Add regional stuff not in jet
0387       tau.addPFCands(RecoTauConstructor::kIsolation,
0388                      RecoTauConstructor::kNeutralHadron,
0389                      boost::make_filter_iterator(pfNeutralJunk, regionalJunk.begin(), regionalJunk.end()),
0390                      boost::make_filter_iterator(pfNeutralJunk, regionalJunk.end(), regionalJunk.end()));
0391
0392       // Build isolation PiZeros
0393       tau.addPiZeros(RecoTauConstructor::kIsolation,
0394                      PiZeroDRFilterIter(isoConePiZeroFilter, cleanPiZeros.begin(), cleanPiZeros.end()),
0395                      PiZeroDRFilterIter(isoConePiZeroFilter, cleanPiZeros.end(), cleanPiZeros.end()));
0396
0397       // Add regional stuff not in jet
0398       tau.addPFCands(RecoTauConstructor::kIsolation,
0399                      RecoTauConstructor::kGamma,
0400                      boost::make_filter_iterator(pfGammaJunk, regionalJunk.begin(), regionalJunk.end()),
0401                      boost::make_filter_iterator(pfGammaJunk, regionalJunk.end(), regionalJunk.end()));
0402
0403       // Put our built tau in the output - 'false' indicates don't build the
0404       // leading candidates, we already did that explicitly above.
0405
0406       std::unique_ptr<reco::PFTau> tauPtr = tau.get(false);
0407
0408       // Set event vertex position for tau
0409       reco::VertexRef primaryVertexRef = primaryVertex(*tauPtr);
0410       if (primaryVertexRef.isNonnull())
0411         tauPtr->setVertex(primaryVertexRef->position());
0412
0413       // Set missing tau quantities
0414       setTauQuantities(*tauPtr, minAbsPhotonSumPt_insideSignalCone_, minRelPhotonSumPt_insideSignalCone_);
0415
0416       output.push_back(std::move(tauPtr));
0417       return output;
0418     }
0419   }  // namespace tau
0420 }  // namespace reco
0421
0422 #include "FWCore/Framework/interface/MakerMacros.h"
0423 DEFINE_EDM_PLUGIN(RecoTauBuilderPluginFactory, reco::tau::RecoTauBuilderConePlugin, "RecoTauBuilderConePlugin");