61 minGammaEt_(cfg.getParameter<double>(
"minGammaEt")),
70 double trackPerr = track.
p() * (track.
ptError() / track.
pt());
71 return trackPerr * trackPerr;
76 double tauPx_modified = tau.
px() + sf * addP4.px();
77 double tauPy_modified = tau.
py() + sf * addP4.py();
78 double tauPz_modified = tau.
pz() + sf * addP4.pz();
79 double tauMass = tau.
mass();
80 double tauEn_modified =
sqrt(tauPx_modified * tauPx_modified + tauPy_modified * tauPy_modified +
81 tauPz_modified * tauPz_modified + tauMass * tauMass);
83 tau.
setP4(tauP4_modified);
88 tau.
setP4(tauP4_modified);
93 template <
class Base,
class Der>
105 std::cout <<
"<PFRecoTauEnergyAlgorithmPlugin::operator()>:" << std::endl;
106 std::cout <<
"tau: Pt = " << tau.
pt() <<
", eta = " << tau.
eta() <<
", phi = " << tau.
phi()
107 <<
" (En = " << tau.
energy() <<
", decayMode = " << tau.
decayMode() <<
")" << std::endl;
111 std::vector<reco::CandidatePtr> addNeutrals;
113 const auto& jetConstituents = tau.
jetRef()->daughterPtrVector();
114 for (
const auto& jetConstituent : jetConstituents) {
115 int jetConstituentPdgId =
std::abs(jetConstituent->pdgId());
117 (jetConstituentPdgId == 22 && jetConstituent->et() >
minGammaEt_)))
120 bool isSignalCand =
false;
122 for (
const auto& signalCand : signalCands) {
129 double dR =
deltaR(jetConstituent->p4(), tau.
p4());
131 if (jetConstituentPdgId == 130)
133 else if (jetConstituentPdgId == 22)
136 addNeutrals.push_back(jetConstituent);
137 addNeutralsSumP4 += jetConstituent->p4();
141 std::cout <<
"addNeutralsSumP4: En = " << addNeutralsSumP4.energy() << std::endl;
144 unsigned numNonPFCandTracks = 0;
145 double nonPFCandTracksSumP = 0.;
146 double nonPFCandTracksSumPerr2 = 0.;
148 for (std::vector<PFRecoTauChargedHadron>::const_iterator
chargedHadron = chargedHadrons.begin();
152 ++numNonPFCandTracks;
154 if (chargedHadronTrack !=
nullptr) {
155 nonPFCandTracksSumP += chargedHadronTrack->
p();
156 nonPFCandTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);
161 std::cout <<
"nonPFCandTracksSumP = " << nonPFCandTracksSumP <<
" +/- " <<
sqrt(nonPFCandTracksSumPerr2)
162 <<
" (numNonPFCandTracks = " << numNonPFCandTracks <<
")" << std::endl;
165 if (numNonPFCandTracks == 0) {
170 std::cout <<
"easy case: all tracks are associated to PFCandidates --> leaving tau momentum untouched." 173 updateTauP4(tau, 1., addNeutralsSumP4);
182 if (nonPFCandTracksSumP < addNeutralsSumP4.energy()) {
183 double scaleFactor = 1. - nonPFCandTracksSumP / addNeutralsSumP4.energy();
184 if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
186 <<
"Failed to compute tau energy --> killing tau candidate !!" << std::endl;
191 std::cout <<
"case (2): addNeutralsSumEn > nonPFCandTracksSumP --> adjusting tau momentum." << std::endl;
193 updateTauP4(tau, scaleFactor, addNeutralsSumP4);
199 std::vector<reco::CandidatePtr> mergedNeutrals;
201 for (std::vector<PFRecoTauChargedHadron>::const_iterator
chargedHadron = chargedHadrons.begin();
205 const std::vector<reco::CandidatePtr>& neutralPFCands =
chargedHadron->getNeutralPFCandidates();
206 for (std::vector<reco::CandidatePtr>::const_iterator neutralPFCand = neutralPFCands.begin();
207 neutralPFCand != neutralPFCands.end();
209 mergedNeutrals.push_back(*neutralPFCand);
210 mergedNeutralsSumP4 += (*neutralPFCand)->p4();
215 std::cout <<
"mergedNeutralsSumP4: En = " << mergedNeutralsSumP4.energy() << std::endl;
220 if (nonPFCandTracksSumP < (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy())) {
221 double scaleFactor = ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy()) - nonPFCandTracksSumP) /
222 mergedNeutralsSumP4.energy();
223 if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
225 <<
"Failed to compute tau energy --> killing tau candidate !!" << std::endl;
230 size_t numChargedHadrons = chargedHadrons.size();
231 for (
size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
237 diffP4 += (chargedHadron.
p4() - chargedHadron_modified.
p4());
241 std::cout <<
"case (3): (addNeutralsSumEn + mergedNeutralsSumEn) > nonPFCandTracksSumP --> adjusting tau " 245 updateTauP4(tau, -1., diffP4);
250 unsigned numChargedHadronNeutrals = 0;
251 std::vector<reco::CandidatePtr> chargedHadronNeutrals;
253 for (std::vector<PFRecoTauChargedHadron>::const_iterator
chargedHadron = chargedHadrons.begin();
257 ++numChargedHadronNeutrals;
258 chargedHadronNeutrals.push_back(
chargedHadron->getChargedPFCandidate());
259 chargedHadronNeutralsSumP4 +=
chargedHadron->getChargedPFCandidate()->p4();
263 std::cout <<
"chargedHadronNeutralsSumP4: En = " << chargedHadronNeutralsSumP4.energy()
264 <<
" (numChargedHadronNeutrals = " << numChargedHadronNeutrals <<
")" << std::endl;
269 if (nonPFCandTracksSumP <
270 (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy())) {
272 ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy()) -
273 nonPFCandTracksSumP) /
274 chargedHadronNeutralsSumP4.energy();
275 if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
277 <<
"Failed to compute tau energy --> killing tau candidate !!" << std::endl;
282 size_t numChargedHadrons = chargedHadrons.size();
283 for (
size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
289 double chargedHadronPx_modified = scaleFactor * chargedPFCand->px();
290 double chargedHadronPy_modified = scaleFactor * chargedPFCand->py();
291 double chargedHadronPz_modified = scaleFactor * chargedPFCand->pz();
293 chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
294 chargedHadron_modified.
setP4(chargedHadronP4_modified);
296 diffP4 += (chargedHadron.
p4() - chargedHadron_modified.
p4());
300 std::cout <<
"case (4): (addNeutralsSumEn + mergedNeutralsSumEn + chargedHadronNeutralsSumEn) > " 301 "nonPFCandTracksSumP --> adjusting momenta of tau and chargedHadrons." 304 updateTauP4(tau, -1., diffP4);
307 double allTracksSumP = 0.;
308 double allTracksSumPerr2 = 0.;
310 for (std::vector<PFRecoTauChargedHadron>::const_iterator
chargedHadron = chargedHadrons.begin();
316 if (chargedHadronTrack !=
nullptr) {
317 allTracksSumP += chargedHadronTrack->
p();
318 allTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);
320 edm::LogInfo(
"PFRecoTauEnergyAlgorithmPlugin::operator()")
321 <<
"PFRecoTauChargedHadron has no associated reco::Track !!";
329 std::cout <<
"allTracksSumP = " << allTracksSumP <<
" +/- " <<
sqrt(allTracksSumPerr2) << std::endl;
331 double allNeutralsSumEn = 0.;
333 for (
const auto& signalCand : signalCands) {
335 std::cout <<
"Candidate #" << signalCand.id() <<
":" << signalCand.key() <<
":" 336 <<
" Pt = " << (signalCand)->
pt() <<
", eta = " << (signalCand)->
eta()
337 <<
", phi = " << (signalCand)->phi() << std::endl;
343 <<
" ECAL = " << (pfCand)->ecalEnergy() <<
"," 344 <<
" HCAL = " << (pfCand)->hcalEnergy() <<
"," 345 <<
" HO = " << (pfCand)->hoEnergy() << std::endl;
353 allNeutralsSumEn += pfCand->
hoEnergy();
356 allNeutralsSumEn += addNeutralsSumP4.energy();
357 if (allNeutralsSumEn < 0.)
358 allNeutralsSumEn = 0.;
360 std::cout <<
"allNeutralsSumEn = " << allNeutralsSumEn << std::endl;
362 if (allNeutralsSumEn > allTracksSumP) {
364 size_t numChargedHadrons = chargedHadrons.size();
365 for (
size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
372 std::cout <<
"chargedHadron #" << iChargedHadron <<
": changing En = " << chargedHadron.
energy()
373 <<
" to " << chargedHadron_modified.
energy() << std::endl;
381 if (chTrack !=
nullptr) {
382 double chargedHadronPx_modified = chTrack->
px();
383 double chargedHadronPy_modified = chTrack->
py();
384 double chargedHadronPz_modified = chTrack->
pz();
386 chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
389 <<
"PFRecoTauChargedHadron has no associated reco::Track !!" << std::endl;
391 chargedHadron.
print();
394 chargedHadron_modified.
setP4(chargedHadronP4_modified);
396 std::cout <<
"chargedHadron #" << iChargedHadron <<
": changing En = " << chargedHadron.
energy()
397 <<
" to " << chargedHadron_modified.
energy() << std::endl;
404 std::cout <<
"case (5): allNeutralsSumEn > allTracksSumP --> adjusting momenta of tau and chargedHadrons." 407 updateTauP4(tau, scaleFactor - 1., tau.
p4());
412 size_t numChargedHadrons = chargedHadrons.size();
413 for (
size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
421 if (chargedHadronTrack !=
nullptr) {
422 double trackP = chargedHadronTrack->
p();
423 double trackPerr2 = getTrackPerr2(*chargedHadronTrack);
425 std::cout <<
"trackP = " << trackP <<
" +/- " <<
sqrt(trackPerr2) << std::endl;
429 double trackP_modified = (trackP * (allTracksSumPerr2 - trackPerr2) +
430 trackPerr2 * (allNeutralsSumEn - (allTracksSumP - trackP))) /
435 if (trackP_modified < 1.
e-1)
436 trackP_modified = 1.e-1;
438 std::cout <<
"trackP (modified) = " << trackP_modified << std::endl;
441 if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
443 <<
"Failed to compute tau energy --> killing tau candidate !!" << std::endl;
447 double chargedHadronPx_modified = scaleFactor * chargedHadronTrack->
px();
448 double chargedHadronPy_modified = scaleFactor * chargedHadronTrack->
py();
449 double chargedHadronPz_modified = scaleFactor * chargedHadronTrack->
pz();
451 chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
453 edm::LogInfo(
"PFRecoTauEnergyAlgorithmPlugin::operator()")
454 <<
"PFRecoTauChargedHadron has no associated reco::Track !!";
456 chargedHadron.
print();
459 chargedHadron_modified.
setP4(chargedHadronP4_modified);
461 std::cout <<
"chargedHadron #" << iChargedHadron <<
": changing En = " << chargedHadron.
energy()
462 <<
" to " << chargedHadron_modified.
energy() << std::endl;
470 <<
"case (6): allNeutralsSumEn < allTracksSumP --> adjusting momenta of tau and chargedHadrons." 473 updateTauP4(tau, scaleFactor - 1., tau.
p4());
480 std::cout <<
"case (7): allNeutralsSumEn < allTracksSumP not compatible with tau decay mode hypothesis " 481 "--> killing tau candidate." 493 std::cout <<
"undefined case: you should never come here !!" << std::endl;
507 "PFRecoTauEnergyAlgorithmPlugin");
double p() const
momentum vector magnitude
double ecalEnergy() const
return corrected Ecal energy
double eta() const final
momentum pseudorapidity
double px() const final
x coordinate of momentum vector
double pt() const final
transverse momentum
double px() const
x coordinate of momentum vector
void beginEvent() override
constexpr bool isFinite(T x)
double minNeutralHadronEt_
hadronicDecayMode decayMode() const
void print(std::ostream &stream=std::cout) const
double pz() const final
z coordinate of momentum vector
const std::vector< RecoTauPiZero > & signalPiZeroCandidates() const
Retrieve the association of signal region gamma candidates into candidate PiZeros.
double pt() const
track transverse momentum
const JetBaseRef & jetRef() const
double ptError() const
error on Pt (set to 1000 TeV if charge==0 for safety)
double energy() const final
energy
Abs< T >::type abs(const T &t)
~PFRecoTauEnergyAlgorithmPlugin() override
const LorentzVector & p4() const final
four-momentum Lorentz vector
std::vector< CandidatePtr > neutralPFCandidates_
double pz() const
z coordinate of momentum vector
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
bool isPtrEqual(const edm::Ptr< Base > &b, const edm::Ptr< Der > &d)
double hoEnergy() const
return corrected Hcal energy
double dRaddNeutralHadron_
double py() const final
y coordinate of momentum vector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
const reco::Track * getTrackFromChargedHadron(const reco::PFRecoTauChargedHadron &chargedHadron)
Particle reconstructed by the particle flow algorithm.
double hcalEnergy() const
return corrected Hcal energy
PFRecoTauEnergyAlgorithmPlugin(const edm::ParameterSet &, edm::ConsumesCollector &&iC)
bool algoIs(PFRecoTauChargedHadronAlgorithm algo) const
Check whether a given algo produced this charged hadron.
#define DEFINE_EDM_PLUGIN(factory, type, name)
std::vector< PFRecoTauChargedHadron > & signalTauChargedHadronCandidatesRestricted()
const std::vector< reco::CandidatePtr > & signalCands() const
Candidates in signal region.
void setStatus(int status) final
set status word
void setChargedHadronP4(reco::PFRecoTauChargedHadron &chargedHadron, double scaleFactor_neutralPFCands=1.0)
const CandidatePtr & getChargedPFCandidate() const
reference to "charged" PFCandidate (either charged PFCandidate or PFNeutralHadron) ...
double phi() const final
momentum azimuthal angle
void setP4(const LorentzVector &p4) final
set 4-momentum
void operator()(reco::PFTau &) const override
double py() const
y coordinate of momentum vector
double mass() const final
mass
const std::vector< CandidatePtr > & getNeutralPFCandidates() const
references to additional neutral PFCandidates
reco::Candidate::LorentzVector compChargedHadronP4fromPxPyPz(double, double, double)