RecoTau/plugins/RecoTauPiZeroStripPlugin2.cc

0001 /*
0002  * RecoTauPiZeroStripPlugin2
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  *                     Christian Veelken (LLR)
0011  *
0012  */
0013
0014 #include <algorithm>
0015 #include <functional>
0016 #include <memory>
0017
0018 #include "RecoTauTag/RecoTau/interface/RecoTauPiZeroPlugins.h"
0019
0020 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0021 #include "DataFormats/Candidate/interface/Candidate.h"
0022 #include "DataFormats/VertexReco/interface/Vertex.h"
0023 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0024 #include "DataFormats/TauReco/interface/RecoTauPiZero.h"
0025 #include "DataFormats/JetReco/interface/PFJet.h"
0026 #include "CommonTools/CandUtils/interface/AddFourMomenta.h"
0027 #include "DataFormats/Math/interface/deltaPhi.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029
0030 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0031 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0032 #include "RecoTauTag/RecoTau/interface/RecoTauVertexAssociator.h"
0033 #include "RecoTauTag/RecoTau/interface/CombinatoricGenerator.h"
0034
0035 //-------------------------------------------------------------------------------
0036 // CV: the following headers are needed only for debug print-out
0037 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0038 #include "DataFormats/TrackReco/interface/Track.h"
0039 //-------------------------------------------------------------------------------
0040
0041 namespace reco {
0042   namespace tau {
0043
0044     namespace {
0045       // Apply a hypothesis on the mass of the strips.
0046       math::XYZTLorentzVector applyMassConstraint(const math::XYZTLorentzVector& vec, double mass) {
0047         double factor = sqrt(vec.energy() * vec.energy() - mass * mass) / vec.P();
0048         return math::XYZTLorentzVector(vec.px() * factor, vec.py() * factor, vec.pz() * factor, vec.energy());
0049       }
0050     }  // namespace
0051
0052     class RecoTauPiZeroStripPlugin2 : public RecoTauPiZeroBuilderPlugin {
0053     public:
0054       explicit RecoTauPiZeroStripPlugin2(const edm::ParameterSet&, edm::ConsumesCollector&& iC);
0055       ~RecoTauPiZeroStripPlugin2() override;
0056       // Return type is unique_ptr<PiZeroVector>
0057       return_type operator()(const reco::Jet&) const override;
0058       // Hook to update PV information
0059       void beginEvent() override;
0060
0061     private:
0062       typedef std::vector<reco::CandidatePtr> CandPtrs;
0063       void addCandsToStrip(RecoTauPiZero&, CandPtrs&, const std::vector<bool>&, std::set<size_t>&, bool&) const;
0064
0065       RecoTauVertexAssociator vertexAssociator_;
0066
0067       std::unique_ptr<RecoTauQualityCuts> qcuts_;
0068       bool applyElecTrackQcuts_;
0069       double minGammaEtStripSeed_;
0070       double minGammaEtStripAdd_;
0071
0072       double minStripEt_;
0073
0074       std::vector<int> inputParticleIds_;  // type of candidates to clusterize
0075       double etaAssociationDistance_;      // size of strip clustering window in eta direction
0076       double phiAssociationDistance_;      // size of strip clustering window in phi direction
0077
0078       bool updateStripAfterEachDaughter_;
0079       int maxStripBuildIterations_;
0080
0081       // Parameters for build strip combinations
0082       bool combineStrips_;
0083       int maxStrips_;
0084       double combinatoricStripMassHypo_;
0085
0086       AddFourMomenta p4Builder_;
0087
0088       int verbosity_;
0089     };
0090
0091     RecoTauPiZeroStripPlugin2::RecoTauPiZeroStripPlugin2(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0092         : RecoTauPiZeroBuilderPlugin(pset, std::move(iC)),
0093           vertexAssociator_(pset.getParameter<edm::ParameterSet>("qualityCuts"), std::move(iC)),
0094           qcuts_(nullptr) {
0095       minGammaEtStripSeed_ = pset.getParameter<double>("minGammaEtStripSeed");
0096       minGammaEtStripAdd_ = pset.getParameter<double>("minGammaEtStripAdd");
0097
0098       minStripEt_ = pset.getParameter<double>("minStripEt");
0099
0100       edm::ParameterSet qcuts_pset = pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts");
0101       //-------------------------------------------------------------------------------
0102       // CV: disable track quality cuts for PFElectronsPFElectron
0103       //       (treat PFElectrons like PFGammas for the purpose of building eta-phi strips)
0104       applyElecTrackQcuts_ = pset.getParameter<bool>("applyElecTrackQcuts");
0105       if (!applyElecTrackQcuts_) {
0106         qcuts_pset.addParameter<double>("minTrackPt", std::min(minGammaEtStripSeed_, minGammaEtStripAdd_));
0107         qcuts_pset.addParameter<double>("maxTrackChi2", 1.e+9);
0108         qcuts_pset.addParameter<double>("maxTransverseImpactParameter", 1.e+9);
0109         qcuts_pset.addParameter<double>("maxDeltaZ", 1.e+9);
0110         qcuts_pset.addParameter<double>("minTrackVertexWeight", -1.);
0111         qcuts_pset.addParameter<unsigned>("minTrackPixelHits", 0);
0112         qcuts_pset.addParameter<unsigned>("minTrackHits", 0);
0113       }
0114       //-------------------------------------------------------------------------------
0115       qcuts_pset.addParameter<double>("minGammaEt", std::min(minGammaEtStripSeed_, minGammaEtStripAdd_));
0116       qcuts_ = std::make_unique<RecoTauQualityCuts>(qcuts_pset);
0117
0118       inputParticleIds_ = pset.getParameter<std::vector<int> >("stripCandidatesParticleIds");
0119       etaAssociationDistance_ = pset.getParameter<double>("stripEtaAssociationDistance");
0120       phiAssociationDistance_ = pset.getParameter<double>("stripPhiAssociationDistance");
0121
0122       updateStripAfterEachDaughter_ = pset.getParameter<bool>("updateStripAfterEachDaughter");
0123       maxStripBuildIterations_ = pset.getParameter<int>("maxStripBuildIterations");
0124
0125       combineStrips_ = pset.getParameter<bool>("makeCombinatoricStrips");
0126       if (combineStrips_) {
0127         maxStrips_ = pset.getParameter<int>("maxInputStrips");
0128         combinatoricStripMassHypo_ = pset.getParameter<double>("stripMassWhenCombining");
0129       }
0130
0131       verbosity_ = pset.getParameter<int>("verbosity");
0132     }
0133
0134     RecoTauPiZeroStripPlugin2::~RecoTauPiZeroStripPlugin2() {}
0135
0136     // Update the primary vertex
0137     void RecoTauPiZeroStripPlugin2::beginEvent() { vertexAssociator_.setEvent(*evt()); }
0138
0139     void RecoTauPiZeroStripPlugin2::addCandsToStrip(RecoTauPiZero& strip,
0140                                                     CandPtrs& cands,
0141                                                     const std::vector<bool>& candFlags,
0142                                                     std::set<size_t>& candIdsCurrentStrip,
0143                                                     bool& isCandAdded) const {
0144       if (verbosity_ >= 1) {
0145         edm::LogPrint("RecoTauPiZeroStripPlugin2") << "<RecoTauPiZeroStripPlugin2::addCandsToStrip>:";
0146       }
0147       size_t numCands = cands.size();
0148       for (size_t candId = 0; candId < numCands; ++candId) {
0149         if ((!candFlags[candId]) &&
0150             candIdsCurrentStrip.find(candId) == candIdsCurrentStrip.end()) {  // do not include same cand twice
0151           reco::CandidatePtr cand = cands[candId];
0152           if (fabs(strip.eta() - cand->eta()) <
0153                   etaAssociationDistance_ &&  // check if cand is within eta-phi window centered on strip
0154               fabs(strip.phi() - cand->phi()) < phiAssociationDistance_) {
0155             if (verbosity_ >= 2) {
0156               edm::LogPrint("RecoTauPiZeroStripPlugin2")
0157                   << "--> adding PFCand #" << candId << " (" << cand.id() << ":" << cand.key()
0158                   << "): Et = " << cand->et() << ", eta = " << cand->eta() << ", phi = " << cand->phi();
0159             }
0160             strip.addDaughter(cand);
0161             if (updateStripAfterEachDaughter_)
0162               p4Builder_.set(strip);
0163             isCandAdded = true;
0164             candIdsCurrentStrip.insert(candId);
0165           }
0166         }
0167       }
0168     }
0169
0170     namespace {
0171       void markCandsInStrip(std::vector<bool>& candFlags, const std::set<size_t>& candIds) {
0172         for (std::set<size_t>::const_iterator candId = candIds.begin(); candId != candIds.end(); ++candId) {
0173           candFlags[*candId] = true;
0174         }
0175       }
0176
0177       inline const reco::TrackBaseRef getTrack(const Candidate& cand) {
0178         const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0179         if (pfCandPtr) {
0180           if (pfCandPtr->trackRef().isNonnull())
0181             return reco::TrackBaseRef(pfCandPtr->trackRef());
0182           else if (pfCandPtr->gsfTrackRef().isNonnull())
0183             return reco::TrackBaseRef(pfCandPtr->gsfTrackRef());
0184           else
0185             return reco::TrackBaseRef();
0186         }
0187
0188         return reco::TrackBaseRef();
0189       }
0190     }  // namespace
0191
0192     RecoTauPiZeroStripPlugin2::return_type RecoTauPiZeroStripPlugin2::operator()(const reco::Jet& jet) const {
0193       if (verbosity_ >= 1) {
0194         edm::LogPrint("RecoTauPiZeroStripPlugin2") << "<RecoTauPiZeroStripPlugin2::operator()>:";
0195         edm::LogPrint("RecoTauPiZeroStripPlugin2") << " minGammaEtStripSeed = " << minGammaEtStripSeed_;
0196         edm::LogPrint("RecoTauPiZeroStripPlugin2") << " minGammaEtStripAdd = " << minGammaEtStripAdd_;
0197         edm::LogPrint("RecoTauPiZeroStripPlugin2") << " minStripEt = " << minStripEt_;
0198       }
0199
0200       PiZeroVector output;
0201
0202       // Get the candidates passing our quality cuts
0203       qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0204       CandPtrs candsVector = qcuts_->filterCandRefs(pfCandidates(jet, inputParticleIds_));
0205
0206       // Convert to stl::list to allow fast deletions
0207       CandPtrs seedCands;
0208       CandPtrs addCands;
0209       int idx = 0;
0210       for (CandPtrs::iterator cand = candsVector.begin(); cand != candsVector.end(); ++cand) {
0211         if (verbosity_ >= 1) {
0212           edm::LogPrint("RecoTauPiZeroStripPlugin2")
0213               << "PFGamma #" << idx << " (" << cand->id() << ":" << cand->key() << "): Et = " << (*cand)->et()
0214               << ", eta = " << (*cand)->eta() << ", phi = " << (*cand)->phi();
0215         }
0216         if ((*cand)->et() > minGammaEtStripSeed_) {
0217           if (verbosity_ >= 2) {
0218             edm::LogPrint("RecoTauPiZeroStripPlugin2") << "--> assigning seedCandId = " << seedCands.size();
0219             const reco::TrackBaseRef candTrack = getTrack(**cand);
0220             if (candTrack.isNonnull()) {
0221               edm::LogPrint("RecoTauPiZeroStripPlugin2")
0222                   << "track: Pt = " << candTrack->pt() << " eta = " << candTrack->eta()
0223                   << ", phi = " << candTrack->phi() << ", charge = " << candTrack->charge();
0224               edm::LogPrint("RecoTauPiZeroStripPlugin2")
0225                   << " (dZ = " << candTrack->dz(vertexAssociator_.associatedVertex(jet)->position())
0226                   << ", dXY = " << candTrack->dxy(vertexAssociator_.associatedVertex(jet)->position()) << ","
0227                   << " numHits = " << candTrack->hitPattern().numberOfValidTrackerHits()
0228                   << ", numPxlHits = " << candTrack->hitPattern().numberOfValidPixelHits() << ","
0229                   << " chi2 = " << candTrack->normalizedChi2()
0230                   << ", dPt/Pt = " << (candTrack->ptError() / candTrack->pt()) << ")";
0231             }
0232           }
0233           seedCands.push_back(*cand);
0234         } else if ((*cand)->et() > minGammaEtStripAdd_) {
0235           if (verbosity_ >= 2) {
0236             edm::LogPrint("RecoTauPiZeroStripPlugin2") << "--> assigning addCandId = " << addCands.size();
0237           }
0238           addCands.push_back(*cand);
0239         }
0240         ++idx;
0241       }
0242
0243       std::vector<bool> seedCandFlags(seedCands.size());  // true/false: seedCand is already/not yet included in strip
0244       std::vector<bool> addCandFlags(addCands.size());    // true/false: addCand  is already/not yet included in strip
0245
0246       std::set<size_t> seedCandIdsCurrentStrip;
0247       std::set<size_t> addCandIdsCurrentStrip;
0248
0249       size_t idxSeed = 0;
0250       while (idxSeed < seedCands.size()) {
0251         if (verbosity_ >= 2)
0252           edm::LogPrint("RecoTauPiZeroStripPlugin2") << "processing seed #" << idxSeed;
0253
0254         seedCandIdsCurrentStrip.clear();
0255         addCandIdsCurrentStrip.clear();
0256
0257         auto strip = std::make_unique<RecoTauPiZero>(*seedCands[idxSeed], RecoTauPiZero::kStrips);
0258         strip->addDaughter(seedCands[idxSeed]);
0259         seedCandIdsCurrentStrip.insert(idxSeed);
0260
0261         bool isCandAdded;
0262         int stripBuildIteration = 0;
0263         do {
0264           isCandAdded = false;
0265
0266           //if ( verbosity_ >= 2 ) edm::LogPrint("RecoTauPiZeroStripPlugin2") << " adding seedCands to strip..." ;
0267           addCandsToStrip(*strip, seedCands, seedCandFlags, seedCandIdsCurrentStrip, isCandAdded);
0268           //if ( verbosity_ >= 2 ) edm::LogPrint("RecoTauPiZeroStripPlugin2") << " adding addCands to strip..." ;
0269           addCandsToStrip(*strip, addCands, addCandFlags, addCandIdsCurrentStrip, isCandAdded);
0270
0271           if (!updateStripAfterEachDaughter_)
0272             p4Builder_.set(*strip);
0273
0274           ++stripBuildIteration;
0275         } while (isCandAdded && (stripBuildIteration < maxStripBuildIterations_ || maxStripBuildIterations_ == -1));
0276
0277         if (strip->et() > minStripEt_) {  // strip passed Et cuts, add it to the event
0278           if (verbosity_ >= 2)
0279             edm::LogPrint("RecoTauPiZeroStripPlugin2")
0280                 << "Building strip: Et = " << strip->et() << ", eta = " << strip->eta() << ", phi = " << strip->phi();
0281
0282           // Update the vertex
0283           if (strip->daughterPtr(0).isNonnull())
0284             strip->setVertex(strip->daughterPtr(0)->vertex());
0285           output.push_back(std::move(strip));
0286
0287           // Mark daughters as being part of this strip
0288           markCandsInStrip(seedCandFlags, seedCandIdsCurrentStrip);
0289           markCandsInStrip(addCandFlags, addCandIdsCurrentStrip);
0290         } else {  // strip failed Et cuts, just skip it
0291           if (verbosity_ >= 2)
0292             edm::LogPrint("RecoTauPiZeroStripPlugin2")
0293                 << "Discarding strip: Et = " << strip->et() << ", eta = " << strip->eta() << ", phi = " << strip->phi();
0294         }
0295
0296         ++idxSeed;
0297         while (idxSeed < seedCands.size() && seedCandFlags[idxSeed]) {
0298           ++idxSeed;  // fast-forward to next seed cand not yet included in any strip
0299         }
0300       }
0301
0302       // Check if we want to combine our strips
0303       if (combineStrips_ && output.size() > 1) {
0304         PiZeroVector stripCombinations;
0305         // Sort the output by descending pt
0306         std::sort(output.begin(), output.end(), [&](auto& arg1, auto& arg2) { return arg1->pt() > arg2->pt(); });
0307         // Get the end of interesting set of strips to try and combine
0308         PiZeroVector::const_iterator end_iter = takeNElements(output.begin(), output.end(), maxStrips_);
0309
0310         // Look at all the combinations
0311         for (PiZeroVector::const_iterator firstIter = output.begin(); firstIter != end_iter - 1; ++firstIter) {
0312           for (PiZeroVector::const_iterator secondIter = firstIter + 1; secondIter != end_iter; ++secondIter) {
0313             auto const& first = *firstIter;
0314             auto const& second = *secondIter;
0315             Candidate::LorentzVector firstP4 = first->p4();
0316             Candidate::LorentzVector secondP4 = second->p4();
0317             // If we assume a certain mass for each strip apply it here.
0318             firstP4 = applyMassConstraint(firstP4, combinatoricStripMassHypo_);
0319             secondP4 = applyMassConstraint(secondP4, combinatoricStripMassHypo_);
0320             Candidate::LorentzVector totalP4 = firstP4 + secondP4;
0321             // Make our new combined strip
0322             auto combinedStrips =
0323                 std::make_unique<RecoTauPiZero>(0,
0324                                                 totalP4,
0325                                                 Candidate::Point(0, 0, 0),
0326                                                 //111, 10001, true, RecoTauPiZero::kCombinatoricStrips));
0327                                                 111,
0328                                                 10001,
0329                                                 true,
0330                                                 RecoTauPiZero::kUndefined);
0331
0332             // Now loop over the strip members
0333             for (auto const& gamma : first->daughterPtrVector()) {
0334               combinedStrips->addDaughter(gamma);
0335             }
0336             for (auto const& gamma : second->daughterPtrVector()) {
0337               combinedStrips->addDaughter(gamma);
0338             }
0339             // Update the vertex
0340             if (combinedStrips->daughterPtr(0).isNonnull())
0341               combinedStrips->setVertex(combinedStrips->daughterPtr(0)->vertex());
0342             // Add to our collection of combined strips
0343             stripCombinations.push_back(std::move(combinedStrips));
0344           }
0345         }
0346         // When done doing all the combinations, add the combined strips to the
0347         // output.
0348         std::move(stripCombinations.begin(), stripCombinations.end(), std::back_inserter(output));
0349       }
0350
0351       return output;
0352     }
0353   }  // namespace tau
0354 }  // namespace reco
0355
0356 #include "FWCore/Framework/interface/MakerMacros.h"
0357 DEFINE_EDM_PLUGIN(RecoTauPiZeroBuilderPluginFactory, reco::tau::RecoTauPiZeroStripPlugin2, "RecoTauPiZeroStripPlugin2");