feat(qudarap): add GPU wavelength shifting (WLS) physics

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
+<define>
+    <matrix coldim="2" name="WLSABSLENGTH" values="1.77e-06 10000.0 ... 4.13e-06 0.01"/>
+    <matrix coldim="2" name="WLSCOMPONENT" values="1.77e-06 0.00 ... 3.10e-06 0.00"/>
+    <matrix coldim="1" name="WLSTIMECONSTANT" values="0.5"/>
+</define>
+<materials>
+    <material name="WLSMaterial">
+        <property name="WLSABSLENGTH" ref="WLSABSLENGTH"/>
+        <property name="WLSCOMPONENT" ref="WLSCOMPONENT"/>
+        <property name="WLSTIMECONSTANT" ref="WLSTIMECONSTANT"/>
+    </material>
+</materials>
+```
+
+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 | ✓ | ✓ | ✓ | ✗ | ✗ |

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
+  }
+}

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

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 ;
 
 

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 ;

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<qbnd>(const qbnd* array, unsigned num_it
 template sevent*        QU::UploadArray<sevent>(const sevent* array, unsigned num_items, const char* label) ;
 template qdebug*        QU::UploadArray<qdebug>(const qdebug* array, unsigned num_items, const char* label) ;
 template qscint*        QU::UploadArray<qscint>(const qscint* array, unsigned num_items, const char* label) ;
+template qwls *QU::UploadArray<qwls>(const qwls *array, unsigned num_items, const char *label);
 template qcerenkov*     QU::UploadArray<qcerenkov>(const qcerenkov* array, unsigned num_items, const char* label) ;
 template qbase*         QU::UploadArray<qbase>(const qbase* array, unsigned num_items, const char* label) ;
 

qudarap/QWls.cc

@@ -0,0 +1,132 @@
+#include <cassert>
+#include <csignal>
+#include <sstream>
+
+#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<qwls>(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<float> *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<float> *tx = new QTex<float>(nx, ny, src->cvalues<float>(), '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<float> *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<void **>(&d_mat_map), mat_map_size));
+    QUDA_CHECK(cudaMemcpy(d_mat_map, mat_map->cvalues<int>(), 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<void **>(&d_tc), tc_size));
+    QUDA_CHECK(cudaMemcpy(d_tc, time_constants->cvalues<float>(), 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();
+}