Skip to content

Commit

Permalink
Support importing optical properties from Geant4 (celeritas-project#1103
Browse files Browse the repository at this point in the history 
)

* Add optical property params and import properties from Geant4

* Add GDML and test for imported optical material properties

* Refactor optical property import

* Build scintillation params from imported data

* Separate imported optical properties from material

* Add helper struct for getting optical properties from the MPT and import rayleigh and absorption properties

* Validate input

* Print optical properties

* Fix test for other units systems and check optical material ID

* Convert Rayleigh and absorption optical property units in ImportDataConverter

* Address review feedback
  • Loading branch information
@amandalund
amandalund committed Feb 13, 2024
1 parent 276c2cd commit 4c469e2
Show file tree
Hide file tree
Showing 26 changed files with 986 additions and 91 deletions.
93 changes: 93 additions & 0 deletions app/celer-dump-data.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -614,6 +614,98 @@ void print_atomic_relaxation_data(
cout << endl;
}

//---------------------------------------------------------------------------//
/*!
* Print optical material properties map.
*/
void print_optical_material_data(ImportData::ImportOpticalMap const& iom)
{
if (iom.empty())
{
CELER_LOG(info) << "Optical material data not available";
return;
}

CELER_LOG(info) << "Loaded optical material data map with size "
<< iom.size();

#define POM_STREAM_SCALAR_COMP(ID, STRUCT, NAME, UNITS, COMP) \
"| " << setw(11) << ID << " | " << setw(20) << #NAME << COMP << " | " \
<< setw(15) << to_cstring(UNITS) << " | " << setprecision(3) \
<< setw(9) << STRUCT.NAME << " | " << setw(52) << "" \
<< " | " << setw(7) << "" \
<< " |\n"
#define POM_STREAM_SCALAR(ID, STRUCT, NAME, UNITS) \
POM_STREAM_SCALAR_COMP(ID, STRUCT, NAME, UNITS, " ")
#define POM_STREAM_VECTOR(ID, STRUCT, NAME, UNITS) \
"| " << setw(11) << ID << " | " << setw(26) << #NAME << " | " << setw(15) \
<< to_cstring(UNITS) << " | " << setw(9) << "" \
<< " | (" << setprecision(3) << setw(10) << STRUCT.NAME.x.front() \
<< ", " << setprecision(3) << setw(10) << STRUCT.NAME.y.front() \
<< ") -> (" << setprecision(3) << setw(10) << STRUCT.NAME.x.back() \
<< ", " << setprecision(3) << setw(10) << STRUCT.NAME.y.back() \
<< ") | " << setw(7) << STRUCT.NAME.x.size() << " |\n";
std::string header = R"gfm(
| Material ID | Property | Units | Scalar | Vector endpoints (MeV, value) | Size |
| ----------- | -------------------------- | --------------- | --------- | ---------------------------------------------------- | ------- |
)gfm";

using IU = ImportUnits;

cout << "\n# Optical properties\n";
cout << "\n## Common properties";
cout << header;
for (auto const& [mid, val] : iom)
{
auto const& prop = val.properties;
cout << POM_STREAM_VECTOR(mid, prop, refractive_index, IU::unitless);
}
cout << "\n## Scintillation";
cout << header;
char const* comp_str[] = {"(fast)", " (mid)", "(slow)"};
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, resolution_scale, IU::unitless);
for (auto i : range(scint.components.size()))
{
auto const& comp = scint.components[i];
cout << POM_STREAM_SCALAR_COMP(
mid, comp, yield, IU::unitless, comp_str[i]);
cout << POM_STREAM_SCALAR_COMP(
mid, comp, lambda_mean, IU::len, comp_str[i]);
cout << POM_STREAM_SCALAR_COMP(
mid, comp, lambda_sigma, IU::len, comp_str[i]);
cout << POM_STREAM_SCALAR_COMP(
mid, comp, rise_time, IU::time, comp_str[i]);
cout << POM_STREAM_SCALAR_COMP(
mid, comp, fall_time, IU::time, comp_str[i]);
}
}
cout << "\n## Rayleigh";
cout << header;
for (auto const& [mid, val] : iom)
{
auto const& rayl = val.rayleigh;
cout << POM_STREAM_SCALAR(mid, rayl, scale_factor, IU::unitless);
cout << POM_STREAM_SCALAR(
mid, rayl, compressibility, IU::len_time_sq_per_mass);
cout << POM_STREAM_VECTOR(mid, rayl, mfp, IU::len);
}
cout << "\n## Absorption";
cout << header;
for (auto const& [mid, val] : iom)
{
auto const& abs = val.absorption;
cout << POM_STREAM_VECTOR(mid, abs, absorption_length, IU::len);
}
cout << endl;
#undef PEP_STREAM_SCALAR
#undef PEP_STREAM_VECTOR
}

//---------------------------------------------------------------------------//
} // namespace
} // namespace app
Expand Down Expand Up @@ -677,6 +769,7 @@ int main(int argc, char* argv[])
print_sb_data(data.sb_data);
print_livermore_pe_data(data.livermore_pe_data);
print_atomic_relaxation_data(data.atomic_relaxation_data);
print_optical_material_data(data.optical);

return EXIT_SUCCESS;
}
1 change: 1 addition & 0 deletions src/celeritas/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@ list(APPEND SOURCES
mat/detail/Utils.cc
optical/CerenkovParams.cc
optical/ScintillationParams.cc
optical/OpticalPropertyParams.cc
phys/CutoffParams.cc
phys/ImportedModelAdapter.cc
phys/ImportedProcessAdapter.cc
Expand Down
138 changes: 138 additions & 0 deletions src/celeritas/ext/GeantImporter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,47 @@ struct ProcessFilter
}
};

//---------------------------------------------------------------------------//
//! Retrieve and store optical material properties, if present.
struct MatPropGetter
{
G4MaterialPropertiesTable const& mpt;

void scalar(double* dst, std::string name, ImportUnits q)
{
if (!mpt.ConstPropertyExists(name))
{
return;
}
*dst = mpt.GetConstProperty(name) * native_value_from_clhep(q);
}

void scalar(double* dst, std::string name, int comp, ImportUnits q)
{
this->scalar(dst, name + std::to_string(comp), q);
}

void vector(ImportPhysicsVector* dst, std::string name, ImportUnits q)
{
auto const* g4vector = mpt.GetProperty(name);
if (!g4vector)
{
return;
}
CELER_ASSERT(g4vector->GetType()
== G4PhysicsVectorType::T_G4PhysicsFreeVector);
double const y_scale = native_value_from_clhep(q);
dst->vector_type = ImportPhysicsVectorType::free;
dst->x.resize(g4vector->GetVectorLength());
dst->y.resize(dst->x.size());
for (auto i : range(dst->x.size()))
{
dst->x[i] = g4vector->Energy(i);
dst->y[i] = (*g4vector)[i] * y_scale;
}
}
};

//---------------------------------------------------------------------------//
/*!
* Safely switch from G4State [G4Material.hh] to ImportMaterialState.
Expand Down Expand Up @@ -327,6 +368,102 @@ std::vector<ImportElement> import_elements()
return elements;
}

//---------------------------------------------------------------------------//
/*!
* Store material-dependent optical properties.
*
* This returns a map of material index to imported optical property data.
*/
ImportData::ImportOpticalMap import_optical()
{
auto const& pct = *G4ProductionCutsTable::GetProductionCutsTable();
auto num_materials = pct.GetTableSize();
CELER_ASSERT(num_materials > 0);

ImportData::ImportOpticalMap result;
for (auto mat_idx : range(num_materials))
{
auto const* mcc = pct.GetMaterialCutsCouple(mat_idx);
CELER_ASSERT(mcc);
CELER_ASSERT(static_cast<std::size_t>(mcc->GetIndex()) == mat_idx);

auto const* material = mcc->GetMaterial();
CELER_ASSERT(material);

// Add optical material properties, if any are present
auto const* mpt = material->GetMaterialPropertiesTable();
if (!mpt)
{
continue;
}
ImportOpticalMaterial optical;
MatPropGetter get_property{*mpt};

// Save common properties
get_property.vector(&optical.properties.refractive_index,
"RINDEX",
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,
"SCINTILLATIONYIELD",
comp_idx,
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.components.push_back(comp);
}
}

// Save Rayleigh properties
get_property.vector(
&optical.rayleigh.mfp, "RAYLEIGH", ImportUnits::len);
get_property.scalar(&optical.rayleigh.scale_factor,
"RS_SCALE_FACTOR",
ImportUnits::unitless);
get_property.scalar(&optical.rayleigh.compressibility,
"ISOTHERMAL_COMPRESSIBILITY",
ImportUnits::len_time_sq_per_mass);

// Save absorption properties
get_property.vector(&optical.absorption.absorption_length,
"ABSLENGTH",
ImportUnits::len);

if (optical)
{
result[mat_idx] = optical;
}
}
CELER_LOG(debug) << "Loaded " << result.size() << " optical materials";
return result;
}

//---------------------------------------------------------------------------//
/*!
* Return a populated \c ImportMaterial vector.
Expand Down Expand Up @@ -749,6 +886,7 @@ ImportData GeantImporter::operator()(DataSelection const& selected)
imported.isotopes = import_isotopes();
imported.elements = import_elements();
imported.materials = import_materials(selected.particles);
imported.optical = import_optical();
}
if (selected.processes != DataSelection::none)
{
Expand Down
6 changes: 6 additions & 0 deletions src/celeritas/ext/RootInterfaceLinkDef.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,12 @@
#pragma link C++ class celeritas::ImportPhysicsTable+;
#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::ImportOpticalRayleigh+;
#pragma link C++ class celeritas::ImportOpticalAbsorption+;
#pragma link C++ class celeritas::ImportOpticalProperty+;
#pragma link C++ class celeritas::ImportOpticalMaterial+;
#pragma link C++ class celeritas::ImportProductionCut+;
#pragma link C++ class celeritas::ImportMatElemComponent+;
#pragma link C++ class celeritas::ImportElement+;
Expand Down
4 changes: 4 additions & 0 deletions src/celeritas/io/ImportData.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "ImportElement.hh"
#include "ImportLivermorePE.hh"
#include "ImportMaterial.hh"
#include "ImportOpticalMaterial.hh"
#include "ImportParameters.hh"
#include "ImportParticle.hh"
#include "ImportProcess.hh"
Expand Down Expand Up @@ -53,9 +54,11 @@ struct ImportData
//!@{
//! \name Type aliases
using ZInt = int;
using MatIdx = int;
using ImportSBMap = std::map<ZInt, ImportSBTable>;
using ImportLivermorePEMap = std::map<ZInt, ImportLivermorePE>;
using ImportAtomicRelaxationMap = std::map<ZInt, ImportAtomicRelaxation>;
using ImportOpticalMap = std::map<MatIdx, ImportOpticalMaterial>;
//!@}

std::vector<ImportParticle> particles;
Expand All @@ -65,6 +68,7 @@ struct ImportData
std::vector<ImportProcess> processes;
std::vector<ImportMscModel> msc_models;
std::vector<ImportVolume> volumes;
ImportOpticalMap optical;
ImportEmParameters em_params;
ImportTransParameters trans_params;
ImportSBMap sb_data;
Expand Down

0 comments on commit 4c469e2

Please sign in to comment.