Skip to content

Commit 281e8da

Browse files
ggalgocziplexoos
andauthored
feat(qudarap): add GPU wavelength shifting (WLS) physics (#269)
Add wavelength shifting (WLS) physics to the GPU optical photon simulation. Implements WLS absorption and re-emission using material properties (WLSABSLENGTH, WLSCOMPONENT, WLSTIMECONSTANT) with ICDF-based wavelength sampling via texture lookup. Adds U4WLS.h for extracting WLS properties from Geant4 materials during geometry conversion. Includes exponential time delay for re-emitted photons matching G4OpWLS behavior. Validated on wls_test.gdml geometry with 990/1000 photon detection rate. --------- Co-authored-by: Dmitri Smirnov <dmixsmi@gmail.com>
1 parent dbfb279 commit 281e8da

18 files changed

Lines changed: 974 additions & 80 deletions

File tree

README.md

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -155,6 +155,7 @@ EIC-Opticks provides several examples demonstrating GPU-accelerated optical phot
155155
| `GPUPhotonSource` | Optical photons (torch) | Any GDML | G4 + GPU side-by-side validation |
156156
| `GPUPhotonSourceMinimal` | Optical photons (torch) | Any GDML | GPU-only test |
157157
| `GPUPhotonFileSource` | Optical photons (text file) | Any GDML | GPU-only, user-defined photons from file |
158+
| WLS test | Wavelength shifting | WLS sphere + detector shell | Validate GPU WLS physics |
158159

159160
### Example 1: GPUCerenkov (Cerenkov Only)
160161

@@ -297,6 +298,55 @@ GPUPhotonFileSource -g tests/geom/opticks_raindrop.gdml -p my_photons.txt -m run
297298

298299
**Source files:** `src/GPUPhotonFileSource.cpp`, `src/GPUPhotonFileSource.h`
299300

301+
### Example 6: Wavelength Shifting (WLS) Test
302+
303+
This test validates the GPU wavelength shifting implementation using a dedicated
304+
geometry with a WLS sphere surrounded by a detector shell:
305+
306+
```
307+
Geometry: wls_test.gdml
308+
├── Air world (r=200 mm)
309+
│ ├── WLS sphere (r=20 mm) ← Absorbs UV, re-emits visible
310+
│ └── Glass detector shell (r=28-30 mm) ← 100% detection efficiency
311+
```
312+
313+
The WLS material absorbs UV photons (350 nm) and re-emits them isotropically at
314+
longer wavelengths (peak ~481 nm) with a 0.5 ns exponential time delay. The test
315+
fires 1000 monochromatic 350 nm photons from the origin into the WLS sphere.
316+
317+
```bash
318+
GPUPhotonSourceMinimal -g tests/geom/wls_test.gdml -c wls_test -m tests/run.mac -s 42
319+
```
320+
321+
**Expected results:**
322+
- ~990/1000 photons detected (10 absorbed after failing energy conservation)
323+
- All hits wavelength-shifted from 350 nm to mean ~487 nm
324+
- Energy conservation: no hits with wavelength < 350 nm
325+
- Isotropic re-emission: mean momentum direction near zero
326+
- Time delay: mean ~0.6 ns (propagation + 0.5 ns exponential WLS decay)
327+
328+
**GDML WLS properties required** (same syntax for G4 10.x and 11.x):
329+
```xml
330+
<define>
331+
<matrix coldim="2" name="WLSABSLENGTH" values="1.77e-06 10000.0 ... 4.13e-06 0.01"/>
332+
<matrix coldim="2" name="WLSCOMPONENT" values="1.77e-06 0.00 ... 3.10e-06 0.00"/>
333+
<matrix coldim="1" name="WLSTIMECONSTANT" values="0.5"/>
334+
</define>
335+
<materials>
336+
<material name="WLSMaterial">
337+
<property name="WLSABSLENGTH" ref="WLSABSLENGTH"/>
338+
<property name="WLSCOMPONENT" ref="WLSCOMPONENT"/>
339+
<property name="WLSTIMECONSTANT" ref="WLSTIMECONSTANT"/>
340+
</material>
341+
</materials>
342+
```
343+
344+
Unlike scintillation properties, WLS property names are the same in both Geant4
345+
10.x and 11.x — no dual-naming is needed.
346+
347+
**Test files:** `tests/geom/wls_test.gdml`, `config/wls_test.json`
348+
**Implementation docs:** `docs/WLS_IMPLEMENTATION.md`
349+
300350
### Torch configuration
301351

302352
`GPUPhotonSource` and `GPUPhotonSourceMinimal` read photon source parameters from a
@@ -317,6 +367,7 @@ JSON config file (default `config/dev.json`). Key fields:
317367
|---------|-------------|-------------|-----------------|----------------------|---------------------|
318368
| Cerenkov genstep collection ||||||
319369
| Scintillation genstep collection ||||||
370+
| Wavelength shifting (WLS) ||||||
320371
| Torch photon generation ||||||
321372
| Photon input from text file ||||||
322373
| G4 optical photon tracking ||||||

config/wls_test.json

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
{
2+
"torch": {
3+
"gentype": "TORCH",
4+
"trackid": 0,
5+
"matline": 0,
6+
"numphoton": 1000,
7+
8+
"pos": [0.0, 0.0, 0.0],
9+
"time": 0.0,
10+
11+
"mom": [0.0, 0.0, 1.0],
12+
"weight": 0.0,
13+
14+
"pol": [1.0, 0.0, 0.0],
15+
"wavelength": 350.0,
16+
17+
"zenith": [0.0, 1.0],
18+
"azimuth": [0.0, 1.0],
19+
20+
"radius": 0.0,
21+
"distance": 0.0,
22+
"mode": 255,
23+
"type": "disc"
24+
},
25+
26+
"event": {
27+
"mode": "DebugLite",
28+
"maxslot": 1000000
29+
}
30+
}

qudarap/CMakeLists.txt

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -48,6 +48,8 @@ set(SOURCES
4848
QScint.cc
4949
QScint.cu
5050

51+
QWls.cc
52+
5153
QCerenkovIntegral.cc
5254
QCerenkov.cc
5355
QCerenkov.cu
@@ -117,6 +119,9 @@ SET(HEADERS
117119
QScint.hh
118120
qscint.h
119121

122+
QWls.hh
123+
qwls.h
124+
120125
QCerenkovIntegral.hh
121126
QCerenkov.hh
122127
qcerenkov.h

qudarap/QSim.cc

Lines changed: 33 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -32,19 +32,19 @@
3232
#include "qdebug.h"
3333

3434
#include "QBase.hh"
35-
#include "QEvt.hh"
36-
#include "QRng.hh"
37-
#include "QTex.hh"
38-
#include "QScint.hh"
39-
#include "QCerenkov.hh"
4035
#include "QBnd.hh"
41-
#include "QProp.hh"
42-
#include "QMultiFilm.hh"
36+
#include "QCerenkov.hh"
37+
#include "QDebug.hh"
4338
#include "QEvt.hh"
39+
#include "QMultiFilm.hh"
4440
#include "QOptical.hh"
45-
#include "QSimLaunch.hh"
46-
#include "QDebug.hh"
4741
#include "QPMT.hh"
42+
#include "QProp.hh"
43+
#include "QRng.hh"
44+
#include "QScint.hh"
45+
#include "QSimLaunch.hh"
46+
#include "QTex.hh"
47+
#include "QWls.hh"
4848

4949
#include "QSim.hh"
5050

@@ -179,6 +179,27 @@ void QSim::UploadComponents( const SSim* ssim )
179179
LOG(LEVEL) << scint->desc();
180180
}
181181

182+
const NP *wls_icdf = ssim->get(snam::WLS_ICDF);
183+
const NP *wls_mat_map = ssim->get(snam::WLS_MAT_MAP);
184+
if (wls_icdf == nullptr || wls_mat_map == nullptr)
185+
{
186+
LOG(LEVEL) << " wls_icdf or wls_mat_map null — no WLS materials in geometry ";
187+
}
188+
else
189+
{
190+
const NP *wls_tc = ssim->get(snam::WLS_TIME_CONSTANTS);
191+
if (wls_tc)
192+
{
193+
unsigned hd_factor = 20u;
194+
QWls *qwls_ = new QWls(wls_icdf, wls_mat_map, wls_tc, hd_factor);
195+
LOG(LEVEL) << qwls_->desc();
196+
}
197+
else
198+
{
199+
LOG(error) << " wls_icdf and wls_mat_map present but wls_time_constants missing ";
200+
}
201+
}
202+
182203
// TODO: make this more like the others : acting on the available inputs rather than the mode
183204
bool is_simtrace = SEventConfig::IsRGModeSimtrace() ;
184205
if(is_simtrace == false )
@@ -258,13 +279,13 @@ singleton components.
258279
259280
**/
260281

261-
QSim::QSim()
262-
:
282+
QSim::QSim() :
263283
base(QBase::Get()),
264284
qev(new QEvt),
265285
sev(qev->sev),
266286
rng(QRng::Get()),
267287
scint(QScint::Get()),
288+
qwls(QWls::Get()),
268289
cerenkov(QCerenkov::Get()),
269290
bnd(QBnd::Get()),
270291
debug_(QDebug::Get()),
@@ -314,6 +335,7 @@ void QSim::init()
314335
sim->multifilm = multifilm ? multifilm->d_multifilm : nullptr ;
315336
sim->cerenkov = cerenkov ? cerenkov->d_cerenkov : nullptr ;
316337
sim->scint = scint ? scint->d_scint : nullptr ;
338+
sim->wls = qwls ? qwls->d_wls : nullptr;
317339
sim->pmt = pmt ? pmt->d_pmt : nullptr ;
318340

319341

qudarap/QSim.hh

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ struct QBase ;
3838
struct QEvt ;
3939
struct QRng ;
4040
struct QScint ;
41+
struct QWls;
4142
struct QCerenkov ;
4243
struct QBnd ;
4344
struct QMultiFilm;
@@ -74,6 +75,7 @@ struct QUDARAP_API QSim
7475

7576
const QRng* rng ;
7677
const QScint* scint ;
78+
const QWls *qwls;
7779
const QCerenkov* cerenkov ;
7880
const QBnd* bnd ;
7981
const QOptical* optical ;

qudarap/QU.cc

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -29,15 +29,15 @@
2929
#include "qsim.h"
3030

3131
#include "qbase.h"
32-
#include "qprop.h"
33-
#include "qpmt.h"
34-
#include "qdebug.h"
35-
#include "qscint.h"
3632
#include "qcerenkov.h"
3733
#include "qcurandwrap.h"
38-
#include "scurandref.h"
34+
#include "qdebug.h"
3935
#include "qmultifilm.h"
40-
36+
#include "qpmt.h"
37+
#include "qprop.h"
38+
#include "qscint.h"
39+
#include "qwls.h"
40+
#include "scurandref.h"
4141

4242
const plog::Severity QU::LEVEL = SLOG::EnvLevel("QU", "DEBUG") ;
4343
bool QU::MEMCHECK = ssys::getenvbool(_MEMCHECK);
@@ -171,6 +171,7 @@ template qbnd* QU::UploadArray<qbnd>(const qbnd* array, unsigned num_it
171171
template sevent* QU::UploadArray<sevent>(const sevent* array, unsigned num_items, const char* label) ;
172172
template qdebug* QU::UploadArray<qdebug>(const qdebug* array, unsigned num_items, const char* label) ;
173173
template qscint* QU::UploadArray<qscint>(const qscint* array, unsigned num_items, const char* label) ;
174+
template qwls *QU::UploadArray<qwls>(const qwls *array, unsigned num_items, const char *label);
174175
template qcerenkov* QU::UploadArray<qcerenkov>(const qcerenkov* array, unsigned num_items, const char* label) ;
175176
template qbase* QU::UploadArray<qbase>(const qbase* array, unsigned num_items, const char* label) ;
176177

qudarap/QWls.cc

Lines changed: 132 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,132 @@
1+
#include <cassert>
2+
#include <csignal>
3+
#include <sstream>
4+
5+
#include "scuda.h"
6+
#include "squad.h"
7+
8+
#include "NP.hh"
9+
#include "SLOG.hh"
10+
#include "ssys.h"
11+
12+
#include "QTex.hh"
13+
#include "QU.hh"
14+
#include "QUDA_CHECK.h"
15+
#include "QWls.hh"
16+
17+
#include "qwls.h"
18+
19+
const plog::Severity QWls::LEVEL = SLOG::EnvLevel("QWls", "DEBUG");
20+
21+
const QWls *QWls::INSTANCE = nullptr;
22+
const QWls *QWls::Get()
23+
{
24+
return INSTANCE;
25+
}
26+
27+
/**
28+
QWls::QWls
29+
------------
30+
31+
1. Narrows ICDF from double to float if needed
32+
2. Uploads ICDF into GPU texture
33+
3. Creates qwls instance with device pointers and uploads it
34+
35+
**/
36+
37+
QWls::QWls(const NP *wls_icdf, const NP *mat_map, const NP *time_constants, unsigned hd_factor) :
38+
dsrc(wls_icdf->ebyte == 8 ? wls_icdf : nullptr),
39+
src(wls_icdf->ebyte == 4 ? wls_icdf : NP::MakeNarrow(dsrc)),
40+
tex(MakeWlsTex(src, hd_factor)),
41+
wls(MakeInstance(tex, mat_map, time_constants, hd_factor, time_constants->shape[0])),
42+
d_wls(QU::UploadArray<qwls>(wls, 1, "QWls::QWls/d_wls"))
43+
{
44+
INSTANCE = this;
45+
}
46+
47+
/**
48+
QWls::MakeWlsTex
49+
-------------------
50+
51+
Creates a 2D CUDA texture from the ICDF array.
52+
Shape: (num_wls*3, 4096, 1) where 3 = HD layers per material.
53+
54+
**/
55+
56+
QTex<float> *QWls::MakeWlsTex(const NP *src, unsigned hd_factor)
57+
{
58+
assert(src);
59+
assert(src->shape.size() == 3);
60+
61+
unsigned ni = src->shape[0]; // height: num_wls * 3
62+
unsigned nj = src->shape[1]; // width: 4096
63+
unsigned nk = src->shape[2]; // 1
64+
65+
assert(nk == 1);
66+
assert(nj == 4096);
67+
assert(ni % 3 == 0); // must be multiple of 3 (3 HD layers per material)
68+
assert(src->uifc == 'f' && src->ebyte == 4);
69+
70+
unsigned ny = ni; // height
71+
unsigned nx = nj; // width
72+
73+
bool normalizedCoords = true;
74+
QTex<float> *tx = new QTex<float>(nx, ny, src->cvalues<float>(), 'L', normalizedCoords, src);
75+
76+
tx->setHDFactor(hd_factor);
77+
tx->uploadMeta();
78+
79+
LOG(LEVEL) << " src " << src->desc() << " nx (width) " << nx << " ny (height) " << ny << " tx.HDFactor "
80+
<< tx->getHDFactor();
81+
82+
return tx;
83+
}
84+
85+
/**
86+
QWls::MakeInstance
87+
---------------------
88+
89+
Creates the host-side qwls struct populated with device pointers.
90+
Uploads material_map and time_constants to device memory.
91+
92+
**/
93+
94+
qwls *QWls::MakeInstance(const QTex<float> *tex, const NP *mat_map, const NP *time_constants, unsigned hd_factor,
95+
unsigned num_wls)
96+
{
97+
assert(mat_map);
98+
assert(time_constants);
99+
assert(mat_map->uifc == 'i' && mat_map->ebyte == 4);
100+
assert(time_constants->uifc == 'f' && time_constants->ebyte == 4);
101+
102+
qwls *w = new qwls;
103+
w->wls_tex = tex->texObj;
104+
w->hd_factor = hd_factor;
105+
w->num_wls = num_wls;
106+
w->tex_height = tex->height;
107+
108+
// Upload material_map to device
109+
unsigned num_mat = mat_map->shape[0];
110+
int *d_mat_map = nullptr;
111+
size_t mat_map_size = num_mat * sizeof(int);
112+
QUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&d_mat_map), mat_map_size));
113+
QUDA_CHECK(cudaMemcpy(d_mat_map, mat_map->cvalues<int>(), mat_map_size, cudaMemcpyHostToDevice));
114+
w->material_map = d_mat_map;
115+
116+
// Upload time_constants to device
117+
float *d_tc = nullptr;
118+
size_t tc_size = num_wls * sizeof(float);
119+
QUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&d_tc), tc_size));
120+
QUDA_CHECK(cudaMemcpy(d_tc, time_constants->cvalues<float>(), tc_size, cudaMemcpyHostToDevice));
121+
w->time_constants = d_tc;
122+
123+
return w;
124+
}
125+
126+
std::string QWls::desc() const
127+
{
128+
std::stringstream ss;
129+
ss << "QWls" << " dsrc " << (dsrc ? dsrc->desc() : "-") << " src " << (src ? src->desc() : "-") << " tex "
130+
<< (tex ? tex->desc() : "-");
131+
return ss.str();
132+
}

0 commit comments

Comments
 (0)