RecoTau/plugins/RecoTauPiZeroStripPlugin3.cc

0001 /*
0002  * RecoTauPiZeroStripPlugin3
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/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0021 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.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 #include "TString.h"
0042 #include "TFormula.h"
0043
0044 namespace reco {
0045   namespace tau {
0046
0047     namespace {
0048       // Apply a hypothesis on the mass of the strips.
0049       math::XYZTLorentzVector applyMassConstraint(const math::XYZTLorentzVector& vec, double mass) {
0050         double factor = sqrt(vec.energy() * vec.energy() - mass * mass) / vec.P();
0051         return math::XYZTLorentzVector(vec.px() * factor, vec.py() * factor, vec.pz() * factor, vec.energy());
0052       }
0053     }  // namespace
0054
0055     class RecoTauPiZeroStripPlugin3 : public RecoTauPiZeroBuilderPlugin {
0056     public:
0057       explicit RecoTauPiZeroStripPlugin3(const edm::ParameterSet&, edm::ConsumesCollector&& iC);
0058       ~RecoTauPiZeroStripPlugin3() override;
0059       // Return type is unique_ptr<PiZeroVector>
0060       return_type operator()(const reco::Jet&) const override;
0061       // Hook to update PV information
0062       void beginEvent() override;
0063
0064     private:
0065       typedef std::vector<reco::CandidatePtr> CandPtrs;
0066       void addCandsToStrip(RecoTauPiZero&, CandPtrs&, const std::vector<bool>&, std::set<size_t>&, bool&) const;
0067
0068       RecoTauVertexAssociator vertexAssociator_;
0069
0070       std::unique_ptr<RecoTauQualityCuts> qcuts_;
0071       bool applyElecTrackQcuts_;
0072       double minGammaEtStripSeed_;
0073       double minGammaEtStripAdd_;
0074
0075       double minStripEt_;
0076
0077       std::vector<int> inputParticleIds_;                       // type of candidates to clusterize
0078       std::unique_ptr<const TFormula> etaAssociationDistance_;  // size of strip clustering window in eta direction
0079       std::unique_ptr<const TFormula> phiAssociationDistance_;  // size of strip clustering window in phi direction
0080
0081       bool updateStripAfterEachDaughter_;
0082       int maxStripBuildIterations_;
0083
0084       // Parameters for build strip combinations
0085       bool combineStrips_;
0086       int maxStrips_;
0087       double combinatoricStripMassHypo_;
0088
0089       AddFourMomenta p4Builder_;
0090
0091       int verbosity_;
0092     };
0093
0094     namespace {
0095       std::unique_ptr<TFormula> makeFunction(const std::string& functionName, const edm::ParameterSet& pset) {
0096         TString formula = pset.getParameter<std::string>("function");
0097         formula = formula.ReplaceAll("pT", "x");
0098         std::unique_ptr<TFormula> function(new TFormula(functionName.data(), formula.Data()));
0099         int numParameter = function->GetNpar();
0100         for (int idxParameter = 0; idxParameter < numParameter; ++idxParameter) {
0101           std::string parameterName = Form("par%i", idxParameter);
0102           double parameter = pset.getParameter<double>(parameterName);
0103           function->SetParameter(idxParameter, parameter);
0104         }
0105         return function;
0106       }
0107     }  // namespace
0108
0109     RecoTauPiZeroStripPlugin3::RecoTauPiZeroStripPlugin3(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC)
0110         : RecoTauPiZeroBuilderPlugin(pset, std::move(iC)),
0111           vertexAssociator_(pset.getParameter<edm::ParameterSet>("qualityCuts"), std::move(iC)),
0112           qcuts_(nullptr),
0113           etaAssociationDistance_(nullptr),
0114           phiAssociationDistance_(nullptr) {
0115       minGammaEtStripSeed_ = pset.getParameter<double>("minGammaEtStripSeed");
0116       minGammaEtStripAdd_ = pset.getParameter<double>("minGammaEtStripAdd");
0117
0118       minStripEt_ = pset.getParameter<double>("minStripEt");
0119
0120       edm::ParameterSet qcuts_pset = pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts");
0121       //-------------------------------------------------------------------------------
0122       // CV: disable track quality cuts for PFElectronsPFElectron
0123       //       (treat PFElectrons like PFGammas for the purpose of building eta-phi strips)
0124       applyElecTrackQcuts_ = pset.getParameter<bool>("applyElecTrackQcuts");
0125       if (!applyElecTrackQcuts_) {
0126         qcuts_pset.addParameter<double>("minTrackPt", std::min(minGammaEtStripSeed_, minGammaEtStripAdd_));
0127         qcuts_pset.addParameter<double>("maxTrackChi2", 1.e+9);
0128         qcuts_pset.addParameter<double>("maxTransverseImpactParameter", 1.e+9);
0129         qcuts_pset.addParameter<double>("maxDeltaZ", 1.e+9);
0130         qcuts_pset.addParameter<double>("minTrackVertexWeight", -1.);
0131         qcuts_pset.addParameter<unsigned>("minTrackPixelHits", 0);
0132         qcuts_pset.addParameter<unsigned>("minTrackHits", 0);
0133       }
0134       //-------------------------------------------------------------------------------
0135       qcuts_pset.addParameter<double>("minGammaEt", std::min(minGammaEtStripSeed_, minGammaEtStripAdd_));
0136       //qcuts_ = new RecoTauQualityCuts(qcuts_pset);
0137       //std::unique_ptr<RecoTauQualityCuts> qcuts_(new RecoTauQualityCuts(qcuts_pset));
0138
0139       qcuts_ = std::make_unique<RecoTauQualityCuts>(qcuts_pset);
0140
0141       inputParticleIds_ = pset.getParameter<std::vector<int> >("stripCandidatesParticleIds");
0142       const edm::ParameterSet& stripSize_eta_pset = pset.getParameterSet("stripEtaAssociationDistanceFunc");
0143       etaAssociationDistance_ = makeFunction("etaAssociationDistance", stripSize_eta_pset);
0144       const edm::ParameterSet& stripSize_phi_pset = pset.getParameterSet("stripPhiAssociationDistanceFunc");
0145       phiAssociationDistance_ = makeFunction("phiAssociationDistance", stripSize_phi_pset);
0146
0147       updateStripAfterEachDaughter_ = pset.getParameter<bool>("updateStripAfterEachDaughter");
0148       maxStripBuildIterations_ = pset.getParameter<int>("maxStripBuildIterations");
0149
0150       combineStrips_ = pset.getParameter<bool>("makeCombinatoricStrips");
0151       if (combineStrips_) {
0152         maxStrips_ = pset.getParameter<int>("maxInputStrips");
0153         combinatoricStripMassHypo_ = pset.getParameter<double>("stripMassWhenCombining");
0154       }
0155
0156       verbosity_ = pset.getParameter<int>("verbosity");
0157     }
0158     RecoTauPiZeroStripPlugin3::~RecoTauPiZeroStripPlugin3() {}
0159
0160     // Update the primary vertex
0161     void RecoTauPiZeroStripPlugin3::beginEvent() { vertexAssociator_.setEvent(*evt()); }
0162
0163     void RecoTauPiZeroStripPlugin3::addCandsToStrip(RecoTauPiZero& strip,
0164                                                     CandPtrs& cands,
0165                                                     const std::vector<bool>& candFlags,
0166                                                     std::set<size_t>& candIdsCurrentStrip,
0167                                                     bool& isCandAdded) const {
0168       if (verbosity_ >= 1) {
0169         edm::LogPrint("RecoTauPiZeroStripPlugin3") << "<RecoTauPiZeroStripPlugin3::addCandsToStrip>:";
0170       }
0171       size_t numCands = cands.size();
0172       for (size_t candId = 0; candId < numCands; ++candId) {
0173         if ((!candFlags[candId]) &&
0174             candIdsCurrentStrip.find(candId) == candIdsCurrentStrip.end()) {  // do not include same cand twice
0175           reco::CandidatePtr cand = cands[candId];
0176           double etaAssociationDistance_value =
0177               etaAssociationDistance_->Eval(strip.pt()) + etaAssociationDistance_->Eval(cand->pt());
0178           double phiAssociationDistance_value =
0179               phiAssociationDistance_->Eval(strip.pt()) + phiAssociationDistance_->Eval(cand->pt());
0180           if (fabs(strip.eta() - cand->eta()) <
0181                   etaAssociationDistance_value &&  // check if cand is within eta-phi window centered on strip
0182               reco::deltaPhi(strip.phi(), cand->phi()) < phiAssociationDistance_value) {
0183             if (verbosity_ >= 2) {
0184               edm::LogPrint("RecoTauPiZeroStripPlugin3")
0185                   << "--> adding PFCand #" << candId << " (" << cand.id() << ":" << cand.key()
0186                   << "): Et = " << cand->et() << ", eta = " << cand->eta() << ", phi = " << cand->phi();
0187             }
0188             strip.addDaughter(cand);
0189             if (updateStripAfterEachDaughter_)
0190               p4Builder_.set(strip);
0191             isCandAdded = true;
0192             candIdsCurrentStrip.insert(candId);
0193           }
0194         }
0195       }
0196     }
0197
0198     namespace {
0199       void markCandsInStrip(std::vector<bool>& candFlags, const std::set<size_t>& candIds) {
0200         for (std::set<size_t>::const_iterator candId = candIds.begin(); candId != candIds.end(); ++candId) {
0201           candFlags[*candId] = true;
0202         }
0203       }
0204
0205       // The following method has not been adapted to work with PackedCandidates,
0206       // however, it is only used for printing/debugging and can be adapted when
0207       // needed.
0208       inline const reco::TrackBaseRef getTrack(const Candidate& cand) {
0209         const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0210         if (pfCandPtr) {
0211           if (pfCandPtr->trackRef().isNonnull())
0212             return reco::TrackBaseRef(pfCandPtr->trackRef());
0213           else if (pfCandPtr->gsfTrackRef().isNonnull())
0214             return reco::TrackBaseRef(pfCandPtr->gsfTrackRef());
0215           else
0216             return reco::TrackBaseRef();
0217         }
0218
0219         return reco::TrackBaseRef();
0220       }
0221     }  // namespace
0222
0223     RecoTauPiZeroStripPlugin3::return_type RecoTauPiZeroStripPlugin3::operator()(const reco::Jet& jet) const {
0224       if (verbosity_ >= 1) {
0225         edm::LogPrint("RecoTauPiZeroStripPlugin3") << "<RecoTauPiZeroStripPlugin3::operator()>:";
0226         edm::LogPrint("RecoTauPiZeroStripPlugin3") << " minGammaEtStripSeed = " << minGammaEtStripSeed_;
0227         edm::LogPrint("RecoTauPiZeroStripPlugin3") << " minGammaEtStripAdd = " << minGammaEtStripAdd_;
0228         edm::LogPrint("RecoTauPiZeroStripPlugin3") << " minStripEt = " << minStripEt_;
0229       }
0230
0231       PiZeroVector output;
0232
0233       // Get the candidates passing our quality cuts
0234       qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0235       CandPtrs candsVector = qcuts_->filterCandRefs(pfCandidates(jet, inputParticleIds_));
0236
0237       // Convert to stl::list to allow fast deletions
0238       CandPtrs seedCands;
0239       CandPtrs addCands;
0240       int idx = 0;
0241       for (CandPtrs::iterator cand = candsVector.begin(); cand != candsVector.end(); ++cand) {
0242         if (verbosity_ >= 1) {
0243           edm::LogPrint("RecoTauPiZeroStripPlugin3")
0244               << "PFGamma #" << idx << " (" << cand->id() << ":" << cand->key() << "): Et = " << (*cand)->et()
0245               << ", eta = " << (*cand)->eta() << ", phi = " << (*cand)->phi();
0246         }
0247         if ((*cand)->et() > minGammaEtStripSeed_) {
0248           if (verbosity_ >= 2) {
0249             edm::LogPrint("RecoTauPiZeroStripPlugin3") << "--> assigning seedCandId = " << seedCands.size();
0250             const reco::TrackBaseRef candTrack = getTrack(**cand);
0251             if (candTrack.isNonnull()) {
0252               edm::LogPrint("RecoTauPiZeroStripPlugin3")
0253                   << "track: Pt = " << candTrack->pt() << " eta = " << candTrack->eta()
0254                   << ", phi = " << candTrack->phi() << ", charge = " << candTrack->charge();
0255               edm::LogPrint("RecoTauPiZeroStripPlugin3")
0256                   << " (dZ = " << candTrack->dz(vertexAssociator_.associatedVertex(jet)->position())
0257                   << ", dXY = " << candTrack->dxy(vertexAssociator_.associatedVertex(jet)->position()) << ","
0258                   << " numHits = " << candTrack->hitPattern().numberOfValidTrackerHits()
0259                   << ", numPxlHits = " << candTrack->hitPattern().numberOfValidPixelHits() << ","
0260                   << " chi2 = " << candTrack->normalizedChi2()
0261                   << ", dPt/Pt = " << (candTrack->ptError() / candTrack->pt()) << ")";
0262             }
0263           }
0264           seedCands.push_back(*cand);
0265         } else if ((*cand)->et() > minGammaEtStripAdd_) {
0266           if (verbosity_ >= 2) {
0267             edm::LogPrint("RecoTauPiZeroStripPlugin3") << "--> assigning addCandId = " << addCands.size();
0268           }
0269           addCands.push_back(*cand);
0270         }
0271         ++idx;
0272       }
0273
0274       std::vector<bool> seedCandFlags(seedCands.size());  // true/false: seedCand is already/not yet included in strip
0275       std::vector<bool> addCandFlags(addCands.size());    // true/false: addCand  is already/not yet included in strip
0276
0277       std::set<size_t> seedCandIdsCurrentStrip;
0278       std::set<size_t> addCandIdsCurrentStrip;
0279
0280       size_t idxSeed = 0;
0281       while (idxSeed < seedCands.size()) {
0282         if (verbosity_ >= 2)
0283           edm::LogPrint("RecoTauPiZeroStripPlugin3") << "processing seed #" << idxSeed;
0284
0285         seedCandIdsCurrentStrip.clear();
0286         addCandIdsCurrentStrip.clear();
0287
0288         std::unique_ptr<RecoTauPiZero> strip(new RecoTauPiZero(*seedCands[idxSeed], RecoTauPiZero::kStrips));
0289         strip->addDaughter(seedCands[idxSeed]);
0290         seedCandIdsCurrentStrip.insert(idxSeed);
0291
0292         bool isCandAdded;
0293         int stripBuildIteration = 0;
0294         do {
0295           isCandAdded = false;
0296
0297           //if ( verbosity_ >= 2 ) edm::LogPrint("RecoTauPiZeroStripPlugin3") << " adding seedCands to strip..." ;
0298           addCandsToStrip(*strip, seedCands, seedCandFlags, seedCandIdsCurrentStrip, isCandAdded);
0299           //if ( verbosity_ >= 2 ) edm::LogPrint("RecoTauPiZeroStripPlugin3") << " adding addCands to strip..." ;
0300           addCandsToStrip(*strip, addCands, addCandFlags, addCandIdsCurrentStrip, isCandAdded);
0301
0302           if (!updateStripAfterEachDaughter_)
0303             p4Builder_.set(*strip);
0304
0305           ++stripBuildIteration;
0306         } while (isCandAdded && (stripBuildIteration < maxStripBuildIterations_ || maxStripBuildIterations_ == -1));
0307
0308         if (strip->et() > minStripEt_) {  // strip passed Et cuts, add it to the event
0309           if (verbosity_ >= 2)
0310             edm::LogPrint("RecoTauPiZeroStripPlugin3")
0311                 << "Building strip: Et = " << strip->et() << ", eta = " << strip->eta() << ", phi = " << strip->phi();
0312
0313           // Update the vertex
0314           if (strip->daughterPtr(0).isNonnull())
0315             strip->setVertex(strip->daughterPtr(0)->vertex());
0316           output.push_back(std::move(strip));
0317
0318           // Mark daughters as being part of this strip
0319           markCandsInStrip(seedCandFlags, seedCandIdsCurrentStrip);
0320           markCandsInStrip(addCandFlags, addCandIdsCurrentStrip);
0321         } else {  // strip failed Et cuts, just skip it
0322           if (verbosity_ >= 2)
0323             edm::LogPrint("RecoTauPiZeroStripPlugin3")
0324                 << "Discarding strip: Et = " << strip->et() << ", eta = " << strip->eta() << ", phi = " << strip->phi();
0325         }
0326
0327         ++idxSeed;
0328         while (idxSeed < seedCands.size() && seedCandFlags[idxSeed]) {
0329           ++idxSeed;  // fast-forward to next seed cand not yet included in any strip
0330         }
0331       }
0332
0333       // Check if we want to combine our strips
0334       if (combineStrips_ && output.size() > 1) {
0335         PiZeroVector stripCombinations;
0336         // Sort the output by descending pt
0337         std::sort(output.begin(), output.end(), [&](auto& arg1, auto& arg2) { return arg1->pt() > arg2->pt(); });
0338         // Get the end of interesting set of strips to try and combine
0339         PiZeroVector::const_iterator end_iter = takeNElements(output.begin(), output.end(), maxStrips_);
0340
0341         // Look at all the combinations
0342         for (PiZeroVector::const_iterator firstIter = output.begin(); firstIter != end_iter - 1; ++firstIter) {
0343           for (PiZeroVector::const_iterator secondIter = firstIter + 1; secondIter != end_iter; ++secondIter) {
0344             auto const& first = *firstIter;
0345             auto const& second = *secondIter;
0346             Candidate::LorentzVector firstP4 = first->p4();
0347             Candidate::LorentzVector secondP4 = second->p4();
0348             // If we assume a certain mass for each strip apply it here.
0349             firstP4 = applyMassConstraint(firstP4, combinatoricStripMassHypo_);
0350             secondP4 = applyMassConstraint(secondP4, combinatoricStripMassHypo_);
0351             Candidate::LorentzVector totalP4 = firstP4 + secondP4;
0352             // Make our new combined strip
0353             auto combinedStrips =
0354                 std::make_unique<RecoTauPiZero>(0,
0355                                                 totalP4,
0356                                                 Candidate::Point(0, 0, 0),
0357                                                 //111, 10001, true, RecoTauPiZero::kCombinatoricStrips));
0358                                                 111,
0359                                                 10001,
0360                                                 true,
0361                                                 RecoTauPiZero::kUndefined);
0362
0363             // Now loop over the strip members
0364             for (auto const& gamma : first->daughterPtrVector()) {
0365               combinedStrips->addDaughter(gamma);
0366             }
0367             for (auto const& gamma : second->daughterPtrVector()) {
0368               combinedStrips->addDaughter(gamma);
0369             }
0370             // Update the vertex
0371             if (combinedStrips->daughterPtr(0).isNonnull()) {
0372               combinedStrips->setVertex(combinedStrips->daughterPtr(0)->vertex());
0373             }
0374
0375             // Add to our collection of combined strips
0376             stripCombinations.push_back(std::move(combinedStrips));
0377           }
0378         }
0379         // When done doing all the combinations, add the combined strips to the
0380         // output.
0381         std::move(stripCombinations.begin(), stripCombinations.end(), std::back_inserter(output));
0382       }
0383
0384       // Compute correction to account for spread of photon energy in eta and phi
0385       // in case charged pions make nuclear interactions or photons convert within the tracking detector
0386       for (auto const& strip : output) {
0387         double bendCorrEta = 0.;
0388         double bendCorrPhi = 0.;
0389         double energySum = 0.;
0390         for (auto const& gamma : strip->daughterPtrVector()) {
0391           bendCorrEta += (gamma->energy() * etaAssociationDistance_->Eval(gamma->pt()));
0392           bendCorrPhi += (gamma->energy() * phiAssociationDistance_->Eval(gamma->pt()));
0393           energySum += gamma->energy();
0394         }
0395         if (energySum > 1.e-2) {
0396           bendCorrEta /= energySum;
0397           bendCorrPhi /= energySum;
0398         }
0399         //std::cout << "stripPt = " << strip->pt() << ": bendCorrEta = " << bendCorrEta << ", bendCorrPhi = " << bendCorrPhi << std::endl;
0400         strip->setBendCorrEta(bendCorrEta);
0401         strip->setBendCorrPhi(bendCorrPhi);
0402       }
0403
0404       return output;
0405     }
0406   }  // namespace tau
0407 }  // namespace reco
0408
0409 #include "FWCore/Framework/interface/MakerMacros.h"
0410 DEFINE_EDM_PLUGIN(RecoTauPiZeroBuilderPluginFactory, reco::tau::RecoTauPiZeroStripPlugin3, "RecoTauPiZeroStripPlugin3");