File indexing completed on 2024-04-06 12:27:51
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 auto function = std::make_unique<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 qcuts_ = std::make_unique<RecoTauQualityCuts>(qcuts_pset);
0137
0138 inputParticleIds_ = pset.getParameter<std::vector<int> >("stripCandidatesParticleIds");
0139 const edm::ParameterSet& stripSize_eta_pset = pset.getParameterSet("stripEtaAssociationDistanceFunc");
0140 etaAssociationDistance_ = makeFunction("etaAssociationDistance", stripSize_eta_pset);
0141 const edm::ParameterSet& stripSize_phi_pset = pset.getParameterSet("stripPhiAssociationDistanceFunc");
0142 phiAssociationDistance_ = makeFunction("phiAssociationDistance", stripSize_phi_pset);
0143
0144 updateStripAfterEachDaughter_ = pset.getParameter<bool>("updateStripAfterEachDaughter");
0145 maxStripBuildIterations_ = pset.getParameter<int>("maxStripBuildIterations");
0146
0147 combineStrips_ = pset.getParameter<bool>("makeCombinatoricStrips");
0148 if (combineStrips_) {
0149 maxStrips_ = pset.getParameter<int>("maxInputStrips");
0150 combinatoricStripMassHypo_ = pset.getParameter<double>("stripMassWhenCombining");
0151 }
0152
0153 verbosity_ = pset.getParameter<int>("verbosity");
0154 }
0155 RecoTauPiZeroStripPlugin3::~RecoTauPiZeroStripPlugin3() {}
0156
0157
0158 void RecoTauPiZeroStripPlugin3::beginEvent() { vertexAssociator_.setEvent(*evt()); }
0159
0160 void RecoTauPiZeroStripPlugin3::addCandsToStrip(RecoTauPiZero& strip,
0161 CandPtrs& cands,
0162 const std::vector<bool>& candFlags,
0163 std::set<size_t>& candIdsCurrentStrip,
0164 bool& isCandAdded) const {
0165 if (verbosity_ >= 1) {
0166 edm::LogPrint("RecoTauPiZeroStripPlugin3") << "<RecoTauPiZeroStripPlugin3::addCandsToStrip>:";
0167 }
0168 size_t numCands = cands.size();
0169 for (size_t candId = 0; candId < numCands; ++candId) {
0170 if ((!candFlags[candId]) &&
0171 candIdsCurrentStrip.find(candId) == candIdsCurrentStrip.end()) {
0172 reco::CandidatePtr cand = cands[candId];
0173 double etaAssociationDistance_value =
0174 etaAssociationDistance_->Eval(strip.pt()) + etaAssociationDistance_->Eval(cand->pt());
0175 double phiAssociationDistance_value =
0176 phiAssociationDistance_->Eval(strip.pt()) + phiAssociationDistance_->Eval(cand->pt());
0177 if (fabs(strip.eta() - cand->eta()) <
0178 etaAssociationDistance_value &&
0179 reco::deltaPhi(strip.phi(), cand->phi()) < phiAssociationDistance_value) {
0180 if (verbosity_ >= 2) {
0181 edm::LogPrint("RecoTauPiZeroStripPlugin3")
0182 << "--> adding PFCand #" << candId << " (" << cand.id() << ":" << cand.key()
0183 << "): Et = " << cand->et() << ", eta = " << cand->eta() << ", phi = " << cand->phi();
0184 }
0185 strip.addDaughter(cand);
0186 if (updateStripAfterEachDaughter_)
0187 p4Builder_.set(strip);
0188 isCandAdded = true;
0189 candIdsCurrentStrip.insert(candId);
0190 }
0191 }
0192 }
0193 }
0194
0195 namespace {
0196 void markCandsInStrip(std::vector<bool>& candFlags, const std::set<size_t>& candIds) {
0197 for (std::set<size_t>::const_iterator candId = candIds.begin(); candId != candIds.end(); ++candId) {
0198 candFlags[*candId] = true;
0199 }
0200 }
0201
0202
0203
0204
0205 inline const reco::TrackBaseRef getTrack(const Candidate& cand) {
0206 const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0207 if (pfCandPtr) {
0208 if (pfCandPtr->trackRef().isNonnull())
0209 return reco::TrackBaseRef(pfCandPtr->trackRef());
0210 else if (pfCandPtr->gsfTrackRef().isNonnull())
0211 return reco::TrackBaseRef(pfCandPtr->gsfTrackRef());
0212 else
0213 return reco::TrackBaseRef();
0214 }
0215
0216 return reco::TrackBaseRef();
0217 }
0218 }
0219
0220 RecoTauPiZeroStripPlugin3::return_type RecoTauPiZeroStripPlugin3::operator()(const reco::Jet& jet) const {
0221 if (verbosity_ >= 1) {
0222 edm::LogPrint("RecoTauPiZeroStripPlugin3") << "<RecoTauPiZeroStripPlugin3::operator()>:";
0223 edm::LogPrint("RecoTauPiZeroStripPlugin3") << " minGammaEtStripSeed = " << minGammaEtStripSeed_;
0224 edm::LogPrint("RecoTauPiZeroStripPlugin3") << " minGammaEtStripAdd = " << minGammaEtStripAdd_;
0225 edm::LogPrint("RecoTauPiZeroStripPlugin3") << " minStripEt = " << minStripEt_;
0226 }
0227
0228 PiZeroVector output;
0229
0230
0231 qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
0232 CandPtrs candsVector = qcuts_->filterCandRefs(pfCandidates(jet, inputParticleIds_));
0233
0234
0235 CandPtrs seedCands;
0236 CandPtrs addCands;
0237 int idx = 0;
0238 for (CandPtrs::iterator cand = candsVector.begin(); cand != candsVector.end(); ++cand) {
0239 if (verbosity_ >= 1) {
0240 edm::LogPrint("RecoTauPiZeroStripPlugin3")
0241 << "PFGamma #" << idx << " (" << cand->id() << ":" << cand->key() << "): Et = " << (*cand)->et()
0242 << ", eta = " << (*cand)->eta() << ", phi = " << (*cand)->phi();
0243 }
0244 if ((*cand)->et() > minGammaEtStripSeed_) {
0245 if (verbosity_ >= 2) {
0246 edm::LogPrint("RecoTauPiZeroStripPlugin3") << "--> assigning seedCandId = " << seedCands.size();
0247 const reco::TrackBaseRef candTrack = getTrack(**cand);
0248 if (candTrack.isNonnull()) {
0249 edm::LogPrint("RecoTauPiZeroStripPlugin3")
0250 << "track: Pt = " << candTrack->pt() << " eta = " << candTrack->eta()
0251 << ", phi = " << candTrack->phi() << ", charge = " << candTrack->charge();
0252 edm::LogPrint("RecoTauPiZeroStripPlugin3")
0253 << " (dZ = " << candTrack->dz(vertexAssociator_.associatedVertex(jet)->position())
0254 << ", dXY = " << candTrack->dxy(vertexAssociator_.associatedVertex(jet)->position()) << ","
0255 << " numHits = " << candTrack->hitPattern().numberOfValidTrackerHits()
0256 << ", numPxlHits = " << candTrack->hitPattern().numberOfValidPixelHits() << ","
0257 << " chi2 = " << candTrack->normalizedChi2()
0258 << ", dPt/Pt = " << (candTrack->ptError() / candTrack->pt()) << ")";
0259 }
0260 }
0261 seedCands.push_back(*cand);
0262 } else if ((*cand)->et() > minGammaEtStripAdd_) {
0263 if (verbosity_ >= 2) {
0264 edm::LogPrint("RecoTauPiZeroStripPlugin3") << "--> assigning addCandId = " << addCands.size();
0265 }
0266 addCands.push_back(*cand);
0267 }
0268 ++idx;
0269 }
0270
0271 std::vector<bool> seedCandFlags(seedCands.size());
0272 std::vector<bool> addCandFlags(addCands.size());
0273
0274 std::set<size_t> seedCandIdsCurrentStrip;
0275 std::set<size_t> addCandIdsCurrentStrip;
0276
0277 size_t idxSeed = 0;
0278 while (idxSeed < seedCands.size()) {
0279 if (verbosity_ >= 2)
0280 edm::LogPrint("RecoTauPiZeroStripPlugin3") << "processing seed #" << idxSeed;
0281
0282 seedCandIdsCurrentStrip.clear();
0283 addCandIdsCurrentStrip.clear();
0284
0285 auto strip = std::make_unique<RecoTauPiZero>(*seedCands[idxSeed], RecoTauPiZero::kStrips);
0286 strip->addDaughter(seedCands[idxSeed]);
0287 seedCandIdsCurrentStrip.insert(idxSeed);
0288
0289 bool isCandAdded;
0290 int stripBuildIteration = 0;
0291 do {
0292 isCandAdded = false;
0293
0294
0295 addCandsToStrip(*strip, seedCands, seedCandFlags, seedCandIdsCurrentStrip, isCandAdded);
0296
0297 addCandsToStrip(*strip, addCands, addCandFlags, addCandIdsCurrentStrip, isCandAdded);
0298
0299 if (!updateStripAfterEachDaughter_)
0300 p4Builder_.set(*strip);
0301
0302 ++stripBuildIteration;
0303 } while (isCandAdded && (stripBuildIteration < maxStripBuildIterations_ || maxStripBuildIterations_ == -1));
0304
0305 if (strip->et() > minStripEt_) {
0306 if (verbosity_ >= 2)
0307 edm::LogPrint("RecoTauPiZeroStripPlugin3")
0308 << "Building strip: Et = " << strip->et() << ", eta = " << strip->eta() << ", phi = " << strip->phi();
0309
0310
0311 if (strip->daughterPtr(0).isNonnull())
0312 strip->setVertex(strip->daughterPtr(0)->vertex());
0313 output.push_back(std::move(strip));
0314
0315
0316 markCandsInStrip(seedCandFlags, seedCandIdsCurrentStrip);
0317 markCandsInStrip(addCandFlags, addCandIdsCurrentStrip);
0318 } else {
0319 if (verbosity_ >= 2)
0320 edm::LogPrint("RecoTauPiZeroStripPlugin3")
0321 << "Discarding strip: Et = " << strip->et() << ", eta = " << strip->eta() << ", phi = " << strip->phi();
0322 }
0323
0324 ++idxSeed;
0325 while (idxSeed < seedCands.size() && seedCandFlags[idxSeed]) {
0326 ++idxSeed;
0327 }
0328 }
0329
0330
0331 if (combineStrips_ && output.size() > 1) {
0332 PiZeroVector stripCombinations;
0333
0334 std::sort(output.begin(), output.end(), [&](auto& arg1, auto& arg2) { return arg1->pt() > arg2->pt(); });
0335
0336 PiZeroVector::const_iterator end_iter = takeNElements(output.begin(), output.end(), maxStrips_);
0337
0338
0339 for (PiZeroVector::const_iterator firstIter = output.begin(); firstIter != end_iter - 1; ++firstIter) {
0340 for (PiZeroVector::const_iterator secondIter = firstIter + 1; secondIter != end_iter; ++secondIter) {
0341 auto const& first = *firstIter;
0342 auto const& second = *secondIter;
0343 Candidate::LorentzVector firstP4 = first->p4();
0344 Candidate::LorentzVector secondP4 = second->p4();
0345
0346 firstP4 = applyMassConstraint(firstP4, combinatoricStripMassHypo_);
0347 secondP4 = applyMassConstraint(secondP4, combinatoricStripMassHypo_);
0348 Candidate::LorentzVector totalP4 = firstP4 + secondP4;
0349
0350 auto combinedStrips =
0351 std::make_unique<RecoTauPiZero>(0,
0352 totalP4,
0353 Candidate::Point(0, 0, 0),
0354
0355 111,
0356 10001,
0357 true,
0358 RecoTauPiZero::kUndefined);
0359
0360
0361 for (auto const& gamma : first->daughterPtrVector()) {
0362 combinedStrips->addDaughter(gamma);
0363 }
0364 for (auto const& gamma : second->daughterPtrVector()) {
0365 combinedStrips->addDaughter(gamma);
0366 }
0367
0368 if (combinedStrips->daughterPtr(0).isNonnull()) {
0369 combinedStrips->setVertex(combinedStrips->daughterPtr(0)->vertex());
0370 }
0371
0372
0373 stripCombinations.push_back(std::move(combinedStrips));
0374 }
0375 }
0376
0377
0378 std::move(stripCombinations.begin(), stripCombinations.end(), std::back_inserter(output));
0379 }
0380
0381
0382
0383 for (auto const& strip : output) {
0384 double bendCorrEta = 0.;
0385 double bendCorrPhi = 0.;
0386 double energySum = 0.;
0387 for (auto const& gamma : strip->daughterPtrVector()) {
0388 bendCorrEta += (gamma->energy() * etaAssociationDistance_->Eval(gamma->pt()));
0389 bendCorrPhi += (gamma->energy() * phiAssociationDistance_->Eval(gamma->pt()));
0390 energySum += gamma->energy();
0391 }
0392 if (energySum > 1.e-2) {
0393 bendCorrEta /= energySum;
0394 bendCorrPhi /= energySum;
0395 }
0396
0397 strip->setBendCorrEta(bendCorrEta);
0398 strip->setBendCorrPhi(bendCorrPhi);
0399 }
0400
0401 return output;
0402 }
0403 }
0404 }
0405
0406 #include "FWCore/Framework/interface/MakerMacros.h"
0407 DEFINE_EDM_PLUGIN(RecoTauPiZeroBuilderPluginFactory, reco::tau::RecoTauPiZeroStripPlugin3, "RecoTauPiZeroStripPlugin3");