GitHub PullRequestReviewCommentEvent at 2014-01-18 14:42:56
Repository: cmssw
Actor: slava77
Payload: {"comment":{"url":"https://api.github.com/repos/cms-sw/cmssw/pulls/comments/8987917","id":8987917,"diff_hunk":"@@ -0,0 +1,463 @@\n+/*\n+ * =============================================================================\n+ *       Filename:  RecoTauEnergyAlgorithmPlugin.cc\n+ *\n+ *    Description:  Determine best estimate for tau energy\n+ *                  for tau candidates reconstructed in different decay modes\n+ *\n+ *        Created:  04/09/2013 11:40:00\n+ *\n+ *         Authors:  Christian Veelken (LLR)\n+ *\n+ * =============================================================================\n+ */\n+\n+#include \"FWCore/MessageLogger/interface/MessageLogger.h\"\n+\n+#include \"RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h\"\n+#include \"DataFormats/ParticleFlowCandidate/interface/PFCandidate.h\"\n+#include \"DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h\"\n+#include \"DataFormats/JetReco/interface/PFJet.h\"\n+#include \"DataFormats/TrackReco/interface/Track.h\"\n+#include \"DataFormats/Common/interface/Ptr.h\"\n+#include \"DataFormats/Common/interface/RefToPtr.h\"\n+#include \"DataFormats/TauReco/interface/PFRecoTauChargedHadron.h\"\n+#include \"DataFormats/TauReco/interface/PFRecoTauChargedHadronFwd.h\"\n+#include \"RecoTauTag/RecoTau/interface/pfRecoTauChargedHadronAuxFunctions.h\"\n+\n+#include \"DataFormats/Math/interface/deltaR.h\"\n+\n+#include <vector>\n+#include <cmath>\n+\n+namespace reco { namespace tau {\n+\n+class PFRecoTauEnergyAlgorithmPlugin : public RecoTauModifierPlugin\n+{\n+ public:\n+\n+  explicit PFRecoTauEnergyAlgorithmPlugin(const edm::ParameterSet&, edm::ConsumesCollector &&iC);\n+  virtual ~PFRecoTauEnergyAlgorithmPlugin();\n+  void operator()(PFTau&) const;\n+  virtual void beginJob(edm::EDProducer*);\n+  virtual void beginEvent();\n+  virtual void endEvent();\n+\n+ private:\n+  \n+  double dRaddNeutralHadron_;\n+  double minNeutralHadronEt_;\n+  double dRaddPhoton_;\n+  double minGammaEt_;\n+\n+  int verbosity_;\n+};\n+\n+  PFRecoTauEnergyAlgorithmPlugin::PFRecoTauEnergyAlgorithmPlugin(const edm::ParameterSet& cfg, edm::ConsumesCollector &&iC)\n+    : RecoTauModifierPlugin(cfg, std::move(iC)),\n+    dRaddNeutralHadron_(cfg.getParameter<double>(\"dRaddNeutralHadron\")),\n+    minNeutralHadronEt_(cfg.getParameter<double>(\"minNeutralHadronEt\")),\n+    dRaddPhoton_(cfg.getParameter<double>(\"dRaddPhoton\")),\n+    minGammaEt_(cfg.getParameter<double>(\"minGammaEt\"))\n+{\n+  verbosity_ = ( cfg.exists(\"verbosity\") ) ?\n+    cfg.getParameter<int>(\"verbosity\") : 0;\n+}\n+\n+PFRecoTauEnergyAlgorithmPlugin::~PFRecoTauEnergyAlgorithmPlugin()\n+{}\n+\n+void PFRecoTauEnergyAlgorithmPlugin::beginJob(edm::EDProducer* producer)\n+{}\n+\n+void PFRecoTauEnergyAlgorithmPlugin::beginEvent()\n+{}\n+\n+namespace\n+{\n+  double getTrackPerr2(const reco::Track& track)\n+  {\n+    double trackPerr = track.p()*(track.ptError()/track.pt());\n+    return trackPerr*trackPerr;\n+  }\n+\n+  void updateTauP4(PFTau& tau, double sf, const reco::Candidate::LorentzVector& addP4)\n+  {\n+    // preserve tau candidate mass when adding extra neutral energy\n+    double tauPx_modified = tau.px() + sf*addP4.px();\n+    double tauPy_modified = tau.py() + sf*addP4.py();\n+    double tauPz_modified = tau.pz() + sf*addP4.pz();\n+    double tauMass = tau.mass();\n+    double tauEn_modified = sqrt(tauPx_modified*tauPx_modified + tauPy_modified*tauPy_modified + tauPz_modified*tauPz_modified + tauMass*tauMass);\n+    reco::Candidate::LorentzVector tauP4_modified(tauPx_modified, tauPy_modified, tauPz_modified, tauEn_modified);\n+    tau.setP4(tauP4_modified);\n+  }\n+\n+  void killTau(PFTau& tau)\n+  {\n+    reco::Candidate::LorentzVector tauP4_modified(0.,0.,0.,0.);\n+    tau.setP4(tauP4_modified);\n+    tau.setStatus(-1);\n+  }\n+}\n+\n+void PFRecoTauEnergyAlgorithmPlugin::operator()(PFTau& tau) const\n+{\n+  if ( verbosity_ ) {\n+    std::cout << \"<PFRecoTauEnergyAlgorithmPlugin::operator()>:\" << std::endl;\n+    std::cout << \"tau: Pt = \" << tau.pt() << \", eta = \" << tau.eta() << \", phi = \" << tau.phi() << \" (En = \" << tau.energy() << \", decayMode = \" << tau.decayMode() << \")\" << std::endl;\n+  }\n+\n+  // Add high Pt PFNeutralHadrons and PFGammas that are not \"used\" by tau decay mode object\n+  std::vector<reco::PFCandidatePtr> addNeutrals;\n+  reco::Candidate::LorentzVector addNeutralsSumP4;\n+  std::vector<reco::PFCandidatePtr> jetConstituents = tau.jetRef()->getPFConstituents();\n+  for ( std::vector<reco::PFCandidatePtr>::const_iterator jetConstituent = jetConstituents.begin();\n+\tjetConstituent != jetConstituents.end(); ++jetConstituent ) {\n+    reco::PFCandidate::ParticleType jetConstituentType = (*jetConstituent)->particleId();\n+    if ( !((jetConstituentType == reco::PFCandidate::h0    && (*jetConstituent)->et() > minNeutralHadronEt_) ||\n+\t   (jetConstituentType == reco::PFCandidate::gamma && (*jetConstituent)->et() > minGammaEt_        )) ) continue;\n+\n+    bool isSignalPFCand = false;\n+    const std::vector<reco::PFCandidatePtr>& signalPFCands = tau.signalPFCands();\n+    for ( std::vector<reco::PFCandidatePtr>::const_iterator signalPFCand = signalPFCands.begin();\n+\t  signalPFCand != signalPFCands.end(); ++signalPFCand ) {\n+      if ( (*jetConstituent) == (*signalPFCand) ) isSignalPFCand = true;\n+    }\n+    if ( isSignalPFCand ) continue;\n+    \n+    double dR = deltaR((*jetConstituent)->p4(), tau.p4());\n+    double dRadd = -1.;      \n+    if      ( jetConstituentType == reco::PFCandidate::h0    ) dRadd = dRaddNeutralHadron_;\n+    else if ( jetConstituentType == reco::PFCandidate::gamma ) dRadd = dRaddPhoton_;\n+    if ( dR < dRadd ) {\n+      addNeutrals.push_back(*jetConstituent);\n+      addNeutralsSumP4 += (*jetConstituent)->p4();\n+    }\n+  }\n+  if ( verbosity_ ) {\n+    std::cout << \"addNeutralsSumP4: En = \" << addNeutralsSumP4.energy() << std::endl;\n+  }\n+  \n+  unsigned numNonPFCandTracks = 0;\n+  double nonPFCandTracksSumP = 0.;\n+  double nonPFCandTracksSumPerr2 = 0.;\n+  const std::vector<PFRecoTauChargedHadron>& chargedHadrons = tau.signalTauChargedHadronCandidates();\n+  for ( std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();\n+\tchargedHadron != chargedHadrons.end(); ++chargedHadron ) {\n+    if ( chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack) ) {\n+      ++numNonPFCandTracks;\n+      const edm::Ptr<Track>& chargedHadronTrack = chargedHadron->getTrack();\n+      nonPFCandTracksSumP += chargedHadronTrack->p();\n+      nonPFCandTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);\n+    }\n+  }\n+  if ( verbosity_ ) {\n+    std::cout << \"nonPFCandTracksSumP = \" << nonPFCandTracksSumP << \" +/- \" << sqrt(nonPFCandTracksSumPerr2) \n+\t      << \" (numNonPFCandTracks = \" << numNonPFCandTracks << \")\" << std::endl;\n+  }\n+ \n+  if ( numNonPFCandTracks == 0 ) {\n+    // This is the easy case: \n+    // All tau energy is taken from PFCandidates reconstructed by PFlow algorithm\n+    // and there is no issue with double-counting of energy.\n+    if ( verbosity_ ) {\n+      std::cout << \"easy case: all tracks are associated to PFCandidates --> leaving tau momentum untouched.\" << std::endl;\n+    }\n+    updateTauP4(tau, 1., addNeutralsSumP4);\n+    return;\n+  } else {\n+    // This is the difficult case: \n+    // The tau energy needs to be computed for an arbitrary mix of charged and neutral PFCandidates plus reco::Tracks.\n+    // We need to make sure not to double-count energy deposited by reco::Track in ECAL and/or HCAL as neutral PFCandidates.\n+    \n+    // Check if we have enough energy in collection of PFNeutralHadrons and PFGammas that are not \"used\" by tau decay mode object\n+    // to balance track momenta:\n+    if ( nonPFCandTracksSumP < addNeutralsSumP4.energy() ) {\n+      double scaleFactor = 1. - nonPFCandTracksSumP/addNeutralsSumP4.energy();\n+      if ( !(scaleFactor >= 0. && scaleFactor <= 1.) ) {\n+\tedm::LogWarning(\"PFRecoTauEnergyAlgorithmPlugin::operator()\") \n+\t  << \"Failed to compute tau energy --> killing tau candidate !!\" << std::endl;\n+\tkillTau(tau);\n+\treturn;\n+      }\n+      if ( verbosity_ ) {\n+\tstd::cout << \"case (2): addNeutralsSumEn > nonPFCandTracksSumP --> adjusting tau momentum.\" << std::endl;\n+      }\n+      updateTauP4(tau, scaleFactor, addNeutralsSumP4);\n+      return;\n+    }\n+\n+    // Determine which neutral PFCandidates are close to PFChargedHadrons\n+    // and have been merged into ChargedHadrons\n+    std::vector<reco::PFCandidatePtr> mergedNeutrals;\n+    reco::Candidate::LorentzVector mergedNeutralsSumP4;\n+    for ( std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();\n+\t  chargedHadron != chargedHadrons.end(); ++chargedHadron ) {\n+      if ( chargedHadron->algoIs(PFRecoTauChargedHadron::kTrack) ) {\n+\tconst std::vector<reco::PFCandidatePtr>& neutralPFCands = chargedHadron->getNeutralPFCandidates();\n+\tfor ( std::vector<reco::PFCandidatePtr>::const_iterator neutralPFCand = neutralPFCands.begin();\n+\t      neutralPFCand != neutralPFCands.end(); ++neutralPFCand ) {\n+\t  mergedNeutrals.push_back(*neutralPFCand);\n+\t  mergedNeutralsSumP4 += (*neutralPFCand)->p4();\n+\t}\n+      }\n+    }\n+    if ( verbosity_ ) {\n+      std::cout << \"mergedNeutralsSumP4: En = \" << mergedNeutralsSumP4.energy() << std::endl;\n+    }\n+\n+    // Check if track momenta are balanced by sum of PFNeutralHadrons and PFGammas that are not \"used\" by tau decay mode object\n+    // plus neutral PFCandidates close to PFChargedHadrons:\n+    if ( nonPFCandTracksSumP < (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy()) ) {\n+      double scaleFactor = ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy()) - nonPFCandTracksSumP)/mergedNeutralsSumP4.energy();\n+      if ( !(scaleFactor >= 0. && scaleFactor <= 1.) ) {\n+      \tedm::LogWarning(\"PFRecoTauEnergyAlgorithmPlugin::operator()\") \n+\t  << \"Failed to compute tau energy --> killing tau candidate !!\" << std::endl;\n+\tkillTau(tau);\n+\treturn;\n+      }\n+      reco::Candidate::LorentzVector diffP4;\n+      size_t numChargedHadrons = chargedHadrons.size();\n+      for ( size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron ) {\n+\tconst PFRecoTauChargedHadron& chargedHadron = chargedHadrons.at(iChargedHadron);\n+\tif ( chargedHadron.getNeutralPFCandidates().size() >= 1 ) {\n+\t  PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;\n+\t  setChargedHadronP4(chargedHadron_modified, scaleFactor);\n+\t  tau.signalTauChargedHadronCandidates_[iChargedHadron] = chargedHadron_modified;\n+\t  diffP4 += (chargedHadron.p4() - chargedHadron_modified.p4());\n+\t}\n+      }\n+      if ( verbosity_ ) {\n+\tstd::cout << \"case (3): (addNeutralsSumEn + mergedNeutralsSumEn) > nonPFCandTracksSumP --> adjusting tau momentum.\" << std::endl;\n+      }\n+      updateTauP4(tau, -1., diffP4);\n+      return;\n+    }\n+\n+    // Determine energy sum of all PFNeutralHadrons interpreted as ChargedHadrons with missing track\n+    unsigned numChargedHadronNeutrals = 0;\n+    std::vector<reco::PFCandidatePtr> chargedHadronNeutrals;\n+    reco::Candidate::LorentzVector chargedHadronNeutralsSumP4;\n+    for ( std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();\n+\t  chargedHadron != chargedHadrons.end(); ++chargedHadron ) {\n+      if ( chargedHadron->algoIs(PFRecoTauChargedHadron::kPFNeutralHadron) ) {\n+\t++numChargedHadronNeutrals;\n+\tchargedHadronNeutrals.push_back(chargedHadron->getChargedPFCandidate());\n+\tchargedHadronNeutralsSumP4 += chargedHadron->getChargedPFCandidate()->p4();\n+      }\n+    }\n+    if ( verbosity_ ) {\n+      std::cout << \"chargedHadronNeutralsSumP4: En = \" << chargedHadronNeutralsSumP4.energy() \n+\t\t<< \" (numChargedHadronNeutrals = \" << numChargedHadronNeutrals << \")\" << std::endl;\n+    }\n+    \n+    // Check if sum of PFNeutralHadrons and PFGammas that are not \"used\" by tau decay mode object\n+    // plus neutral PFCandidates close to PFChargedHadrons plus PFNeutralHadrons interpreted as ChargedHadrons with missing track balances track momenta\n+    if ( nonPFCandTracksSumP < (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy()) ) {\n+      double scaleFactor = ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy()) - nonPFCandTracksSumP)/chargedHadronNeutralsSumP4.energy();\n+      if ( !(scaleFactor >= 0. && scaleFactor <= 1.) ) {\n+      \tedm::LogWarning(\"PFRecoTauEnergyAlgorithmPlugin::operator()\") \n+\t  << \"Failed to compute tau energy --> killing tau candidate !!\" << std::endl;\n+\tkillTau(tau);\n+\treturn;\n+      }\n+      reco::Candidate::LorentzVector diffP4;\n+      size_t numChargedHadrons = chargedHadrons.size();\n+      for ( size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron ) {\n+\tconst PFRecoTauChargedHadron& chargedHadron = chargedHadrons.at(iChargedHadron);\n+\tif ( chargedHadron.algoIs(PFRecoTauChargedHadron::kPFNeutralHadron) ) {\n+\t  PFRecoTauChargedHadron chargedHadron_modified = chargedHadron;\n+\t  chargedHadron_modified.neutralPFCandidates_.clear();\n+\t  const PFCandidatePtr& chargedPFCand = chargedHadron.getChargedPFCandidate();\n+\t  double chargedHadronP_modified     = scaleFactor*chargedPFCand->p();\n+\t  double chargedHadronTheta_modified = chargedPFCand->theta();\n+\t  double chargedHadronPhi_modified   = chargedPFCand->phi();\t\n+\t  reco::Candidate::LorentzVector chargedHadronP4_modified = compChargedHadronP4(chargedHadronP_modified, chargedHadronTheta_modified, chargedHadronPhi_modified);\t  ","path":"RecoTauTag/RecoTau/plugins/PFRecoTauEnergyAlgorithmPlugin.cc","position":276,"original_position":276,"commit_id":"3f03da5f073c5ddb07fbcb8877d55c8e2bd3e8dd","original_commit_id":"3f03da5f073c5ddb07fbcb8877d55c8e2bd3e8dd","user":{"login":"slava77","id":4676718,"avatar_url":"https://gravatar.com/avatar/c5a167808fdc334783250d6c3e861364?d=https%3A%2F%2Fidenticons.github.com%2F3a641d68584204b82a8b9449d36c83f9.png&r=x","gravatar_id":"c5a167808fdc334783250d6c3e861364","url":"https://api.github.com/users/slava77","html_url":"https://github.com/slava77","followers_url":"https://api.github.com/users/slava77/followers","following_url":"https://api.github.com/users/slava77/following{/other_user}","gists_url":"https://api.github.com/users/slava77/gists{/gist_id}","starred_url":"https://api.github.com/users/slava77/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/slava77/subscriptions","organizations_url":"https://api.github.com/users/slava77/orgs","repos_url":"https://api.github.com/users/slava77/repos","events_url":"https://api.github.com/users/slava77/events{/privacy}","received_events_url":"https://api.github.com/users/slava77/received_events","type":"User","site_admin":false},"body":"add a new method to set from px,py,pz. polar-cartesian ping pong is not good","created_at":"2014-01-18T14:42:56Z","updated_at":"2014-01-18T14:42:56Z","html_url":"https://github.com/cms-sw/cmssw/pull/1817#discussion_r8987917","pull_request_url":"https://api.github.com/repos/cms-sw/cmssw/pulls/1817","_links":{"self":{"href":"https://api.github.com/repos/cms-sw/cmssw/pulls/comments/8987917"},"html":{"href":"https://github.com/cms-sw/cmssw/pull/1817#discussion_r8987917"},"pull_request":{"href":"https://api.github.com/repos/cms-sw/cmssw/pulls/1817"}}}}