From 2e412a326342cdc4d6fc59311a4afcac56b3da43 Mon Sep 17 00:00:00 2001
From: Jason Veatch <Jason.Veatch@cern.ch>
Date: Thu, 12 Mar 2020 19:30:25 +0100
Subject: [PATCH] Adding C4 and some protections

---
 .../EnergyCorrelatorRatiosTool.h              |  1 +
 .../EnergyCorrelatorTool.h                    |  1 +
 .../EnergyCorrelatorGeneralizedRatiosTool.cxx |  4 ++--
 .../Root/EnergyCorrelatorGeneralizedTool.cxx  |  1 +
 .../Root/EnergyCorrelatorRatiosTool.cxx       | 24 ++++++++++++++++++-
 .../Root/EnergyCorrelatorTool.cxx             | 19 ++++++++++++---
 6 files changed, 44 insertions(+), 6 deletions(-)

diff --git a/Reconstruction/Jet/JetSubStructureMomentTools/JetSubStructureMomentTools/EnergyCorrelatorRatiosTool.h b/Reconstruction/Jet/JetSubStructureMomentTools/JetSubStructureMomentTools/EnergyCorrelatorRatiosTool.h
index cdd4e9b70705..4cecb2b18673 100644
--- a/Reconstruction/Jet/JetSubStructureMomentTools/JetSubStructureMomentTools/EnergyCorrelatorRatiosTool.h
+++ b/Reconstruction/Jet/JetSubStructureMomentTools/JetSubStructureMomentTools/EnergyCorrelatorRatiosTool.h
@@ -36,6 +36,7 @@ class EnergyCorrelatorRatiosTool :
 
     private:
       bool m_doC3;
+      bool m_doC4;
       std::vector<float> m_rawBetaVals; /// Vector of input values before cleaning
       std::vector<float> m_betaVals; /// Local vector for cleaned up inputs
       bool m_doDichroic;
diff --git a/Reconstruction/Jet/JetSubStructureMomentTools/JetSubStructureMomentTools/EnergyCorrelatorTool.h b/Reconstruction/Jet/JetSubStructureMomentTools/JetSubStructureMomentTools/EnergyCorrelatorTool.h
index 7a1cd2af3fe7..1182b5c61900 100644
--- a/Reconstruction/Jet/JetSubStructureMomentTools/JetSubStructureMomentTools/EnergyCorrelatorTool.h
+++ b/Reconstruction/Jet/JetSubStructureMomentTools/JetSubStructureMomentTools/EnergyCorrelatorTool.h
@@ -37,6 +37,7 @@ class EnergyCorrelatorTool :
     private:
       float m_Beta;
       bool m_doC3;
+      bool m_doC4;
       std::vector<float> m_rawBetaVals; /// Vector of input values before cleaning
       std::vector<float> m_betaVals; /// Local vector for cleaned up inputs
       bool m_doDichroic;
diff --git a/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorGeneralizedRatiosTool.cxx b/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorGeneralizedRatiosTool.cxx
index 8102c3d015f4..c0438a37449a 100644
--- a/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorGeneralizedRatiosTool.cxx
+++ b/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorGeneralizedRatiosTool.cxx
@@ -15,6 +15,7 @@ EnergyCorrelatorGeneralizedRatiosTool::EnergyCorrelatorGeneralizedRatiosTool(std
 }
 
 StatusCode EnergyCorrelatorGeneralizedRatiosTool::initialize() {
+
   // Add beta = 1.0 by default
   m_betaVals.push_back(1.0);
 
@@ -106,7 +107,6 @@ int EnergyCorrelatorGeneralizedRatiosTool::modifyJet(xAOD::Jet &jet) const {
     }
 
     // N2
-
     if(ecfg_2_1 > 1e-8) { // Prevent div-0
       jet.setAttribute(m_prefix+"N2"+suffix, ecfg_3_2 / (pow(ecfg_2_1, 2.0)));
 
@@ -121,7 +121,7 @@ int EnergyCorrelatorGeneralizedRatiosTool::modifyJet(xAOD::Jet &jet) const {
     }
 
     // N3
-    if(ecfg_3_1 > 1e-8) // Prevent div-0
+    if(m_doN3 && ecfg_3_1 > 1e-8) // Prevent div-0
       jet.setAttribute(m_prefix+"N3"+suffix, ecfg_4_2 / (pow(ecfg_3_1, 2.0)));
     else
       jet.setAttribute(m_prefix+"N3"+suffix, -999.0);
diff --git a/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorGeneralizedTool.cxx b/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorGeneralizedTool.cxx
index 5df6927d5356..8bff896e130c 100644
--- a/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorGeneralizedTool.cxx
+++ b/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorGeneralizedTool.cxx
@@ -17,6 +17,7 @@ EnergyCorrelatorGeneralizedTool::EnergyCorrelatorGeneralizedTool(std::string nam
 }
 
 StatusCode EnergyCorrelatorGeneralizedTool::initialize() {
+
   // Add beta = 1.0 by default
   m_betaVals.push_back(1.0);
 
diff --git a/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorRatiosTool.cxx b/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorRatiosTool.cxx
index 4e97b88f3f4a..67c6b03e30ee 100644
--- a/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorRatiosTool.cxx
+++ b/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorRatiosTool.cxx
@@ -10,10 +10,12 @@ EnergyCorrelatorRatiosTool::EnergyCorrelatorRatiosTool(std::string name) :
 {
   declareProperty("BetaList", m_rawBetaVals = {});
   declareProperty("DoC3",  m_doC3 = false);
+  declareProperty("DoC4",  m_doC4 = false);
   declareProperty("DoDichroic", m_doDichroic = false);
 }
 
 StatusCode EnergyCorrelatorRatiosTool::initialize() {
+
   // Add beta = 1.0 by default
   m_betaVals.push_back(1.0);
 
@@ -38,6 +40,10 @@ StatusCode EnergyCorrelatorRatiosTool::initialize() {
     ATH_MSG_DEBUG("Including beta = " << beta);
   }
 
+  // If DoC4 is set to true, set DoC3 to true by default since it won't
+  // add any additional computational overhead
+  if(m_doC4) m_doC3 = true;
+
   ATH_CHECK(JetSubStructureMomentToolsBase::initialize());
 
   return StatusCode::SUCCESS;
@@ -68,6 +74,11 @@ int EnergyCorrelatorRatiosTool::modifyJet(xAOD::Jet &jet) const {
       return 1;
     }
 
+    if (m_doC4 && !jet.isAvailable<float>(m_prefix+"ECF5"+suffix)) {
+      ATH_MSG_WARNING("Energy correlation function " << m_prefix << "ECF5" << suffix << " is not available. Exiting..");
+      return 1;
+    }
+
     if(m_doDichroic) {
       if (!jet.isAvailable<float>(m_prefix+"ECF1_ungroomed"+suffix)) {
         ATH_MSG_WARNING("Energy correlation function " << m_prefix << "ECF1_ungroomed" << suffix << " is not available. Exiting..");
@@ -94,6 +105,11 @@ int EnergyCorrelatorRatiosTool::modifyJet(xAOD::Jet &jet) const {
       ecf4 = jet.getAttribute<float>(m_prefix+"ECF4"+suffix);
     }
 
+    float ecf5 = -999.0;
+    if(m_doC4) {
+      ecf5 = jet.getAttribute<float>(m_prefix+"ECF5"+suffix);
+    }
+
     float ecf1_ungroomed = -999.0;
     float ecf2_ungroomed = -999.0;
     float ecf3_ungroomed = -999.0;
@@ -132,10 +148,16 @@ int EnergyCorrelatorRatiosTool::modifyJet(xAOD::Jet &jet) const {
       jet.setAttribute(m_prefix+"C2"+suffix, -999.0);
 
     // C3
-    if(ecf3 > 1e-8) // Prevent div-0
+    if(m_doC3 && ecf3 > 1e-8) // Prevent div-0
       jet.setAttribute(m_prefix+"C3"+suffix, ecf4 * ecf2 / pow(ecf3, 2.0));
     else
       jet.setAttribute(m_prefix+"C3"+suffix, -999.0);
+
+    // C4
+    if(m_doC4 && ecf4 > 1e-8) // Prevent div-0
+      jet.setAttribute(m_prefix+"C4"+suffix, ecf5 * ecf3 / pow(ecf4, 2.0));
+    else
+      jet.setAttribute(m_prefix+"C4"+suffix, -999.0);
   }
 
   return 0;
diff --git a/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorTool.cxx b/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorTool.cxx
index 0f5cafb6f300..23cc78c2d22f 100644
--- a/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorTool.cxx
+++ b/Reconstruction/Jet/JetSubStructureMomentTools/Root/EnergyCorrelatorTool.cxx
@@ -11,10 +11,12 @@ EnergyCorrelatorTool::EnergyCorrelatorTool(std::string name) :
   declareProperty("Beta", m_Beta = 1.0);
   declareProperty("BetaList", m_rawBetaVals = {});
   declareProperty("DoC3", m_doC3 = false);
+  declareProperty("DoC4", m_doC4 = false);
   declareProperty("DoDichroic", m_doDichroic = false);
 }
 
 StatusCode EnergyCorrelatorTool::initialize() {
+
   // Add beta = 1.0 by default
   m_betaVals.push_back(1.0);
 
@@ -42,13 +44,17 @@ StatusCode EnergyCorrelatorTool::initialize() {
     ATH_MSG_DEBUG("Including beta = " << beta);
   }
 
+  // If DoC4 is set to true, set DoC3 to true by default since it won't
+  // add any additional computational overhead
+  if(m_doC4) m_doC3 = true;
+
   ATH_CHECK(JetSubStructureMomentToolsBase::initialize());
 
   return StatusCode::SUCCESS;
 }
 
 int EnergyCorrelatorTool::modifyJet(xAOD::Jet &injet) const {
-  
+
   fastjet::PseudoJet jet;
   fastjet::PseudoJet jet_ungroomed;
 
@@ -81,6 +87,7 @@ int EnergyCorrelatorTool::modifyJet(xAOD::Jet &injet) const {
     float result_ECF2 = -999;
     float result_ECF3 = -999;
     float result_ECF4 = -999;
+    float result_ECF5 = -999;
 
     float result_ECF1_ungroomed = -999;
     float result_ECF2_ungroomed = -999;
@@ -100,7 +107,12 @@ int EnergyCorrelatorTool::modifyJet(xAOD::Jet &injet) const {
         JetSubStructureUtils::EnergyCorrelator ECF4(4, beta, JetSubStructureUtils::EnergyCorrelator::pt_R);
         result_ECF4 = ECF4.result(jet);
       }
-    
+
+      if(m_doC4) {
+        JetSubStructureUtils::EnergyCorrelator ECF5(5, beta, JetSubStructureUtils::EnergyCorrelator::pt_R);
+        result_ECF5 = ECF5.result(jet);
+      }
+
       if(decorate_ungroomed) {
         result_ECF1_ungroomed = ECF1.result(jet_ungroomed);
         result_ECF2_ungroomed = ECF2.result(jet_ungroomed);
@@ -113,7 +125,8 @@ int EnergyCorrelatorTool::modifyJet(xAOD::Jet &injet) const {
     injet.setAttribute(m_prefix+"ECF2"+suffix, result_ECF2);
     injet.setAttribute(m_prefix+"ECF3"+suffix, result_ECF3);
     injet.setAttribute(m_prefix+"ECF4"+suffix, result_ECF4);
-    
+    injet.setAttribute(m_prefix+"ECF5"+suffix, result_ECF5);
+
     injet.setAttribute(m_prefix+"ECF1_ungroomed"+suffix, result_ECF1_ungroomed);
     injet.setAttribute(m_prefix+"ECF2_ungroomed"+suffix, result_ECF2_ungroomed);
     injet.setAttribute(m_prefix+"ECF3_ungroomed"+suffix, result_ECF3_ungroomed);
-- 
GitLab