Import scintillation particle data (celeritas-project#1151)

* Update `ImportOpticalMaterial` to include particle data
* Import particle scintillation data
* Update `GeantImporter` test
* Fix celer-dump-data
* Temporary fix before ScintillationData is fully overhauled
* Fix bug; Add particle scintillation import test
* Fix root update command
* Move new enum to end of list

app/celer-dump-data.cc

@@ -667,13 +667,13 @@ 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, yield, IU::unitless);
+        cout << POM_STREAM_SCALAR(mid, scint, material.yield, IU::inv_mev);
         cout << POM_STREAM_SCALAR(mid, scint, resolution_scale, IU::unitless);
-        for (auto i : range(scint.components.size()))
+        for (auto i : range(scint.material.components.size()))
         {
-            auto const& comp = scint.components[i];
+            auto const& comp = scint.material.components[i];
             cout << POM_STREAM_SCALAR_COMP(
-                mid, comp, yield, IU::unitless, comp_str[i]);
+                mid, comp, yield, 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/ext/GeantImporter.cc

@@ -192,6 +192,71 @@ struct MatPropGetter
     }
 };
 
+//---------------------------------------------------------------------------//
+//! Map particles defined in \c G4MaterialConstPropertyIndex .
+auto& optical_particles_map()
+{
+    static std::unordered_map<std::string, PDGNumber> const map
+        = {{"PROTON", pdg::proton()},
+           {"DEUTERON", pdg::deuteron()},
+           {"TRITON", pdg::triton()},
+           {"ALPHA", pdg::alpha()},
+           {"ION", pdg::ion()},
+           {"ELECTRON", pdg::electron()}};
+    return map;
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Populate an \c ImportScintComponent .
+ * To retrieve a material-only component simply do not use particle name.
+ */
+std::vector<ImportScintComponent>
+fill_vec_import_scint_comp(MatPropGetter& get_property,
+                           std::string particle_name = "")
+{
+    CELER_EXPECT(particle_name.empty()
+                 || optical_particles_map().count(particle_name));
+
+    std::vector<ImportScintComponent> components;
+    for (int comp_idx : range(1, 4))
+    {
+        ImportScintComponent comp;
+        get_property.scalar(&comp.yield,
+                            particle_name + "SCINTILLATIONYIELD",
+                            comp_idx,
+                            ImportUnits::inv_mev);
+
+        // Custom-defined properties not available in G4MaterialPropertyIndex
+        {
+            get_property.scalar(&comp.lambda_mean,
+                                particle_name + "SCINTILLATIONLAMBDAMEAN",
+                                comp_idx,
+                                ImportUnits::len);
+            get_property.scalar(&comp.lambda_sigma,
+                                particle_name + "SCINTILLATIONLAMBDASIGMA",
+                                comp_idx,
+                                ImportUnits::len);
+        }
+
+        // Rise time is not defined for particle type in Geant4
+        get_property.scalar(&comp.rise_time,
+                            particle_name + "SCINTILLATIONRISETIME",
+                            comp_idx,
+                            ImportUnits::time);
+
+        get_property.scalar(&comp.fall_time,
+                            particle_name + "SCINTILLATIONTIMECONSTANT",
+                            comp_idx,
+                            ImportUnits::time);
+        if (comp)
+        {
+            components.push_back(std::move(comp));
+        }
+    }
+    return components;
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Safely switch from G4State [G4Material.hh] to ImportMaterialState.
@@ -404,7 +469,10 @@ ImportData::ImportOpticalMap import_optical()
     auto num_materials = pct.GetTableSize();
     CELER_ASSERT(num_materials > 0);
 
+    auto const& particle_map = optical_particles_map();
     ImportData::ImportOpticalMap result;
+
+    // Loop over optical materials
     for (auto mat_idx : range(num_materials))
     {
         auto const* mcc = pct.GetMaterialCutsCouple(mat_idx);
@@ -429,38 +497,34 @@ ImportData::ImportOpticalMap import_optical()
                             ImportUnits::unitless);
 
         // Save scintillation properties
-        get_property.scalar(&optical.scintillation.resolution_scale,
-                            "RESOLUTIONSCALE",
-                            ImportUnits::unitless);
-        get_property.scalar(&optical.scintillation.yield,
-                            "SCINTILLATIONYIELD",
-                            ImportUnits::unitless);
-        for (int comp_idx : range(1, 4))
         {
-            ImportScintComponent comp;
-            get_property.scalar(&comp.yield,
+            // Material scintillation properties
+            get_property.scalar(&optical.scintillation.material.yield,
                                 "SCINTILLATIONYIELD",
-                                comp_idx,
+                                ImportUnits::inv_mev);
+            get_property.scalar(&optical.scintillation.resolution_scale,
+                                "RESOLUTIONSCALE",
                                 ImportUnits::unitless);
-            get_property.scalar(&comp.lambda_mean,
-                                "SCINTILLATIONLAMBDAMEAN",
-                                comp_idx,
-                                ImportUnits::len);
-            get_property.scalar(&comp.lambda_sigma,
-                                "SCINTILLATIONLAMBDASIGMA",
-                                comp_idx,
-                                ImportUnits::len);
-            get_property.scalar(&comp.rise_time,
-                                "SCINTILLATIONRISETIME",
-                                comp_idx,
-                                ImportUnits::time);
-            get_property.scalar(&comp.fall_time,
-                                "SCINTILLATIONTIMECONSTANT",
-                                comp_idx,
-                                ImportUnits::time);
-            if (comp)
+            optical.scintillation.material.components
+                = fill_vec_import_scint_comp(get_property);
+
+            // Particle scintillation properties
+            for (auto const& iter : particle_map)
             {
-                optical.scintillation.components.push_back(comp);
+                auto const& particle_name = iter.first;
+
+                ImportScintData::IPSS scint_part_spec;
+                get_property.vector(&scint_part_spec.yield_vector,
+                                    particle_name + "SCINTILLATIONYIELD",
+                                    ImportUnits::inv_mev);
+                scint_part_spec.components
+                    = fill_vec_import_scint_comp(get_property, particle_name);
+
+                if (scint_part_spec)
+                {
+                    optical.scintillation.particles.insert(
+                        {iter.second.get(), std::move(scint_part_spec)});
+                }
             }
         }
 
@@ -484,6 +548,7 @@ ImportData::ImportOpticalMap import_optical()
             result[mat_idx] = optical;
         }
     }
+
     CELER_LOG(debug) << "Loaded " << result.size() << " optical materials";
     return result;
 }
@@ -507,7 +572,8 @@ import_materials(GeantImporter::DataSelection::Flags particle_flags)
     std::vector<ImportMaterial> materials;
     materials.resize(g4production_cuts_table.GetTableSize());
     CELER_VALIDATE(!materials.empty(),
-                   << "no Geant4 production cuts are defined (you may need "
+                   << "no Geant4 production cuts are defined (you may "
+                      "need "
                       "to call G4RunManager::RunInitialization)");
 
     using CutRange = std::pair<G4ProductionCutsIndex,

src/celeritas/ext/RootInterfaceLinkDef.h

@@ -19,7 +19,9 @@
 #pragma link C++ class celeritas::ImportPhysicsVector+;
 #pragma link C++ class celeritas::ImportMaterial+;
 #pragma link C++ class celeritas::ImportScintComponent+;
-#pragma link C++ class celeritas::ImportScintSpectrum+;
+#pragma link C++ class celeritas::ImportMaterialScintSpectrum+;
+#pragma link C++ class celeritas::ImportParticleScintSpectrum+;
+#pragma link C++ class celeritas::ImportScintData+;
 #pragma link C++ class celeritas::ImportOpticalRayleigh+;
 #pragma link C++ class celeritas::ImportOpticalAbsorption+;
 #pragma link C++ class celeritas::ImportOpticalProperty+;

src/celeritas/io/ImportOpticalMaterial.hh

@@ -9,22 +9,28 @@
 
 #include <vector>
 
+#include "celeritas/phys/PDGNumber.hh"
+
 #include "ImportPhysicsVector.hh"
 
 namespace celeritas
 {
 //---------------------------------------------------------------------------//
 /*!
- * Store material-dependent properties for scintillation spectrum.
+ * Store basic properties for different scintillation component types.
+ *
+ * Fast/intermediate/slow/etc scintillation components can be used for both
+ * particle- and material-dependent spectra, as well as material-only spectra.
  */
 struct ImportScintComponent
 {
-    double yield{};  //!< Yield for this component
-    double lambda_mean{};  //!< Mean wavelength
+    double yield{};  //!< 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
-    double fall_time{};  //!< Decay time
+    double rise_time{};  //!< Rise time [time]
+    double fall_time{};  //!< Decay time [time]
 
+    //! Whether all data are assigned and valid
     explicit operator bool() const
     {
         return yield > 0 && lambda_mean > 0 && lambda_sigma > 0
@@ -34,17 +40,59 @@ struct ImportScintComponent
 
 //---------------------------------------------------------------------------//
 /*!
- * Store optical material properties for scintillation.
+ * Store material-only scintillation spectrum information.
+ * TODO: Components are not necessary in Geant4, but are in our generator.
  */
-struct ImportScintSpectrum
+struct ImportMaterialScintSpectrum
 {
-    double yield{};  //!< Characteristic light yield of the material
+    double yield{};  //!< Characteristic light yields 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(); }
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Store per-particle material scintillation spectrum information.
+ *
+ * The yield vector is the only necessary element, needed to calculate the
+ * yield based on the particle energy-loss during the stepping loop.
+ * Components may not be assigned---they are the equivalent of
+ * \c k[Particle]ScintillationYield[i] in \c G4MaterialPropertiesIndex.hh
+ */
+struct ImportParticleScintSpectrum
+{
+    ImportPhysicsVector yield_vector;  //!< Particle yield vector
+    std::vector<ImportScintComponent> components;  //!< Scintillation
+                                                   //!< components
+
+    //! Whether all data are assigned and valid
+    explicit operator bool() const
+    {
+        return yield_vector
+               && yield_vector.vector_type == ImportPhysicsVectorType::free;
+    }
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Store optical properties for scintillation.
+ */
+struct ImportScintData
+{
+    using PDGint = int;
+    using IPSS = ImportParticleScintSpectrum;
+
+    ImportMaterialScintSpectrum material;  //!< Material scintillation data
+    std::map<PDGint, IPSS> particles;  //!< Particle scintillation data
     double resolution_scale{};  //!< Scales the stdev of photon distribution
-    std::vector<ImportScintComponent> components;
 
+    //! Whether all data are assigned and valid
     explicit operator bool() const
     {
-        return yield > 0 && resolution_scale >= 0 && !components.empty();
+        return static_cast<bool>(material) && resolution_scale >= 0;
     }
 };
 
@@ -61,6 +109,7 @@ struct ImportOpticalRayleigh
     double compressibility{};  //!< Isothermal compressibility
     ImportPhysicsVector mfp;  //!< Rayleigh mean free path
 
+    //! Whether all data are assigned and valid
     explicit operator bool() const
     {
         return scale_factor >= 0
@@ -76,6 +125,7 @@ struct ImportOpticalAbsorption
 {
     ImportPhysicsVector absorption_length;
 
+    //! Whether all data are assigned and valid
     explicit operator bool() const
     {
         return static_cast<bool>(absorption_length);
@@ -90,6 +140,7 @@ struct ImportOpticalProperty
 {
     ImportPhysicsVector refractive_index;
 
+    //! Whether all data are assigned and valid
     explicit operator bool() const
     {
         return static_cast<bool>(refractive_index);
@@ -102,11 +153,12 @@ struct ImportOpticalProperty
  */
 struct ImportOpticalMaterial
 {
-    ImportScintSpectrum scintillation;
+    ImportScintData scintillation;
     ImportOpticalRayleigh rayleigh;
     ImportOpticalAbsorption absorption;
     ImportOpticalProperty properties;
 
+    //! Whether all data are assigned and valid
     explicit operator bool() const
     {
         return static_cast<bool>(scintillation) || static_cast<bool>(rayleigh)

src/celeritas/io/ImportUnits.cc

@@ -83,6 +83,7 @@ char const* to_cstring(ImportUnits value)
         "time",
         "1/len^3",
         "len-time^2/mass",
+        "1/MeV"
     };
     return to_cstring_impl(value);
 }
@@ -103,6 +104,8 @@ double native_value_from(UnitSystem sys, ImportUnits q)
             return 1;
         case ImportUnits::mev:
             return mev;
+        case ImportUnits::inv_mev:
+            return 1 / mev;
         case ImportUnits::mev_per_len:
             return mev / len;
         case ImportUnits::len:

src/celeritas/io/ImportUnits.hh

@@ -33,6 +33,7 @@ enum class ImportUnits
     time,  //!< Time [time]
     inv_len_cb,  //!< Number density [1/len^3]
     len_time_sq_per_mass,  //!< Inverse pressure [len-time^2/mass]
+    inv_mev,  //!< Inverse energy [1/MeV]
     size_,
     // Deprecated aliases
     none = unitless,  //!< Deprecated

src/celeritas/io/detail/ImportDataConverter.cc

@@ -89,13 +89,23 @@ void ImportDataConverter::operator()(ImportOpticalMaterial* data)
 {
     CELER_EXPECT(data);
 
-    for (auto& comp : data->scintillation.components)
+    for (auto& comp : data->scintillation.material.components)
     {
         comp.lambda_mean *= len_;
         comp.lambda_sigma *= len_;
         comp.rise_time *= time_;
         comp.fall_time *= time_;
     }
+    for (auto& iter : data->scintillation.particles)
+    {
+        for (auto& comp : iter.second.components)
+        {
+            comp.lambda_mean *= len_;
+            comp.lambda_sigma *= len_;
+            comp.rise_time *= time_;
+            comp.fall_time *= time_;
+        }
+    }
     for (auto& mfp : data->rayleigh.mfp.y)
     {
         mfp *= len_;

src/celeritas/optical/ScintillationParams.cc

@@ -58,7 +58,7 @@ ScintillationParams::ScintillationParams(Input const& input)
     for (auto const& spec : input.data)
     {
         // Check validity of scintillation data
-        auto const& comp_inp = spec.components;
+        auto const& comp_inp = spec.material.components;
         CELER_ASSERT(!comp_inp.empty());
         std::vector<ScintillationComponent> comp(comp_inp.size());
         real_type norm{0};
@@ -95,7 +95,7 @@ ScintillationParams::ScintillationParams(Input const& input)
                            << " (should be positive)");
         }
         ScintillationSpectrum spectrum;
-        spectrum.yield = spec.yield;
+        spectrum.yield = spec.material.yield;
         CELER_VALIDATE(spectrum.yield > 0,
                        << "invalid yield=" << spectrum.yield
                        << " for scintillation (should be positive)");