Import WLS data for optical photons (celeritas-project#1165)

* Add and import `ImportWavelengthShifting` data
* Update GeantImporter test harness
* Import wls mean number of photons and time constant
* Expand GeantImporter test
* Update celer-dump-data; add imported data to the RootInterface linkdef

app/celer-dump-data.cc

@@ -702,6 +702,17 @@ void print_optical_material_data(ImportData::ImportOpticalMap const& iom)
         cout << POM_STREAM_VECTOR(mid, abs, absorption_length, IU::len);
     }
     cout << endl;
+    cout << "\n## WLS";
+    cout << header;
+    for (auto const& [mid, val] : iom)
+    {
+        auto const& wls = val.wls;
+        cout << POM_STREAM_SCALAR(mid, wls, mean_num_photons, IU::unitless);
+        cout << POM_STREAM_SCALAR(mid, wls, time_constant, IU::time);
+        cout << POM_STREAM_VECTOR(mid, wls, absorption_length, IU::len);
+        cout << POM_STREAM_VECTOR(mid, wls, component, IU::unitless);
+    }
+    cout << endl;
 #undef PEP_STREAM_SCALAR
 #undef PEP_STREAM_VECTOR
 }

src/celeritas/ext/GeantImporter.cc

@@ -548,6 +548,17 @@ ImportData::ImportOpticalMap import_optical()
                             "ABSLENGTH",
                             ImportUnits::len);
 
+        // Save WLS properties
+        get_property.scalar(&optical.wls.mean_num_photons,
+                            "WLSMEANNUMBERPHOTONS",
+                            ImportUnits::unitless);
+        get_property.scalar(
+            &optical.wls.time_constant, "WLSTIMECONSTANT", ImportUnits::time);
+        get_property.vector(
+            &optical.wls.absorption_length, "WLSABSLENGTH", ImportUnits::len);
+        get_property.vector(
+            &optical.wls.component, "WLSCOMPONENT", ImportUnits::unitless);
+
         if (optical)
         {
             result[mat_idx] = optical;

src/celeritas/ext/RootInterfaceLinkDef.h

@@ -24,6 +24,7 @@
 #pragma link C++ class celeritas::ImportScintData+;
 #pragma link C++ class celeritas::ImportOpticalRayleigh+;
 #pragma link C++ class celeritas::ImportOpticalAbsorption+;
+#pragma link C++ class celeritas::ImportWavelengthShift+;
 #pragma link C++ class celeritas::ImportOpticalProperty+;
 #pragma link C++ class celeritas::ImportOpticalMaterial+;
 #pragma link C++ class celeritas::ImportProductionCut+;

src/celeritas/io/ImportOpticalMaterial.hh

@@ -147,6 +147,32 @@ struct ImportOpticalProperty
     }
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * Store optical photon wavelength shifting properties.
+ *
+ * The component vector represents the relative population as a function of the
+ * re-emission energy. It is used to define an inverse CDF needed to sample the
+ * re-emitted optical photon energy.
+ */
+struct ImportWavelengthShift
+{
+    double mean_num_photons;  //!< Mean number of re-emitted photons
+    double time_constant;  //!< Time delay between absorption and re-emission
+    ImportPhysicsVector absorption_length;  //!< Absorption length [MeV, len]
+    ImportPhysicsVector component;  //!< Re-emission population [MeV, unitless]
+
+    //! Whether all data are assigned and valid
+    explicit operator bool() const
+    {
+        return mean_num_photons > 0 && time_constant > 0
+               && static_cast<bool>(absorption_length)
+               && static_cast<bool>(component)
+               && absorption_length.vector_type == ImportPhysicsVectorType::free
+               && component.vector_type == absorption_length.vector_type;
+    }
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * Store optical material properties.
@@ -157,13 +183,14 @@ struct ImportOpticalMaterial
     ImportOpticalRayleigh rayleigh;
     ImportOpticalAbsorption absorption;
     ImportOpticalProperty properties;
+    ImportWavelengthShift wls;
 
     //! Whether all data are assigned and valid
     explicit operator bool() const
     {
         return static_cast<bool>(scintillation) || static_cast<bool>(rayleigh)
                || static_cast<bool>(absorption)
-               || static_cast<bool>(properties);
+               || static_cast<bool>(properties) || static_cast<bool>(wls);
     }
 };
 //---------------------------------------------------------------------------//

test/celeritas/ext/GeantImporter.test.cc

@@ -1458,6 +1458,31 @@ TEST_F(LarSphere, optical)
     EXPECT_REAL_EQ(86.4473, to_cm(abs.absorption_length.y.front()));
     EXPECT_REAL_EQ(0.000296154, to_cm(abs.absorption_length.y.back()));
 
+    // Check WLS optical properties
+    auto const& wls = optical.wls;
+    EXPECT_TRUE(wls);
+    EXPECT_REAL_EQ(3, wls.mean_num_photons);
+    EXPECT_REAL_EQ(6e-9, wls.time_constant);
+    EXPECT_EQ(2, wls.absorption_length.x.size());
+    EXPECT_EQ(wls.absorption_length.x.size(), wls.absorption_length.y.size());
+    EXPECT_EQ(ImportPhysicsVectorType::free, wls.absorption_length.vector_type);
+    EXPECT_EQ(wls.component.vector_type, wls.absorption_length.vector_type);
+
+    std::vector<double> abslen_grid, comp_grid;
+    for (auto i : range(wls.absorption_length.x.size()))
+    {
+        abslen_grid.push_back(wls.absorption_length.x[i]);
+        abslen_grid.push_back(wls.absorption_length.y[i]);
+        comp_grid.push_back(wls.component.x[i]);
+        comp_grid.push_back(wls.component.y[i]);
+    }
+
+    static double const expected_abslen_grid[]
+        = {1.3778e-06, 86.4473, 1.55e-05, 0.000296154};
+    static double const expected_comp_grid[] = {1.3778e-06, 10, 1.55e-05, 20};
+    EXPECT_VEC_SOFT_EQ(expected_abslen_grid, abslen_grid);
+    EXPECT_VEC_SOFT_EQ(expected_comp_grid, comp_grid);
+
     // Check common optical properties
     // Refractive index data in the geometry comes from the refractive index
     // database https://refractiveindex.info and was calculating using the

test/geocel/data/lar-sphere.gdml

@@ -36,6 +36,8 @@
     <matrix name="RAY" coldim="2" values="1.55e-06 321429 1.7714e-06 192857 2.102e-06 109286 2.255e-06 64285.7 2.531e-06 39857.1 2.884e-06 27000 3.024e-06 19285.7 4.133e-06 4885.71 6.2e-06 546.429 1.033e-05 546.429 1.55e-05 546.429"/>
     <matrix name="IC" coldim="1" values="3.95306e-5"/>
     <matrix name="ABS" coldim="2" values="1.3778e-06 864.473 1.55e-05 0.00296154"/>
+    <matrix name="WLSC" coldim="2" values="1.3778e-06 10 1.55e-05 20"/>
+    <matrix name="WLSMNP" coldim="1" values="3"/>
   </define>
 
   <materials>
@@ -107,6 +109,10 @@
       <property name="RAYLEIGH" ref="RAY"/>
       <property name="ISOTHERMAL_COMPRESSIBILITY" ref="IC"/>
       <property name="ABSLENGTH" ref="ABS"/>
+      <property name="WLSCOMPONENT" ref="WLSC"/> 
+      <property name="WLSMEANNUMBERPHOTONS" ref="WLSMNP"/> 
+      <property name="WLSABSLENGTH" ref="ABS"/> <!-- Reused for simplicity -->
+      <property name="WLSTIMECONSTANT" ref="TC1"/> <!-- Reused for simplicity -->
     </material>
   </materials>