Back to home page

Project CMSSW displayed by LXR
 
 

    

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     // Initialize tau
0030     tau_ = std::make_unique<PFTau>();
0031 
0032     copyGammas_ = copyGammasFromPiZeros;
0033     // Initialize our Accessors
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     // Build our temporary sorted collections, since you can't use stl sorts on
0045     // RefVectors
0046     for (auto const& colkey : collections_) {
0047       // Build an empty list for each collection
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       // Keep track of the four vector of the signal vector products added so far.
0059       // If a photon add it if we are not using PiZeros to build the gammas
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     // Add to global collection
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     // Reserve global collection as well
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   }  // namespace
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     // CV: need to make sure that PFGammas merged with ChargedHadrons are not part of PiZeros
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       // If we are building the gammas with the pizeros, resize that
0176       // vector as well
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       // Copy the daughter gammas into the gamma collection if desired
0192       if (copyGammas_) {
0193         // If we are using the pizeros to build the gammas, make sure we update
0194         // the four vector correctly.
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   // Trivial converter needed for polymorphism
0215   CandidatePtr RecoTauConstructor::convertToPtr(const CandidatePtr& pfPtr) const { return pfPtr; }
0216 
0217   namespace {
0218     // Make sure the two products come from the same EDM source
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   }  // namespace
0230 
0231   // Convert from a CandidateRef to a Ptr
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   }  // namespace
0250 
0251   void RecoTauConstructor::sortAndCopyIntoTau() {
0252     // The charged hadrons and pizeros are a special case, as we can sort them in situ
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     // Sort each of our sortable collections, and copy them into the final
0267     // tau RefVector.
0268     for (auto const& colkey : collections_) {
0269       SortedListPtr sortedCollection = sortedCollections_[colkey.first];
0270       std::sort(sortedCollection->begin(), sortedCollection->end(), ptDescendingPtr<CandidatePtr>);
0271       // Copy into the real tau collection
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       // If no tracks exist, this is definitely not a tau!
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       // Find the maximum number of PiZeros our parameterization can hold
0315       const unsigned int maxPiZeros = PFTau::kOneProngNPiZero;
0316 
0317       // Determine our track index
0318       unsigned int trackIndex = (nCharged - 1) * (maxPiZeros + 1);
0319 
0320       // Check if we handle the given number of tracks
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   }  // namespace
0329 
0330   std::unique_ptr<reco::PFTau> RecoTauConstructor::get(bool setupLeadingObjects) {
0331     LogDebug("TauConstructorGet") << "Start getting";
0332 
0333     // Copy the sorted collections into the interal tau refvectors
0334     sortAndCopyIntoTau();
0335 
0336     // Setup all the important member variables of the tau
0337     // Set charge of tau
0338     //  tau_->setCharge(
0339     //      sumPFCandCharge(getCollection(kSignal, kChargedHadron)->begin(),
0340     //                      getCollection(kSignal, kChargedHadron)->end()));
0341     // CV: take charge of highest pT charged hadron as charge of tau,
0342     //     in case tau does not have three reconstructed tracks
0343     //    (either because tau is reconstructed as 2prong or because PFRecoTauChargedHadron is built from a PFNeutralHadron)
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     // Set PDG id
0368     tau_->setPdgId(tau_->charge() < 0 ? 15 : -15);
0369 
0370     // Set P4
0371     tau_->setP4(p4_);
0372 
0373     // Set Decay Mode
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     // Set charged isolation quantities
0388     tau_->setisolationPFChargedHadrCandsPtSum(sumPFCandPt(getCollection(kIsolation, kChargedHadron)->begin(),
0389                                                           getCollection(kIsolation, kChargedHadron)->end()));
0390 
0391     // Set gamma isolation quantities
0392     tau_->setisolationPFGammaCandsEtSum(
0393         sumPFCandPt(getCollection(kIsolation, kGamma)->begin(), getCollection(kIsolation, kGamma)->end()));
0394 
0395     // Set em fraction
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       // Find the highest PT object in the signal cone
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       // Hardest charged object in signal cone
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       // Hardest gamma object in signal cone
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 }  // end namespace reco::tau
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!