File indexing completed on 2021-10-25 04:55:30
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
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
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
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 }
0054
0055 class RecoTauPiZeroStripPlugin3 : public RecoTauPiZeroBuilderPlugin {
0056 public:
0057 explicit RecoTauPiZeroStripPlugin3(const edm::ParameterSet&, edm::ConsumesCollector&& iC);
0058 ~RecoTauPiZeroStripPlugin3() override;
0059
0060 return_type operator()(const reco::Jet&) const override;
0061
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_;
0078 std::unique_ptr<const TFormula> etaAssociationDistance_;
0079 std::unique_ptr<const TFormula> phiAssociationDistance_;
0080
0081 bool updateStripAfterEachDaughter_;
0082 int maxStripBuildIterations_;
0083
0084
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 }
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
0123
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
0137
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
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()) {
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 &&
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
0206
0207
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 }
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
0234 qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0235 CandPtrs candsVector = qcuts_->filterCandRefs(pfCandidates(jet, inputParticleIds_));
0236
0237
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());
0275 std::vector<bool> addCandFlags(addCands.size());
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
0298 addCandsToStrip(*strip, seedCands, seedCandFlags, seedCandIdsCurrentStrip, isCandAdded);
0299
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_) {
0309 if (verbosity_ >= 2)
0310 edm::LogPrint("RecoTauPiZeroStripPlugin3")
0311 << "Building strip: Et = " << strip->et() << ", eta = " << strip->eta() << ", phi = " << strip->phi();
0312
0313
0314 if (strip->daughterPtr(0).isNonnull())
0315 strip->setVertex(strip->daughterPtr(0)->vertex());
0316 output.push_back(std::move(strip));
0317
0318
0319 markCandsInStrip(seedCandFlags, seedCandIdsCurrentStrip);
0320 markCandsInStrip(addCandFlags, addCandIdsCurrentStrip);
0321 } else {
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;
0330 }
0331 }
0332
0333
0334 if (combineStrips_ && output.size() > 1) {
0335 PiZeroVector stripCombinations;
0336
0337 std::sort(output.begin(), output.end(), [&](auto& arg1, auto& arg2) { return arg1->pt() > arg2->pt(); });
0338
0339 PiZeroVector::const_iterator end_iter = takeNElements(output.begin(), output.end(), maxStrips_);
0340
0341
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
0349 firstP4 = applyMassConstraint(firstP4, combinatoricStripMassHypo_);
0350 secondP4 = applyMassConstraint(secondP4, combinatoricStripMassHypo_);
0351 Candidate::LorentzVector totalP4 = firstP4 + secondP4;
0352
0353 auto combinedStrips =
0354 std::make_unique<RecoTauPiZero>(0,
0355 totalP4,
0356 Candidate::Point(0, 0, 0),
0357
0358 111,
0359 10001,
0360 true,
0361 RecoTauPiZero::kUndefined);
0362
0363
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
0371 if (combinedStrips->daughterPtr(0).isNonnull()) {
0372 combinedStrips->setVertex(combinedStrips->daughterPtr(0)->vertex());
0373 }
0374
0375
0376 stripCombinations.push_back(std::move(combinedStrips));
0377 }
0378 }
0379
0380
0381 std::move(stripCombinations.begin(), stripCombinations.end(), std::back_inserter(output));
0382 }
0383
0384
0385
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
0400 strip->setBendCorrEta(bendCorrEta);
0401 strip->setBendCorrPhi(bendCorrPhi);
0402 }
0403
0404 return output;
0405 }
0406 }
0407 }
0408
0409 #include "FWCore/Framework/interface/MakerMacros.h"
0410 DEFINE_EDM_PLUGIN(RecoTauPiZeroBuilderPluginFactory, reco::tau::RecoTauPiZeroStripPlugin3, "RecoTauPiZeroStripPlugin3");