Add particle scintillation data to ScintillationParams and add `Sci…

…ntillationPreGenerator` (celeritas-project#1153)

app/celer-dump-data.cc

@@ -667,13 +667,14 @@ void print_optical_material_data(ImportData::ImportOpticalMap const& iom)
     for (auto const& [mid, val] : iom)
     {
         auto const& scint = val.scintillation;
-        cout << POM_STREAM_SCALAR(mid, scint, material.yield, IU::inv_mev);
+        cout << POM_STREAM_SCALAR(
+            mid, scint, material.yield_per_energy, IU::inv_mev);
         cout << POM_STREAM_SCALAR(mid, scint, resolution_scale, IU::unitless);
         for (auto i : range(scint.material.components.size()))
         {
             auto const& comp = scint.material.components[i];
             cout << POM_STREAM_SCALAR_COMP(
-                mid, comp, yield, IU::inv_mev, comp_str[i]);
+                mid, comp, yield_per_energy, IU::inv_mev, comp_str[i]);
             cout << POM_STREAM_SCALAR_COMP(
                 mid, comp, lambda_mean, IU::len, comp_str[i]);
             cout << POM_STREAM_SCALAR_COMP(

src/celeritas/Types.hh

@@ -63,9 +63,6 @@ using ElementComponentId = OpaqueId<struct MatElementComponent>;
 //! Opaque index to one isotopic component datum in a particular element
 using IsotopeComponentId = OpaqueId<struct ElIsotopeComponent>;
 
-//! Opaque index to a material with optical properties
-using OpticalMaterialId = OpaqueId<struct OpticalMaterial_>;
-
 //! Opaque index of a process applicable to a single particle type
 using ParticleProcessId = OpaqueId<ProcessId>;
 

src/celeritas/ext/GeantImporter.cc

@@ -229,7 +229,7 @@ fill_vec_import_scint_comp(MatPropGetter& get_property,
     for (int comp_idx : range(1, 4))
     {
         ImportScintComponent comp;
-        get_property.scalar(&comp.yield,
+        get_property.scalar(&comp.yield_per_energy,
                             particle_name + "SCINTILLATIONYIELD",
                             comp_idx,
                             ImportUnits::inv_mev);
@@ -504,9 +504,10 @@ ImportData::ImportOpticalMap import_optical()
         // Save scintillation properties
         {
             // Material scintillation properties
-            get_property.scalar(&optical.scintillation.material.yield,
-                                "SCINTILLATIONYIELD",
-                                ImportUnits::inv_mev);
+            get_property.scalar(
+                &optical.scintillation.material.yield_per_energy,
+                "SCINTILLATIONYIELD",
+                ImportUnits::inv_mev);
             get_property.scalar(&optical.scintillation.resolution_scale,
                                 "RESOLUTIONSCALE",
                                 ImportUnits::unitless);

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::ImportOpticalParameters+;
 #pragma link C++ class celeritas::ImportWavelengthShift+;
 #pragma link C++ class celeritas::ImportOpticalProperty+;
 #pragma link C++ class celeritas::ImportOpticalMaterial+;

src/celeritas/io/ImportData.hh

@@ -72,6 +72,7 @@ struct ImportData
     ImportOpticalMap optical;
     ImportEmParameters em_params;
     ImportTransParameters trans_params;
+    ImportOpticalParameters optical_params;
     ImportSBMap sb_data;
     ImportLivermorePEMap livermore_pe_data;
     ImportNeutronElasticMap neutron_elastic_data;

src/celeritas/io/ImportOpticalMaterial.hh

@@ -12,6 +12,7 @@
 #include "celeritas/phys/PDGNumber.hh"
 
 #include "ImportPhysicsVector.hh"
+#include "ImportUnits.hh"
 
 namespace celeritas
 {
@@ -24,7 +25,7 @@ namespace celeritas
  */
 struct ImportScintComponent
 {
-    double yield{};  //!< Yield for this material component [1/MeV]
+    double yield_per_energy{};  //!< Yield for this material component [1/MeV]
     double lambda_mean{};  //!< Mean wavelength [len]
     double lambda_sigma{};  //!< Standard deviation of wavelength
     double rise_time{};  //!< Rise time [time]
@@ -33,7 +34,7 @@ struct ImportScintComponent
     //! Whether all data are assigned and valid
     explicit operator bool() const
     {
-        return yield > 0 && lambda_mean > 0 && lambda_sigma > 0
+        return yield_per_energy > 0 && lambda_mean > 0 && lambda_sigma > 0
                && rise_time >= 0 && fall_time > 0;
     }
 };
@@ -45,12 +46,15 @@ struct ImportScintComponent
  */
 struct ImportMaterialScintSpectrum
 {
-    double yield{};  //!< Characteristic light yields of the material [1/MeV]
+    double yield_per_energy{};  //!< Light yield of the material [1/MeV]
     std::vector<ImportScintComponent> components;  //!< Scintillation
                                                    //!< components
 
     //! Whether all data are assigned and valid
-    explicit operator bool() const { return yield > 0 && !components.empty(); }
+    explicit operator bool() const
+    {
+        return yield_per_energy > 0 && !components.empty();
+    }
 };
 
 //---------------------------------------------------------------------------//
@@ -64,7 +68,12 @@ struct ImportMaterialScintSpectrum
  */
 struct ImportParticleScintSpectrum
 {
-    ImportPhysicsVector yield_vector;  //!< Particle yield vector
+#ifndef SWIG
+    static constexpr auto x_units{ImportUnits::mev};
+    static constexpr auto y_units{ImportUnits::unitless};
+#endif
+
+    ImportPhysicsVector yield_vector;  //!< Particle yield per energy bin
     std::vector<ImportScintComponent> components;  //!< Scintillation
                                                    //!< components
 
@@ -92,7 +101,8 @@ struct ImportScintData
     //! Whether all data are assigned and valid
     explicit operator bool() const
     {
-        return static_cast<bool>(material) && resolution_scale >= 0;
+        return (static_cast<bool>(material) || !particles.empty())
+               && resolution_scale >= 0;
     }
 };
 

src/celeritas/io/ImportParameters.hh

@@ -117,5 +117,15 @@ struct ImportTransParameters
     }
 };
 
+//---------------------------------------------------------------------------//
+/*!
+ * TODO: Placeholder for optical parameter data.
+ * See \c G4OpticalParameters .
+ */
+struct ImportOpticalParameters
+{
+    bool scintillation_by_particle{false};
+};
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/mat/MaterialData.hh

@@ -13,6 +13,7 @@
 #include "corecel/data/CollectionBuilder.hh"
 #include "celeritas/Quantities.hh"
 #include "celeritas/Types.hh"
+#include "celeritas/optical/Types.hh"
 #include "celeritas/phys/AtomicNumber.hh"
 
 namespace celeritas

src/celeritas/optical/CerenkovData.hh

@@ -12,6 +12,7 @@
 #include "corecel/data/Collection.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/grid/GenericGridData.hh"
+#include "celeritas/optical/Types.hh"
 
 namespace celeritas
 {

src/celeritas/optical/CerenkovPreGenerator.hh

@@ -48,6 +48,7 @@ class CerenkovPreGenerator
 {
   public:
     // Placeholder for data that is not available through Views
+    // TODO: Merge Cerenkov and Scintillation pre-gen data into one struct
     struct OpticalPreGenStepData
     {
         real_type time{};  //!< Pre-step time

src/celeritas/optical/OpticalDistributionData.hh

@@ -12,6 +12,7 @@
 #include "corecel/cont/EnumArray.hh"
 #include "celeritas/Quantities.hh"
 #include "celeritas/Types.hh"
+#include "celeritas/optical/Types.hh"
 
 namespace celeritas
 {

src/celeritas/optical/OpticalPropertyData.hh

@@ -12,6 +12,7 @@
 #include "corecel/data/Collection.hh"
 #include "celeritas/Types.hh"
 #include "celeritas/grid/GenericGridData.hh"
+#include "celeritas/optical/Types.hh"
 
 namespace celeritas
 {

src/celeritas/optical/OpticalPropertyParams.cc

@@ -12,7 +12,6 @@
 
 #include "corecel/cont/Range.hh"
 #include "corecel/data/CollectionBuilder.hh"
-#include "corecel/data/DedupeCollectionBuilder.hh"
 #include "corecel/math/Algorithms.hh"
 #include "celeritas/Quantities.hh"
 #include "celeritas/Types.hh"

src/celeritas/optical/ScintillationData.hh

@@ -11,16 +11,22 @@
 #include "corecel/Types.hh"
 #include "corecel/data/Collection.hh"
 #include "celeritas/Types.hh"
+#include "celeritas/grid/GenericGridData.hh"
+#include "celeritas/optical/Types.hh"
 
 namespace celeritas
 {
 //---------------------------------------------------------------------------//
 /*!
  * Material dependent scintillation property.
+ *
+ * Components represent different scintillation emissions, such as
+ * prompt/fast, intermediate, and slow. They can be material-only or depend on
+ * the incident particle type.
  */
 struct ScintillationComponent
 {
-    real_type yield_prob{};  //!< Probability of the yield
+    real_type yield_frac{};  //!< Fraction of total yield (yield/sum(yields))
     real_type lambda_mean{};  //!< Mean wavelength
     real_type lambda_sigma{};  //!< Standard dev. of wavelength
     real_type rise_time{};  //!< Rise time
@@ -29,35 +35,77 @@ struct ScintillationComponent
     //! Whether all data are assigned and valid
     explicit CELER_FUNCTION operator bool() const
     {
-        return yield_prob > 0 && lambda_mean > 0 && lambda_sigma > 0
-               && rise_time >= 0 && fall_time > 0;
+        return yield_frac > 0 && yield_frac <= 1 && lambda_mean > 0
+               && lambda_sigma > 0 && rise_time >= 0 && fall_time > 0;
     }
 };
 
 //---------------------------------------------------------------------------//
 /*!
- * Data characterizing the scintillation spectrum.
+ * Data characterizing material-only scintillation spectrum information.
  *
- * \c yield is the characteristic light yield of the material.
- * \c resolution_scale scales the standard deviation of the distribution of the
- * number of photons generated.
+ * - \c yield_per_energy is the characteristic light yield of the material in
+ *   [1/MeV] units. The total light yield per step is then
+ *   `yield_per_energy * energy_dep`, which results in a (unitless) number of
+ *   photons.
+ * - \c resolution_scale scales the standard deviation of the distribution of
+ *   the number of photons generated.
+ * - \c components stores the different scintillation components
+ *   (fast/slow/etc) for this material.
  */
-struct ScintillationSpectrum
+struct MaterialScintillationSpectrum
 {
-    real_type yield{};
-    real_type resolution_scale{};
+    real_type yield_per_energy{};  //!< [1/MeV]
     ItemRange<ScintillationComponent> components;
 
     //! Whether all data are assigned and valid
     explicit CELER_FUNCTION operator bool() const
     {
-        return yield > 0 && resolution_scale >= 0 && !components.empty();
+        return yield_per_energy > 0 && !components.empty();
     }
 };
 
 //---------------------------------------------------------------------------//
 /*!
- * Scintillation data tabulated with the optical material id.
+ * Data characterizing the scintillation spectrum for a given particle in a
+ * given material.
+ *
+ * \c yield_vector is the characteristic light yield for different energies.
+ * \c components stores the fast/slow/etc scintillation components for this
+ * particle type.
+ */
+struct ParticleScintillationSpectrum
+{
+    GenericGridData yield_vector;
+    ItemRange<ScintillationComponent> components;
+
+    //! Whether all data are assigned and valid
+    explicit CELER_FUNCTION operator bool() const
+    {
+        return static_cast<bool>(yield_vector);
+    }
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Data characterizing the scintillation spectrum for all particles and
+ * materials.
+ *
+ * Sampling using material-only data or particle- and material-dependent data
+ * are mutually exclusive. Therefore, either \c materials or \c particles are
+ * loaded at the beginning of the simulation, but *never* both at the same
+ * time. The \c scintillation_by_particle() function can be used to check that.
+ *
+ * - \c pid_to_scintpid returns a \c ScintillationParticleId given a
+ *   \c ParticleId .
+ * - \c resolution_scale is indexed by \c OpticalMaterialId .
+ * - \c materials stores material-only scintillation data. Indexed by
+ *   \c OpticalMaterialId
+ * - \c particles stores scintillation spectrum for each particle type for each
+ *   material, being a grid of size `num_particles * num_materials`. Therefore
+ *   it is indexed by \c ParticleScintSpectrumId , which combines
+ *   \c ScintillationParticleId and \c OpticalMaterialId . Use the
+ *   \c spectrum_index() function to retrieve the correct index.
  */
 template<Ownership W, MemSpace M>
 struct ScintillationData
@@ -66,30 +114,73 @@ struct ScintillationData
     using Items = Collection<T, W, M>;
     template<class T>
     using OpticalMaterialItems = Collection<T, W, M, OpticalMaterialId>;
+    using ParticleScintillationSpectra
+        = Collection<ParticleScintillationSpectrum, W, M, ParticleScintSpectrumId>;
 
     //// MEMBER DATA ////
 
+    //! Resolution scale for each material [OpticalMaterialid]
+    OpticalMaterialItems<real_type> resolution_scale;
+
+    //! Material-only scintillation spectrum data [OpticalMaterialid]
+    OpticalMaterialItems<MaterialScintillationSpectrum> materials;
+
+    //! Index between ScintillationParticleId and ParticleId
+    Collection<ScintillationParticleId, W, M, ParticleId> pid_to_scintpid;
+    //! Cache number of scintillation particles; Used by this->spectrum_index
+    size_type num_scint_particles{};
+    //! Particle/material scintillation spectrum data [ParticleScintSpectrumId]
+    ParticleScintillationSpectra particles;
+    //! Backend storage for ParticleScintillationSpectrum::yield_vector
+    Items<real_type> reals;
+
+    //! Components for either material or particle items
     Items<ScintillationComponent> components;
-    OpticalMaterialItems<ScintillationSpectrum> spectra;
 
     //// MEMBER FUNCTIONS ////
 
     //! Whether all data are assigned and valid
     explicit CELER_FUNCTION operator bool() const
     {
-        return !components.empty() && !spectra.empty();
+        return !resolution_scale.empty()
+               && (materials.empty() != particles.empty())
+               && (!pid_to_scintpid.empty() == !particles.empty())
+               && (!pid_to_scintpid.empty() == (num_scint_particles > 0));
+    }
+
+    //! Whether sampling must happen by particle type
+    CELER_FUNCTION bool scintillation_by_particle() const
+    {
+        return !particles.empty();
+    }
+
+    //! Retrieve spectrum index given optical particle and material ids
+    ParticleScintSpectrumId
+    spectrum_index(ScintillationParticleId pid, OpticalMaterialId mat_id) const
+    {
+        // Resolution scale exists independent of material-only data and it's
+        // indexed by optical material id
+        CELER_EXPECT(pid < num_scint_particles
+                     && mat_id < resolution_scale.size());
+        return ParticleScintSpectrumId{resolution_scale.size() * pid.get()
+                                       + mat_id.get()};
     }
 
     //! Assign from another set of data
     template<Ownership W2, MemSpace M2>
     ScintillationData& operator=(ScintillationData<W2, M2> const& other)
     {
         CELER_EXPECT(other);
+        resolution_scale = other.resolution_scale;
+        materials = other.materials;
+        pid_to_scintpid = other.pid_to_scintpid;
+        num_scint_particles = other.num_scint_particles;
+        particles = other.particles;
+        reals = other.reals;
         components = other.components;
-        spectra = other.spectra;
         return *this;
     }
 };
 
 //---------------------------------------------------------------------------//
-}  // namespace celeritas
+}  // namespace celeritas