CMS 3D CMS Logo
PFRecoTauEnergyAlgorithmPlugin.cc
Go to the documentation of this file.
1 /*
2  * =============================================================================
3  * Filename: RecoTauEnergyAlgorithmPlugin.cc
4  *
5  * Description: Determine best estimate for tau energy
6  * for tau candidates reconstructed in different decay modes
7  *
8  * Created: 04/09/2013 11:40:00
9  *
10  * Authors: Christian Veelken (LLR)
11  *
12  * =============================================================================
13  */
14 
15 #include "FWCore/MessageLogger/interface/MessageLogger.h"
16 
17 #include "RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h"
18 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
19 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
20 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
21 #include "DataFormats/JetReco/interface/PFJet.h"
22 #include "DataFormats/TrackReco/interface/Track.h"
23 #include "DataFormats/Common/interface/Ptr.h"
24 #include "DataFormats/Common/interface/RefToPtr.h"
25 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
26 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadronFwd.h"
27 #include "RecoTauTag/RecoTau/interface/pfRecoTauChargedHadronAuxFunctions.h"
28 
29 #include "DataFormats/Math/interface/deltaR.h"
30 #include "FWCore/Utilities/interface/isFinite.h"
31 
32 #include <vector>
33 #include <cmath>
34 
35 namespace reco {
36  namespace tau {
37 
38  class PFRecoTauEnergyAlgorithmPlugin : public RecoTauModifierPlugin {
39  public:
40  explicit PFRecoTauEnergyAlgorithmPlugin(const edm::ParameterSet&, edm::ConsumesCollector&& iC);
41  ~PFRecoTauEnergyAlgorithmPlugin() override;
42  void operator()(reco::PFTau&) const override;
43  void beginEvent() override;
44  void endEvent() override;
45 
46  private:
47  double dRaddNeutralHadron_;
48  double minNeutralHadronEt_;
49  double dRaddPhoton_;
50  double minGammaEt_;
51 
52  int verbosity_;
53  };
54 
55  PFRecoTauEnergyAlgorithmPlugin::PFRecoTauEnergyAlgorithmPlugin(const edm::ParameterSet& cfg,
56  edm::ConsumesCollector&& iC)
57  : RecoTauModifierPlugin(cfg, std::move(iC)),
58  dRaddNeutralHadron_(cfg.getParameter<double>("dRaddNeutralHadron")),
59  minNeutralHadronEt_(cfg.getParameter<double>("minNeutralHadronEt")),
60  dRaddPhoton_(cfg.getParameter<double>("dRaddPhoton")),
61  minGammaEt_(cfg.getParameter<double>("minGammaEt")),
62  verbosity_(cfg.getParameter<int>("verbosity")) {}
63 
64  PFRecoTauEnergyAlgorithmPlugin::~PFRecoTauEnergyAlgorithmPlugin() {}
65 
66  void PFRecoTauEnergyAlgorithmPlugin::beginEvent() {}
67 
68  namespace {
69  double getTrackPerr2(const reco::Track& track) {
70  double trackPerr = track.p() * (track.ptError() / track.pt());
71  return trackPerr * trackPerr;
72  }
73 
74  void updateTauP4(reco::PFTau& tau, double sf, const reco::Candidate::LorentzVector& addP4) {
75  // preserve tau candidate mass when adding extra neutral energy
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);
82  reco::Candidate::LorentzVector tauP4_modified(tauPx_modified, tauPy_modified, tauPz_modified, tauEn_modified);
83  tau.setP4(tauP4_modified);
84  }
85 
86  void killTau(reco::PFTau& tau) {
87  reco::Candidate::LorentzVector tauP4_modified(0., 0., 0., 0.);
88  tau.setP4(tauP4_modified);
89  tau.setStatus(-1);
90  }
91  } // namespace
92 
93  template <class Base, class Der>
94  bool isPtrEqual(const edm::Ptr<Base>& b, const edm::Ptr<Der>& d) {
95  return edm::Ptr<Der>(b) == d;
96  }
97 
98  template <class Base>
99  bool isPtrEqual(const edm::Ptr<Base>& b, const edm::Ptr<Base>& d) {
100  return b == d;
101  }
102 
103  void PFRecoTauEnergyAlgorithmPlugin::operator()(reco::PFTau& tau) const {
104  if (verbosity_) {
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;
108  }
109 
110  // Add high Pt PFNeutralHadrons and PFGammas that are not "used" by tau decay mode object
111  std::vector<reco::CandidatePtr> addNeutrals;
112  reco::Candidate::LorentzVector addNeutralsSumP4;
113  const auto& jetConstituents = tau.jetRef()->daughterPtrVector();
114  for (const auto& jetConstituent : jetConstituents) {
115  int jetConstituentPdgId = std::abs(jetConstituent->pdgId());
116  if (!((jetConstituentPdgId == 130 && jetConstituent->et() > minNeutralHadronEt_) ||
117  (jetConstituentPdgId == 22 && jetConstituent->et() > minGammaEt_)))
118  continue;
119 
120  bool isSignalCand = false;
121  const auto& signalCands = tau.signalCands();
122  for (const auto& signalCand : signalCands) {
123  if (isPtrEqual(jetConstituent, signalCand))
124  isSignalCand = true;
125  }
126  if (isSignalCand)
127  continue;
128 
129  double dR = deltaR(jetConstituent->p4(), tau.p4());
130  double dRadd = -1.;
131  if (jetConstituentPdgId == 130)
132  dRadd = dRaddNeutralHadron_;
133  else if (jetConstituentPdgId == 22)
134  dRadd = dRaddPhoton_;
135  if (dR < dRadd) {
136  addNeutrals.push_back(jetConstituent);
137  addNeutralsSumP4 += jetConstituent->p4();
138  }
139  }
140  if (verbosity_) {
141  std::cout << "addNeutralsSumP4: En = " << addNeutralsSumP4.energy() << std::endl;
142  }
143 
144  unsigned numNonPFCandTracks = 0;
145  double nonPFCandTracksSumP = 0.;
146  double nonPFCandTracksSumPerr2 = 0.;
147  const std::vector<PFRecoTauChargedHadron>& chargedHadrons = tau.signalTauChargedHadronCandidatesRestricted();
148  for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
149  chargedHadron != chargedHadrons.end();
150  ++chargedHadron) {
151  if (chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack)) {
152  ++numNonPFCandTracks;
153  const reco::Track* chargedHadronTrack = getTrackFromChargedHadron(*chargedHadron);
154  if (chargedHadronTrack != nullptr) {
155  nonPFCandTracksSumP += chargedHadronTrack->p();
156  nonPFCandTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);
157  }
158  }
159  }
160  if (verbosity_) {
161  std::cout << "nonPFCandTracksSumP = " << nonPFCandTracksSumP << " +/- " << sqrt(nonPFCandTracksSumPerr2)
162  << " (numNonPFCandTracks = " << numNonPFCandTracks << ")" << std::endl;
163  }
164 
165  if (numNonPFCandTracks == 0) {
166  // This is the easy case:
167  // All tau energy is taken from PFCandidates reconstructed by PFlow algorithm
168  // and there is no issue with double-counting of energy.
169  if (verbosity_) {
170  std::cout << "easy case: all tracks are associated to PFCandidates --> leaving tau momentum untouched."
171  << std::endl;
172  }
173  updateTauP4(tau, 1., addNeutralsSumP4);
174  return;
175  } else {
176  // This is the difficult case:
177  // The tau energy needs to be computed for an arbitrary mix of charged and neutral PFCandidates plus reco::Tracks.
178  // We need to make sure not to double-count energy deposited by reco::Track in ECAL and/or HCAL as neutral PFCandidates.
179 
180  // Check if we have enough energy in collection of PFNeutralHadrons and PFGammas that are not "used" by tau decay mode object
181  // to balance track momenta:
182  if (nonPFCandTracksSumP < addNeutralsSumP4.energy()) {
183  double scaleFactor = 1. - nonPFCandTracksSumP / addNeutralsSumP4.energy();
184  if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
185  edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
186  << "Failed to compute tau energy --> killing tau candidate !!" << std::endl;
187  killTau(tau);
188  return;
189  }
190  if (verbosity_) {
191  std::cout << "case (2): addNeutralsSumEn > nonPFCandTracksSumP --> adjusting tau momentum." << std::endl;
192  }
193  updateTauP4(tau, scaleFactor, addNeutralsSumP4);
194  return;
195  }
196 
197  // Determine which neutral PFCandidates are close to PFChargedHadrons
198  // and have been merged into ChargedHadrons
199  std::vector<reco::CandidatePtr> mergedNeutrals;
200  reco::Candidate::LorentzVector mergedNeutralsSumP4;
201  for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
202  chargedHadron != chargedHadrons.end();
203  ++chargedHadron) {
204  if (chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack)) {
205  const std::vector<reco::CandidatePtr>& neutralPFCands = chargedHadron->getNeutralPFCandidates();
206  for (std::vector<reco::CandidatePtr>::const_iterator neutralPFCand = neutralPFCands.begin();
207  neutralPFCand != neutralPFCands.end();
208  ++neutralPFCand) {
209  mergedNeutrals.push_back(*neutralPFCand);
210  mergedNeutralsSumP4 += (*neutralPFCand)->p4();
211  }
212  }
213  }
214  if (verbosity_) {
215  std::cout << "mergedNeutralsSumP4: En = " << mergedNeutralsSumP4.energy() << std::endl;
216  }
217 
218  // Check if track momenta are balanced by sum of PFNeutralHadrons and PFGammas that are not "used" by tau decay mode object
219  // plus neutral PFCandidates close to PFChargedHadrons:
220  if (nonPFCandTracksSumP < (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy())) {
221  double scaleFactor = ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy()) - nonPFCandTracksSumP) /
222  mergedNeutralsSumP4.energy();
223  if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
224  edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
225  << "Failed to compute tau energy --> killing tau candidate !!" << std::endl;
226  killTau(tau);
227  return;
228  }
229  reco::Candidate::LorentzVector diffP4;
230  size_t numChargedHadrons = chargedHadrons.size();
231  for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
232  const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
233  if (!chargedHadron.getNeutralPFCandidates().empty()) {
234  PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
235  setChargedHadronP4(chargedHadron_modified, scaleFactor);
236  tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
237  diffP4 += (chargedHadron.p4() - chargedHadron_modified.p4());
238  }
239  }
240  if (verbosity_) {
241  std::cout << "case (3): (addNeutralsSumEn + mergedNeutralsSumEn) > nonPFCandTracksSumP --> adjusting tau "
242  "momentum."
243  << std::endl;
244  }
245  updateTauP4(tau, -1., diffP4);
246  return;
247  }
248 
249  // Determine energy sum of all PFNeutralHadrons interpreted as ChargedHadrons with missing track
250  unsigned numChargedHadronNeutrals = 0;
251  std::vector<reco::CandidatePtr> chargedHadronNeutrals;
252  reco::Candidate::LorentzVector chargedHadronNeutralsSumP4;
253  for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
254  chargedHadron != chargedHadrons.end();
255  ++chargedHadron) {
256  if (chargedHadron->algoIs(PFRecoTauChargedHadron::kPFNeutralHadron)) {
257  ++numChargedHadronNeutrals;
258  chargedHadronNeutrals.push_back(chargedHadron->getChargedPFCandidate());
259  chargedHadronNeutralsSumP4 += chargedHadron->getChargedPFCandidate()->p4();
260  }
261  }
262  if (verbosity_) {
263  std::cout << "chargedHadronNeutralsSumP4: En = " << chargedHadronNeutralsSumP4.energy()
264  << " (numChargedHadronNeutrals = " << numChargedHadronNeutrals << ")" << std::endl;
265  }
266 
267  // Check if sum of PFNeutralHadrons and PFGammas that are not "used" by tau decay mode object
268  // plus neutral PFCandidates close to PFChargedHadrons plus PFNeutralHadrons interpreted as ChargedHadrons with missing track balances track momenta
269  if (nonPFCandTracksSumP <
270  (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy())) {
271  double scaleFactor =
272  ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy()) -
273  nonPFCandTracksSumP) /
274  chargedHadronNeutralsSumP4.energy();
275  if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
276  edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
277  << "Failed to compute tau energy --> killing tau candidate !!" << std::endl;
278  killTau(tau);
279  return;
280  }
281  reco::Candidate::LorentzVector diffP4;
282  size_t numChargedHadrons = chargedHadrons.size();
283  for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
284  const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
285  if (chargedHadron.algoIs(PFRecoTauChargedHadron::kPFNeutralHadron)) {
286  PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
287  chargedHadron_modified.neutralPFCandidates_.clear();
288  const CandidatePtr& chargedPFCand = chargedHadron.getChargedPFCandidate();
289  double chargedHadronPx_modified = scaleFactor * chargedPFCand->px();
290  double chargedHadronPy_modified = scaleFactor * chargedPFCand->py();
291  double chargedHadronPz_modified = scaleFactor * chargedPFCand->pz();
292  reco::Candidate::LorentzVector chargedHadronP4_modified = compChargedHadronP4fromPxPyPz(
293  chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
294  chargedHadron_modified.setP4(chargedHadronP4_modified);
295  tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
296  diffP4 += (chargedHadron.p4() - chargedHadron_modified.p4());
297  }
298  }
299  if (verbosity_) {
300  std::cout << "case (4): (addNeutralsSumEn + mergedNeutralsSumEn + chargedHadronNeutralsSumEn) > "
301  "nonPFCandTracksSumP --> adjusting momenta of tau and chargedHadrons."
302  << std::endl;
303  }
304  updateTauP4(tau, -1., diffP4);
305  return;
306  } else {
307  double allTracksSumP = 0.;
308  double allTracksSumPerr2 = 0.;
309  const std::vector<PFRecoTauChargedHadron> chargedHadrons = tau.signalTauChargedHadronCandidatesRestricted();
310  for (std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
311  chargedHadron != chargedHadrons.end();
312  ++chargedHadron) {
313  if (chargedHadron->algoIs(PFRecoTauChargedHadron::kChargedPFCandidate) ||
314  chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack)) {
315  const reco::Track* chargedHadronTrack = getTrackFromChargedHadron(*chargedHadron);
316  if (chargedHadronTrack != nullptr) {
317  allTracksSumP += chargedHadronTrack->p();
318  allTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);
319  } else {
320  // don't print warning as it is caused by missing detailed track inforamtion in lostTrack in miniAOD
321  if (!(chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack) &&
322  chargedHadron->getLostTrackCandidate().isNonnull() &&
323  chargedHadron->getLostTrackCandidate()->bestTrack() == nullptr)) {
324  edm::LogInfo("PFRecoTauEnergyAlgorithmPlugin::operator()")
325  << "PFRecoTauChargedHadron has no associated reco::Track !!";
326  }
327  if (verbosity_) {
328  chargedHadron->print();
329  }
330  }
331  }
332  }
333  if (verbosity_) {
334  std::cout << "allTracksSumP = " << allTracksSumP << " +/- " << sqrt(allTracksSumPerr2) << std::endl;
335  }
336  double allNeutralsSumEn = 0.;
337  const auto& signalCands = tau.signalCands();
338  for (const auto& signalCand : signalCands) {
339  if (verbosity_) {
340  std::cout << "Candidate #" << signalCand.id() << ":" << signalCand.key() << ":"
341  << " Pt = " << (signalCand)->pt() << ", eta = " << (signalCand)->eta()
342  << ", phi = " << (signalCand)->phi() << std::endl;
343  }
344  const PFCandidate* pfCand = dynamic_cast<const PFCandidate*>(&*signalCand);
345  if (pfCand) {
346  if (verbosity_) {
347  std::cout << "calorimeter energy:"
348  << " ECAL = " << (pfCand)->ecalEnergy() << ","
349  << " HCAL = " << (pfCand)->hcalEnergy() << ","
350  << " HO = " << (pfCand)->hoEnergy() << std::endl;
351  }
352  if (edm::isFinite(pfCand->ecalEnergy()))
353  allNeutralsSumEn += pfCand->ecalEnergy();
354  if (edm::isFinite(pfCand->hcalEnergy()))
355  allNeutralsSumEn += pfCand->hcalEnergy();
356  if (edm::isFinite(pfCand->hoEnergy()))
357  allNeutralsSumEn += pfCand->hoEnergy();
358  } else {
359  // TauReco@MiniAOD: individual ECAL and HCAL energies recovered from fractions. Info on HO energy is not (yet?) available in miniAOD.
360  const pat::PackedCandidate* packedCand = dynamic_cast<const pat::PackedCandidate*>(&*signalCand);
361  assert(packedCand); // Taus are built either from reco::PFCandidates or pat::PackedCandidates
362  double caloEn = packedCand->caloFraction() * packedCand->energy();
363  double hcalEn = caloEn * packedCand->hcalFraction();
364  double ecalEn = caloEn - hcalEn;
365  if (verbosity_) {
366  std::cout << "calorimeter energy:"
367  << " ECAL = " << ecalEn << ","
368  << " HCAL = " << hcalEn << ","
369  << " HO unknown for PackedCand" << std::endl;
370  }
371  if (edm::isFinite(ecalEn))
372  allNeutralsSumEn += ecalEn;
373  if (edm::isFinite(hcalEn))
374  allNeutralsSumEn += hcalEn;
375  }
376  }
377  allNeutralsSumEn += addNeutralsSumP4.energy();
378  if (allNeutralsSumEn < 0.)
379  allNeutralsSumEn = 0.;
380  if (verbosity_) {
381  std::cout << "allNeutralsSumEn = " << allNeutralsSumEn << std::endl;
382  }
383  if (allNeutralsSumEn > allTracksSumP) {
384  // Adjust momenta of neutral PFCandidates merged into ChargedHadrons
385  size_t numChargedHadrons = chargedHadrons.size();
386  for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
387  const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
388  if (chargedHadron.algoIs(PFRecoTauChargedHadron::kChargedPFCandidate)) {
389  PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
390  chargedHadron_modified.neutralPFCandidates_.clear();
391  chargedHadron_modified.setP4(chargedHadron.getChargedPFCandidate()->p4());
392  if (verbosity_) {
393  std::cout << "chargedHadron #" << iChargedHadron << ": changing En = " << chargedHadron.energy()
394  << " to " << chargedHadron_modified.energy() << std::endl;
395  }
396  tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
397  } else if (chargedHadron.algoIs(PFRecoTauChargedHadron::kTrack)) {
398  PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
399  chargedHadron_modified.neutralPFCandidates_.clear();
400  reco::Candidate::LorentzVector chargedHadronP4_modified(0., 0., 0., 0.);
401  const reco::Track* chTrack = getTrackFromChargedHadron(chargedHadron);
402  if (chTrack != nullptr) {
403  double chargedHadronPx_modified = chTrack->px();
404  double chargedHadronPy_modified = chTrack->py();
405  double chargedHadronPz_modified = chTrack->pz();
406  chargedHadronP4_modified = compChargedHadronP4fromPxPyPz(
407  chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
408  } else {
409  // don't print warning as it is caused by missing detailed track inforamtion in lostTrack in miniAOD
410  if (!(chargedHadron.algoIs(PFRecoTauChargedHadron::kTrack) &&
411  chargedHadron.getLostTrackCandidate().isNonnull() &&
412  chargedHadron.getLostTrackCandidate()->bestTrack() == nullptr)) {
413  edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
414  << "PFRecoTauChargedHadron has no associated reco::Track !!" << std::endl;
415  }
416  if (verbosity_) {
417  chargedHadron.print();
418  }
419  }
420  chargedHadron_modified.setP4(chargedHadronP4_modified);
421  if (verbosity_) {
422  std::cout << "chargedHadron #" << iChargedHadron << ": changing En = " << chargedHadron.energy()
423  << " to " << chargedHadron_modified.energy() << std::endl;
424  }
425  tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
426  }
427  }
428  double scaleFactor = allNeutralsSumEn / tau.energy();
429  if (verbosity_) {
430  std::cout << "case (5): allNeutralsSumEn > allTracksSumP --> adjusting momenta of tau and chargedHadrons."
431  << std::endl;
432  }
433  updateTauP4(tau, scaleFactor - 1., tau.p4());
434  return;
435  } else {
436  if (numChargedHadronNeutrals == 0 && tau.signalPiZeroCandidates().empty()) {
437  // Adjust momenta of ChargedHadrons build from reco::Tracks to match sum of energy deposits in ECAL + HCAL + HO
438  size_t numChargedHadrons = chargedHadrons.size();
439  for (size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron) {
440  const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
441  if (chargedHadron.algoIs(PFRecoTauChargedHadron::kChargedPFCandidate) ||
442  chargedHadron.algoIs(PFRecoTauChargedHadron::kTrack)) {
443  PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;
444  chargedHadron_modified.neutralPFCandidates_.clear();
445  reco::Candidate::LorentzVector chargedHadronP4_modified(0., 0., 0., 0.);
446  const reco::Track* chargedHadronTrack = getTrackFromChargedHadron(chargedHadron);
447  if (chargedHadronTrack != nullptr) {
448  double trackP = chargedHadronTrack->p();
449  double trackPerr2 = getTrackPerr2(*chargedHadronTrack);
450  if (verbosity_) {
451  std::cout << "trackP = " << trackP << " +/- " << sqrt(trackPerr2) << std::endl;
452  }
453  // CV: adjust track momenta such that difference beeen (measuredTrackP - adjustedTrackP)/sigmaMeasuredTrackP is minimal
454  // (expression derived using Mathematica)
455  double trackP_modified = (trackP * (allTracksSumPerr2 - trackPerr2) +
456  trackPerr2 * (allNeutralsSumEn - (allTracksSumP - trackP))) /
457  allTracksSumPerr2;
458  // CV: trackP_modified may actually become negative in case sum of energy deposits in ECAL + HCAL + HO is small
459  // and one of the tracks has a significantly larger momentum uncertainty than the other tracks.
460  // In this case set track momentum to small positive value.
461  if (trackP_modified < 1.e-1)
462  trackP_modified = 1.e-1;
463  if (verbosity_) {
464  std::cout << "trackP (modified) = " << trackP_modified << std::endl;
465  }
466  double scaleFactor = trackP_modified / trackP;
467  if (!(scaleFactor >= 0. && scaleFactor <= 1.)) {
468  edm::LogWarning("PFRecoTauEnergyAlgorithmPlugin::operator()")
469  << "Failed to compute tau energy --> killing tau candidate !!" << std::endl;
470  killTau(tau);
471  return;
472  }
473  double chargedHadronPx_modified = scaleFactor * chargedHadronTrack->px();
474  double chargedHadronPy_modified = scaleFactor * chargedHadronTrack->py();
475  double chargedHadronPz_modified = scaleFactor * chargedHadronTrack->pz();
476  chargedHadronP4_modified = compChargedHadronP4fromPxPyPz(
477  chargedHadronPx_modified, chargedHadronPy_modified, chargedHadronPz_modified);
478  } else {
479  // don't print warning as it is caused by missing detailed track inforamtion in lostTrack in miniAOD
480  if (!(chargedHadron.algoIs(PFRecoTauChargedHadron::kTrack) &&
481  chargedHadron.getLostTrackCandidate().isNonnull() &&
482  chargedHadron.getLostTrackCandidate()->bestTrack() == nullptr)) {
483  edm::LogInfo("PFRecoTauEnergyAlgorithmPlugin::operator()")
484  << "PFRecoTauChargedHadron has no associated reco::Track !!";
485  }
486  if (verbosity_) {
487  chargedHadron.print();
488  }
489  }
490  chargedHadron_modified.setP4(chargedHadronP4_modified);
491  if (verbosity_) {
492  std::cout << "chargedHadron #" << iChargedHadron << ": changing En = " << chargedHadron.energy()
493  << " to " << chargedHadron_modified.energy() << std::endl;
494  }
495  tau.signalTauChargedHadronCandidatesRestricted()[iChargedHadron] = chargedHadron_modified;
496  }
497  }
498  double scaleFactor = allNeutralsSumEn / tau.energy();
499  if (verbosity_) {
500  std::cout
501  << "case (6): allNeutralsSumEn < allTracksSumP --> adjusting momenta of tau and chargedHadrons."
502  << std::endl;
503  }
504  updateTauP4(tau, scaleFactor - 1., tau.p4());
505  return;
506  } else {
507  // Interpretation of PFNeutralHadrons as ChargedHadrons with missing track and/or reconstruction of extra PiZeros
508  // is not compatible with the fact that sum of reco::Track momenta exceeds sum of energy deposits in ECAL + HCAL + HO:
509  // kill tau candidate (by setting its four-vector to zero)
510  if (verbosity_) {
511  std::cout << "case (7): allNeutralsSumEn < allTracksSumP not compatible with tau decay mode hypothesis "
512  "--> killing tau candidate."
513  << std::endl;
514  }
515  killTau(tau);
516  return;
517  }
518  }
519  }
520  }
521 
522  // CV: You should never come here.
523  if (verbosity_) {
524  std::cout << "undefined case: you should never come here !!" << std::endl;
525  }
526  assert(0);
527  }
528 
529  void PFRecoTauEnergyAlgorithmPlugin::endEvent() {}
530 
531  } // namespace tau
532 } // namespace reco
533 
534 #include "FWCore/Framework/interface/MakerMacros.h"
535 
536 DEFINE_EDM_PLUGIN(RecoTauModifierPluginFactory,
537  reco::tau::PFRecoTauEnergyAlgorithmPlugin,
538  "PFRecoTauEnergyAlgorithmPlugin");
reco::tau::PFRecoTauEnergyAlgorithmPlugin::endEvent
void endEvent() override
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:529
HLT_FULL_cff.track
track
Definition: HLT_FULL_cff.py:11724
MessageLogger.h
reco::tau::PFRecoTauEnergyAlgorithmPlugin::beginEvent
void beginEvent() override
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:66
metsig::tau
Definition: SignAlgoResolutions.h:49
PFCandidate.h
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:631
pat::PackedCandidate::caloFraction
float caloFraction() const
Set the fraction of ECAL+HCAL energy over candidate energy.
Definition: PackedCandidate.h:929
gather_cfg.cout
cout
Definition: gather_cfg.py:144
PFJet.h
reco::tau::isPtrEqual
bool isPtrEqual(const edm::Ptr< Base > &b, const edm::Ptr< Der > &d)
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:94
cms::cuda::assert
assert(be >=bs)
reco::PFTau
Definition: PFTau.h:36
reco::TrackBase::px
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:640
reco::tau::getTrackFromChargedHadron
const reco::Track * getTrackFromChargedHadron(const reco::PFRecoTauChargedHadron &chargedHadron)
Definition: pfRecoTauChargedHadronAuxFunctions.cc:12
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
reco::tau::PFRecoTauEnergyAlgorithmPlugin::dRaddNeutralHadron_
double dRaddNeutralHadron_
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:47
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:122
reco::tau::PFRecoTauEnergyAlgorithmPlugin::~PFRecoTauEnergyAlgorithmPlugin
~PFRecoTauEnergyAlgorithmPlugin() override
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:64
HPSPFTauProducerPuppi_cfi.chargedHadron
chargedHadron
Definition: HPSPFTauProducerPuppi_cfi.py:7
reco::PFRecoTauChargedHadron
Definition: PFRecoTauChargedHadron.h:23
MakerMacros.h
pat::PackedCandidate::energy
double energy() const override
energy
Definition: PackedCandidate.h:469
reco::tau::setChargedHadronP4
void setChargedHadronP4(reco::PFRecoTauChargedHadron &chargedHadron, double scaleFactor_neutralPFCands=1.0)
Definition: pfRecoTauChargedHadronAuxFunctions.cc:31
Track.h
reco::TrackBase::py
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:643
PVValHelper::eta
Definition: PVValidationHelpers.h:70
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
reco::Track
Definition: Track.h:27
reco::tau::PFRecoTauEnergyAlgorithmPlugin::PFRecoTauEnergyAlgorithmPlugin
PFRecoTauEnergyAlgorithmPlugin(const edm::ParameterSet &, edm::ConsumesCollector &&iC)
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:55
DEFINE_EDM_PLUGIN
#define DEFINE_EDM_PLUGIN(factory, type, name)
Definition: PluginFactory.h:124
b
double b
Definition: hdecay.h:118
edm::ParameterSet
Definition: ParameterSet.h:47
reco::tau::compChargedHadronP4fromPxPyPz
reco::Candidate::LorentzVector compChargedHadronP4fromPxPyPz(double, double, double)
Definition: pfRecoTauChargedHadronAuxFunctions.cc:80
deltaR.h
pat::PackedCandidate
Definition: PackedCandidate.h:22
edmplugin::PluginFactory
Definition: PluginFactory.h:34
PackedCandidate.h
reco::tau::PFRecoTauEnergyAlgorithmPlugin::minGammaEt_
double minGammaEt_
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:50
reco::PFRecoTauChargedHadron::kChargedPFCandidate
Definition: PFRecoTauChargedHadron.h:30
Ptr.h
createfilelist.int
int
Definition: createfilelist.py:10
RefToPtr.h
reco::LeafCandidate::p4
const LorentzVector & p4() const final
four-momentum Lorentz vector
Definition: LeafCandidate.h:114
reco::tau::PFRecoTauEnergyAlgorithmPlugin::operator()
void operator()(reco::PFTau &) const override
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:103
reco::tau::PFRecoTauEnergyAlgorithmPlugin
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:38
reco::PFCandidate::ecalEnergy
double ecalEnergy() const
return corrected Ecal energy
Definition: PFCandidate.h:221
edm::Ptr
Definition: AssociationVector.h:31
looper.cfg
cfg
Definition: looper.py:296
RecoTauBuilderPlugins.h
reco::PFRecoTauChargedHadron::neutralPFCandidates_
std::vector< CandidatePtr > neutralPFCandidates_
Definition: PFRecoTauChargedHadron.h:88
reco::LeafCandidate::setP4
void setP4(const LorentzVector &p4) final
set 4-momentum
Definition: LeafCandidate.h:158
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
pat::PackedCandidate::hcalFraction
float hcalFraction() const
Definition: PackedCandidate.h:934
reco::PFRecoTauChargedHadron::kPFNeutralHadron
Definition: PFRecoTauChargedHadron.h:32
isFinite.h
reco::tau::RecoTauModifierPlugin
Definition: RecoTauBuilderPlugins.h:104
reco::PFCandidate::hcalEnergy
double hcalEnergy() const
return corrected Hcal energy
Definition: PFCandidate.h:233
reco::tau::PFRecoTauEnergyAlgorithmPlugin::dRaddPhoton_
double dRaddPhoton_
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:49
multiplicitycorr_cfi.scaleFactor
scaleFactor
Definition: multiplicitycorr_cfi.py:7
reco::deltaR
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
reco::LeafCandidate::energy
double energy() const final
energy
Definition: LeafCandidate.h:125
pfRecoTauChargedHadronAuxFunctions.h
PFRecoTauChargedHadron.h
reco::PFCandidate
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:41
ztail.d
d
Definition: ztail.py:151
reco::TrackBase::pz
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:646
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
edm::isFinite
constexpr bool isFinite(T x)
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
reco::tau::PFRecoTauEnergyAlgorithmPlugin::minNeutralHadronEt_
double minNeutralHadronEt_
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:48
reco::Candidate::LorentzVector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:36
reco::tau::PFRecoTauEnergyAlgorithmPlugin::verbosity_
int verbosity_
Definition: PFRecoTauEnergyAlgorithmPlugin.cc:52
edm::ConsumesCollector
Definition: ConsumesCollector.h:45
PFCandidateFwd.h
PFRecoTauChargedHadronFwd.h
reco::PFCandidate::hoEnergy
double hoEnergy() const
return corrected Hcal energy
Definition: PFCandidate.h:245
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
reco::PFRecoTauChargedHadron::kTrack
Definition: PFRecoTauChargedHadron.h:31
Generated for CMSSW Reference Manual by doxygen 1.8.16