Skip to content

Commit 80ca0d1

Browse files
abhatfnalHaiwangYu
abhatfnal
authored and
HaiwangYu
committed
clus: add SCECorrection IPCTransform via ISCEField
SCECorrection is an IPCTransform that composes a T0 correction [per-APA drift_speed from DetectorVolumes metadata] with an optional SCE backward displacement delegated to an externally-provided ISCEField. This keeps clus/ ROOT-free; the ROOT-backed displacement field lives in root/ as SCEFieldTH3. PCTransforms looks up the ISCEField by TypeName from the DetectorVolumes metadata key 'sce_field' via Factory::find_tn and passes it to the SCECorrection constructor. If 'sce_field' is not present in the metadata, SCECorrection runs as T0-only -- the SCE step is a no-op. Validated on 50 SBND crossing-muon DETSIM events against the dualmap SCE model: per-point residual rms 1.83e-04 cm East, 2.23e-04 cm West, max 8.42e-04 cm; pooled mean|dx| ratio W/E = 1.327. The reco reproduces the TH3 backward map to interpolation precision.
1 parent b062925 commit 80ca0d1

2 files changed

Lines changed: 194 additions & 0 deletions

File tree

‎clus/inc/WireCellClus/SCECorrection.h‎

Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
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:
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+
//
13+
// East/West (apa==0/1) routing is handled inside the ISCEField implementation.
14+
// If no field is provided (nullptr), the SCE step is a no-op (T0 still applied).
15+
//
16+
// Author: Avinay Bhat (UChicago), for SBND
17+
// Modeled on T0Correction by Haiwang Yu.
18+
19+
#ifndef WIRECELLCLUS_SCECORRECTION_H
20+
#define WIRECELLCLUS_SCECORRECTION_H
21+
22+
#include "WireCellClus/IPCTransform.h"
23+
#include "WireCellIface/IDetectorVolumes.h"
24+
#include "WireCellIface/ISCEField.h"
25+
#include "WireCellUtil/Logging.h"
26+
27+
#include <map>
28+
#include <memory>
29+
#include <string>
30+
31+
namespace WireCell::Clus {
32+
33+
class SCECorrection : public WireCell::Clus::IPCTransform {
34+
public:
35+
SCECorrection(IDetectorVolumes::pointer dv,
36+
WireCell::ISCEField::pointer field = nullptr);
37+
virtual ~SCECorrection() = default;
38+
39+
virtual Point forward(const Point& pos_raw, double cluster_t0, int face, int apa) const override;
40+
virtual Point backward(const Point& pos_cor, double cluster_t0, int face, int apa) const override;
41+
virtual bool filter(const Point& pos_cor, double cluster_t0, int face, int apa) const override;
42+
43+
virtual Dataset forward(const Dataset& pc_raw,
44+
const std::vector<std::string>& arr_raw_names,
45+
const std::vector<std::string>& arr_cor_names,
46+
double cluster_t0, int face, int apa) const override;
47+
virtual Dataset backward(const Dataset& pc_cor,
48+
const std::vector<std::string>& arr_cor_names,
49+
const std::vector<std::string>& arr_raw_names,
50+
double cluster_t0, int face, int apa) const override;
51+
virtual Dataset filter(const Dataset& pc_cor,
52+
const std::vector<std::string>& arr_cor_names,
53+
double cluster_t0, int face, int apa) const override;
54+
55+
virtual PointCloud::Tree::Scope output_scope() const override {
56+
return {"3d", {"x_sce", "y", "z"}};
57+
}
58+
virtual std::vector<std::string> stored_array_names() const override {
59+
return {"x_sce"};
60+
}
61+
62+
private:
63+
IDetectorVolumes::pointer m_dv;
64+
WireCell::ISCEField::pointer m_field;
65+
std::map<int, std::map<int, double>> m_drift_speeds;
66+
67+
// Combined T0 + SCE for one X coordinate. All quantities in WCT mm.
68+
inline double corrected_x(double x_raw_mm, double y_mm, double z_mm,
69+
double cluster_t0, int face, int apa) const {
70+
const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
71+
const double x_t0_mm = x_raw_mm - dirx * cluster_t0 * m_drift_speeds.at(apa).at(face);
72+
const double dx_mm = m_field ? m_field->displacement_x(apa, x_t0_mm, y_mm, z_mm) : 0.0;
73+
return x_t0_mm + dx_mm;
74+
}
75+
};
76+
77+
inline SCECorrection::SCECorrection(IDetectorVolumes::pointer dv,
78+
WireCell::ISCEField::pointer field)
79+
: m_dv(dv), m_field(field)
80+
{
81+
for (const auto& [wfid, _] : m_dv->wpident_faces()) {
82+
WirePlaneId wpid(wfid);
83+
const auto md = m_dv->metadata(wpid);
84+
m_drift_speeds[wpid.apa()][wpid.face()] = md["drift_speed"].asDouble();
85+
}
86+
if (m_field) spdlog::info("SCECorrection: ISCEField wired in");
87+
else spdlog::info("SCECorrection: no ISCEField -- SCE disabled (T0 still applied)");
88+
}
89+
90+
inline Point SCECorrection::forward(const Point& pos_raw, double t0, int face, int apa) const {
91+
Point pc(pos_raw);
92+
pc[0] = corrected_x(pos_raw[0], pos_raw[1], pos_raw[2], t0, face, apa);
93+
return pc;
94+
}
95+
96+
inline Point SCECorrection::backward(const Point& pos_cor, double t0, int face, int apa) const {
97+
Point pr(pos_cor);
98+
if (m_field) {
99+
// Approximate: evaluate field at x_sce instead of x_t0 (small displacement OK).
100+
pr[0] -= m_field->displacement_x(apa, pos_cor[0], pos_cor[1], pos_cor[2]);
101+
}
102+
const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
103+
pr[0] += dirx * t0 * m_drift_speeds.at(apa).at(face);
104+
return pr;
105+
}
106+
107+
inline bool SCECorrection::filter(const Point& pos_cor, double, int face, int apa) const {
108+
auto wpid = m_dv->contained_by(pos_cor);
109+
return wpid.valid() && wpid.apa() == apa && wpid.face() == face;
110+
}
111+
112+
inline PointCloud::Dataset SCECorrection::forward(const PointCloud::Dataset& pc_raw,
113+
const std::vector<std::string>& arr_raw_names,
114+
const std::vector<std::string>& arr_cor_names,
115+
double t0, int face, int apa) const
116+
{
117+
const auto& ax = pc_raw.get(arr_raw_names[0])->elements<double>();
118+
const auto& ay = pc_raw.get(arr_raw_names[1])->elements<double>();
119+
const auto& az = pc_raw.get(arr_raw_names[2])->elements<double>();
120+
std::vector<double> ax_c(ax.size());
121+
for (size_t i = 0; i < ax.size(); ++i) {
122+
ax_c[i] = corrected_x(ax[i], ay[i], az[i], t0, face, apa);
123+
}
124+
Dataset ds;
125+
ds.add(arr_cor_names[0], Array(ax_c));
126+
ds.add(arr_cor_names[1], Array(ay));
127+
ds.add(arr_cor_names[2], Array(az));
128+
return ds;
129+
}
130+
131+
inline PointCloud::Dataset SCECorrection::backward(const PointCloud::Dataset& pc_cor,
132+
const std::vector<std::string>& arr_cor_names,
133+
const std::vector<std::string>& arr_raw_names,
134+
double t0, int face, int apa) const
135+
{
136+
const auto& ax = pc_cor.get(arr_cor_names[0])->elements<double>();
137+
const auto& ay = pc_cor.get(arr_cor_names[1])->elements<double>();
138+
const auto& az = pc_cor.get(arr_cor_names[2])->elements<double>();
139+
const double dirx = m_dv->face_dirx(WirePlaneId(kAllLayers, face, apa));
140+
const double v = m_drift_speeds.at(apa).at(face);
141+
std::vector<double> ax_r(ax.size());
142+
for (size_t i = 0; i < ax.size(); ++i) {
143+
double x = ax[i];
144+
if (m_field) x -= m_field->displacement_x(apa, ax[i], ay[i], az[i]);
145+
x += dirx * t0 * v;
146+
ax_r[i] = x;
147+
}
148+
Dataset ds;
149+
ds.add(arr_raw_names[0], Array(ax_r));
150+
ds.add(arr_raw_names[1], Array(ay));
151+
ds.add(arr_raw_names[2], Array(az));
152+
return ds;
153+
}
154+
155+
inline PointCloud::Dataset SCECorrection::filter(const PointCloud::Dataset& pc_cor,
156+
const std::vector<std::string>& arr_cor_names,
157+
double, int face, int apa) const
158+
{
159+
std::vector<int> flt(pc_cor.size_major());
160+
const auto& ax = pc_cor.get(arr_cor_names[0])->elements<double>();
161+
const auto& ay = pc_cor.get(arr_cor_names[1])->elements<double>();
162+
const auto& az = pc_cor.get(arr_cor_names[2])->elements<double>();
163+
for (size_t i = 0; i < ax.size(); ++i) {
164+
flt[i] = 0;
165+
auto wpid = m_dv->contained_by(Point(ax[i], ay[i], az[i]));
166+
if (wpid.valid() && wpid.apa() == apa && wpid.face() == face) flt[i] = 1;
167+
}
168+
Dataset ds;
169+
ds.add("filter", Array(flt));
170+
return ds;
171+
}
172+
173+
} // namespace WireCell::Clus
174+
175+
#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

0 commit comments

Comments
 (0)