Add theta-energy calibration parameters

k4GaudiPandora/include/DDPandoraPFANewAlgorithm.h

@@ -87,6 +87,20 @@ struct DDPandoraPFANewAlgorithm final
     std::vector<float> m_ecalOutputEnergyCorrectionPoints{}; ///< The input energy points for non-linearity energy
                                                              ///< correction in the ECAL
 
+    // Theta-energy (2D) calibration, EM branch
+    bool m_electromagneticThetaEnergyCorrectionEnabled = false;
+    std::string m_electromagneticThetaEnergyCorrectionPluginName = "PhotonEMNonLinearity";
+    std::vector<float> m_electromagneticThetaEnergyCorrectionThetaBinEdges{};
+    std::vector<float> m_electromagneticThetaEnergyCorrectionEnergyBinEdges{};
+    std::vector<float> m_electromagneticThetaEnergyCorrectionScaleFactors{};
+
+    // Theta-energy (2D) calibration, HAD branch
+    bool m_hadronicThetaEnergyCorrectionEnabled = false;
+    std::string m_hadronicThetaEnergyCorrectionPluginName = "HadronicThetaEnergyBinned";
+    std::vector<float> m_hadronicThetaEnergyCorrectionThetaBinEdges{};
+    std::vector<float> m_hadronicThetaEnergyCorrectionEnergyBinEdges{};
+    std::vector<float> m_hadronicThetaEnergyCorrectionScaleFactors{};
+
     // Software compensation parameters
     std::vector<float> m_softCompParameters{};
     std::vector<float> m_softCompEnergyDensityBins{};
@@ -159,6 +173,9 @@ struct DDPandoraPFANewAlgorithm final
    */
   void finaliseSteeringParameters();
 
+  /// Check the theta-energy calibration tables are well formed. Logs and returns false if not.
+  bool validateThetaEnergyCorrectionSettings() const;
+
   /**
    *  @brief  Reset the pandora pfa new processor
    */
@@ -359,6 +376,40 @@ struct DDPandoraPFANewAlgorithm final
       this, "ECALInputEnergyCorrectionPoints", {}, "The input energy points for electromagnetic energy correction"};
   Gaudi::Property<std::vector<float>> m_ecalOutputEnergyCorrectionPoints{
       this, "ECALOutputEnergyCorrectionPoints", {}, "The output energy points for electromagnetic energy correction"};
+
+  // Theta-energy (2D) calibration, EM branch
+  Gaudi::Property<bool> m_emThetaEnergyEnabled{
+      this, "ElectromagneticThetaEnergyCorrectionEnabled", false,
+      "Whether to enable theta-energy correction on the Pandora EM branch"};
+  Gaudi::Property<std::string> m_emThetaEnergyPluginName{
+      this, "ElectromagneticThetaEnergyCorrectionPluginName", std::string("PhotonEMNonLinearity"),
+      "The name of the Pandora EM theta-energy correction plugin"};
+  Gaudi::Property<std::vector<float>> m_emThetaEnergyThetaBinEdges{
+      this, "ElectromagneticThetaEnergyCorrectionThetaBinEdges", {},
+      "The theta bin edges for Pandora EM theta-energy correction"};
+  Gaudi::Property<std::vector<float>> m_emThetaEnergyEnergyBinEdges{
+      this, "ElectromagneticThetaEnergyCorrectionEnergyBinEdges", {},
+      "The energy bin edges for Pandora EM theta-energy correction"};
+  Gaudi::Property<std::vector<float>> m_emThetaEnergyScaleFactors{
+      this, "ElectromagneticThetaEnergyCorrectionScaleFactors", {},
+      "The row-major theta-energy scale factors for Pandora EM theta-energy correction"};
+
+  // Theta-energy (2D) calibration, HAD branch
+  Gaudi::Property<bool> m_hadThetaEnergyEnabled{
+      this, "HadronicThetaEnergyCorrectionEnabled", false,
+      "Whether to enable theta-energy correction on the Pandora HAD branch"};
+  Gaudi::Property<std::string> m_hadThetaEnergyPluginName{
+      this, "HadronicThetaEnergyCorrectionPluginName", std::string("HadronicThetaEnergyBinned"),
+      "The name of the Pandora HAD theta-energy correction plugin"};
+  Gaudi::Property<std::vector<float>> m_hadThetaEnergyThetaBinEdges{
+      this, "HadronicThetaEnergyCorrectionThetaBinEdges", {},
+      "The theta bin edges for Pandora HAD theta-energy correction"};
+  Gaudi::Property<std::vector<float>> m_hadThetaEnergyEnergyBinEdges{
+      this, "HadronicThetaEnergyCorrectionEnergyBinEdges", {},
+      "The energy bin edges for Pandora HAD theta-energy correction"};
+  Gaudi::Property<std::vector<float>> m_hadThetaEnergyScaleFactors{
+      this, "HadronicThetaEnergyCorrectionScaleFactors", {},
+      "The row-major theta-energy scale factors for Pandora HAD theta-energy correction"};
   // EXTRA PARAMETERS FROM NIKIFOROS m_caloEncodingString
   Gaudi::Property<std::string> m_trackCreatorName{this, "TrackCreatorName", "DDTrackCreatorCLIC",
                                                   "The name of the DDTrackCreator implementation"};

k4GaudiPandora/src/DDPandoraPFANewAlgorithm.cpp

@@ -39,6 +39,30 @@
 #include <LCContent.h>
 #include <LCPlugins/LCSoftwareCompensation.h>
 
+namespace {
+
+bool isStrictlyIncreasing(const std::vector<float>& values) {
+  if (values.size() < 2)
+    return false;
+  for (std::size_t i = 1; i < values.size(); ++i) {
+    if (values[i] <= values[i - 1])
+      return false;
+  }
+  return true;
+}
+
+bool isConsistentThetaEnergyTable(const std::vector<float>& thetaBinEdges,
+                                  const std::vector<float>& energyBinEdges,
+                                  const std::vector<float>& scaleFactors) {
+  if (!isStrictlyIncreasing(thetaBinEdges) || !isStrictlyIncreasing(energyBinEdges))
+    return false;
+  const std::size_t nThetaBins = thetaBinEdges.size() - 1;
+  const std::size_t nEnergyBins = energyBinEdges.size() - 1;
+  return (nThetaBins * nEnergyBins == scaleFactors.size());
+}
+
+} // namespace
+
 #include <DD4hep/DD4hepUnits.h>
 #include <DD4hep/DetType.h>
 #include <DD4hep/Detector.h>
@@ -131,6 +155,9 @@ StatusCode DDPandoraPFANewAlgorithm::initialize() {
 
   finaliseSteeringParameters();
 
+  if (!validateThetaEnergyCorrectionSettings())
+    return StatusCode::FAILURE;
+
   if (m_settings.m_detectorName == "ALLEGRO") {
     m_settings.m_trackCreatorName = "DDTrackCreatorALLEGRO";
     m_geometryCreator = std::make_unique<DDGeometryCreatorALLEGRO>(m_geometryCreatorSettings, m_pPandora, this);
@@ -353,6 +380,41 @@ pandora::StatusCode DDPandoraPFANewAlgorithm::registerUserComponents() const {
                                m_pPandora, "NonLinearity", pandora::HADRONIC, m_settings.m_inputEnergyCorrectionPoints,
                                m_settings.m_outputEnergyCorrectionPoints))
 
+  // Theta-energy (2D) calibration. Both plugin names are registered unconditionally so that a
+  // Pandora settings XML naming them is always valid: without a payload they fall back to the
+  // 1D overload with empty points, which is the identity.
+  if (m_settings.m_electromagneticThetaEnergyCorrectionEnabled) {
+    PANDORA_RETURN_RESULT_IF(
+        pandora::STATUS_CODE_SUCCESS, !=,
+        LCContent::RegisterNonLinearityEnergyCorrection(
+            m_pPandora, m_settings.m_electromagneticThetaEnergyCorrectionPluginName, pandora::ELECTROMAGNETIC,
+            m_settings.m_electromagneticThetaEnergyCorrectionThetaBinEdges,
+            m_settings.m_electromagneticThetaEnergyCorrectionEnergyBinEdges,
+            m_settings.m_electromagneticThetaEnergyCorrectionScaleFactors))
+  } else {
+    PANDORA_RETURN_RESULT_IF(
+        pandora::STATUS_CODE_SUCCESS, !=,
+        LCContent::RegisterNonLinearityEnergyCorrection(
+            m_pPandora, m_settings.m_electromagneticThetaEnergyCorrectionPluginName, pandora::ELECTROMAGNETIC,
+            std::vector<float>(), std::vector<float>()))
+  }
+
+  if (m_settings.m_hadronicThetaEnergyCorrectionEnabled) {
+    PANDORA_RETURN_RESULT_IF(
+        pandora::STATUS_CODE_SUCCESS, !=,
+        LCContent::RegisterNonLinearityEnergyCorrection(
+            m_pPandora, m_settings.m_hadronicThetaEnergyCorrectionPluginName, pandora::HADRONIC,
+            m_settings.m_hadronicThetaEnergyCorrectionThetaBinEdges,
+            m_settings.m_hadronicThetaEnergyCorrectionEnergyBinEdges,
+            m_settings.m_hadronicThetaEnergyCorrectionScaleFactors))
+  } else {
+    PANDORA_RETURN_RESULT_IF(
+        pandora::STATUS_CODE_SUCCESS, !=,
+        LCContent::RegisterNonLinearityEnergyCorrection(
+            m_pPandora, m_settings.m_hadronicThetaEnergyCorrectionPluginName, pandora::HADRONIC,
+            std::vector<float>(), std::vector<float>()))
+  }
+
   lc_content::LCSoftwareCompensationParameters softwareCompensationParameters;
   softwareCompensationParameters.m_softCompParameters = m_settings.m_softCompParameters;
   softwareCompensationParameters.m_softCompEnergyDensityBins = m_settings.m_softCompEnergyDensityBins;
@@ -369,6 +431,34 @@ pandora::StatusCode DDPandoraPFANewAlgorithm::registerUserComponents() const {
   return pandora::STATUS_CODE_SUCCESS;
 }
 
+bool DDPandoraPFANewAlgorithm::validateThetaEnergyCorrectionSettings() const {
+  if (m_settings.m_electromagneticThetaEnergyCorrectionEnabled &&
+      !isConsistentThetaEnergyTable(m_settings.m_electromagneticThetaEnergyCorrectionThetaBinEdges,
+                                    m_settings.m_electromagneticThetaEnergyCorrectionEnergyBinEdges,
+                                    m_settings.m_electromagneticThetaEnergyCorrectionScaleFactors)) {
+    error() << "ElectromagneticThetaEnergyCorrection: bin edges must be strictly increasing and "
+            << "ScaleFactors must hold (nTheta-1)*(nEnergy-1) entries; got "
+            << m_settings.m_electromagneticThetaEnergyCorrectionThetaBinEdges.size() << " theta edges, "
+            << m_settings.m_electromagneticThetaEnergyCorrectionEnergyBinEdges.size() << " energy edges, "
+            << m_settings.m_electromagneticThetaEnergyCorrectionScaleFactors.size() << " scale factors" << endmsg;
+    return false;
+  }
+
+  if (m_settings.m_hadronicThetaEnergyCorrectionEnabled &&
+      !isConsistentThetaEnergyTable(m_settings.m_hadronicThetaEnergyCorrectionThetaBinEdges,
+                                    m_settings.m_hadronicThetaEnergyCorrectionEnergyBinEdges,
+                                    m_settings.m_hadronicThetaEnergyCorrectionScaleFactors)) {
+    error() << "HadronicThetaEnergyCorrection: bin edges must be strictly increasing and "
+            << "ScaleFactors must hold (nTheta-1)*(nEnergy-1) entries; got "
+            << m_settings.m_hadronicThetaEnergyCorrectionThetaBinEdges.size() << " theta edges, "
+            << m_settings.m_hadronicThetaEnergyCorrectionEnergyBinEdges.size() << " energy edges, "
+            << m_settings.m_hadronicThetaEnergyCorrectionScaleFactors.size() << " scale factors" << endmsg;
+    return false;
+  }
+
+  return true;
+}
+
 void DDPandoraPFANewAlgorithm::finaliseSteeringParameters() {
   // ATTN: This function seems to be necessary for operations that cannot easily be performed at construction of the
   // processor, when the steering file is parsed e.g. the call to GEAR to get the inner bfield
@@ -448,6 +538,16 @@ void DDPandoraPFANewAlgorithm::finaliseSteeringParameters() {
   m_settings.m_outputEnergyCorrectionPoints = m_outputEnergyCorrectionPoints;
   m_settings.m_ecalInputEnergyCorrectionPoints = m_ecalInputEnergyCorrectionPoints;
   m_settings.m_ecalOutputEnergyCorrectionPoints = m_ecalOutputEnergyCorrectionPoints;
+  m_settings.m_electromagneticThetaEnergyCorrectionEnabled = m_emThetaEnergyEnabled.value();
+  m_settings.m_electromagneticThetaEnergyCorrectionPluginName = m_emThetaEnergyPluginName.value();
+  m_settings.m_electromagneticThetaEnergyCorrectionThetaBinEdges = m_emThetaEnergyThetaBinEdges.value();
+  m_settings.m_electromagneticThetaEnergyCorrectionEnergyBinEdges = m_emThetaEnergyEnergyBinEdges.value();
+  m_settings.m_electromagneticThetaEnergyCorrectionScaleFactors = m_emThetaEnergyScaleFactors.value();
+  m_settings.m_hadronicThetaEnergyCorrectionEnabled = m_hadThetaEnergyEnabled.value();
+  m_settings.m_hadronicThetaEnergyCorrectionPluginName = m_hadThetaEnergyPluginName.value();
+  m_settings.m_hadronicThetaEnergyCorrectionThetaBinEdges = m_hadThetaEnergyThetaBinEdges.value();
+  m_settings.m_hadronicThetaEnergyCorrectionEnergyBinEdges = m_hadThetaEnergyEnergyBinEdges.value();
+  m_settings.m_hadronicThetaEnergyCorrectionScaleFactors = m_hadThetaEnergyScaleFactors.value();
   m_settings.m_trackCreatorName = m_trackCreatorName;
   m_settings.m_detectorName = m_detectorName;
   m_caloHitCreatorSettings.m_eCalBarrelNormalVector = m_eCalBarrelNormalVector;