File indexing completed on 2024-04-06 12:27:54
0001 #include <memory>
0002
0003 #include "RecoTauTag/RecoTau/interface/RecoTauConstructor.h"
0004
0005 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
0006 #include "DataFormats/Common/interface/RefToPtr.h"
0007
0008 #include "DataFormats/Math/interface/deltaR.h"
0009
0010 #include "RecoTauTag/RecoTau/interface/pfRecoTauChargedHadronAuxFunctions.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 namespace reco::tau {
0014
0015 RecoTauConstructor::RecoTauConstructor(const JetBaseRef& jet,
0016 const edm::Handle<edm::View<reco::Candidate> >& pfCands,
0017 bool copyGammasFromPiZeros,
0018 const StringObjectFunction<reco::PFTau>* signalConeSize,
0019 double minAbsPhotonSumPt_insideSignalCone,
0020 double minRelPhotonSumPt_insideSignalCone,
0021 double minAbsPhotonSumPt_outsideSignalCone,
0022 double minRelPhotonSumPt_outsideSignalCone)
0023 : signalConeSize_(signalConeSize),
0024 minAbsPhotonSumPt_insideSignalCone_(minAbsPhotonSumPt_insideSignalCone),
0025 minRelPhotonSumPt_insideSignalCone_(minRelPhotonSumPt_insideSignalCone),
0026 minAbsPhotonSumPt_outsideSignalCone_(minAbsPhotonSumPt_outsideSignalCone),
0027 minRelPhotonSumPt_outsideSignalCone_(minRelPhotonSumPt_outsideSignalCone),
0028 pfCands_(pfCands) {
0029
0030 tau_ = std::make_unique<PFTau>();
0031
0032 copyGammas_ = copyGammasFromPiZeros;
0033
0034 collections_[std::make_pair(kSignal, kChargedHadron)] = &tau_->selectedSignalChargedHadrCands_;
0035 collections_[std::make_pair(kSignal, kGamma)] = &tau_->selectedSignalGammaCands_;
0036 collections_[std::make_pair(kSignal, kNeutralHadron)] = &tau_->selectedSignalNeutrHadrCands_;
0037 collections_[std::make_pair(kSignal, kAll)] = &tau_->selectedSignalCands_;
0038
0039 collections_[std::make_pair(kIsolation, kChargedHadron)] = &tau_->selectedIsolationChargedHadrCands_;
0040 collections_[std::make_pair(kIsolation, kGamma)] = &tau_->selectedIsolationGammaCands_;
0041 collections_[std::make_pair(kIsolation, kNeutralHadron)] = &tau_->selectedIsolationNeutrHadrCands_;
0042 collections_[std::make_pair(kIsolation, kAll)] = &tau_->selectedIsolationCands_;
0043
0044
0045
0046 for (auto const& colkey : collections_) {
0047
0048 sortedCollections_[colkey.first] = std::make_shared<SortedListPtr::element_type>();
0049 }
0050
0051 tau_->setjetRef(jet);
0052 }
0053
0054 void RecoTauConstructor::addPFCand(Region region, ParticleType type, const CandidatePtr& ptr, bool skipAddToP4) {
0055 LogDebug("TauConstructorAddPFCand") << " region = " << region << ", type = " << type << ": Pt = " << ptr->pt()
0056 << ", eta = " << ptr->eta() << ", phi = " << ptr->phi();
0057 if (region == kSignal) {
0058
0059
0060 if (((type != kGamma) || !copyGammas_) && !skipAddToP4) {
0061 LogDebug("TauConstructorAddPFCand") << "--> adding PFCand to tauP4.";
0062 p4_ += ptr->p4();
0063 }
0064 }
0065 getSortedCollection(region, type)->push_back(ptr);
0066
0067 getSortedCollection(region, kAll)->push_back(ptr);
0068 }
0069
0070 void RecoTauConstructor::reserve(Region region, ParticleType type, size_t size) {
0071 getSortedCollection(region, type)->reserve(size);
0072 getCollection(region, type)->reserve(size);
0073
0074 getSortedCollection(region, kAll)->reserve(getSortedCollection(region, kAll)->size() + size);
0075 getCollection(region, kAll)->reserve(getCollection(region, kAll)->size() + size);
0076 }
0077
0078 void RecoTauConstructor::reserveTauChargedHadron(Region region, size_t size) {
0079 if (region == kSignal) {
0080 tau_->signalTauChargedHadronCandidatesRestricted().reserve(size);
0081 tau_->selectedSignalChargedHadrCands_.reserve(size);
0082 } else {
0083 tau_->isolationTauChargedHadronCandidatesRestricted().reserve(size);
0084 tau_->selectedIsolationChargedHadrCands_.reserve(size);
0085 }
0086 }
0087
0088 namespace {
0089 void checkOverlap(const CandidatePtr& neutral, const std::vector<CandidatePtr>& pfGammas, bool& isUnique) {
0090 LogDebug("TauConstructorCheckOverlap") << " pfGammas: #entries = " << pfGammas.size();
0091 for (std::vector<CandidatePtr>::const_iterator pfGamma = pfGammas.begin(); pfGamma != pfGammas.end(); ++pfGamma) {
0092 LogDebug("TauConstructorCheckOverlap") << "pfGamma = " << pfGamma->id() << ":" << pfGamma->key();
0093 if ((*pfGamma).refCore() == neutral.refCore() && (*pfGamma).key() == neutral.key())
0094 isUnique = false;
0095 }
0096 }
0097
0098 void checkOverlap(const CandidatePtr& neutral, const std::vector<reco::RecoTauPiZero>& piZeros, bool& isUnique) {
0099 LogDebug("TauConstructorCheckOverlap") << " piZeros: #entries = " << piZeros.size();
0100 for (std::vector<reco::RecoTauPiZero>::const_iterator piZero = piZeros.begin(); piZero != piZeros.end();
0101 ++piZero) {
0102 size_t numPFGammas = piZero->numberOfDaughters();
0103 for (size_t iPFGamma = 0; iPFGamma < numPFGammas; ++iPFGamma) {
0104 reco::CandidatePtr pfGamma = piZero->daughterPtr(iPFGamma);
0105 LogDebug("TauConstructorCheckOverlap") << "pfGamma = " << pfGamma.id() << ":" << pfGamma.key();
0106 if (pfGamma.id() == neutral.id() && pfGamma.key() == neutral.key())
0107 isUnique = false;
0108 }
0109 }
0110 }
0111 }
0112
0113 void RecoTauConstructor::addTauChargedHadron(Region region, const PFRecoTauChargedHadron& chargedHadron) {
0114 LogDebug("TauConstructorAddChH") << " region = " << region << ": Pt = " << chargedHadron.pt()
0115 << ", eta = " << chargedHadron.eta() << ", phi = " << chargedHadron.phi();
0116
0117 const std::vector<CandidatePtr>& neutrals = chargedHadron.getNeutralPFCandidates();
0118 std::vector<CandidatePtr> neutrals_cleaned;
0119 for (std::vector<CandidatePtr>::const_iterator neutral = neutrals.begin(); neutral != neutrals.end(); ++neutral) {
0120 LogDebug("TauConstructorAddChH") << "neutral = " << neutral->id() << ":" << neutral->key();
0121 bool isUnique = true;
0122 if (copyGammas_)
0123 checkOverlap(*neutral, *getSortedCollection(kSignal, kGamma), isUnique);
0124 else
0125 checkOverlap(*neutral, tau_->signalPiZeroCandidatesRestricted(), isUnique);
0126 if (region == kIsolation) {
0127 if (copyGammas_)
0128 checkOverlap(*neutral, *getSortedCollection(kIsolation, kGamma), isUnique);
0129 else
0130 checkOverlap(*neutral, tau_->isolationPiZeroCandidatesRestricted(), isUnique);
0131 }
0132 LogDebug("TauConstructorAddChH") << "--> isUnique = " << isUnique;
0133 if (isUnique)
0134 neutrals_cleaned.push_back(*neutral);
0135 }
0136 PFRecoTauChargedHadron chargedHadron_cleaned = chargedHadron;
0137 if (neutrals_cleaned.size() != neutrals.size()) {
0138 chargedHadron_cleaned.neutralPFCandidates_ = neutrals_cleaned;
0139 setChargedHadronP4(chargedHadron_cleaned);
0140 }
0141 if (region == kSignal) {
0142 tau_->signalTauChargedHadronCandidatesRestricted().push_back(chargedHadron_cleaned);
0143 p4_ += chargedHadron_cleaned.p4();
0144 if (chargedHadron_cleaned.getChargedPFCandidate().isNonnull()) {
0145 addPFCand(kSignal, kChargedHadron, convertToPtr(chargedHadron_cleaned.getChargedPFCandidate()), true);
0146 }
0147 const std::vector<CandidatePtr>& neutrals = chargedHadron_cleaned.getNeutralPFCandidates();
0148 for (std::vector<CandidatePtr>::const_iterator neutral = neutrals.begin(); neutral != neutrals.end(); ++neutral) {
0149 if (std::abs((*neutral)->pdgId()) == 22)
0150 addPFCand(kSignal, kGamma, convertToPtr(*neutral), true);
0151 else if (std::abs((*neutral)->pdgId()) == 130)
0152 addPFCand(kSignal, kNeutralHadron, convertToPtr(*neutral), true);
0153 };
0154 } else {
0155 tau_->isolationTauChargedHadronCandidatesRestricted().push_back(chargedHadron_cleaned);
0156 if (chargedHadron_cleaned.getChargedPFCandidate().isNonnull()) {
0157 if (std::abs(chargedHadron_cleaned.getChargedPFCandidate()->pdgId()) == 211)
0158 addPFCand(kIsolation, kChargedHadron, convertToPtr(chargedHadron_cleaned.getChargedPFCandidate()));
0159 else if (std::abs(chargedHadron_cleaned.getChargedPFCandidate()->pdgId()) == 130)
0160 addPFCand(kIsolation, kNeutralHadron, convertToPtr(chargedHadron_cleaned.getChargedPFCandidate()));
0161 }
0162 const std::vector<CandidatePtr>& neutrals = chargedHadron_cleaned.getNeutralPFCandidates();
0163 for (std::vector<CandidatePtr>::const_iterator neutral = neutrals.begin(); neutral != neutrals.end(); ++neutral) {
0164 if (std::abs((*neutral)->pdgId()) == 22)
0165 addPFCand(kIsolation, kGamma, convertToPtr(*neutral));
0166 else if (std::abs((*neutral)->pdgId()) == 130)
0167 addPFCand(kIsolation, kNeutralHadron, convertToPtr(*neutral));
0168 };
0169 }
0170 }
0171
0172 void RecoTauConstructor::reservePiZero(Region region, size_t size) {
0173 if (region == kSignal) {
0174 tau_->signalPiZeroCandidatesRestricted().reserve(size);
0175
0176
0177 if (copyGammas_)
0178 reserve(kSignal, kGamma, 2 * size);
0179 } else {
0180 tau_->isolationPiZeroCandidatesRestricted().reserve(size);
0181 if (copyGammas_)
0182 reserve(kIsolation, kGamma, 2 * size);
0183 }
0184 }
0185
0186 void RecoTauConstructor::addPiZero(Region region, const RecoTauPiZero& piZero) {
0187 LogDebug("TauConstructorAddPi0") << " region = " << region << ": Pt = " << piZero.pt() << ", eta = " << piZero.eta()
0188 << ", phi = " << piZero.phi();
0189 if (region == kSignal) {
0190 tau_->signalPiZeroCandidatesRestricted().push_back(piZero);
0191
0192 if (copyGammas_) {
0193
0194
0195 p4_ += piZero.p4();
0196 addPFCands(kSignal, kGamma, piZero.daughterPtrVector().begin(), piZero.daughterPtrVector().end());
0197 }
0198 } else {
0199 tau_->isolationPiZeroCandidatesRestricted().push_back(piZero);
0200 if (copyGammas_) {
0201 addPFCands(kIsolation, kGamma, piZero.daughterPtrVector().begin(), piZero.daughterPtrVector().end());
0202 }
0203 }
0204 }
0205
0206 std::vector<CandidatePtr>* RecoTauConstructor::getCollection(Region region, ParticleType type) {
0207 return collections_[std::make_pair(region, type)];
0208 }
0209
0210 RecoTauConstructor::SortedListPtr RecoTauConstructor::getSortedCollection(Region region, ParticleType type) {
0211 return sortedCollections_[std::make_pair(region, type)];
0212 }
0213
0214
0215 CandidatePtr RecoTauConstructor::convertToPtr(const CandidatePtr& pfPtr) const { return pfPtr; }
0216
0217 namespace {
0218
0219 template <typename T1, typename T2>
0220 void checkMatchedProductIds(const T1& t1, const T2& t2) {
0221 if (t1.id() != t2.id()) {
0222 throw cms::Exception("MismatchedPFCandSrc")
0223 << "Error: the input tag"
0224 << " for the PF candidate collection provided to the RecoTauBuilder "
0225 << " does not match the one that was used to build the source jets."
0226 << " Please update the pfCandSrc paramters for the PFTau builders.";
0227 }
0228 }
0229 }
0230
0231
0232 CandidatePtr RecoTauConstructor::convertToPtr(const PFCandidatePtr& candPtr) const {
0233 if (candPtr.isNonnull()) {
0234 checkMatchedProductIds(candPtr, pfCands_);
0235 return CandidatePtr(pfCands_, candPtr.key());
0236 } else
0237 return PFCandidatePtr();
0238 }
0239
0240 namespace {
0241 template <typename T>
0242 bool ptDescending(const T& a, const T& b) {
0243 return a.pt() > b.pt();
0244 }
0245 template <typename T>
0246 bool ptDescendingPtr(const T& a, const T& b) {
0247 return a->pt() > b->pt();
0248 }
0249 }
0250
0251 void RecoTauConstructor::sortAndCopyIntoTau() {
0252
0253 std::sort(tau_->signalTauChargedHadronCandidatesRestricted().begin(),
0254 tau_->signalTauChargedHadronCandidatesRestricted().end(),
0255 ptDescending<PFRecoTauChargedHadron>);
0256 std::sort(tau_->isolationTauChargedHadronCandidatesRestricted().begin(),
0257 tau_->isolationTauChargedHadronCandidatesRestricted().end(),
0258 ptDescending<PFRecoTauChargedHadron>);
0259 std::sort(tau_->signalPiZeroCandidatesRestricted().begin(),
0260 tau_->signalPiZeroCandidatesRestricted().end(),
0261 ptDescending<RecoTauPiZero>);
0262 std::sort(tau_->isolationPiZeroCandidatesRestricted().begin(),
0263 tau_->isolationPiZeroCandidatesRestricted().end(),
0264 ptDescending<RecoTauPiZero>);
0265
0266
0267
0268 for (auto const& colkey : collections_) {
0269 SortedListPtr sortedCollection = sortedCollections_[colkey.first];
0270 std::sort(sortedCollection->begin(), sortedCollection->end(), ptDescendingPtr<CandidatePtr>);
0271
0272 for (std::vector<CandidatePtr>::const_iterator particle = sortedCollection->begin();
0273 particle != sortedCollection->end();
0274 ++particle) {
0275 colkey.second->push_back(*particle);
0276 }
0277 }
0278 }
0279
0280 namespace {
0281 PFTau::hadronicDecayMode calculateDecayMode(const reco::PFTau& tau,
0282 double dRsignalCone,
0283 double minAbsPhotonSumPt_insideSignalCone,
0284 double minRelPhotonSumPt_insideSignalCone,
0285 double minAbsPhotonSumPt_outsideSignalCone,
0286 double minRelPhotonSumPt_outsideSignalCone) {
0287 unsigned int nCharged = tau.signalTauChargedHadronCandidates().size();
0288
0289 if (!nCharged)
0290 return PFTau::kNull;
0291
0292 unsigned int nPiZeros = 0;
0293 const std::vector<RecoTauPiZero>& piZeros = tau.signalPiZeroCandidates();
0294 for (std::vector<RecoTauPiZero>::const_iterator piZero = piZeros.begin(); piZero != piZeros.end(); ++piZero) {
0295 double photonSumPt_insideSignalCone = 0.;
0296 double photonSumPt_outsideSignalCone = 0.;
0297 int numPhotons = piZero->numberOfDaughters();
0298 for (int idxPhoton = 0; idxPhoton < numPhotons; ++idxPhoton) {
0299 const reco::Candidate* photon = piZero->daughter(idxPhoton);
0300 double dR = deltaR(photon->p4(), tau.p4());
0301 if (dR < dRsignalCone) {
0302 photonSumPt_insideSignalCone += photon->pt();
0303 } else {
0304 photonSumPt_outsideSignalCone += photon->pt();
0305 }
0306 }
0307 if (photonSumPt_insideSignalCone > minAbsPhotonSumPt_insideSignalCone ||
0308 photonSumPt_insideSignalCone > (minRelPhotonSumPt_insideSignalCone * tau.pt()) ||
0309 photonSumPt_outsideSignalCone > minAbsPhotonSumPt_outsideSignalCone ||
0310 photonSumPt_outsideSignalCone > (minRelPhotonSumPt_outsideSignalCone * tau.pt()))
0311 ++nPiZeros;
0312 }
0313
0314
0315 const unsigned int maxPiZeros = PFTau::kOneProngNPiZero;
0316
0317
0318 unsigned int trackIndex = (nCharged - 1) * (maxPiZeros + 1);
0319
0320
0321 if (trackIndex >= PFTau::kRareDecayMode)
0322 return PFTau::kRareDecayMode;
0323
0324 if (nPiZeros > maxPiZeros)
0325 nPiZeros = maxPiZeros;
0326 return static_cast<PFTau::hadronicDecayMode>(trackIndex + nPiZeros);
0327 }
0328 }
0329
0330 std::unique_ptr<reco::PFTau> RecoTauConstructor::get(bool setupLeadingObjects) {
0331 LogDebug("TauConstructorGet") << "Start getting";
0332
0333
0334 sortAndCopyIntoTau();
0335
0336
0337
0338
0339
0340
0341
0342
0343
0344 unsigned int nCharged = 0;
0345 int charge = 0;
0346 double leadChargedHadronPt = 0.;
0347 int leadChargedHadronCharge = 0;
0348 for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron =
0349 tau_->signalTauChargedHadronCandidatesRestricted().begin();
0350 chargedHadron != tau_->signalTauChargedHadronCandidatesRestricted().end();
0351 ++chargedHadron) {
0352 if (chargedHadron->algoIs(PFRecoTauChargedHadron::kChargedPFCandidate) ||
0353 chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack)) {
0354 ++nCharged;
0355 charge += chargedHadron->charge();
0356 if (chargedHadron->pt() > leadChargedHadronPt) {
0357 leadChargedHadronPt = chargedHadron->pt();
0358 leadChargedHadronCharge = chargedHadron->charge();
0359 }
0360 }
0361 }
0362 if (nCharged == 3)
0363 tau_->setCharge(charge);
0364 else
0365 tau_->setCharge(leadChargedHadronCharge);
0366
0367
0368 tau_->setPdgId(tau_->charge() < 0 ? 15 : -15);
0369
0370
0371 tau_->setP4(p4_);
0372
0373
0374 double dRsignalCone = (signalConeSize_) ? (*signalConeSize_)(*tau_) : 0.5;
0375 tau_->setSignalConeSize(dRsignalCone);
0376 PFTau::hadronicDecayMode dm = calculateDecayMode(*tau_,
0377 dRsignalCone,
0378 minAbsPhotonSumPt_insideSignalCone_,
0379 minRelPhotonSumPt_insideSignalCone_,
0380 minAbsPhotonSumPt_outsideSignalCone_,
0381 minRelPhotonSumPt_outsideSignalCone_);
0382 tau_->setDecayMode(dm);
0383
0384 LogDebug("TauConstructorGet") << "Pt = " << tau_->pt() << ", eta = " << tau_->eta() << ", phi = " << tau_->phi()
0385 << ", mass = " << tau_->mass() << ", dm = " << tau_->decayMode();
0386
0387
0388 tau_->setisolationPFChargedHadrCandsPtSum(sumPFCandPt(getCollection(kIsolation, kChargedHadron)->begin(),
0389 getCollection(kIsolation, kChargedHadron)->end()));
0390
0391
0392 tau_->setisolationPFGammaCandsEtSum(
0393 sumPFCandPt(getCollection(kIsolation, kGamma)->begin(), getCollection(kIsolation, kGamma)->end()));
0394
0395
0396 tau_->setemFraction(sumPFCandPt(getCollection(kSignal, kGamma)->begin(), getCollection(kSignal, kGamma)->end()) /
0397 tau_->pt());
0398
0399 if (setupLeadingObjects) {
0400 typedef std::vector<CandidatePtr>::const_iterator Iter;
0401
0402 Iter leadingCand = leadCand(getCollection(kSignal, kAll)->begin(), getCollection(kSignal, kAll)->end());
0403
0404 if (leadingCand != getCollection(kSignal, kAll)->end())
0405 tau_->setleadCand(*leadingCand);
0406
0407
0408 Iter leadingChargedCand =
0409 leadCand(getCollection(kSignal, kChargedHadron)->begin(), getCollection(kSignal, kChargedHadron)->end());
0410
0411 if (leadingChargedCand != getCollection(kSignal, kChargedHadron)->end())
0412 tau_->setleadChargedHadrCand(*leadingChargedCand);
0413
0414
0415 Iter leadingGammaCand = leadCand(getCollection(kSignal, kGamma)->begin(), getCollection(kSignal, kGamma)->end());
0416
0417 if (leadingGammaCand != getCollection(kSignal, kGamma)->end())
0418 tau_->setleadNeutralCand(*leadingGammaCand);
0419 }
0420 return std::move(tau_);
0421 }
0422 }