Merge branch 'HEP-FCC:main' into gmarchio-main-20240719

RecCaloCommon/include/RecCaloCommon/ClusterJet.h

@@ -3,25 +3,23 @@
 
 // std
 #include <map>
+#include <string>
 #include <vector>
 
-// Gaudi
-#include "GaudiAlg/GaudiAlgorithm.h"
-
+// FastJet
+#include "fastjet/ClusterSequence.hh"
 #include "fastjet/JetDefinition.hh"
-
-// EDM4hep
-#include "edm4hep/ClusterCollection.h"
+#include "fastjet/PseudoJet.hh"
 
 namespace k4::recCalo {
 
 class ClusterInfo : public fastjet::PseudoJet::UserInfoBase {
 public:
-  ClusterInfo(const int &index) : _index(index) {}
+  explicit ClusterInfo(const int &index) : _index(index) {}
   int index() const { return _index; }
 
 protected:
-  int _index;
+  int _index = 0;
 };
 
 /** @class ClusterJet
@@ -46,17 +44,17 @@ class ClusterInfo : public fastjet::PseudoJet::UserInfoBase {
 
 class ClusterJet {
 public:
-  ClusterJet(std::string jetAlg, double jetRadius,
+  ClusterJet(const std::string &jetAlg, double jetRadius,
              int isExclusiveClustering = 0, double minPt = 0,
              int clusterArgs = 0);
-  StatusCode initialize();
+  bool initialize();
   ~ClusterJet() {
-    if (m_clustSeq)
-      delete m_clustSeq;
+    delete m_jetDef;
+    delete m_clustSeq;
   }
 
   std::vector<fastjet::PseudoJet>
-  cluster(const std::vector<fastjet::PseudoJet> clustersPJ);
+  cluster(const std::vector<fastjet::PseudoJet> &clustersPJ);
 
 private:
   std::map<std::string, fastjet::JetAlgorithm> m_jetAlgMap = {
@@ -74,7 +72,8 @@ class ClusterJet {
   double m_minPt = 10;             // Only relevant for inclusive clustering
   int m_njets = 0;                 // Only relevant for exclusive clustering
 
-  fastjet::ClusterSequence *m_clustSeq;
+  fastjet::JetDefinition *m_jetDef = nullptr;
+  fastjet::ClusterSequence *m_clustSeq = nullptr;
 };
 
 } /* namespace k4::recCalo */

RecCaloCommon/src/ClusterJet.cpp

@@ -1,41 +1,44 @@
 #include "RecCaloCommon/ClusterJet.h"
 
+// std
+#include <iostream>
+
 namespace k4::recCalo {
 
-ClusterJet::ClusterJet(std::string jetAlg, double jetRadius,
+ClusterJet::ClusterJet(const std::string &jetAlg, double jetRadius,
                        int isExclusiveClustering, double minPt, int njets)
     : m_jetAlg(jetAlg), m_jetRadius(jetRadius),
       m_isExclusiveClustering(isExclusiveClustering), m_minPt(minPt),
-      m_njets(njets) {
-  m_clustSeq = nullptr;
-}
+      m_njets(njets) {}
 
-StatusCode ClusterJet::initialize() {
+bool ClusterJet::initialize() {
   if (m_jetAlgMap.find(m_jetAlg) == m_jetAlgMap.end()) {
     std::cout << "ERROR: "
               << " is not in the list of supported jet algorithms" << std::endl;
-    ;
-    return StatusCode::FAILURE;
+    return false;
   }
+
+  m_jetDef = new fastjet::JetDefinition(m_jetAlgMap.at(m_jetAlg), m_jetRadius);
+
   if (m_isExclusiveClustering > 1) {
     std::cout << "ERROR: "
               << "Clustering of " << m_isExclusiveClustering
               << " is currently not supported" << std::endl;
     ;
-    return StatusCode::FAILURE;
+    return false;
   }
-  return StatusCode::SUCCESS;
+
+  return true;
 }
 
 std::vector<fastjet::PseudoJet>
-ClusterJet::cluster(const std::vector<fastjet::PseudoJet> clustersPJ) {
+ClusterJet::cluster(const std::vector<fastjet::PseudoJet> &clustersPJ) {
+  // Deleting cluster sequence from previous event
+  delete m_clustSeq;
+
   std::vector<fastjet::PseudoJet> jets;
 
-  fastjet::JetDefinition *jetDef =
-      new fastjet::JetDefinition(m_jetAlgMap.at(m_jetAlg), m_jetRadius);
-  if (m_clustSeq)
-    delete m_clustSeq;
-  m_clustSeq = new fastjet::ClusterSequence(clustersPJ, *jetDef);
+  m_clustSeq = new fastjet::ClusterSequence(clustersPJ, *m_jetDef);
 
   // Note: initialize has already checked if m_isExclusiveClustering has the
   // right range

RecCalorimeter/src/components/CreateCaloJet.cpp

@@ -1,20 +1,22 @@
 // Gaudi
 #include "Gaudi/Property.h"
+#include <GaudiKernel/StatusCode.h>
 
 // k4FWCore
 #include "k4FWCore/Transformer.h"
 
-// std
-#include <math.h>
-#include <vector>
-
 // EDM4hep
 #include "edm4hep/ClusterCollection.h"
 #include "edm4hep/ReconstructedParticleCollection.h"
 
 // k4RecCalorimeter
 #include "RecCaloCommon/ClusterJet.h"
 
+// std
+#include <math.h>
+#include <string>
+#include <vector>
+
 /** @class CreateCaloJet
  * k4RecCalorimeter/RecFCCeeCalorimeter/src/components/CreateCaloJet.h
  *
@@ -41,11 +43,15 @@ struct CreateCaloJet final
             {KeyValues("InputClusterCollection", {"CorrectedCaloClusters"})},
             {KeyValues("OutputJetCollection", {"Jets"})}) {}
 
-  StatusCode initialize() {
-    clusterer = new k4::recCalo::ClusterJet(m_jetAlg, m_jetRadius,
-                                            m_isExclusive, m_minPt);
+  StatusCode initialize() override {
+    m_clusterer = new k4::recCalo::ClusterJet(m_jetAlg, m_jetRadius,
+                                              m_isExclusive, m_minPt);
+
+    if (!m_clusterer->initialize()) {
+      return StatusCode::FAILURE;
+    }
 
-    return clusterer->initialize();
+    return StatusCode::SUCCESS;
   }
 
   edm4hep::ReconstructedParticleCollection
@@ -75,7 +81,7 @@ struct CreateCaloJet final
     }
 
     std::vector<fastjet::PseudoJet> inclusiveJets =
-        clusterer->cluster(clustersPJ);
+        m_clusterer->cluster(clustersPJ);
 
     edm4hep::ReconstructedParticleCollection edmJets =
         edm4hep::ReconstructedParticleCollection();
@@ -98,6 +104,12 @@ struct CreateCaloJet final
     return edmJets;
   }
 
+  StatusCode finalize() override {
+    delete m_clusterer;
+
+    return StatusCode::SUCCESS;
+  }
+
 private:
   Gaudi::Property<std::string> m_jetAlg{this, "JetAlg", "antikt",
                                         "Name of jet clustering algorithm"};
@@ -108,7 +120,7 @@ struct CreateCaloJet final
   Gaudi::Property<int> m_isExclusive{this, "IsExclusiveClustering", 0,
                                      "1 if exclusive, 0 if inclusive"};
 
-  k4::recCalo::ClusterJet *clusterer;
+  k4::recCalo::ClusterJet *m_clusterer = nullptr;
 };
 
 DECLARE_COMPONENT(CreateCaloJet)

RecCalorimeter/src/components/CreateTruthJet.cpp

@@ -1,5 +1,6 @@
 // Gaudi
 #include "Gaudi/Property.h"
+#include <GaudiKernel/StatusCode.h>
 
 // k4FWCore
 #include "k4FWCore/Transformer.h"
@@ -49,10 +50,14 @@ struct CreateTruthJet final
              KeyValues("OutputAssociationsCollection",
                        {"TruthJetParticleAssociations"})}) {}
 
-  StatusCode initialize() {
-    clusterer = new k4::recCalo::ClusterJet(m_jetAlg, m_jetRadius,
-                                            m_isExclusive, m_minPt);
-    return clusterer->initialize();
+  StatusCode initialize() override {
+    m_clusterer = new k4::recCalo::ClusterJet(m_jetAlg, m_jetRadius,
+                                              m_isExclusive, m_minPt);
+    if (!m_clusterer->initialize()) {
+      return StatusCode::FAILURE;
+    }
+
+    return StatusCode::SUCCESS;
   }
 
   std::tuple<edm4hep::ReconstructedParticleCollection,
@@ -71,7 +76,7 @@ struct CreateTruthJet final
     }
 
     std::vector<fastjet::PseudoJet> inclusiveJets =
-        clusterer->cluster(clustersPJ);
+        m_clusterer->cluster(clustersPJ);
 
     auto edmJets = edm4hep::ReconstructedParticleCollection();
     auto assoc = edm4hep::MCRecoParticleAssociationCollection();
@@ -99,8 +104,13 @@ struct CreateTruthJet final
     return std::make_tuple(std::move(edmJets), std::move(assoc));
   }
 
-  private:
+  StatusCode finalize() override {
+    delete m_clusterer;
+
+    return StatusCode::SUCCESS;
+  }
 
+private:
   Gaudi::Property<std::string> m_jetAlg{this, "JetAlg", "antikt",
                                         "Name of jet clustering algorithm"};
   Gaudi::Property<double> m_jetRadius{this, "JetRadius", 0.4,
@@ -109,7 +119,7 @@ struct CreateTruthJet final
                                   "Minimum pT for saved jets"};
   Gaudi::Property<bool> m_isExclusive{this, "isExclusiveClustering", false,
                                       "true if exclusive, false if inclusive"};
-  k4::recCalo::ClusterJet *clusterer;
+  k4::recCalo::ClusterJet *m_clusterer = nullptr;
 };
 
 DECLARE_COMPONENT(CreateTruthJet)

RecFCCeeCalorimeter/src/components/AugmentClustersFCCee.cpp

@@ -633,16 +633,33 @@ StatusCode AugmentClustersFCCee::execute([[maybe_unused]] const EventContext &ev
 
         // do pi0/photon shape var only for EMB
         if (m_do_photon_shapeVar && systemID == systemID_EMB) {
-	  double w_theta;
           if (m_do_widthTheta_logE_weights) {
-            w_theta = (sumWeightLayer[layer+startPositionToFill] != 0.) ? sqrt(theta2_E_layer[layer+startPositionToFill] / sumWeightLayer[layer+startPositionToFill] - std::pow(theta_E_layer[layer+startPositionToFill] / sumWeightLayer[layer+startPositionToFill], 2)) : 0. ;
+            double w_theta2 = theta2_E_layer[layer+startPositionToFill] / sumWeightLayer[layer+startPositionToFill] - std::pow(theta_E_layer[layer+startPositionToFill] / sumWeightLayer[layer+startPositionToFill], 2);
+            // Very small but negative values caused by computational precision are allowed,
+            // otherwise crash.
+            if (w_theta2 < -1e-9) {
+              error() << "w_theta2 in theta width calculation is negative: " << w_theta2 << endmsg;
+              return StatusCode::FAILURE;
+            }
+            width_theta[layer+startPositionToFill] = (sumWeightLayer[layer+startPositionToFill] > 0.) ? std::sqrt(std::fabs(w_theta2)) : 0. ;
           } else {
-            w_theta = (sumEnLayer[layer+startPositionToFill] > 0.) ? sqrt(theta2_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill] - std::pow(theta_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill], 2)) : 0. ;
+            double w_theta2 = theta2_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill] - std::pow(theta_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill], 2);
+            // Very small but negative values caused by computational precision are allowed,
+            // otherwise crash.
+            if (w_theta2 < -1e-9) {
+              error() << "w_theta2 in theta width calculation is negative: " << w_theta2 << endmsg;
+              return StatusCode::FAILURE;
+            }
+            width_theta[layer+startPositionToFill] = (sumEnLayer[layer+startPositionToFill] > 0.) ? std::sqrt(std::fabs(w_theta2)) : 0. ;
 	  }
-	  width_theta[layer+startPositionToFill] = w_theta;
-
-          double w_module = (sumEnLayer[layer+startPositionToFill] > 0.) ? sqrt(module2_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill] - std::pow(module_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill], 2)) : 0. ;
-          width_module[layer+startPositionToFill] = w_module;
+          double w_module2 = module2_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill] - std::pow(module_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill], 2);
+          // Very small but negative values caused by computational precision are allowed,
+          // otherwise crash.
+          if (w_module2 < -1e-9) {
+            error() << "w_module2 in module width calculation is negative: " << w_module2 << endmsg;
+            return StatusCode::FAILURE;
+          }
+          width_module[layer+startPositionToFill] = (sumEnLayer[layer+startPositionToFill] > 0.) ? std::sqrt(std::fabs(w_module2)) : 0. ;
 
           double Ratio_E = (E_cell_Max[layer+startPositionToFill] - E_cell_secMax[layer+startPositionToFill]) /
                            (E_cell_Max[layer+startPositionToFill] + E_cell_secMax[layer+startPositionToFill]);
@@ -766,15 +783,34 @@ StatusCode AugmentClustersFCCee::execute([[maybe_unused]] const EventContext &ev
               double theta2_E_9Bin = theta2_E_7Bin + theta_m4 * theta_m4 * weightLog_m4 + theta_p4 * theta_p4 * weightLog_p4;
               double theta_E_9Bin = theta_E_7Bin + theta_m4 * weightLog_m4 + theta_p4 * weightLog_p4;
 
-              double _w_theta_3Bin = sqrt(theta2_E_3Bin / sum_weightLog_3Bin - std::pow(theta_E_3Bin / sum_weightLog_3Bin, 2));
-              double _w_theta_5Bin = sqrt(theta2_E_5Bin / sum_weightLog_5Bin - std::pow(theta_E_5Bin / sum_weightLog_5Bin, 2));
-              double _w_theta_7Bin = sqrt(theta2_E_7Bin / sum_weightLog_7Bin - std::pow(theta_E_7Bin / sum_weightLog_7Bin, 2));
-              double _w_theta_9Bin = sqrt(theta2_E_9Bin / sum_weightLog_9Bin - std::pow(theta_E_9Bin / sum_weightLog_9Bin, 2));
-
-              width_theta_3Bin[layer+startPositionToFill] = (sum_weightLog_3Bin > 0.) ? _w_theta_3Bin : 0. ;
-              width_theta_5Bin[layer+startPositionToFill] = (sum_weightLog_5Bin > 0.) ? _w_theta_5Bin : 0. ;
-              width_theta_7Bin[layer+startPositionToFill] = (sum_weightLog_7Bin > 0.) ? _w_theta_7Bin : 0. ;
-              width_theta_9Bin[layer+startPositionToFill] = (sum_weightLog_9Bin > 0.) ? _w_theta_9Bin : 0. ;
+              double _w_theta_3Bin2 = theta2_E_3Bin / sum_weightLog_3Bin - std::pow(theta_E_3Bin / sum_weightLog_3Bin, 2);
+              double _w_theta_5Bin2 = theta2_E_5Bin / sum_weightLog_5Bin - std::pow(theta_E_5Bin / sum_weightLog_5Bin, 2);
+              double _w_theta_7Bin2 = theta2_E_7Bin / sum_weightLog_7Bin - std::pow(theta_E_7Bin / sum_weightLog_7Bin, 2);
+              double _w_theta_9Bin2 = theta2_E_9Bin / sum_weightLog_9Bin - std::pow(theta_E_9Bin / sum_weightLog_9Bin, 2);
+
+              // Very small but negative values caused by computational precision are allowed,
+              // otherwise crash.
+              if (_w_theta_3Bin2 < -1e-9) {
+                error() << "_w_theta_3Bin2 in theta width calculation is negative: " << _w_theta_3Bin2 << endmsg;
+                return StatusCode::FAILURE;
+              }
+              if (_w_theta_5Bin2 < -1e-9) {
+                error() << "_w_theta_5Bin2 in theta width calculation is negative: " << _w_theta_5Bin2 << endmsg;
+                return StatusCode::FAILURE;
+              }
+              if (_w_theta_7Bin2 < -1e-9) {
+                error() << "_w_theta_7Bin2 in theta width calculation is negative: " << _w_theta_7Bin2 << endmsg;
+                return StatusCode::FAILURE;
+              }
+              if (_w_theta_9Bin2 < -1e-9) {
+                error() << "_w_theta_9Bin2 in theta width calculation is negative: " << _w_theta_9Bin2 << endmsg;
+                return StatusCode::FAILURE;
+              }
+
+              width_theta_3Bin[layer+startPositionToFill] = (sum_weightLog_3Bin > 0.) ? std::sqrt(std::fabs(_w_theta_3Bin2)) : 0. ;
+              width_theta_5Bin[layer+startPositionToFill] = (sum_weightLog_5Bin > 0.) ? std::sqrt(std::fabs(_w_theta_5Bin2)) : 0. ;
+              width_theta_7Bin[layer+startPositionToFill] = (sum_weightLog_7Bin > 0.) ? std::sqrt(std::fabs(_w_theta_7Bin2)) : 0. ;
+              width_theta_9Bin[layer+startPositionToFill] = (sum_weightLog_9Bin > 0.) ? std::sqrt(std::fabs(_w_theta_9Bin2)) : 0. ;
             } else {
               double theta2_E_3Bin = theta_m1 * theta_m1 * E_m1 + E_cell_Max_theta[layer+startPositionToFill] * E_cell_Max_theta[layer+startPositionToFill] * E_cell_Max[layer+startPositionToFill] + theta_p1 * theta_p1 * E_p1;
               double theta_E_3Bin = theta_m1 * E_m1 + E_cell_Max_theta[layer+startPositionToFill] * E_cell_Max[layer+startPositionToFill] + theta_p1 * E_p1;
@@ -785,15 +821,34 @@ StatusCode AugmentClustersFCCee::execute([[maybe_unused]] const EventContext &ev
               double theta2_E_9Bin = theta2_E_7Bin + theta_m4 * theta_m4 * E_m4 + theta_p4 * theta_p4 * E_p4;
               double theta_E_9Bin = theta_E_7Bin + theta_m4 * E_m4 + theta_p4 * E_p4;
 
-              double _w_theta_3Bin = sqrt(theta2_E_3Bin / sum_E_3Bin - std::pow(theta_E_3Bin / sum_E_3Bin, 2));
-              double _w_theta_5Bin = sqrt(theta2_E_5Bin / sum_E_5Bin - std::pow(theta_E_5Bin / sum_E_5Bin, 2));
-              double _w_theta_7Bin = sqrt(theta2_E_7Bin / sum_E_7Bin - std::pow(theta_E_7Bin / sum_E_7Bin, 2));
-              double _w_theta_9Bin = sqrt(theta2_E_9Bin / sum_E_9Bin - std::pow(theta_E_9Bin / sum_E_9Bin, 2));
-
-              width_theta_3Bin[layer+startPositionToFill] = (sum_E_3Bin > 0.) ? _w_theta_3Bin : 0. ;
-              width_theta_5Bin[layer+startPositionToFill] = (sum_E_5Bin > 0.) ? _w_theta_5Bin : 0. ;
-              width_theta_7Bin[layer+startPositionToFill] = (sum_E_7Bin > 0.) ? _w_theta_7Bin : 0. ;
-              width_theta_9Bin[layer+startPositionToFill] = (sum_E_9Bin > 0.) ? _w_theta_9Bin : 0. ;
+              double _w_theta_3Bin2 = theta2_E_3Bin / sum_E_3Bin - std::pow(theta_E_3Bin / sum_E_3Bin, 2);
+              double _w_theta_5Bin2 = theta2_E_5Bin / sum_E_5Bin - std::pow(theta_E_5Bin / sum_E_5Bin, 2);
+              double _w_theta_7Bin2 = theta2_E_7Bin / sum_E_7Bin - std::pow(theta_E_7Bin / sum_E_7Bin, 2);
+              double _w_theta_9Bin2 = theta2_E_9Bin / sum_E_9Bin - std::pow(theta_E_9Bin / sum_E_9Bin, 2);
+
+              // Very small but negative values caused by computational precision are allowed,
+              // otherwise crash.
+              if (_w_theta_3Bin2 < -1e-9) {
+                error() << "_w_theta_3Bin2 in theta width calculation is negative: " << _w_theta_3Bin2 << endmsg;
+                return StatusCode::FAILURE;
+              }
+              if (_w_theta_5Bin2 < -1e-9) {
+                error() << "_w_theta_5Bin2 in theta width calculation is negative: " << _w_theta_5Bin2 << endmsg;
+                return StatusCode::FAILURE;
+              }
+              if (_w_theta_7Bin2 < -1e-9) {
+                error() << "_w_theta_7Bin2 in theta width calculation is negative: " << _w_theta_7Bin2 << endmsg;
+                return StatusCode::FAILURE;
+              }
+              if (_w_theta_9Bin2 < -1e-9) {
+                error() << "_w_theta_9Bin2 in theta width calculation is negative: " << _w_theta_9Bin2 << endmsg;
+                return StatusCode::FAILURE;
+              }
+
+              width_theta_3Bin[layer+startPositionToFill] = (sum_E_3Bin > 0.) ? std::sqrt(std::fabs(_w_theta_3Bin2)) : 0. ;
+              width_theta_5Bin[layer+startPositionToFill] = (sum_E_5Bin > 0.) ? std::sqrt(std::fabs(_w_theta_5Bin2)) : 0. ;
+              width_theta_7Bin[layer+startPositionToFill] = (sum_E_7Bin > 0.) ? std::sqrt(std::fabs(_w_theta_7Bin2)) : 0. ;
+              width_theta_9Bin[layer+startPositionToFill] = (sum_E_9Bin > 0.) ? std::sqrt(std::fabs(_w_theta_9Bin2)) : 0. ;
             }
             E_fr_side_pm2[layer+startPositionToFill] = (sum_E_3Bin > 0.) ? (sum_E_5Bin / sum_E_3Bin - 1.) : 0. ;
             E_fr_side_pm3[layer+startPositionToFill] = (sum_E_3Bin > 0.) ? (sum_E_7Bin / sum_E_3Bin - 1.) : 0. ;