Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-25 04:55:29

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