Implement H4, H3 and S2 staggering in HEXPLIT

src/algorithms/calorimetry/CalorimeterHitReco.cc

@@ -270,7 +270,8 @@ void CalorimeterHitReco::process(const CalorimeterHitReco::Input& input,
       segmentation_type = segmentation->type();
     }
 
-    if (segmentation_type == "CartesianGridXY" || segmentation_type == "HexGridXY") {
+    if (segmentation_type == "CartesianGridXY" || segmentation_type == "HexGridXY" ||
+        segmentation_type == "CartesianGridXYStaggered") {
       auto cell_dim = m_converter->cellDimensions(cellID);
       cdim.resize(3);
       cdim[0] = cell_dim[0];

src/algorithms/calorimetry/HEXPLIT.cc

@@ -15,8 +15,9 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdlib>
-#include <gsl/pointers> // for not_null
+#include <iostream>
 #include <numbers>
+#include <tuple>
 #include <vector>
 
 #include "HEXPLIT.h"
@@ -25,7 +26,7 @@
 namespace eicrecon {
 
 //positions where the overlapping cells are relative to a given cell (in units of hexagon side length)
-const std::vector<double> HEXPLIT::neighbor_offsets_x = []() {
+const std::vector<double> HEXPLIT::neighbor_offsets_x_H4 = []() {
   std::vector<double> x;
   double rs[2]      = {1.5, std::numbers::sqrt3 / 2.};
   double offsets[2] = {0, M_PI / 2};
@@ -37,7 +38,7 @@ const std::vector<double> HEXPLIT::neighbor_offsets_x = []() {
   return x;
 }();
 
-const std::vector<double> HEXPLIT::neighbor_offsets_y = []() {
+const std::vector<double> HEXPLIT::neighbor_offsets_y_H4 = []() {
   std::vector<double> y;
   double rs[2]      = {1.5, std::numbers::sqrt3 / 2.};
   double offsets[2] = {0, M_PI / 2};
@@ -50,12 +51,12 @@ const std::vector<double> HEXPLIT::neighbor_offsets_y = []() {
 }();
 
 //indices of the neighboring cells which overlap to produce a given subcell
-const int HEXPLIT::neighbor_indices[SUBCELLS][OVERLAP] = {
+const std::vector<std::vector<int>> HEXPLIT::neighbor_indices_H4 = {
     {0, 11, 10}, {1, 6, 11}, {2, 7, 6}, {3, 8, 7},  {4, 9, 8},   {5, 10, 9},
     {6, 11, 7},  {7, 6, 8},  {8, 7, 9}, {9, 8, 10}, {10, 9, 11}, {11, 10, 6}};
 
 //positions of the centers of subcells
-const std::vector<double> HEXPLIT::subcell_offsets_x = []() {
+const std::vector<double> HEXPLIT::subcell_offsets_x_H4 = []() {
   std::vector<double> x;
   double rs[2]      = {0.75, std::numbers::sqrt3 / 4.};
   double offsets[2] = {0, M_PI / 2};
@@ -67,7 +68,7 @@ const std::vector<double> HEXPLIT::subcell_offsets_x = []() {
   return x;
 }();
 
-const std::vector<double> HEXPLIT::subcell_offsets_y = []() {
+const std::vector<double> HEXPLIT::subcell_offsets_y_H4 = []() {
   std::vector<double> y;
   double rs[2]      = {0.75, std::numbers::sqrt3 / 4.};
   double offsets[2] = {0, M_PI / 2};
@@ -79,7 +80,23 @@ const std::vector<double> HEXPLIT::subcell_offsets_y = []() {
   return y;
 }();
 
-void HEXPLIT::init() {}
+const std::vector<double> HEXPLIT::neighbor_offsets_x_S2 = {1, -1, -1, 1};
+const std::vector<double> HEXPLIT::neighbor_offsets_y_S2 = {1, 1, -1, -1};
+
+const std::vector<std::vector<int>> HEXPLIT::neighbor_indices_S2 = {{0}, {1}, {2}, {3}};
+
+const std::vector<double> HEXPLIT::subcell_offsets_x_S2 = {0.5, -0.5, -0.5, 0.5};
+const std::vector<double> HEXPLIT::subcell_offsets_y_S2 = {0.5, 0.5, -0.5, -0.5};
+
+void HEXPLIT::init() {
+  if (m_cfg.stag_type == HEXPLITConfig::StaggerType::H4) {
+    stag = stag_H4;
+  } else if (m_cfg.stag_type == HEXPLITConfig::StaggerType::S2) {
+    stag = stag_S2;
+  } else if (m_cfg.stag_type == HEXPLITConfig::StaggerType::H3) {
+    std::cout << "error   H3 staggering not implemented yet in EICrecon" << std::endl;
+  }
+}
 
 void HEXPLIT::process(const HEXPLIT::Input& input, const HEXPLIT::Output& output) const {
 
@@ -100,7 +117,7 @@ void HEXPLIT::process(const HEXPLIT::Input& input, const HEXPLIT::Output& output
     }
 
     //keep track of the energy in each neighboring cell
-    std::vector<double> Eneighbors(NEIGHBORS, 0.0);
+    std::vector<double> Eneighbors(stag.NEIGHBORS, 0.0);
 
     double sl = hit.getDimension().x / 2.;
     for (const auto& other_hit : *hits) {
@@ -126,30 +143,45 @@ void HEXPLIT::process(const HEXPLIT::Input& input, const HEXPLIT::Output& output
 
       //loop over locations of the neighboring cells
       //and check if the jth hit matches this location
-      for (int k = 0; k < NEIGHBORS; k++) {
-        if (std::abs(dx - neighbor_offsets_x[k]) < tol &&
-            std::abs(dy - neighbor_offsets_y[k]) < tol) {
+      for (int k = 0; k < stag.NEIGHBORS; k++) {
+        if (std::abs(dx - stag.neighbor_offsets_x[k]) < tol &&
+            std::abs(dy - stag.neighbor_offsets_y[k]) < tol) {
           Eneighbors[k] += other_hit.getEnergy();
           break;
         }
       }
     }
-    double weights[SUBCELLS];
-    for (int k = 0; k < NEIGHBORS; k++) {
+
+    double weights[12];
+    for (int k = 0; k < stag.NEIGHBORS; k++) {
       Eneighbors[k] = std::max(Eneighbors[k], delta);
     }
     double sum_weights = 0;
-    for (int k = 0; k < SUBCELLS; k++) {
-      weights[k] = Eneighbors[neighbor_indices[k][0]] * Eneighbors[neighbor_indices[k][1]] *
-                   Eneighbors[neighbor_indices[k][2]];
-      sum_weights += weights[k];
+    if (m_cfg.stag_type == HEXPLITConfig::StaggerType::H4)
+      for (int k = 0; k < stag.SUBCELLS; k++) {
+        weights[k] = Eneighbors[stag.neighbor_indices[k][0]] *
+                     Eneighbors[stag.neighbor_indices[k][1]] *
+                     Eneighbors[stag.neighbor_indices[k][2]];
+        sum_weights += weights[k];
+      }
+    else if (m_cfg.stag_type == HEXPLITConfig::StaggerType::S2) {
+      for (int k = 0; k < stag.SUBCELLS; k++) {
+        weights[k] = Eneighbors[stag.neighbor_indices[k][0]];
+        sum_weights += weights[k];
+      }
+    } else if (m_cfg.stag_type == HEXPLITConfig::StaggerType::H3) {
+      for (int k = 0; k < stag.SUBCELLS; k++) {
+        weights[k] =
+            Eneighbors[stag.neighbor_indices[k][0]] * Eneighbors[stag.neighbor_indices[k][1]];
+        sum_weights += weights[k];
+      }
     }
-    for (int k = 0; k < SUBCELLS; k++) {
+    for (int k = 0; k < stag.SUBCELLS; k++) {
 
       //create the subcell hits.  First determine their positions in local coordinates.
       const decltype(edm4eic::CalorimeterHitData::local) local(
-          hit.getLocal().x + subcell_offsets_x[k] * sl,
-          hit.getLocal().y + subcell_offsets_y[k] * sl, hit.getLocal().z);
+          hit.getLocal().x + stag.subcell_offsets_x[k] * sl,
+          hit.getLocal().y + stag.subcell_offsets_y[k] * sl, hit.getLocal().z);
 
       //convert this to a position object so that the global position can be determined
       dd4hep::Position local_position;

src/algorithms/calorimetry/HEXPLIT.h

@@ -38,23 +38,72 @@ class HEXPLIT : public HEXPLITAlgorithm, public WithPodConfig<HEXPLITConfig> {
   void process(const Input&, const Output&) const final;
 
 private:
+  typedef struct stagger_pattern {
+    int SUBCELLS;
+    int NEIGHBORS;
+    int OVERLAP;
+    std::vector<double> neighbor_offsets_x;
+    std::vector<double> neighbor_offsets_y;
+    std::vector<std::vector<int>> neighbor_indices;
+    std::vector<double> subcell_offsets_x;
+    std::vector<double> subcell_offsets_y;
+  } stagger_pattern;
+
   // number of subcells that a single cell is divided into
-  static const int SUBCELLS = 12;
+  static const int SUBCELLS_H4 = 12;
+  // number of neighboring positions whose overlap define the subcells
+  static const int NEIGHBORS_H4 = 12;
+  // number of neighboring cells that overlap to obtain a subcell
+  static const int OVERLAP_H4 = 3;
+  //positions where the overlapping cells are relative to a given cell (in units of hexagon side length)
+  static const std::vector<double> neighbor_offsets_x_H4;
+  static const std::vector<double> neighbor_offsets_y_H4;
+  //indices of the neighboring cells which overlap to produce a given subcell
+  static const std::vector<std::vector<int>> neighbor_indices_H4;
+  //positions of the centers of subcells
+  static const std::vector<double> subcell_offsets_x_H4;
+  static const std::vector<double> subcell_offsets_y_H4;
+
+  const stagger_pattern stag_H4 = {
+      .SUBCELLS           = SUBCELLS_H4,
+      .NEIGHBORS          = NEIGHBORS_H4,
+      .OVERLAP            = OVERLAP_H4,
+      .neighbor_offsets_x = neighbor_offsets_x_H4,
+      .neighbor_offsets_y = neighbor_offsets_y_H4,
+      .neighbor_indices   = neighbor_indices_H4,
+      .subcell_offsets_x  = subcell_offsets_x_H4,
+      .subcell_offsets_y  = subcell_offsets_y_H4,
+  };
+
+  static const int SUBCELLS_S2 = 4;
   // number of neighboring positions whose overlap define the subcells
-  static const int NEIGHBORS = 12;
+  static const int NEIGHBORS_S2 = 4;
   // number of neighboring cells that overlap to obtain a subcell
-  static const int OVERLAP = 3;
+  static const int OVERLAP_S2 = 1;
   //positions where the overlapping cells are relative to a given cell (in units of hexagon side length)
-  static const std::vector<double> neighbor_offsets_x;
-  static const std::vector<double> neighbor_offsets_y;
+  static const std::vector<double> neighbor_offsets_x_S2;
+  static const std::vector<double> neighbor_offsets_y_S2;
   //indices of the neighboring cells which overlap to produce a given subcell
-  static const int neighbor_indices[SUBCELLS][OVERLAP];
+  static const std::vector<std::vector<int>> neighbor_indices_S2;
   //positions of the centers of subcells
-  static const std::vector<double> subcell_offsets_x;
-  static const std::vector<double> subcell_offsets_y;
+  static const std::vector<double> subcell_offsets_x_S2;
+  static const std::vector<double> subcell_offsets_y_S2;
+
+  const stagger_pattern stag_S2 = {
+      .SUBCELLS           = SUBCELLS_S2,
+      .NEIGHBORS          = NEIGHBORS_S2,
+      .OVERLAP            = OVERLAP_S2,
+      .neighbor_offsets_x = neighbor_offsets_x_S2,
+      .neighbor_offsets_y = neighbor_offsets_y_S2,
+      .neighbor_indices   = neighbor_indices_S2,
+      .subcell_offsets_x  = subcell_offsets_x_S2,
+      .subcell_offsets_y  = subcell_offsets_y_S2,
+  };
+
+  stagger_pattern stag = stag_H4;
 
 private:
   const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()};
-};
+}; // namespace eicrecon
 
 } // namespace eicrecon

src/algorithms/calorimetry/HEXPLITConfig.h

@@ -10,6 +10,7 @@ struct HEXPLITConfig {
   double Emin_in_MIPs{0.1};
   double delta_in_MIPs{0.01};
   double tmax{325 * dd4hep::ns};
+  enum StaggerType { H4 = 0, H3 = 1, S2 = 2 } stag_type = H4;
 };
 
 } // namespace eicrecon

src/detectors/ZDC/ZDC.cc

@@ -3,12 +3,13 @@
 
 #include <Evaluator/DD4hepUnits.h>
 #include <JANA/JApplicationFwd.h>
-#include <edm4eic/EDM4eicVersion.h>
 #include <JANA/Utils/JTypeInfo.h>
+#include <edm4eic/EDM4eicVersion.h>
 #include <string>
 #include <variant>
 #include <vector>
 
+#include "algorithms/calorimetry/HEXPLITConfig.h"
 #include "algorithms/calorimetry/ImagingTopoClusterConfig.h"
 #include "extensions/jana/JOmniFactoryGeneratorT.h"
 #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h"
@@ -194,36 +195,36 @@ void InitPlugin(JApplication* app) {
       app // TODO: Remove me once fixed
       ));
 
-  app->Add(new JOmniFactoryGeneratorT<HEXPLIT_factory>("HcalFarForwardZDCSubcellHits",
-                                                       {"HcalFarForwardZDCRecHits"},
-                                                       {"HcalFarForwardZDCSubcellHits"},
-                                                       {
-                                                           .MIP           = 472. * dd4hep::keV,
-                                                           .Emin_in_MIPs  = 0.5,
-                                                           .delta_in_MIPs = 0.01,
-                                                           .tmax          = 269 * dd4hep::ns,
-                                                       },
-                                                       app // TODO: Remove me once fixed
-                                                       ));
+  app->Add(new JOmniFactoryGeneratorT<HEXPLIT_factory>(
+      "HcalFarForwardZDCSubcellHits", {"HcalFarForwardZDCRecHits"},
+      {"HcalFarForwardZDCSubcellHits"},
+      {
+          .MIP           = 630. * dd4hep::keV,
+          .Emin_in_MIPs  = 0.5,
+          .delta_in_MIPs = 0.01,
+          .tmax          = 269 * dd4hep::ns,
+          .stag_type     = HEXPLITConfig::StaggerType::S2,
+      },
+      app // TODO: Remove me once fixed
+      ));
 
+  double side_length = 48.8 * dd4hep::mm;
   app->Add(new JOmniFactoryGeneratorT<ImagingTopoCluster_factory>(
       "HcalFarForwardZDCImagingProtoClusters", {"HcalFarForwardZDCSubcellHits"},
       {"HcalFarForwardZDCImagingProtoClusters"},
       {
           .neighbourLayersRange = 1,
-          .sameLayerDistXY   = {"0.5 * HcalFarForwardZDC_SiPMonTile_HexSideLength",
-                                "0.5 * HcalFarForwardZDC_SiPMonTile_HexSideLength * sin(pi / 3)"},
-          .diffLayerDistXY   = {"0.5 * HcalFarForwardZDC_SiPMonTile_HexSideLength",
-                                "0.5 * HcalFarForwardZDC_SiPMonTile_HexSideLength * sin(pi / 3)"},
-          .sameLayerMode     = eicrecon::ImagingTopoClusterConfig::ELayerMode::xy,
-          .sectorDist        = 10.0 * dd4hep::cm,
-          .minClusterHitEdep = 50.0 * dd4hep::keV,
-          .minClusterCenterEdep = 3.0 * dd4hep::MeV,
-          .minClusterEdep       = 11.0 * dd4hep::MeV,
-          .minClusterNhits      = 30,
+          .sameLayerDistXY      = {side_length * 0.75, side_length * 0.75},
+          .diffLayerDistXY      = {side_length * 0.75, side_length * 0.75},
+          .sameLayerMode        = eicrecon::ImagingTopoClusterConfig::ELayerMode::xy,
+          .diffLayerMode        = eicrecon::ImagingTopoClusterConfig::ELayerMode::xy,
+          .sectorDist           = 10.0 * dd4hep::cm,
+          .minClusterHitEdep    = 315.0 * dd4hep::keV,
+          .minClusterCenterEdep = 25 * dd4hep::MeV,
+          .minClusterEdep       = 50.0 * dd4hep::MeV,
+          .minClusterNhits      = 10,
       },
-      app // TODO: Remove me once fixed
-      ));
+      app));
 
   app->Add(new JOmniFactoryGeneratorT<CalorimeterIslandCluster_factory>(
       "HcalFarForwardZDCIslandProtoClusters", {"HcalFarForwardZDCSubcellHits"},
@@ -232,8 +233,8 @@ void InitPlugin(JApplication* app) {
        .peakNeighbourhoodMatrix{},
        .readout{},
        .sectorDist  = 1.5 * dd4hep::cm,
-       .localDistXY = {"0.9 * HcalFarForwardZDC_SiPMonTile_HexSideLength",
-                       "0.76 * HcalFarForwardZDC_SiPMonTile_HexSideLength * sin(pi / 3)"},
+       .localDistXY = {"0.55 * HcalFarForwardZDC_SiPMonTile_SquareSideLength",
+                       "0.55 * HcalFarForwardZDC_SiPMonTile_SquareSideLength"},
        .localDistXZ{},
        .localDistYZ{},
        .globalDistRPhi{},
@@ -261,8 +262,9 @@ void InitPlugin(JApplication* app) {
        "HcalFarForwardZDCClusterLinksWithoutShapes",
 #endif
        "HcalFarForwardZDCClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
-      {.energyWeight        = "log",
-       .sampFrac            = 0.0203,
+      {.energyWeight = "log",
+       .sampFrac     = dd4hep::_toDouble(
+           "0.0273 * HcalFarForwardZDC_SiPMonTile_PolystyreneThickness / ( 4 * mm )"),
        .logWeightBaseCoeffs = {5.8, 0.65, 0.31},
        .logWeightBase_Eref  = 50 * dd4hep::GeV},
       app // TODO: Remove me once fixed
@@ -279,9 +281,10 @@ void InitPlugin(JApplication* app) {
        "HcalFarForwardZDCClusterAssociations"},
       {.longitudinalShowerInfoAvailable = true,
        .energyWeight                    = "log",
-       .sampFrac                        = 0.0203,
-       .logWeightBaseCoeffs             = {5.8, 0.65, 0.31},
-       .logWeightBase_Eref              = 50 * dd4hep::GeV},
+       .sampFrac                        = dd4hep::_toDouble(
+           "0.0273 * HcalFarForwardZDC_SiPMonTile_PolystyreneThickness / ( 4 * mm )"),
+       .logWeightBaseCoeffs = {5.8, 0.65, 0.31},
+       .logWeightBase_Eref  = 50 * dd4hep::GeV},
       app));
 
   app->Add(new JOmniFactoryGeneratorT<CalorimeterTruthClustering_factory>(
@@ -357,8 +360,9 @@ void InitPlugin(JApplication* app) {
 #endif
        "HcalFarForwardZDCClusterAssociationsBaselineWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
       {
-          .energyWeight    = "log",
-          .sampFrac        = 0.0203,
+          .energyWeight = "log",
+          .sampFrac     = dd4hep::_toDouble(
+              "0.0273 * HcalFarForwardZDC_SiPMonTile_PolystyreneThickness / ( 4 * mm )"),
           .logWeightBase   = 6.2,
           .enableEtaBounds = false,
       },
@@ -376,8 +380,9 @@ void InitPlugin(JApplication* app) {
        "HcalFarForwardZDCClusterAssociationsBaseline"},
       {.longitudinalShowerInfoAvailable = true,
        .energyWeight                    = "log",
-       .sampFrac                        = 0.0203,
-       .logWeightBase                   = 6.2},
+       .sampFrac                        = dd4hep::_toDouble(
+           "0.0273 * HcalFarForwardZDC_SiPMonTile_PolystyreneThickness / ( 4 * mm )"),
+       .logWeightBase = 6.2},
       app));
 }
 }

src/global/reco/reco.cc

@@ -225,7 +225,7 @@ void InitPlugin(JApplication* app) {
        .globalToProtonRotation    = -0.025,
        .gammaZMaxOffset           = 300 * dd4hep::mm,
        .gammaMaxLength            = 100 * dd4hep::mm,
-       .gammaMaxWidth             = 12 * dd4hep::mm},
+       .gammaMaxWidth             = 27 * dd4hep::mm},
       app // TODO: Remove me once fixed
       ));