diff --git a/README.md b/README.md
index 5806392a4..ac65f82ee 100644
--- a/README.md
+++ b/README.md
@@ -155,6 +155,7 @@ EIC-Opticks provides several examples demonstrating GPU-accelerated optical phot
| `GPUPhotonSource` | Optical photons (torch) | Any GDML | G4 + GPU side-by-side validation |
| `GPUPhotonSourceMinimal` | Optical photons (torch) | Any GDML | GPU-only test |
| `GPUPhotonFileSource` | Optical photons (text file) | Any GDML | GPU-only, user-defined photons from file |
+| WLS test | Wavelength shifting | WLS sphere + detector shell | Validate GPU WLS physics |
### Example 1: GPUCerenkov (Cerenkov Only)
@@ -297,6 +298,55 @@ GPUPhotonFileSource -g tests/geom/opticks_raindrop.gdml -p my_photons.txt -m run
**Source files:** `src/GPUPhotonFileSource.cpp`, `src/GPUPhotonFileSource.h`
+### Example 6: Wavelength Shifting (WLS) Test
+
+This test validates the GPU wavelength shifting implementation using a dedicated
+geometry with a WLS sphere surrounded by a detector shell:
+
+```
+Geometry: wls_test.gdml
+├── Air world (r=200 mm)
+│ ├── WLS sphere (r=20 mm) ← Absorbs UV, re-emits visible
+│ └── Glass detector shell (r=28-30 mm) ← 100% detection efficiency
+```
+
+The WLS material absorbs UV photons (350 nm) and re-emits them isotropically at
+longer wavelengths (peak ~481 nm) with a 0.5 ns exponential time delay. The test
+fires 1000 monochromatic 350 nm photons from the origin into the WLS sphere.
+
+```bash
+GPUPhotonSourceMinimal -g tests/geom/wls_test.gdml -c wls_test -m tests/run.mac -s 42
+```
+
+**Expected results:**
+- ~990/1000 photons detected (10 absorbed after failing energy conservation)
+- All hits wavelength-shifted from 350 nm to mean ~487 nm
+- Energy conservation: no hits with wavelength < 350 nm
+- Isotropic re-emission: mean momentum direction near zero
+- Time delay: mean ~0.6 ns (propagation + 0.5 ns exponential WLS decay)
+
+**GDML WLS properties required** (same syntax for G4 10.x and 11.x):
+```xml
+
+
+
+
+
+
+
+
+
+
+
+
+```
+
+Unlike scintillation properties, WLS property names are the same in both Geant4
+10.x and 11.x — no dual-naming is needed.
+
+**Test files:** `tests/geom/wls_test.gdml`, `config/wls_test.json`
+**Implementation docs:** `docs/WLS_IMPLEMENTATION.md`
+
### Torch configuration
`GPUPhotonSource` and `GPUPhotonSourceMinimal` read photon source parameters from a
@@ -317,6 +367,7 @@ JSON config file (default `config/dev.json`). Key fields:
|---------|-------------|-------------|-----------------|----------------------|---------------------|
| Cerenkov genstep collection | ✓ | ✓ | ✗ | ✗ | ✗ |
| Scintillation genstep collection | ✗ | ✓ | ✗ | ✗ | ✗ |
+| Wavelength shifting (WLS) | ✓ | ✓ | ✓ | ✓ | ✓ |
| Torch photon generation | ✗ | ✗ | ✓ | ✓ | ✗ |
| Photon input from text file | ✗ | ✗ | ✗ | ✗ | ✓ |
| G4 optical photon tracking | ✓ | ✓ | ✓ | ✗ | ✗ |
diff --git a/config/wls_test.json b/config/wls_test.json
new file mode 100644
index 000000000..b8572b5b9
--- /dev/null
+++ b/config/wls_test.json
@@ -0,0 +1,30 @@
+{
+ "torch": {
+ "gentype": "TORCH",
+ "trackid": 0,
+ "matline": 0,
+ "numphoton": 1000,
+
+ "pos": [0.0, 0.0, 0.0],
+ "time": 0.0,
+
+ "mom": [0.0, 0.0, 1.0],
+ "weight": 0.0,
+
+ "pol": [1.0, 0.0, 0.0],
+ "wavelength": 350.0,
+
+ "zenith": [0.0, 1.0],
+ "azimuth": [0.0, 1.0],
+
+ "radius": 0.0,
+ "distance": 0.0,
+ "mode": 255,
+ "type": "disc"
+ },
+
+ "event": {
+ "mode": "DebugLite",
+ "maxslot": 1000000
+ }
+}
diff --git a/qudarap/CMakeLists.txt b/qudarap/CMakeLists.txt
index 72c76bd41..9a7e34673 100644
--- a/qudarap/CMakeLists.txt
+++ b/qudarap/CMakeLists.txt
@@ -48,6 +48,8 @@ set(SOURCES
QScint.cc
QScint.cu
+ QWls.cc
+
QCerenkovIntegral.cc
QCerenkov.cc
QCerenkov.cu
@@ -117,6 +119,9 @@ SET(HEADERS
QScint.hh
qscint.h
+ QWls.hh
+ qwls.h
+
QCerenkovIntegral.hh
QCerenkov.hh
qcerenkov.h
diff --git a/qudarap/QSim.cc b/qudarap/QSim.cc
index 3731e664a..da5110486 100644
--- a/qudarap/QSim.cc
+++ b/qudarap/QSim.cc
@@ -32,19 +32,19 @@
#include "qdebug.h"
#include "QBase.hh"
-#include "QEvt.hh"
-#include "QRng.hh"
-#include "QTex.hh"
-#include "QScint.hh"
-#include "QCerenkov.hh"
#include "QBnd.hh"
-#include "QProp.hh"
-#include "QMultiFilm.hh"
+#include "QCerenkov.hh"
+#include "QDebug.hh"
#include "QEvt.hh"
+#include "QMultiFilm.hh"
#include "QOptical.hh"
-#include "QSimLaunch.hh"
-#include "QDebug.hh"
#include "QPMT.hh"
+#include "QProp.hh"
+#include "QRng.hh"
+#include "QScint.hh"
+#include "QSimLaunch.hh"
+#include "QTex.hh"
+#include "QWls.hh"
#include "QSim.hh"
@@ -179,6 +179,27 @@ void QSim::UploadComponents( const SSim* ssim )
LOG(LEVEL) << scint->desc();
}
+ const NP *wls_icdf = ssim->get(snam::WLS_ICDF);
+ const NP *wls_mat_map = ssim->get(snam::WLS_MAT_MAP);
+ if (wls_icdf == nullptr || wls_mat_map == nullptr)
+ {
+ LOG(LEVEL) << " wls_icdf or wls_mat_map null — no WLS materials in geometry ";
+ }
+ else
+ {
+ const NP *wls_tc = ssim->get(snam::WLS_TIME_CONSTANTS);
+ if (wls_tc)
+ {
+ unsigned hd_factor = 20u;
+ QWls *qwls_ = new QWls(wls_icdf, wls_mat_map, wls_tc, hd_factor);
+ LOG(LEVEL) << qwls_->desc();
+ }
+ else
+ {
+ LOG(error) << " wls_icdf and wls_mat_map present but wls_time_constants missing ";
+ }
+ }
+
// TODO: make this more like the others : acting on the available inputs rather than the mode
bool is_simtrace = SEventConfig::IsRGModeSimtrace() ;
if(is_simtrace == false )
@@ -258,13 +279,13 @@ singleton components.
**/
-QSim::QSim()
- :
+QSim::QSim() :
base(QBase::Get()),
qev(new QEvt),
sev(qev->sev),
rng(QRng::Get()),
scint(QScint::Get()),
+ qwls(QWls::Get()),
cerenkov(QCerenkov::Get()),
bnd(QBnd::Get()),
debug_(QDebug::Get()),
@@ -314,6 +335,7 @@ void QSim::init()
sim->multifilm = multifilm ? multifilm->d_multifilm : nullptr ;
sim->cerenkov = cerenkov ? cerenkov->d_cerenkov : nullptr ;
sim->scint = scint ? scint->d_scint : nullptr ;
+ sim->wls = qwls ? qwls->d_wls : nullptr;
sim->pmt = pmt ? pmt->d_pmt : nullptr ;
diff --git a/qudarap/QSim.hh b/qudarap/QSim.hh
index de92f7c23..3e59a9c9c 100644
--- a/qudarap/QSim.hh
+++ b/qudarap/QSim.hh
@@ -38,6 +38,7 @@ struct QBase ;
struct QEvt ;
struct QRng ;
struct QScint ;
+struct QWls;
struct QCerenkov ;
struct QBnd ;
struct QMultiFilm;
@@ -74,6 +75,7 @@ struct QUDARAP_API QSim
const QRng* rng ;
const QScint* scint ;
+ const QWls *qwls;
const QCerenkov* cerenkov ;
const QBnd* bnd ;
const QOptical* optical ;
diff --git a/qudarap/QU.cc b/qudarap/QU.cc
index 97aacf985..3ba65bcf5 100644
--- a/qudarap/QU.cc
+++ b/qudarap/QU.cc
@@ -29,15 +29,15 @@
#include "qsim.h"
#include "qbase.h"
-#include "qprop.h"
-#include "qpmt.h"
-#include "qdebug.h"
-#include "qscint.h"
#include "qcerenkov.h"
#include "qcurandwrap.h"
-#include "scurandref.h"
+#include "qdebug.h"
#include "qmultifilm.h"
-
+#include "qpmt.h"
+#include "qprop.h"
+#include "qscint.h"
+#include "qwls.h"
+#include "scurandref.h"
const plog::Severity QU::LEVEL = SLOG::EnvLevel("QU", "DEBUG") ;
bool QU::MEMCHECK = ssys::getenvbool(_MEMCHECK);
@@ -171,6 +171,7 @@ template qbnd* QU::UploadArray(const qbnd* array, unsigned num_it
template sevent* QU::UploadArray(const sevent* array, unsigned num_items, const char* label) ;
template qdebug* QU::UploadArray(const qdebug* array, unsigned num_items, const char* label) ;
template qscint* QU::UploadArray(const qscint* array, unsigned num_items, const char* label) ;
+template qwls *QU::UploadArray(const qwls *array, unsigned num_items, const char *label);
template qcerenkov* QU::UploadArray(const qcerenkov* array, unsigned num_items, const char* label) ;
template qbase* QU::UploadArray(const qbase* array, unsigned num_items, const char* label) ;
diff --git a/qudarap/QWls.cc b/qudarap/QWls.cc
new file mode 100644
index 000000000..62a83e404
--- /dev/null
+++ b/qudarap/QWls.cc
@@ -0,0 +1,132 @@
+#include
+#include
+#include
+
+#include "scuda.h"
+#include "squad.h"
+
+#include "NP.hh"
+#include "SLOG.hh"
+#include "ssys.h"
+
+#include "QTex.hh"
+#include "QU.hh"
+#include "QUDA_CHECK.h"
+#include "QWls.hh"
+
+#include "qwls.h"
+
+const plog::Severity QWls::LEVEL = SLOG::EnvLevel("QWls", "DEBUG");
+
+const QWls *QWls::INSTANCE = nullptr;
+const QWls *QWls::Get()
+{
+ return INSTANCE;
+}
+
+/**
+QWls::QWls
+------------
+
+1. Narrows ICDF from double to float if needed
+2. Uploads ICDF into GPU texture
+3. Creates qwls instance with device pointers and uploads it
+
+**/
+
+QWls::QWls(const NP *wls_icdf, const NP *mat_map, const NP *time_constants, unsigned hd_factor) :
+ dsrc(wls_icdf->ebyte == 8 ? wls_icdf : nullptr),
+ src(wls_icdf->ebyte == 4 ? wls_icdf : NP::MakeNarrow(dsrc)),
+ tex(MakeWlsTex(src, hd_factor)),
+ wls(MakeInstance(tex, mat_map, time_constants, hd_factor, time_constants->shape[0])),
+ d_wls(QU::UploadArray(wls, 1, "QWls::QWls/d_wls"))
+{
+ INSTANCE = this;
+}
+
+/**
+QWls::MakeWlsTex
+-------------------
+
+Creates a 2D CUDA texture from the ICDF array.
+Shape: (num_wls*3, 4096, 1) where 3 = HD layers per material.
+
+**/
+
+QTex *QWls::MakeWlsTex(const NP *src, unsigned hd_factor)
+{
+ assert(src);
+ assert(src->shape.size() == 3);
+
+ unsigned ni = src->shape[0]; // height: num_wls * 3
+ unsigned nj = src->shape[1]; // width: 4096
+ unsigned nk = src->shape[2]; // 1
+
+ assert(nk == 1);
+ assert(nj == 4096);
+ assert(ni % 3 == 0); // must be multiple of 3 (3 HD layers per material)
+ assert(src->uifc == 'f' && src->ebyte == 4);
+
+ unsigned ny = ni; // height
+ unsigned nx = nj; // width
+
+ bool normalizedCoords = true;
+ QTex *tx = new QTex(nx, ny, src->cvalues(), 'L', normalizedCoords, src);
+
+ tx->setHDFactor(hd_factor);
+ tx->uploadMeta();
+
+ LOG(LEVEL) << " src " << src->desc() << " nx (width) " << nx << " ny (height) " << ny << " tx.HDFactor "
+ << tx->getHDFactor();
+
+ return tx;
+}
+
+/**
+QWls::MakeInstance
+---------------------
+
+Creates the host-side qwls struct populated with device pointers.
+Uploads material_map and time_constants to device memory.
+
+**/
+
+qwls *QWls::MakeInstance(const QTex *tex, const NP *mat_map, const NP *time_constants, unsigned hd_factor,
+ unsigned num_wls)
+{
+ assert(mat_map);
+ assert(time_constants);
+ assert(mat_map->uifc == 'i' && mat_map->ebyte == 4);
+ assert(time_constants->uifc == 'f' && time_constants->ebyte == 4);
+
+ qwls *w = new qwls;
+ w->wls_tex = tex->texObj;
+ w->hd_factor = hd_factor;
+ w->num_wls = num_wls;
+ w->tex_height = tex->height;
+
+ // Upload material_map to device
+ unsigned num_mat = mat_map->shape[0];
+ int *d_mat_map = nullptr;
+ size_t mat_map_size = num_mat * sizeof(int);
+ QUDA_CHECK(cudaMalloc(reinterpret_cast(&d_mat_map), mat_map_size));
+ QUDA_CHECK(cudaMemcpy(d_mat_map, mat_map->cvalues(), mat_map_size, cudaMemcpyHostToDevice));
+ w->material_map = d_mat_map;
+
+ // Upload time_constants to device
+ float *d_tc = nullptr;
+ size_t tc_size = num_wls * sizeof(float);
+ QUDA_CHECK(cudaMalloc(reinterpret_cast(&d_tc), tc_size));
+ QUDA_CHECK(cudaMemcpy(d_tc, time_constants->cvalues(), tc_size, cudaMemcpyHostToDevice));
+ w->time_constants = d_tc;
+
+ return w;
+}
+
+std::string QWls::desc() const
+{
+ std::stringstream ss;
+ ss << "QWls" << " dsrc " << (dsrc ? dsrc->desc() : "-") << " src " << (src ? src->desc() : "-") << " tex "
+ << (tex ? tex->desc() : "-");
+ return ss.str();
+}
diff --git a/qudarap/QWls.hh b/qudarap/QWls.hh
new file mode 100644
index 000000000..3134eba2b
--- /dev/null
+++ b/qudarap/QWls.hh
@@ -0,0 +1,41 @@
+#pragma once
+
+#include "QUDARAP_API_EXPORT.hh"
+#include "plog/Severity.h"
+#include
+
+struct NP;
+template struct QTex;
+struct qwls;
+
+/**
+QWls : Host-side WLS ICDF Texture Upload
+============================================
+
+Uploads the WLS inverse CDF array into a GPU texture and creates
+the device-side qwls struct with material mapping and time constants.
+
+Follows the same pattern as QScint for scintillation ICDF textures.
+
+**/
+
+struct QUDARAP_API QWls
+{
+ static const plog::Severity LEVEL;
+ static const QWls *INSTANCE;
+ static const QWls *Get();
+
+ static QTex *MakeWlsTex(const NP *src, unsigned hd_factor);
+ static qwls *MakeInstance(const QTex *tex, const NP *mat_map, const NP *time_constants, unsigned hd_factor,
+ unsigned num_wls);
+
+ const NP *dsrc; // original double-precision ICDF
+ const NP *src; // narrowed float ICDF
+ QTex *tex; // GPU texture
+ qwls *wls; // host-side instance (with device pointers)
+ qwls *d_wls; // device copy of qwls struct
+
+ QWls(const NP *wls_icdf, const NP *mat_map, const NP *time_constants, unsigned hd_factor);
+
+ std::string desc() const;
+};
diff --git a/qudarap/qsim.h b/qudarap/qsim.h
index 39ecd6caf..a4b7143d6 100644
--- a/qudarap/qsim.h
+++ b/qudarap/qsim.h
@@ -55,17 +55,17 @@ Canonical use is from CSGOptiX/CSGOptiX7.cu:simulate
#include "sctx.h"
-#include "qrng.h"
#include "qbase.h"
-#include "qprop.h"
-#include "qmultifilm.h"
#include "qbnd.h"
-#include "qscint.h"
#include "qcerenkov.h"
+#include "qmultifilm.h"
#include "qpmt.h"
+#include "qprop.h"
+#include "qrng.h"
+#include "qscint.h"
+#include "qwls.h"
#include "tcomplex.h"
-
struct qcerenkov ;
struct qsim
@@ -77,6 +77,7 @@ struct qsim
qmultifilm* multifilm;
qcerenkov* cerenkov ;
qscint* scint ;
+ qwls *wls;
qpmt* pmt ;
#if defined(__CUDACC__) || defined(__CUDABE__)
@@ -136,18 +137,19 @@ struct qsim
// CTOR
#if defined(__CUDACC__) || defined(__CUDABE__)
#else
-inline qsim::qsim() // instanciated on CPU (see QSim::init_sim) and copied to device so no ctor in device code
- :
- base(nullptr),
- evt(nullptr),
- rng(nullptr),
- bnd(nullptr),
- multifilm(nullptr),
- cerenkov(nullptr),
- scint(nullptr),
- pmt(nullptr)
- {
- }
+inline qsim::qsim() // instanciated on CPU (see QSim::init_sim) and copied to device so no ctor in device code
+ :
+ base(nullptr),
+ evt(nullptr),
+ rng(nullptr),
+ bnd(nullptr),
+ multifilm(nullptr),
+ cerenkov(nullptr),
+ scint(nullptr),
+ wls(nullptr),
+ pmt(nullptr)
+{
+}
#endif
inline QSIM_METHOD void qsim::generate_photon_dummy(sphoton& p_, RNG& rng, const quad6& gs, unsigned long long photon_id, unsigned genstep_id ) const
@@ -215,7 +217,6 @@ inline QSIM_METHOD float qsim::RandGaussQ_shoot( RNG& rng, float mean, float std
return v ;
}
-
/**
qsim::SmearNormal_SigmaAlpha
------------------------------
@@ -719,6 +720,7 @@ inline QSIM_METHOD int qsim::propagate_to_boundary(unsigned& flag, RNG& rng, sct
const float& scattering_length = s.material1.z ;
const float& reemission_prob = s.material1.w ;
const float& group_velocity = ctx.current_group_velocity;
+ const float &wls_absorption_length = s.m1group2.z;
const float& distance_to_boundary = ctx.prd->q0.f.w ;
@@ -728,6 +730,7 @@ inline QSIM_METHOD int qsim::propagate_to_boundary(unsigned& flag, RNG& rng, sct
#endif
float u_scattering = curand_uniform(&rng) ;
float u_absorption = curand_uniform(&rng) ;
+ float u_wls_absorption = (wls != nullptr) ? curand_uniform(&rng) : 2.f;
#if !defined(PRODUCTION) && defined(DEBUG_TAG)
stagr& tagr = ctx.tagr ;
@@ -742,9 +745,11 @@ inline QSIM_METHOD int qsim::propagate_to_boundary(unsigned& flag, RNG& rng, sct
// see notes/issues/U4LogTest_maybe_replacing_G4Log_G4UniformRand_in_Absorption_and_Scattering_with_float_version_will_avoid_deviations.rst
float scattering_distance = -scattering_length*KLUDGE_FASTMATH_LOGF(u_scattering);
float absorption_distance = -absorption_length*KLUDGE_FASTMATH_LOGF(u_absorption);
+ float wls_absorption_distance = -wls_absorption_length * KLUDGE_FASTMATH_LOGF(u_wls_absorption);
#else
float scattering_distance = -scattering_length*logf(u_scattering);
float absorption_distance = -absorption_length*logf(u_absorption);
+ float wls_absorption_distance = -wls_absorption_length * logf(u_wls_absorption);
#endif
#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
@@ -766,11 +771,82 @@ inline QSIM_METHOD int qsim::propagate_to_boundary(unsigned& flag, RNG& rng, sct
}
#endif
+ // WLS absorption only competes in materials that actually have WLS properties.
+ // Without this gate, non-WLS materials would create an extra absorption channel
+ // whenever the WLS distance happened to be the shortest of the three.
+ unsigned mat_idx = s.index.x - 1u; // 0-based material index from 1-based optical index
+ bool has_wls_material = (wls != nullptr) && wls->has_wls(mat_idx);
+
+ bool wls_wins = has_wls_material && wls_absorption_distance <= absorption_distance &&
+ wls_absorption_distance <= scattering_distance;
+ if (wls_wins && wls_absorption_distance <= distance_to_boundary)
+ {
+ // WLS ABSORPTION: photon absorbed by wavelength shifting material
+ p.time += wls_absorption_distance / group_velocity;
+ p.pos += wls_absorption_distance * (p.mom);
+
+ // Sample re-emitted wavelength from WLS emission spectrum ICDF
+ float u_wls_wl = curand_uniform(&rng);
+ float new_wavelength = wls->wavelength(mat_idx, u_wls_wl);
+
+ // Energy conservation: re-emitted photon must have lower energy (longer wavelength).
+ // Matches G4OpWLS algorithm: retry up to 100 times.
+ int attempts = 0;
+ while (new_wavelength < p.wavelength && attempts < 100)
+ {
+ u_wls_wl = curand_uniform(&rng);
+ new_wavelength = wls->wavelength(mat_idx, u_wls_wl);
+ attempts++;
+ }
+ if (new_wavelength < p.wavelength)
+ {
+ // Failed energy conservation after 100 attempts — absorb without re-emission
+ flag = BULK_ABSORB;
+ return BREAK;
+ }
+ p.wavelength = new_wavelength;
+
+ // Isotropic re-emission direction and polarization perpendicular by construction.
+ // Same pattern as qscint::mom_pol_wavelength to avoid the near-parallel
+ // random-cross-product instability.
+ float u_wls_cost = curand_uniform(&rng);
+ float u_wls_phi = curand_uniform(&rng);
+ float u_wls_pol = curand_uniform(&rng);
+
+ float cost = 1.f - 2.f * u_wls_cost;
+ float sint = sqrtf((1.f - cost) * (1.f + cost));
+ float phi = 2.f * M_PIf * u_wls_phi;
+ float sinp = sinf(phi);
+ float cosp = cosf(phi);
+
+ p.mom.x = sint * cosp;
+ p.mom.y = sint * sinp;
+ p.mom.z = cost;
+
+ p.pol.x = cost * cosp;
+ p.pol.y = cost * sinp;
+ p.pol.z = -sint;
+
+ phi = 2.f * M_PIf * u_wls_pol;
+ sinp = sinf(phi);
+ cosp = cosf(phi);
+ p.pol = normalize(cosp * p.pol + sinp * cross(p.mom, p.pol));
+
+ // Apply WLS time delay (exponential decay)
+ float tc = wls->time_constant(mat_idx);
+ if (tc > 0.f)
+ {
+ float u_wls_time = curand_uniform(&rng);
+ p.time += -tc * logf(u_wls_time);
+ }
- if (absorption_distance <= scattering_distance)
+ flag = BULK_REEMIT;
+ return CONTINUE;
+ }
+ else if (absorption_distance <= scattering_distance)
{
if (absorption_distance <= distance_to_boundary)
{
@@ -2139,7 +2215,6 @@ inline QSIM_METHOD int qsim::propagate(const int bounce, RNG& rng, sctx& ctx )
// Q: Does flag need to be single bit at this point OR can multiple "flags" be OR-ed together here ?
// A: Decided to keep the flag as single bitted, and directly set EFFICENCY_COLLECT/CULL into ctx.p.flagmask
-
#if !defined(PRODUCTION) && defined(DEBUG_PIDX)
if( ctx.pidx == base->pidx )
printf("//qsim.propagate.tail pidx %7lld bounce %d command %d flag %d ctx.s.optical.y(ems) %d \n",
diff --git a/qudarap/qwls.h b/qudarap/qwls.h
new file mode 100644
index 000000000..0e4100726
--- /dev/null
+++ b/qudarap/qwls.h
@@ -0,0 +1,148 @@
+#pragma once
+/**
+qwls.h : GPU-side Wavelength Shifting
+=========================================
+
+Device-side struct for WLS wavelength sampling via ICDF texture lookup.
+Supports multiple WLS materials indexed by material ID.
+
+The ICDF texture layout:
+- Each WLS material occupies 3 rows (standard, LHS HD, RHS HD)
+- material_map[mat_idx] gives the base row for that material (-1 = no WLS)
+- time_constants[wls_idx] gives the re-emission time constant in ns
+
+Wavelength sampling uses the same HD (high-definition) technique as qscint.h:
+- hd_factor=20: 20x resolution at extremes (u < 0.05 or u > 0.95)
+- Normalized texture coordinates with linear interpolation
+
+**/
+
+#if defined(__CUDACC__) || defined(__CUDABE__)
+#define QWLS_METHOD __device__
+#else
+#define QWLS_METHOD
+#endif
+
+struct qwls
+{
+ cudaTextureObject_t wls_tex; // ICDF texture: (num_wls*3, 4096, 1)
+ unsigned hd_factor; // 0, 10, or 20
+ int *material_map; // device ptr: mat_idx -> base ICDF row (-1 = no WLS)
+ float *time_constants; // device ptr: per-WLS-material time constant (ns)
+ unsigned num_wls; // number of WLS materials
+ unsigned tex_height; // total rows in texture = num_wls * 3
+
+#if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
+
+ QWLS_METHOD bool has_wls(unsigned mat_idx) const;
+ QWLS_METHOD float wavelength(unsigned mat_idx, const float &u0) const;
+ QWLS_METHOD float wavelength_at_row(unsigned base_row, const float &u0) const;
+ QWLS_METHOD float time_constant(unsigned mat_idx) const;
+
+#endif
+};
+
+#if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_CURAND) || defined(MOCK_CUDA)
+
+/**
+qwls::has_wls
+---------------
+
+Returns true if material at mat_idx has WLS properties.
+The material_map holds -1 for non-WLS materials.
+
+**/
+
+inline QWLS_METHOD bool qwls::has_wls(unsigned mat_idx) const
+{
+ return material_map[mat_idx] >= 0;
+}
+
+/**
+qwls::time_constant
+---------------------
+
+Returns the WLS re-emission time constant in ns for the given material.
+Returns 0.f if material has no WLS (instant re-emission / no delay).
+
+**/
+
+inline QWLS_METHOD float qwls::time_constant(unsigned mat_idx) const
+{
+ int base_row = material_map[mat_idx];
+ if (base_row < 0)
+ return 0.f;
+ unsigned wls_idx = base_row / 3;
+ return time_constants[wls_idx];
+}
+
+/**
+qwls::wavelength
+-------------------
+
+Sample a re-emitted wavelength from the WLS emission spectrum ICDF
+for the material at mat_idx, using uniform random u0 in [0,1).
+
+Returns 0.f if material has no WLS (should not happen in practice
+as callers check has_wls first).
+
+**/
+
+inline QWLS_METHOD float qwls::wavelength(unsigned mat_idx, const float &u0) const
+{
+ int base_row = material_map[mat_idx];
+ if (base_row < 0)
+ return 0.f;
+ return wavelength_at_row(unsigned(base_row), u0);
+}
+
+/**
+qwls::wavelength_at_row
+--------------------------
+
+ICDF texture lookup with HD (high-definition) support.
+base_row is the first of 3 rows for this WLS material:
+ row 0: standard resolution (full CDF range)
+ row 1: LHS high-res (0.00 -> 0.05 for hd_factor=20)
+ row 2: RHS high-res (0.95 -> 1.00 for hd_factor=20)
+
+Uses normalized texture coordinates with linear interpolation,
+matching the qscint.h implementation.
+
+**/
+
+inline QWLS_METHOD float qwls::wavelength_at_row(unsigned base_row, const float &u0) const
+{
+ float y0 = (float(base_row) + 0.5f) / float(tex_height);
+ float y1 = (float(base_row + 1) + 0.5f) / float(tex_height);
+ float y2 = (float(base_row + 2) + 0.5f) / float(tex_height);
+
+ float wl;
+
+ if (hd_factor == 0)
+ {
+ wl = tex2D(wls_tex, u0, y0);
+ }
+ else if (hd_factor == 10)
+ {
+ if (u0 < 0.1f)
+ wl = tex2D(wls_tex, u0 * 10.f, y1);
+ else if (u0 > 0.9f)
+ wl = tex2D(wls_tex, (u0 - 0.9f) * 10.f, y2);
+ else
+ wl = tex2D(wls_tex, u0, y0);
+ }
+ else // hd_factor == 20
+ {
+ if (u0 < 0.05f)
+ wl = tex2D(wls_tex, u0 * 20.f, y1);
+ else if (u0 > 0.95f)
+ wl = tex2D(wls_tex, (u0 - 0.95f) * 20.f, y2);
+ else
+ wl = tex2D(wls_tex, u0, y0);
+ }
+
+ return wl;
+}
+
+#endif
diff --git a/sysrap/snam.h b/sysrap/snam.h
index 1f2e49f57..4dd713e9f 100644
--- a/sysrap/snam.h
+++ b/sysrap/snam.h
@@ -16,6 +16,10 @@ struct snam
static constexpr const char* OPTICAL = "optical.npy" ;
static constexpr const char* ICDF = "icdf.npy" ;
+ static constexpr const char *WLS_ICDF = "wls_icdf.npy";
+ static constexpr const char *WLS_MAT_MAP = "wls_mat_map.npy";
+ static constexpr const char *WLS_TIME_CONSTANTS = "wls_time_constants.npy";
+
static constexpr const char* MULTIFILM = "multifilm.npy" ;
static constexpr const char* PROPCOM = "propcom.npy" ;
diff --git a/sysrap/sproplist.h b/sysrap/sproplist.h
index deeffa02a..a6421ef8c 100644
--- a/sysrap/sproplist.h
+++ b/sysrap/sproplist.h
@@ -29,16 +29,16 @@ not to define ABSLENGTH and RAYLEIGH properties.
struct sproplist
{
- static constexpr const char* MATERIAL = R"(
+ static constexpr const char *MATERIAL = R"(
0 0 RINDEX 1
0 1 ABSLENGTH 1e12
0 2 RAYLEIGH 1e12
0 3 REEMISSIONPROB 0.
1 0 GROUPVEL 299.792458
1 1 SPARE11 0.
- 1 2 SPARE12 0.
+ 1 2 WLSABSLENGTH 1e12
1 3 SPARE13 0.
- )" ;
+ )";
// default GROUPVEL set to c_light_mm_per_ns, see U4PhysicalConstants.h
static constexpr const char* SURFACE = R"(
diff --git a/sysrap/sstandard.h b/sysrap/sstandard.h
index 1c62d6f0c..cefe1cc51 100644
--- a/sysrap/sstandard.h
+++ b/sysrap/sstandard.h
@@ -105,6 +105,9 @@ struct sstandard
const NP* icdf ;
+ const NP *wls_icdf;
+ const NP *wls_mat_map;
+ const NP *wls_time_constants;
sstandard();
@@ -146,9 +149,7 @@ struct sstandard
static NP* unused_create(const sproplist* pl, const std::vector& names, const NPFold* fold );
};
-
-inline sstandard::sstandard()
- :
+inline sstandard::sstandard() :
dom(nullptr),
wavelength(nullptr),
energy(nullptr),
@@ -158,7 +159,10 @@ inline sstandard::sstandard()
bd(nullptr),
bnd(nullptr),
optical(nullptr),
- icdf(nullptr)
+ icdf(nullptr),
+ wls_icdf(nullptr),
+ wls_mat_map(nullptr),
+ wls_time_constants(nullptr)
{
}
@@ -211,6 +215,10 @@ inline NPFold* sstandard::serialize() const
fold->add(snam::ICDF, icdf) ;
+ fold->add(snam::WLS_ICDF, wls_icdf);
+ fold->add(snam::WLS_MAT_MAP, wls_mat_map);
+ fold->add(snam::WLS_TIME_CONSTANTS, wls_time_constants);
+
return fold ;
}
@@ -228,6 +236,10 @@ inline void sstandard::import(const NPFold* fold )
optical = fold->get(snam::OPTICAL);
icdf = fold->get(snam::ICDF);
+
+ wls_icdf = fold->get(snam::WLS_ICDF);
+ wls_mat_map = fold->get(snam::WLS_MAT_MAP);
+ wls_time_constants = fold->get(snam::WLS_TIME_CONSTANTS);
}
inline void sstandard::save(const char* base, const char* rel )
diff --git a/sysrap/sstate.h b/sysrap/sstate.h
index 5795fcd70..b74e7c75b 100644
--- a/sysrap/sstate.h
+++ b/sysrap/sstate.h
@@ -3,26 +3,26 @@
#if defined(__CUDACC__) || defined(__CUDABE__)
#define SSTATE_METHOD __device__
#else
- #define SSTATE_METHOD
-#endif
+ #define SSTATE_METHOD
+#endif
/**
sstate.h
=========
-This was formerly qstate.h but as no CUDA specifics it
-belongs down in sysrap not up in QUDARap.
+This was formerly qstate.h but as no CUDA specifics it
+belongs down in sysrap not up in QUDARap.
-Populated by qsim::fill_state from texture and buffer lookups
-using photon wavelength and the boundary obtained from geometry
-intersect.
+Populated by qsim::fill_state from texture and buffer lookups
+using photon wavelength and the boundary obtained from geometry
+intersect.
**/
struct sstate
{
float4 material1 ; // refractive_index/absorption_length/scattering_length/reemission_prob
- float4 m1group2 ; // material1_group_velocity/material2_group_velocity/spare2/spare3
+ float4 m1group2; // material1_group_velocity/material2_group_velocity/wls_absorption_length/spare3
float4 material2 ;
float4 surface ; // detect/absorb/reflect_specular/reflect_diffuse
@@ -65,27 +65,16 @@ inline void sstate::save(const char* dir) const
#else
inline std::ostream& operator<<(std::ostream& os, const sstate& s )
{
- os << "sstate"
- << std::endl
- << " material1 " << s.material1
- << " (refractive_index/absorption_length/scattering_length/reemission_prob) "
- << std::endl
+ os << "sstate" << std::endl
+ << " material1 " << s.material1 << " (refractive_index/absorption_length/scattering_length/reemission_prob) "
+ << std::endl
<< " m1group2 " << s.m1group2
- << " (material1_group_velocity/material2_group_velocity/spare2/spare3) "
- << std::endl
- << " material2 " << s.material2
- << " (refractive_index/absorption_length/scattering_length/reemission_prob) "
- << std::endl
- << " surface " << s.surface
- << " (detect/absorb/reflect_specular/reflect_diffuse) "
- << std::endl
- << " optical " << s.optical
- << " (x/y/z/w index/type/finish/value) "
- << std::endl
- << " index " << s.index
- << " (indices of m1/m2/surf/sensor) "
- << std::endl
- ;
+ << " (material1_group_velocity/material2_group_velocity/wls_absorption_length/spare3) " << std::endl
+ << " material2 " << s.material2 << " (refractive_index/absorption_length/scattering_length/reemission_prob) "
+ << std::endl
+ << " surface " << s.surface << " (detect/absorb/reflect_specular/reflect_diffuse) " << std::endl
+ << " optical " << s.optical << " (x/y/z/w index/type/finish/value) " << std::endl
+ << " index " << s.index << " (indices of m1/m2/surf/sensor) " << std::endl;
return os;
}
#endif
diff --git a/tests/geom/wls_test.gdml b/tests/geom/wls_test.gdml
new file mode 100644
index 000000000..984f99bf4
--- /dev/null
+++ b/tests/geom/wls_test.gdml
@@ -0,0 +1,144 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/u4/CMakeLists.txt b/u4/CMakeLists.txt
index dba03458d..8290c5942 100644
--- a/u4/CMakeLists.txt
+++ b/u4/CMakeLists.txt
@@ -61,6 +61,7 @@ set(HEADERS
U4Material.hh
U4Mat.h
U4Scint.h
+ U4WLS.h
U4Volume.h
U4Surface.h
diff --git a/u4/U4Tree.h b/u4/U4Tree.h
index c7495afad..c9222900a 100644
--- a/u4/U4Tree.h
+++ b/u4/U4Tree.h
@@ -82,6 +82,7 @@ controlled via envvar::
#include "U4Mesh.h"
#include "U4Scint.h"
+#include "U4WLS.h"
#include "U4Solid.h"
#include "U4PhysicsTable.h"
@@ -112,6 +113,7 @@ struct U4Tree
std::vector solids ;
U4PhysicsTable* rayleigh_table ;
U4Scint* scint ;
+ U4WLS *wls;
// disable the below with settings with by defining the below envvar
static constexpr const char* __DISABLE_OSUR_IMPLICIT = "U4Tree__DISABLE_OSUR_IMPLICIT" ;
@@ -152,6 +154,7 @@ struct U4Tree
void initMaterial(const G4Material* const mt);
void initScint();
+ void initWLS();
void initSurfaces();
void initSolids();
@@ -249,12 +252,7 @@ inline U4Tree* U4Tree::Create(
return tree ;
}
-inline U4Tree::U4Tree(
- stree* st_,
- const G4VPhysicalVolume* const top_,
- U4SensorIdentifier* sid_
- )
- :
+inline U4Tree::U4Tree(stree *st_, const G4VPhysicalVolume *const top_, U4SensorIdentifier *sid_) :
st(st_),
top(top_),
sid(sid_ ? sid_ : new U4SensorIdentifierDefault),
@@ -262,10 +260,11 @@ inline U4Tree::U4Tree(
num_surface_standard(-1),
rayleigh_table(CreateRayleighTable()),
scint(nullptr),
+ wls(nullptr),
enable_osur(!ssys::getenvbool(__DISABLE_OSUR_IMPLICIT)),
enable_isur(!ssys::getenvbool(__DISABLE_ISUR_IMPLICIT)),
- material_debug(ssys::getenvint(__MATERIAL_DEBUG,0)),
- solid_debug(ssys::getenvint(__SOLID_DEBUG,0))
+ material_debug(ssys::getenvint(__MATERIAL_DEBUG, 0)),
+ solid_debug(ssys::getenvint(__SOLID_DEBUG, 0))
{
init();
}
@@ -292,6 +291,9 @@ inline void U4Tree::init()
LOG(LEVEL) << "-initScint" ;
initScint();
+ LOG(LEVEL) << "-initWLS";
+ initWLS();
+
LOG(LEVEL) << "-initSurfaces" ;
initSurfaces();
@@ -384,6 +386,28 @@ inline void U4Tree::initScint()
}
}
+/**
+U4Tree::initWLS
+------------------
+
+Scans all G4 materials for WLS properties (WLSCOMPONENT, WLSTIMECONSTANT).
+Creates inverse CDF texture data and material mapping for GPU-side WLS
+wavelength sampling. Stored in st->standard for serialization and upload.
+
+**/
+
+inline void U4Tree::initWLS()
+{
+ wls = U4WLS::Create(materials);
+ if (wls)
+ {
+ st->standard->wls_icdf = wls->icdf;
+ st->standard->wls_mat_map = wls->mat_map;
+ st->standard->wls_time_constants = wls->time_constants;
+ LOG(LEVEL) << wls->desc();
+ }
+}
+
/**
U4Tree::CreateRayleighTable
----------------------------
diff --git a/u4/U4WLS.h b/u4/U4WLS.h
new file mode 100644
index 000000000..409c2ff57
--- /dev/null
+++ b/u4/U4WLS.h
@@ -0,0 +1,213 @@
+#pragma once
+/**
+U4WLS.h : Wavelength Shifting ICDF Creation
+===============================================
+
+Creates inverse CDF textures for wavelength shifting (WLS) materials,
+analogous to U4Scint.h for scintillation. Supports multiple WLS materials
+by stacking ICDF rows into a single texture.
+
+For each material with a WLSCOMPONENT property:
+1. Integrates the emission spectrum to get a CDF
+2. Inverts it at 4096 uniformly-spaced CDF values (3 resolutions for HD)
+3. Extracts WLSTIMECONSTANT from the material properties table
+
+The output arrays:
+- icdf: shape (num_wls_mat*3, 4096, 1) — stacked HD ICDF rows
+- mat_map: shape (num_total_mat,) int — maps material index to WLS row (-1 = no WLS)
+- time_constants: shape (num_wls_mat,) float — per-WLS-material time constant
+
+The G4 WLS process (G4OpWLS) uses these material properties:
+- WLSABSLENGTH: absorption length as f(energy) — handled via boundary texture
+- WLSCOMPONENT: emission spectrum as f(energy) — converted to ICDF here
+- WLSTIMECONSTANT: re-emission time delay (scalar) — extracted here
+
+**/
+
+#include
+#include
+#include
+#include
+#include
+
+#include "G4Material.hh"
+#include "G4MaterialPropertiesTable.hh"
+#include "G4MaterialPropertyVector.hh"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+
+#include "NP.hh"
+#include "NPFold.h"
+#include "SLOG.hh"
+#include "U4MaterialPropertyVector.h"
+#include "U4Scint.h" // reuse Integral and CreateGeant4InterpolatedInverseCDF
+
+struct U4WLS
+{
+ static constexpr const char *WLSCOMPONENT_KEY = "WLSCOMPONENT";
+ static constexpr const char *WLSTIMECONSTANT_KEY = "WLSTIMECONSTANT";
+
+ static U4WLS *Create(const std::vector &mats);
+
+ const NP *icdf; // (num_wls*3, 4096, 1) stacked HD ICDF for all WLS materials
+ const NP *mat_map; // (num_total_mat,) int: material idx -> base ICDF row, or -1
+ const NP *time_constants; // (num_wls,) float: time constant per WLS material
+
+ unsigned num_wls;
+ unsigned num_mat;
+
+ U4WLS(const std::vector &mats, const std::vector &wls_indices,
+ const std::vector &wls_components,
+ const std::vector &wls_time_consts);
+
+ std::string desc() const;
+};
+
+/**
+U4WLS::Create
+---------------
+
+Scans all materials for WLSCOMPONENT property. For each material
+that has it, extracts the emission spectrum and time constant.
+
+Returns nullptr if no WLS materials are found.
+
+**/
+
+inline U4WLS *U4WLS::Create(const std::vector &mats)
+{
+ std::vector wls_indices;
+ std::vector wls_components;
+ std::vector wls_time_consts;
+
+ for (unsigned i = 0; i < mats.size(); i++)
+ {
+ const G4Material *mat = mats[i];
+ G4MaterialPropertiesTable *mpt = mat->GetMaterialPropertiesTable();
+ if (mpt == nullptr)
+ continue;
+
+ G4MaterialPropertyVector *wlscomp = mpt->GetProperty(WLSCOMPONENT_KEY);
+ if (wlscomp == nullptr)
+ continue;
+
+ // Found a WLS material
+ wls_indices.push_back(i);
+ wls_components.push_back(wlscomp);
+
+ // Extract time constant (scalar property, default 0 = instant re-emission)
+ double tc = 0.0;
+ if (mpt->ConstPropertyExists(WLSTIMECONSTANT_KEY))
+ {
+ tc = mpt->GetConstProperty(WLSTIMECONSTANT_KEY) / ns; // convert to ns
+ }
+ wls_time_consts.push_back(tc);
+ }
+
+ if (wls_indices.empty())
+ return nullptr;
+
+ return new U4WLS(mats, wls_indices, wls_components, wls_time_consts);
+}
+
+/**
+U4WLS::U4WLS
+--------------
+
+Builds the ICDF texture data and material mapping arrays.
+
+For each WLS material:
+1. Integrate WLSCOMPONENT to get CDF (reuses U4Scint::Integral)
+2. Build 3-layer HD ICDF (reuses U4Scint::CreateGeant4InterpolatedInverseCDF)
+3. Stack into combined ICDF array
+
+**/
+
+inline U4WLS::U4WLS(const std::vector &mats, const std::vector &wls_indices,
+ const std::vector &wls_components,
+ const std::vector &wls_time_consts) :
+ icdf(nullptr),
+ mat_map(nullptr),
+ time_constants(nullptr),
+ num_wls(wls_indices.size()),
+ num_mat(mats.size())
+{
+ assert(num_wls > 0);
+ assert(wls_components.size() == num_wls);
+ assert(wls_time_consts.size() == num_wls);
+
+ int num_bins = 4096;
+ int hd_factor = 20;
+
+ // Build per-material ICDFs and stack them
+ std::vector icdfs;
+ for (unsigned w = 0; w < num_wls; w++)
+ {
+ const G4MaterialPropertyVector *comp = wls_components[w];
+ const G4Material *mat = mats[wls_indices[w]];
+ const char *matname = mat->GetName().c_str();
+
+ // Integrate emission spectrum to get CDF
+ G4MaterialPropertyVector *integral = U4Scint::Integral(comp);
+
+ // Build 3-layer HD ICDF (wavelength values in nm)
+ NP *one_icdf = U4Scint::CreateGeant4InterpolatedInverseCDF(integral, num_bins, hd_factor, matname,
+ false /*energy_not_wavelength*/
+ );
+
+ assert(one_icdf);
+ assert(one_icdf->has_shape(3, num_bins, 1));
+ icdfs.push_back(one_icdf);
+ }
+
+ // Stack all ICDFs into a single array: (num_wls*3, 4096, 1)
+ {
+ NP *stacked = NP::Make(num_wls * 3, num_bins, 1);
+ double *dst = stacked->values();
+ for (unsigned w = 0; w < num_wls; w++)
+ {
+ const double *src = icdfs[w]->cvalues();
+ unsigned row_size = 3 * num_bins * 1;
+ memcpy(dst + w * row_size, src, row_size * sizeof(double));
+ }
+ stacked->set_meta("hd_factor", hd_factor);
+ stacked->set_meta("num_bins", num_bins);
+ stacked->set_meta("num_wls", num_wls);
+ icdf = stacked;
+ }
+
+ // Build material index -> ICDF row mapping
+ // For material i, mat_map[i] = base row in ICDF texture (0, 3, 6, ...)
+ // or -1 if material has no WLS
+ {
+ NP *mm = NP::Make(num_mat);
+ int *mm_v = mm->values();
+ for (unsigned i = 0; i < num_mat; i++)
+ mm_v[i] = -1;
+ for (unsigned w = 0; w < num_wls; w++)
+ {
+ mm_v[wls_indices[w]] = w * 3; // base row for this WLS material's 3 HD layers
+ }
+ mat_map = mm;
+ }
+
+ // Build time constants array (in ns)
+ {
+ NP *tc = NP::Make(num_wls);
+ float *tc_v = tc->values();
+ for (unsigned w = 0; w < num_wls; w++)
+ {
+ tc_v[w] = float(wls_time_consts[w]);
+ }
+ time_constants = tc;
+ }
+}
+
+inline std::string U4WLS::desc() const
+{
+ std::stringstream ss;
+ ss << "U4WLS::desc" << " num_wls " << num_wls << " num_mat " << num_mat << " icdf " << (icdf ? icdf->sstr() : "-")
+ << " mat_map " << (mat_map ? mat_map->sstr() : "-") << " time_constants "
+ << (time_constants ? time_constants->sstr() : "-");
+ return ss.str();
+}