Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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_(new RecoTauQualityCuts(pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts"))),
0071           vertexAssociator_(pset.getParameter<edm::ParameterSet>("qualityCuts"), std::move(iC)) {
0072       inputParticleIds_ = pset.getParameter<std::vector<int> >("stripCandidatesParticleIds");
0073       etaAssociationDistance_ = pset.getParameter<double>("stripEtaAssociationDistance");
0074       phiAssociationDistance_ = pset.getParameter<double>("stripPhiAssociationDistance");
0075       combineStrips_ = pset.getParameter<bool>("makeCombinatoricStrips");
0076       if (combineStrips_) {
0077         maxStrips_ = pset.getParameter<int>("maxInputStrips");
0078         combinatoricStripMassHypo_ = pset.getParameter<double>("stripMassWhenCombining");
0079       }
0080     }
0081 
0082     // Update the primary vertex
0083     void RecoTauPiZeroStripPlugin::beginEvent() { vertexAssociator_.setEvent(*evt()); }
0084 
0085     RecoTauPiZeroStripPlugin::return_type RecoTauPiZeroStripPlugin::operator()(const reco::Jet& jet) const {
0086       // Get list of gamma candidates
0087       typedef std::vector<reco::CandidatePtr> CandPtrs;
0088       typedef CandPtrs::iterator CandIter;
0089       PiZeroVector output;
0090 
0091       // Get the candidates passing our quality cuts
0092       qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0093       CandPtrs candsVector = qcuts_->filterCandRefs(pfCandidates(jet, inputParticleIds_));
0094       //PFCandPtrs candsVector = qcuts_->filterCandRefs(pfGammas(jet));
0095 
0096       // Convert to stl::list to allow fast deletions
0097       std::list<reco::CandidatePtr> cands;
0098       cands.insert(cands.end(), candsVector.begin(), candsVector.end());
0099 
0100       while (!cands.empty()) {
0101         // Seed this new strip, and delete it from future strips
0102         CandidatePtr seed = cands.front();
0103         cands.pop_front();
0104 
0105         // Add a new candidate to our collection using this seed
0106         auto strip = std::make_unique<RecoTauPiZero>(*seed, RecoTauPiZero::kStrips);
0107         strip->addDaughter(seed);
0108 
0109         // Find all other objects in the strip
0110         auto stripCand = cands.begin();
0111         while (stripCand != cands.end()) {
0112           if (fabs(strip->eta() - (*stripCand)->eta()) < etaAssociationDistance_ &&
0113               fabs(deltaPhi(*strip, **stripCand)) < phiAssociationDistance_) {
0114             // Add candidate to strip
0115             strip->addDaughter(*stripCand);
0116             // Update the strips four momenta
0117             p4Builder_.set(*strip);
0118             // Delete this candidate from future strips and move on to
0119             // the next potential candidate
0120             stripCand = cands.erase(stripCand);
0121           } else {
0122             // This candidate isn't compatabile - just move to the next candidate
0123             ++stripCand;
0124           }
0125         }
0126         // Update the vertex
0127         if (strip->daughterPtr(0).isNonnull())
0128           strip->setVertex(strip->daughterPtr(0)->vertex());
0129         output.emplace_back(strip.release());
0130       }
0131 
0132       // Check if we want to combine our strips
0133       if (combineStrips_ && output.size() > 1) {
0134         PiZeroVector stripCombinations;
0135         // Sort the output by descending pt
0136         std::sort(output.begin(), output.end(), [&](auto& arg1, auto& arg2) { return arg1->pt() > arg2->pt(); });
0137         // Get the end of interesting set of strips to try and combine
0138         PiZeroVector::const_iterator end_iter = takeNElements(output.begin(), output.end(), maxStrips_);
0139 
0140         // Look at all the combinations
0141         for (PiZeroVector::const_iterator firstIter = output.begin(); firstIter != end_iter - 1; ++firstIter) {
0142           for (PiZeroVector::const_iterator secondIter = firstIter + 1; secondIter != end_iter; ++secondIter) {
0143             auto const& first = *firstIter;
0144             auto const& second = *secondIter;
0145             Candidate::LorentzVector firstP4 = first->p4();
0146             Candidate::LorentzVector secondP4 = second->p4();
0147             // If we assume a certain mass for each strip apply it here.
0148             firstP4 = applyMassConstraint(firstP4, combinatoricStripMassHypo_);
0149             secondP4 = applyMassConstraint(secondP4, combinatoricStripMassHypo_);
0150             Candidate::LorentzVector totalP4 = firstP4 + secondP4;
0151             // Make our new combined strip
0152             auto combinedStrips =
0153                 std::make_unique<RecoTauPiZero>(0,
0154                                                 totalP4,
0155                                                 Candidate::Point(0, 0, 0),
0156                                                 //111, 10001, true, RecoTauPiZero::kCombinatoricStrips));
0157                                                 111,
0158                                                 10001,
0159                                                 true,
0160                                                 RecoTauPiZero::kUndefined);
0161 
0162             // Now loop over the strip members
0163             for (auto const& gamma : first->daughterPtrVector()) {
0164               combinedStrips->addDaughter(gamma);
0165             }
0166             for (auto const& gamma : second->daughterPtrVector()) {
0167               combinedStrips->addDaughter(gamma);
0168             }
0169             // Update the vertex
0170             if (combinedStrips->daughterPtr(0).isNonnull())
0171               combinedStrips->setVertex(combinedStrips->daughterPtr(0)->vertex());
0172             // Add to our collection of combined strips
0173             stripCombinations.emplace_back(combinedStrips.release());
0174           }
0175         }
0176         // When done doing all the combinations, add the combined strips to the output.
0177         std::move(stripCombinations.begin(), stripCombinations.end(), std::back_inserter(output));
0178       }
0179 
0180       return output;
0181     }
0182   }  // namespace tau
0183 }  // namespace reco
0184 
0185 #include "FWCore/Framework/interface/MakerMacros.h"
0186 DEFINE_EDM_PLUGIN(RecoTauPiZeroBuilderPluginFactory, reco::tau::RecoTauPiZeroStripPlugin, "RecoTauPiZeroStripPlugin");