eic-opticks/u4/U4WLS.h

0001 #pragma once
0002 /**
0003 U4WLS.h : Wavelength Shifting ICDF Creation
0004 ===============================================
0005
0006 Creates inverse CDF textures for wavelength shifting (WLS) materials,
0007 analogous to U4Scint.h for scintillation. Supports multiple WLS materials
0008 by stacking ICDF rows into a single texture.
0009
0010 For each material with a WLSCOMPONENT property:
0011 1. Integrates the emission spectrum to get a CDF
0012 2. Inverts it at 4096 uniformly-spaced CDF values (3 resolutions for HD)
0013 3. Extracts WLSTIMECONSTANT from the material properties table
0014
0015 The output arrays:
0016 - icdf: shape (num_wls_mat*3, 4096, 1) — stacked HD ICDF rows
0017 - mat_map: shape (num_total_mat,) int — maps material index to WLS row (-1 = no WLS)
0018 - time_constants: shape (num_wls_mat,) float — per-WLS-material time constant
0019
0020 The G4 WLS process (G4OpWLS) uses these material properties:
0021 - WLSABSLENGTH: absorption length as f(energy) — handled via boundary texture
0022 - WLSCOMPONENT: emission spectrum as f(energy) — converted to ICDF here
0023 - WLSTIMECONSTANT: re-emission time delay (scalar) — extracted here
0024
0025 **/
0026
0027 #include <cassert>
0028 #include <iomanip>
0029 #include <sstream>
0030 #include <string>
0031 #include <vector>
0032
0033 #include "G4Material.hh"
0034 #include "G4MaterialPropertiesTable.hh"
0035 #include "G4MaterialPropertyVector.hh"
0036 #include "G4PhysicalConstants.hh"
0037 #include "G4SystemOfUnits.hh"
0038
0039 #include "NP.hh"
0040 #include "NPFold.h"
0041 #include "SLOG.hh"
0042 #include "U4MaterialPropertyVector.h"
0043 #include "U4Scint.h" // reuse Integral and CreateGeant4InterpolatedInverseCDF
0044
0045 struct U4WLS
0046 {
0047     static constexpr const char *WLSCOMPONENT_KEY = "WLSCOMPONENT";
0048     static constexpr const char *WLSTIMECONSTANT_KEY = "WLSTIMECONSTANT";
0049
0050     static U4WLS *Create(const std::vector<const G4Material *> &mats);
0051
0052     const NP *icdf;           // (num_wls*3, 4096, 1) stacked HD ICDF for all WLS materials
0053     const NP *mat_map;        // (num_total_mat,) int: material idx -> base ICDF row, or -1
0054     const NP *time_constants; // (num_wls,) float: time constant per WLS material
0055
0056     unsigned num_wls;
0057     unsigned num_mat;
0058
0059     U4WLS(const std::vector<const G4Material *> &mats, const std::vector<int> &wls_indices,
0060           const std::vector<const G4MaterialPropertyVector *> &wls_components,
0061           const std::vector<double> &wls_time_consts);
0062
0063     std::string desc() const;
0064 };
0065
0066 /**
0067 U4WLS::Create
0068 ---------------
0069
0070 Scans all materials for WLSCOMPONENT property. For each material
0071 that has it, extracts the emission spectrum and time constant.
0072
0073 Returns nullptr if no WLS materials are found.
0074
0075 **/
0076
0077 inline U4WLS *U4WLS::Create(const std::vector<const G4Material *> &mats)
0078 {
0079     std::vector<int> wls_indices;
0080     std::vector<const G4MaterialPropertyVector *> wls_components;
0081     std::vector<double> wls_time_consts;
0082
0083     for (unsigned i = 0; i < mats.size(); i++)
0084     {
0085         const G4Material *mat = mats[i];
0086         G4MaterialPropertiesTable *mpt = mat->GetMaterialPropertiesTable();
0087         if (mpt == nullptr)
0088             continue;
0089
0090         G4MaterialPropertyVector *wlscomp = mpt->GetProperty(WLSCOMPONENT_KEY);
0091         if (wlscomp == nullptr)
0092             continue;
0093
0094         // Found a WLS material
0095         wls_indices.push_back(i);
0096         wls_components.push_back(wlscomp);
0097
0098         // Extract time constant (scalar property, default 0 = instant re-emission)
0099         double tc = 0.0;
0100         if (mpt->ConstPropertyExists(WLSTIMECONSTANT_KEY))
0101         {
0102             tc = mpt->GetConstProperty(WLSTIMECONSTANT_KEY) / ns; // convert to ns
0103         }
0104         wls_time_consts.push_back(tc);
0105     }
0106
0107     if (wls_indices.empty())
0108         return nullptr;
0109
0110     return new U4WLS(mats, wls_indices, wls_components, wls_time_consts);
0111 }
0112
0113 /**
0114 U4WLS::U4WLS
0115 --------------
0116
0117 Builds the ICDF texture data and material mapping arrays.
0118
0119 For each WLS material:
0120 1. Integrate WLSCOMPONENT to get CDF (reuses U4Scint::Integral)
0121 2. Build 3-layer HD ICDF (reuses U4Scint::CreateGeant4InterpolatedInverseCDF)
0122 3. Stack into combined ICDF array
0123
0124 **/
0125
0126 inline U4WLS::U4WLS(const std::vector<const G4Material *> &mats, const std::vector<int> &wls_indices,
0127                     const std::vector<const G4MaterialPropertyVector *> &wls_components,
0128                     const std::vector<double> &wls_time_consts) :
0129     icdf(nullptr),
0130     mat_map(nullptr),
0131     time_constants(nullptr),
0132     num_wls(wls_indices.size()),
0133     num_mat(mats.size())
0134 {
0135     assert(num_wls > 0);
0136     assert(wls_components.size() == num_wls);
0137     assert(wls_time_consts.size() == num_wls);
0138
0139     int num_bins = 4096;
0140     int hd_factor = 20;
0141
0142     // Build per-material ICDFs and stack them
0143     std::vector<const NP *> icdfs;
0144     for (unsigned w = 0; w < num_wls; w++)
0145     {
0146         const G4MaterialPropertyVector *comp = wls_components[w];
0147         const G4Material *mat = mats[wls_indices[w]];
0148         const char *matname = mat->GetName().c_str();
0149
0150         // Integrate emission spectrum to get CDF
0151         G4MaterialPropertyVector *integral = U4Scint::Integral(comp);
0152
0153         // Build 3-layer HD ICDF (wavelength values in nm)
0154         NP *one_icdf = U4Scint::CreateGeant4InterpolatedInverseCDF(integral, num_bins, hd_factor, matname,
0155                                                                    false /*energy_not_wavelength*/
0156         );
0157
0158         assert(one_icdf);
0159         assert(one_icdf->has_shape(3, num_bins, 1));
0160         icdfs.push_back(one_icdf);
0161     }
0162
0163     // Stack all ICDFs into a single array: (num_wls*3, 4096, 1)
0164     {
0165         NP *stacked = NP::Make<double>(num_wls * 3, num_bins, 1);
0166         double *dst = stacked->values<double>();
0167         for (unsigned w = 0; w < num_wls; w++)
0168         {
0169             const double *src = icdfs[w]->cvalues<double>();
0170             unsigned row_size = 3 * num_bins * 1;
0171             memcpy(dst + w * row_size, src, row_size * sizeof(double));
0172         }
0173         stacked->set_meta<int>("hd_factor", hd_factor);
0174         stacked->set_meta<int>("num_bins", num_bins);
0175         stacked->set_meta<int>("num_wls", num_wls);
0176         icdf = stacked;
0177     }
0178
0179     // Build material index -> ICDF row mapping
0180     // For material i, mat_map[i] = base row in ICDF texture (0, 3, 6, ...)
0181     // or -1 if material has no WLS
0182     {
0183         NP *mm = NP::Make<int>(num_mat);
0184         int *mm_v = mm->values<int>();
0185         for (unsigned i = 0; i < num_mat; i++)
0186             mm_v[i] = -1;
0187         for (unsigned w = 0; w < num_wls; w++)
0188         {
0189             mm_v[wls_indices[w]] = w * 3; // base row for this WLS material's 3 HD layers
0190         }
0191         mat_map = mm;
0192     }
0193
0194     // Build time constants array (in ns)
0195     {
0196         NP *tc = NP::Make<float>(num_wls);
0197         float *tc_v = tc->values<float>();
0198         for (unsigned w = 0; w < num_wls; w++)
0199         {
0200             tc_v[w] = float(wls_time_consts[w]);
0201         }
0202         time_constants = tc;
0203     }
0204 }
0205
0206 inline std::string U4WLS::desc() const
0207 {
0208     std::stringstream ss;
0209     ss << "U4WLS::desc" << " num_wls " << num_wls << " num_mat " << num_mat << " icdf " << (icdf ? icdf->sstr() : "-")
0210        << " mat_map " << (mat_map ? mat_map->sstr() : "-") << " time_constants "
0211        << (time_constants ? time_constants->sstr() : "-");
0212     return ss.str();
0213 }