RecoTauBuilderCombinatoricPlugin.cc

Go to the documentation of this file.

#include <vector>

#include "RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h"
#include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"

#include "RecoTauTag/RecoTau/interface/CombinatoricGenerator.h"
#include "RecoTauTag/RecoTau/interface/RecoTauCrossCleaning.h"
#include "RecoTauTag/RecoTau/interface/ConeTools.h"

#include "CommonTools/Utils/interface/StringObjectFunction.h"

#include "DataFormats/TauReco/interface/PFTau.h"
#include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
#include "DataFormats/TauReco/interface/RecoTauPiZero.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "RecoTauTag/RecoTau/interface/RecoTauConstructor.h"
#include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"

#include <algorithm> 

namespace reco { namespace tau {

typedef std::vector<reco::PFRecoTauChargedHadron> ChargedHadronList;
typedef tau::CombinatoricGenerator<ChargedHadronList> ChargedHadronCombo;
typedef std::vector<RecoTauPiZero> PiZeroList;
typedef tau::CombinatoricGenerator<PiZeroList> PiZeroCombo;

class RecoTauBuilderCombinatoricPlugin : public RecoTauBuilderPlugin 
{
 public:
  explicit RecoTauBuilderCombinatoricPlugin(const edm::ParameterSet& pset, edm::ConsumesCollector && iC);
  ~RecoTauBuilderCombinatoricPlugin() override {}

  return_type operator()(
      const reco::JetBaseRef&, 
      const std::vector<reco::PFRecoTauChargedHadron>&, 
      const std::vector<RecoTauPiZero>&, 
      const std::vector<CandidatePtr>&) const override;

 private:
  RecoTauQualityCuts qcuts_;

  double isolationConeSize_;

  struct decayModeInfo 
  {
    uint32_t maxPiZeros_;
    uint32_t maxPFCHs_;
    uint32_t nCharged_;
    uint32_t nPiZeros_;
  };
  std::vector<decayModeInfo> decayModesToBuild_;

  StringObjectFunction<reco::PFTau> signalConeSize_;
  double minAbsPhotonSumPt_insideSignalCone_;
  double minRelPhotonSumPt_insideSignalCone_;
  double minAbsPhotonSumPt_outsideSignalCone_;
  double minRelPhotonSumPt_outsideSignalCone_;

  int verbosity_;
};

RecoTauBuilderCombinatoricPlugin::RecoTauBuilderCombinatoricPlugin(const edm::ParameterSet& pset, edm::ConsumesCollector && iC)
  : RecoTauBuilderPlugin(pset, std::move(iC)),
    qcuts_(pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts")),
    isolationConeSize_(pset.getParameter<double>("isolationConeSize")),
    signalConeSize_(pset.getParameter<std::string>("signalConeSize")),
    minAbsPhotonSumPt_insideSignalCone_(pset.getParameter<double>("minAbsPhotonSumPt_insideSignalCone")),
    minRelPhotonSumPt_insideSignalCone_(pset.getParameter<double>("minRelPhotonSumPt_insideSignalCone")),
    minAbsPhotonSumPt_outsideSignalCone_(pset.getParameter<double>("minAbsPhotonSumPt_outsideSignalCone")),
    minRelPhotonSumPt_outsideSignalCone_(pset.getParameter<double>("minRelPhotonSumPt_outsideSignalCone"))
{
  typedef std::vector<edm::ParameterSet> VPSet;
  const VPSet& decayModes = pset.getParameter<VPSet>("decayModes");
  for ( VPSet::const_iterator decayMode = decayModes.begin();
    decayMode != decayModes.end(); ++decayMode ) {
    decayModeInfo info;
    info.nCharged_ = decayMode->getParameter<uint32_t>("nCharged");
    info.nPiZeros_ = decayMode->getParameter<uint32_t>("nPiZeros");
    info.maxPFCHs_ = decayMode->getParameter<uint32_t>("maxTracks");
    info.maxPiZeros_ = decayMode->getParameter<uint32_t>("maxPiZeros");
    decayModesToBuild_.push_back(info);
  }

  verbosity_ = pset.getParameter<int>("verbosity");
}

// define template specialization for cross-cleaning 
namespace xclean
{
  template<>
  inline void CrossCleanPiZeros<ChargedHadronCombo::combo_iterator>::initialize(const ChargedHadronCombo::combo_iterator& chargedHadronsBegin, const ChargedHadronCombo::combo_iterator& chargedHadronsEnd) 
  {
    // Get the list of objects we need to clean
    for ( ChargedHadronCombo::combo_iterator chargedHadron = chargedHadronsBegin; chargedHadron != chargedHadronsEnd; ++chargedHadron ) {
      // CV: Remove PFGammas that are merged into TauChargedHadrons from isolation PiZeros, but not from signal PiZeros.
      //     The overlap between PFGammas contained in signal PiZeros and merged into TauChargedHadrons
      //     is resolved by RecoTauConstructor::addTauChargedHadron,
      //     which gives preference to PFGammas that are within PiZeros and removes those PFGammas from TauChargedHadrons.
      if ( mode_ == kRemoveChargedDaughterOverlaps ) {
    if ( chargedHadron->getChargedPFCandidate().isNonnull() ) toRemove_.insert(reco::CandidatePtr(chargedHadron->getChargedPFCandidate()));
      } else if ( mode_ == kRemoveChargedAndNeutralDaughterOverlaps ) {
    const reco::CompositePtrCandidate::daughters& daughters = chargedHadron->daughterPtrVector();
    for ( reco::CompositePtrCandidate::daughters::const_iterator daughter = daughters.begin();
          daughter != daughters.end(); ++daughter ) {
      toRemove_.insert(reco::CandidatePtr(*daughter));
    }
      } else assert(0);
    }
  }

  template<>
  inline void CrossCleanPtrs<PiZeroList::const_iterator>::initialize(const PiZeroList::const_iterator& piZerosBegin, const PiZeroList::const_iterator& piZerosEnd) 
  {
    for(auto const& ptr : flattenPiZeros(piZerosBegin, piZerosEnd) ) {
      toRemove_.insert(CandidatePtr(ptr));
    }
  }

  template<>
  inline void CrossCleanPtrs<ChargedHadronCombo::combo_iterator>::initialize(const ChargedHadronCombo::combo_iterator& chargedHadronsBegin, const ChargedHadronCombo::combo_iterator& chargedHadronsEnd) 
  {
    //std::cout << "<CrossCleanPtrs<ChargedHadronCombo>::initialize>:" << std::endl;
    for ( ChargedHadronCombo::combo_iterator chargedHadron = chargedHadronsBegin; chargedHadron != chargedHadronsEnd; ++chargedHadron ) {
      const reco::CompositePtrCandidate::daughters& daughters = chargedHadron->daughterPtrVector();
      for ( reco::CompositePtrCandidate::daughters::const_iterator daughter = daughters.begin();
        daughter != daughters.end(); ++daughter ) { 
    //std::cout << " adding PFCandidate = " << daughter->id() << ":" << daughter->key() << std::endl;
    toRemove_.insert(reco::CandidatePtr(*daughter));
      }
    }
  }

  template<>
  inline void CrossCleanPtrs<ChargedHadronList::const_iterator>::initialize(const ChargedHadronList::const_iterator& chargedHadronsBegin, const ChargedHadronList::const_iterator& chargedHadronsEnd) 
  {
    //std::cout << "<CrossCleanPtrs<ChargedHadronList>::initialize>:" << std::endl;
    for ( ChargedHadronList::const_iterator chargedHadron = chargedHadronsBegin; chargedHadron != chargedHadronsEnd; ++chargedHadron ) {
      const reco::CompositePtrCandidate::daughters& daughters = chargedHadron->daughterPtrVector();
      for ( reco::CompositePtrCandidate::daughters::const_iterator daughter = daughters.begin();
        daughter != daughters.end(); ++daughter ) {
    //std::cout << " adding PFCandidate = " << daughter->id() << ":" << daughter->key() << std::endl;
    toRemove_.insert(reco::CandidatePtr(*daughter));
      }
    }
  }
}

namespace
{
  // auxiliary class for sorting pizeros by descending transverse momentum
  class SortPi0sDescendingPt 
  {
   public:
    bool operator()(const RecoTauPiZero& a, const RecoTauPiZero& b) const 
    {
      return a.pt() > b.pt();
    }  
  };

  double square(double x)
  {
    return x*x;
  }
}

RecoTauBuilderCombinatoricPlugin::return_type
RecoTauBuilderCombinatoricPlugin::operator()(
    const reco::JetBaseRef& jet, 
    const std::vector<reco::PFRecoTauChargedHadron>& chargedHadrons, 
    const std::vector<RecoTauPiZero>& piZeros, 
    const std::vector<CandidatePtr>& regionalExtras) const 
{
  if ( verbosity_ ) {
    std::cout << "<RecoTauBuilderCombinatoricPlugin::operator()>:" << std::endl;
    std::cout << " processing jet: Pt = " << jet->pt() << ", eta = " << jet->eta() << ", phi = " << jet->eta() << ","
          << " mass = " << jet->mass() << ", area = " << jet->jetArea() << std::endl;
  }
  
  // Define output.  
  output_type output;
  
  // Update the primary vertex used by the quality cuts.  The PV is supplied by
  // the base class.
  qcuts_.setPV( primaryVertex(jet) );
  
  typedef std::vector<CandidatePtr> CandPtrs;
  
  if ( verbosity_ ) {
    std::cout << "#chargedHadrons = " << chargedHadrons.size() << std::endl;
    int idx = 0;
    for ( ChargedHadronList::const_iterator chargedHadron = chargedHadrons.begin();
      chargedHadron != chargedHadrons.end(); ++chargedHadron ) {
      std::cout << "chargedHadron #" << idx << ":" << std::endl;
      chargedHadron->print(std::cout);
      ++idx;
    }
    std::cout << "#piZeros = " << piZeros.size() << std::endl;
    idx = 0;
    for ( PiZeroList::const_iterator piZero = piZeros.begin();
      piZero != piZeros.end(); ++piZero ) {
      std::cout << "piZero #" << idx << ":" << std::endl;
      piZero->print(std::cout);
      ++idx;
    }
  }

  CandPtrs pfchs = qcuts_.filterCandRefs(pfChargedCands(*jet));
  CandPtrs pfnhs = qcuts_.filterCandRefs(pfCandidatesByPdgId(*jet, 130));
  CandPtrs pfgammas = qcuts_.filterCandRefs(pfCandidatesByPdgId(*jet, 22));

  CandPtrs regionalJunk = qcuts_.filterCandRefs(regionalExtras);
    
  // Loop over the decay modes we want to build
  for ( std::vector<decayModeInfo>::const_iterator decayMode = decayModesToBuild_.begin();
    decayMode != decayModesToBuild_.end(); ++decayMode ) {
    // Find how many piZeros are in this decay mode
    size_t piZerosToBuild = decayMode->nPiZeros_;
    // Find how many tracks are in this decay mode
    size_t tracksToBuild = decayMode->nCharged_;
    if ( verbosity_ ) {
      std::cout << "piZerosToBuild = " << piZerosToBuild << std::endl;
      std::cout << "#piZeros = " << piZeros.size() << std::endl;
      std::cout << "tracksToBuild = " << tracksToBuild << std::endl;
      std::cout << "#chargedHadrons = " << chargedHadrons.size() << std::endl;      
    }
    
    // Skip decay mode if jet doesn't have the multiplicity to support it
    if ( chargedHadrons.size() < tracksToBuild ) continue;

    // Find the start and end of potential signal tracks
    ChargedHadronList::const_iterator chargedHadron_begin = chargedHadrons.begin();
    ChargedHadronList::const_iterator chargedHadron_end = chargedHadrons.end();
    chargedHadron_end = takeNElements(chargedHadron_begin, chargedHadron_end, decayMode->maxPFCHs_);

    // Build our track combo generator
    ChargedHadronCombo trackCombos(chargedHadron_begin, chargedHadron_end, tracksToBuild);

    CandPtrs::iterator pfch_end = pfchs.end();
    pfch_end = takeNElements(pfchs.begin(), pfch_end, decayMode->maxPFCHs_);

    //-------------------------------------------------------
    // Begin combinatoric loop for this decay mode
    //-------------------------------------------------------
    
    // Loop over the different combinations of tracks
    for ( ChargedHadronCombo::iterator trackCombo = trackCombos.begin();
      trackCombo != trackCombos.end(); ++trackCombo ) {
      xclean::CrossCleanPiZeros<ChargedHadronCombo::combo_iterator> signalPiZeroXCleaner(
          trackCombo->combo_begin(), trackCombo->combo_end(), 
      xclean::CrossCleanPiZeros<ChargedHadronCombo::combo_iterator>::kRemoveChargedDaughterOverlaps);

      PiZeroList cleanSignalPiZeros = signalPiZeroXCleaner(piZeros);
      
      // CV: sort collection of cross-cleaned pi0s by descending Pt
      std::sort(cleanSignalPiZeros.begin(), cleanSignalPiZeros.end(), SortPi0sDescendingPt());
      
      // Skip decay mode if we don't have enough remaining clean pizeros to
      // build it.
      if ( cleanSignalPiZeros.size() < piZerosToBuild ) continue;
      
      // Find the start and end of potential signal tracks
      PiZeroList::iterator signalPiZero_begin = cleanSignalPiZeros.begin();
      PiZeroList::iterator signalPiZero_end = cleanSignalPiZeros.end();
      signalPiZero_end = takeNElements(signalPiZero_begin, signalPiZero_end, decayMode->maxPiZeros_);
      
      // Build our piZero combo generator     
      PiZeroCombo piZeroCombos(signalPiZero_begin, signalPiZero_end, piZerosToBuild);
      // Loop over the different combinations of PiZeros
      for ( PiZeroCombo::iterator piZeroCombo = piZeroCombos.begin();
            piZeroCombo != piZeroCombos.end(); ++piZeroCombo ) {
        // Output tau
        RecoTauConstructor tau(
          jet, getPFCands(), true, 
      &signalConeSize_, 
      minAbsPhotonSumPt_insideSignalCone_, minRelPhotonSumPt_insideSignalCone_, minAbsPhotonSumPt_outsideSignalCone_, minRelPhotonSumPt_outsideSignalCone_);
        // Reserve space in our collections
        tau.reserve(
        RecoTauConstructor::kSignal,
        RecoTauConstructor::kChargedHadron, tracksToBuild);
        tau.reserve(
            RecoTauConstructor::kSignal,
            RecoTauConstructor::kGamma, 2*piZerosToBuild); // k-factor = 2
        tau.reservePiZero(RecoTauConstructor::kSignal, piZerosToBuild);
    
    xclean::CrossCleanPiZeros<ChargedHadronCombo::combo_iterator> isolationPiZeroXCleaner(
          trackCombo->combo_begin(), trackCombo->combo_end(), 
      xclean::CrossCleanPiZeros<ChargedHadronCombo::combo_iterator>::kRemoveChargedAndNeutralDaughterOverlaps);

    PiZeroList precleanedIsolationPiZeros = isolationPiZeroXCleaner(piZeros);
    std::set<reco::CandidatePtr> toRemove;
    for ( PiZeroCombo::combo_iterator signalPiZero = piZeroCombo->combo_begin();
          signalPiZero != piZeroCombo->combo_end(); ++signalPiZero ) {
      toRemove.insert(signalPiZero->daughterPtrVector().begin(), signalPiZero->daughterPtrVector().end());
    }
    PiZeroList cleanIsolationPiZeros;
    for(auto const& precleanedPiZero : precleanedIsolationPiZeros ) {
      std::set<reco::CandidatePtr> toCheck(precleanedPiZero.daughterPtrVector().begin(), precleanedPiZero.daughterPtrVector().end());
      std::vector<reco::CandidatePtr> cleanDaughters;
      std::set_difference(toCheck.begin(), toCheck.end(), toRemove.begin(), toRemove.end(), std::back_inserter(cleanDaughters));
      // CV: piZero is signal piZero if at least one daughter overlaps
      if ( cleanDaughters.size() == precleanedPiZero.daughterPtrVector().size() ) {
        cleanIsolationPiZeros.push_back(precleanedPiZero);
      }
    }
    if ( verbosity_ ) {
      std::cout << "#cleanIsolationPiZeros = " << cleanIsolationPiZeros.size() << std::endl;
      int idx = 0;
      for ( PiZeroList::const_iterator piZero = cleanIsolationPiZeros.begin();
        piZero != cleanIsolationPiZeros.end(); ++piZero ) {
        std::cout << "piZero #" << idx << ":" << std::endl;
        piZero->print(std::cout);
        ++idx;
      }
    }

        // FIXME - are all these reserves okay?  will they get propagated to the
        // dataformat size if they are wrong?
        tau.reserve(
            RecoTauConstructor::kIsolation,
            RecoTauConstructor::kChargedHadron, chargedHadrons.size() - tracksToBuild);
        tau.reserve(
            RecoTauConstructor::kIsolation,
        RecoTauConstructor::kGamma,
        (piZeros.size() - piZerosToBuild)*2);
        tau.reservePiZero(
        RecoTauConstructor::kIsolation,
        (piZeros.size() - piZerosToBuild));

        // Get signal PiZero constituents and add them to the tau.
        // The sub-gammas are automatically added.
        tau.addPiZeros(
            RecoTauConstructor::kSignal,
            piZeroCombo->combo_begin(), piZeroCombo->combo_end());

    // Set signal and isolation components for charged hadrons, after
        // converting them to a PFCandidateRefVector
    //
    // NOTE: signal ChargedHadrons need to be added **after** signal PiZeros
    //       to avoid double-counting PFGammas as part of PiZero and merged with ChargedHadron
    //
        tau.addTauChargedHadrons(
            RecoTauConstructor::kSignal, 
            trackCombo->combo_begin(), trackCombo->combo_end());

        // Now build isolation collections
        // Load our isolation tools
        using namespace reco::tau::cone;
        CandPtrDRFilter isolationConeFilter(tau.p4(), -0.1, isolationConeSize_);

        // Cross cleaning predicate: Remove any PFCandidatePtrs that are contained within existing ChargedHadrons or PiZeros.  
    // The predicate will return false for any object that overlaps with chargedHadrons or cleanPiZeros.
    //  1.) to select charged PFCandidates within jet that are not signalPFChargedHadrons 
    typedef xclean::CrossCleanPtrs<ChargedHadronCombo::combo_iterator> pfChargedHadronXCleanerType;
    pfChargedHadronXCleanerType pfChargedHadronXCleaner_comboChargedHadrons(trackCombo->combo_begin(), trackCombo->combo_end());
    // And this cleaning filter predicate with our Iso cone filter
        xclean::PredicateAND<CandPtrDRFilter, pfChargedHadronXCleanerType> pfCandFilter_comboChargedHadrons(isolationConeFilter, pfChargedHadronXCleaner_comboChargedHadrons);
    //  2.) to select neutral PFCandidates within jet
    xclean::CrossCleanPtrs<ChargedHadronList::const_iterator> pfChargedHadronXCleaner_allChargedHadrons(chargedHadrons.begin(), chargedHadrons.end());
    xclean::CrossCleanPtrs<PiZeroList::const_iterator> piZeroXCleaner(piZeros.begin(), piZeros.end());
    typedef xclean::PredicateAND<xclean::CrossCleanPtrs<ChargedHadronList::const_iterator>, xclean::CrossCleanPtrs<PiZeroList::const_iterator> > pfCandXCleanerType;
        pfCandXCleanerType pfCandXCleaner_allChargedHadrons(pfChargedHadronXCleaner_allChargedHadrons, piZeroXCleaner);
        // And this cleaning filter predicate with our Iso cone filter
        xclean::PredicateAND<CandPtrDRFilter, pfCandXCleanerType> pfCandFilter_allChargedHadrons(isolationConeFilter, pfCandXCleaner_allChargedHadrons);

    ChargedHadronDRFilter isolationConeFilterChargedHadron(tau.p4(), -0.1, isolationConeSize_);
        PiZeroDRFilter isolationConeFilterPiZero(tau.p4(), -0.1, isolationConeSize_);

        // Additionally make predicates to select the different PF object types
        // of the regional junk objects to add
        typedef xclean::PredicateAND<xclean::FilterCandByAbsPdgId,
        CandPtrDRFilter> RegionalJunkConeAndIdFilter;

        xclean::FilterCandByAbsPdgId
          pfchCandSelector(211);
        xclean::FilterCandByAbsPdgId
          pfgammaCandSelector(22);
        xclean::FilterCandByAbsPdgId
          pfnhCandSelector(130);

        RegionalJunkConeAndIdFilter pfChargedJunk(
            pfchCandSelector, // select charged stuff from junk
            isolationConeFilter); // only take those in iso cone

        RegionalJunkConeAndIdFilter pfGammaJunk(
            pfgammaCandSelector, // select gammas from junk
            isolationConeFilter); // only take those in iso cone

        RegionalJunkConeAndIdFilter pfNeutralJunk(
            pfnhCandSelector, // select neutral stuff from junk
            isolationConeFilter); // select stuff in iso cone
    
        tau.addPiZeros(
            RecoTauConstructor::kIsolation,
            boost::make_filter_iterator(
                isolationConeFilterPiZero,
                cleanIsolationPiZeros.begin(), cleanIsolationPiZeros.end()),
            boost::make_filter_iterator(
                isolationConeFilterPiZero,
                cleanIsolationPiZeros.end(), cleanIsolationPiZeros.end()));

        // Filter the isolation candidates in a DR cone
    //
    // NOTE: isolation ChargedHadrons need to be added **after** signal and isolation PiZeros
    //       to avoid double-counting PFGammas as part of PiZero and merged with ChargedHadron
    //
    if ( verbosity_ >= 2 ) {
      std::cout << "adding isolation PFChargedHadrons from trackCombo:" << std::endl;
    }
        tau.addTauChargedHadrons(
            RecoTauConstructor::kIsolation,
            boost::make_filter_iterator(
                isolationConeFilterChargedHadron,
                trackCombo->remainder_begin(), trackCombo->remainder_end()),
            boost::make_filter_iterator(
                isolationConeFilterChargedHadron,
                trackCombo->remainder_end(), trackCombo->remainder_end()));

        // Add all the candidates that weren't included in the combinatoric
        // generation
    if ( verbosity_ >= 2 ) {
      std::cout << "adding isolation PFChargedHadrons not considered in trackCombo:" << std::endl;
    }
        tau.addPFCands(
            RecoTauConstructor::kIsolation, RecoTauConstructor::kChargedHadron,
            boost::make_filter_iterator(
                pfCandFilter_comboChargedHadrons,
                pfch_end, pfchs.end()),
            boost::make_filter_iterator(
                pfCandFilter_comboChargedHadrons,
                pfchs.end(), pfchs.end()));
        // Add all charged candidates that are in the iso cone but weren't in the
        // original PFJet
    if ( verbosity_ >= 2 ) {
      std::cout << "adding isolation PFChargedHadrons from 'regional junk':" << std::endl;
    }
        tau.addPFCands(
            RecoTauConstructor::kIsolation, RecoTauConstructor::kChargedHadron,
            boost::make_filter_iterator(
                pfChargedJunk, regionalJunk.begin(), regionalJunk.end()),
            boost::make_filter_iterator(
                pfChargedJunk, regionalJunk.end(), regionalJunk.end()));
    
    // Add all PFGamma constituents of the jet that are not part of a PiZero
    if ( verbosity_ >= 2 ) {
      std::cout << "adding isolation PFGammas not considered in PiZeros:" << std::endl;
    }
        tau.addPFCands(
        RecoTauConstructor::kIsolation, RecoTauConstructor::kGamma,
            boost::make_filter_iterator(
                pfCandFilter_allChargedHadrons,
        pfgammas.begin(), pfgammas.end()),
            boost::make_filter_iterator(
                pfCandFilter_allChargedHadrons,
        pfgammas.end(), pfgammas.end()));
        // Add all gammas that are in the iso cone but weren't in the
        // orginal PFJet
        tau.addPFCands(
            RecoTauConstructor::kIsolation, RecoTauConstructor::kGamma,
            boost::make_filter_iterator(
                pfGammaJunk, regionalJunk.begin(), regionalJunk.end()),
            boost::make_filter_iterator(
                pfGammaJunk, regionalJunk.end(), regionalJunk.end()));

        // Add all the neutral hadron candidates to the isolation collection
        tau.addPFCands(
            RecoTauConstructor::kIsolation, RecoTauConstructor::kNeutralHadron,
            boost::make_filter_iterator(
                pfCandFilter_allChargedHadrons,
                pfnhs.begin(), pfnhs.end()),
            boost::make_filter_iterator(
                pfCandFilter_allChargedHadrons,
                pfnhs.end(), pfnhs.end()));
        // Add all the neutral hadrons from the region collection that are in
        // the iso cone to the tau
        tau.addPFCands(
            RecoTauConstructor::kIsolation,  RecoTauConstructor::kNeutralHadron,
            boost::make_filter_iterator(
              pfNeutralJunk, regionalJunk.begin(), regionalJunk.end()),
            boost::make_filter_iterator(
              pfNeutralJunk, regionalJunk.end(), regionalJunk.end()));

        std::auto_ptr<reco::PFTau> tauPtr = tau.get(true);
    
    // Set event vertex position for tau
    reco::VertexRef primaryVertexRef = primaryVertex(*tauPtr);
    if ( primaryVertexRef.isNonnull() ) {
      tauPtr->setVertex(primaryVertexRef->position());
    }

    double tauEn = tauPtr->energy();
    double tauPz = tauPtr->pz();
        const double chargedPionMass = 0.13957; // GeV
    double tauMass = std::max(tauPtr->mass(), chargedPionMass);
    double bendCorrMass2 = 0.;
    const std::vector<RecoTauPiZero>& piZeros = tauPtr->signalPiZeroCandidates();
    for (auto const& piZero : piZeros ) {
      double piZeroEn = piZero.energy();
      double piZeroPx = piZero.px();
      double piZeroPy = piZero.py();
      double piZeroPz = piZero.pz();      
      double tau_wo_piZeroPx = tauPtr->px() - piZeroPx;
      double tau_wo_piZeroPy = tauPtr->py() - piZeroPy;
      // CV: Compute effect of varying strip four-vector by eta and phi correction on tau mass
      //    (derrivative of tau mass by strip eta, phi has been computed using Mathematica) 
      bendCorrMass2 += square(((piZeroPz*tauEn - piZeroEn*tauPz)/tauMass)*piZero.bendCorrEta());
      bendCorrMass2 += square(((piZeroPy*tau_wo_piZeroPx - piZeroPx*tau_wo_piZeroPy)/tauMass)*piZero.bendCorrPhi());
    }
    //edm::LogPrint("RecoTauBuilderCombinatoricPlugin") << "bendCorrMass2 = " << sqrt(bendCorrMass2) << std::endl;
    tauPtr->setBendCorrMass(sqrt(bendCorrMass2));

        output.push_back(tauPtr);
      }
    }
  }

  return output.release();
}

}}  // end namespace reco::tau

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_EDM_PLUGIN(RecoTauBuilderPluginFactory,
                  reco::tau::RecoTauBuilderCombinatoricPlugin,
                  "RecoTauBuilderCombinatoricPlugin");