RecoTau/plugins/PFRecoTauEnergyAlgorithmPlugin.cc

0001 /*
0002  * =============================================================================
0003  *       Filename:  RecoTauEnergyAlgorithmPlugin.cc
0004  *
0005  *    Description:  Determine best estimate for tau energy
0006  *                  for tau candidates reconstructed in different decay modes
0007  *
0008  *        Created:  04/09/2013 11:40:00
0009  *
0010  *         Authors:  Christian Veelken (LLR)
0011  *
0012  * =============================================================================
0013  */
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h"
0018 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0019 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0020 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0021 #include "DataFormats/JetReco/interface/PFJet.h"
0022 #include "DataFormats/TrackReco/interface/Track.h"
0023 #include "DataFormats/Common/interface/Ptr.h"
0024 #include "DataFormats/Common/interface/RefToPtr.h"
0025 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
0026 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadronFwd.h"
0027 #include "RecoTauTag/RecoTau/interface/pfRecoTauChargedHadronAuxFunctions.h"
0028
0029 #include "DataFormats/Math/interface/deltaR.h"
0030 #include "FWCore/Utilities/interface/isFinite.h"
0031
0032 #include <vector>
0033 #include <cmath>
0034
0035 namespace reco {
0036   namespace tau {
0037
0038     class PFRecoTauEnergyAlgorithmPlugin : public RecoTauModifierPlugin {
0039     public:
0040       explicit PFRecoTauEnergyAlgorithmPlugin(const edm::ParameterSet&, edm::ConsumesCollector&& iC);
0041       ~PFRecoTauEnergyAlgorithmPlugin() override;
0042       void operator()(reco::PFTau&) const override;
0043       void beginEvent() override;
0044       void endEvent() override;
0045
0046     private:
0047       double dRaddNeutralHadron_;
0048       double minNeutralHadronEt_;
0049       double dRaddPhoton_;
0050       double minGammaEt_;
0051
0052       int verbosity_;
0053     };
0054
0055     PFRecoTauEnergyAlgorithmPlugin::PFRecoTauEnergyAlgorithmPlugin(const edm::ParameterSet& cfg,
0056                                                                    edm::ConsumesCollector&& iC)
0057         : RecoTauModifierPlugin(cfg, std::move(iC)),
0058           dRaddNeutralHadron_(cfg.getParameter<double>("dRaddNeutralHadron")),
0059           minNeutralHadronEt_(cfg.getParameter<double>("minNeutralHadronEt")),
0060           dRaddPhoton_(cfg.getParameter<double>("dRaddPhoton")),
0061           minGammaEt_(cfg.getParameter<double>("minGammaEt")),
0062           verbosity_(cfg.getParameter<int>("verbosity")) {}
0063
0064     PFRecoTauEnergyAlgorithmPlugin::~PFRecoTauEnergyAlgorithmPlugin() {}
0065
0066     void PFRecoTauEnergyAlgorithmPlugin::beginEvent() {}
0067
0068     namespace {
0069       double getTrackPerr2(const reco::Track& track) {
0070         double trackPerr = track.p() * (track.ptError() / track.pt());
0071         return trackPerr * trackPerr;
0072       }
0073
0074       void updateTauP4(reco::PFTau& tau, double sf, const reco::Candidate::LorentzVector& addP4) {
0075         // preserve tau candidate mass when adding extra neutral energy
0076         double tauPx_modified = tau.px() + sf * addP4.px();
0077         double tauPy_modified = tau.py() + sf * addP4.py();
0078         double tauPz_modified = tau.pz() + sf * addP4.pz();
0079         double tauMass = tau.mass();
0080         double tauEn_modified = sqrt(tauPx_modified * tauPx_modified + tauPy_modified * tauPy_modified +
0081                                      tauPz_modified * tauPz_modified + tauMass * tauMass);
0082         reco::Candidate::LorentzVector tauP4_modified(tauPx_modified, tauPy_modified, tauPz_modified, tauEn_modified);
0083         tau.setP4(tauP4_modified);
0084       }
0085
0086       void killTau(reco::PFTau& tau) {
0087         reco::Candidate::LorentzVector tauP4_modified(0., 0., 0., 0.);
0088         tau.setP4(tauP4_modified);
0089         tau.setStatus(-1);
0090       }
0091     }  // namespace
0092
0093     template <class Base, class Der>
0094     bool isPtrEqual(const edm::Ptr<Base>& b, const edm::Ptr<Der>& d) {
0095       return edm::Ptr<Der>(b) == d;
0096     }
0097
0098     template <class Base>
0099     bool isPtrEqual(const edm::Ptr<Base>& b, const edm::Ptr<Base>& d) {
0100       return b == d;
0101     }
0102
0103     void PFRecoTauEnergyAlgorithmPlugin::operator()(reco::PFTau& tau) const {
0104       if (verbosity_) {
0105         std::cout << "<PFRecoTauEnergyAlgorithmPlugin::operator()>:" << std::endl;
0106         std::cout << "tau: Pt = " << tau.pt() << ", eta = " << tau.eta() << ", phi = " << tau.phi()
0107                   << " (En = " << tau.energy() << ", decayMode = " << tau.decayMode() << ")" << std::endl;
0108       }
0109
0110       // Add high Pt PFNeutralHadrons and PFGammas that are not "used" by tau decay mode object
0111       std::vector<reco::CandidatePtr> addNeutrals;
0112       reco::Candidate::LorentzVector addNeutralsSumP4;
0113       const auto& jetConstituents = tau.jetRef()->daughterPtrVector();
0114       for (const auto& jetConstituent : jetConstituents) {
0115         int jetConstituentPdgId = std::abs(jetConstituent->pdgId());
0116         if (!((jetConstituentPdgId == 130 && jetConstituent->et() > minNeutralHadronEt_) ||
0117               (jetConstituentPdgId == 22 && jetConstituent->et() > minGammaEt_)))
0118           continue;
0119
0120         bool isSignalCand = false;
0121         const auto& signalCands = tau.signalCands();
0122         for (const auto& signalCand : signalCands) {
0123           if (isPtrEqual(jetConstituent, signalCand))
0124             isSignalCand = true;
0125         }
0126         if (isSignalCand)
0127           continue;
0128
0129         double dR = deltaR(jetConstituent->p4(), tau.p4());
0130         double dRadd = -1.;
0131         if (jetConstituentPdgId == 130)
0132           dRadd = dRaddNeutralHadron_;
0133         else if (jetConstituentPdgId == 22)
0134           dRadd = dRaddPhoton_;
0135         if (dR < dRadd) {
0136           addNeutrals.push_back(jetConstituent);
0137           addNeutralsSumP4 += jetConstituent->p4();
0138         }
0139       }
0140       if (verbosity_) {
0141         std::cout << "addNeutralsSumP4: En = " << addNeutralsSumP4.energy() << std::endl;
0142       }
0143
0144       unsigned numNonPFCandTracks = 0;
0145       double nonPFCandTracksSumP = 0.;
0146       double nonPFCandTracksSumPerr2 = 0.;
0147       const std::vector<PFRecoTauChargedHadron>& chargedHadrons = tau.signalTauChargedHadronCandidatesRestricted();
0148       for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
0149            chargedHadron != chargedHadrons.end();
0150            ++chargedHadron) {
0151         if (chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack)) {
0152           ++numNonPFCandTracks;
0153           const reco::Track* chargedHadronTrack = getTrackFromChargedHadron(*chargedHadron);
0154           if (chargedHadronTrack != nullptr) {
0155             nonPFCandTracksSumP += chargedHadronTrack->p();
0156             nonPFCandTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);
0157           }
0158         }
0159       }
0160       if (verbosity_) {
0161         std::cout << "nonPFCandTracksSumP = " << nonPFCandTracksSumP << " +/- " << sqrt(nonPFCandTracksSumPerr2)
0162                   << " (numNonPFCandTracks = " << numNonPFCandTracks << ")" << std::endl;
0163       }
0164
0165       if (numNonPFCandTracks == 0) {
0166         // This is the easy case:
0167         // All tau energy is taken from PFCandidates reconstructed by PFlow algorithm
0168         // and there is no issue with double-counting of energy.
0169         if (verbosity_) {
0170           std::cout << "easy case: all tracks are associated to PFCandidates --> leaving tau momentum untouched."
0171                     << std::endl;
0172         }
0173         updateTauP4(tau, 1., addNeutralsSumP4);
0174         return;
0175       } else {
0176         // This is the difficult case:
0177         // The tau energy needs to be computed for an arbitrary mix of charged and neutral PFCandidates plus reco::Tracks.
0178         // We need to make sure not to double-count energy deposited by reco::Track in ECAL and/or HCAL as neutral PFCandidates.
0179
0180         // Check if we have enough energy in collection of PFNeutralHadrons and PFGammas that are not "used" by tau decay mode object
0181         // to balance track momenta:
0182         if (nonPFCandTracksSumP < addNeutralsSumP4.energy()) {
0183           double scaleFactor = 1. - nonPFCandTracksSumP / addNeutralsSumP4.energy();
0184           if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
0185             edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
0186                 << "Failed to compute tau energy --> killing tau candidate !!" << std::endl;
0187             killTau(tau);
0188             return;
0189           }
0190           if (verbosity_) {
0191             std::cout << "case (2): addNeutralsSumEn > nonPFCandTracksSumP --> adjusting tau momentum." << std::endl;
0192           }
0193           updateTauP4(tau, scaleFactor, addNeutralsSumP4);
0194           return;
0195         }
0196
0197         // Determine which neutral PFCandidates are close to PFChargedHadrons
0198         // and have been merged into ChargedHadrons
0199         std::vector<reco::CandidatePtr> mergedNeutrals;
0200         reco::Candidate::LorentzVector mergedNeutralsSumP4;
0201         for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
0202              chargedHadron != chargedHadrons.end();
0203              ++chargedHadron) {
0204           if (chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack)) {
0205             const std::vector<reco::CandidatePtr>& neutralPFCands = chargedHadron->getNeutralPFCandidates();
0206             for (std::vector<reco::CandidatePtr>::const_iterator neutralPFCand = neutralPFCands.begin();
0207                  neutralPFCand != neutralPFCands.end();
0208                  ++neutralPFCand) {
0209               mergedNeutrals.push_back(*neutralPFCand);
0210               mergedNeutralsSumP4 += (*neutralPFCand)->p4();
0211             }
0212           }
0213         }
0214         if (verbosity_) {
0215           std::cout << "mergedNeutralsSumP4: En = " << mergedNeutralsSumP4.energy() << std::endl;
0216         }
0217
0218         // Check if track momenta are balanced by sum of PFNeutralHadrons and PFGammas that are not "used" by tau decay mode object
0219         // plus neutral PFCandidates close to PFChargedHadrons:
0220         if (nonPFCandTracksSumP < (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy())) {
0221           double scaleFactor = ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy()) - nonPFCandTracksSumP) /
0222                                mergedNeutralsSumP4.energy();
0223           if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
0224             edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
0225                 << "Failed to compute tau energy --> killing tau candidate !!" << std::endl;
0226             killTau(tau);
0227             return;
0228           }
0229           reco::Candidate::LorentzVector diffP4;
0230           size_t numChargedHadrons = chargedHadrons.size();
0231           for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
0232             const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
0233             if (!chargedHadron.getNeutralPFCandidates().empty()) {
0234               PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
0235               setChargedHadronP4(chargedHadron_modified, scaleFactor);
0236               tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
0237               diffP4 += (chargedHadron.p4() - chargedHadron_modified.p4());
0238             }
0239           }
0240           if (verbosity_) {
0241             std::cout << "case (3): (addNeutralsSumEn + mergedNeutralsSumEn) > nonPFCandTracksSumP --> adjusting tau "
0242                          "momentum."
0243                       << std::endl;
0244           }
0245           updateTauP4(tau, -1., diffP4);
0246           return;
0247         }
0248
0249         // Determine energy sum of all PFNeutralHadrons interpreted as ChargedHadrons with missing track
0250         unsigned numChargedHadronNeutrals = 0;
0251         std::vector<reco::CandidatePtr> chargedHadronNeutrals;
0252         reco::Candidate::LorentzVector chargedHadronNeutralsSumP4;
0253         for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
0254              chargedHadron != chargedHadrons.end();
0255              ++chargedHadron) {
0256           if (chargedHadron->algoIs(PFRecoTauChargedHadron::kPFNeutralHadron)) {
0257             ++numChargedHadronNeutrals;
0258             chargedHadronNeutrals.push_back(chargedHadron->getChargedPFCandidate());
0259             chargedHadronNeutralsSumP4 += chargedHadron->getChargedPFCandidate()->p4();
0260           }
0261         }
0262         if (verbosity_) {
0263           std::cout << "chargedHadronNeutralsSumP4: En = " << chargedHadronNeutralsSumP4.energy()
0264                     << " (numChargedHadronNeutrals = " << numChargedHadronNeutrals << ")" << std::endl;
0265         }
0266
0267         // Check if sum of PFNeutralHadrons and PFGammas that are not "used" by tau decay mode object
0268         // plus neutral PFCandidates close to PFChargedHadrons plus PFNeutralHadrons interpreted as ChargedHadrons with missing track balances track momenta
0269         if (nonPFCandTracksSumP <
0270             (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy())) {
0271           double scaleFactor =
0272               ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy()) -
0273                nonPFCandTracksSumP) /
0274               chargedHadronNeutralsSumP4.energy();
0275           if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
0276             edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
0277                 << "Failed to compute tau energy --> killing tau candidate !!" << std::endl;
0278             killTau(tau);
0279             return;
0280           }
0281           reco::Candidate::LorentzVector diffP4;
0282           size_t numChargedHadrons = chargedHadrons.size();
0283           for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
0284             const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
0285             if (chargedHadron.algoIs(PFRecoTauChargedHadron::kPFNeutralHadron)) {
0286               PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
0287               chargedHadron_modified.neutralPFCandidates_.clear();
0288               const CandidatePtr& chargedPFCand = chargedHadron.getChargedPFCandidate();
0289               double chargedHadronPx_modified = scaleFactor * chargedPFCand->px();
0290               double chargedHadronPy_modified = scaleFactor * chargedPFCand->py();
0291               double chargedHadronPz_modified = scaleFactor * chargedPFCand->pz();
0292               reco::Candidate::LorentzVector chargedHadronP4_modified = compChargedHadronP4fromPxPyPz(
0293                   chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
0294               chargedHadron_modified.setP4(chargedHadronP4_modified);
0295               tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
0296               diffP4 += (chargedHadron.p4() - chargedHadron_modified.p4());
0297             }
0298           }
0299           if (verbosity_) {
0300             std::cout << "case (4): (addNeutralsSumEn + mergedNeutralsSumEn + chargedHadronNeutralsSumEn) > "
0301                          "nonPFCandTracksSumP --> adjusting momenta of tau and chargedHadrons."
0302                       << std::endl;
0303           }
0304           updateTauP4(tau, -1., diffP4);
0305           return;
0306         } else {
0307           double allTracksSumP = 0.;
0308           double allTracksSumPerr2 = 0.;
0309           const std::vector<PFRecoTauChargedHadron> chargedHadrons = tau.signalTauChargedHadronCandidatesRestricted();
0310           for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
0311                chargedHadron != chargedHadrons.end();
0312                ++chargedHadron) {
0313             if (chargedHadron->algoIs(PFRecoTauChargedHadron::kChargedPFCandidate) ||
0314                 chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack)) {
0315               const reco::Track* chargedHadronTrack = getTrackFromChargedHadron(*chargedHadron);
0316               if (chargedHadronTrack != nullptr) {
0317                 allTracksSumP += chargedHadronTrack->p();
0318                 allTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);
0319               } else {
0320                 // don't print warning as it is caused by missing detailed track inforamtion in lostTrack in miniAOD
0321                 if (!(chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack) &&
0322                       chargedHadron->getLostTrackCandidate().isNonnull() &&
0323                       chargedHadron->getLostTrackCandidate()->bestTrack() == nullptr)) {
0324                   edm::LogInfo("PFRecoTauEnergyAlgorithmPlugin::operator()")
0325                       << "PFRecoTauChargedHadron has no associated reco::Track !!";
0326                 }
0327                 if (verbosity_) {
0328                   chargedHadron->print();
0329                 }
0330               }
0331             }
0332           }
0333           if (verbosity_) {
0334             std::cout << "allTracksSumP = " << allTracksSumP << " +/- " << sqrt(allTracksSumPerr2) << std::endl;
0335           }
0336           double allNeutralsSumEn = 0.;
0337           const auto& signalCands = tau.signalCands();
0338           for (const auto& signalCand : signalCands) {
0339             if (verbosity_) {
0340               std::cout << "Candidate #" << signalCand.id() << ":" << signalCand.key() << ":"
0341                         << " Pt = " << (signalCand)->pt() << ", eta = " << (signalCand)->eta()
0342                         << ", phi = " << (signalCand)->phi() << std::endl;
0343             }
0344             const PFCandidate* pfCand = dynamic_cast<const PFCandidate*>(&*signalCand);
0345             if (pfCand) {
0346               if (verbosity_) {
0347                 std::cout << "calorimeter energy:"
0348                           << " ECAL = " << (pfCand)->ecalEnergy() << ","
0349                           << " HCAL = " << (pfCand)->hcalEnergy() << ","
0350                           << " HO = " << (pfCand)->hoEnergy() << std::endl;
0351               }
0352               if (edm::isFinite(pfCand->ecalEnergy()))
0353                 allNeutralsSumEn += pfCand->ecalEnergy();
0354               if (edm::isFinite(pfCand->hcalEnergy()))
0355                 allNeutralsSumEn += pfCand->hcalEnergy();
0356               if (edm::isFinite(pfCand->hoEnergy()))
0357                 allNeutralsSumEn += pfCand->hoEnergy();
0358             } else {
0359               // TauReco@MiniAOD: individual ECAL and HCAL energies recovered from fractions. Info on HO energy is not (yet?) available in miniAOD.
0360               const pat::PackedCandidate* packedCand = dynamic_cast<const pat::PackedCandidate*>(&*signalCand);
0361               assert(packedCand);  // Taus are built either from reco::PFCandidates or pat::PackedCandidates
0362               double caloEn = packedCand->caloFraction() * packedCand->energy();
0363               double hcalEn = caloEn * packedCand->hcalFraction();
0364               double ecalEn = caloEn - hcalEn;
0365               if (verbosity_) {
0366                 std::cout << "calorimeter energy:"
0367                           << " ECAL = " << ecalEn << ","
0368                           << " HCAL = " << hcalEn << ","
0369                           << " HO unknown for PackedCand" << std::endl;
0370               }
0371               if (edm::isFinite(ecalEn))
0372                 allNeutralsSumEn += ecalEn;
0373               if (edm::isFinite(hcalEn))
0374                 allNeutralsSumEn += hcalEn;
0375             }
0376           }
0377           allNeutralsSumEn += addNeutralsSumP4.energy();
0378           if (allNeutralsSumEn < 0.)
0379             allNeutralsSumEn = 0.;
0380           if (verbosity_) {
0381             std::cout << "allNeutralsSumEn = " << allNeutralsSumEn << std::endl;
0382           }
0383           if (allNeutralsSumEn > allTracksSumP) {
0384             // Adjust momenta of neutral PFCandidates merged into ChargedHadrons
0385             size_t numChargedHadrons = chargedHadrons.size();
0386             for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
0387               const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
0388               if (chargedHadron.algoIs(PFRecoTauChargedHadron::kChargedPFCandidate)) {
0389                 PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
0390                 chargedHadron_modified.neutralPFCandidates_.clear();
0391                 chargedHadron_modified.setP4(chargedHadron.getChargedPFCandidate()->p4());
0392                 if (verbosity_) {
0393                   std::cout << "chargedHadron #" << iChargedHadron << ": changing En = " << chargedHadron.energy()
0394                             << " to " << chargedHadron_modified.energy() << std::endl;
0395                 }
0396                 tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
0397               } else if (chargedHadron.algoIs(PFRecoTauChargedHadron::kTrack)) {
0398                 PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
0399                 chargedHadron_modified.neutralPFCandidates_.clear();
0400                 reco::Candidate::LorentzVector chargedHadronP4_modified(0., 0., 0., 0.);
0401                 const reco::Track* chTrack = getTrackFromChargedHadron(chargedHadron);
0402                 if (chTrack != nullptr) {
0403                   double chargedHadronPx_modified = chTrack->px();
0404                   double chargedHadronPy_modified = chTrack->py();
0405                   double chargedHadronPz_modified = chTrack->pz();
0406                   chargedHadronP4_modified = compChargedHadronP4fromPxPyPz(
0407                       chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
0408                 } else {
0409                   // don't print warning as it is caused by missing detailed track inforamtion in lostTrack in miniAOD
0410                   if (!(chargedHadron.algoIs(PFRecoTauChargedHadron::kTrack) &&
0411                         chargedHadron.getLostTrackCandidate().isNonnull() &&
0412                         chargedHadron.getLostTrackCandidate()->bestTrack() == nullptr)) {
0413                     edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
0414                         << "PFRecoTauChargedHadron has no associated reco::Track !!" << std::endl;
0415                   }
0416                   if (verbosity_) {
0417                     chargedHadron.print();
0418                   }
0419                 }
0420                 chargedHadron_modified.setP4(chargedHadronP4_modified);
0421                 if (verbosity_) {
0422                   std::cout << "chargedHadron #" << iChargedHadron << ": changing En = " << chargedHadron.energy()
0423                             << " to " << chargedHadron_modified.energy() << std::endl;
0424                 }
0425                 tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
0426               }
0427             }
0428             double scaleFactor = allNeutralsSumEn / tau.energy();
0429             if (verbosity_) {
0430               std::cout << "case (5): allNeutralsSumEn > allTracksSumP --> adjusting momenta of tau and chargedHadrons."
0431                         << std::endl;
0432             }
0433             updateTauP4(tau, scaleFactor - 1., tau.p4());
0434             return;
0435           } else {
0436             if (numChargedHadronNeutrals == 0 && tau.signalPiZeroCandidates().empty()) {
0437               // Adjust momenta of ChargedHadrons build from reco::Tracks to match sum of energy deposits in ECAL + HCAL + HO
0438               size_t numChargedHadrons = chargedHadrons.size();
0439               for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
0440                 const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
0441                 if (chargedHadron.algoIs(PFRecoTauChargedHadron::kChargedPFCandidate) ||
0442                     chargedHadron.algoIs(PFRecoTauChargedHadron::kTrack)) {
0443                   PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
0444                   chargedHadron_modified.neutralPFCandidates_.clear();
0445                   reco::Candidate::LorentzVector chargedHadronP4_modified(0., 0., 0., 0.);
0446                   const reco::Track* chargedHadronTrack = getTrackFromChargedHadron(chargedHadron);
0447                   if (chargedHadronTrack != nullptr) {
0448                     double trackP = chargedHadronTrack->p();
0449                     double trackPerr2 = getTrackPerr2(*chargedHadronTrack);
0450                     if (verbosity_) {
0451                       std::cout << "trackP = " << trackP << " +/- " << sqrt(trackPerr2) << std::endl;
0452                     }
0453                     // CV: adjust track momenta such that difference beeen (measuredTrackP - adjustedTrackP)/sigmaMeasuredTrackP is minimal
0454                     //    (expression derived using Mathematica)
0455                     double trackP_modified = (trackP * (allTracksSumPerr2 - trackPerr2) +
0456                                               trackPerr2 * (allNeutralsSumEn - (allTracksSumP - trackP))) /
0457                                              allTracksSumPerr2;
0458                     // CV: trackP_modified may actually become negative in case sum of energy deposits in ECAL + HCAL + HO is small
0459                     //     and one of the tracks has a significantly larger momentum uncertainty than the other tracks.
0460                     //     In this case set track momentum to small positive value.
0461                     if (trackP_modified < 1.e-1)
0462                       trackP_modified = 1.e-1;
0463                     if (verbosity_) {
0464                       std::cout << "trackP (modified) = " << trackP_modified << std::endl;
0465                     }
0466                     double scaleFactor = trackP_modified / trackP;
0467                     if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
0468                       edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
0469                           << "Failed to compute tau energy --> killing tau candidate !!" << std::endl;
0470                       killTau(tau);
0471                       return;
0472                     }
0473                     double chargedHadronPx_modified = scaleFactor * chargedHadronTrack->px();
0474                     double chargedHadronPy_modified = scaleFactor * chargedHadronTrack->py();
0475                     double chargedHadronPz_modified = scaleFactor * chargedHadronTrack->pz();
0476                     chargedHadronP4_modified = compChargedHadronP4fromPxPyPz(
0477                         chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
0478                   } else {
0479                     // don't print warning as it is caused by missing detailed track inforamtion in lostTrack in miniAOD
0480                     if (!(chargedHadron.algoIs(PFRecoTauChargedHadron::kTrack) &&
0481                           chargedHadron.getLostTrackCandidate().isNonnull() &&
0482                           chargedHadron.getLostTrackCandidate()->bestTrack() == nullptr)) {
0483                       edm::LogInfo("PFRecoTauEnergyAlgorithmPlugin::operator()")
0484                           << "PFRecoTauChargedHadron has no associated reco::Track !!";
0485                     }
0486                     if (verbosity_) {
0487                       chargedHadron.print();
0488                     }
0489                   }
0490                   chargedHadron_modified.setP4(chargedHadronP4_modified);
0491                   if (verbosity_) {
0492                     std::cout << "chargedHadron #" << iChargedHadron << ": changing En = " << chargedHadron.energy()
0493                               << " to " << chargedHadron_modified.energy() << std::endl;
0494                   }
0495                   tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
0496                 }
0497               }
0498               double scaleFactor = allNeutralsSumEn / tau.energy();
0499               if (verbosity_) {
0500                 std::cout
0501                     << "case (6): allNeutralsSumEn < allTracksSumP --> adjusting momenta of tau and chargedHadrons."
0502                     << std::endl;
0503               }
0504               updateTauP4(tau, scaleFactor - 1., tau.p4());
0505               return;
0506             } else {
0507               // Interpretation of PFNeutralHadrons as ChargedHadrons with missing track and/or reconstruction of extra PiZeros
0508               // is not compatible with the fact that sum of reco::Track momenta exceeds sum of energy deposits in ECAL + HCAL + HO:
0509               // kill tau candidate (by setting its four-vector to zero)
0510               if (verbosity_) {
0511                 std::cout << "case (7): allNeutralsSumEn < allTracksSumP not compatible with tau decay mode hypothesis "
0512                              "--> killing tau candidate."
0513                           << std::endl;
0514               }
0515               killTau(tau);
0516               return;
0517             }
0518           }
0519         }
0520       }
0521
0522       // CV: You should never come here.
0523       if (verbosity_) {
0524         std::cout << "undefined case: you should never come here !!" << std::endl;
0525       }
0526       assert(0);
0527     }
0528
0529     void PFRecoTauEnergyAlgorithmPlugin::endEvent() {}
0530
0531   }  // namespace tau
0532 }  // namespace reco
0533
0534 #include "FWCore/Framework/interface/MakerMacros.h"
0535
0536 DEFINE_EDM_PLUGIN(RecoTauModifierPluginFactory,
0537                   reco::tau::PFRecoTauEnergyAlgorithmPlugin,
0538                   "PFRecoTauEnergyAlgorithmPlugin");