clus, iface, root: extend SCECorrection to transverse (y, z) displacements

abhatfnal · HaiwangYu · commit b3df3d3d30be · 2026-06-16T09:47:37.000-05:00
Adds the transverse components to the SCECorrection chain introduced
in this PR.  ISCEField gains optional displacement_y / displacement_z
[default no-op, so maps without transverse components are unaffected].
SCEFieldTH3 soft-loads the TrueBkwd_Displacement_{Y,Z}_{E,W}
histograms and routes all three components through a single shared
interp() helper.  SCECorrection now outputs the fully-corrected triple.

Behavior change [intentional]: SCECorrection output_scope goes from
{x_sce, y, z} to {x_sce, y_sce, z_sce}; stored_array_names likewise.
Downstream configs that cluster on the SCE scope therefore cluster in
fully-corrected 3D.

Validated end-to-end on 50 SBND crossing-muon DETSIM events: per-point
Dx residual vs the TH3 map is unchanged from the prior baseline
[1.83 / 2.23 um E/W], and the new Dy / Dz residuals close at the same
interpolation precision -- Dy 2.63 / 1.65 um, Dz 3.35 / 2.99 um;
max &lt;= 9.91 um across all three components.  Evidence in the
companion wcp-porting-validation PR (sbnd-sce-036-stage3a-yz branch).
diff --git a/clus/inc/WireCellClus/SCECorrection.h b/clus/inc/WireCellClus/SCECorrection.h
@@ -6,12 +6,15 @@
 // (typically WireCell::Root::SCEFieldTH3) lives in the root/ subpackage and
 // is looked up by jsonnet TypeName.
 //
-// Applies:
+// Applies (all in WCT mm):
 //   (1) T0:  x_t0  = x_raw - dirx * cluster_t0 * drift_speed
-//   (2) SCE: x_sce = x_t0  + field->displacement_x(apa, x_t0, y, z)
+//   (2) SCE: x_sce = x_t0 + field->displacement_x(apa, x_t0, y, z)
+//            y_sce = y    + field->displacement_y(apa, x_t0, y, z)
+//            z_sce = z    + field->displacement_z(apa, x_t0, y, z)
 //
 // East/West (apa==0/1) routing is handled inside the ISCEField implementation.
-// If no field is provided (nullptr), the SCE step is a no-op (T0 still applied).
+// If no field is provided (nullptr), the SCE step is a no-op (T0 still applied
+// to x; y, z pass through).
 //
 // Author: Avinay Bhat (UChicago), for SBND
 // Modeled on T0Correction by Haiwang Yu.
@@ -53,24 +56,30 @@ namespace WireCell::Clus {
                                double cluster_t0, int face, int apa) const override;
 
         virtual PointCloud::Tree::Scope output_scope() const override {
-            return {"3d", {"x_sce", "y", "z"}};
+            return {"3d", {"x_sce", "y_sce", "z_sce"}};
         }
         virtual std::vector<std::string> stored_array_names() const override {
-            return {"x_sce"};
+            return {"x_sce", "y_sce", "z_sce"};
         }
 
     private:
         IDetectorVolumes::pointer m_dv;
         WireCell::ISCEField::pointer m_field;
         std::map<int, std::map<int, double>> m_drift_speeds;
 
-        // Combined T0 + SCE for one X coordinate.  All quantities in WCT mm.
-        inline double corrected_x(double x_raw_mm, double y_mm, double z_mm,
-                                  double cluster_t0, int face, int apa) const {
+        // Combined T0 + full 3D SCE for one point.  All quantities in WCT mm.
+        inline void corrected_xyz(double x_raw_mm, double y_mm, double z_mm,
+                                  double cluster_t0, int face, int apa,
+                                  double& xs, double& ys, double& zs) const {
             const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
             const double x_t0_mm = x_raw_mm - dirx * cluster_t0 * m_drift_speeds.at(apa).at(face);
-            const double dx_mm = m_field ? m_field->displacement_x(apa, x_t0_mm, y_mm, z_mm) : 0.0;
-            return x_t0_mm + dx_mm;
+            if (m_field) {
+                xs = x_t0_mm + m_field->displacement_x(apa, x_t0_mm, y_mm, z_mm);
+                ys = y_mm    + m_field->displacement_y(apa, x_t0_mm, y_mm, z_mm);
+                zs = z_mm    + m_field->displacement_z(apa, x_t0_mm, y_mm, z_mm);
+            } else {
+                xs = x_t0_mm; ys = y_mm; zs = z_mm;
+            }
         }
     };
 
@@ -83,21 +92,23 @@ namespace WireCell::Clus {
             const auto md = m_dv->metadata(wpid);
             m_drift_speeds[wpid.apa()][wpid.face()] = md["drift_speed"].asDouble();
         }
-        if (m_field) spdlog::info("SCECorrection: ISCEField wired in");
+        if (m_field) spdlog::info("SCECorrection: ISCEField wired in (x,y,z)");
         else         spdlog::info("SCECorrection: no ISCEField -- SCE disabled (T0 still applied)");
     }
 
     inline Point SCECorrection::forward(const Point& pos_raw, double t0, int face, int apa) const {
         Point pc(pos_raw);
-        pc[0] = corrected_x(pos_raw[0], pos_raw[1], pos_raw[2], t0, face, apa);
+        corrected_xyz(pos_raw[0], pos_raw[1], pos_raw[2], t0, face, apa, pc[0], pc[1], pc[2]);
         return pc;
     }
 
     inline Point SCECorrection::backward(const Point& pos_cor, double t0, int face, int apa) const {
         Point pr(pos_cor);
         if (m_field) {
-            // Approximate: evaluate field at x_sce instead of x_t0 (small displacement OK).
+            // Approximate: evaluate field at the corrected position.
             pr[0] -= m_field->displacement_x(apa, pos_cor[0], pos_cor[1], pos_cor[2]);
+            pr[1] -= m_field->displacement_y(apa, pos_cor[0], pos_cor[1], pos_cor[2]);
+            pr[2] -= m_field->displacement_z(apa, pos_cor[0], pos_cor[1], pos_cor[2]);
         }
         const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
         pr[0] += dirx * t0 * m_drift_speeds.at(apa).at(face);
@@ -117,14 +128,15 @@ namespace WireCell::Clus {
         const auto& ax = pc_raw.get(arr_raw_names[0])->elements<double>();
         const auto& ay = pc_raw.get(arr_raw_names[1])->elements<double>();
         const auto& az = pc_raw.get(arr_raw_names[2])->elements<double>();
-        std::vector<double> ax_c(ax.size());
+        std::vector<double> ax_c(ax.size()), ay_c(ax.size()), az_c(ax.size());
         for (size_t i = 0; i < ax.size(); ++i) {
-            ax_c[i] = corrected_x(ax[i], ay[i], az[i], t0, face, apa);
+            corrected_xyz(ax[i], ay[i], az[i], t0, face, apa,
+                          ax_c[i], ay_c[i], az_c[i]);
         }
         Dataset ds;
-        ds.add(arr_cor_names[0], Array(ax_c));
-        ds.add(arr_cor_names[1], Array(ay));
-        ds.add(arr_cor_names[2], Array(az));
+        ds.add(arr_cor_names[0], Array(ax_c));   // x_sce
+        ds.add(arr_cor_names[1], Array(ay_c));   // y_sce
+        ds.add(arr_cor_names[2], Array(az_c));   // z_sce
         return ds;
     }
 
@@ -138,17 +150,21 @@ namespace WireCell::Clus {
         const auto& az = pc_cor.get(arr_cor_names[2])->elements<double>();
         const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
         const double v = m_drift_speeds.at(apa).at(face);
-        std::vector<double> ax_r(ax.size());
+        std::vector<double> ax_r(ax.size()), ay_r(ax.size()), az_r(ax.size());
         for (size_t i = 0; i < ax.size(); ++i) {
-            double x = ax[i];
-            if (m_field) x -= m_field->displacement_x(apa, ax[i], ay[i], az[i]);
+            double x = ax[i], y = ay[i], z = az[i];
+            if (m_field) {
+                x -= m_field->displacement_x(apa, ax[i], ay[i], az[i]);
+                y -= m_field->displacement_y(apa, ax[i], ay[i], az[i]);
+                z -= m_field->displacement_z(apa, ax[i], ay[i], az[i]);
+            }
             x += dirx * t0 * v;
-            ax_r[i] = x;
+            ax_r[i] = x; ay_r[i] = y; az_r[i] = z;
         }
         Dataset ds;
         ds.add(arr_raw_names[0], Array(ax_r));
-        ds.add(arr_raw_names[1], Array(ay));
-        ds.add(arr_raw_names[2], Array(az));
+        ds.add(arr_raw_names[1], Array(ay_r));
+        ds.add(arr_raw_names[2], Array(az_r));
         return ds;
     }
 
diff --git a/iface/inc/WireCellIface/ISCEField.h b/iface/inc/WireCellIface/ISCEField.h
@@ -8,16 +8,17 @@ namespace WireCell {
     /** Per-TPC backward (reco -> true) Space-Charge-Effect displacement field.
      *
      *  Implementations expose, for a given APA/TPC index and 3D position,
-     *  the SIGNED backward X displacement to be added directly to a t0-corrected
-     *  X coordinate to recover the true X.  Out-of-domain queries are clamped
-     *  internally by the implementation.  Implementations that carry no map
-     *  should return 0 (no-op safe).
+     *  the SIGNED backward displacement to be added directly to a t0-corrected
+     *  coordinate to recover the true coordinate.  Out-of-domain queries are
+     *  clamped internally by the implementation.  Implementations that carry no
+     *  map should return 0 (no-op safe).
      *
      *  Units: WCT internal (mm) in, (mm) out.  Implementations handle any
      *  conversion required by their underlying map.
      *
-     *  Future revisions may add displacement_y / displacement_z if the
-     *  underlying map provides those components.
+     *  displacement_x is mandatory.  displacement_y / displacement_z are
+     *  optional: implementations whose underlying map lacks transverse
+     *  components inherit the default no-op (return 0).
      */
     class ISCEField : public IComponent<ISCEField> {
        public:
@@ -26,6 +27,12 @@ namespace WireCell {
         /// Backward X-displacement at (x,y,z) for APA index `apa`.
         /// SBND convention: apa=0 East, apa=1 West.  Units: mm.
         virtual double displacement_x(int apa, double x, double y, double z) const = 0;
+
+        /// Backward Y-displacement at (x,y,z).  Units: mm.  Optional (default no-op).
+        virtual double displacement_y(int /*apa*/, double /*x*/, double /*y*/, double /*z*/) const { return 0.0; }
+
+        /// Backward Z-displacement at (x,y,z).  Units: mm.  Optional (default no-op).
+        virtual double displacement_z(int /*apa*/, double /*x*/, double /*y*/, double /*z*/) const { return 0.0; }
     };
 
 }  // namespace WireCell
diff --git a/root/inc/WireCellRoot/SCEFieldTH3.h b/root/inc/WireCellRoot/SCEFieldTH3.h
@@ -1,8 +1,8 @@
 /** ISCEField implementation backed by per-TPC TH3F histograms in a
- *  ROOT file (SBND dualmap layout: TrueBkwd_Displacement_X_E,
- *  TrueBkwd_Displacement_X_W).  IConfigurable so jsonnet wiring works
- *  the usual way; the actual SCE-correction transform in clus/
- *  receives an ISCEField::pointer and never sees ROOT.
+ *  ROOT file (SBND dualmap layout: TrueBkwd_Displacement_{X,Y,Z}_{E,W}).
+ *  IConfigurable so jsonnet wiring works the usual way; the actual
+ *  SCE-correction transform in clus/ receives an ISCEField::pointer and
+ *  never sees ROOT.
  */
 
 #ifndef WIRECELLROOT_SCEFIELDTH3
@@ -30,6 +30,8 @@ namespace WireCell {
 
             // ISCEField
             virtual double displacement_x(int apa, double x, double y, double z) const override;
+            virtual double displacement_y(int apa, double x, double y, double z) const override;
+            virtual double displacement_z(int apa, double x, double y, double z) const override;
 
             // IConfigurable
             virtual void configure(const WireCell::Configuration& cfg) override;
@@ -38,19 +40,37 @@ namespace WireCell {
            private:
             // configuration
             std::string m_sce_map_file;
-            std::string m_th3_name_E;
+            std::string m_th3_name_E;     // X (mandatory)
             std::string m_th3_name_W;
+            std::string m_th3_name_E_y;   // Y (optional transverse)
+            std::string m_th3_name_W_y;
+            std::string m_th3_name_E_z;   // Z (optional transverse)
+            std::string m_th3_name_W_z;
             bool   m_axis_unit_mm{false};
             double m_sign{-1.0};
 
-            // loaded state
+            // loaded state (all owned by m_file)
             std::unique_ptr<TFile> m_file;
-            TH3F* m_hE{nullptr};  // owned by m_file
+            TH3F* m_hE{nullptr};
             TH3F* m_hW{nullptr};
+            TH3F* m_hE_y{nullptr};
+            TH3F* m_hW_y{nullptr};
+            TH3F* m_hE_z{nullptr};
+            TH3F* m_hW_z{nullptr};
 
             Log::logptr_t l;
 
+            // Shared: pick E/W histo by apa, clamp, convert units, apply sign.
+            // Returns 0 if the requested histo is absent (no-op component).
+            double interp(TH3F* hE, TH3F* hW, int apa,
+                          double x, double y, double z) const;
+
             static double clamp_axis(const TAxis* a, double v);
+
+            // Load one (E,W) histo pair; hard=true throws on missing,
+            // hard=false warns and leaves pointers null.
+            void load_pair(const std::string& nameE, const std::string& nameW,
+                           TH3F*& hE, TH3F*& hW, bool hard);
         };
 
     }  // namespace Root
diff --git a/root/src/SCEFieldTH3.cxx b/root/src/SCEFieldTH3.cxx
@@ -16,6 +16,10 @@ using namespace WireCell::Root;
 SCEFieldTH3::SCEFieldTH3()
     : m_th3_name_E("TrueBkwd_Displacement_X_E"),
       m_th3_name_W("TrueBkwd_Displacement_X_W"),
+      m_th3_name_E_y("TrueBkwd_Displacement_Y_E"),
+      m_th3_name_W_y("TrueBkwd_Displacement_Y_W"),
+      m_th3_name_E_z("TrueBkwd_Displacement_Z_E"),
+      m_th3_name_W_z("TrueBkwd_Displacement_Z_W"),
       l(Log::logger("wct"))
 {
 }
@@ -25,19 +29,46 @@ SCEFieldTH3::~SCEFieldTH3() = default;
 Configuration SCEFieldTH3::default_configuration() const
 {
     Configuration cfg;
-    cfg["sce_map_file"] = m_sce_map_file;
-    cfg["th3_name_E"]   = m_th3_name_E;
-    cfg["th3_name_W"]   = m_th3_name_W;
-    cfg["axis_unit_mm"] = m_axis_unit_mm;
-    cfg["sign"]         = m_sign;
+    cfg["sce_map_file"]  = m_sce_map_file;
+    cfg["th3_name_E"]    = m_th3_name_E;
+    cfg["th3_name_W"]    = m_th3_name_W;
+    cfg["th3_name_E_y"]  = m_th3_name_E_y;
+    cfg["th3_name_W_y"]  = m_th3_name_W_y;
+    cfg["th3_name_E_z"]  = m_th3_name_E_z;
+    cfg["th3_name_W_z"]  = m_th3_name_W_z;
+    cfg["axis_unit_mm"]  = m_axis_unit_mm;
+    cfg["sign"]          = m_sign;
     return cfg;
 }
 
+void SCEFieldTH3::load_pair(const std::string& nameE, const std::string& nameW,
+                            TH3F*& hE, TH3F*& hW, bool hard)
+{
+    hE = dynamic_cast<TH3F*>(m_file->Get(nameE.c_str()));
+    hW = dynamic_cast<TH3F*>(m_file->Get(nameW.c_str()));
+    if (!hE || !hW) {
+        if (hard) {
+            throw std::runtime_error("SCEFieldTH3: missing TH3F '" + nameE
+                                     + "' or '" + nameW + "' in " + m_sce_map_file);
+        }
+        l->warn("SCEFieldTH3: optional TH3 '{}' / '{}' absent -- that component is a no-op",
+                nameE, nameW);
+        hE = nullptr; hW = nullptr;
+        return;
+    }
+    hE->SetDirectory(nullptr);
+    hW->SetDirectory(nullptr);
+}
+
 void SCEFieldTH3::configure(const Configuration& cfg)
 {
     m_sce_map_file = get<std::string>(cfg, "sce_map_file", m_sce_map_file);
     m_th3_name_E   = get<std::string>(cfg, "th3_name_E",   m_th3_name_E);
     m_th3_name_W   = get<std::string>(cfg, "th3_name_W",   m_th3_name_W);
+    m_th3_name_E_y = get<std::string>(cfg, "th3_name_E_y", m_th3_name_E_y);
+    m_th3_name_W_y = get<std::string>(cfg, "th3_name_W_y", m_th3_name_W_y);
+    m_th3_name_E_z = get<std::string>(cfg, "th3_name_E_z", m_th3_name_E_z);
+    m_th3_name_W_z = get<std::string>(cfg, "th3_name_W_z", m_th3_name_W_z);
     m_axis_unit_mm = get<bool>       (cfg, "axis_unit_mm", m_axis_unit_mm);
     m_sign         = get<double>     (cfg, "sign",         m_sign);
 
@@ -51,16 +82,15 @@ void SCEFieldTH3::configure(const Configuration& cfg)
     if (!m_file || m_file->IsZombie()) {
         throw std::runtime_error("SCEFieldTH3: cannot open " + m_sce_map_file);
     }
-    m_hE = dynamic_cast<TH3F*>(m_file->Get(m_th3_name_E.c_str()));
-    m_hW = dynamic_cast<TH3F*>(m_file->Get(m_th3_name_W.c_str()));
-    if (!m_hE || !m_hW) {
-        throw std::runtime_error("SCEFieldTH3: missing TH3F '" + m_th3_name_E
-                                 + "' or '" + m_th3_name_W + "' in " + m_sce_map_file);
-    }
-    m_hE->SetDirectory(nullptr);
-    m_hW->SetDirectory(nullptr);
-    l->info("SCEFieldTH3: loaded TH3 E='{}' W='{}' axis_mm={} sign={}",
-            m_th3_name_E, m_th3_name_W, m_axis_unit_mm, m_sign);
+
+    load_pair(m_th3_name_E,   m_th3_name_W,   m_hE,   m_hW,   /*hard=*/true);
+    load_pair(m_th3_name_E_y, m_th3_name_W_y, m_hE_y, m_hW_y, /*hard=*/false);
+    load_pair(m_th3_name_E_z, m_th3_name_W_z, m_hE_z, m_hW_z, /*hard=*/false);
+
+    l->info("SCEFieldTH3: loaded X[E='{}' W='{}'] Y[{}] Z[{}] axis_mm={} sign={}",
+            m_th3_name_E, m_th3_name_W,
+            (m_hE_y ? "ok" : "absent"), (m_hE_z ? "ok" : "absent"),
+            m_axis_unit_mm, m_sign);
 }
 
 double SCEFieldTH3::clamp_axis(const TAxis* a, double v)
@@ -73,9 +103,10 @@ double SCEFieldTH3::clamp_axis(const TAxis* a, double v)
     return v;
 }
 
-double SCEFieldTH3::displacement_x(int apa, double x, double y, double z) const
+double SCEFieldTH3::interp(TH3F* hE, TH3F* hW, int apa,
+                           double x, double y, double z) const
 {
-    TH3F* h = (apa == 0) ? m_hE : (apa == 1) ? m_hW : nullptr;
+    TH3F* h = (apa == 0) ? hE : (apa == 1) ? hW : nullptr;
     if (!h) return 0.0;
 
     // WCT input is mm; map axes are cm by default (mm if m_axis_unit_mm).
@@ -85,7 +116,22 @@ double SCEFieldTH3::displacement_x(int apa, double x, double y, double z) const
     const double xq = clamp_axis(h->GetXaxis(), x * s_in);
     const double yq = clamp_axis(h->GetYaxis(), y * s_in);
     const double zq = clamp_axis(h->GetZaxis(), z * s_in);
-    const double dx_map = h->Interpolate(xq, yq, zq);
+    const double d_map = h->Interpolate(xq, yq, zq);
 
-    return m_sign * dx_map * s_out;  // signed displacement in mm
+    return m_sign * d_map * s_out;  // signed displacement in mm
+}
+
+double SCEFieldTH3::displacement_x(int apa, double x, double y, double z) const
+{
+    return interp(m_hE, m_hW, apa, x, y, z);
+}
+
+double SCEFieldTH3::displacement_y(int apa, double x, double y, double z) const
+{
+    return interp(m_hE_y, m_hW_y, apa, x, y, z);
+}
+
+double SCEFieldTH3::displacement_z(int apa, double x, double y, double z) const
+{
+    return interp(m_hE_z, m_hW_z, apa, x, y, z);
 }
