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