optical/detail/MatScintSpecInserter.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file celeritas/optical/detail/MatScintSpecInserter.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include "corecel/Assert.hh"
0010 #include "corecel/Types.hh"
0011 #include "corecel/data/CollectionBuilder.hh"
0012 #include "corecel/data/DedupeCollectionBuilder.hh"
0013 #include "corecel/math/SoftEqual.hh"
0014 #include "celeritas/io/ImportOpticalMaterial.hh"
0015
0016 #include "../ScintillationData.hh"
0017
0018 namespace celeritas
0019 {
0020 namespace optical
0021 {
0022 namespace detail
0023 {
0024 //---------------------------------------------------------------------------//
0025 /*!
0026  * Build scintillation spectrum data.
0027  */
0028 class MatScintSpecInserter
0029 {
0030   public:
0031     //!@{
0032     //! \name Type aliases
0033     using Data = HostVal<ScintillationData>;
0034     //!@}
0035
0036   public:
0037     // Construct with data to insert into
0038     explicit inline MatScintSpecInserter(Data* data);
0039
0040     // Add scintillation data for a single material
0041     auto operator()(ImportMaterialScintSpectrum const& mat);
0042
0043   private:
0044     using MatId = OptMatId;
0045
0046     CollectionBuilder<MatScintSpectrumRecord, MemSpace::host, MatId> materials_;
0047     DedupeCollectionBuilder<real_type> reals_;
0048     CollectionBuilder<ScintRecord> scint_records_;
0049 };
0050
0051 //---------------------------------------------------------------------------//
0052 // INLINE DEFINITIONS
0053 //---------------------------------------------------------------------------//
0054 /*!
0055  * Construct with defaults.
0056  */
0057 MatScintSpecInserter::MatScintSpecInserter(Data* data)
0058     : materials_{&data->materials}
0059     , reals_{&data->reals}
0060     , scint_records_{&data->scint_records}
0061 {
0062     CELER_EXPECT(data);
0063 }
0064
0065 //---------------------------------------------------------------------------//
0066 /*!
0067  * Add scintillation data for a single material.
0068  */
0069 auto MatScintSpecInserter::operator()(ImportMaterialScintSpectrum const& mat)
0070 {
0071     CELER_EXPECT(mat);
0072
0073     CELER_VALIDATE(mat.yield_per_energy > 0,
0074                    << "invalid yield=" << mat.yield_per_energy
0075                    << " for scintillation (should be positive)");
0076
0077     double total_yield{0};
0078     std::vector<double> yield_pdf;
0079     auto const begin_components = scint_records_.size_id();
0080     for (ImportScintComponent const& comp : mat.components)
0081     {
0082         CELER_VALIDATE(comp.lambda_mean > 0,
0083                        << "invalid lambda_mean=" << comp.lambda_mean
0084                        << " for scintillation component (should be positive)");
0085         CELER_VALIDATE(comp.lambda_sigma > 0,
0086                        << "invalid lambda_sigma=" << comp.lambda_sigma
0087                        << " (should be positive)");
0088         CELER_VALIDATE(comp.rise_time >= 0,
0089                        << "invalid rise_time=" << comp.rise_time
0090                        << " (should be nonnegative)");
0091         CELER_VALIDATE(comp.fall_time > 0,
0092                        << "invalid fall_time=" << comp.fall_time
0093                        << " (should be positive)");
0094         CELER_VALIDATE(comp.yield_frac > 0,
0095                        << "invalid yield=" << comp.yield_frac);
0096
0097         ScintRecord scint;
0098         scint.lambda_mean = comp.lambda_mean;
0099         scint.lambda_sigma = comp.lambda_sigma;
0100         scint.rise_time = comp.rise_time;
0101         scint.fall_time = comp.fall_time;
0102         scint_records_.push_back(scint);
0103
0104         yield_pdf.push_back(comp.yield_frac);
0105         total_yield += comp.yield_frac;
0106     }
0107
0108     // Normalize yield PDF by total yield
0109     for (auto& y : yield_pdf)
0110     {
0111         y /= total_yield;
0112     }
0113
0114     MatScintSpectrumRecord spectrum;
0115     spectrum.yield_per_energy = mat.yield_per_energy;
0116     spectrum.components = {begin_components, scint_records_.size_id()};
0117     spectrum.yield_pdf = reals_.insert_back(yield_pdf.begin(), yield_pdf.end());
0118
0119     CELER_ENSURE(spectrum.components.size() == mat.components.size());
0120     return materials_.push_back(std::move(spectrum));
0121 }
0122
0123 //---------------------------------------------------------------------------//
0124 }  // namespace detail
0125 }  // namespace optical
0126 }  // namespace celeritas