RecoTau/plugins/RecoTauBuilderCombinatoricPlugin.cc

0001 #include <vector>
0002
0003 #include "RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h"
0004 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0005
0006 #include "RecoTauTag/RecoTau/interface/CombinatoricGenerator.h"
0007 #include "RecoTauTag/RecoTau/interface/RecoTauCrossCleaning.h"
0008 #include "RecoTauTag/RecoTau/interface/ConeTools.h"
0009
0010 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0011
0012 #include "DataFormats/TauReco/interface/PFTau.h"
0013 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
0014 #include "DataFormats/TauReco/interface/RecoTauPiZero.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016 #include "RecoTauTag/RecoTau/interface/RecoTauConstructor.h"
0017 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0018
0019 #include <algorithm>
0020
0021 namespace reco {
0022   namespace tau {
0023
0024     typedef std::vector<reco::PFRecoTauChargedHadron> ChargedHadronList;
0025     typedef tau::CombinatoricGenerator<ChargedHadronList> ChargedHadronCombo;
0026     typedef std::vector<RecoTauPiZero> PiZeroList;
0027     typedef tau::CombinatoricGenerator<PiZeroList> PiZeroCombo;
0028
0029     class RecoTauBuilderCombinatoricPlugin : public RecoTauBuilderPlugin {
0030     public:
0031       explicit RecoTauBuilderCombinatoricPlugin(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC);
0032       ~RecoTauBuilderCombinatoricPlugin() override {}
0033
0034       return_type operator()(const reco::JetBaseRef&,
0035                              const std::vector<reco::PFRecoTauChargedHadron>&,
0036                              const std::vector<RecoTauPiZero>&,
0037                              const std::vector<CandidatePtr>&) const override;
0038
0039     private:
0040       std::unique_ptr<RecoTauQualityCuts> qcuts_;
0041
0042       double isolationConeSize_;
0043
0044       struct decayModeInfo {
0045         uint32_t maxPiZeros_;
0046         uint32_t maxPFCHs_;
0047         uint32_t nCharged_;
0048         uint32_t nPiZeros_;
0049       };
0050       std::vector<decayModeInfo> decayModesToBuild_;
0051
0052       StringObjectFunction<reco::PFTau> signalConeSize_;
0053       double minAbsPhotonSumPt_insideSignalCone_;
0054       double minRelPhotonSumPt_insideSignalCone_;
0055       double minAbsPhotonSumPt_outsideSignalCone_;
0056       double minRelPhotonSumPt_outsideSignalCone_;
0057
0058       int verbosity_;
0059     };
0060
0061     RecoTauBuilderCombinatoricPlugin::RecoTauBuilderCombinatoricPlugin(const edm::ParameterSet& pset,
0062                                                                        edm::ConsumesCollector&& iC)
0063         : RecoTauBuilderPlugin(pset, std::move(iC)),
0064           qcuts_(std::make_unique<RecoTauQualityCuts>(
0065               pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts"))),
0066           isolationConeSize_(pset.getParameter<double>("isolationConeSize")),
0067           signalConeSize_(pset.getParameter<std::string>("signalConeSize")),
0068           minAbsPhotonSumPt_insideSignalCone_(pset.getParameter<double>("minAbsPhotonSumPt_insideSignalCone")),
0069           minRelPhotonSumPt_insideSignalCone_(pset.getParameter<double>("minRelPhotonSumPt_insideSignalCone")),
0070           minAbsPhotonSumPt_outsideSignalCone_(pset.getParameter<double>("minAbsPhotonSumPt_outsideSignalCone")),
0071           minRelPhotonSumPt_outsideSignalCone_(pset.getParameter<double>("minRelPhotonSumPt_outsideSignalCone")) {
0072       typedef std::vector<edm::ParameterSet> VPSet;
0073       const VPSet& decayModes = pset.getParameter<VPSet>("decayModes");
0074       for (VPSet::const_iterator decayMode = decayModes.begin(); decayMode != decayModes.end(); ++decayMode) {
0075         decayModeInfo info;
0076         info.nCharged_ = decayMode->getParameter<uint32_t>("nCharged");
0077         info.nPiZeros_ = decayMode->getParameter<uint32_t>("nPiZeros");
0078         info.maxPFCHs_ = decayMode->getParameter<uint32_t>("maxTracks");
0079         info.maxPiZeros_ = decayMode->getParameter<uint32_t>("maxPiZeros");
0080         decayModesToBuild_.push_back(info);
0081       }
0082
0083       verbosity_ = pset.getParameter<int>("verbosity");
0084     }
0085
0086     // define template specialization for cross-cleaning
0087     namespace xclean {
0088       template <>
0089       inline void CrossCleanPiZeros<ChargedHadronCombo::combo_iterator>::initialize(
0090           const ChargedHadronCombo::combo_iterator& chargedHadronsBegin,
0091           const ChargedHadronCombo::combo_iterator& chargedHadronsEnd) {
0092         // Get the list of objects we need to clean
0093         for (ChargedHadronCombo::combo_iterator chargedHadron = chargedHadronsBegin; chargedHadron != chargedHadronsEnd;
0094              ++chargedHadron) {
0095           // CV: Remove PFGammas that are merged into TauChargedHadrons from isolation PiZeros, but not from signal PiZeros.
0096           //     The overlap between PFGammas contained in signal PiZeros and merged into TauChargedHadrons
0097           //     is resolved by RecoTauConstructor::addTauChargedHadron,
0098           //     which gives preference to PFGammas that are within PiZeros and removes those PFGammas from TauChargedHadrons.
0099           if (mode_ == kRemoveChargedDaughterOverlaps) {
0100             if (chargedHadron->getChargedPFCandidate().isNonnull())
0101               toRemove_.insert(reco::CandidatePtr(chargedHadron->getChargedPFCandidate()));
0102           } else if (mode_ == kRemoveChargedAndNeutralDaughterOverlaps) {
0103             const reco::CompositePtrCandidate::daughters& daughters = chargedHadron->daughterPtrVector();
0104             for (reco::CompositePtrCandidate::daughters::const_iterator daughter = daughters.begin();
0105                  daughter != daughters.end();
0106                  ++daughter) {
0107               toRemove_.insert(reco::CandidatePtr(*daughter));
0108             }
0109           } else
0110             assert(0);
0111         }
0112       }
0113
0114       template <>
0115       inline void CrossCleanPtrs<PiZeroList::const_iterator>::initialize(const PiZeroList::const_iterator& piZerosBegin,
0116                                                                          const PiZeroList::const_iterator& piZerosEnd) {
0117         for (auto const& ptr : flattenPiZeros(piZerosBegin, piZerosEnd)) {
0118           toRemove_.insert(CandidatePtr(ptr));
0119         }
0120       }
0121
0122       template <>
0123       inline void CrossCleanPtrs<ChargedHadronCombo::combo_iterator>::initialize(
0124           const ChargedHadronCombo::combo_iterator& chargedHadronsBegin,
0125           const ChargedHadronCombo::combo_iterator& chargedHadronsEnd) {
0126         //std::cout << "<CrossCleanPtrs<ChargedHadronCombo>::initialize>:" << std::endl;
0127         for (ChargedHadronCombo::combo_iterator chargedHadron = chargedHadronsBegin; chargedHadron != chargedHadronsEnd;
0128              ++chargedHadron) {
0129           const reco::CompositePtrCandidate::daughters& daughters = chargedHadron->daughterPtrVector();
0130           for (reco::CompositePtrCandidate::daughters::const_iterator daughter = daughters.begin();
0131                daughter != daughters.end();
0132                ++daughter) {
0133             //std::cout << " adding PFCandidate = " << daughter->id() << ":" << daughter->key() << std::endl;
0134             toRemove_.insert(reco::CandidatePtr(*daughter));
0135           }
0136         }
0137       }
0138
0139       template <>
0140       inline void CrossCleanPtrs<ChargedHadronList::const_iterator>::initialize(
0141           const ChargedHadronList::const_iterator& chargedHadronsBegin,
0142           const ChargedHadronList::const_iterator& chargedHadronsEnd) {
0143         //std::cout << "<CrossCleanPtrs<ChargedHadronList>::initialize>:" << std::endl;
0144         for (ChargedHadronList::const_iterator chargedHadron = chargedHadronsBegin; chargedHadron != chargedHadronsEnd;
0145              ++chargedHadron) {
0146           const reco::CompositePtrCandidate::daughters& daughters = chargedHadron->daughterPtrVector();
0147           for (reco::CompositePtrCandidate::daughters::const_iterator daughter = daughters.begin();
0148                daughter != daughters.end();
0149                ++daughter) {
0150             //std::cout << " adding PFCandidate = " << daughter->id() << ":" << daughter->key() << std::endl;
0151             toRemove_.insert(reco::CandidatePtr(*daughter));
0152           }
0153         }
0154       }
0155     }  // namespace xclean
0156
0157     namespace {
0158       // auxiliary class for sorting pizeros by descending transverse momentum
0159       class SortPi0sDescendingPt {
0160       public:
0161         bool operator()(const RecoTauPiZero& a, const RecoTauPiZero& b) const { return a.pt() > b.pt(); }
0162       };
0163
0164       double square(double x) { return x * x; }
0165     }  // namespace
0166
0167     RecoTauBuilderCombinatoricPlugin::return_type RecoTauBuilderCombinatoricPlugin::operator()(
0168         const reco::JetBaseRef& jet,
0169         const std::vector<reco::PFRecoTauChargedHadron>& chargedHadrons,
0170         const std::vector<RecoTauPiZero>& piZeros,
0171         const std::vector<CandidatePtr>& regionalExtras) const {
0172       if (verbosity_) {
0173         std::cout << "<RecoTauBuilderCombinatoricPlugin::operator()>:" << std::endl;
0174         std::cout << " processing jet: Pt = " << jet->pt() << ", eta = " << jet->eta() << ", phi = " << jet->eta()
0175                   << ","
0176                   << " mass = " << jet->mass() << ", area = " << jet->jetArea() << std::endl;
0177       }
0178
0179       // Define output.
0180       output_type output;
0181
0182       // Update the primary vertex used by the quality cuts.  The PV is supplied by
0183       // the base class.
0184       qcuts_->setPV(primaryVertex(jet));
0185
0186       typedef std::vector<CandidatePtr> CandPtrs;
0187
0188       if (verbosity_) {
0189         std::cout << "#chargedHadrons = " << chargedHadrons.size() << std::endl;
0190         int idx = 0;
0191         for (ChargedHadronList::const_iterator chargedHadron = chargedHadrons.begin();
0192              chargedHadron != chargedHadrons.end();
0193              ++chargedHadron) {
0194           std::cout << "chargedHadron #" << idx << ":" << std::endl;
0195           chargedHadron->print(std::cout);
0196           ++idx;
0197         }
0198         std::cout << "#piZeros = " << piZeros.size() << std::endl;
0199         idx = 0;
0200         for (PiZeroList::const_iterator piZero = piZeros.begin(); piZero != piZeros.end(); ++piZero) {
0201           std::cout << "piZero #" << idx << ":" << std::endl;
0202           piZero->print(std::cout);
0203           ++idx;
0204         }
0205       }
0206
0207       CandPtrs pfchs = qcuts_->filterCandRefs(pfChargedCands(*jet));
0208       CandPtrs pfnhs = qcuts_->filterCandRefs(pfCandidatesByPdgId(*jet, 130));
0209       CandPtrs pfgammas = qcuts_->filterCandRefs(pfCandidatesByPdgId(*jet, 22));
0210
0211       /// Apply quality cuts to the regional junk around the jet.  Note that the
0212       /// particle contents of the junk is exclusive to the jet content.
0213       CandPtrs regionalJunk = qcuts_->filterCandRefs(regionalExtras);
0214
0215       // Loop over the decay modes we want to build
0216       for (std::vector<decayModeInfo>::const_iterator decayMode = decayModesToBuild_.begin();
0217            decayMode != decayModesToBuild_.end();
0218            ++decayMode) {
0219         // Find how many piZeros are in this decay mode
0220         size_t piZerosToBuild = decayMode->nPiZeros_;
0221         // Find how many tracks are in this decay mode
0222         size_t tracksToBuild = decayMode->nCharged_;
0223         if (verbosity_) {
0224           std::cout << "piZerosToBuild = " << piZerosToBuild << std::endl;
0225           std::cout << "#piZeros = " << piZeros.size() << std::endl;
0226           std::cout << "tracksToBuild = " << tracksToBuild << std::endl;
0227           std::cout << "#chargedHadrons = " << chargedHadrons.size() << std::endl;
0228         }
0229
0230         // Skip decay mode if jet doesn't have the multiplicity to support it
0231         if (chargedHadrons.size() < tracksToBuild)
0232           continue;
0233
0234         // Find the start and end of potential signal tracks
0235         ChargedHadronList::const_iterator chargedHadron_begin = chargedHadrons.begin();
0236         ChargedHadronList::const_iterator chargedHadron_end = chargedHadrons.end();
0237         chargedHadron_end = takeNElements(chargedHadron_begin, chargedHadron_end, decayMode->maxPFCHs_);
0238
0239         // Build our track combo generator
0240         ChargedHadronCombo trackCombos(chargedHadron_begin, chargedHadron_end, tracksToBuild);
0241
0242         CandPtrs::iterator pfch_end = pfchs.end();
0243         pfch_end = takeNElements(pfchs.begin(), pfch_end, decayMode->maxPFCHs_);
0244
0245         //-------------------------------------------------------
0246         // Begin combinatoric loop for this decay mode
0247         //-------------------------------------------------------
0248
0249         // Loop over the different combinations of tracks
0250         for (ChargedHadronCombo::iterator trackCombo = trackCombos.begin(); trackCombo != trackCombos.end();
0251              ++trackCombo) {
0252           xclean::CrossCleanPiZeros<ChargedHadronCombo::combo_iterator> signalPiZeroXCleaner(
0253               trackCombo->combo_begin(),
0254               trackCombo->combo_end(),
0255               xclean::CrossCleanPiZeros<ChargedHadronCombo::combo_iterator>::kRemoveChargedDaughterOverlaps);
0256
0257           PiZeroList cleanSignalPiZeros = signalPiZeroXCleaner(piZeros);
0258
0259           // CV: sort collection of cross-cleaned pi0s by descending Pt
0260           std::sort(cleanSignalPiZeros.begin(), cleanSignalPiZeros.end(), SortPi0sDescendingPt());
0261
0262           // Skip decay mode if we don't have enough remaining clean pizeros to
0263           // build it.
0264           if (cleanSignalPiZeros.size() < piZerosToBuild)
0265             continue;
0266
0267           // Find the start and end of potential signal tracks
0268           PiZeroList::iterator signalPiZero_begin = cleanSignalPiZeros.begin();
0269           PiZeroList::iterator signalPiZero_end = cleanSignalPiZeros.end();
0270           signalPiZero_end = takeNElements(signalPiZero_begin, signalPiZero_end, decayMode->maxPiZeros_);
0271
0272           // Build our piZero combo generator
0273           PiZeroCombo piZeroCombos(signalPiZero_begin, signalPiZero_end, piZerosToBuild);
0274           // Loop over the different combinations of PiZeros
0275           for (PiZeroCombo::iterator piZeroCombo = piZeroCombos.begin(); piZeroCombo != piZeroCombos.end();
0276                ++piZeroCombo) {
0277             // Output tau
0278             RecoTauConstructor tau(jet,
0279                                    getPFCands(),
0280                                    true,
0281                                    &signalConeSize_,
0282                                    minAbsPhotonSumPt_insideSignalCone_,
0283                                    minRelPhotonSumPt_insideSignalCone_,
0284                                    minAbsPhotonSumPt_outsideSignalCone_,
0285                                    minRelPhotonSumPt_outsideSignalCone_);
0286             // Reserve space in our collections
0287             tau.reserve(RecoTauConstructor::kSignal, RecoTauConstructor::kChargedHadron, tracksToBuild);
0288             tau.reserve(RecoTauConstructor::kSignal, RecoTauConstructor::kGamma, 2 * piZerosToBuild);  // k-factor = 2
0289             tau.reservePiZero(RecoTauConstructor::kSignal, piZerosToBuild);
0290
0291             xclean::CrossCleanPiZeros<ChargedHadronCombo::combo_iterator> isolationPiZeroXCleaner(
0292                 trackCombo->combo_begin(),
0293                 trackCombo->combo_end(),
0294                 xclean::CrossCleanPiZeros<ChargedHadronCombo::combo_iterator>::kRemoveChargedAndNeutralDaughterOverlaps);
0295
0296             PiZeroList precleanedIsolationPiZeros = isolationPiZeroXCleaner(piZeros);
0297             std::set<reco::CandidatePtr> toRemove;
0298             for (PiZeroCombo::combo_iterator signalPiZero = piZeroCombo->combo_begin();
0299                  signalPiZero != piZeroCombo->combo_end();
0300                  ++signalPiZero) {
0301               toRemove.insert(signalPiZero->daughterPtrVector().begin(), signalPiZero->daughterPtrVector().end());
0302             }
0303             PiZeroList cleanIsolationPiZeros;
0304             for (auto const& precleanedPiZero : precleanedIsolationPiZeros) {
0305               std::set<reco::CandidatePtr> toCheck(precleanedPiZero.daughterPtrVector().begin(),
0306                                                    precleanedPiZero.daughterPtrVector().end());
0307               std::vector<reco::CandidatePtr> cleanDaughters;
0308               std::set_difference(
0309                   toCheck.begin(), toCheck.end(), toRemove.begin(), toRemove.end(), std::back_inserter(cleanDaughters));
0310               // CV: piZero is signal piZero if at least one daughter overlaps
0311               if (cleanDaughters.size() == precleanedPiZero.daughterPtrVector().size()) {
0312                 cleanIsolationPiZeros.push_back(precleanedPiZero);
0313               }
0314             }
0315             if (verbosity_) {
0316               std::cout << "#cleanIsolationPiZeros = " << cleanIsolationPiZeros.size() << std::endl;
0317               int idx = 0;
0318               for (PiZeroList::const_iterator piZero = cleanIsolationPiZeros.begin();
0319                    piZero != cleanIsolationPiZeros.end();
0320                    ++piZero) {
0321                 std::cout << "piZero #" << idx << ":" << std::endl;
0322                 piZero->print(std::cout);
0323                 ++idx;
0324               }
0325             }
0326
0327             // FIXME - are all these reserves okay?  will they get propagated to the
0328             // dataformat size if they are wrong?
0329             tau.reserve(RecoTauConstructor::kIsolation,
0330                         RecoTauConstructor::kChargedHadron,
0331                         chargedHadrons.size() - tracksToBuild);
0332             tau.reserve(
0333                 RecoTauConstructor::kIsolation, RecoTauConstructor::kGamma, (piZeros.size() - piZerosToBuild) * 2);
0334             tau.reservePiZero(RecoTauConstructor::kIsolation, (piZeros.size() - piZerosToBuild));
0335
0336             // Get signal PiZero constituents and add them to the tau.
0337             // The sub-gammas are automatically added.
0338             tau.addPiZeros(RecoTauConstructor::kSignal, piZeroCombo->combo_begin(), piZeroCombo->combo_end());
0339
0340             // Set signal and isolation components for charged hadrons, after
0341             // converting them to a PFCandidateRefVector
0342             //
0343             // NOTE: signal ChargedHadrons need to be added **after** signal PiZeros
0344             //       to avoid double-counting PFGammas as part of PiZero and merged with ChargedHadron
0345             //
0346             tau.addTauChargedHadrons(RecoTauConstructor::kSignal, trackCombo->combo_begin(), trackCombo->combo_end());
0347
0348             // Now build isolation collections
0349             // Load our isolation tools
0350             using namespace reco::tau::cone;
0351             CandPtrDRFilter isolationConeFilter(tau.p4(), -0.1, isolationConeSize_);
0352
0353             // Cross cleaning predicate: Remove any PFCandidatePtrs that are contained within existing ChargedHadrons or PiZeros.
0354             // The predicate will return false for any object that overlaps with chargedHadrons or cleanPiZeros.
0355             //  1.) to select charged PFCandidates within jet that are not signalPFChargedHadrons
0356             typedef xclean::CrossCleanPtrs<ChargedHadronCombo::combo_iterator> pfChargedHadronXCleanerType;
0357             pfChargedHadronXCleanerType pfChargedHadronXCleaner_comboChargedHadrons(trackCombo->combo_begin(),
0358                                                                                     trackCombo->combo_end());
0359             // And this cleaning filter predicate with our Iso cone filter
0360             xclean::PredicateAND<CandPtrDRFilter, pfChargedHadronXCleanerType> pfCandFilter_comboChargedHadrons(
0361                 isolationConeFilter, pfChargedHadronXCleaner_comboChargedHadrons);
0362             //  2.) to select neutral PFCandidates within jet
0363             xclean::CrossCleanPtrs<ChargedHadronList::const_iterator> pfChargedHadronXCleaner_allChargedHadrons(
0364                 chargedHadrons.begin(), chargedHadrons.end());
0365             xclean::CrossCleanPtrs<PiZeroList::const_iterator> piZeroXCleaner(piZeros.begin(), piZeros.end());
0366             typedef xclean::PredicateAND<xclean::CrossCleanPtrs<ChargedHadronList::const_iterator>,
0367                                          xclean::CrossCleanPtrs<PiZeroList::const_iterator> >
0368                 pfCandXCleanerType;
0369             pfCandXCleanerType pfCandXCleaner_allChargedHadrons(pfChargedHadronXCleaner_allChargedHadrons,
0370                                                                 piZeroXCleaner);
0371             // And this cleaning filter predicate with our Iso cone filter
0372             xclean::PredicateAND<CandPtrDRFilter, pfCandXCleanerType> pfCandFilter_allChargedHadrons(
0373                 isolationConeFilter, pfCandXCleaner_allChargedHadrons);
0374
0375             ChargedHadronDRFilter isolationConeFilterChargedHadron(tau.p4(), -0.1, isolationConeSize_);
0376             PiZeroDRFilter isolationConeFilterPiZero(tau.p4(), -0.1, isolationConeSize_);
0377
0378             // Additionally make predicates to select the different PF object types
0379             // of the regional junk objects to add
0380             typedef xclean::PredicateAND<xclean::FilterCandByAbsPdgId, CandPtrDRFilter> RegionalJunkConeAndIdFilter;
0381
0382             xclean::FilterCandByAbsPdgId pfchCandSelector(211);
0383             xclean::FilterCandByAbsPdgId pfgammaCandSelector(22);
0384             xclean::FilterCandByAbsPdgId pfnhCandSelector(130);
0385
0386             RegionalJunkConeAndIdFilter pfChargedJunk(pfchCandSelector,      // select charged stuff from junk
0387                                                       isolationConeFilter);  // only take those in iso cone
0388
0389             RegionalJunkConeAndIdFilter pfGammaJunk(pfgammaCandSelector,   // select gammas from junk
0390                                                     isolationConeFilter);  // only take those in iso cone
0391
0392             RegionalJunkConeAndIdFilter pfNeutralJunk(pfnhCandSelector,      // select neutral stuff from junk
0393                                                       isolationConeFilter);  // select stuff in iso cone
0394
0395             tau.addPiZeros(RecoTauConstructor::kIsolation,
0396                            boost::make_filter_iterator(
0397                                isolationConeFilterPiZero, cleanIsolationPiZeros.begin(), cleanIsolationPiZeros.end()),
0398                            boost::make_filter_iterator(
0399                                isolationConeFilterPiZero, cleanIsolationPiZeros.end(), cleanIsolationPiZeros.end()));
0400
0401             // Filter the isolation candidates in a DR cone
0402             //
0403             // NOTE: isolation ChargedHadrons need to be added **after** signal and isolation PiZeros
0404             //       to avoid double-counting PFGammas as part of PiZero and merged with ChargedHadron
0405             //
0406             if (verbosity_ >= 2) {
0407               std::cout << "adding isolation PFChargedHadrons from trackCombo:" << std::endl;
0408             }
0409             tau.addTauChargedHadrons(
0410                 RecoTauConstructor::kIsolation,
0411                 boost::make_filter_iterator(
0412                     isolationConeFilterChargedHadron, trackCombo->remainder_begin(), trackCombo->remainder_end()),
0413                 boost::make_filter_iterator(
0414                     isolationConeFilterChargedHadron, trackCombo->remainder_end(), trackCombo->remainder_end()));
0415
0416             // Add all the candidates that weren't included in the combinatoric
0417             // generation
0418             if (verbosity_ >= 2) {
0419               std::cout << "adding isolation PFChargedHadrons not considered in trackCombo:" << std::endl;
0420             }
0421             tau.addPFCands(RecoTauConstructor::kIsolation,
0422                            RecoTauConstructor::kChargedHadron,
0423                            boost::make_filter_iterator(pfCandFilter_comboChargedHadrons, pfch_end, pfchs.end()),
0424                            boost::make_filter_iterator(pfCandFilter_comboChargedHadrons, pfchs.end(), pfchs.end()));
0425             // Add all charged candidates that are in the iso cone but weren't in the
0426             // original PFJet
0427             if (verbosity_ >= 2) {
0428               std::cout << "adding isolation PFChargedHadrons from 'regional junk':" << std::endl;
0429             }
0430             tau.addPFCands(RecoTauConstructor::kIsolation,
0431                            RecoTauConstructor::kChargedHadron,
0432                            boost::make_filter_iterator(pfChargedJunk, regionalJunk.begin(), regionalJunk.end()),
0433                            boost::make_filter_iterator(pfChargedJunk, regionalJunk.end(), regionalJunk.end()));
0434
0435             // Add all PFGamma constituents of the jet that are not part of a PiZero
0436             if (verbosity_ >= 2) {
0437               std::cout << "adding isolation PFGammas not considered in PiZeros:" << std::endl;
0438             }
0439             tau.addPFCands(
0440                 RecoTauConstructor::kIsolation,
0441                 RecoTauConstructor::kGamma,
0442                 boost::make_filter_iterator(pfCandFilter_allChargedHadrons, pfgammas.begin(), pfgammas.end()),
0443                 boost::make_filter_iterator(pfCandFilter_allChargedHadrons, pfgammas.end(), pfgammas.end()));
0444             // Add all gammas that are in the iso cone but weren't in the
0445             // orginal PFJet
0446             tau.addPFCands(RecoTauConstructor::kIsolation,
0447                            RecoTauConstructor::kGamma,
0448                            boost::make_filter_iterator(pfGammaJunk, regionalJunk.begin(), regionalJunk.end()),
0449                            boost::make_filter_iterator(pfGammaJunk, regionalJunk.end(), regionalJunk.end()));
0450
0451             // Add all the neutral hadron candidates to the isolation collection
0452             tau.addPFCands(RecoTauConstructor::kIsolation,
0453                            RecoTauConstructor::kNeutralHadron,
0454                            boost::make_filter_iterator(pfCandFilter_allChargedHadrons, pfnhs.begin(), pfnhs.end()),
0455                            boost::make_filter_iterator(pfCandFilter_allChargedHadrons, pfnhs.end(), pfnhs.end()));
0456             // Add all the neutral hadrons from the region collection that are in
0457             // the iso cone to the tau
0458             tau.addPFCands(RecoTauConstructor::kIsolation,
0459                            RecoTauConstructor::kNeutralHadron,
0460                            boost::make_filter_iterator(pfNeutralJunk, regionalJunk.begin(), regionalJunk.end()),
0461                            boost::make_filter_iterator(pfNeutralJunk, regionalJunk.end(), regionalJunk.end()));
0462
0463             std::unique_ptr<reco::PFTau> tauPtr = tau.get(true);
0464
0465             // Set event vertex position for tau
0466             reco::VertexRef primaryVertexRef = primaryVertex(*tauPtr);
0467             if (primaryVertexRef.isNonnull()) {
0468               tauPtr->setVertex(primaryVertexRef->position());
0469             }
0470
0471             double tauEn = tauPtr->energy();
0472             double tauPz = tauPtr->pz();
0473             const double chargedPionMass = 0.13957;  // GeV
0474             double tauMass = std::max(tauPtr->mass(), chargedPionMass);
0475             double bendCorrMass2 = 0.;
0476             const std::vector<RecoTauPiZero>& piZeros = tauPtr->signalPiZeroCandidates();
0477             for (auto const& piZero : piZeros) {
0478               double piZeroEn = piZero.energy();
0479               double piZeroPx = piZero.px();
0480               double piZeroPy = piZero.py();
0481               double piZeroPz = piZero.pz();
0482               double tau_wo_piZeroPx = tauPtr->px() - piZeroPx;
0483               double tau_wo_piZeroPy = tauPtr->py() - piZeroPy;
0484               // CV: Compute effect of varying strip four-vector by eta and phi correction on tau mass
0485               //    (derrivative of tau mass by strip eta, phi has been computed using Mathematica)
0486               bendCorrMass2 += square(((piZeroPz * tauEn - piZeroEn * tauPz) / tauMass) * piZero.bendCorrEta());
0487               bendCorrMass2 +=
0488                   square(((piZeroPy * tau_wo_piZeroPx - piZeroPx * tau_wo_piZeroPy) / tauMass) * piZero.bendCorrPhi());
0489             }
0490             //edm::LogPrint("RecoTauBuilderCombinatoricPlugin") << "bendCorrMass2 = " << sqrt(bendCorrMass2) << std::endl;
0491             tauPtr->setBendCorrMass(sqrt(bendCorrMass2));
0492
0493             output.push_back(std::move(tauPtr));
0494           }
0495         }
0496       }
0497
0498       return output;
0499     }
0500
0501   }  // namespace tau
0502 }  // namespace reco
0503
0504 #include "FWCore/Framework/interface/MakerMacros.h"
0505 DEFINE_EDM_PLUGIN(RecoTauBuilderPluginFactory,
0506                   reco::tau::RecoTauBuilderCombinatoricPlugin,
0507                   "RecoTauBuilderCombinatoricPlugin");