RecoTau/plugins/RecoTauPiZeroStripPlugin.cc

0001 /*
0002  * RecoTauPiZeroStripPlugin
0003  *
0004  * Merges PFGammas in a PFJet into Candidate piZeros defined as
0005  * strips in eta-phi.
0006  *
0007  * Author: Michail Bachtis (University of Wisconsin)
0008  *
0009  * Code modifications: Evan Friis (UC Davis)
0010  *
0011  */
0012 #include <algorithm>
0013 #include <functional>
0014 #include <memory>
0015
0016 #include "RecoTauTag/RecoTau/interface/RecoTauPiZeroPlugins.h"
0017
0018 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0019 #include "DataFormats/Candidate/interface/Candidate.h"
0020 #include "DataFormats/VertexReco/interface/Vertex.h"
0021 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0022 #include "DataFormats/TauReco/interface/RecoTauPiZero.h"
0023 #include "DataFormats/JetReco/interface/PFJet.h"
0024 #include "CommonTools/CandUtils/interface/AddFourMomenta.h"
0025 #include "DataFormats/Math/interface/deltaPhi.h"
0026
0027 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0028 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0029 #include "RecoTauTag/RecoTau/interface/RecoTauVertexAssociator.h"
0030 #include "RecoTauTag/RecoTau/interface/CombinatoricGenerator.h"
0031
0032 namespace reco {
0033   namespace tau {
0034
0035     namespace {
0036       // Apply a hypothesis on the mass of the strips.
0037       math::XYZTLorentzVector applyMassConstraint(const math::XYZTLorentzVector& vec, double mass) {
0038         double factor = sqrt(vec.energy() * vec.energy() - mass * mass) / vec.P();
0039         return math::XYZTLorentzVector(vec.px() * factor, vec.py() * factor, vec.pz() * factor, vec.energy());
0040       }
0041     }  // namespace
0042
0043     class RecoTauPiZeroStripPlugin : public RecoTauPiZeroBuilderPlugin {
0044     public:
0045       explicit RecoTauPiZeroStripPlugin(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC);
0046       ~RecoTauPiZeroStripPlugin() override {}
0047       // Return type is unique_ptr<PiZeroVector>
0048       return_type operator()(const reco::Jet& jet) const override;
0049       // Hook to update PV information
0050       void beginEvent() override;
0051
0052     private:
0053       std::unique_ptr<RecoTauQualityCuts> qcuts_;
0054       RecoTauVertexAssociator vertexAssociator_;
0055
0056       std::vector<int> inputParticleIds_;  //type of candidates to clusterize
0057       double etaAssociationDistance_;      //eta Clustering Association Distance
0058       double phiAssociationDistance_;      //phi Clustering Association Distance
0059
0060       // Parameters for build strip combinations
0061       bool combineStrips_;
0062       int maxStrips_;
0063       double combinatoricStripMassHypo_;
0064
0065       AddFourMomenta p4Builder_;
0066     };
0067
0068     RecoTauPiZeroStripPlugin::RecoTauPiZeroStripPlugin(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0069         : RecoTauPiZeroBuilderPlugin(pset, std::move(iC)),
0070           qcuts_(std::make_unique<RecoTauQualityCuts>(
0071               pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts"))),
0072           vertexAssociator_(pset.getParameter<edm::ParameterSet>("qualityCuts"), std::move(iC)) {
0073       inputParticleIds_ = pset.getParameter<std::vector<int> >("stripCandidatesParticleIds");
0074       etaAssociationDistance_ = pset.getParameter<double>("stripEtaAssociationDistance");
0075       phiAssociationDistance_ = pset.getParameter<double>("stripPhiAssociationDistance");
0076       combineStrips_ = pset.getParameter<bool>("makeCombinatoricStrips");
0077       if (combineStrips_) {
0078         maxStrips_ = pset.getParameter<int>("maxInputStrips");
0079         combinatoricStripMassHypo_ = pset.getParameter<double>("stripMassWhenCombining");
0080       }
0081     }
0082
0083     // Update the primary vertex
0084     void RecoTauPiZeroStripPlugin::beginEvent() { vertexAssociator_.setEvent(*evt()); }
0085
0086     RecoTauPiZeroStripPlugin::return_type RecoTauPiZeroStripPlugin::operator()(const reco::Jet& jet) const {
0087       // Get list of gamma candidates
0088       typedef std::vector<reco::CandidatePtr> CandPtrs;
0089       typedef CandPtrs::iterator CandIter;
0090       PiZeroVector output;
0091
0092       // Get the candidates passing our quality cuts
0093       qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0094       CandPtrs candsVector = qcuts_->filterCandRefs(pfCandidates(jet, inputParticleIds_));
0095       //PFCandPtrs candsVector = qcuts_->filterCandRefs(pfGammas(jet));
0096
0097       // Convert to stl::list to allow fast deletions
0098       std::list<reco::CandidatePtr> cands;
0099       cands.insert(cands.end(), candsVector.begin(), candsVector.end());
0100
0101       while (!cands.empty()) {
0102         // Seed this new strip, and delete it from future strips
0103         CandidatePtr seed = cands.front();
0104         cands.pop_front();
0105
0106         // Add a new candidate to our collection using this seed
0107         auto strip = std::make_unique<RecoTauPiZero>(*seed, RecoTauPiZero::kStrips);
0108         strip->addDaughter(seed);
0109
0110         // Find all other objects in the strip
0111         auto stripCand = cands.begin();
0112         while (stripCand != cands.end()) {
0113           if (fabs(strip->eta() - (*stripCand)->eta()) < etaAssociationDistance_ &&
0114               fabs(deltaPhi(*strip, **stripCand)) < phiAssociationDistance_) {
0115             // Add candidate to strip
0116             strip->addDaughter(*stripCand);
0117             // Update the strips four momenta
0118             p4Builder_.set(*strip);
0119             // Delete this candidate from future strips and move on to
0120             // the next potential candidate
0121             stripCand = cands.erase(stripCand);
0122           } else {
0123             // This candidate isn't compatabile - just move to the next candidate
0124             ++stripCand;
0125           }
0126         }
0127         // Update the vertex
0128         if (strip->daughterPtr(0).isNonnull())
0129           strip->setVertex(strip->daughterPtr(0)->vertex());
0130         output.emplace_back(strip.release());
0131       }
0132
0133       // Check if we want to combine our strips
0134       if (combineStrips_ && output.size() > 1) {
0135         PiZeroVector stripCombinations;
0136         // Sort the output by descending pt
0137         std::sort(output.begin(), output.end(), [&](auto& arg1, auto& arg2) { return arg1->pt() > arg2->pt(); });
0138         // Get the end of interesting set of strips to try and combine
0139         PiZeroVector::const_iterator end_iter = takeNElements(output.begin(), output.end(), maxStrips_);
0140
0141         // Look at all the combinations
0142         for (PiZeroVector::const_iterator firstIter = output.begin(); firstIter != end_iter - 1; ++firstIter) {
0143           for (PiZeroVector::const_iterator secondIter = firstIter + 1; secondIter != end_iter; ++secondIter) {
0144             auto const& first = *firstIter;
0145             auto const& second = *secondIter;
0146             Candidate::LorentzVector firstP4 = first->p4();
0147             Candidate::LorentzVector secondP4 = second->p4();
0148             // If we assume a certain mass for each strip apply it here.
0149             firstP4 = applyMassConstraint(firstP4, combinatoricStripMassHypo_);
0150             secondP4 = applyMassConstraint(secondP4, combinatoricStripMassHypo_);
0151             Candidate::LorentzVector totalP4 = firstP4 + secondP4;
0152             // Make our new combined strip
0153             auto combinedStrips =
0154                 std::make_unique<RecoTauPiZero>(0,
0155                                                 totalP4,
0156                                                 Candidate::Point(0, 0, 0),
0157                                                 //111, 10001, true, RecoTauPiZero::kCombinatoricStrips));
0158                                                 111,
0159                                                 10001,
0160                                                 true,
0161                                                 RecoTauPiZero::kUndefined);
0162
0163             // Now loop over the strip members
0164             for (auto const& gamma : first->daughterPtrVector()) {
0165               combinedStrips->addDaughter(gamma);
0166             }
0167             for (auto const& gamma : second->daughterPtrVector()) {
0168               combinedStrips->addDaughter(gamma);
0169             }
0170             // Update the vertex
0171             if (combinedStrips->daughterPtr(0).isNonnull())
0172               combinedStrips->setVertex(combinedStrips->daughterPtr(0)->vertex());
0173             // Add to our collection of combined strips
0174             stripCombinations.emplace_back(combinedStrips.release());
0175           }
0176         }
0177         // When done doing all the combinations, add the combined strips to the output.
0178         std::move(stripCombinations.begin(), stripCombinations.end(), std::back_inserter(output));
0179       }
0180
0181       return output;
0182     }
0183   }  // namespace tau
0184 }  // namespace reco
0185
0186 #include "FWCore/Framework/interface/MakerMacros.h"
0187 DEFINE_EDM_PLUGIN(RecoTauPiZeroBuilderPluginFactory, reco::tau::RecoTauPiZeroStripPlugin, "RecoTauPiZeroStripPlugin");