Skip to content

Commit 3d1ebf9

Browse files
HaiwangYuHaiwangYu
authored
Merge pull request #482 from HaiwangYu/sce
Cherry-pick PR479 to master instead of match
2 parents 527dfe3 + b3df3d3 commit 3d1ebf9

5 files changed

Lines changed: 466 additions & 0 deletions

File tree

‎clus/inc/WireCellClus/SCECorrection.h‎

Lines changed: 191 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
1+
// SCECorrection.h
2+
//
3+
// A Clus::IPCTransform that performs combined T0 + SCE (Space Charge Effect)
4+
// correction for SBND, delegating the SCE displacement field to an externally
5+
// provided ISCEField. This keeps clus/ ROOT-free; the field implementation
6+
// (typically WireCell::Root::SCEFieldTH3) lives in the root/ subpackage and
7+
// is looked up by jsonnet TypeName.
8+
//
9+
// Applies (all in WCT mm):
10+
// (1) T0: x_t0 = x_raw - dirx * cluster_t0 * drift_speed
11+
// (2) SCE: x_sce = x_t0 + field->displacement_x(apa, x_t0, y, z)
12+
// y_sce = y + field->displacement_y(apa, x_t0, y, z)
13+
// z_sce = z + field->displacement_z(apa, x_t0, y, z)
14+
//
15+
// East/West (apa==0/1) routing is handled inside the ISCEField implementation.
16+
// If no field is provided (nullptr), the SCE step is a no-op (T0 still applied
17+
// to x; y, z pass through).
18+
//
19+
// Author: Avinay Bhat (UChicago), for SBND
20+
// Modeled on T0Correction by Haiwang Yu.
21+
22+
#ifndef WIRECELLCLUS_SCECORRECTION_H
23+
#define WIRECELLCLUS_SCECORRECTION_H
24+
25+
#include "WireCellClus/IPCTransform.h"
26+
#include "WireCellIface/IDetectorVolumes.h"
27+
#include "WireCellIface/ISCEField.h"
28+
#include "WireCellUtil/Logging.h"
29+
30+
#include <map>
31+
#include <memory>
32+
#include <string>
33+
34+
namespace WireCell::Clus {
35+
36+
class SCECorrection : public WireCell::Clus::IPCTransform {
37+
public:
38+
SCECorrection(IDetectorVolumes::pointer dv,
39+
WireCell::ISCEField::pointer field = nullptr);
40+
virtual ~SCECorrection() = default;
41+
42+
virtual Point forward(const Point& pos_raw, double cluster_t0, int face, int apa) const override;
43+
virtual Point backward(const Point& pos_cor, double cluster_t0, int face, int apa) const override;
44+
virtual bool filter(const Point& pos_cor, double cluster_t0, int face, int apa) const override;
45+
46+
virtual Dataset forward(const Dataset& pc_raw,
47+
const std::vector<std::string>& arr_raw_names,
48+
const std::vector<std::string>& arr_cor_names,
49+
double cluster_t0, int face, int apa) const override;
50+
virtual Dataset backward(const Dataset& pc_cor,
51+
const std::vector<std::string>& arr_cor_names,
52+
const std::vector<std::string>& arr_raw_names,
53+
double cluster_t0, int face, int apa) const override;
54+
virtual Dataset filter(const Dataset& pc_cor,
55+
const std::vector<std::string>& arr_cor_names,
56+
double cluster_t0, int face, int apa) const override;
57+
58+
virtual PointCloud::Tree::Scope output_scope() const override {
59+
return {"3d", {"x_sce", "y_sce", "z_sce"}};
60+
}
61+
virtual std::vector<std::string> stored_array_names() const override {
62+
return {"x_sce", "y_sce", "z_sce"};
63+
}
64+
65+
private:
66+
IDetectorVolumes::pointer m_dv;
67+
WireCell::ISCEField::pointer m_field;
68+
std::map<int, std::map<int, double>> m_drift_speeds;
69+
70+
// Combined T0 + full 3D SCE for one point. All quantities in WCT mm.
71+
inline void corrected_xyz(double x_raw_mm, double y_mm, double z_mm,
72+
double cluster_t0, int face, int apa,
73+
double& xs, double& ys, double& zs) const {
74+
const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
75+
const double x_t0_mm = x_raw_mm - dirx * cluster_t0 * m_drift_speeds.at(apa).at(face);
76+
if (m_field) {
77+
xs = x_t0_mm + m_field->displacement_x(apa, x_t0_mm, y_mm, z_mm);
78+
ys = y_mm + m_field->displacement_y(apa, x_t0_mm, y_mm, z_mm);
79+
zs = z_mm + m_field->displacement_z(apa, x_t0_mm, y_mm, z_mm);
80+
} else {
81+
xs = x_t0_mm; ys = y_mm; zs = z_mm;
82+
}
83+
}
84+
};
85+
86+
inline SCECorrection::SCECorrection(IDetectorVolumes::pointer dv,
87+
WireCell::ISCEField::pointer field)
88+
: m_dv(dv), m_field(field)
89+
{
90+
for (const auto& [wfid, _] : m_dv->wpident_faces()) {
91+
WirePlaneId wpid(wfid);
92+
const auto md = m_dv->metadata(wpid);
93+
m_drift_speeds[wpid.apa()][wpid.face()] = md["drift_speed"].asDouble();
94+
}
95+
if (m_field) spdlog::info("SCECorrection: ISCEField wired in (x,y,z)");
96+
else spdlog::info("SCECorrection: no ISCEField -- SCE disabled (T0 still applied)");
97+
}
98+
99+
inline Point SCECorrection::forward(const Point& pos_raw, double t0, int face, int apa) const {
100+
Point pc(pos_raw);
101+
corrected_xyz(pos_raw[0], pos_raw[1], pos_raw[2], t0, face, apa, pc[0], pc[1], pc[2]);
102+
return pc;
103+
}
104+
105+
inline Point SCECorrection::backward(const Point& pos_cor, double t0, int face, int apa) const {
106+
Point pr(pos_cor);
107+
if (m_field) {
108+
// Approximate: evaluate field at the corrected position.
109+
pr[0] -= m_field->displacement_x(apa, pos_cor[0], pos_cor[1], pos_cor[2]);
110+
pr[1] -= m_field->displacement_y(apa, pos_cor[0], pos_cor[1], pos_cor[2]);
111+
pr[2] -= m_field->displacement_z(apa, pos_cor[0], pos_cor[1], pos_cor[2]);
112+
}
113+
const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
114+
pr[0] += dirx * t0 * m_drift_speeds.at(apa).at(face);
115+
return pr;
116+
}
117+
118+
inline bool SCECorrection::filter(const Point& pos_cor, double, int face, int apa) const {
119+
auto wpid = m_dv->contained_by(pos_cor);
120+
return wpid.valid() && wpid.apa() == apa && wpid.face() == face;
121+
}
122+
123+
inline PointCloud::Dataset SCECorrection::forward(const PointCloud::Dataset& pc_raw,
124+
const std::vector<std::string>& arr_raw_names,
125+
const std::vector<std::string>& arr_cor_names,
126+
double t0, int face, int apa) const
127+
{
128+
const auto& ax = pc_raw.get(arr_raw_names[0])->elements<double>();
129+
const auto& ay = pc_raw.get(arr_raw_names[1])->elements<double>();
130+
const auto& az = pc_raw.get(arr_raw_names[2])->elements<double>();
131+
std::vector<double> ax_c(ax.size()), ay_c(ax.size()), az_c(ax.size());
132+
for (size_t i = 0; i < ax.size(); ++i) {
133+
corrected_xyz(ax[i], ay[i], az[i], t0, face, apa,
134+
ax_c[i], ay_c[i], az_c[i]);
135+
}
136+
Dataset ds;
137+
ds.add(arr_cor_names[0], Array(ax_c)); // x_sce
138+
ds.add(arr_cor_names[1], Array(ay_c)); // y_sce
139+
ds.add(arr_cor_names[2], Array(az_c)); // z_sce
140+
return ds;
141+
}
142+
143+
inline PointCloud::Dataset SCECorrection::backward(const PointCloud::Dataset& pc_cor,
144+
const std::vector<std::string>& arr_cor_names,
145+
const std::vector<std::string>& arr_raw_names,
146+
double t0, int face, int apa) const
147+
{
148+
const auto& ax = pc_cor.get(arr_cor_names[0])->elements<double>();
149+
const auto& ay = pc_cor.get(arr_cor_names[1])->elements<double>();
150+
const auto& az = pc_cor.get(arr_cor_names[2])->elements<double>();
151+
const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
152+
const double v = m_drift_speeds.at(apa).at(face);
153+
std::vector<double> ax_r(ax.size()), ay_r(ax.size()), az_r(ax.size());
154+
for (size_t i = 0; i < ax.size(); ++i) {
155+
double x = ax[i], y = ay[i], z = az[i];
156+
if (m_field) {
157+
x -= m_field->displacement_x(apa, ax[i], ay[i], az[i]);
158+
y -= m_field->displacement_y(apa, ax[i], ay[i], az[i]);
159+
z -= m_field->displacement_z(apa, ax[i], ay[i], az[i]);
160+
}
161+
x += dirx * t0 * v;
162+
ax_r[i] = x; ay_r[i] = y; az_r[i] = z;
163+
}
164+
Dataset ds;
165+
ds.add(arr_raw_names[0], Array(ax_r));
166+
ds.add(arr_raw_names[1], Array(ay_r));
167+
ds.add(arr_raw_names[2], Array(az_r));
168+
return ds;
169+
}
170+
171+
inline PointCloud::Dataset SCECorrection::filter(const PointCloud::Dataset& pc_cor,
172+
const std::vector<std::string>& arr_cor_names,
173+
double, int face, int apa) const
174+
{
175+
std::vector<int> flt(pc_cor.size_major());
176+
const auto& ax = pc_cor.get(arr_cor_names[0])->elements<double>();
177+
const auto& ay = pc_cor.get(arr_cor_names[1])->elements<double>();
178+
const auto& az = pc_cor.get(arr_cor_names[2])->elements<double>();
179+
for (size_t i = 0; i < ax.size(); ++i) {
180+
flt[i] = 0;
181+
auto wpid = m_dv->contained_by(Point(ax[i], ay[i], az[i]));
182+
if (wpid.valid() && wpid.apa() == apa && wpid.face() == face) flt[i] = 1;
183+
}
184+
Dataset ds;
185+
ds.add("filter", Array(flt));
186+
return ds;
187+
}
188+
189+
} // namespace WireCell::Clus
190+
191+
#endif // WIRECELLCLUS_SCECORRECTION_H

‎clus/src/PCTransforms.cxx‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,8 @@
66
//
77

88
#include "WireCellClus/IPCTransform.h"
9+
#include "WireCellClus/SCECorrection.h"
10+
#include "WireCellIface/ISCEField.h"
911
#include "WireCellIface/IDetectorVolumes.h"
1012

1113
#include "WireCellIface/IConfigurable.h"
@@ -239,6 +241,23 @@ class PCTransformSet : public WireCell::Clus::IPCTransformSet,
239241
auto dv = Factory::find_tn<WireCell::IDetectorVolumes>(dvtn);
240242
const bool relax = get(cfg, "relax_containment_filter", false);
241243
m_pcts["T0Correction"] = std::make_shared<T0Correction>(dv, relax);
244+
// SBND SCE (Space Charge Effect) -- per-TPC TH3 displacement field
245+
// provided externally as an ISCEField (typically WireCell::Root::SCEFieldTH3),
246+
// looked up by TypeName from DetectorVolumes per-APA metadata key "sce_field".
247+
// Falls back to a no-op (T0 only) when no field is configured.
248+
WireCell::ISCEField::pointer sce_field;
249+
for (const auto& [wfid, _] : dv->wpident_faces()) {
250+
WirePlaneId wpid(wfid);
251+
const auto md = dv->metadata(wpid);
252+
if (md.isMember("sce_field")) {
253+
std::string fld_tn = md["sce_field"].asString();
254+
if (!fld_tn.empty()) {
255+
sce_field = Factory::find_tn<WireCell::ISCEField>(fld_tn);
256+
}
257+
break;
258+
}
259+
}
260+
m_pcts["SCECorrection"] = std::make_shared<SCECorrection>(dv, sce_field);
242261
// ...
243262
}
244263

‎iface/inc/WireCellIface/ISCEField.h‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
#ifndef WIRECELLIFACE_ISCEFIELD
2+
#define WIRECELLIFACE_ISCEFIELD
3+
4+
#include "WireCellUtil/IComponent.h"
5+
6+
namespace WireCell {
7+
8+
/** Per-TPC backward (reco -> true) Space-Charge-Effect displacement field.
9+
*
10+
* Implementations expose, for a given APA/TPC index and 3D position,
11+
* the SIGNED backward displacement to be added directly to a t0-corrected
12+
* coordinate to recover the true coordinate. Out-of-domain queries are
13+
* clamped internally by the implementation. Implementations that carry no
14+
* map should return 0 (no-op safe).
15+
*
16+
* Units: WCT internal (mm) in, (mm) out. Implementations handle any
17+
* conversion required by their underlying map.
18+
*
19+
* displacement_x is mandatory. displacement_y / displacement_z are
20+
* optional: implementations whose underlying map lacks transverse
21+
* components inherit the default no-op (return 0).
22+
*/
23+
class ISCEField : public IComponent<ISCEField> {
24+
public:
25+
virtual ~ISCEField() = default;
26+
27+
/// Backward X-displacement at (x,y,z) for APA index `apa`.
28+
/// SBND convention: apa=0 East, apa=1 West. Units: mm.
29+
virtual double displacement_x(int apa, double x, double y, double z) const = 0;
30+
31+
/// Backward Y-displacement at (x,y,z). Units: mm. Optional (default no-op).
32+
virtual double displacement_y(int /*apa*/, double /*x*/, double /*y*/, double /*z*/) const { return 0.0; }
33+
34+
/// Backward Z-displacement at (x,y,z). Units: mm. Optional (default no-op).
35+
virtual double displacement_z(int /*apa*/, double /*x*/, double /*y*/, double /*z*/) const { return 0.0; }
36+
};
37+
38+
} // namespace WireCell
39+
40+
#endif

‎root/inc/WireCellRoot/SCEFieldTH3.h‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
/** ISCEField implementation backed by per-TPC TH3F histograms in a
2+
* ROOT file (SBND dualmap layout: TrueBkwd_Displacement_{X,Y,Z}_{E,W}).
3+
* IConfigurable so jsonnet wiring works the usual way; the actual
4+
* SCE-correction transform in clus/ receives an ISCEField::pointer and
5+
* never sees ROOT.
6+
*/
7+
8+
#ifndef WIRECELLROOT_SCEFIELDTH3
9+
#define WIRECELLROOT_SCEFIELDTH3
10+
11+
#include "WireCellIface/ISCEField.h"
12+
#include "WireCellIface/IConfigurable.h"
13+
#include "WireCellUtil/Logging.h"
14+
15+
#include <memory>
16+
#include <string>
17+
18+
class TFile;
19+
class TH3F;
20+
class TAxis;
21+
22+
namespace WireCell {
23+
namespace Root {
24+
25+
class SCEFieldTH3 : public WireCell::ISCEField,
26+
public WireCell::IConfigurable {
27+
public:
28+
SCEFieldTH3();
29+
virtual ~SCEFieldTH3();
30+
31+
// ISCEField
32+
virtual double displacement_x(int apa, double x, double y, double z) const override;
33+
virtual double displacement_y(int apa, double x, double y, double z) const override;
34+
virtual double displacement_z(int apa, double x, double y, double z) const override;
35+
36+
// IConfigurable
37+
virtual void configure(const WireCell::Configuration& cfg) override;
38+
virtual WireCell::Configuration default_configuration() const override;
39+
40+
private:
41+
// configuration
42+
std::string m_sce_map_file;
43+
std::string m_th3_name_E; // X (mandatory)
44+
std::string m_th3_name_W;
45+
std::string m_th3_name_E_y; // Y (optional transverse)
46+
std::string m_th3_name_W_y;
47+
std::string m_th3_name_E_z; // Z (optional transverse)
48+
std::string m_th3_name_W_z;
49+
bool m_axis_unit_mm{false};
50+
double m_sign{-1.0};
51+
52+
// loaded state (all owned by m_file)
53+
std::unique_ptr<TFile> m_file;
54+
TH3F* m_hE{nullptr};
55+
TH3F* m_hW{nullptr};
56+
TH3F* m_hE_y{nullptr};
57+
TH3F* m_hW_y{nullptr};
58+
TH3F* m_hE_z{nullptr};
59+
TH3F* m_hW_z{nullptr};
60+
61+
Log::logptr_t l;
62+
63+
// Shared: pick E/W histo by apa, clamp, convert units, apply sign.
64+
// Returns 0 if the requested histo is absent (no-op component).
65+
double interp(TH3F* hE, TH3F* hW, int apa,
66+
double x, double y, double z) const;
67+
68+
static double clamp_axis(const TAxis* a, double v);
69+
70+
// Load one (E,W) histo pair; hard=true throws on missing,
71+
// hard=false warns and leaves pointers null.
72+
void load_pair(const std::string& nameE, const std::string& nameW,
73+
TH3F*& hE, TH3F*& hW, bool hard);
74+
};
75+
76+
} // namespace Root
77+
} // namespace WireCell
78+
79+
#endif

0 commit comments

Comments
 (0)