From 8f5849a5098a8792ef2fdf1ddf014eee901d2f01 Mon Sep 17 00:00:00 2001 From: "Seth R. Johnson" Date: Wed, 31 Jan 2024 20:11:47 -0500 Subject: [PATCH] Support native CLHEP unit system (#1085) --- .github/workflows/build-full.yml | 5 + CMakeLists.txt | 9 + app/celer-g4/DetectorConstruction.cc | 8 +- app/celer-g4/PGPrimaryGeneratorAction.cc | 2 +- app/celer-g4/test-harness.py | 2 +- interface/celeritas.i | 7 + scripts/cmake-presets/ci-ubuntu-cuda.json | 9 + scripts/cmake-presets/goldfinger.json | 9 + scripts/cmake-presets/wildstyle.json | 9 + src/accel/AlongStepFactory.cc | 8 +- src/accel/ExceptionConverter.cc | 2 +- src/accel/LocalTransporter.cc | 4 +- src/accel/RZMapMagneticField.hh | 3 +- src/accel/detail/HitProcessor.cc | 4 +- src/accel/detail/TouchableUpdater.cc | 8 +- src/celeritas/CMakeLists.txt | 4 + src/celeritas/Quantities.hh | 1 + src/celeritas/UnitTypes.hh | 6 + src/celeritas/Units.hh | 6 + src/celeritas/em/data/RayleighData.hh | 10 +- .../em/interactor/RayleighInteractor.hh | 7 - src/celeritas/em/model/RayleighModel.cc | 3 + src/celeritas/em/msc/detail/MscStepToGeo.hh | 4 +- src/celeritas/em/msc/detail/UrbanMscHelper.hh | 2 + src/celeritas/em/xs/LPMCalculator.hh | 2 +- src/celeritas/ext/Convert.geant.hh | 22 ++- src/celeritas/ext/GeantGeoData.hh | 1 + src/celeritas/ext/GeantGeoParams.cc | 4 +- src/celeritas/ext/GeantGeoTrackView.hh | 27 +-- src/celeritas/ext/GeantImporter.cc | 39 ++-- src/celeritas/ext/GeantPhysicsOptions.hh | 2 +- src/celeritas/ext/RootImporter.cc | 3 + src/celeritas/ext/RootImporter.hh | 3 +- src/celeritas/ext/VecgeomParams.cc | 9 +- .../ext/detail/GeantMicroXsCalculator.cc | 43 ++--- .../ext/detail/GeantMicroXsCalculator.hh | 5 +- src/celeritas/ext/detail/GeantPhysicsList.cc | 8 +- .../ext/detail/GeantProcessImporter.cc | 50 +----- src/celeritas/ext/g4vg/Scaler.hh | 4 +- src/celeritas/grid/EnergyLossCalculator.hh | 6 +- src/celeritas/io/EventWriter.cc | 8 +- src/celeritas/io/ImportData.cc | 53 ++++++ src/celeritas/io/ImportData.hh | 15 ++ src/celeritas/io/ImportModel.hh | 2 + src/celeritas/io/ImportParameters.hh | 10 ++ src/celeritas/io/ImportPhysicsTable.cc | 19 -- src/celeritas/io/ImportPhysicsTable.hh | 34 +--- src/celeritas/io/ImportPhysicsVector.cc | 8 +- src/celeritas/io/ImportPhysicsVector.hh | 1 + src/celeritas/io/ImportUnits.cc | 125 +++++++++++++ src/celeritas/io/ImportUnits.hh | 60 +++++++ .../io/detail/ImportDataConverter.cc | 170 ++++++++++++++++++ .../io/detail/ImportDataConverter.hh | 61 +++++++ src/celeritas/mat/ElementSelector.hh | 4 +- src/celeritas/phys/Interaction.hh | 6 +- src/celeritas/phys/ParticleTrackView.hh | 4 +- src/celeritas/phys/PhysicsStepUtils.hh | 2 +- src/celeritas/phys/ProcessBuilder.cc | 4 + test/CMakeLists.txt | 22 ++- test/Test.cc | 6 +- test/accel/HepMC3PrimaryGenerator.test.cc | 2 + test/accel/detail/HitProcessor.test.cc | 30 ++-- test/accel/detail/TouchableUpdater.test.cc | 59 +++--- test/celeritas/AllGeoTypedTestBase.hh | 5 + test/celeritas/GeantTestBase.cc | 1 + test/celeritas/GenericGeoTestBase.cc | 38 ++-- test/celeritas/GenericGeoTestBase.hh | 13 +- test/celeritas/MockTestBase.cc | 16 +- test/celeritas/SimpleTestBase.cc | 14 +- test/celeritas/em/CombinedBrem.test.cc | 2 +- test/celeritas/em/EPlusGG.test.cc | 4 +- test/celeritas/em/Fluctuation.test.cc | 11 +- test/celeritas/em/LivermorePE.test.cc | 5 +- test/celeritas/em/MollerBhabha.test.cc | 2 +- test/celeritas/em/Rayleigh.test.cc | 2 +- test/celeritas/em/RelativisticBrem.test.cc | 6 +- test/celeritas/em/SeltzerBerger.test.cc | 20 ++- test/celeritas/em/UrbanMsc.test.cc | 92 +++++++--- test/celeritas/em/Wentzel.test.cc | 53 +++--- test/celeritas/ext/GeantGeo.test.cc | 77 ++++---- test/celeritas/ext/GeantImporter.test.cc | 57 ++++-- test/celeritas/ext/Vecgeom.test.cc | 154 +++++++++------- test/celeritas/field/CMSParameterizedField.hh | 1 - test/celeritas/field/LinearPropagator.test.cc | 63 ++++--- test/celeritas/field/MagFieldEquation.test.cc | 41 +++-- test/celeritas/geo/CheckedGeoTrackView.hh | 5 +- test/celeritas/geo/GeoMaterial.test.cc | 5 +- test/celeritas/global/AlongStep.test.cc | 16 +- test/celeritas/global/AlongStepTestBase.cc | 13 +- test/celeritas/global/AlongStepTestBase.hh | 4 +- .../global/KernelContextException.test.cc | 7 +- test/celeritas/global/Stepper.test.cc | 5 +- test/celeritas/global/StepperTestBase.cc | 1 + test/celeritas/io/EventIOTestBase.cc | 13 +- test/celeritas/io/ImportUnits.test.cc | 91 ++++++++++ test/celeritas/mat/ElementSelector.test.cc | 8 +- test/celeritas/mat/Material.test.cc | 47 +++-- test/celeritas/mat/Material.test.cu | 10 +- test/celeritas/mat/MaterialTestBase.hh | 9 +- test/celeritas/optical/Cerenkov.test.cc | 18 +- .../optical/ScintillationGenerator.test.cc | 2 +- test/celeritas/phys/CutoffParams.test.cc | 7 +- test/celeritas/phys/InteractorHostTestBase.cc | 10 +- test/celeritas/phys/MockProcess.cc | 33 ++-- test/celeritas/phys/MockProcess.hh | 25 ++- test/celeritas/phys/Particle.test.cc | 9 +- test/celeritas/phys/Particle.test.cu | 7 +- test/celeritas/phys/Physics.test.cc | 46 +++-- test/celeritas/phys/Physics.test.cu | 5 +- test/celeritas/phys/PhysicsStepUtils.test.cc | 55 +++--- test/celeritas/user/Diagnostic.test.cc | 5 +- test/celeritas/user/MctruthTestBase.cc | 6 +- test/celeritas/user/StepCollector.test.cc | 3 +- 113 files changed, 1557 insertions(+), 650 deletions(-) create mode 100644 src/celeritas/io/ImportData.cc create mode 100644 src/celeritas/io/ImportUnits.cc create mode 100644 src/celeritas/io/ImportUnits.hh create mode 100644 src/celeritas/io/detail/ImportDataConverter.cc create mode 100644 src/celeritas/io/detail/ImportDataConverter.hh create mode 100644 test/celeritas/io/ImportUnits.test.cc diff --git a/.github/workflows/build-full.yml b/.github/workflows/build-full.yml index 238461b4f6..9ac6225430 100644 --- a/.github/workflows/build-full.yml +++ b/.github/workflows/build-full.yml @@ -13,6 +13,7 @@ jobs: gpu: name: gpu strategy: + fail-fast: false matrix: special: [null] geometry: ['orange', 'vecgeom'] @@ -46,6 +47,10 @@ jobs: - geometry: 'vecgeom' buildtype: 'reldeb' image: 'ubuntu-cuda' + - special: 'clhep' + geometry: 'vecgeom' + buildtype: 'debug' + image: 'ubuntu-cuda' env: ASAN_OPTIONS: "detect_leaks=0" CELER_TEST_STRICT: 1 diff --git a/CMakeLists.txt b/CMakeLists.txt index d77fbcc7e4..5f11f6dbfa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -136,6 +136,15 @@ celeritas_setup_option(CELERITAS_REAL_TYPE float) celeritas_define_options(CELERITAS_REAL_TYPE "Global runtime precision for real numbers") +if((CELERITAS_CORE_GEO STREQUAL "ORANGE") + AND (NOT CELERITAS_UNITS STREQUAL "CGS")) + celeritas_error_incompatible_option( + "ORANGE currently requires CGS units" + CELERITAS_UNITS + CGS + ) +endif() + #----------------------------------------------------------------------------# # CMAKE VERSION CHECKS #----------------------------------------------------------------------------# diff --git a/app/celer-g4/DetectorConstruction.cc b/app/celer-g4/DetectorConstruction.cc index a1a36fa20e..66a0c31e13 100644 --- a/app/celer-g4/DetectorConstruction.cc +++ b/app/celer-g4/DetectorConstruction.cc @@ -205,7 +205,7 @@ auto DetectorConstruction::construct_field() const -> FieldData if (norm(field_val) > 0) { CELER_LOG_LOCAL(info) - << "Using a uniform field " << field_val << " [tesla]"; + << "Using a uniform field " << field_val << " [T]"; g4field = std::make_shared( convert_to_geant(field_val, CLHEP::tesla)); } @@ -239,9 +239,9 @@ void DetectorConstruction::ConstructSDandField() auto const& field_options = GlobalSetup::Instance()->GetFieldOptions(); auto chord_finder = std::make_unique( mag_field_.get(), - convert_to_geant(field_options.minimum_step, CLHEP::cm)); + convert_to_geant(field_options.minimum_step, clhep_length)); chord_finder->SetDeltaChord( - convert_to_geant(field_options.delta_chord, CLHEP::cm)); + convert_to_geant(field_options.delta_chord, clhep_length)); // Construct the magnetic field G4FieldManager* field_manager @@ -251,7 +251,7 @@ void DetectorConstruction::ConstructSDandField() field_manager->SetChordFinder(chord_finder.release()); field_manager->SetMinimumEpsilonStep(field_options.epsilon_step); field_manager->SetDeltaIntersection( - convert_to_geant(field_options.delta_intersection, CLHEP::cm)); + convert_to_geant(field_options.delta_intersection, clhep_length)); } auto sd_type = GlobalSetup::Instance()->input().sd_type; diff --git a/app/celer-g4/PGPrimaryGeneratorAction.cc b/app/celer-g4/PGPrimaryGeneratorAction.cc index 9322731be7..1869201709 100644 --- a/app/celer-g4/PGPrimaryGeneratorAction.cc +++ b/app/celer-g4/PGPrimaryGeneratorAction.cc @@ -74,7 +74,7 @@ void PGPrimaryGeneratorAction::GeneratePrimaries(G4Event* event) { gun_.SetParticleDefinition(particle_def_[i % particle_def_.size()]); gun_.SetParticlePosition( - convert_to_geant(sample_pos_(rng_), CLHEP::cm)); + convert_to_geant(sample_pos_(rng_), clhep_length)); gun_.SetParticleMomentumDirection( convert_to_geant(sample_dir_(rng_), 1)); gun_.SetParticleEnergy( diff --git a/app/celer-g4/test-harness.py b/app/celer-g4/test-harness.py index 71b40eb553..8dfae72834 100755 --- a/app/celer-g4/test-harness.py +++ b/app/celer-g4/test-harness.py @@ -132,7 +132,7 @@ def strtobool(text): except: pass else: - outfilename = f'{run_name}.out.failed.json' + outfilename = f'{problem_name}.out.failed.json' with open(outfilename, 'w') as f: json.dump(j, f, indent=1) print("Failure written to", outfilename, file=stderr) diff --git a/interface/celeritas.i b/interface/celeritas.i index b79b84b365..787037cb4f 100644 --- a/interface/celeritas.i +++ b/interface/celeritas.i @@ -69,6 +69,11 @@ #include "celeritas/Constants.hh" %} +namespace celeritas +{ +enum class UnitSystem; +} + %include "celeritas/Units.hh" %include "celeritas/Constants.hh" @@ -111,6 +116,8 @@ namespace celeritas %include "celeritas/io/ImportPhysicsVector.hh" %template(VecImportPhysicsVector) std::vector; +%include "celeritas/io/ImportUnits.hh" + %include "celeritas/io/ImportPhysicsTable.hh" %template(VecImportPhysicsTable) std::vector; diff --git a/scripts/cmake-presets/ci-ubuntu-cuda.json b/scripts/cmake-presets/ci-ubuntu-cuda.json index 4862c6e8a8..c5aa55efec 100644 --- a/scripts/cmake-presets/ci-ubuntu-cuda.json +++ b/scripts/cmake-presets/ci-ubuntu-cuda.json @@ -73,6 +73,15 @@ "CELERITAS_USE_OpenMP": {"type": "BOOL", "value": "OFF"} } }, + { + "name": "debug-vecgeom-clhep", + "description": "Build with VecGeom and CLHEP units", + "inherits": [".vecgeom", ".debug", "base"], + "cacheVariables": { + "CELERITAS_USE_CUDA": {"type": "BOOL", "value": "OFF"}, + "CELERITAS_UNITS": "CLHEP" + } + }, { "name": "reldeb-vecgeom", "description": "Build with RelWithDebInfo, assertions, and VecGeom", diff --git a/scripts/cmake-presets/goldfinger.json b/scripts/cmake-presets/goldfinger.json index 3c871338ed..d6580d8037 100644 --- a/scripts/cmake-presets/goldfinger.json +++ b/scripts/cmake-presets/goldfinger.json @@ -58,6 +58,15 @@ "CELERITAS_REAL_TYPE": "float" } }, + { + "name": "clhep", + "displayName": "With CLHEP units", + "inherits": [".base", ".debug", "default"], + "cacheVariables": { + "CELERITAS_UNITS": "CLHEP", + "CELERITAS_USE_VecGeom": {"type": "BOOL", "value": "ON"} + } + }, { "name": "vecgeom", "displayName": "With vecgeom", diff --git a/scripts/cmake-presets/wildstyle.json b/scripts/cmake-presets/wildstyle.json index f4ad2086ba..50e5f0631d 100644 --- a/scripts/cmake-presets/wildstyle.json +++ b/scripts/cmake-presets/wildstyle.json @@ -72,6 +72,15 @@ "CELERITAS_USE_VecGeom": {"type": "BOOL", "value": "ON"} } }, + { + "name": "clhep", + "displayName": "Debug with VecGeom and clhep units", + "inherits": [".debug", ".base"], + "cacheVariables": { + "CELERITAS_UNITS": "CLHEP", + "CELERITAS_USE_VecGeom": {"type": "BOOL", "value": "ON"} + } + }, { "name": "reldeb", "displayName": "Everything but vecgeom in release mode with debug symbols and assertions", diff --git a/src/accel/AlongStepFactory.cc b/src/accel/AlongStepFactory.cc index 88672fe7df..5cc3f36ce4 100644 --- a/src/accel/AlongStepFactory.cc +++ b/src/accel/AlongStepFactory.cc @@ -11,6 +11,7 @@ #include "corecel/io/Logger.hh" #include "corecel/math/ArrayUtils.hh" +#include "corecel/math/QuantityIO.hh" #include "celeritas/em/UrbanMscParams.hh" #include "celeritas/ext/Convert.geant.hh" #include "celeritas/field/RZMapFieldInput.hh" @@ -47,13 +48,14 @@ auto UniformAlongStepFactory::operator()(AlongStepFactoryInput const& input) con { // Get the field strength in tesla (or zero if accessor is undefined) auto field_params = get_field_ ? get_field_() : UniformFieldParams{}; - real_type magnitude_tesla = norm(field_params.field) / units::tesla; + auto magnitude + = native_value_to(norm(field_params.field)); - if (magnitude_tesla > 0) + if (magnitude > zero_quantity()) { // Create a uniform field CELER_LOG(info) << "Creating along-step action with field strength " - << magnitude_tesla << "T"; + << magnitude; return celeritas::AlongStepUniformMscAction::from_params( input.action_id, *input.material, diff --git a/src/accel/ExceptionConverter.cc b/src/accel/ExceptionConverter.cc index 104902a630..34a6c69384 100644 --- a/src/accel/ExceptionConverter.cc +++ b/src/accel/ExceptionConverter.cc @@ -102,7 +102,7 @@ void log_state(Logger::Message& msg, msg << "\n- Particle type ID: " << kce.particle(); } msg << "\n- Energy: " << kce.energy() << "\n- Position: " << kce.pos() - << " (cm)" + << " [" << units::NativeTraits::Length::label() << "]" << "\n- Direction: " << kce.dir(); if (core_params && kce.volume()) diff --git a/src/accel/LocalTransporter.cc b/src/accel/LocalTransporter.cc index af530d460f..db5bfa8a9b 100644 --- a/src/accel/LocalTransporter.cc +++ b/src/accel/LocalTransporter.cc @@ -156,9 +156,9 @@ void LocalTransporter::Push(G4Track const& g4track) << g4track.GetDefinition()->GetParticleName() << "' particles"); - track.position = convert_from_geant(g4track.GetPosition(), CLHEP::cm); + track.position = convert_from_geant(g4track.GetPosition(), clhep_length); track.direction = convert_from_geant(g4track.GetMomentumDirection(), 1); - track.time = convert_from_geant(g4track.GetGlobalTime(), CLHEP::s); + track.time = convert_from_geant(g4track.GetGlobalTime(), clhep_time); // TODO: Celeritas track IDs are independent from Geant4 track IDs, since // they must be sequential from zero for a given event. We may need to save diff --git a/src/accel/RZMapMagneticField.hh b/src/accel/RZMapMagneticField.hh index 289db8ea7a..43198c94b3 100644 --- a/src/accel/RZMapMagneticField.hh +++ b/src/accel/RZMapMagneticField.hh @@ -62,8 +62,7 @@ RZMapMagneticField::RZMapMagneticField(SPConstFieldParams params) void RZMapMagneticField::GetFieldValue(double const pos[3], double* field) const { // Calculate the magnetic field value in the native Celeritas unit system - Real3 result - = calc_field_(convert_from_geant(pos, CLHEP::cm) * units::centimeter); + Real3 result = calc_field_(convert_from_geant(pos, clhep_length)); for (auto i = 0; i < 3; ++i) { // Return values of the field vector in CLHEP::tesla for Geant4 diff --git a/src/accel/detail/HitProcessor.cc b/src/accel/detail/HitProcessor.cc index bd63b57897..0728f51eb7 100644 --- a/src/accel/detail/HitProcessor.cc +++ b/src/accel/detail/HitProcessor.cc @@ -196,8 +196,8 @@ void HitProcessor::operator()(DetectorStepOutput const& out) const { continue; } - HP_SET(points[sp]->SetGlobalTime, out.points[sp].time, CLHEP::s); - HP_SET(points[sp]->SetPosition, out.points[sp].pos, CLHEP::cm); + HP_SET(points[sp]->SetGlobalTime, out.points[sp].time, clhep_time); + HP_SET(points[sp]->SetPosition, out.points[sp].pos, clhep_length); HP_SET(points[sp]->SetKineticEnergy, out.points[sp].energy, CLHEP::MeV); diff --git a/src/accel/detail/TouchableUpdater.cc b/src/accel/detail/TouchableUpdater.cc index b859fd4e56..5753da23d1 100644 --- a/src/accel/detail/TouchableUpdater.cc +++ b/src/accel/detail/TouchableUpdater.cc @@ -57,7 +57,7 @@ bool TouchableUpdater::operator()(Real3 const& pos, Real3 const& dir, G4LogicalVolume const* lv) { - auto g4pos = convert_to_geant(pos, CLHEP::cm); + auto g4pos = convert_to_geant(pos, clhep_length); auto g4dir = convert_to_geant(dir, 1); // Locate pre-step point @@ -74,9 +74,9 @@ bool TouchableUpdater::operator()(Real3 const& pos, return true; } - constexpr double g4max_step = convert_to_geant(max_step(), CLHEP::cm); + constexpr double g4max_step = convert_to_geant(max_step(), clhep_length); constexpr double g4max_quiet_step - = convert_to_geant(max_quiet_step(), CLHEP::cm); + = convert_to_geant(max_quiet_step(), clhep_length); double g4safety{-1}; double g4step{-1}; @@ -157,7 +157,7 @@ bool TouchableUpdater::operator()(Real3 const& pos, } // Reset the position and flip the direction - g4pos = convert_to_geant(pos, CLHEP::cm); + g4pos = convert_to_geant(pos, clhep_length); g4dir *= -1; find_next_step(); if (try_cross_boundary()) diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index dbffe585a8..0485797260 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -46,13 +46,16 @@ list(APPEND SOURCES grid/ValueGridInserter.cc grid/VectorUtils.cc io/AtomicRelaxationReader.cc + io/ImportData.cc io/ImportMaterial.cc io/ImportModel.cc io/ImportPhysicsTable.cc io/ImportPhysicsVector.cc io/ImportProcess.cc + io/ImportUnits.cc io/LivermorePEReader.cc io/SeltzerBergerReader.cc + io/detail/ImportDataConverter.cc mat/MaterialParams.cc mat/MaterialParamsOutput.cc mat/detail/Utils.cc @@ -194,6 +197,7 @@ if(CELERITAS_USE_ROOT) "celeritas/io/ImportPhysicsTable.hh" "celeritas/io/ImportPhysicsVector.hh" "celeritas/io/ImportProcess.hh" + "celeritas/io/ImportUnits.hh" "celeritas/io/ImportVolume.hh" "celeritas/io/EventData.hh" NOINSTALL diff --git a/src/celeritas/Quantities.hh b/src/celeritas/Quantities.hh index 5a43de9ce6..85064c81ef 100644 --- a/src/celeritas/Quantities.hh +++ b/src/celeritas/Quantities.hh @@ -37,6 +37,7 @@ using CmLength = Quantity; using InvCmXs = Quantity>; using InvCcDensity = Quantity; using MolCcDensity = Quantity; +using GramCcDensity = Quantity; using FieldTesla = Quantity; //!@} //---------------------------------------------------------------------------// diff --git a/src/celeritas/UnitTypes.hh b/src/celeritas/UnitTypes.hh index 2d6fc27be9..4028a272ef 100644 --- a/src/celeritas/UnitTypes.hh +++ b/src/celeritas/UnitTypes.hh @@ -151,6 +151,12 @@ struct MolPerCentimeterCubed : UnitProduct static char const* label() { return "mol/cm^3"; } }; +//! Mass density +struct GramPerCentimeterCubed : UnitProduct +{ + static char const* label() { return "g/cm^3"; } +}; + //!@} //---------------------------------------------------------------------------// //!@{ diff --git a/src/celeritas/Units.hh b/src/celeritas/Units.hh index 8eb9af4284..1cdf9d8bee 100644 --- a/src/celeritas/Units.hh +++ b/src/celeritas/Units.hh @@ -43,6 +43,12 @@ namespace units * Additionally: * - radians are used for measures of angle (unitless) * - steradians are used for measures of solid angle (unitless) + * + * TODO: if we're serious about supporting single-precision arithmetic, we + * should define a helper class that stores the constant as full precision but + * when multiplied by a single/double is truncated to that precision. + * Otherwise, if \c real_type is single-precision, then we lose accuracy in + * places like the GeantImporter where everything is double precision. */ #define CELER_ICRT inline constexpr real_type diff --git a/src/celeritas/em/data/RayleighData.hh b/src/celeritas/em/data/RayleighData.hh index a0d05dab60..4028221fff 100644 --- a/src/celeritas/em/data/RayleighData.hh +++ b/src/celeritas/em/data/RayleighData.hh @@ -23,9 +23,9 @@ namespace celeritas * \f[ * FF(E,cos)^2 = \Sigma_{j} \frac{a_j}{[1 + b_j x]^{n}} * \f] - * where \f$ x= E^{2}(1-cos\theta) \f$ and \em n is the high energy slope of - * the form factor and \em a and \em b are free parameters to obtain the best - * fit to the form factor. The unit for the energy (\em E) is in MeV. + * where \f$ x = E^{2}(1 - \cos\theta) \f$ and \em n is the high energy slope + * of the form factor and \em a and \em b are free parameters to obtain the + * best fit to the form factor. The unit for the energy (\em E) is in MeV. */ struct RayleighParameters { @@ -79,8 +79,6 @@ struct RayleighData } }; -using RayleighDeviceRef = DeviceCRef; -using RayleighHostRef = HostCRef; using RayleighRef = NativeCRef; - +//---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/em/interactor/RayleighInteractor.hh b/src/celeritas/em/interactor/RayleighInteractor.hh index 4481cca2ee..87d4de7b2c 100644 --- a/src/celeritas/em/interactor/RayleighInteractor.hh +++ b/src/celeritas/em/interactor/RayleighInteractor.hh @@ -58,13 +58,6 @@ class RayleighInteractor //// CONSTANTS //// - //! cm/hc in the MeV energy unit - static CELER_CONSTEXPR_FUNCTION real_type hc_factor() - { - return units::centimeter * native_value_from(units::MevEnergy{1.0}) - / (constants::c_light * constants::h_planck); - } - //! A point where the functional form of the form factor fit changes static CELER_CONSTEXPR_FUNCTION real_type fit_slice() { return 0.02; } diff --git a/src/celeritas/em/model/RayleighModel.cc b/src/celeritas/em/model/RayleighModel.cc index 7613c97d53..dc7159b1ce 100644 --- a/src/celeritas/em/model/RayleighModel.cc +++ b/src/celeritas/em/model/RayleighModel.cc @@ -153,6 +153,9 @@ void RayleighModel::build_data(HostValue* data, MaterialParams const& materials) * Parameters for Z = 0 are dropped as they are zeros and not used. * Reshaped as [el][param][3] with params.T.reshape((100, 3, 3)), * then updated 'n' with params[:,2,:] -= 1 + * + * The fit data embeds centimeters as a unit system: this is accounted for in + * the interactor. */ auto RayleighModel::get_el_parameters(AtomicNumber z) -> ElScatParams const& { diff --git a/src/celeritas/em/msc/detail/MscStepToGeo.hh b/src/celeritas/em/msc/detail/MscStepToGeo.hh index 87099ed671..d0615238cd 100644 --- a/src/celeritas/em/msc/detail/MscStepToGeo.hh +++ b/src/celeritas/em/msc/detail/MscStepToGeo.hh @@ -77,8 +77,8 @@ class MscStepToGeo struct result_type { - real_type step{}; //!< Geometrical step length - real_type alpha{0}; //!< Scaled MFP slope + real_type step{}; //!< Geometrical step length [len] + real_type alpha{0}; //!< Scaled MFP slope [1/len] }; public: diff --git a/src/celeritas/em/msc/detail/UrbanMscHelper.hh b/src/celeritas/em/msc/detail/UrbanMscHelper.hh index d5f0ae4226..dac7bf4342 100644 --- a/src/celeritas/em/msc/detail/UrbanMscHelper.hh +++ b/src/celeritas/em/msc/detail/UrbanMscHelper.hh @@ -30,6 +30,8 @@ namespace detail /*! * This is a helper class for the UrbanMscStepLimit and UrbanMscScatter. * + * NOTE: units are "native" units, listed here as CGS. + * * \todo Refactor to UrbanMscTrackView . */ class UrbanMscHelper diff --git a/src/celeritas/em/xs/LPMCalculator.hh b/src/celeritas/em/xs/LPMCalculator.hh index 6f7902ced4..1769eda578 100644 --- a/src/celeritas/em/xs/LPMCalculator.hh +++ b/src/celeritas/em/xs/LPMCalculator.hh @@ -62,7 +62,7 @@ class LPMCalculator // Current element ElementView const& element_; - // Electron density of the current material [1/cm^3] + // Electron density of the current material [1/len^3] real_type const electron_density_; // Characteristic energy for the LPM effect for this material [MeV] real_type const lpm_energy_; diff --git a/src/celeritas/ext/Convert.geant.hh b/src/celeritas/ext/Convert.geant.hh index 17b26eaff3..d869f0c491 100644 --- a/src/celeritas/ext/Convert.geant.hh +++ b/src/celeritas/ext/Convert.geant.hh @@ -7,17 +7,31 @@ //---------------------------------------------------------------------------// #pragma once -#include #include #include "corecel/Types.hh" #include "corecel/cont/Array.hh" #include "orange/Types.hh" #include "celeritas/Quantities.hh" +#include "celeritas/UnitTypes.hh" namespace celeritas { //---------------------------------------------------------------------------// +// CONSTANTS +//---------------------------------------------------------------------------// +//! Value of a unit CLHEP length in the native Celeritas system +inline constexpr real_type clhep_length + = 1 / units::ClhepTraits::Length::value(); +//! Value of a unit CLHEP field in the native Celeritas system +inline constexpr real_type clhep_field = 1 + / units::ClhepTraits::BField::value(); +//! Value of a unit CLHEP time in the native Celeritas system +inline constexpr real_type clhep_time = 1 / units::ClhepTraits::Time::value(); + +//---------------------------------------------------------------------------// +// FREE FUNCTIONS +//---------------------------------------------------------------------------// /*! * Convert a value from Geant4/CLHEP to Celeritas native units. */ @@ -90,12 +104,14 @@ inline G4ThreeVector convert_to_geant(Array const& arr, double units) //---------------------------------------------------------------------------// /*! * Convert Celeritas energy quantities to Geant4. + * + * The unit value should always be CLHEP::MeV which is defined to be unity. */ inline constexpr double convert_to_geant(units::MevEnergy const& energy, double units) { - CELER_EXPECT(units == CLHEP::MeV); - return energy.value() * CLHEP::MeV; + CELER_EXPECT(units == 1); + return energy.value(); } //---------------------------------------------------------------------------// diff --git a/src/celeritas/ext/GeantGeoData.hh b/src/celeritas/ext/GeantGeoData.hh index 46d5a26da8..9a5fe1df6b 100644 --- a/src/celeritas/ext/GeantGeoData.hh +++ b/src/celeritas/ext/GeantGeoData.hh @@ -14,6 +14,7 @@ #include "corecel/data/CollectionBuilder.hh" #include "corecel/sys/ThreadId.hh" #include "celeritas/Types.hh" +#include "celeritas/UnitTypes.hh" #include "detail/GeantGeoNavCollection.hh" diff --git a/src/celeritas/ext/GeantGeoParams.cc b/src/celeritas/ext/GeantGeoParams.cc index 800cdd5430..c3d68d1ba4 100644 --- a/src/celeritas/ext/GeantGeoParams.cc +++ b/src/celeritas/ext/GeantGeoParams.cc @@ -267,11 +267,11 @@ void GeantGeoParams::build_metadata() return BBox({convert_from_geant(G4ThreeVector(extent.GetXmin(), extent.GetYmin(), extent.GetZmin()), - CLHEP::cm), + clhep_length), convert_from_geant(G4ThreeVector(extent.GetXmax(), extent.GetYmax(), extent.GetZmax()), - CLHEP::cm)}); + clhep_length)}); }(); } diff --git a/src/celeritas/ext/GeantGeoTrackView.hh b/src/celeritas/ext/GeantGeoTrackView.hh index 204393a9e3..2affebaceb 100644 --- a/src/celeritas/ext/GeantGeoTrackView.hh +++ b/src/celeritas/ext/GeantGeoTrackView.hh @@ -72,7 +72,10 @@ class GeantGeoTrackView //// STATIC ACCESSORS //// //! A tiny push to make sure tracks do not get stuck at boundaries - static constexpr real_type extra_push() { return 1e-13; } + static constexpr real_type extra_push() + { + return 1e-13 * units::centimeter; + } //// ACCESSORS //// @@ -169,9 +172,9 @@ GeantGeoTrackView::GeantGeoTrackView(ParamsRef const&, , touch_handle_(states.nav_state.touch_handle(tid)) , navi_(states.nav_state.navigator(tid)) { - g4pos_ = convert_to_geant(pos_, CLHEP::cm); + g4pos_ = convert_to_geant(pos_, clhep_length); g4dir_ = convert_to_geant(dir_, 1); - g4safety_ = convert_to_geant(safety_radius_, CLHEP::cm); + g4safety_ = convert_to_geant(safety_radius_, clhep_length); } //---------------------------------------------------------------------------// @@ -188,7 +191,7 @@ GeantGeoTrackView& GeantGeoTrackView::operator=(Initializer_t const& init) next_step_ = 0; safety_radius_ = -1; // Assume *not* on a boundary - g4pos_ = convert_to_geant(pos_, CLHEP::cm); + g4pos_ = convert_to_geant(pos_, clhep_length); g4dir_ = convert_to_geant(dir_, 1); g4safety_ = -1; @@ -302,20 +305,20 @@ Propagation GeantGeoTrackView::find_next_step(real_type max_step) CELER_EXPECT(max_step > 0); // Compute the step - real_type g4step = convert_to_geant(max_step, CLHEP::cm); + real_type g4step = convert_to_geant(max_step, clhep_length); g4step = navi_.ComputeStep(g4pos_, g4dir_, g4step, g4safety_); if (g4safety_ != 0 && !this->is_on_boundary()) { // Save the resulting safety distance if computed: allow to be // "negative" to prevent accidentally changing the boundary state - safety_radius_ = convert_from_geant(g4safety_, CLHEP::cm); + safety_radius_ = convert_from_geant(g4safety_, clhep_length); CELER_ASSERT(!this->is_on_boundary()); } // Update result Propagation result; - result.distance = convert_from_geant(g4step, CLHEP::cm); + result.distance = convert_from_geant(g4step, clhep_length); if (result.distance <= max_step) { result.boundary = true; @@ -362,9 +365,9 @@ auto GeantGeoTrackView::find_safety(real_type max_step) -> real_type CELER_EXPECT(max_step > 0); if (!this->is_on_boundary() && (safety_radius_ < max_step)) { - real_type g4step = convert_to_geant(max_step, CLHEP::cm); + real_type g4step = convert_to_geant(max_step, clhep_length); g4safety_ = navi_.ComputeSafety(g4pos_, g4step); - safety_radius_ = max(convert_from_geant(g4safety_, CLHEP::cm), 0.0); + safety_radius_ = max(convert_from_geant(g4safety_, clhep_length), 0.0); } return safety_radius_; @@ -380,7 +383,7 @@ void GeantGeoTrackView::move_to_boundary() // Move next step axpy(next_step_, dir_, &pos_); - axpy(convert_to_geant(next_step_, CLHEP::cm), g4dir_, &g4pos_); + axpy(convert_to_geant(next_step_, clhep_length), g4dir_, &g4pos_); next_step_ = 0; safety_radius_ = 0; g4safety_ = 0; @@ -422,7 +425,7 @@ void GeantGeoTrackView::move_internal(real_type dist) // Move and update next_step axpy(dist, dir_, &pos_); - axpy(convert_to_geant(dist, CLHEP::cm), g4dir_, &g4pos_); + axpy(convert_to_geant(dist, clhep_length), g4dir_, &g4pos_); next_step_ -= dist; navi_.LocateGlobalPointWithinVolume(g4pos_); @@ -440,7 +443,7 @@ void GeantGeoTrackView::move_internal(real_type dist) void GeantGeoTrackView::move_internal(Real3 const& pos) { pos_ = pos; - g4pos_ = convert_to_geant(pos_, CLHEP::cm); + g4pos_ = convert_to_geant(pos_, clhep_length); next_step_ = 0; navi_.LocateGlobalPointWithinVolume(g4pos_); diff --git a/src/celeritas/ext/GeantImporter.cc b/src/celeritas/ext/GeantImporter.cc index 7a3a2ffb7a..bd95669c65 100644 --- a/src/celeritas/ext/GeantImporter.cc +++ b/src/celeritas/ext/GeantImporter.cc @@ -50,6 +50,7 @@ #include #include +#include "celeritas_config.h" #include "corecel/Assert.hh" #include "corecel/cont/Range.hh" #include "corecel/io/Logger.hh" @@ -70,12 +71,7 @@ #include "detail/GeantProcessImporter.hh" #include "detail/GeantVolumeVisitor.hh" -using CLHEP::cm; -using CLHEP::cm2; -using CLHEP::cm3; -using CLHEP::g; -using CLHEP::MeV; -using CLHEP::s; +inline constexpr double mev_scale = 1 / CLHEP::MeV; namespace celeritas { @@ -207,6 +203,8 @@ import_particles(GeantImporter::DataSelection::Flags particle_flags) std::vector particles; + double const time_scale = native_value_from_clhep(ImportUnits::time); + ParticleFilter include_particle{particle_flags}; while (particle_iterator()) { @@ -231,7 +229,7 @@ import_particles(GeantImporter::DataSelection::Flags particle_flags) if (!particle.is_stable) { // Convert lifetime of unstable particles to seconds - particle.lifetime /= s; + particle.lifetime *= time_scale; } particles.push_back(particle); @@ -331,10 +329,10 @@ std::vector import_elements() /*! * Return a populated \c ImportMaterial vector. * - * TODO: there seems to be an inconsitency between "materials" (index in the + * TODO: there is an inconsitency between "materials" (index in the * global material table) and "material cut couple" (which is what we're - * defining here?) Maybe we need another level of indirection for material and - * material+cutoff values? + * defining here?). We need another level of indirection for material and + * material+cutoff ("physics material"). */ std::vector import_materials(GeantImporter::DataSelection::Flags particle_flags) @@ -383,6 +381,10 @@ import_materials(GeantImporter::DataSelection::Flags particle_flags) cut_converters.emplace_back(gi, std::move(converter)); } + double const numdens_scale + = native_value_from_clhep(ImportUnits::inv_len_cb); + double const len_scale = native_value_from_clhep(ImportUnits::len); + // Loop over material data for (auto i : range(materials.size())) { @@ -406,7 +408,7 @@ import_materials(GeantImporter::DataSelection::Flags particle_flags) material.state = to_material_state(g4material->GetState()); material.temperature = g4material->GetTemperature(); // [K] material.number_density = g4material->GetTotNbOfAtomsPerVolume() - / (1. / cm3); + * numdens_scale; // Populate material production cut values for (auto const& idx_convert : cut_converters) @@ -418,8 +420,8 @@ import_materials(GeantImporter::DataSelection::Flags particle_flags) double const energy = converter.Convert(range, g4material); ImportProductionCut cutoffs; - cutoffs.energy = energy / MeV; - cutoffs.range = range / cm; + cutoffs.energy = energy * mev_scale; + cutoffs.range = range * len_scale; material.pdg_cutoffs.insert({to_pdg(g4i).get(), cutoffs}); } @@ -433,7 +435,7 @@ import_materials(GeantImporter::DataSelection::Flags particle_flags) ImportMatElemComponent elem_comp; elem_comp.element_id = g4element->GetIndex(); double elem_num_density = g4material->GetVecNbOfAtomsPerVolume()[j] - / (1. / cm3); + * numdens_scale; elem_comp.number_fraction = elem_num_density / material.number_density; @@ -637,7 +639,8 @@ import_trans_parameters(GeantImporter::DataSelection::Flags particle_flags) // Get the threshold values for killing looping tracks ImportLoopingThreshold looping; looping.threshold_trials = trans->GetThresholdTrials(); - looping.important_energy = trans->GetThresholdImportantEnergy() / MeV; + looping.important_energy = trans->GetThresholdImportantEnergy() + * mev_scale; CELER_ASSERT(looping); result.looping.insert({particle->GetPDGEncoding(), looping}); } @@ -655,17 +658,18 @@ ImportEmParameters import_em_parameters() ImportEmParameters import; auto const& g4 = *G4EmParameters::Instance(); + double const len_scale = native_value_from_clhep(ImportUnits::len); import.energy_loss_fluct = g4.LossFluctuation(); import.lpm = g4.LPM(); import.integral_approach = g4.Integral(); import.linear_loss_limit = g4.LinearLossLimit(); - import.lowest_electron_energy = g4.LowestElectronEnergy() / MeV; + import.lowest_electron_energy = g4.LowestElectronEnergy() * mev_scale; import.auger = g4.Auger(); import.msc_range_factor = g4.MscRangeFactor(); #if G4VERSION_NUMBER >= 1060 import.msc_safety_factor = g4.MscSafetyFactor(); - import.msc_lambda_limit = g4.MscLambdaLimit() / cm; + import.msc_lambda_limit = g4.MscLambdaLimit() * len_scale; #endif import.apply_cuts = g4.ApplyCuts(); import.screening_factor = g4.ScreeningFactor(); @@ -799,6 +803,7 @@ ImportData GeantImporter::operator()(DataSelection const& selected) } } + imported.units = units::NativeTraits::label(); return imported; } diff --git a/src/celeritas/ext/GeantPhysicsOptions.hh b/src/celeritas/ext/GeantPhysicsOptions.hh index 183a528d86..d9596e9cce 100644 --- a/src/celeritas/ext/GeantPhysicsOptions.hh +++ b/src/celeritas/ext/GeantPhysicsOptions.hh @@ -112,7 +112,7 @@ struct GeantPhysicsOptions //! Kill secondaries below the production cut bool apply_cuts{false}; //! Set the default production cut for all particle types [len] - double default_cutoff{0.1}; + double default_cutoff{0.1 * units::centimeter}; //!@} //!@{ diff --git a/src/celeritas/ext/RootImporter.cc b/src/celeritas/ext/RootImporter.cc index 9deec87470..4df979feec 100644 --- a/src/celeritas/ext/RootImporter.cc +++ b/src/celeritas/ext/RootImporter.cc @@ -61,6 +61,9 @@ ImportData RootImporter::operator()() CELER_ASSERT(err_code >= 0); tree_data->GetEntry(0); + // Convert (if necessary) the resulting data to the native unit system + convert_to_native(&import_data); + return import_data; } diff --git a/src/celeritas/ext/RootImporter.hh b/src/celeritas/ext/RootImporter.hh index 8c77f233d6..6fd2cdf901 100644 --- a/src/celeritas/ext/RootImporter.hh +++ b/src/celeritas/ext/RootImporter.hh @@ -26,7 +26,8 @@ namespace celeritas * * RootImporter loads particle, element, material, process, and volume * information from a ROOT file that contains an \c ImportData object. - * Currently, said ROOT file is created by the \c RootExporter class. + * Currently, said ROOT file is created by the \c RootExporter class. The + * imported data will be converted to the native unit system. * * \c RootImporter , along with all \c Import[Class] type of classes, are the * link between Geant4 and Celeritas. Every Celeritas' host/device class that diff --git a/src/celeritas/ext/VecgeomParams.cc b/src/celeritas/ext/VecgeomParams.cc index 5f71553fc5..f4df3e85ff 100644 --- a/src/celeritas/ext/VecgeomParams.cc +++ b/src/celeritas/ext/VecgeomParams.cc @@ -41,10 +41,11 @@ #include "corecel/sys/ScopedLimitSaver.hh" #include "corecel/sys/ScopedMem.hh" #include "corecel/sys/ScopedProfiling.hh" -#include "celeritas/ext/detail/VecgeomCompatibility.hh" +#include "celeritas/Units.hh" #include "GeantGeoUtils.hh" #include "VecgeomData.hh" // IWYU pragma: associated +#include "detail/VecgeomCompatibility.hh" #include "g4vg/Converter.hh" #ifdef VECGEOM_USE_SURF @@ -291,8 +292,12 @@ void VecgeomParams::build_volumes_vgdml(std::string const& filename) ScopedProfiling profile_this{"load-vecgeom"}; ScopedMem record_mem("VecgeomParams.load_geant_geometry"); ScopedTimeAndRedirect time_and_output_("vgdml::Frontend"); + #ifdef VECGEOM_GDML - vgdml::Frontend::Load(filename, /* validate_xml_schema = */ false); + vgdml::Frontend::Load(filename, + /* validate_xml_schema = */ false, + /* mm_unit = */ units::millimeter, + /* verbose = */ vecgeom_verbosity()); #else CELER_DISCARD(filename); CELER_NOT_CONFIGURED("VGDML"); diff --git a/src/celeritas/ext/detail/GeantMicroXsCalculator.cc b/src/celeritas/ext/detail/GeantMicroXsCalculator.cc index 99e979408d..f2fb9acffe 100644 --- a/src/celeritas/ext/detail/GeantMicroXsCalculator.cc +++ b/src/celeritas/ext/detail/GeantMicroXsCalculator.cc @@ -7,11 +7,12 @@ //---------------------------------------------------------------------------// #include "GeantMicroXsCalculator.hh" -#include #include #include #include "corecel/cont/Range.hh" +#include "corecel/math/Algorithms.hh" +#include "celeritas/io/ImportUnits.hh" namespace celeritas { @@ -54,17 +55,31 @@ void GeantMicroXsCalculator::operator()(VecDouble const& energy_grid, mxs_vec.resize(energy_grid.size()); } + auto calc_element_xs + = [this, &elements](std::size_t elcomp_idx, double energy) { + // Calculate microscopic cross section + double xs = model_.ComputeCrossSectionPerAtom( + &particle_, + elements[elcomp_idx], + energy, + secondary_cut_, + /* max_energy = */ std::numeric_limits::max()); + return clamp_to_nonneg(xs); + }; + double const xs_scaling = native_value_from_clhep(ImportUnits::len_sq); + // Outer loop over energy to reduce material setup calls for (auto energy_idx : range(energy_grid.size())) { double const energy = energy_grid[energy_idx]; + CELER_ASSERT(energy > 0); model_.SetupForMaterial(&particle_, &material_, energy); // Inner loop over elements for (auto elcomp_idx : range(elements.size())) { (*result_xs)[elcomp_idx][energy_idx] - = this->calc_element_xs(*elements[elcomp_idx], energy); + = calc_element_xs(elcomp_idx, energy) * xs_scaling; } } @@ -74,13 +89,13 @@ void GeantMicroXsCalculator::operator()(VecDouble const& energy_grid, // Geant4 simply uses the next/previous bin value when the vector's // front/back are zero. This probably isn't correct but it replicates // Geant4's behavior. - if (xs[0] == 0.0) + if (xs[0] == 0) { xs[0] = xs[1]; } auto const last_idx = xs.size() - 1; - if (xs[last_idx] == 0.0) + if (xs[last_idx] == 0) { // Cross-section ends with zero, use previous bin value xs[last_idx] = xs[last_idx - 1]; @@ -88,26 +103,6 @@ void GeantMicroXsCalculator::operator()(VecDouble const& energy_grid, } } -//---------------------------------------------------------------------------// -/*! - * Calculate after setting up for material and energy. - */ -double GeantMicroXsCalculator::calc_element_xs(G4Element const& el, - double energy) const -{ - CELER_EXPECT(energy >= 0); - // Calculate microscopic cross-section - double xs = model_.ComputeCrossSectionPerAtom( - &particle_, - &el, - energy, - secondary_cut_, - /* max_energy = */ std::numeric_limits::max()); - - // Convert to celeritas units and clamp to zero - return std::max(0, xs / CLHEP::cm2); -} - //---------------------------------------------------------------------------// } // namespace detail } // namespace celeritas diff --git a/src/celeritas/ext/detail/GeantMicroXsCalculator.hh b/src/celeritas/ext/detail/GeantMicroXsCalculator.hh index b88499d77f..044a526eba 100644 --- a/src/celeritas/ext/detail/GeantMicroXsCalculator.hh +++ b/src/celeritas/ext/detail/GeantMicroXsCalculator.hh @@ -33,7 +33,7 @@ class GeantMicroXsCalculator //!@{ //! \name Type aliases using EnergyUnits = units::Mev; - using XsUnits = units::Native; // cm^2 + using XsUnits = units::Native; // len^2 using VecDouble = std::vector; using VecVecDouble = std::vector>; //!@} @@ -52,9 +52,6 @@ class GeantMicroXsCalculator G4ParticleDefinition const& particle_; G4Material const& material_; double secondary_cut_; - - // Calculate after setting up for material and energy - double calc_element_xs(G4Element const& g4el, double energy) const; }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/ext/detail/GeantPhysicsList.cc b/src/celeritas/ext/detail/GeantPhysicsList.cc index d8e3ad1269..624ea71a1f 100644 --- a/src/celeritas/ext/detail/GeantPhysicsList.cc +++ b/src/celeritas/ext/detail/GeantPhysicsList.cc @@ -59,6 +59,8 @@ GeantPhysicsList::GeantPhysicsList(Options const& options) : options_(options) << "number of EM bins per decade=" << options.em_bins_per_decade << " (must be at least 5)"); + using ClhepLen = Quantity; + em_parameters.SetNumberOfBinsPerDecade(options.em_bins_per_decade); em_parameters.SetLossFluctuations(options.eloss_fluctuation); em_parameters.SetMinEnergy(value_as(options.min_energy) @@ -75,7 +77,8 @@ GeantPhysicsList::GeantPhysicsList(Options const& options) : options_(options) // Customizable MSC safety factor/lambda limit were added in // emutils-V10-05-18 em_parameters.SetMscSafetyFactor(options.msc_safety_factor); - em_parameters.SetMscLambdaLimit(options.msc_lambda_limit * CLHEP::cm); + em_parameters.SetMscLambdaLimit( + native_value_to(options.msc_lambda_limit).value()); #endif em_parameters.SetLowestElectronEnergy( value_as(options.lowest_electron_energy) @@ -86,7 +89,8 @@ GeantPhysicsList::GeantPhysicsList(Options const& options) : options_(options) em_parameters.SetMscEnergyLimit(100 * CLHEP::TeV); } em_parameters.SetApplyCuts(options.apply_cuts); - this->SetDefaultCutValue(options.default_cutoff * CLHEP::cm); + this->SetDefaultCutValue( + native_value_to(options.default_cutoff).value()); int verb = options_.verbose ? 1 : 0; this->SetVerboseLevel(verb); diff --git a/src/celeritas/ext/detail/GeantProcessImporter.cc b/src/celeritas/ext/detail/GeantProcessImporter.cc index 0d4c929f86..e095d629bb 100644 --- a/src/celeritas/ext/detail/GeantProcessImporter.cc +++ b/src/celeritas/ext/detail/GeantProcessImporter.cc @@ -32,14 +32,12 @@ #include "corecel/Assert.hh" #include "corecel/cont/Range.hh" #include "corecel/io/Logger.hh" +#include "celeritas/UnitTypes.hh" +#include "celeritas/io/ImportUnits.hh" #include "celeritas/phys/PDGNumber.hh" #include "GeantModelImporter.hh" -using CLHEP::cm; -using CLHEP::cm2; -using CLHEP::MeV; - namespace celeritas { namespace detail @@ -142,36 +140,6 @@ int get_secondary_pdg(T const& process) return 0; } -//---------------------------------------------------------------------------// -/*! - * Get a multiplicative geant4-natural-units constant to convert the units. - */ -double units_to_scaling(ImportUnits units) -{ - switch (units) - { - case ImportUnits::none: - return 1; - case ImportUnits::mev: - return 1 / MeV; - case ImportUnits::mev_per_cm: - return cm / MeV; - case ImportUnits::cm: - return 1 / cm; - case ImportUnits::cm_inv: - return cm; - case ImportUnits::cm_mev_inv: - return cm * MeV; - case ImportUnits::mev_2_per_cm: - return cm / (MeV * MeV); - case ImportUnits::cm_2: - return 1 / (cm * cm); - case ImportUnits::size_: - CELER_ASSERT_UNREACHABLE(); - } - CELER_ASSERT_UNREACHABLE(); -} - //---------------------------------------------------------------------------// /*! * Convert physics vector type from Geant4 to Celeritas IO. @@ -225,31 +193,31 @@ void append_table(G4PhysicsTable const* g4table, { case ImportTableType::dedx: table.x_units = ImportUnits::mev; - table.y_units = ImportUnits::mev_per_cm; + table.y_units = ImportUnits::mev_per_len; break; case ImportTableType::range: table.x_units = ImportUnits::mev; - table.y_units = ImportUnits::cm; + table.y_units = ImportUnits::len; break; case ImportTableType::lambda: table.x_units = ImportUnits::mev; - table.y_units = ImportUnits::cm_inv; + table.y_units = ImportUnits::len_inv; break; case ImportTableType::lambda_prim: table.x_units = ImportUnits::mev; - table.y_units = ImportUnits::cm_mev_inv; + table.y_units = ImportUnits::len_mev_inv; break; case ImportTableType::msc_xs: table.x_units = ImportUnits::mev; - table.y_units = ImportUnits::mev_2_per_cm; + table.y_units = ImportUnits::mev_sq_per_len; break; default: CELER_ASSERT_UNREACHABLE(); }; // Convert units - double x_scaling = units_to_scaling(table.x_units); - double y_scaling = units_to_scaling(table.y_units); + double const x_scaling = native_value_from_clhep(table.x_units); + double const y_scaling = native_value_from_clhep(table.y_units); // Save physics vectors for (auto const* g4vector : *g4table) diff --git a/src/celeritas/ext/g4vg/Scaler.hh b/src/celeritas/ext/g4vg/Scaler.hh index 6d73775272..ae05f9f1c8 100644 --- a/src/celeritas/ext/g4vg/Scaler.hh +++ b/src/celeritas/ext/g4vg/Scaler.hh @@ -20,8 +20,8 @@ namespace g4vg /*! * Convert a unit from Geant4 scale to another. * - * Currently the scale is hardcoded as 1/mm but could easily be a class - * attribute. + * Currently the scale is hardcoded as mm (i.e., CLHEP units) but could easily + * be a class attribute. */ class Scaler { diff --git a/src/celeritas/grid/EnergyLossCalculator.hh b/src/celeritas/grid/EnergyLossCalculator.hh index be4797900e..699775a54a 100644 --- a/src/celeritas/grid/EnergyLossCalculator.hh +++ b/src/celeritas/grid/EnergyLossCalculator.hh @@ -12,7 +12,11 @@ namespace celeritas { //---------------------------------------------------------------------------// -//! For now, energy loss calculation has the same behavior as cross sections +/*! + * For now, energy loss calculation has the same behavior as cross sections. + * + * The return value is [MeV / len] but isn't wrapped with a Quantity. + */ using EnergyLossCalculator = XsCalculator; //---------------------------------------------------------------------------// diff --git a/src/celeritas/io/EventWriter.cc b/src/celeritas/io/EventWriter.cc index 4e504a664a..16dd31e7ca 100644 --- a/src/celeritas/io/EventWriter.cc +++ b/src/celeritas/io/EventWriter.cc @@ -133,10 +133,10 @@ void EventWriter::operator()(VecPrimary const& primaries) vtx_primary = &p; HepMC3::FourVector pos; - pos.set_x(p.position[0]); - pos.set_y(p.position[1]); - pos.set_z(p.position[2]); - pos.set_t(p.time / units::centimeter * constants::c_light); + pos.set_x(p.position[0] / units::centimeter); + pos.set_y(p.position[1] / units::centimeter); + pos.set_z(p.position[2] / units::centimeter); + pos.set_t(p.time * (constants::c_light / units::centimeter)); // Need to create a new virtual particle for each vertex auto temp_par = std::make_shared(); diff --git a/src/celeritas/io/ImportData.cc b/src/celeritas/io/ImportData.cc new file mode 100644 index 0000000000..798c2565e8 --- /dev/null +++ b/src/celeritas/io/ImportData.cc @@ -0,0 +1,53 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/io/ImportData.cc +//---------------------------------------------------------------------------// +#include "ImportData.hh" + +#include "corecel/Assert.hh" +#include "corecel/io/Logger.hh" +#include "celeritas/UnitTypes.hh" + +#include "detail/ImportDataConverter.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Recursively convert imported data to the native unit type. + */ +void convert_to_native(ImportData* data) +{ + CELER_EXPECT(data); + + // Convert data + if (data->units.empty()) + { + CELER_LOG(debug) << "Unit system missing from import data: assuming " + "CGS"; + data->units = to_cstring(UnitSystem::cgs); + } + + // Convert string to unit system enum + UnitSystem const usys = to_unit_system(data->units); + CELER_VALIDATE(usys != UnitSystem::none, << "invalid unit system"); + + if (usys == UnitSystem::native) + { + // No unit conversion needed + return; + } + CELER_LOG(info) << "Converting imported units from '" << to_cstring(usys) + << "' to '" << to_cstring(UnitSystem::native) << "'"; + + detail::ImportDataConverter convert{usys}; + convert(data); + + CELER_ENSURE(data->units == units::NativeTraits::label()); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/io/ImportData.hh b/src/celeritas/io/ImportData.hh index 4a466dae97..29f70568d5 100644 --- a/src/celeritas/io/ImportData.hh +++ b/src/celeritas/io/ImportData.hh @@ -40,6 +40,12 @@ namespace celeritas * atomic numbers, and thus are stored in maps. To retrieve specific data use * \c find(atomic_number) . * + * The unit system of the data is stored in the "units" string. If empty + * (backward compatibility) or "CGS" the embedded contents are in CGS. If + * "CLHEP" the units are CLHEP. The \c convert_to_native function will + * convert a data structure in place and update the units label. Refer to \c + * base/Units.hh for further information on unit systems. + * * The "processes" field may be empty for testing applications. */ struct ImportData @@ -64,7 +70,16 @@ struct ImportData ImportSBMap sb_data; ImportLivermorePEMap livermore_pe_data; ImportAtomicRelaxationMap atomic_relaxation_data; + + std::string units; //!< "cgs", "clhep", or "si" }; +//---------------------------------------------------------------------------// +// FREE FUNCTIONS +//---------------------------------------------------------------------------// + +// Recursively convert imported data to the native unit type +void convert_to_native(ImportData* data); + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/io/ImportModel.hh b/src/celeritas/io/ImportModel.hh index 9796c07ed9..57070ae52a 100644 --- a/src/celeritas/io/ImportModel.hh +++ b/src/celeritas/io/ImportModel.hh @@ -11,6 +11,7 @@ #include #include "ImportPhysicsTable.hh" +#include "ImportUnits.hh" namespace celeritas { @@ -131,6 +132,7 @@ char const* to_cstring(ImportModelClass value); // Get the default Geant4 process name char const* to_geant_name(ImportModelClass value); + // Convert a Geant4 process name to an IMC (throw RuntimeError if unsupported) ImportModelClass geant_name_to_import_model_class(std::string_view s); diff --git a/src/celeritas/io/ImportParameters.hh b/src/celeritas/io/ImportParameters.hh index b23fecf7da..6a400c108c 100644 --- a/src/celeritas/io/ImportParameters.hh +++ b/src/celeritas/io/ImportParameters.hh @@ -11,6 +11,8 @@ #include "celeritas/Units.hh" +#include "ImportUnits.hh" + namespace celeritas { //---------------------------------------------------------------------------// @@ -25,6 +27,10 @@ namespace celeritas */ struct ImportEmParameters { +#ifndef SWIG + static constexpr auto energy_units{ImportUnits::mev}; +#endif + //! Energy loss fluctuation bool energy_loss_fluct{false}; //! LPM effect for bremsstrahlung and pair production @@ -64,6 +70,10 @@ struct ImportEmParameters */ struct ImportLoopingThreshold { +#ifndef SWIG + static constexpr auto energy_units{ImportUnits::mev}; +#endif + //! Number of steps a higher-energy looping track takes before it's killed int threshold_trials{10}; //! Energy below which looping tracks are immediately killed [MeV] diff --git a/src/celeritas/io/ImportPhysicsTable.cc b/src/celeritas/io/ImportPhysicsTable.cc index 980bdbba40..475f9b8acf 100644 --- a/src/celeritas/io/ImportPhysicsTable.cc +++ b/src/celeritas/io/ImportPhysicsTable.cc @@ -27,24 +27,5 @@ char const* to_cstring(ImportTableType value) return to_cstring_impl(value); } -//---------------------------------------------------------------------------// -/*! - * Get a printable label for units. - */ -char const* to_cstring(ImportUnits value) -{ - static EnumStringMapper const to_cstring_impl{ - "unitless", - "MeV", - "MeV/len", - "len", - "1/len", - "1/len-MeV", - "MeV^2/len", - "len^2", - }; - return to_cstring_impl(value); -} - //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/io/ImportPhysicsTable.hh b/src/celeritas/io/ImportPhysicsTable.hh index a72bb4e2e5..13931de8d4 100644 --- a/src/celeritas/io/ImportPhysicsTable.hh +++ b/src/celeritas/io/ImportPhysicsTable.hh @@ -9,7 +9,10 @@ #include +#include "celeritas/Units.hh" + #include "ImportPhysicsVector.hh" +#include "ImportUnits.hh" namespace celeritas { @@ -38,31 +41,6 @@ enum class ImportTableType size_ }; -//---------------------------------------------------------------------------// -/*! - * Units of a physics table. - * - * These depend on the unit selection in the main import data - */ -enum class ImportUnits -{ - none, //!< Unitless - mev, //!< Energy [MeV] - mev_per_len, //!< Energy loss [MeV/len] - mev_per_cm = mev_per_len, //!< Deprecated - len, //!< Range [len] - cm = len, //!< Deprecated - len_inv, //!< Macroscopic xs [1/len] - cm_inv = len_inv, //!< Deprecated - len_mev_inv, //!< Scaled [1/E] macroscopic xs [1/len-MeV] - cm_mev_inv = len_mev_inv, //!< Deprecated - mev_sq_per_len, //!< Scaled [E^2] macroscopic xs [MeV^2/len] - mev_2_per_cm = mev_sq_per_len, //!< Deprecated - len_sq, //!< Microscopic cross section [len^2] - cm_2 = len_sq, //!< Deprecated - size_ -}; - //---------------------------------------------------------------------------// /*! * Imported physics table. Each table stores physics vectors for all materials. @@ -70,8 +48,8 @@ enum class ImportUnits struct ImportPhysicsTable { ImportTableType table_type{ImportTableType::size_}; - ImportUnits x_units{ImportUnits::none}; - ImportUnits y_units{ImportUnits::none}; + ImportUnits x_units{ImportUnits::unitless}; + ImportUnits y_units{ImportUnits::unitless}; std::vector physics_vectors; explicit operator bool() const @@ -84,8 +62,8 @@ struct ImportPhysicsTable // FREE FUNCTIONS //---------------------------------------------------------------------------// +// Get the string value for a table type char const* to_cstring(ImportTableType value); -char const* to_cstring(ImportUnits value); //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/io/ImportPhysicsVector.cc b/src/celeritas/io/ImportPhysicsVector.cc index 64c835eea3..b03dda4a36 100644 --- a/src/celeritas/io/ImportPhysicsVector.cc +++ b/src/celeritas/io/ImportPhysicsVector.cc @@ -8,6 +8,7 @@ #include "ImportPhysicsVector.hh" #include "corecel/Assert.hh" +#include "corecel/io/EnumStringMapper.hh" namespace celeritas { @@ -17,10 +18,9 @@ namespace celeritas */ char const* to_cstring(ImportPhysicsVectorType value) { - static char const* const strings[] = {"unknown", "linear", "log", "free"}; - CELER_EXPECT(static_cast(value) * sizeof(char const*) - < sizeof(strings)); - return strings[static_cast(value)]; + static EnumStringMapper const to_cstring_impl{ + "unknown", "linear", "log", "free"}; + return to_cstring_impl(value); } //---------------------------------------------------------------------------// diff --git a/src/celeritas/io/ImportPhysicsVector.hh b/src/celeritas/io/ImportPhysicsVector.hh index 655f1f2e9f..8dcaa505ad 100644 --- a/src/celeritas/io/ImportPhysicsVector.hh +++ b/src/celeritas/io/ImportPhysicsVector.hh @@ -23,6 +23,7 @@ enum class ImportPhysicsVectorType linear, //!< Uniform and linear in x log, //!< Uniform and logarithmic in x free, //!< Nonuniform in x + size_ }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/io/ImportUnits.cc b/src/celeritas/io/ImportUnits.cc new file mode 100644 index 0000000000..7ca4d81a9e --- /dev/null +++ b/src/celeritas/io/ImportUnits.cc @@ -0,0 +1,125 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/io/ImportUnits.cc +//---------------------------------------------------------------------------// +#include "ImportUnits.hh" + +#include + +#include "corecel/io/EnumStringMapper.hh" +#include "corecel/math/Algorithms.hh" +#include "corecel/math/Quantity.hh" +#include "celeritas/UnitTypes.hh" + +using celeritas::units::visit_unit_system; + +namespace celeritas +{ +namespace +{ +//---------------------------------------------------------------------------// +/*! + * Get the length scale of a unit system. + */ +double length_scale(UnitSystem sys) +{ + return visit_unit_system( + [](auto traits) { + using Unit = typename decltype(traits)::Length; + return native_value_from(Quantity{1}); + }, + sys); +} + +//---------------------------------------------------------------------------// +/*! + * Get the time scale of a unit system. + */ +double time_scale(UnitSystem sys) +{ + return visit_unit_system( + [](auto traits) { + using Unit = typename decltype(traits)::Time; + return native_value_from(Quantity{1}); + }, + sys); +} + +//---------------------------------------------------------------------------// +} // namespace + +//---------------------------------------------------------------------------// +/*! + * Get a printable label for units. + */ +char const* to_cstring(ImportUnits value) +{ + static EnumStringMapper const to_cstring_impl{ + "unitless", + "MeV", + "MeV/len", + "len", + "1/len", + "1/len-MeV", + "MeV^2/len", + "len^2", + "time", + "1/len^3", + }; + return to_cstring_impl(value); +} + +//---------------------------------------------------------------------------// +/*! + * Get the native value from a quantity of this type. + */ +double native_value_from(UnitSystem sys, ImportUnits q) +{ + constexpr double mev = 1; + double const len = length_scale(sys); + + switch (q) + { + case ImportUnits::none: + return 1; + case ImportUnits::mev: + return mev; + case ImportUnits::mev_per_len: + return mev / len; + case ImportUnits::len: + return len; + case ImportUnits::len_inv: + return 1 / len; + case ImportUnits::len_mev_inv: + return 1 / (len * mev); + case ImportUnits::mev_sq_per_len: + return ipow<2>(mev) / len; + case ImportUnits::len_sq: + return ipow<2>(len); + case ImportUnits::time: + return time_scale(sys); + case ImportUnits::inv_len_cb: + return 1 / ipow<3>(len); + case ImportUnits::size_: + CELER_ASSERT_UNREACHABLE(); + } + CELER_ASSERT_UNREACHABLE(); +} + +//---------------------------------------------------------------------------// +/*! + * Get the native value from a unit CLHEP quantity of this type. + * + * Multiply a Geant4 quantity by the result to convert to the native unit + * system. + */ +double native_value_from_clhep(ImportUnits q) +{ + return native_value_from(celeritas::UnitSystem::clhep, q); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/io/ImportUnits.hh b/src/celeritas/io/ImportUnits.hh new file mode 100644 index 0000000000..d62736bd50 --- /dev/null +++ b/src/celeritas/io/ImportUnits.hh @@ -0,0 +1,60 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/io/ImportUnits.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/Types.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Quantity of measure enumeration for imported data. + * + * These enumeration values are used to scale values between the Celeritas + * native unit system and the CLHEP/Geant4 values. + * + * \todo Rename to ImportUnit?? + */ +enum class ImportUnits +{ + unitless, //!< No dimension + mev, //!< Energy [MeV] + mev_per_len, //!< Energy loss [MeV/len] + len, //!< Range [len] + len_inv, //!< Macroscopic xs [1/len] + len_mev_inv, //!< Scaled (1/E) macroscopic xs [1/len-MeV] + mev_sq_per_len, //!< Scaled [E^2] macroscopic xs [MeV^2/len] + len_sq, //!< Microscopic cross section [len^2] + time, //!< Time [time] + inv_len_cb, //!< Number density [1/len^3] + size_, + // Deprecated aliases + none = unitless, //!< Deprecated + mev_per_cm = mev_per_len, //!< Deprecated + cm = len, //!< Deprecated + cm_inv = len_inv, //!< Deprecated + cm_mev_inv = len_mev_inv, //!< Deprecated + mev_2_per_cm = mev_sq_per_len, //!< Deprecated + cm_2 = len_sq, //!< Deprecated +}; + +//---------------------------------------------------------------------------// +// FREE FUNCTIONS +//---------------------------------------------------------------------------// + +// Get the string label for units +char const* to_cstring(ImportUnits q); + +// Get the multiplier to turn this quantity to a native value +double native_value_from(UnitSystem sys, ImportUnits q); + +// Get the multiplier to turn a unit Geant4 value to a native value +double native_value_from_clhep(ImportUnits q); + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/io/detail/ImportDataConverter.cc b/src/celeritas/io/detail/ImportDataConverter.cc new file mode 100644 index 0000000000..8c1b9c7980 --- /dev/null +++ b/src/celeritas/io/detail/ImportDataConverter.cc @@ -0,0 +1,170 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/io/detail/ImportDataConverter.cc +//---------------------------------------------------------------------------// +#include "ImportDataConverter.hh" + +#include "corecel/Assert.hh" +#include "celeritas/UnitTypes.hh" + +#include "../ImportData.hh" +#include "../ImportUnits.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Construct with a unit system. + */ +ImportDataConverter::ImportDataConverter(UnitSystem usys) : usys_{usys} +{ + len_ = native_value_from(usys_, ImportUnits::len); + numdens_ = native_value_from(usys_, ImportUnits::inv_len_cb); + time_ = native_value_from(usys_, ImportUnits::time); + xs_ = native_value_from(usys_, ImportModelMaterial::xs_units); +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportData* data) +{ + for (auto& m : data->materials) + { + (*this)(&m); + } + + for (auto& p : data->processes) + { + (*this)(&p); + } + + for (auto& m : data->msc_models) + { + (*this)(&m); + } + + (*this)(&data->em_params); + + data->units = units::NativeTraits::label(); +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportElement* data) +{ + CELER_EXPECT(data); +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportEmParameters* data) +{ + CELER_EXPECT(data); + data->msc_lambda_limit *= len_; +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportMaterial* data) +{ + CELER_EXPECT(data); + + data->number_density *= numdens_; + + for (auto& [pdg, cut] : data->pdg_cutoffs) + { + cut.range *= len_; + } +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportModelMaterial* data) +{ + CELER_EXPECT(data); + + for (auto& xs_vec : data->micro_xs) + { + for (double& xs : xs_vec) + { + xs *= xs_; + } + } +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportModel* data) +{ + CELER_EXPECT(data); + + for (auto& mm : data->materials) + { + (*this)(&mm); + } +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportMscModel* data) +{ + CELER_EXPECT(data); + + (*this)(&data->xs_table); +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportParticle* data) +{ + CELER_EXPECT(data); + + data->lifetime *= time_; +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportPhysicsTable* data) +{ + CELER_EXPECT(data); + + if (double const units = native_value_from(usys_, data->x_units); + units != 1) + { + for (auto& v : data->physics_vectors) + { + for (auto& xval : v.x) + { + xval *= units; + } + } + } + + if (double const units = native_value_from(usys_, data->y_units); + units != 1) + { + for (auto& v : data->physics_vectors) + { + for (auto& yval : v.y) + { + yval *= units; + } + } + } +} + +//---------------------------------------------------------------------------// +void ImportDataConverter::operator()(ImportProcess* data) +{ + CELER_EXPECT(data); + + for (auto& m : data->models) + { + (*this)(&m); + } + + for (auto& t : data->tables) + { + (*this)(&t); + } +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/io/detail/ImportDataConverter.hh b/src/celeritas/io/detail/ImportDataConverter.hh new file mode 100644 index 0000000000..e798d4265c --- /dev/null +++ b/src/celeritas/io/detail/ImportDataConverter.hh @@ -0,0 +1,61 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/io/detail/ImportDataConverter.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/Types.hh" + +namespace celeritas +{ +struct ImportData; +struct ImportElement; +struct ImportEmParameters; +struct ImportMaterial; +struct ImportModel; +struct ImportModelMaterial; +struct ImportMscModel; +struct ImportParticle; +struct ImportPhysicsTable; +struct ImportProcess; + +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Convert imported data from one unit system to another. + */ +class ImportDataConverter +{ + public: + // Construct with a unit system + explicit ImportDataConverter(UnitSystem usys); + + //!@{ + //! Convert imported data to the native unit type + void operator()(ImportData* data); + void operator()(ImportElement* data); + void operator()(ImportEmParameters* data); + void operator()(ImportMaterial* data); + void operator()(ImportModel* data); + void operator()(ImportModelMaterial* data); + void operator()(ImportMscModel* data); + void operator()(ImportParticle* data); + void operator()(ImportPhysicsTable* data); + void operator()(ImportProcess* data); + //!@} + + private: + UnitSystem usys_; + double len_; + double numdens_; + double time_; + double xs_; +}; + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mat/ElementSelector.hh b/src/celeritas/mat/ElementSelector.hh index c94ddb911c..981777868a 100644 --- a/src/celeritas/mat/ElementSelector.hh +++ b/src/celeritas/mat/ElementSelector.hh @@ -49,8 +49,8 @@ namespace celeritas * * Note that the units of the calculated microscopic cross section will be * identical to the units returned by `calc_micro_xs`. The macroscopic cross - * section units (micro times \c mat.number_density() ) will be 1/cm if and - * only if calc_micro units are cm^2. + * section units (micro times \c mat.number_density() ) will be 1/len if and + * only if calc_micro units are len^2. * * \todo Refactor to use Selector. */ diff --git a/src/celeritas/phys/Interaction.hh b/src/celeritas/phys/Interaction.hh index 83e9e0d71a..1caa23e891 100644 --- a/src/celeritas/phys/Interaction.hh +++ b/src/celeritas/phys/Interaction.hh @@ -89,7 +89,7 @@ struct MscStep bool is_displaced{true}; //!< Flag for the lateral displacement real_type true_path{}; //!< True path length due to the msc [len] real_type geom_path{}; //!< Geometrical path length [len] - real_type alpha = small_step_alpha(); //!< Scaled MFP slope [len^-1] + real_type alpha = small_step_alpha(); //!< Scaled MFP slope [1/len] }; //---------------------------------------------------------------------------// @@ -101,9 +101,9 @@ struct MscStep */ struct MscRange { - real_type range_init{}; //!< Initial msc range + real_type range_init{}; //!< Initial msc range [len] real_type range_fact{}; //!< Scale factor for the msc range - real_type limit_min{}; //!< Minimum of the true path limit + real_type limit_min{}; //!< Minimum of the true path limit [len] explicit CELER_FUNCTION operator bool() const { diff --git a/src/celeritas/phys/ParticleTrackView.hh b/src/celeritas/phys/ParticleTrackView.hh index 269a7061fd..8491ceefde 100644 --- a/src/celeritas/phys/ParticleTrackView.hh +++ b/src/celeritas/phys/ParticleTrackView.hh @@ -80,7 +80,7 @@ class ParticleTrackView // Charge [elemental charge e+] CELER_FORCEINLINE_FUNCTION units::ElementaryCharge charge() const; - // Decay constant [1/s] + // Decay constant [1 / time] CELER_FORCEINLINE_FUNCTION real_type decay_constant() const; // Whether it is an antiparticle @@ -228,7 +228,7 @@ CELER_FUNCTION units::ElementaryCharge ParticleTrackView::charge() const //---------------------------------------------------------------------------// /*! - * Decay constant. + * Decay constant in native units. */ CELER_FUNCTION real_type ParticleTrackView::decay_constant() const { diff --git a/src/celeritas/phys/PhysicsStepUtils.hh b/src/celeritas/phys/PhysicsStepUtils.hh index 038c158a02..b0cf12469a 100644 --- a/src/celeritas/phys/PhysicsStepUtils.hh +++ b/src/celeritas/phys/PhysicsStepUtils.hh @@ -45,7 +45,7 @@ calc_physics_step_limit(MaterialTrackView const& material, using VGT = ValueGridType; // TODO: for particles with decay, macro XS calculation will incorporate - // decay probability, dividing decay constant by speed to become 1/cm to + // decay probability, dividing decay constant by speed to become 1/len to // compete with interactions // Loop over all processes that apply to this track (based on particle diff --git a/src/celeritas/phys/ProcessBuilder.cc b/src/celeritas/phys/ProcessBuilder.cc index 0750093b6f..16cb336893 100644 --- a/src/celeritas/phys/ProcessBuilder.cc +++ b/src/celeritas/phys/ProcessBuilder.cc @@ -49,6 +49,9 @@ auto ProcessBuilder::get_all_process_classes( /*! * Construct imported process data. * + * \pre The import data must have already been converted to the native unit + * system. + * * \warning If Livermore and SB data is present in the import data, their * lifetime must extend beyond the \c ProcessBuilder instance. */ @@ -67,6 +70,7 @@ ProcessBuilder::ProcessBuilder(ImportData const& data, { CELER_EXPECT(input_.material); CELER_EXPECT(input_.particle); + CELER_EXPECT(std::string(data.units) == units::NativeTraits::label()); input_.imported = std::make_shared(data.processes); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3431dd66e9..4da5bdf638 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -104,6 +104,13 @@ else() set(_fixme_single DISABLE) endif() +if(CELERITAS_UNITS STREQUAL "CGS") + set(_fixme_cgs) +else() + # Unit tests haven't yet been updated to include the correct units + set(_fixme_cgs DISABLE) +endif() + #-----------------------------------------------------------------------------# # GOOGLETEST EXTENSION TESTS #-----------------------------------------------------------------------------# @@ -159,8 +166,8 @@ celeritas_add_test(corecel/grid/UniformGrid.test.cc) # IO set(CELERITASTEST_PREFIX corecel/io) celeritas_add_test(corecel/io/EnumStringMapper.test.cc) -celeritas_add_test(corecel/io/Label.test.cc) celeritas_add_test(corecel/io/Join.test.cc) +celeritas_add_test(corecel/io/Label.test.cc) celeritas_add_test(corecel/io/Logger.test.cc) celeritas_add_test(corecel/io/OutputRegistry.test.cc LINK_LIBRARIES ${nlohmann_json_LIBRARIES}) @@ -401,6 +408,7 @@ if(CELERITAS_USE_VecGeom) celeritas_add_device_test(celeritas/ext/Vecgeom FILTER ${_vecgeom_tests} + ${_fixme_cgs} ) endif() if(CELERITAS_USE_Geant4 AND CELERITAS_REAL_TYPE STREQUAL "double") @@ -413,7 +421,8 @@ if(CELERITAS_USE_Geant4 AND CELERITAS_REAL_TYPE STREQUAL "double") endif() celeritas_add_test(celeritas/ext/GeantImporter.test.cc - ${_needs_geant4} LINK_LIBRARIES ${nlohmann_json_LIBRARIES} + ${_needs_geant4} + LINK_LIBRARIES ${nlohmann_json_LIBRARIES} FILTER "FourSteelSlabs*" "TestEm3*" @@ -430,10 +439,10 @@ celeritas_add_test(celeritas/ext/RootImporter.test.cc ${_needs_root}) set(CELERITASTEST_PREFIX celeritas/field) celeritas_add_test(celeritas/field/Fields.test.cc) -celeritas_add_test(celeritas/field/Steppers.test.cc) -celeritas_add_test(celeritas/field/FieldDriver.test.cc) +celeritas_add_test(celeritas/field/Steppers.test.cc ${_fixme_cgs}) +celeritas_add_test(celeritas/field/FieldDriver.test.cc ${_fixme_cgs}) celeritas_add_test(celeritas/field/FieldPropagator.test.cc - ${_needs_geo}) + ${_needs_geo} ${_fixme_cgs}) celeritas_add_test(celeritas/field/LinearPropagator.test.cc ${_needs_geo} LINK_LIBRARIES ${_geo_libs}) celeritas_add_test(celeritas/field/MagFieldEquation.test.cc) @@ -448,7 +457,7 @@ if(CELERITAS_USE_CUDA OR CELERITAS_USE_HIP) list(APPEND _geo_args celeritas/geo/HeuristicGeoTestBase.cu) endif() celeritas_add_test(celeritas/geo/Geometry.test.cc - ${_geo_args}) + ${_geo_args} ${_fixme_cgs}) if(NOT (CELERITAS_USE_Geant4 OR CELERITAS_USE_ROOT)) set(_needs_geant_or_root DISABLE) @@ -536,6 +545,7 @@ celeritas_add_test(celeritas/grid/XsCalculator.test.cc) set(CELERITASTEST_PREFIX celeritas/io) celeritas_add_test(celeritas/io/EventIO.test.cc ${_needs_hepmc} LINK_LIBRARIES ${HepMC3_LIBRARIES}) +celeritas_add_test(celeritas/io/ImportUnits.test.cc) celeritas_add_test(celeritas/io/RootEventIO.test.cc ${_needs_root}) celeritas_add_test(celeritas/io/SeltzerBergerReader.test.cc ${_needs_geant4}) diff --git a/test/Test.cc b/test/Test.cc index d189d21c0b..6462b8dc82 100644 --- a/test/Test.cc +++ b/test/Test.cc @@ -134,12 +134,16 @@ bool Test::strict_testing() // Disable strict testing for Geant4 return false; } - if (CELERITAS_REAL_TYPE != CELERITAS_REAL_TYPE_DOUBLE) { // Disable strict testing for single precision return false; } + if (CELERITAS_UNITS != CELERITAS_UNITS_CGS) + { + // Disable strict testing for non-CLHEP units + return false; + } return !envstr.empty(); } diff --git a/test/accel/HepMC3PrimaryGenerator.test.cc b/test/accel/HepMC3PrimaryGenerator.test.cc index 94821846b9..0fee4b318a 100644 --- a/test/accel/HepMC3PrimaryGenerator.test.cc +++ b/test/accel/HepMC3PrimaryGenerator.test.cc @@ -7,6 +7,8 @@ //---------------------------------------------------------------------------// #include "accel/HepMC3PrimaryGenerator.hh" +#include + #include "celeritas/GeantTestBase.hh" #include "celeritas/SimpleCmsTestBase.hh" #include "celeritas/ext/Convert.geant.hh" diff --git a/test/accel/detail/HitProcessor.test.cc b/test/accel/detail/HitProcessor.test.cc index a00912bcce..5901124858 100644 --- a/test/accel/detail/HitProcessor.test.cc +++ b/test/accel/detail/HitProcessor.test.cc @@ -10,6 +10,7 @@ #include #include "celeritas/SimpleCmsTestBase.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/phys/PDGNumber.hh" #include "celeritas/user/DetectorSteps.hh" #include "celeritas/user/StepData.hh" @@ -18,7 +19,9 @@ #include "celeritas_test.hh" -using namespace celeritas::units; +using celeritas::test::from_cm; +using celeritas::test::SimpleHitsResult; +using celeritas::units::MevEnergy; namespace celeritas { @@ -26,7 +29,6 @@ namespace detail { namespace test { -using ::celeritas::test::SimpleHitsResult; //---------------------------------------------------------------------------// class SimpleCmsTest : public ::celeritas::test::SDTestBase, @@ -145,6 +147,8 @@ DetectorStepOutput SimpleCmsTest::make_dso() const } if (selection_.points[StepPoint::post].time) { + using celeritas::units::second; + dso.points[StepPoint::post].time = { 1e-9 * second, 2e-10 * second, @@ -155,9 +159,9 @@ DetectorStepOutput SimpleCmsTest::make_dso() const { // note: points must correspond to detector volumes! dso.points[StepPoint::pre].pos = { - {100, 0, 0}, - {0, 150, 10}, - {0, 200, -20}, + from_cm(Real3{100, 0, 0}), + from_cm(Real3{0, 150, 10}), + from_cm(Real3{0, 200, -20}), }; } if (selection_.points[StepPoint::pre].dir) @@ -275,19 +279,19 @@ TEST_F(SimpleCmsTest, touchable_edgecase) auto& pos = dso_hits.points[StepPoint::pre].pos; auto& dir = dso_hits.points[StepPoint::pre].dir; pos = { - {30, 0, 0}, - {0, 125, 10}, - {0, 175, -20}, + from_cm(Real3{30, 0, 0}), + from_cm(Real3{0, 125, 10}), + from_cm(Real3{0, 175, -20}), }; process_hits(dso_hits); pos = { - {-120.20472398905, 34.290294993135, -58.348475076307}, - {-58.042349740868, -165.09417202481, -315.41125902053}, - {0, 275, -20}, + from_cm(Real3{-120.20472398905, 34.290294993135, -58.348475076307}), + from_cm(Real3{-58.042349740868, -165.09417202481, -315.41125902053}), + from_cm(Real3{0, 275, -20}), }; - EXPECT_SOFT_EQ(125.0, std::hypot(pos[0][0], pos[0][1])); - EXPECT_SOFT_EQ(175.0, std::hypot(pos[1][0], pos[1][1])); + EXPECT_SOFT_EQ(from_cm(125.0), std::hypot(pos[0][0], pos[0][1])); + EXPECT_SOFT_EQ(from_cm(175.0), std::hypot(pos[1][0], pos[1][1])); dir = { {0.39117837162751, -0.78376148752334, -0.48238720157779}, {0.031769215780742, 0.6378450322959, -0.76950921482729}, diff --git a/test/accel/detail/TouchableUpdater.test.cc b/test/accel/detail/TouchableUpdater.test.cc index 1ddf350964..7bf49eb22b 100644 --- a/test/accel/detail/TouchableUpdater.test.cc +++ b/test/accel/detail/TouchableUpdater.test.cc @@ -16,11 +16,13 @@ #include "corecel/math/ArrayOperators.hh" #include "corecel/math/ArrayUtils.hh" #include "celeritas/GenericGeoTestBase.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/Units.hh" #include "celeritas/ext/GeantGeoParams.hh" #include "celeritas_test.hh" +using celeritas::test::from_cm; using celeritas::test::ScopedLogStorer; namespace celeritas @@ -82,11 +84,13 @@ class TouchableUpdaterTest : public ::celeritas::test::GenericGeantGeoTestBase TEST_F(TouchableUpdaterTest, correct) { TouchableUpdater update = this->make_touchable_updater(); - Real3 const dir{1, 0, 0}; + auto update_cm = [&](Real3 const& pos_cm, std ::string lv_name) { + return update(from_cm(pos_cm), Real3{1, 0, 0}, this->find_lv(lv_name)); + }; - EXPECT_TRUE(update({15, 0, 0}, dir, this->find_lv("vacuum_tube"))); - EXPECT_TRUE(update({100, 0, 0}, dir, this->find_lv("si_tracker"))); - EXPECT_TRUE(update({150, 0, 0}, dir, this->find_lv("em_calorimeter"))); + EXPECT_TRUE(update_cm({15, 0, 0}, "vacuum_tube")); + EXPECT_TRUE(update_cm({100, 0, 0}, "si_tracker")); + EXPECT_TRUE(update_cm({150, 0, 0}, "em_calorimeter")); } TEST_F(TouchableUpdaterTest, just_inside) @@ -99,11 +103,11 @@ TEST_F(TouchableUpdaterTest, just_inside) ScopedLogStorer scoped_log_{&celeritas::self_logger(), LogLevel::diagnostic}; - EXPECT_TRUE(update({30 + eps, 0, 0}, {1, 0, 0}, tracker_lv)); - EXPECT_TRUE(update({125 - eps, 0, 0}, {1, 0, 0}, tracker_lv)); + EXPECT_TRUE(update({from_cm(30) + eps, 0, 0}, {1, 0, 0}, tracker_lv)); + EXPECT_TRUE(update({from_cm(125) - eps, 0, 0}, {1, 0, 0}, tracker_lv)); - EXPECT_TRUE(update({125 + eps, 0, 0}, {-1, 0, 0}, calo_lv)); - EXPECT_TRUE(update({175 - eps, 0, 0}, {-1, 0, 0}, calo_lv)); + EXPECT_TRUE(update({from_cm(125) + eps, 0, 0}, {-1, 0, 0}, calo_lv)); + EXPECT_TRUE(update({from_cm(175) - eps, 0, 0}, {-1, 0, 0}, calo_lv)); EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; } @@ -115,10 +119,11 @@ TEST_F(TouchableUpdaterTest, coincident) LogLevel::diagnostic}; // Coincident point should work in either volume, in or out + real_type const r = from_cm(125); for (char const* lv : {"si_tracker", "em_calorimeter"}) { - EXPECT_TRUE(update({125, 0, 0}, {1, 0, 0}, this->find_lv(lv))); - EXPECT_TRUE(update({125, 0, 0}, {-1, 0, 0}, this->find_lv(lv))); + EXPECT_TRUE(update({r, 0, 0}, {1, 0, 0}, this->find_lv(lv))); + EXPECT_TRUE(update({r, 0, 0}, {-1, 0, 0}, this->find_lv(lv))); } EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; @@ -133,9 +138,9 @@ TEST_F(TouchableUpdaterTest, coincident_tangent) // TODO: we can't seem to test the volume on the other side of an exact // coincident surface on a tangent - EXPECT_FALSE(update({125, 0, 0}, {0, 1, 0}, this->find_lv("si_tracker"))); - EXPECT_TRUE( - update({125, 0, 0}, {0, 1, 0}, this->find_lv("em_calorimeter"))); + real_type const r = from_cm(125); + EXPECT_FALSE(update({r, 0, 0}, {0, 1, 0}, this->find_lv("si_tracker"))); + EXPECT_TRUE(update({r, 0, 0}, {0, 1, 0}, this->find_lv("em_calorimeter"))); static char const* const expected_log_messages[] = { "Failed to bump navigation state up to a distance of 1 [mm] at {1250, " @@ -157,8 +162,10 @@ TEST_F(TouchableUpdaterTest, just_outside_nowarn) for (real_type xdir : {1.0, -1.0}) { - EXPECT_TRUE(update({30 - eps, 0, 0}, {xdir, 0, 0}, tracker_lv)); - EXPECT_TRUE(update({125 + 2 * eps, 0, 0}, {-xdir, 0, 0}, tracker_lv)); + EXPECT_TRUE( + update({from_cm(30) - eps, 0, 0}, {xdir, 0, 0}, tracker_lv)); + EXPECT_TRUE( + update({from_cm(125) + 2 * eps, 0, 0}, {-xdir, 0, 0}, tracker_lv)); } EXPECT_TRUE(scoped_log_.empty()) << scoped_log_; @@ -175,8 +182,10 @@ TEST_F(TouchableUpdaterTest, just_outside_warn) for (real_type xdir : {1.0, -1.0}) { - EXPECT_TRUE(update({30 - eps, 0, 0}, {xdir, 0, 0}, tracker_lv)); - EXPECT_TRUE(update({125 + eps, 0, 0}, {-xdir, 0, 0}, tracker_lv)); + EXPECT_TRUE( + update({from_cm(30) - eps, 0, 0}, {xdir, 0, 0}, tracker_lv)); + EXPECT_TRUE( + update({from_cm(125) + eps, 0, 0}, {-xdir, 0, 0}, tracker_lv)); } static char const* const expected_log_messages[] @@ -196,7 +205,10 @@ TEST_F(TouchableUpdaterTest, just_outside_warn) "0} [mm] along {-1, -0, -0} from {{pv='em_calorimeter_pv', " "lv=2='em_calorimeter'}} to try to reach \"si_tracker\"@0x0 (ID=1)", "...bumped to {{pv='si_tracker_pv', lv=1='si_tracker'}}"}; - EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()); + if (CELERITAS_UNITS == CELERITAS_UNITS_CGS) + { + EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()); + } static char const* const expected_log_levels[] = {"warning", "diagnostic", "warning", @@ -219,8 +231,10 @@ TEST_F(TouchableUpdaterTest, too_far) for (real_type xdir : {1.0, -1.0}) { - EXPECT_FALSE(update({30 - eps, 0, 0}, {xdir, 0, 0}, tracker_lv)); - EXPECT_FALSE(update({125 + eps, 0, 0}, {-xdir, 0, 0}, tracker_lv)); + EXPECT_FALSE( + update({from_cm(30) - eps, 0, 0}, {xdir, 0, 0}, tracker_lv)); + EXPECT_FALSE( + update({from_cm(125) + eps, 0, 0}, {-xdir, 0, 0}, tracker_lv)); } static char const* const expected_log_messages[] = { @@ -323,7 +337,10 @@ TEST_F(TouchableUpdaterTest, regression) "{{pv='vacuum_tube_pv', lv=0='vacuum_tube'}} to try to reach " "\"si_tracker\"@0x0 (ID=1)", "...bumped to {{pv='si_tracker_pv', lv=1='si_tracker'}}"}; - EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()); + if (CELERITAS_UNITS == CELERITAS_UNITS_CGS) + { + EXPECT_VEC_EQ(expected_log_messages, scoped_log_.messages()); + } static char const* const expected_log_levels[] = {"warning", "diagnostic", "warning", diff --git a/test/celeritas/AllGeoTypedTestBase.hh b/test/celeritas/AllGeoTypedTestBase.hh index 2c5aa6c69d..c75761ef66 100644 --- a/test/celeritas/AllGeoTypedTestBase.hh +++ b/test/celeritas/AllGeoTypedTestBase.hh @@ -68,6 +68,11 @@ class AllGeoTypedTestBase : public GenericGeoTestBase public: using SPConstGeo = typename GenericGeoTestBase::SPConstGeo; + static std::string geo_name() + { + return testdetail::GenericGeoTraits::name; + } + SPConstGeo build_geometry() override { return this->build_geometry_from_basename(); diff --git a/test/celeritas/GeantTestBase.cc b/test/celeritas/GeantTestBase.cc index f3ab08b5d4..8093038a43 100644 --- a/test/celeritas/GeantTestBase.cc +++ b/test/celeritas/GeantTestBase.cc @@ -72,6 +72,7 @@ bool GeantTestBase::is_ci_build() { return CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE && CELERITAS_CORE_GEO != CELERITAS_CORE_GEO_GEANT4 + && CELERITAS_UNITS == CELERITAS_UNITS_CGS && cstring_equal(celeritas_core_rng, "xorwow") && (cstring_equal(celeritas_clhep_version, "2.4.6.0") || cstring_equal(celeritas_clhep_version, "2.4.6.4")) diff --git a/test/celeritas/GenericGeoTestBase.cc b/test/celeritas/GenericGeoTestBase.cc index 411bcd1041..16757560f6 100644 --- a/test/celeritas/GenericGeoTestBase.cc +++ b/test/celeritas/GenericGeoTestBase.cc @@ -23,6 +23,7 @@ #include "celeritas/ext/GeantImporter.hh" #include "celeritas/io/ImportVolume.hh" +#include "UnitUtils.hh" #include "geo/CheckedGeoTrackView.hh" #if CELERITAS_USE_VECGEOM @@ -251,51 +252,40 @@ auto GenericGeoTestBase::make_geo_track_view() -> GeoTrackView //---------------------------------------------------------------------------// // Get and initialize a single-thread host track view template -auto GenericGeoTestBase::make_geo_track_view(Real3 const& pos, Real3 dir) +auto GenericGeoTestBase::make_geo_track_view(Real3 const& pos_cm, Real3 dir) -> GeoTrackView { auto geo = this->make_geo_track_view(); - geo = GeoTrackInitializer{pos, make_unit_vector(dir)}; + geo = GeoTrackInitializer{from_cm(pos_cm), make_unit_vector(dir)}; return geo; } -//---------------------------------------------------------------------------// -// Get and initialize a single-thread host track view -template -auto GenericGeoTestBase::calc_bump_pos(GeoTrackView const& geo, - real_type delta) const -> Real3 -{ - CELER_EXPECT(delta > 0); - auto pos = geo.pos(); - axpy(delta, geo.dir(), &pos); - return pos; -} - //---------------------------------------------------------------------------// template -auto GenericGeoTestBase::track(Real3 const& pos, Real3 const& dir) +auto GenericGeoTestBase::track(Real3 const& pos_cm, Real3 const& dir) -> TrackingResult { - return this->track(pos, dir, std::numeric_limits::max()); + return this->track(pos_cm, dir, std::numeric_limits::max()); } //---------------------------------------------------------------------------// template -auto GenericGeoTestBase::track(Real3 const& pos, +auto GenericGeoTestBase::track(Real3 const& pos_cm, Real3 const& dir, int max_step) -> TrackingResult { CELER_EXPECT(max_step > 0); TrackingResult result; - GeoTrackView geo = CheckedGeoTrackView{this->make_geo_track_view(pos, dir)}; + GeoTrackView geo + = CheckedGeoTrackView{this->make_geo_track_view(pos_cm, dir)}; if (geo.is_outside()) { - // Initial step is outside but may approach insidfe + // Initial step is outside but may approach inside result.volumes.push_back("[OUTSIDE]"); auto next = geo.find_next_step(); - result.distances.push_back(next.distance); + result.distances.push_back(to_cm(next.distance)); if (next.boundary) { geo.move_to_boundary(); @@ -309,7 +299,7 @@ auto GenericGeoTestBase::track(Real3 const& pos, { result.volumes.push_back(this->volume_name(geo)); auto next = geo.find_next_step(); - result.distances.push_back(next.distance); + result.distances.push_back(to_cm(next.distance)); if (!next.boundary) { // Failure to find the next boundary while inside the geometry @@ -317,11 +307,11 @@ auto GenericGeoTestBase::track(Real3 const& pos, result.volumes.push_back("[NO INTERCEPT]"); break; } - if (next.distance > real_type(1e-7)) + if (next.distance > real_type(from_cm(1e-7))) { geo.move_internal(next.distance / 2); geo.find_next_step(); - result.halfway_safeties.push_back(geo.find_safety()); + result.halfway_safeties.push_back(to_cm(geo.find_safety())); if (result.halfway_safeties.back() > 0) { @@ -343,7 +333,7 @@ auto GenericGeoTestBase::track(Real3 const& pos, << init.pos << " along " << init.dir << " from " << result.volumes.back() << " to " << this->volume_name(geo) << " (alleged safety: " - << result.halfway_safeties.back() << ")"; + << to_cm(result.halfway_safeties.back()) << ")"; result.volumes.back() += "/" + this->volume_name(geo); } auto new_next = geo.find_next_step(); diff --git a/test/celeritas/GenericGeoTestBase.hh b/test/celeritas/GenericGeoTestBase.hh index 2ba4dbc1c8..9cb09204c6 100644 --- a/test/celeritas/GenericGeoTestBase.hh +++ b/test/celeritas/GenericGeoTestBase.hh @@ -38,8 +38,8 @@ namespace test struct GenericGeoTrackingResult { std::vector volumes; - std::vector distances; - std::vector halfway_safeties; + std::vector distances; //!< [cm] + std::vector halfway_safeties; //!< [cm] void print_expected(); }; @@ -152,15 +152,12 @@ class GenericGeoTestBase : virtual public Test, private LazyGeoManager //! Get a single-thread host track view GeoTrackView make_geo_track_view(); //! Get and initialize a single-thread host track view - GeoTrackView make_geo_track_view(Real3 const& pos, Real3 dir); - - //! Calculate a "bumped" position based on the geo's state - Real3 calc_bump_pos(GeoTrackView const& geo, real_type delta) const; + GeoTrackView make_geo_track_view(Real3 const& pos_cm, Real3 dir); //! Find linear segments until outside - TrackingResult track(Real3 const& pos, Real3 const& dir); + TrackingResult track(Real3 const& pos_cm, Real3 const& dir); //! Find linear segments until outside (maximum count - TrackingResult track(Real3 const& pos, Real3 const& dir, int max_step); + TrackingResult track(Real3 const& pos_cm, Real3 const& dir, int max_step); //! Try to map Geant4 volumes using ImportVolume and name GeantVolResult diff --git a/test/celeritas/MockTestBase.cc b/test/celeritas/MockTestBase.cc index d7d5695d4c..509f3e35a2 100644 --- a/test/celeritas/MockTestBase.cc +++ b/test/celeritas/MockTestBase.cc @@ -7,6 +7,7 @@ //---------------------------------------------------------------------------// #include "MockTestBase.hh" +#include "corecel/math/Algorithms.hh" #include "celeritas/geo/GeoMaterialParams.hh" #include "celeritas/global/ActionRegistry.hh" #include "celeritas/global/alongstep/AlongStepGeneralLinearAction.hh" @@ -62,23 +63,23 @@ auto MockTestBase::build_material() -> SPConstMaterial inp.elements = {{AtomicNumber{1}, AmuMass{1.0}, {}, "celerogen"}, {AtomicNumber{4}, AmuMass{10.0}, {}, "celerinium"}, {AtomicNumber{15}, AmuMass{30.0}, {}, "celeron"}}; - inp.materials.push_back({1e20, + inp.materials.push_back({native_value_from(InvCcDensity{1e20}), 300, MatterState::gas, {{ElementId{0}, 1.0}}, "lo density celerogen"}); - inp.materials.push_back({1e21, + inp.materials.push_back({native_value_from(InvCcDensity{1e21}), 300, MatterState::liquid, {{ElementId{0}, 1.0}}, "hi density celerogen"}); inp.materials.push_back( - {1e23, + {native_value_from(InvCcDensity{1e23}), 300, MatterState::solid, {{ElementId{0}, 0.1}, {ElementId{1}, 0.3}, {ElementId{2}, 0.6}}, "celer composite"}); - inp.materials.push_back({1.0, + inp.materials.push_back({native_value_from(InvCcDensity{1.0}), 2.7, MatterState::gas, {{ElementId{0}, 1.0}}, @@ -181,7 +182,8 @@ auto MockTestBase::build_physics() -> SPConstPhysics make_applicability("celeriton", 1, 10), make_applicability("celeriton", 10, 100)}; inp.xs = {Barn{3.0}, Barn{3.0}}; - inp.energy_loss = 0.6 * 1e-20; // 0.6 MeV/cm in celerogen + inp.energy_loss = MevCmSqLossDens{0.6 * 1e-20}; // 0.6 MeV/cm in + // celerogen physics_inp.processes.push_back(std::make_shared(inp)); } { @@ -191,7 +193,7 @@ auto MockTestBase::build_physics() -> SPConstPhysics inp.applic = {make_applicability("anti-celeriton", 1e-3, 1), make_applicability("anti-celeriton", 1, 100)}; inp.xs = {Barn{4.0}, Barn{4.0}}; - inp.energy_loss = 0.7 * 1e-20; + inp.energy_loss = MevCmSqLossDens{0.7 * 1e-20}; physics_inp.processes.push_back(std::make_shared(inp)); } { @@ -209,7 +211,7 @@ auto MockTestBase::build_physics() -> SPConstPhysics inp.use_integral_xs = true; inp.applic = {make_applicability("electron", 1e-5, 10)}; inp.xs = {Barn{0}, Barn{6.0}, Barn{12.0}, Barn{6.0}}; - inp.energy_loss = 0.5 * 1e-20; + inp.energy_loss = MevCmSqLossDens{0.5 * 1e-20}; physics_inp.processes.push_back(std::make_shared(inp)); } return std::make_shared(std::move(physics_inp)); diff --git a/test/celeritas/SimpleTestBase.cc b/test/celeritas/SimpleTestBase.cc index d90609427e..335e6c81ea 100644 --- a/test/celeritas/SimpleTestBase.cc +++ b/test/celeritas/SimpleTestBase.cc @@ -13,6 +13,7 @@ #include "celeritas/global/ActionRegistry.hh" #include "celeritas/global/alongstep/AlongStepNeutralAction.hh" #include "celeritas/io/ImportProcess.hh" +#include "celeritas/io/detail/ImportDataConverter.hh" #include "celeritas/mat/MaterialParams.hh" #include "celeritas/phys/CutoffParams.hh" #include "celeritas/phys/ImportedProcessAdapter.hh" @@ -33,7 +34,7 @@ auto SimpleTestBase::build_material() -> SPConstMaterial MaterialParams::Input inp; inp.elements = {{AtomicNumber{13}, AmuMass{27}, {}, "Al"}}; - inp.materials = {{2.7 * constants::na_avogadro / 27, + inp.materials = {{native_value_from(MolCcDensity{0.1}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, @@ -116,7 +117,7 @@ auto SimpleTestBase::build_physics() -> SPConstPhysics ImportPhysicsTable lambda; lambda.table_type = ImportTableType::lambda; lambda.x_units = ImportUnits::mev; - lambda.y_units = ImportUnits::cm_inv; + lambda.y_units = ImportUnits::len_inv; lambda.physics_vectors = { {ImportPhysicsVectorType::log, {1e-4, 1.0}, // energy @@ -131,7 +132,7 @@ auto SimpleTestBase::build_physics() -> SPConstPhysics ImportPhysicsTable lambdap; lambdap.table_type = ImportTableType::lambda_prim; lambdap.x_units = ImportUnits::mev; - lambdap.y_units = ImportUnits::cm_mev_inv; + lambdap.y_units = ImportUnits::len_mev_inv; lambdap.physics_vectors = { {ImportPhysicsVectorType::log, {1.0, 1e4, 1e8}, // energy @@ -143,6 +144,13 @@ auto SimpleTestBase::build_physics() -> SPConstPhysics compton_data.tables.push_back(std::move(lambdap)); } + // Update data values from CGS + { + celeritas::detail::ImportDataConverter convert{ + celeritas::UnitSystem::cgs}; + convert(&compton_data); + } + auto process_data = std::make_shared( std::vector{std::move(compton_data)}); diff --git a/test/celeritas/em/CombinedBrem.test.cc b/test/celeritas/em/CombinedBrem.test.cc index 586c89f134..0db4a1d922 100644 --- a/test/celeritas/em/CombinedBrem.test.cc +++ b/test/celeritas/em/CombinedBrem.test.cc @@ -47,7 +47,7 @@ class CombinedBremTest : public InteractorHostTestBase mat_inp.elements = {{AtomicNumber{29}, units::AmuMass{63.546}, {}, "Cu"}}; mat_inp.materials = { - {0.141 * na_avogadro, + {native_value_from(MolCcDensity{0.141}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, diff --git a/test/celeritas/em/EPlusGG.test.cc b/test/celeritas/em/EPlusGG.test.cc index 2759736739..70488a6681 100644 --- a/test/celeritas/em/EPlusGG.test.cc +++ b/test/celeritas/em/EPlusGG.test.cc @@ -246,7 +246,8 @@ TEST_F(EPlusGGInteractorTest, macro_xs) { real_type e = std::exp(loge); energy.push_back(e); - macro_xs.push_back(calc_macro_xs(MevEnergy(e))); + macro_xs.push_back( + native_value_to(calc_macro_xs(MevEnergy{e})).value()); loge += delta; } real_type const expected_macro_xs[] @@ -259,6 +260,7 @@ TEST_F(EPlusGGInteractorTest, macro_xs) 6.355265134801e-10, 2.068312058021e-10}; EXPECT_VEC_SOFT_EQ(expected_macro_xs, macro_xs); } + //---------------------------------------------------------------------------// } // namespace test } // namespace celeritas diff --git a/test/celeritas/em/Fluctuation.test.cc b/test/celeritas/em/Fluctuation.test.cc index 9b710c2710..1925b27879 100644 --- a/test/celeritas/em/Fluctuation.test.cc +++ b/test/celeritas/em/Fluctuation.test.cc @@ -7,6 +7,7 @@ //---------------------------------------------------------------------------// #include "corecel/data/CollectionStateStore.hh" #include "celeritas/MockTestBase.hh" +#include "celeritas/Quantities.hh" #include "celeritas/em/FluctuationParams.hh" #include "celeritas/em/distribution/EnergyLossDeltaDistribution.hh" #include "celeritas/em/distribution/EnergyLossGammaDistribution.hh" @@ -66,7 +67,7 @@ class EnergyLossDistributionTest : public Test MaterialParams::Input mat_inp; mat_inp.elements = {{AtomicNumber{18}, AmuMass{39.948}, {}, "Ar"}}; mat_inp.materials = { - {1.0 * na_avogadro, + {native_value_from(MolCcDensity{1.0}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, @@ -150,7 +151,7 @@ TEST_F(EnergyLossDistributionTest, none) MevEnergy mean_loss{2e-6}; // Tiny step, little energy loss - real_type step = 1e-6; + real_type step = 1e-6 * units::centimeter; EnergyLossHelper helper( fluct->host_ref(), cutoff, material, particle, mean_loss, step); EXPECT_EQ(EnergyLossFluctuationModel::none, helper.model()); @@ -181,7 +182,7 @@ TEST_F(EnergyLossDistributionTest, gaussian) // Larger step samples from gamma distribution, smaller step from Gaussian { real_type sum = 0; - real_type step = 5e-2; + real_type step = 5e-2 * units::centimeter; EnergyLossHelper helper( fluct->host_ref(), cutoff, material, particle, mean_loss, step); EXPECT_EQ(EnergyLossFluctuationModel::gamma, helper.model()); @@ -206,7 +207,7 @@ TEST_F(EnergyLossDistributionTest, gaussian) } { real_type sum = 0; - real_type step = 5e-4; + real_type step = 5e-4 * units::centimeter; EnergyLossHelper helper( fluct->host_ref(), cutoff, material, particle, mean_loss, step); EXPECT_SOFT_EQ(0.00019160444039613, @@ -244,7 +245,7 @@ TEST_F(EnergyLossDistributionTest, urban) material = {MaterialId{0}}; CutoffView cutoff(cutoffs->host_ref(), MaterialId{0}); MevEnergy mean_loss{0.01}; - real_type step = 0.01; + real_type step = 0.01 * units::centimeter; int num_samples = 10000; std::vector counts(20); diff --git a/test/celeritas/em/LivermorePE.test.cc b/test/celeritas/em/LivermorePE.test.cc index f430f970d2..ec04497a6f 100644 --- a/test/celeritas/em/LivermorePE.test.cc +++ b/test/celeritas/em/LivermorePE.test.cc @@ -61,7 +61,7 @@ class LivermorePETest : public InteractorHostTestBase // Set up shared material data MaterialParams::Input mi; mi.elements = {{AtomicNumber{19}, AmuMass{39.0983}, {}, "K"}}; - mi.materials = {{1e-5 * na_avogadro, + mi.materials = {{native_value_from(MolCcDensity{1e-5}), 293., MatterState::solid, {{ElementId{0}, 1.0}}, @@ -499,7 +499,8 @@ TEST_F(LivermorePETest, macro_xs) { real_type e = std::exp(loge); energy.push_back(e); - macro_xs.push_back(calc_macro_xs(MevEnergy{e})); + macro_xs.push_back( + native_value_to(calc_macro_xs(MevEnergy{e})).value()); loge += delta; } real_type const expected_macro_xs[] diff --git a/test/celeritas/em/MollerBhabha.test.cc b/test/celeritas/em/MollerBhabha.test.cc index ff9bf8c993..7bad9fe834 100644 --- a/test/celeritas/em/MollerBhabha.test.cc +++ b/test/celeritas/em/MollerBhabha.test.cc @@ -37,7 +37,7 @@ class MollerBhabhaInteractorTest : public InteractorHostTestBase MaterialParams::Input inp; inp.elements = {{AtomicNumber{29}, units::AmuMass{63.546}, {}, "Cu"}}; inp.materials = { - {1.0 * constants::na_avogadro, + {native_value_from(MolCcDensity{1.0}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, diff --git a/test/celeritas/em/Rayleigh.test.cc b/test/celeritas/em/Rayleigh.test.cc index 2b1dd94ae5..661296eaaf 100644 --- a/test/celeritas/em/Rayleigh.test.cc +++ b/test/celeritas/em/Rayleigh.test.cc @@ -49,7 +49,7 @@ class RayleighInteractorTest : public InteractorHostTestBase {AtomicNumber{74}, units::AmuMass{183.84}, {}, "W"}, {AtomicNumber{82}, units::AmuMass{207.2}, {}, "Pb"}}; inp.materials = { - {1.0 * constants::na_avogadro, + {native_value_from(MolCcDensity{1.0}), 293.0, MatterState::solid, {{ElementId{0}, 0.5}, {ElementId{1}, 0.3}, {ElementId{2}, 0.2}}, diff --git a/test/celeritas/em/RelativisticBrem.test.cc b/test/celeritas/em/RelativisticBrem.test.cc index 11e766cf54..9e7964af5f 100644 --- a/test/celeritas/em/RelativisticBrem.test.cc +++ b/test/celeritas/em/RelativisticBrem.test.cc @@ -35,10 +35,12 @@ class RelativisticBremTest : public InteractorHostTestBase protected: void SetUp() override { + using namespace units; + // Set up shared material data MaterialParams::Input mi; - mi.elements = {{AtomicNumber{82}, units::AmuMass{207.2}, {}, "Pb"}}; - mi.materials = {{0.05477 * constants::na_avogadro, + mi.elements = {{AtomicNumber{82}, AmuMass{207.2}, {}, "Pb"}}; + mi.materials = {{native_value_from(MolCcDensity{0.05477}), 293.15, MatterState::solid, {{ElementId{0}, 1.0}}, diff --git a/test/celeritas/em/SeltzerBerger.test.cc b/test/celeritas/em/SeltzerBerger.test.cc index d02b0791ab..d8fbe94a34 100644 --- a/test/celeritas/em/SeltzerBerger.test.cc +++ b/test/celeritas/em/SeltzerBerger.test.cc @@ -43,6 +43,8 @@ class SeltzerBergerTest : public InteractorHostTestBase void SetUp() override { + using namespace celeritas::units; + auto const& particles = *this->particle_params(); data_.ids.electron = particles.find(pdg::electron()); data_.ids.positron = particles.find(pdg::positron()); @@ -51,10 +53,9 @@ class SeltzerBergerTest : public InteractorHostTestBase // Set up shared material data MaterialParams::Input mat_inp; - mat_inp.elements - = {{AtomicNumber{29}, units::AmuMass{63.546}, {}, "Cu"}}; + mat_inp.elements = {{AtomicNumber{29}, AmuMass{63.546}, {}, "Cu"}}; mat_inp.materials = { - {0.141 * constants::na_avogadro, + {native_value_from(MolCcDensity{0.141}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, @@ -442,12 +443,13 @@ TEST_F(SeltzerBergerTest, stress_test) TEST_F(SeltzerBergerTest, positron_xs_corrector_edge_case) { + using namespace celeritas::units; // See https://github.com/celeritas-project/celeritas/issues/617 // Set up material data (only value used in this test is the atomic number) MaterialParams::Input mat_inp; - mat_inp.elements = {{AtomicNumber{26}, units::AmuMass{55.845}, {}, "Fe"}}; - mat_inp.materials = {{0.128 * constants::na_avogadro, + mat_inp.elements = {{AtomicNumber{26}, AmuMass{55.845}, {}, "Fe"}}; + mat_inp.materials = {{native_value_from(MolCcDensity{0.128}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, @@ -456,10 +458,10 @@ TEST_F(SeltzerBergerTest, positron_xs_corrector_edge_case) auto const material_params = std::make_shared(std::move(mat_inp)); - units::MevMass const positron_mass{0.51099890999999997}; - units::MevEnergy const min_gamma_energy{0.020822442086622296}; - units::MevEnergy const inc_energy{241.06427050865221}; - units::MevEnergy const sampled_energy{0.020822442732819097}; + MevMass const positron_mass{0.51099890999999997}; + MevEnergy const min_gamma_energy{0.020822442086622296}; + MevEnergy const inc_energy{241.06427050865221}; + MevEnergy const sampled_energy{0.020822442732819097}; SBPositronXsCorrector xs_corrector(positron_mass, material_params->get(ElementId{0}), min_gamma_energy, diff --git a/test/celeritas/em/UrbanMsc.test.cc b/test/celeritas/em/UrbanMsc.test.cc index 479428461b..f62ade10f7 100644 --- a/test/celeritas/em/UrbanMsc.test.cc +++ b/test/celeritas/em/UrbanMsc.test.cc @@ -13,6 +13,7 @@ #include "corecel/data/CollectionStateStore.hh" #include "corecel/grid/Interpolator.hh" #include "celeritas/RootTestBase.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/em/UrbanMscParams.hh" #include "celeritas/em/msc/detail/MscStepFromGeo.hh" #include "celeritas/em/msc/detail/MscStepToGeo.hh" @@ -43,6 +44,21 @@ namespace detail { namespace test { +//---------------------------------------------------------------------------// +struct InvCentimeter +{ + static CELER_CONSTEXPR_FUNCTION real_type value() + { + return 1 / units::centimeter; + } + static char const* label() { return "1/cm"; } +}; + +using InvCmAlpha = Quantity; +using celeritas::test::from_cm; +using celeritas::test::to_cm; +using units::MevEnergy; + //---------------------------------------------------------------------------// TEST(UrbanPositronCorrectorTest, all) { @@ -247,13 +263,16 @@ TEST_F(UrbanMscTest, helper) auto phys = this->make_phys_view(par, "G4_STAINLESS-STEEL"); UrbanMscHelper helper(msc_params_->host_ref(), par, phys); - EXPECT_SOFT_EQ(0.90681578657668238, phys.dedx_range()); - EXPECT_SOFT_EQ(1.0897296072933604, helper.calc_msc_mfp(MevEnergy{10.01})); - EXPECT_SOFT_EQ(0.90820266262324023, helper.calc_msc_mfp(MevEnergy{9.01})); - EXPECT_SOFT_EQ(11.039692548085707, - value_as(helper.calc_inverse_range(1.0))); + EXPECT_SOFT_EQ(0.90681578657668238, to_cm(phys.dedx_range())); + EXPECT_SOFT_EQ(1.0897296072933604, + to_cm(helper.calc_msc_mfp(MevEnergy{10.01}))); + EXPECT_SOFT_EQ(0.90820266262324023, + to_cm(helper.calc_msc_mfp(MevEnergy{9.01}))); + EXPECT_SOFT_EQ( + 11.039692548085707, + value_as(helper.calc_inverse_range(from_cm(1.0)))); EXPECT_SOFT_EQ(4.5491422239586035, - value_as(helper.calc_end_energy(0.5))); + value_as(helper.calc_end_energy(from_cm(0.5)))); } TEST_F(UrbanMscTest, step_conversion) @@ -381,7 +400,7 @@ TEST_F(UrbanMscTest, TEST_IF_CELERITAS_DOUBLE(msc_scattering)) // are bit-for-bit identical when range limits the step. The first three // steps are not limited by range constexpr double step_is_range = -1; - std::vector step = {0.00279169, 0.412343, 0.0376414}; + std::vector step = {0.00279169, 0.412343, 0.0376414}; // [cm] step.resize(nsamples, step_is_range); ASSERT_EQ(nsamples, step.size()); @@ -395,18 +414,19 @@ TEST_F(UrbanMscTest, TEST_IF_CELERITAS_DOUBLE(msc_scattering)) RandomEngine rng; // Input and helper data - std::vector pstep; - std::vector range; - std::vector lambda; + std::vector pstep; // [cm] + std::vector range; // [cm] + std::vector lambda; // [cm] // Step limit - std::vector tstep; - std::vector gstep; + std::vector tstep; // [cm] + std::vector gstep; // [cm] std::vector alpha; + std::vector msc_range_limit; // [cm] // Scatter std::vector angle; - std::vector displace; + std::vector displace; // [cm] std::vector action; // Total RNG count (we only sample once per particle/energy so count and @@ -416,30 +436,40 @@ TEST_F(UrbanMscTest, TEST_IF_CELERITAS_DOUBLE(msc_scattering)) auto const& msc_params = msc_params_->host_ref(); auto sample_one = [&](PDGNumber ptype, int i) { + real_type radius = from_cm(i * 2 - real_type(1e-4)); + auto par = this->make_par_view(ptype, MevEnergy{energy[i]}); auto phys = this->make_phys_view(par, "G4_STAINLESS-STEEL"); - auto geo = this->make_geo_view(/* r = */ i * 2 - real_type(1e-4)); + auto geo = this->make_geo_view(radius); MaterialView mat = this->material()->get(phys.material_id()); rng.reset_count(); UrbanMscHelper helper(msc_params, par, phys); - range.push_back(phys.dedx_range()); - lambda.push_back(helper.msc_mfp()); + range.push_back(to_cm(phys.dedx_range())); + lambda.push_back(to_cm(helper.msc_mfp())); - real_type this_pstep = step[i]; - if (this_pstep == step_is_range) - { - this_pstep = phys.dedx_range(); - } - pstep.push_back(this_pstep); + real_type const this_pstep = [i, &phys, &step] { + if (step[i] == step_is_range) + { + return phys.dedx_range(); + } + real_type const pstep = from_cm(step[i]); + CELER_VALIDATE(pstep <= phys.dedx_range(), + << "unit test input pstep=" << pstep + << " exceeds physics range " << phys.dedx_range()); + return pstep; + }(); + pstep.push_back(to_cm(this_pstep)); // Calculate physical step limit due to MSC real_type safety = geo.find_safety(); auto [true_path, displaced] = [&]() -> std::pair { + EXPECT_FALSE(phys.msc_range()); if (this_pstep < msc_params.params.limit_min_fix() || safety >= helper.max_step()) { // Small step or far from boundary + msc_range_limit.push_back(-1); return {this_pstep, false}; } UrbanMscStepLimit calc_limit(msc_params, @@ -451,9 +481,13 @@ TEST_F(UrbanMscTest, TEST_IF_CELERITAS_DOUBLE(msc_scattering)) safety, this_pstep); + // MSC range should be updated during construction of the limit + // calculator + msc_range_limit.push_back(to_cm(phys.msc_range().limit_min)); + return {calc_limit(rng), true}; }(); - tstep.push_back(true_path); + tstep.push_back(to_cm(true_path)); // Convert physical step limit to geometrical step MscStepToGeo calc_geom_path(msc_params, @@ -462,8 +496,8 @@ TEST_F(UrbanMscTest, TEST_IF_CELERITAS_DOUBLE(msc_scattering)) helper.msc_mfp(), phys.dedx_range()); auto gp = calc_geom_path(true_path); - gstep.push_back(gp.step); - alpha.push_back(gp.alpha); + gstep.push_back(to_cm(gp.step)); + alpha.push_back(native_value_to(gp.alpha).value()); MscStep step_result; step_result.true_path = true_path; @@ -481,7 +515,7 @@ TEST_F(UrbanMscTest, TEST_IF_CELERITAS_DOUBLE(msc_scattering)) ? sample_result.direction[0] : 0); displace.push_back(sample_result.action == Action::displaced - ? sample_result.displacement[0] + ? to_cm(sample_result.displacement[0]) : 0); action.push_back(sample_result.action == Action::displaced ? 'd' : sample_result.action == Action::scattered ? 's' @@ -535,6 +569,11 @@ TEST_F(UrbanMscTest, TEST_IF_CELERITAS_DOUBLE(msc_scattering)) 3.4500251375335, 12.793635146437, 0, 0, 1347.4250715939, 32089.171921536, 0, 1.7061753085921, 3.3439279763905, 12.693879333563, 0, 0, 1479.3077264327, 38389.697977001}; + static double const expected_msc_range_limit[] = {3.5686548881735e-05, + 4.0510166112733e-05, 4.0073894011965e-05, -1, 2.5270068437103e-05, + 1.0695397351015e-05, -1, -1, 1.6703727150203e-05, 2.0782727059196e-05, + 2.0774322573966e-05, -1, 1.3419878391705e-05, 5.6685934033174e-06, -1, + -1}; static double const expected_angle[] = {0.00031474130607916, -0.0179600098014372, -0.14560882721751, 0, -0.32650640360665, 0.013072020086723, 0, 0, 0.003112817663327, 0.385707466417037, @@ -556,6 +595,7 @@ TEST_F(UrbanMscTest, TEST_IF_CELERITAS_DOUBLE(msc_scattering)) EXPECT_VEC_SOFT_EQ(expected_tstep, tstep); EXPECT_VEC_SOFT_EQ(expected_gstep, gstep); EXPECT_VEC_SOFT_EQ(expected_alpha, alpha); + EXPECT_VEC_SOFT_EQ(expected_msc_range_limit, msc_range_limit); EXPECT_VEC_NEAR(expected_angle, angle, 2e-12); EXPECT_VEC_NEAR(expected_displace, displace, 1e-11); EXPECT_VEC_EQ(expected_action, action); diff --git a/test/celeritas/em/Wentzel.test.cc b/test/celeritas/em/Wentzel.test.cc index 034937b906..270a7130e4 100644 --- a/test/celeritas/em/Wentzel.test.cc +++ b/test/celeritas/em/Wentzel.test.cc @@ -29,44 +29,41 @@ class CoulombScatteringTest : public InteractorHostTestBase protected: void SetUp() override { - using namespace constants; + using namespace celeritas::units; + using constants::stable_decay_constant; + // Need to include protons constexpr units::MevMass emass{0.5109989461}; - ParticleParams::Input par_inp - = {{"electron", - pdg::electron(), - emass, - celeritas::units::ElementaryCharge{-1}, - stable_decay_constant}, - {"positron", - pdg::positron(), - emass, - celeritas::units::ElementaryCharge{1}, - stable_decay_constant}, - {"proton", - pdg::proton(), - units::MevMass{938.28}, - celeritas::units::ElementaryCharge{1}, - stable_decay_constant}}; + ParticleParams::Input par_inp = { + {"electron", + pdg::electron(), + emass, + ElementaryCharge{-1}, + stable_decay_constant}, + {"positron", + pdg::positron(), + emass, + ElementaryCharge{1}, + stable_decay_constant}, + {"proton", + pdg::proton(), + units::MevMass{938.28}, + ElementaryCharge{1}, + stable_decay_constant}, + }; this->set_particle_params(std::move(par_inp)); // Set up shared material data - // TODO: Use multiple elements to test elements are picked correctly MaterialParams::Input mat_inp; - mat_inp.isotopes = {{AtomicNumber{29}, - AtomicNumber{63}, - units::MevMass{58618.5}, - "63Cu"}, - {AtomicNumber{29}, - AtomicNumber{65}, - units::MevMass{60479.8}, - "65Cu"}}; + mat_inp.isotopes + = {{AtomicNumber{29}, AtomicNumber{63}, MevMass{58618.5}, "63Cu"}, + {AtomicNumber{29}, AtomicNumber{65}, MevMass{60479.8}, "65Cu"}}; mat_inp.elements = {{AtomicNumber{29}, - units::AmuMass{63.546}, + AmuMass{63.546}, {{IsotopeId{0}, 0.692}, {IsotopeId{1}, 0.308}}, "Cu"}}; mat_inp.materials = { - {0.141 * constants::na_avogadro, + {native_value_from(MolCcDensity{0.141}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, diff --git a/test/celeritas/ext/GeantGeo.test.cc b/test/celeritas/ext/GeantGeo.test.cc index c3ab3f9b6d..58087b3f7f 100644 --- a/test/celeritas/ext/GeantGeo.test.cc +++ b/test/celeritas/ext/GeantGeo.test.cc @@ -14,6 +14,7 @@ #include "corecel/io/StringUtils.hh" #include "celeritas/GenericGeoTestBase.hh" #include "celeritas/LazyGeoManager.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/ext/GeantGeoData.hh" #include "celeritas/ext/GeantGeoParams.hh" #include "celeritas/ext/GeantGeoTrackView.hh" @@ -79,8 +80,8 @@ TEST_F(FourLevelsTest, accessors) { auto const& geom = *this->geometry(); auto const& bbox = geom.bbox(); - EXPECT_VEC_SOFT_EQ((Real3{-24., -24., -24.}), bbox.lower()); - EXPECT_VEC_SOFT_EQ((Real3{24., 24., 24.}), bbox.upper()); + EXPECT_VEC_SOFT_EQ((Real3{-24., -24., -24.}), to_cm(bbox.lower())); + EXPECT_VEC_SOFT_EQ((Real3{24., 24., 24.}), to_cm(bbox.upper())); ASSERT_EQ(4, geom.num_volumes()); EXPECT_EQ("Shape2", geom.id_to_label(VolumeId{0}).name); @@ -103,13 +104,13 @@ TEST_F(FourLevelsTest, consecutive_compute) EXPECT_EQ(VolumeId{0}, geo.volume_id()); EXPECT_FALSE(geo.is_on_boundary()); - auto next = geo.find_next_step(10.0); - EXPECT_SOFT_EQ(4.0, next.distance); - EXPECT_SOFT_EQ(4.0, geo.find_safety()); + auto next = geo.find_next_step(from_cm(10.0)); + EXPECT_SOFT_EQ(4.0, to_cm(next.distance)); + EXPECT_SOFT_EQ(4.0, to_cm(geo.find_safety())); - next = geo.find_next_step(10.0); - EXPECT_SOFT_EQ(4.0, next.distance); - EXPECT_SOFT_EQ(4.0, geo.find_safety()); + next = geo.find_next_step(from_cm(10.0)); + EXPECT_SOFT_EQ(4.0, to_cm(next.distance)); + EXPECT_SOFT_EQ(4.0, to_cm(geo.find_safety())); } //---------------------------------------------------------------------------// @@ -124,18 +125,18 @@ TEST_F(FourLevelsTest, detailed_track) EXPECT_FALSE(geo.is_on_boundary()); // Check for surfaces up to a distance of 4 units away - auto next = geo.find_next_step(4.0); - EXPECT_SOFT_EQ(4.0, next.distance); + auto next = geo.find_next_step(from_cm(4.0)); + EXPECT_SOFT_EQ(4.0, to_cm(next.distance)); EXPECT_FALSE(next.boundary); - next = geo.find_next_step(4.0); - EXPECT_SOFT_EQ(4.0, next.distance); + next = geo.find_next_step(from_cm(4.0)); + EXPECT_SOFT_EQ(4.0, to_cm(next.distance)); EXPECT_FALSE(next.boundary); - geo.move_internal(3.5); + geo.move_internal(from_cm(3.5)); EXPECT_FALSE(geo.is_on_boundary()); // Find one a bit further, then cross it - next = geo.find_next_step(4.0); - EXPECT_SOFT_EQ(1.5, next.distance); + next = geo.find_next_step(from_cm(4.0)); + EXPECT_SOFT_EQ(1.5, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); EXPECT_EQ(VolumeId{0}, geo.volume_id()); @@ -146,11 +147,11 @@ TEST_F(FourLevelsTest, detailed_track) // Find the next boundary and make sure that nearer distances aren't // accepted next = geo.find_next_step(); - EXPECT_SOFT_EQ(1.0, next.distance); + EXPECT_SOFT_EQ(1.0, to_cm(next.distance)); EXPECT_TRUE(next.boundary); EXPECT_TRUE(geo.is_on_boundary()); - next = geo.find_next_step(0.5); - EXPECT_SOFT_EQ(0.5, next.distance); + next = geo.find_next_step(from_cm(0.5)); + EXPECT_SOFT_EQ(0.5, to_cm(next.distance)); EXPECT_FALSE(next.boundary); } { @@ -159,8 +160,8 @@ TEST_F(FourLevelsTest, detailed_track) EXPECT_FALSE(geo.is_outside()); EXPECT_EQ(VolumeId{3}, geo.volume_id()); - auto next = geo.find_next_step(2); - EXPECT_SOFT_EQ(0.5, next.distance); + auto next = geo.find_next_step(from_cm(2)); + EXPECT_SOFT_EQ(0.5, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); @@ -180,8 +181,8 @@ TEST_F(FourLevelsTest, reentrant_boundary) EXPECT_FALSE(geo.is_on_boundary()); // Check for surfaces: we should hit the outside of the sphere Shape2 - auto next = geo.find_next_step(1.0); - EXPECT_SOFT_EQ(0.5, next.distance); + auto next = geo.find_next_step(from_cm(1.0)); + EXPECT_SOFT_EQ(0.5, to_cm(next.distance)); // Move to the boundary but scatter perpendicularly, away from the sphere geo.move_to_boundary(); EXPECT_TRUE(geo.is_on_boundary()); @@ -190,14 +191,14 @@ TEST_F(FourLevelsTest, reentrant_boundary) EXPECT_EQ(VolumeId{1}, geo.volume_id()); // Move a bit internally, then scatter back toward the sphere - next = geo.find_next_step(10.0); - EXPECT_SOFT_EQ(6, next.distance); + next = geo.find_next_step(from_cm(10.0)); + EXPECT_SOFT_EQ(6, to_cm(next.distance)); geo.set_dir({-1, 0, 0}); EXPECT_EQ(VolumeId{1}, geo.volume_id()); // Move to the sphere boundary then scatter still into the sphere - next = geo.find_next_step(10.0); - EXPECT_SOFT_EQ(1e-13, next.distance); + next = geo.find_next_step(from_cm(10.0)); + EXPECT_SOFT_EQ(1e-13, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); EXPECT_TRUE(geo.is_on_boundary()); @@ -209,8 +210,8 @@ TEST_F(FourLevelsTest, reentrant_boundary) // Travel nearly tangent to the right edge of the sphere, then scatter to // still outside - next = geo.find_next_step(1.0); - EXPECT_SOFT_EQ(9.9794624025613538e-07, next.distance); + next = geo.find_next_step(from_cm(1.0)); + EXPECT_SOFT_EQ(9.9794624025613538e-07, to_cm(next.distance)); geo.move_to_boundary(); EXPECT_TRUE(geo.is_on_boundary()); geo.set_dir({1, 0, 0}); @@ -218,7 +219,7 @@ TEST_F(FourLevelsTest, reentrant_boundary) geo.cross_boundary(); EXPECT_EQ(VolumeId{1}, geo.volume_id()); EXPECT_TRUE(geo.is_on_boundary()); - next = geo.find_next_step(10.0); + next = geo.find_next_step(from_cm(10.0)); } //---------------------------------------------------------------------------// @@ -313,13 +314,13 @@ TEST_F(FourLevelsTest, safety) for (auto i : range(11)) { - real_type r = 2.0 * i + 0.1; + real_type r = from_cm(2.0 * i + 0.1); geo = {{r, r, r}, {1, 0, 0}}; if (!geo.is_outside()) { geo.find_next_step(); - safeties.push_back(geo.find_safety()); - lim_safeties.push_back(geo.find_safety(1.5)); + safeties.push_back(to_cm(geo.find_safety())); + lim_safeties.push_back(to_cm(geo.find_safety(from_cm(1.5)))); } } @@ -356,8 +357,8 @@ TEST_F(SolidsTest, accessors) { auto const& geom = *this->geometry(); auto const& bbox = geom.bbox(); - EXPECT_VEC_SOFT_EQ((Real3{-600., -300., -75.}), bbox.lower()); - EXPECT_VEC_SOFT_EQ((Real3{600., 300., 75.}), bbox.upper()); + EXPECT_VEC_SOFT_EQ((Real3{-600., -300., -75.}), to_cm(bbox.lower())); + EXPECT_VEC_SOFT_EQ((Real3{600., 300., 75.}), to_cm(bbox.upper())); // NOTE: because SolidsTest gets loaded after FourLevelsTest, the existing // volumes still have incremented the volume ID counter, so there is an @@ -377,13 +378,11 @@ TEST_F(SolidsTest, output) GeoParamsOutput out(this->geometry()); EXPECT_EQ("geometry", out.label()); - if (CELERITAS_USE_JSON) + if (CELERITAS_USE_JSON && CELERITAS_UNITS == CELERITAS_UNITS_CGS) { - EXPECT_EQ( + EXPECT_JSON_EQ( R"json({"bbox":[[-600.0,-300.0,-75.0],[600.0,300.0,75.0]],"supports_safety":true,"volumes":{"label":["","","","","box500","cone1","para1","sphere1","parabol1","trap1","trd1","trd2","","trd3_refl","tube100","boolean1","polycone1","genPocone1","ellipsoid1","tetrah1","orb1","polyhedr1","hype1","elltube1","ellcone1","arb8b","arb8a","xtru1","World","trd3_refl"]}})json", - to_string(out)) - << "\n/*** REPLACE ***/\nR\"json(" << to_string(out) - << ")json\"\n/******/"; + to_string(out)); } } diff --git a/test/celeritas/ext/GeantImporter.test.cc b/test/celeritas/ext/GeantImporter.test.cc index bdfc430b52..5773cf2dd7 100644 --- a/test/celeritas/ext/GeantImporter.test.cc +++ b/test/celeritas/ext/GeantImporter.test.cc @@ -14,6 +14,7 @@ #include "corecel/io/Repr.hh" #include "corecel/io/StringUtils.hh" #include "corecel/sys/Version.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/ext/GeantSetup.hh" #include "celeritas/io/ImportData.hh" #include "celeritas/phys/PDGNumber.hh" @@ -24,6 +25,8 @@ # include "celeritas/ext/GeantPhysicsOptionsIO.json.hh" #endif +using namespace celeritas::units; + namespace celeritas { namespace test @@ -49,6 +52,11 @@ std::string sub_pointer_string(std::string const& s) return std::regex_replace(s, r, "0x0"); } +double to_inv_cm(double v) +{ + return native_value_to(v).value(); +} + auto const geant4_version = Version::from_string(celeritas_geant4_version); } // namespace @@ -119,10 +127,21 @@ class GeantImporterTest : public GeantTestBase << " for particle PDG=" << pdg.get()); return *result; } + real_type comparison_tolerance() const { - // Some values change substantially between geant versions - return geant4_version == Version(11, 0, 3) ? 1e-12 : 5e-3; + if (geant4_version != Version(11, 0, 3)) + { + // Some values change substantially between geant versions + return 5e-3; + } + if (CELERITAS_REAL_TYPE != CELERITAS_REAL_TYPE_DOUBLE) + { + // Single-precision unit constants cause single-precision + // differences from reference + return 1e-6; + } + return 1e-12; } protected: @@ -233,6 +252,7 @@ class FourSteelSlabsEmStandard : public GeantImporterTest opts.relaxation = RelaxationSelection::all; opts.verbose = true; #if CELERITAS_USE_JSON + if (CELERITAS_UNITS == CELERITAS_UNITS_CGS) { nlohmann::json out = opts; static char const expected[] @@ -521,14 +541,15 @@ TEST_F(FourSteelSlabsEmStandard, materials) { names.push_back(material.name); states.push_back(static_cast(material.state)); - num_densities.push_back(material.number_density); + num_densities.push_back( + native_value_to(material.number_density).value()); temperatures.push_back(material.temperature); for (auto const& key : material.pdg_cutoffs) { pdgs.push_back(key.first); cutoff_energies.push_back(key.second.energy); - cutoff_ranges.push_back(key.second.range); + cutoff_ranges.push_back(to_cm(key.second.range)); } for (auto const& el_comp : material.elements) @@ -538,6 +559,8 @@ TEST_F(FourSteelSlabsEmStandard, materials) } } + real_type const tol = this->comparison_tolerance(); + static char const* expected_names[] = {"G4_Galactic", "G4_STAINLESS-STEEL"}; EXPECT_VEC_EQ(expected_names, names); static int const expected_states[] = {3, 1}; @@ -555,10 +578,10 @@ TEST_F(FourSteelSlabsEmStandard, materials) geant4_version.major() == 10 ? 1e-12 : 0.02); static double const expected_cutoff_ranges[] = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1}; - EXPECT_VEC_SOFT_EQ(expected_cutoff_ranges, cutoff_ranges); + EXPECT_VEC_NEAR(expected_cutoff_ranges, cutoff_ranges, tol); static double const expected_num_densities[] = {0.05974697167543, 8.699348925899e+22}; - EXPECT_VEC_SOFT_EQ(expected_num_densities, num_densities); + EXPECT_VEC_NEAR(expected_num_densities, num_densities, tol); static double const expected_temperatures[] = {2.73, 293.15}; EXPECT_VEC_SOFT_EQ(expected_temperatures, temperatures); static double const expected_el_comps_ids[] = {3, 0, 1, 2}; @@ -605,8 +628,8 @@ TEST_F(FourSteelSlabsEmStandard, eioni) ASSERT_EQ(85, steel.x.size()); EXPECT_SOFT_EQ(1e-4, steel.x.front()); EXPECT_SOFT_EQ(1e8, steel.x.back()); - EXPECT_SOFT_NEAR(839.66835335480653, steel.y.front(), tol); - EXPECT_SOFT_NEAR(11.378226755591747, steel.y.back(), tol); + EXPECT_SOFT_NEAR(839.66835335480653, to_inv_cm(steel.y.front()), tol); + EXPECT_SOFT_NEAR(11.378226755591747, to_inv_cm(steel.y.back()), tol); } { // Test range table @@ -622,8 +645,8 @@ TEST_F(FourSteelSlabsEmStandard, eioni) ASSERT_EQ(85, steel.x.size()); EXPECT_SOFT_EQ(1e-4, steel.x.front()); EXPECT_SOFT_EQ(1e8, steel.x.back()); - EXPECT_SOFT_NEAR(2.3818927937550707e-07, steel.y.front(), tol); - EXPECT_SOFT_NEAR(8788715.7877501156, steel.y.back(), tol); + EXPECT_SOFT_NEAR(2.3818927937550707e-07, to_cm(steel.y.front()), tol); + EXPECT_SOFT_NEAR(8788715.7877501156, to_cm(steel.y.back()), tol); } { // Test cross-section table @@ -640,8 +663,8 @@ TEST_F(FourSteelSlabsEmStandard, eioni) EXPECT_SOFT_NEAR(2.616556310615175, steel.x.front(), tol); EXPECT_SOFT_EQ(1e8, steel.x.back()); EXPECT_SOFT_EQ(0, steel.y.front()); - EXPECT_SOFT_NEAR(0.1905939505829807, steel.y[1], tol); - EXPECT_SOFT_NEAR(0.4373910150880348, steel.y.back(), tol); + EXPECT_SOFT_NEAR(0.1905939505829807, to_inv_cm(steel.y[1]), tol); + EXPECT_SOFT_NEAR(0.4373910150880348, to_inv_cm(steel.y.back()), tol); } } @@ -774,9 +797,9 @@ TEST_F(FourSteelSlabsEmStandard, em_parameters) EXPECT_DOUBLE_EQ(0.01, em_params.linear_loss_limit); EXPECT_DOUBLE_EQ(0.001, em_params.lowest_electron_energy); EXPECT_EQ(true, em_params.auger); - EXPECT_EQ(0.04, em_params.msc_range_factor); - EXPECT_EQ(0.6, em_params.msc_safety_factor); - EXPECT_EQ(0.1, em_params.msc_lambda_limit); + EXPECT_DOUBLE_EQ(0.04, em_params.msc_range_factor); + EXPECT_DOUBLE_EQ(0.6, em_params.msc_safety_factor); + EXPECT_REAL_EQ(0.1, to_cm(em_params.msc_lambda_limit)); } //---------------------------------------------------------------------------// @@ -1111,14 +1134,14 @@ TEST_F(OneSteelSphere, cutoffs) for (auto const& cut : mat.pdg_cutoffs) { pdg.push_back(cut.first); - range_cut.push_back(cut.second.range); + range_cut.push_back(to_cm(cut.second.range)); energy_cut.push_back(cut.second.energy); } } static int const expected_pdg[] = {-11, 11, 22, -11, 11, 22}; EXPECT_VEC_EQ(expected_pdg, pdg); // 1 mm range cut in vacuum, 50 m range cut in steel - static double const expected_range_cut[] + static real_type const expected_range_cut[] = {0.1, 0.1, 0.1, 5000, 5000, 5000}; EXPECT_VEC_SOFT_EQ(expected_range_cut, range_cut); static double const expected_energy_cut[] = {0.00099, diff --git a/test/celeritas/ext/Vecgeom.test.cc b/test/celeritas/ext/Vecgeom.test.cc index 951267b3bd..c4d10366b9 100644 --- a/test/celeritas/ext/Vecgeom.test.cc +++ b/test/celeritas/ext/Vecgeom.test.cc @@ -20,6 +20,7 @@ #include "corecel/sys/Environment.hh" #include "corecel/sys/Version.hh" #include "celeritas/GenericGeoTestBase.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/ext/GeantGeoUtils.hh" #include "celeritas/ext/GeantSetup.hh" #include "celeritas/ext/VecgeomData.hh" @@ -156,39 +157,39 @@ TEST_F(SimpleCmsTest, track) auto geo = this->make_geo_track_view({0, 0, 0}, {0, 0, 1}); EXPECT_EQ("vacuum_tube", this->volume_name(geo)); - auto next = geo.find_next_step(100); - EXPECT_SOFT_EQ(100, next.distance); + auto next = geo.find_next_step(from_cm(100)); + EXPECT_SOFT_EQ(100, to_cm(next.distance)); EXPECT_FALSE(next.boundary); - geo.move_internal(20); - EXPECT_SOFT_EQ(30, geo.find_safety()); + geo.move_internal(from_cm(20)); + EXPECT_SOFT_EQ(30, to_cm(geo.find_safety())); geo.set_dir({1, 0, 0}); - next = geo.find_next_step(50); - EXPECT_SOFT_EQ(30, next.distance); + next = geo.find_next_step(from_cm(50)); + EXPECT_SOFT_EQ(30, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); EXPECT_FALSE(geo.is_outside()); geo.cross_boundary(); EXPECT_EQ("si_tracker", this->volume_name(geo)); - EXPECT_VEC_SOFT_EQ(Real3({30, 0, 20}), geo.pos()); + EXPECT_VEC_SOFT_EQ(Real3({30, 0, 20}), to_cm(geo.pos())); // Scatter to tangent geo.set_dir({0, 1, 0}); - next = geo.find_next_step(1000); - EXPECT_SOFT_EQ(121.34661099511597, next.distance); + next = geo.find_next_step(from_cm(1000)); + EXPECT_SOFT_EQ(121.34661099511597, to_cm(next.distance)); EXPECT_TRUE(next.boundary); - geo.move_internal(10); - EXPECT_SOFT_EQ(1.6227766016837926, geo.find_safety()); + geo.move_internal(from_cm(10)); + EXPECT_SOFT_EQ(1.6227766016837926, to_cm(geo.find_safety())); // Move to boundary and scatter back inside - next = geo.find_next_step(1000); - EXPECT_SOFT_EQ(111.34661099511597, next.distance); + next = geo.find_next_step(from_cm(1000)); + EXPECT_SOFT_EQ(111.34661099511597, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); geo.set_dir({-1, 0, 0}); - next = geo.find_next_step(1000); - EXPECT_SOFT_EQ(60., next.distance); + next = geo.find_next_step(from_cm(1000)); + EXPECT_SOFT_EQ(60., to_cm(next.distance)); } //---------------------------------------------------------------------------// @@ -301,8 +302,8 @@ TEST_F(FourLevelsTest, accessors) EXPECT_EQ(4, geom.max_depth()); auto const& bbox = geom.bbox(); - EXPECT_VEC_SOFT_EQ((Real3{-24.001, -24.001, -24.001}), bbox.lower()); - EXPECT_VEC_SOFT_EQ((Real3{24.001, 24.001, 24.001}), bbox.upper()); + EXPECT_VEC_SOFT_EQ((Real3{-24.001, -24.001, -24.001}), to_cm(bbox.lower())); + EXPECT_VEC_SOFT_EQ((Real3{24.001, 24.001, 24.001}), to_cm(bbox.upper())); ASSERT_EQ(4, geom.num_volumes()); EXPECT_EQ("Shape2", geom.id_to_label(VolumeId{0}).name); @@ -321,13 +322,13 @@ TEST_F(FourLevelsTest, consecutive_compute) EXPECT_EQ(VolumeId{0}, geo.volume_id()); EXPECT_FALSE(geo.is_on_boundary()); - auto next = geo.find_next_step(10.0); - EXPECT_SOFT_EQ(4.0, next.distance); - EXPECT_SOFT_EQ(4.0, geo.find_safety()); + auto next = geo.find_next_step(from_cm(10.0)); + EXPECT_SOFT_EQ(4.0, to_cm(next.distance)); + EXPECT_SOFT_EQ(4.0, to_cm(geo.find_safety())); - next = geo.find_next_step(10.0); - EXPECT_SOFT_EQ(4.0, next.distance); - EXPECT_SOFT_EQ(4.0, geo.find_safety()); + next = geo.find_next_step(from_cm(10.0)); + EXPECT_SOFT_EQ(4.0, to_cm(next.distance)); + EXPECT_SOFT_EQ(4.0, to_cm(geo.find_safety())); } //---------------------------------------------------------------------------// @@ -342,18 +343,18 @@ TEST_F(FourLevelsTest, detailed_track) EXPECT_FALSE(geo.is_on_boundary()); // Check for surfaces up to a distance of 4 units away - auto next = geo.find_next_step(4.0); - EXPECT_SOFT_EQ(4.0, next.distance); + auto next = geo.find_next_step(from_cm(4.0)); + EXPECT_SOFT_EQ(4.0, to_cm(next.distance)); EXPECT_FALSE(next.boundary); - next = geo.find_next_step(4.0); - EXPECT_SOFT_EQ(4.0, next.distance); + next = geo.find_next_step(from_cm(4.0)); + EXPECT_SOFT_EQ(4.0, to_cm(next.distance)); EXPECT_FALSE(next.boundary); - geo.move_internal(3.5); + geo.move_internal(from_cm(3.5)); EXPECT_FALSE(geo.is_on_boundary()); // Find one a bit further, then cross it - next = geo.find_next_step(4.0); - EXPECT_SOFT_EQ(1.5, next.distance); + next = geo.find_next_step(from_cm(4.0)); + EXPECT_SOFT_EQ(1.5, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); EXPECT_EQ(VolumeId{0}, geo.volume_id()); @@ -364,11 +365,11 @@ TEST_F(FourLevelsTest, detailed_track) // Find the next boundary and make sure that nearer distances aren't // accepted next = geo.find_next_step(); - EXPECT_SOFT_EQ(1.0, next.distance); + EXPECT_SOFT_EQ(1.0, to_cm(next.distance)); EXPECT_TRUE(next.boundary); EXPECT_TRUE(geo.is_on_boundary()); - next = geo.find_next_step(0.5); - EXPECT_SOFT_EQ(0.5, next.distance); + next = geo.find_next_step(from_cm(0.5)); + EXPECT_SOFT_EQ(0.5, to_cm(next.distance)); EXPECT_FALSE(next.boundary); } { @@ -377,7 +378,7 @@ TEST_F(FourLevelsTest, detailed_track) EXPECT_TRUE(geo.is_outside()); auto next = geo.find_next_step(); - EXPECT_SOFT_EQ(1.0, next.distance); + EXPECT_SOFT_EQ(1.0, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); @@ -392,8 +393,8 @@ TEST_F(FourLevelsTest, detailed_track) EXPECT_FALSE(geo.is_outside()); EXPECT_EQ(VolumeId{3}, geo.volume_id()); - auto next = geo.find_next_step(2); - EXPECT_SOFT_EQ(0.5, next.distance); + auto next = geo.find_next_step(from_cm(2)); + EXPECT_SOFT_EQ(0.5, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); @@ -411,52 +412,52 @@ TEST_F(FourLevelsTest, detailed_track) TEST_F(FourLevelsTest, reentrant_boundary) { - auto geo = this->make_geo_track_view(); - geo = GeoTrackInitializer{{15.5, 10, 10}, {-1, 0, 0}}; + auto geo = this->make_geo_track_view({15.5, 10, 10}, {-1, 0, 0}); ASSERT_FALSE(geo.is_outside()); - EXPECT_EQ(VolumeId{1}, geo.volume_id()); + EXPECT_EQ("Shape1", this->volume_name(geo)); EXPECT_FALSE(geo.is_on_boundary()); // Check for surfaces: we should hit the outside of the sphere Shape2 - auto next = geo.find_next_step(1.0); - EXPECT_SOFT_EQ(0.5, next.distance); + auto next = geo.find_next_step(from_cm(1.0)); + EXPECT_SOFT_EQ(0.5, to_cm(next.distance)); // Move to the boundary but scatter perpendicularly, away from the sphere geo.move_to_boundary(); EXPECT_TRUE(geo.is_on_boundary()); geo.set_dir({0, 1, 0}); EXPECT_TRUE(geo.is_on_boundary()); - EXPECT_EQ(VolumeId{1}, geo.volume_id()); + EXPECT_EQ("Shape1", this->volume_name(geo)); // Move a bit internally, then scatter back toward the sphere - next = geo.find_next_step(10.0); - EXPECT_SOFT_EQ(6, next.distance); + next = geo.find_next_step(from_cm(10.0)); + EXPECT_SOFT_EQ(6, to_cm(next.distance)); geo.set_dir({-1, 0, 0}); - EXPECT_EQ(VolumeId{1}, geo.volume_id()); + EXPECT_EQ("Shape1", this->volume_name(geo)); // Move to the sphere boundary then scatter still into the sphere - next = geo.find_next_step(10.0); - EXPECT_SOFT_EQ(1e-8, next.distance); + next = geo.find_next_step(from_cm(10.0)); + EXPECT_SOFT_EQ(1e-8, to_cm(next.distance)); EXPECT_TRUE(next.boundary); geo.move_to_boundary(); EXPECT_TRUE(geo.is_on_boundary()); geo.set_dir({0, -1, 0}); EXPECT_TRUE(geo.is_on_boundary()); geo.cross_boundary(); - EXPECT_EQ(VolumeId{0}, geo.volume_id()); + EXPECT_EQ("Shape2", this->volume_name(geo)); EXPECT_TRUE(geo.is_on_boundary()); // Travel nearly tangent to the right edge of the sphere, then scatter to // still outside - next = geo.find_next_step(1.0); - EXPECT_SOFT_EQ(0.00031622777925735285, next.distance); + next = geo.find_next_step(from_cm(1.0)); + EXPECT_SOFT_EQ(0.00031622777925735285, to_cm(next.distance)); geo.move_to_boundary(); EXPECT_TRUE(geo.is_on_boundary()); geo.set_dir({1, 0, 0}); EXPECT_TRUE(geo.is_on_boundary()); geo.cross_boundary(); - EXPECT_EQ(VolumeId{1}, geo.volume_id()); + EXPECT_EQ("Shape1", this->volume_name(geo)); + EXPECT_TRUE(geo.is_on_boundary()); - next = geo.find_next_step(10.0); + next = geo.find_next_step(from_cm(10.0)); } //---------------------------------------------------------------------------// @@ -552,8 +553,8 @@ TEST_F(FourLevelsTest, safety) if (!geo.is_outside()) { geo.find_next_step(); - safeties.push_back(geo.find_safety()); - lim_safeties.push_back(geo.find_safety(1.5)); + safeties.push_back(to_cm(geo.find_safety())); + lim_safeties.push_back(to_cm(geo.find_safety(from_cm(1.5)))); } } @@ -671,8 +672,9 @@ TEST_F(SolidsTest, accessors) } auto const& bbox = geom.bbox(); - EXPECT_VEC_SOFT_EQ((Real3{-600.001, -300.001, -75.001}), bbox.lower()); - EXPECT_VEC_SOFT_EQ((Real3{600.001, 300.001, 75.001}), bbox.upper()); + EXPECT_VEC_SOFT_EQ((Real3{-600.001, -300.001, -75.001}), + to_cm(bbox.lower())); + EXPECT_VEC_SOFT_EQ((Real3{600.001, 300.001, 75.001}), to_cm(bbox.upper())); } //---------------------------------------------------------------------------// @@ -703,7 +705,7 @@ TEST_F(SolidsTest, output) GeoParamsOutput out(this->geometry()); EXPECT_EQ("geometry", out.label()); - if (CELERITAS_USE_JSON) + if (CELERITAS_USE_JSON && CELERITAS_UNITS == CELERITAS_UNITS_CGS) { auto out_str = simplify_pointers(to_string(out)); @@ -1080,8 +1082,8 @@ TEST_F(FourLevelsGeantTest, accessors) EXPECT_EQ(4, geom.max_depth()); auto const& bbox = geom.bbox(); - EXPECT_VEC_SOFT_EQ((Real3{-24.001, -24.001, -24.001}), bbox.lower()); - EXPECT_VEC_SOFT_EQ((Real3{24.001, 24.001, 24.001}), bbox.upper()); + EXPECT_VEC_SOFT_EQ((Real3{-24.001, -24.001, -24.001}), to_cm(bbox.lower())); + EXPECT_VEC_SOFT_EQ((Real3{24.001, 24.001, 24.001}), to_cm(bbox.upper())); ASSERT_EQ(4, geom.num_volumes()); EXPECT_EQ("Shape2", geom.id_to_label(VolumeId{0}).name); @@ -1236,8 +1238,9 @@ TEST_F(SolidsGeantTest, accessors) } auto const& bbox = geom.bbox(); - EXPECT_VEC_SOFT_EQ((Real3{-600.001, -300.001, -75.001}), bbox.lower()); - EXPECT_VEC_SOFT_EQ((Real3{600.001, 300.001, 75.001}), bbox.upper()); + EXPECT_VEC_SOFT_EQ((Real3{-600.001, -300.001, -75.001}), + to_cm(bbox.lower())); + EXPECT_VEC_SOFT_EQ((Real3{600.001, 300.001, 75.001}), to_cm(bbox.upper())); } //---------------------------------------------------------------------------// @@ -1269,7 +1272,7 @@ TEST_F(SolidsGeantTest, output) GeoParamsOutput out(this->geometry()); EXPECT_EQ("geometry", out.label()); - if (CELERITAS_USE_JSON) + if (CELERITAS_USE_JSON && CELERITAS_UNITS == CELERITAS_UNITS_CGS) { auto out_str = simplify_pointers(to_string(out)); @@ -1497,6 +1500,33 @@ TEST_F(ZnenvGeantTest, trace) } } +//---------------------------------------------------------------------------// +// ARBITRARY INPUT TESTS +//---------------------------------------------------------------------------// + +#define ArbitraryVgdmlTest DISABLED_ArbitraryVgdmlTest +class ArbitraryVgdmlTest : public VecgeomTestBase +{ + public: + SPConstGeo build_geometry() final + { + auto filename = celeritas::getenv("GDML"); + CELER_VALIDATE(!filename.empty(), + << "Set the 'GDML' environment variable and run this " + "test with " + "--gtest_filter=*ArbitraryVgdmlTest* " + "--gtest_also_run_disabled_tests"); + return std::make_shared(filename); + } +}; + +TEST_F(ArbitraryVgdmlTest, dump) +{ + this->geometry(); + auto const* world = vecgeom::GeoManager::Instance().GetWorld(); + world->PrintContent(); +} + //---------------------------------------------------------------------------// #define ArbitraryGeantTest DISABLED_ArbitraryGeantTest diff --git a/test/celeritas/field/CMSParameterizedField.hh b/test/celeritas/field/CMSParameterizedField.hh index 2c83028e82..81bd4ccfa4 100644 --- a/test/celeritas/field/CMSParameterizedField.hh +++ b/test/celeritas/field/CMSParameterizedField.hh @@ -33,7 +33,6 @@ class CMSParameterizedField public: //!@{ //! \name Type aliases - using Real3 = Array; using Real4 = Array; //!@} diff --git a/test/celeritas/field/LinearPropagator.test.cc b/test/celeritas/field/LinearPropagator.test.cc index 691d71a904..209c074f51 100644 --- a/test/celeritas/field/LinearPropagator.test.cc +++ b/test/celeritas/field/LinearPropagator.test.cc @@ -12,6 +12,7 @@ #include "corecel/io/Logger.hh" #include "corecel/sys/Device.hh" #include "celeritas/AllGeoTypedTestBase.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas_test.hh" @@ -30,6 +31,15 @@ class LinearPropagatorTest : public AllGeoTypedTestBase using SPConstGeo = typename GenericGeoTestBase::SPConstGeo; using GeoTrackView = typename GenericGeoTestBase::GeoTrackView; + void SetUp() override + { + if (CELERITAS_UNITS != CELERITAS_UNITS_CGS + && this->geo_name() == "ORANGE") + { + GTEST_SKIP() << "ORANGE currently requires CGS [cm]"; + } + } + std::string geometry_basename() const final { return "simple-cms"; } }; @@ -49,11 +59,12 @@ TYPED_TEST(LinearPropagatorTest, rvalue_type) this->make_geo_track_view({0, 0, 0}, {0, 0, 1})); EXPECT_TRUE(( std::is_same_v>)); - Propagation result = propagate(10); - EXPECT_SOFT_EQ(10, result.distance); + Propagation result = propagate(from_cm(10)); + EXPECT_SOFT_EQ(10, to_cm(result.distance)); EXPECT_FALSE(result.boundary); } - EXPECT_VEC_SOFT_EQ(Real3({0, 0, 10}), this->make_geo_track_view().pos()); + EXPECT_VEC_SOFT_EQ(Real3({0, 0, 10}), + to_cm(this->make_geo_track_view().pos())); } TYPED_TEST(LinearPropagatorTest, simple_cms) @@ -69,64 +80,64 @@ TYPED_TEST(LinearPropagatorTest, simple_cms) LinearPropagator>)); // Move up to a small distance - Propagation result = propagate(20); - EXPECT_SOFT_EQ(20, result.distance); + Propagation result = propagate(from_cm(20)); + EXPECT_SOFT_EQ(20, to_cm(result.distance)); EXPECT_FALSE(result.boundary); } // Check state and scatter - EXPECT_VEC_SOFT_EQ(Real3({0, 0, 20}), geo.pos()); + EXPECT_VEC_SOFT_EQ(Real3({0, 0, 20}), to_cm(geo.pos())); EXPECT_EQ("vacuum_tube", this->volume_name(geo)); geo.set_dir({1, 0, 0}); { LinearPropagator propagate(geo); - // Move to the next layer - Propagation result = propagate(1e20); - EXPECT_SOFT_EQ(30, result.distance); + // Move to the result layer + Propagation result = propagate(from_cm(1e20)); + EXPECT_SOFT_EQ(30, to_cm(result.distance)); EXPECT_TRUE(result.boundary); geo.cross_boundary(); } // Check state - EXPECT_VEC_SOFT_EQ(Real3({30, 0, 20}), geo.pos()); + EXPECT_VEC_SOFT_EQ(Real3({30, 0, 20}), to_cm(geo.pos())); EXPECT_EQ("si_tracker", this->volume_name(geo)); { LinearPropagator propagate(geo); // Move two steps internally - Propagation result = propagate(35); - EXPECT_SOFT_EQ(35, result.distance); + Propagation result = propagate(from_cm(35)); + EXPECT_SOFT_EQ(35, to_cm(result.distance)); EXPECT_FALSE(result.boundary); - result = propagate(40); - EXPECT_SOFT_EQ(40, result.distance); + result = propagate(from_cm(40)); + EXPECT_SOFT_EQ(40, to_cm(result.distance)); EXPECT_FALSE(result.boundary); } // Check state - EXPECT_VEC_SOFT_EQ(Real3({105, 0, 20}), geo.pos()); + EXPECT_VEC_SOFT_EQ(Real3({105, 0, 20}), to_cm(geo.pos())); EXPECT_EQ("si_tracker", this->volume_name(geo)); { LinearPropagator propagate(geo); - // Move to next boundary (infinite max distance) + // Move to result boundary (infinite max distance) Propagation result = propagate(); - EXPECT_SOFT_EQ(20, result.distance); + EXPECT_SOFT_EQ(20, to_cm(result.distance)); EXPECT_TRUE(result.boundary); geo.cross_boundary(); - // Move slightly inside before next scatter - result = propagate(0.1); - EXPECT_SOFT_EQ(0.1, result.distance); + // Move slightly inside before result scatter + result = propagate(from_cm(0.1)); + EXPECT_SOFT_EQ(0.1, to_cm(result.distance)); EXPECT_FALSE(result.boundary); } // Check state and scatter - EXPECT_VEC_SOFT_EQ(Real3({125.1, 0, 20}), geo.pos()); + EXPECT_VEC_SOFT_EQ(Real3({125.1, 0, 20}), to_cm(geo.pos())); EXPECT_EQ("em_calorimeter", this->volume_name(geo)); geo.set_dir({0, 0, -1}); @@ -134,19 +145,19 @@ TYPED_TEST(LinearPropagatorTest, simple_cms) LinearPropagator propagate(geo); // Move to world volume - Propagation result = propagate(10000); - EXPECT_SOFT_EQ(720, result.distance); + Propagation result = propagate(from_cm(10000)); + EXPECT_SOFT_EQ(720, to_cm(result.distance)); EXPECT_TRUE(result.boundary); geo.cross_boundary(); // Move outside - result = propagate(10000); - EXPECT_SOFT_EQ(1300, result.distance); + result = propagate(from_cm(10000)); + EXPECT_SOFT_EQ(1300, to_cm(result.distance)); EXPECT_TRUE(result.boundary); geo.cross_boundary(); } - EXPECT_VEC_SOFT_EQ(Real3({125.1, 0, -2000}), geo.pos()); + EXPECT_VEC_SOFT_EQ(Real3({125.1, 0, -2000}), to_cm(geo.pos())); EXPECT_EQ("[OUTSIDE]", this->volume_name(geo)); } diff --git a/test/celeritas/field/MagFieldEquation.test.cc b/test/celeritas/field/MagFieldEquation.test.cc index a749c389a8..4541abec36 100644 --- a/test/celeritas/field/MagFieldEquation.test.cc +++ b/test/celeritas/field/MagFieldEquation.test.cc @@ -7,6 +7,8 @@ //---------------------------------------------------------------------------// #include "celeritas/field/MagFieldEquation.hh" +#include "celeritas/UnitUtils.hh" + #include "celeritas_test.hh" using celeritas::units::ElementaryCharge; @@ -21,9 +23,9 @@ Real3 dummy_field(Real3 const& pos) { // Rotate and scale Real3 result; - result[0] = real_type(0.5) * pos[1]; - result[1] = real_type(1.0) * pos[2]; - result[2] = real_type(2.0) * pos[0]; + result[0] = real_type(0.5) * to_cm(pos[1]) * units::gauss; + result[1] = real_type(1.0) * to_cm(pos[2]) * units::gauss; + result[2] = real_type(2.0) * to_cm(pos[0]) * units::gauss; return result; } @@ -35,6 +37,14 @@ make_equation(FieldT&& field, ElementaryCharge charge) return Equation_t{::celeritas::forward(field), charge}; } +OdeState make_state(Real3 const& pos, Real3 const& mom) +{ + OdeState result; + result.pos = from_cm(pos); + result.mom = mom; + return result; +} + void print_expected(OdeState const& s) { cout << "/*** BEGIN CODE ***/\n" @@ -47,26 +57,32 @@ void print_expected(OdeState const& s) TEST(MagFieldEquationTest, charged) { + constexpr auto inv_cm = 1 / units::centimeter; + + // NOTE: result.pos (dx / ds) is unitless, + // and result.dir (dp / ds) has units of 1/length auto eval = make_equation(dummy_field, ElementaryCharge{3}); { - OdeState result = eval({{1, 2, 3}, {0, 0, 1}}); + OdeState result = eval(make_state({1, 2, 3}, {0, 0, 1})); EXPECT_VEC_SOFT_EQ(Real3({0, 0, 1}), result.pos); - EXPECT_VEC_SOFT_EQ(Real3({-0.002698132122, 0.000899377374, 0}), - result.mom); + EXPECT_VEC_SOFT_EQ( + Real3({-0.002698132122 * inv_cm, 0.000899377374 * inv_cm, 0}), + result.mom); } { - OdeState result = eval({{0.5, -2, -1}, {1, 2, 3}}); + OdeState result = eval(make_state({0.5, -2, -1}, {1, 2, 3})); EXPECT_VEC_SOFT_EQ( Real3({0.26726124191242, 0.53452248382485, 0.80178372573727}), result.pos); - EXPECT_VEC_SOFT_EQ(Real3({0.0012018435696159, - -0.0009614748556927, - 0.00024036871392318}), + EXPECT_VEC_SOFT_EQ(Real3({0.0012018435696159 * inv_cm, + -0.0009614748556927 * inv_cm, + 0.00024036871392318 * inv_cm}), result.mom); } if (CELERITAS_DEBUG) { - EXPECT_THROW(eval({{0.5, -2, -1}, {0, 0, 0}}), DebugError); + EXPECT_THROW(eval(make_state(Real3{0.5, -2, -1}, {0, 0, 0})), + DebugError); } } @@ -74,8 +90,7 @@ TEST(MagFieldEquationTest, neutral) { auto eval = make_equation(dummy_field, ElementaryCharge{0}); { - OdeState s{{0.5, -2, -1}, {0, 0, 1}}; - OdeState result = eval(s); + OdeState result = eval(make_state({0.5, -2, -1}, {0, 0, 1})); EXPECT_VEC_SOFT_EQ(Real3({0, 0, 1}), result.pos); EXPECT_VEC_SOFT_EQ(Real3({0, 0, 0}), result.mom); } diff --git a/test/celeritas/geo/CheckedGeoTrackView.hh b/test/celeritas/geo/CheckedGeoTrackView.hh index b360543199..9d4e7a97ce 100644 --- a/test/celeritas/geo/CheckedGeoTrackView.hh +++ b/test/celeritas/geo/CheckedGeoTrackView.hh @@ -60,7 +60,10 @@ class CheckedGeoTrackView : public GTV } // Check volume consistency this far from the boundary - static constexpr real_type safety_tol() { return 1e-4; } + static constexpr real_type safety_tol() + { + return 1e-4 * units::centimeter; + } // Calculate or return the safety up to an infinite distance inline real_type find_safety(); diff --git a/test/celeritas/geo/GeoMaterial.test.cc b/test/celeritas/geo/GeoMaterial.test.cc index bb25790901..c950bb44cd 100644 --- a/test/celeritas/geo/GeoMaterial.test.cc +++ b/test/celeritas/geo/GeoMaterial.test.cc @@ -6,6 +6,7 @@ //! \file celeritas/geo/GeoMaterial.test.cc //---------------------------------------------------------------------------// #include "corecel/data/CollectionStateStore.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/geo/GeoData.hh" #include "celeritas/geo/GeoMaterialParams.hh" #include "celeritas/geo/GeoMaterialView.hh" @@ -41,7 +42,7 @@ class GeoMaterialTestBase : virtual public GlobalTestBase VecString trace_materials(Real3 const& pos, Real3 dir); }; -auto GeoMaterialTestBase::trace_materials(Real3 const& pos, Real3 dir) +auto GeoMaterialTestBase::trace_materials(Real3 const& pos_cm, Real3 dir) -> VecString { CollectionStateStore host_state{ @@ -55,7 +56,7 @@ auto GeoMaterialTestBase::trace_materials(Real3 const& pos, Real3 dir) // comparison of material IDs encountered. VecString result; - geo = {pos, make_unit_vector(dir)}; + geo = {from_cm(pos_cm), make_unit_vector(dir)}; while (!geo.is_outside()) { result.push_back( diff --git a/test/celeritas/global/AlongStep.test.cc b/test/celeritas/global/AlongStep.test.cc index f4dd8491fe..e54b49591a 100644 --- a/test/celeritas/global/AlongStep.test.cc +++ b/test/celeritas/global/AlongStep.test.cc @@ -9,6 +9,7 @@ #include "celeritas/SimpleCmsTestBase.hh" #include "celeritas/TestEm3Base.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/em/UrbanMscParams.hh" #include "celeritas/ext/GeantPhysicsOptions.hh" #include "celeritas/field/RZMapFieldInput.hh" @@ -518,7 +519,9 @@ TEST_F(SimpleCmsAlongStepTest, msc_field_finegrid) inp.energy = MevEnergy{1.76660104663773580e-3}; // The track is taking its second step in the EM calorimeter, so uses // the cached MSC range values from the previous step - inp.msc_range = {8.43525996595540601e-4, 0.04, 1.34976131122020193e-5}; + inp.msc_range = {from_cm(8.43525996595540601e-4), + 0.04, + from_cm(1.34976131122020193e-5)}; inp.position = { 59.3935490766840459, -109.988210668881749, -81.7228237502843484}; inp.direction = { @@ -539,8 +542,13 @@ TEST_F(SimpleCmsAlongStepTest, msc_field_finegrid) } // Test nearly tangent value nearly on the boundary -TEST_F(SimpleCmsRZFieldAlongStepTest, TEST_IF_CELERITAS_DOUBLE(msc_rzfield)) +TEST_F(SimpleCmsRZFieldAlongStepTest, msc_rzfield) { + if (CELERITAS_REAL_TYPE != CELERITAS_REAL_TYPE_DOUBLE) + { + GTEST_SKIP() << "This edge case only occurs with double"; + } + size_type num_tracks = 128; Input inp; { @@ -569,7 +577,9 @@ TEST_F(SimpleCmsRZFieldAlongStepTest, msc_rzfield_finegrid) inp.energy = MevEnergy{1.76660104663773580e-3}; // The track is taking its second step in the EM calorimeter, so uses // the cached MSC range values from the previous step - inp.msc_range = {8.43525996595540601e-4, 0.04, 1.34976131122020193e-5}; + inp.msc_range = {from_cm(8.43525996595540601e-4), + 0.04, + from_cm(1.34976131122020193e-5)}; inp.position = { 59.3935490766840459, -109.988210668881749, -81.7228237502843484}; inp.direction = { diff --git a/test/celeritas/global/AlongStepTestBase.cc b/test/celeritas/global/AlongStepTestBase.cc index c5c7a08a1b..58529c0cc3 100644 --- a/test/celeritas/global/AlongStepTestBase.cc +++ b/test/celeritas/global/AlongStepTestBase.cc @@ -13,6 +13,7 @@ #include "corecel/io/LogContextException.hh" #include "corecel/io/Repr.hh" #include "corecel/math/ArrayUtils.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/global/ActionRegistry.hh" #include "celeritas/global/CoreParams.hh" #include "celeritas/global/CoreState.hh" @@ -23,6 +24,8 @@ #include "celeritas/track/TrackInitData.hh" #include "celeritas/track/TrackInitParams.hh" +using TimeSecond = celeritas::Quantity; + namespace celeritas { namespace test @@ -44,9 +47,9 @@ auto AlongStepTestBase::run(Input const& inp, size_type num_tracks) -> RunResult Primary p; p.particle_id = inp.particle_id; p.energy = inp.energy; - p.position = inp.position; + p.position = from_cm(inp.position); p.direction = inp.direction; - p.time = inp.time; + p.time = native_value_from(TimeSecond{inp.time}); std::vector primaries(num_tracks, p); for (auto i : range(num_tracks)) @@ -95,10 +98,10 @@ auto AlongStepTestBase::run(Input const& inp, size_type num_tracks) -> RunResult result.eloss += value_as(inp.energy) - value_as(particle.energy()); - result.displacement += distance(geo.pos(), inp.position); + result.displacement += distance(to_cm(geo.pos()), inp.position); result.angle += dot_product(geo.dir(), inp.direction); - result.time += sim.time(); - result.step += sim.step_length(); + result.time += native_value_to(sim.time()).value(); + result.step += to_cm(sim.step_length()); result.mfp += inp.phys_mfp - phys.interaction_mfp(); result.alive += sim.status() == TrackStatus::alive ? 1 : 0; actions[sim.post_step_action()] += 1; diff --git a/test/celeritas/global/AlongStepTestBase.hh b/test/celeritas/global/AlongStepTestBase.hh index ed2c5e8a90..d54fec8f40 100644 --- a/test/celeritas/global/AlongStepTestBase.hh +++ b/test/celeritas/global/AlongStepTestBase.hh @@ -42,9 +42,9 @@ class AlongStepTestBase : virtual public GlobalTestBase { ParticleId particle_id; MevEnergy energy{0}; - Real3 position{0, 0, 0}; + Real3 position{0, 0, 0}; // [cm] Real3 direction{0, 0, 1}; - real_type time{0}; + real_type time{0}; // [s] real_type phys_mfp{1}; //!< Number of MFP to collision MscRange msc_range{0, 0, 0}; diff --git a/test/celeritas/global/KernelContextException.test.cc b/test/celeritas/global/KernelContextException.test.cc index 9fb25a65eb..3722ae0af5 100644 --- a/test/celeritas/global/KernelContextException.test.cc +++ b/test/celeritas/global/KernelContextException.test.cc @@ -15,6 +15,7 @@ #include "corecel/Assert.hh" #include "corecel/Macros.hh" #include "corecel/io/JsonPimpl.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/global/CoreParams.hh" #include "celeritas/global/Stepper.hh" #include "celeritas/phys/PDGNumber.hh" @@ -65,7 +66,7 @@ class KernelContextExceptionTest : public SimpleTestBase, public StepperTestBase Primary p; p.particle_id = this->particle()->find(pdg::gamma()); p.energy = MevEnergy{10}; - p.position = {0, 1, 0}; + p.position = from_cm(Real3{0, 1, 0}); p.direction = {0, 0, 1}; p.time = 0; @@ -137,14 +138,14 @@ TEST_F(KernelContextExceptionTest, typical) EXPECT_EQ(1, e.num_steps()); EXPECT_EQ(ParticleId{0}, e.particle()); EXPECT_EQ(10, e.energy().value()); - EXPECT_VEC_SOFT_EQ((Real3{0, 1, 5}), e.pos()); + EXPECT_VEC_SOFT_EQ(from_cm(Real3{0, 1, 5}), e.pos()); EXPECT_VEC_SOFT_EQ((Real3{0, 0, 1}), e.dir()); if (CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_ORANGE) { EXPECT_EQ(VolumeId{2}, e.volume()); EXPECT_EQ(SurfaceId{11}, e.surface()); } - if (CELERITAS_USE_JSON + if (CELERITAS_USE_JSON && CELERITAS_UNITS == CELERITAS_UNITS_CGS && CELERITAS_CORE_GEO == CELERITAS_CORE_GEO_ORANGE) { std::stringstream ss; diff --git a/test/celeritas/global/Stepper.test.cc b/test/celeritas/global/Stepper.test.cc index a0b4f8a1e7..4903f9106c 100644 --- a/test/celeritas/global/Stepper.test.cc +++ b/test/celeritas/global/Stepper.test.cc @@ -12,6 +12,7 @@ #include "corecel/Types.hh" #include "corecel/cont/Range.hh" #include "corecel/cont/Span.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/em/UrbanMscParams.hh" #include "celeritas/ext/GeantPhysicsOptions.hh" #include "celeritas/field/UniformFieldData.hh" @@ -50,7 +51,7 @@ class SimpleComptonTest : public SimpleTestBase, public StepperTestBase CELER_ASSERT(p.particle_id); p.energy = units::MevEnergy{100}; p.track_id = TrackId{0}; - p.position = {-22, 0, 0}; + p.position = from_cm(Real3{-22, 0, 0}); p.direction = {1, 0, 0}; p.time = 0; @@ -80,7 +81,7 @@ class TestEm3StepperTestBase : public TestEm3Base, public StepperTestBase CELER_ASSERT(p.particle_id); p.energy = energy; p.track_id = TrackId{0}; - p.position = {-22, 0, 0}; + p.position = from_cm(Real3{-22, 0, 0}); p.direction = {1, 0, 0}; p.time = 0; diff --git a/test/celeritas/global/StepperTestBase.cc b/test/celeritas/global/StepperTestBase.cc index e6def736f7..ce6b28c882 100644 --- a/test/celeritas/global/StepperTestBase.cc +++ b/test/celeritas/global/StepperTestBase.cc @@ -43,6 +43,7 @@ StepperTestBase::StepperTestBase() bool StepperTestBase::is_default_build() { return CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE + && CELERITAS_UNITS == CELERITAS_UNITS_CGS && CELERITAS_CORE_RNG == CELERITAS_CORE_RNG_XORWOW; } diff --git a/test/celeritas/io/EventIOTestBase.cc b/test/celeritas/io/EventIOTestBase.cc index 52b896d476..9ea7e008e9 100644 --- a/test/celeritas/io/EventIOTestBase.cc +++ b/test/celeritas/io/EventIOTestBase.cc @@ -12,6 +12,7 @@ #include "corecel/io/Logger.hh" #include "corecel/io/Repr.hh" #include "celeritas/Quantities.hh" +#include "celeritas/UnitUtils.hh" #include "TestMacros.hh" @@ -108,11 +109,11 @@ auto EventIOTestBase::read_all(Reader& read_event) const -> ReadAllResult result.pdg.push_back( particles_->id_to_pdg(p.particle_id).unchecked_get()); result.energy.push_back(p.energy.value()); - result.pos.insert( - result.pos.end(), p.position.begin(), p.position.end()); + auto pos = to_cm(p.position); + result.pos.insert(result.pos.end(), pos.begin(), pos.end()); result.dir.insert( result.dir.end(), p.direction.begin(), p.direction.end()); - result.time.push_back(p.time); + result.time.push_back(p.time / units::second); result.event.push_back(p.event_id.unchecked_get()); result.track.push_back(p.track_id.unchecked_get()); } @@ -131,20 +132,20 @@ void EventIOTestBase::write_test_event(Writer& write_event) const auto gamma_id = particles.find(pdg::gamma()); Primary gamma{gamma_id, MevEnergy{1.23}, - Real3{2, 4, 5}, + from_cm(Real3{2, 4, 5}), Real3{1, 0, 0}, 5.67e-9 * units::second, EventId{0}, TrackId{}}; Primary proton{proton_id, MevEnergy{2.34}, - Real3{3, 5, 8}, + from_cm(Real3{3, 5, 8}), Real3{0, 1, 0}, 5.78e-9 * units::second, EventId{0}, TrackId{}}; std::vector primaries{gamma, proton, gamma, proton}; - primaries[1].position = {-3, -4, 5}; + primaries[1].position = from_cm(Real3{-3, -4, 5}); primaries[3].position = primaries[2].position; primaries[3].time = primaries[2].time; for (auto i : range(primaries.size())) diff --git a/test/celeritas/io/ImportUnits.test.cc b/test/celeritas/io/ImportUnits.test.cc new file mode 100644 index 0000000000..378e99e92f --- /dev/null +++ b/test/celeritas/io/ImportUnits.test.cc @@ -0,0 +1,91 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/io/ImportUnits.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/io/ImportUnits.hh" + +#include "corecel/math/Algorithms.hh" +#include "celeritas/Units.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// + +class ImportUnitsTest : public ::celeritas::test::Test +{ +}; + +TEST_F(ImportUnitsTest, native) +{ + auto from_native = [](ImportUnits iu) { + return native_value_from(UnitSystem::native, iu); + }; + + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::unitless)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::mev)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::mev_per_len)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::len)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::len_inv)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::len_mev_inv)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::mev_sq_per_len)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::len_sq)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::time)); + EXPECT_SOFT_EQ(1.0, from_native(ImportUnits::inv_len_cb)); +} + +TEST_F(ImportUnitsTest, csg) +{ + constexpr auto cm = units::centimeter; + constexpr auto cm_sq = ipow<2>(units::centimeter); + constexpr auto cm_cb = ipow<3>(units::centimeter); + constexpr auto s = units::second; + + auto from_cgs = [](ImportUnits iu) { + return native_value_from(UnitSystem::cgs, iu); + }; + + EXPECT_SOFT_EQ(1.0, from_cgs(ImportUnits::unitless)); + EXPECT_SOFT_EQ(1.0, from_cgs(ImportUnits::mev)); + EXPECT_SOFT_EQ(1.0 / cm, from_cgs(ImportUnits::mev_per_len)); + EXPECT_SOFT_EQ(cm, from_cgs(ImportUnits::len)); + EXPECT_SOFT_EQ(1 / cm, from_cgs(ImportUnits::len_inv)); + EXPECT_SOFT_EQ(1 / cm, from_cgs(ImportUnits::len_mev_inv)); + EXPECT_SOFT_EQ(1 / cm, from_cgs(ImportUnits::mev_sq_per_len)); + EXPECT_SOFT_EQ(cm_sq, from_cgs(ImportUnits::len_sq)); + EXPECT_SOFT_EQ(s, from_cgs(ImportUnits::time)); + EXPECT_SOFT_EQ(1 / cm_cb, from_cgs(ImportUnits::inv_len_cb)); +} + +TEST_F(ImportUnitsTest, clhep) +{ + constexpr auto mm = units::millimeter; + constexpr auto mm_sq = ipow<2>(units::millimeter); + constexpr auto mm_cb = ipow<3>(units::millimeter); + constexpr auto ns = units::nanosecond; + + auto from_clhep = [](ImportUnits iu) { + return native_value_from(UnitSystem::clhep, iu); + }; + + EXPECT_SOFT_EQ(1.0, from_clhep(ImportUnits::unitless)); + EXPECT_SOFT_EQ(1.0, from_clhep(ImportUnits::mev)); + EXPECT_SOFT_EQ(1.0 / mm, from_clhep(ImportUnits::mev_per_len)); + EXPECT_SOFT_EQ(mm, from_clhep(ImportUnits::len)); + EXPECT_SOFT_EQ(1 / mm, from_clhep(ImportUnits::len_inv)); + EXPECT_SOFT_EQ(1 / mm, from_clhep(ImportUnits::len_mev_inv)); + EXPECT_SOFT_EQ(1 / mm, from_clhep(ImportUnits::mev_sq_per_len)); + EXPECT_SOFT_EQ(mm_sq, from_clhep(ImportUnits::len_sq)); + EXPECT_SOFT_EQ(ns, from_clhep(ImportUnits::time)); + EXPECT_SOFT_EQ(1 / mm_cb, from_clhep(ImportUnits::inv_len_cb)); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/mat/ElementSelector.test.cc b/test/celeritas/mat/ElementSelector.test.cc index 88c80727b6..7d4ddc81d3 100644 --- a/test/celeritas/mat/ElementSelector.test.cc +++ b/test/celeritas/mat/ElementSelector.test.cc @@ -38,7 +38,7 @@ class ElementSelectorTest : public Test protected: void SetUp() override { - using units::AmuMass; + using namespace units; MaterialParams::Input inp; inp.elements = { @@ -49,12 +49,12 @@ class ElementSelectorTest : public Test }; inp.materials = { {0.0, 0.0, MatterState::unspecified, {}, "hard_vacuum"}, - {0.1 * constants::na_avogadro, + {native_value_from(MolCcDensity{0.1}), 293.0, MatterState::gas, {{ElementId{2}, 1.0}}, "Al"}, - {0.05 * constants::na_avogadro, + {native_value_from(MolCcDensity{0.05}), 293.0, MatterState::solid, {{ElementId{0}, 0.25}, @@ -62,7 +62,7 @@ class ElementSelectorTest : public Test {ElementId{2}, 0.25}, {ElementId{3}, 0.25}}, "everything_even"}, - {1 * constants::na_avogadro, + {native_value_from(MolCcDensity{1}), 293.0, MatterState::solid, {{ElementId{0}, 0.48}, diff --git a/test/celeritas/mat/Material.test.cc b/test/celeritas/mat/Material.test.cc index 1b93980f71..0abfb1c2d4 100644 --- a/test/celeritas/mat/Material.test.cc +++ b/test/celeritas/mat/Material.test.cc @@ -12,6 +12,7 @@ #include "corecel/data/CollectionStateStore.hh" #include "celeritas/Quantities.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/ext/RootImporter.hh" #include "celeritas/ext/ScopedRootErrorHandler.hh" #include "celeritas/io/ImportData.hh" @@ -27,6 +28,11 @@ #include "MaterialTestBase.hh" #include "celeritas_test.hh" +using celeritas::units::GramCcDensity; +using celeritas::units::InvCcDensity; +using celeritas::units::MevEnergy; +using celeritas::units::MolCcDensity; + namespace celeritas { //---------------------------------------------------------------------------// @@ -79,7 +85,8 @@ TEST(MaterialUtils, radiation_length) el.atomic_mass = units::AmuMass{amu_mass}; el.coulomb_correction = calc_coulomb_correction(el.atomic_number); - return 1 / detail::calc_mass_rad_coeff(el); + return 1 / detail::calc_mass_rad_coeff(el) + * (ipow<2>(units::centimeter) / units::gram); }; // Hydrogen @@ -157,6 +164,11 @@ TEST_F(MaterialTest, params) TEST_F(MaterialTest, material_view) { + if (CELERITAS_UNITS != CELERITAS_UNITS_CGS) + { + GTEST_SKIP() << "Test needs updating to include units"; + } + { // NaI MaterialView mat = params->get(MaterialId{0}); @@ -220,14 +232,20 @@ TEST_F(MaterialTest, material_view) { // H2_3 MaterialView mat = params->get(MaterialId{3}); - EXPECT_SOFT_EQ(1.072e+20, mat.number_density()); + EXPECT_SOFT_EQ( + 1.072e+20, + native_value_to(mat.number_density()).value()); EXPECT_SOFT_EQ(110, mat.temperature()); EXPECT_EQ(MatterState::gas, mat.matter_state()); EXPECT_SOFT_EQ(1.0, mat.zeff()); - EXPECT_SOFT_EQ(0.00017943386624303615, mat.density()); - EXPECT_SOFT_EQ(1.072e+20, mat.electron_density()); - EXPECT_SOFT_EQ(351367.47504673258, mat.radiation_length()); - EXPECT_SOFT_EQ(19.2e-6, mat.mean_excitation_energy().value()); + EXPECT_SOFT_EQ(0.00017943386624303615, + native_value_to(mat.density()).value()); + EXPECT_SOFT_EQ( + 1.072e+20, + native_value_to(mat.electron_density()).value()); + EXPECT_SOFT_EQ(351367.47504673258, to_cm(mat.radiation_length())); + EXPECT_SOFT_EQ(19.2e-6, + value_as(mat.mean_excitation_energy())); EXPECT_SOFT_EQ(std::log(19.2e-6), mat.log_mean_excitation_energy().value()); @@ -248,7 +266,9 @@ TEST_F(MaterialTest, element_view) EXPECT_SOFT_EQ(std::pow(13.0 * 14.0, 1.0 / 3), el.cbrt_zzp()); EXPECT_SOFT_EQ(std::log(13.0), el.log_z()); EXPECT_SOFT_EQ(0.010734632775699565, el.coulomb_correction()); - EXPECT_SOFT_EQ(0.04164723292591279, el.mass_radiation_coeff()); + EXPECT_SOFT_EQ( + 0.04164723292591279 * ipow<2>(units::centimeter) / units::gram, + el.mass_radiation_coeff()); // Test its isotopes EXPECT_EQ(2, el.num_isotopes()); @@ -292,7 +312,7 @@ TEST_F(MaterialTest, TEST_IF_CELERITAS_DOUBLE(output)) MaterialParamsOutput out(params); EXPECT_EQ("material", out.label()); - if (CELERITAS_USE_JSON) + if (CELERITAS_USE_JSON && CELERITAS_UNITS == CELERITAS_UNITS_CGS) { EXPECT_JSON_EQ( R"json({"_units":{"atomic_mass":"amu","mean_excitation_energy":"MeV","nuclear_mass":"MeV/c^2"},"elements":{"atomic_mass":[1.008,26.9815385,22.98976928,126.90447],"atomic_number":[1,13,11,53],"coulomb_correction":[6.400821803338426e-05,0.010734632775699565,0.00770256745342534,0.15954439947436763],"isotope_fractions":[[0.9,0.1],[0.7,0.3],[1.0],[0.05,0.15,0.8]],"isotope_ids":[[0,1],[2,3],[4],[5,6,7]],"label":["H","Al","Na","I"],"mass_radiation_coeff":[0.0158611264432063,0.04164723292591279,0.03605392839455309,0.11791841505608874]},"isotopes":{"atomic_mass_number":[1,2,27,28,23,125,126,127],"atomic_number":[1,1,13,13,11,53,53,53],"label":["1H","2H","27Al","28Al","23Na","125I","126I","127I"],"nuclear_mass":[938.272,1875.61,25126.5,26058.3,21409.2,116321.0,117253.0,118184.0]},"materials":{"density":[3.6700020622594716,0.0,0.00017976000000000003,0.00017943386624303615],"electron_density":[9.4365282069664e+23,0.0,1.073948435904467e+20,1.072e+20],"element_frac":[[0.5,0.5],[],[1.0],[1.0]],"element_id":[[2,3],[],[0],[0]],"label":["NaI","hard vacuum","H2@1","H2@2"],"matter_state":["solid","unspecified","gas","gas"],"mean_excitation_energy":[0.00040000760709482647,0.0,1.9199999999999986e-05,1.9199999999999986e-05],"number_density":[2.948915064677e+22,0.0,1.073948435904467e+20,1.072e+20],"radiation_length":[3.5393292693170424,null,350729.99844063615,351367.4750467326],"temperature":[293.0,0.0,100.0,110.0],"zeff":[32.0,0.0,1.0,1.0]}})json", @@ -339,10 +359,13 @@ TEST_F(MaterialParamsImportTest, TEST_IF_CELERITAS_USE_ROOT(import_materials)) EXPECT_EQ(MatterState::solid, mat.matter_state()); EXPECT_SOFT_EQ(293.15, mat.temperature()); // [K] - EXPECT_SOFT_EQ(7.9999999972353661, mat.density()); // [g/cm^3] - EXPECT_SOFT_EQ(2.2444320228819809e+24, - mat.electron_density()); // [1/cm^3] - EXPECT_SOFT_EQ(8.6993489258991514e+22, mat.number_density()); // [1/cm^3] + EXPECT_SOFT_EQ(7.9999999972353661, + native_value_to(mat.density()).value()); + EXPECT_SOFT_EQ( + 2.2444320228819809e+24, + native_value_to(mat.electron_density()).value()); + EXPECT_SOFT_EQ(8.6993489258991514e+22, + native_value_to(mat.number_density()).value()); // Test elements by unpacking them std::vector els; diff --git a/test/celeritas/mat/Material.test.cu b/test/celeritas/mat/Material.test.cu index daa7154546..9dc4c1fd3b 100644 --- a/test/celeritas/mat/Material.test.cu +++ b/test/celeritas/mat/Material.test.cu @@ -11,8 +11,10 @@ #include "corecel/device_runtime_api.h" #include "corecel/cont/Range.hh" +#include "corecel/math/Quantity.hh" #include "corecel/sys/Device.hh" #include "corecel/sys/KernelParamCalculator.device.hh" +#include "celeritas/Quantities.hh" #include "celeritas/mat/MaterialTrackView.hh" using thrust::raw_pointer_cast; @@ -48,7 +50,8 @@ __global__ void m_test_kernel(unsigned int const size, // Get material properties auto const& mat = mat_track.make_material_view(); temperatures[tid.get()] = mat.temperature(); - rad_len[tid.get()] = mat.radiation_length(); + rad_len[tid.get()] + = native_value_to(mat.radiation_length()).value(); // Fill elements with finctional cross sections Span scratch = mat_track.element_scratch(); @@ -65,7 +68,10 @@ __global__ void m_test_kernel(unsigned int const size, for (auto ec : range(mat.num_elements())) { // Get its atomic number weighted by its fractional number density - tz += scratch[ec] * mat.get_element_density(ElementComponentId{ec}); + tz += scratch[ec] + * native_value_to( + mat.get_element_density(ElementComponentId{ec})) + .value(); } tot_z[tid.get()] = tz; } diff --git a/test/celeritas/mat/MaterialTestBase.hh b/test/celeritas/mat/MaterialTestBase.hh index b27740df6c..91d6ebb11f 100644 --- a/test/celeritas/mat/MaterialTestBase.hh +++ b/test/celeritas/mat/MaterialTestBase.hh @@ -22,8 +22,7 @@ class MaterialTestBase SPConstMaterial build_material() { - using units::AmuMass; - using units::MevMass; + using namespace celeritas::units; MaterialParams::Input inp; @@ -64,7 +63,7 @@ class MaterialTestBase inp.materials = { // Sodium iodide - {2.948915064677e+22, + {native_value_from(InvCcDensity{2.948915064677e+22}), 293.0, MatterState::solid, {{ElementId{2}, 0.5}, {ElementId{3}, 0.5}}, @@ -72,13 +71,13 @@ class MaterialTestBase // Void {0, 0, MatterState::unspecified, {}, "hard vacuum"}, // Diatomic hydrogen - {1.0739484359044669e+20, + {native_value_from(InvCcDensity{1.0739484359044669e+20}), 100.0, MatterState::gas, {{ElementId{0}, 1.0}}, Label{"H2", "1"}}, // Diatomic hydrogen with the same name and different properties - {1.072e+20, + {native_value_from(InvCcDensity{1.072e+20}), 110.0, MatterState::gas, {{ElementId{0}, 1.0}}, diff --git a/test/celeritas/optical/Cerenkov.test.cc b/test/celeritas/optical/Cerenkov.test.cc index 5a85b6fbc3..68f9c1c7e1 100644 --- a/test/celeritas/optical/Cerenkov.test.cc +++ b/test/celeritas/optical/Cerenkov.test.cc @@ -13,6 +13,7 @@ #include "corecel/math/Algorithms.hh" #include "corecel/math/ArrayOperators.hh" #include "corecel/math/ArrayUtils.hh" +#include "corecel/math/Quantity.hh" #include "celeritas/Constants.hh" #include "celeritas/Units.hh" #include "celeritas/grid/GenericGridData.hh" @@ -29,6 +30,16 @@ namespace celeritas { namespace test { +struct InvCentimeter +{ + static CELER_CONSTEXPR_FUNCTION real_type value() + { + return 1 / units::centimeter; + } + static char const* label() { return "1/cm"; } +}; + +using InvCmDnDx = Quantity; //---------------------------------------------------------------------------// /*! * Tabulated refractive index in water as a function of photon wavelength [μm]. @@ -173,6 +184,7 @@ TEST_F(CerenkovTest, dndx) { EXPECT_SOFT_NEAR(369.81e6, constants::alpha_fine_structure * units::Mev::value() + * units::centimeter / (constants::hbar_planck * constants::c_light), 1e-6); @@ -182,7 +194,9 @@ TEST_F(CerenkovTest, dndx) for (real_type beta : {0.5, 0.6813, 0.69, 0.71, 0.73, 0.752, 0.756, 0.8, 0.9, 0.999}) { - dndx.push_back(calc_dndx(units::LightSpeed(beta))); + dndx.push_back( + native_value_to(calc_dndx(units::LightSpeed(beta))) + .value()); } if (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) { @@ -297,7 +311,7 @@ TEST_F(CerenkovTest, TEST_IF_CELERITAS_DOUBLE(generator)) } avg_costheta /= total_num_photons; avg_energy /= total_num_photons; - avg_displacement /= total_num_photons; + avg_displacement /= (units::centimeter * total_num_photons); avg_engine_samples = real_type(rng.count()) / num_samples; }; diff --git a/test/celeritas/optical/ScintillationGenerator.test.cc b/test/celeritas/optical/ScintillationGenerator.test.cc index 9eb6f1a8b5..026f66822e 100644 --- a/test/celeritas/optical/ScintillationGenerator.test.cc +++ b/test/celeritas/optical/ScintillationGenerator.test.cc @@ -128,7 +128,7 @@ TEST_F(ScintillationGeneratorTest, basic) for (auto i : range(dist_.num_photons)) { energy.push_back(photons[i].energy.value()); - time.push_back(photons[i].time); + time.push_back(photons[i].time / units::second); cos_theta.push_back( dot_product(photons[i].direction, dist_.points[StepPoint::post].pos diff --git a/test/celeritas/phys/CutoffParams.test.cc b/test/celeritas/phys/CutoffParams.test.cc index ddfc514b76..548f10c903 100644 --- a/test/celeritas/phys/CutoffParams.test.cc +++ b/test/celeritas/phys/CutoffParams.test.cc @@ -11,6 +11,7 @@ #include "celeritas/Quantities.hh" #include "celeritas/RootTestBase.hh" #include "celeritas/Types.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/mat/ElementView.hh" #include "celeritas/mat/MaterialData.hh" #include "celeritas/mat/MaterialParams.hh" @@ -48,7 +49,7 @@ class CutoffParamsTest : public Test }; m_input.materials = { // Sodium iodide - {2.948915064677e+22, + {native_value_from(InvCcDensity{2.948915064677e+22}), 293.0, MatterState::solid, {{ElementId{2}, 0.5}, {ElementId{3}, 0.5}}, @@ -56,7 +57,7 @@ class CutoffParamsTest : public Test // Void {0, 0, MatterState::unspecified, {}, "hard vacuum"}, // Diatomic hydrogen - {1.0739484359044669e+20, + {native_value_from(InvCcDensity{1.0739484359044669e+20}), 100.0, MatterState::gas, {{ElementId{0}, 1.0}}, @@ -225,7 +226,7 @@ TEST_F(CutoffParamsImportTest, import_cutoffs) { CutoffView cutoffs(this->cutoff()->host_ref(), mid); energies.push_back(cutoffs.energy(pid).value()); - ranges.push_back(cutoffs.range(pid)); + ranges.push_back(to_cm(cutoffs.range(pid))); EXPECT_FALSE(cutoffs.apply_post_interaction()); } } diff --git a/test/celeritas/phys/InteractorHostTestBase.cc b/test/celeritas/phys/InteractorHostTestBase.cc index f01bb53453..9aa8e2fcb9 100644 --- a/test/celeritas/phys/InteractorHostTestBase.cc +++ b/test/celeritas/phys/InteractorHostTestBase.cc @@ -67,27 +67,27 @@ InteractorHostTestBase::InteractorHostTestBase() {AtomicNumber{74}, AmuMass{183.84}, {}, Label{"W"}}, {AtomicNumber{82}, AmuMass{207.2}, {}, Label{"Pb"}}}; mat_inp.materials = { - {0.141 * na_avogadro, + {native_value_from(MolCcDensity{0.141}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, Label{"Cu"}}, - {0.05477 * na_avogadro, + {native_value_from(MolCcDensity{0.05477}), 293.15, MatterState::solid, {{ElementId{0}, 1.0}}, Label{"Pb"}}, - {1e-5 * na_avogadro, + {native_value_from(MolCcDensity{1e-5}), 293., MatterState::solid, {{ElementId{1}, 1.0}}, Label{"K"}}, - {1.0 * na_avogadro, + {native_value_from(MolCcDensity{1.0}), 293.0, MatterState::solid, {{ElementId{0}, 1.0}}, Label{"Cu-1.0"}}, - {1.0 * constants::na_avogadro, + {native_value_from(MolCcDensity{1.0}), 293.0, MatterState::solid, {{ElementId{2}, 0.5}, {ElementId{3}, 0.3}, {ElementId{4}, 0.2}}, diff --git a/test/celeritas/phys/MockProcess.cc b/test/celeritas/phys/MockProcess.cc index c7a9698076..0c6052c543 100644 --- a/test/celeritas/phys/MockProcess.cc +++ b/test/celeritas/phys/MockProcess.cc @@ -30,7 +30,7 @@ MockProcess::MockProcess(Input data) : data_(std::move(data)) || std::any_of(data_.xs.begin(), data_.xs.end(), [](BarnMicroXs x) { return x > zero_quantity(); })); - CELER_EXPECT(data_.energy_loss >= 0); + CELER_EXPECT(data_.energy_loss >= zero_quantity()); } //---------------------------------------------------------------------------// @@ -52,14 +52,14 @@ auto MockProcess::build_models(ActionIdIter start_id) const -> VecModel } //---------------------------------------------------------------------------// -auto MockProcess::step_limits(Applicability range) const -> StepLimitBuilders +auto MockProcess::step_limits(Applicability applic) const -> StepLimitBuilders { - CELER_EXPECT(range.material); - CELER_EXPECT(range.particle); + CELER_EXPECT(applic.material); + CELER_EXPECT(applic.particle); using VecDbl = std::vector; - MaterialView mat(data_.materials->host_ref(), range.material); + MaterialView mat(data_.materials->host_ref(), applic.material); real_type numdens = mat.number_density(); StepLimitBuilders builders; @@ -67,25 +67,28 @@ auto MockProcess::step_limits(Applicability range) const -> StepLimitBuilders { VecDbl xs_grid; for (auto xs : data_.xs) + { xs_grid.push_back(native_value_from(xs) * numdens); + } builders[ValueGridType::macro_xs] = std::make_unique( - range.lower.value(), range.upper.value(), xs_grid); + applic.lower.value(), applic.upper.value(), xs_grid); } - if (data_.energy_loss > 0) + if (data_.energy_loss > zero_quantity()) { - real_type eloss_rate = data_.energy_loss * numdens; + auto eloss_rate = native_value_to( + native_value_from(data_.energy_loss) * numdens); builders[ValueGridType::energy_loss] = std::make_unique( - range.lower.value(), - range.upper.value(), - VecDbl{eloss_rate, eloss_rate}); + applic.lower.value(), + applic.upper.value(), + VecDbl{eloss_rate.value(), eloss_rate.value()}); builders[ValueGridType::range] = std::make_unique( - range.lower.value(), - range.upper.value(), - VecDbl{range.lower.value() / eloss_rate, - range.upper.value() / eloss_rate}); + applic.lower.value(), + applic.upper.value(), + VecDbl{applic.lower.value() / eloss_rate.value(), + applic.upper.value() / eloss_rate.value()}); } return builders; diff --git a/test/celeritas/phys/MockProcess.hh b/test/celeritas/phys/MockProcess.hh index b4bdab38c4..91b706d044 100644 --- a/test/celeritas/phys/MockProcess.hh +++ b/test/celeritas/phys/MockProcess.hh @@ -19,6 +19,29 @@ namespace celeritas { namespace test { +//---------------------------------------------------------------------------// +//! Energy loss rate [MeV/cm] per volume [cm^-3] -> [MeV * cm^2] +struct MevCmSq +{ + static CELER_CONSTEXPR_FUNCTION real_type value() + { + return units::Mev::value() * ipow<2>(units::centimeter); + } +}; + +using MevCmSqLossDens = Quantity; + +//! Energy loss rate +struct MevPerCm +{ + static CELER_CONSTEXPR_FUNCTION real_type value() + { + return units::Mev::value() / units::centimeter; + } +}; + +using MevPerCmLoss = Quantity; + //---------------------------------------------------------------------------// /*! * Mock process. @@ -55,7 +78,7 @@ class MockProcess : public Process VecApplicability applic; //!< Applicablity per model ModelCallback interact; //!< MockModel::interact callback VecMicroXs xs; //!< Constant per atom [bn] - real_type energy_loss{}; //!< Constant per atom [MeV/cm / cm^-3] + MevCmSqLossDens energy_loss{}; //!< Constant per atom }; public: diff --git a/test/celeritas/phys/Particle.test.cc b/test/celeritas/phys/Particle.test.cc index f279d43b07..7630db08f6 100644 --- a/test/celeritas/phys/Particle.test.cc +++ b/test/celeritas/phys/Particle.test.cc @@ -91,7 +91,7 @@ TEST_F(ParticleTest, TEST_IF_CELERITAS_DOUBLE(output)) ParticleParamsOutput out(this->particle_params); EXPECT_EQ("particle", out.label()); - if (CELERITAS_USE_JSON) + if (CELERITAS_USE_JSON && CELERITAS_UNITS == CELERITAS_UNITS_CGS) { EXPECT_EQ( R"json({"_units":{"charge":"e","mass":"MeV/c^2"},"charge":[-1.0,0.0,0.0,1.0],"decay_constant":[0.0,0.0,0.0011371389583807142,0.0],"is_antiparticle":[false,false,false,true],"label":["electron","gamma","neutron","positron"],"mass":[0.5109989461,0.0,939.565413,0.5109989461],"pdg":[11,22,2112,-11]})json", @@ -184,8 +184,8 @@ TEST_F(ParticleTestHost, electron) EXPECT_FALSE(particle.is_antiparticle()); EXPECT_TRUE(particle.is_stable()); EXPECT_SOFT_EQ(0.74453076757415848, particle.beta_sq()); - EXPECT_SOFT_EQ(0.86286196322132447, particle.speed().value()); - EXPECT_SOFT_EQ(25867950886.882648, native_value_from(particle.speed())); + EXPECT_SOFT_EQ(0.86286196322132447, + value_as(particle.speed())); EXPECT_SOFT_EQ(1.9784755992474248, particle.lorentz_factor()); EXPECT_SOFT_EQ(0.87235253544653601, particle.momentum().value()); EXPECT_SOFT_EQ(0.7609989461, particle.momentum_sq().value()); @@ -235,7 +235,8 @@ TEST_F(ParticleTestHost, neutron) particle = Initializer_t{ParticleId{2}, MevEnergy{20}}; EXPECT_REAL_EQ(20, particle.energy().value()); - EXPECT_REAL_EQ(1.0 / 879.4, particle.decay_constant()); + EXPECT_REAL_EQ(1.0 / 879.4 * (1 / units::second), + particle.decay_constant()); EXPECT_FALSE(particle.is_antiparticle()); EXPECT_FALSE(particle.is_stable()); } diff --git a/test/celeritas/phys/Particle.test.cu b/test/celeritas/phys/Particle.test.cu index e1a1276800..d99bef8548 100644 --- a/test/celeritas/phys/Particle.test.cu +++ b/test/celeritas/phys/Particle.test.cu @@ -10,8 +10,11 @@ #include #include "corecel/device_runtime_api.h" +#include "corecel/math/Quantity.hh" +#include "corecel/math/UnitUtils.hh" #include "corecel/sys/Device.hh" #include "corecel/sys/KernelParamCalculator.device.hh" +#include "celeritas/UnitTypes.hh" #include "celeritas/phys/ParticleTrackView.hh" using thrust::raw_pointer_cast; @@ -32,6 +35,8 @@ __global__ void ptv_test_kernel(unsigned int size, ParticleTrackInitializer const* init, double* result) { + using InvSecDecay = Quantity>; + auto local_tid = TrackSlotId{KernelParamCalculator::thread_id().unchecked_get()}; if (!(local_tid < size)) @@ -49,7 +54,7 @@ __global__ void ptv_test_kernel(unsigned int size, *result++ = p.energy().value(); *result++ = p.mass().value(); *result++ = p.charge().value(); - *result++ = p.decay_constant(); + *result++ = native_value_to(p.decay_constant()).value(); *result++ = p.speed().value(); *result++ = (p.mass() > zero_quantity() ? p.lorentz_factor() : -1); *result++ = p.momentum().value(); diff --git a/test/celeritas/phys/Physics.test.cc b/test/celeritas/phys/Physics.test.cc index 1829e41249..8148135948 100644 --- a/test/celeritas/phys/Physics.test.cc +++ b/test/celeritas/phys/Physics.test.cc @@ -12,6 +12,7 @@ #include "corecel/cont/Range.hh" #include "corecel/data/CollectionStateStore.hh" #include "celeritas/MockTestBase.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/em/process/EPlusAnnihilationProcess.hh" #include "celeritas/grid/EnergyLossCalculator.hh" #include "celeritas/grid/RangeCalculator.hh" @@ -32,6 +33,15 @@ namespace test //---------------------------------------------------------------------------// using MevEnergy = units::MevEnergy; +namespace +{ +real_type to_inv_cm(real_type xs_native) +{ + return native_value_to(xs_native).value(); +} + +} // namespace + //---------------------------------------------------------------------------// // PHYSICS BASE CLASS //---------------------------------------------------------------------------// @@ -130,6 +140,10 @@ TEST_F(PhysicsParamsTest, output) { GTEST_SKIP() << "Test results are based on double-precision data"; } + if (CELERITAS_UNITS != CELERITAS_UNITS_CGS) + { + GTEST_SKIP() << "Test results are based on CGS units"; + } if (!CELERITAS_USE_JSON) { GTEST_SKIP() << "JSON required to test output"; @@ -290,7 +304,7 @@ TEST_F(PhysicsTrackViewHostTest, track_view) auto s = celer.range_to_step(r); EXPECT_GT(s, real_type(0)); EXPECT_LE(s, r) << "s - r == " << s - r; - step.push_back(s); + step.push_back(to_cm(s)); } if (CELERITAS_REAL_TYPE != CELERITAS_REAL_TYPE_DOUBLE) @@ -451,7 +465,7 @@ TEST_F(PhysicsTrackViewHostTest, calc_xs) auto id = phys.value_grid(ValueGridType::macro_xs, scat_ppid); ASSERT_TRUE(id); auto calc_xs = phys.make_calculator(id); - xs.push_back(calc_xs(MevEnergy{1.0})); + xs.push_back(to_inv_cm(calc_xs(MevEnergy{1.0}))); } } @@ -487,9 +501,11 @@ TEST_F(PhysicsTrackViewHostTest, calc_eloss_range) auto calc_range = phys.make_calculator(range_id); for (real_type energy : {1e-6, 0.01, 1.0, 1e2}) { - eloss.push_back(calc_eloss(MevEnergy{energy})); - range.push_back(calc_range(MevEnergy{energy})); - step.push_back(phys.range_to_step(range.back())); + // Energy loss rate per unit length (MeV / len) + eloss.push_back(calc_eloss(MevEnergy{energy}) * units::centimeter); + auto r = calc_range(MevEnergy{energy}); + range.push_back(to_cm(r)); + step.push_back(to_cm(phys.range_to_step(r))); } } @@ -526,7 +542,8 @@ TEST_F(PhysicsTrackViewHostTest, use_integral) EXPECT_FALSE(phys.integral_xs_process(ppid)); MaterialView material = this->material()->get(MaterialId{2}); - EXPECT_SOFT_EQ(0.1, phys.calc_xs(ppid, material, MevEnergy{1.0})); + EXPECT_SOFT_EQ( + 0.1, to_inv_cm(phys.calc_xs(ppid, material, MevEnergy{1.0}))); } { // Energy loss tables and energy-dependent macro xs @@ -540,9 +557,10 @@ TEST_F(PhysicsTrackViewHostTest, use_integral) MaterialView material = this->material()->get(MaterialId{2}); for (real_type energy : {0.001, 0.01, 0.1, 0.11, 10.0}) { - xs.push_back(phys.calc_xs(ppid, material, MevEnergy{energy})); - max_xs.push_back(phys.calc_max_xs( - integral_proc, ppid, material, MevEnergy{energy})); + xs.push_back( + to_inv_cm(phys.calc_xs(ppid, material, MevEnergy{energy}))); + max_xs.push_back(to_inv_cm(phys.calc_max_xs( + integral_proc, ppid, material, MevEnergy{energy}))); } double const expected_xs[] = {0.6, 36. / 55, 1.2, 1979. / 1650, 0.6}; double const expected_max_xs[] = {0.6, 36. / 55, 1.2, 1.2, 357. / 495}; @@ -620,7 +638,8 @@ TEST_F(PhysicsTrackViewHostTest, cuda_surrogate) for (real_type energy : {1e-5, 1e-3, 1., 100., 1e5}) { - step.push_back(test::calc_step(phys, pstep, MevEnergy{energy})); + step.push_back( + to_cm(test::calc_step(phys, pstep, MevEnergy{energy}))); } } @@ -733,7 +752,7 @@ auto EPlusAnnihilationTest::build_material() -> SPConstMaterial MaterialParams::Input mi; mi.elements = {{AtomicNumber{19}, AmuMass{39.0983}, {}, "K"}}; - mi.materials = {{1e-5 * na_avogadro, + mi.materials = {{native_value_from(MolCcDensity{1e-5}), 293., MatterState::solid, {{ElementId{0}, 1.0}}, @@ -806,8 +825,9 @@ TEST_F(EPlusAnnihilationTest, host_track_view) // Check cross section MaterialView material_view = this->material()->get(MaterialId{0}); - EXPECT_SOFT_EQ(5.1172452607412999e-05, - phys.calc_xs(ppid, material_view, MevEnergy{0.1})); + EXPECT_SOFT_EQ( + 5.1172452607412999e-05, + to_inv_cm(phys.calc_xs(ppid, material_view, MevEnergy{0.1}))); } //---------------------------------------------------------------------------// } // namespace test diff --git a/test/celeritas/phys/Physics.test.cu b/test/celeritas/phys/Physics.test.cu index 4349e3f2af..c4b3787f6f 100644 --- a/test/celeritas/phys/Physics.test.cu +++ b/test/celeritas/phys/Physics.test.cu @@ -10,6 +10,7 @@ #include "corecel/device_runtime_api.h" #include "corecel/sys/Device.hh" #include "corecel/sys/KernelParamCalculator.device.hh" +#include "celeritas/Quantities.hh" #include "celeritas/phys/PhysicsStepView.hh" #include "celeritas/phys/PhysicsTrackView.hh" @@ -36,7 +37,9 @@ __global__ void phys_test_kernel(PTestInput const inp) PhysicsStepView step(inp.params, inp.states, tid); phys = PhysicsTrackInitializer{}; - inp.result[tid.get()] = calc_step(phys, step, init.energy); + inp.result[tid.get()] + = native_value_to(calc_step(phys, step, init.energy)) + .value(); } } // namespace diff --git a/test/celeritas/phys/PhysicsStepUtils.test.cc b/test/celeritas/phys/PhysicsStepUtils.test.cc index 265e1280a2..3cc1d7187b 100644 --- a/test/celeritas/phys/PhysicsStepUtils.test.cc +++ b/test/celeritas/phys/PhysicsStepUtils.test.cc @@ -9,20 +9,20 @@ #include "corecel/data/CollectionStateStore.hh" #include "celeritas/MockTestBase.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/phys/CutoffParams.hh" #include "celeritas/phys/ParticleParams.hh" #include "celeritas/phys/PhysicsParams.hh" #include "DiagnosticRngEngine.hh" +#include "MockProcess.hh" #include "celeritas_test.hh" namespace celeritas { namespace test { -//---------------------------------------------------------------------------// -using units::MevEnergy; - //---------------------------------------------------------------------------// // TEST HARNESS //---------------------------------------------------------------------------// @@ -40,6 +40,8 @@ class PhysicsStepUtilsTest : public MockTestBase using PhysicsStateStore = CollectionStateStore; + using MevEnergy = units::MevEnergy; + PhysicsOptions build_physics_options() const override { return PhysicsOptions{}; @@ -123,7 +125,7 @@ TEST_F(PhysicsStepUtilsTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(discrete_action, step.action); - EXPECT_SOFT_EQ(1. / 3.e-4, step.step); + EXPECT_SOFT_EQ(1. / 3.e-4, to_cm(step.step)); } { PhysicsTrackView phys = this->init_track( @@ -132,13 +134,13 @@ TEST_F(PhysicsStepUtilsTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(discrete_action, step.action); - EXPECT_SOFT_EQ(1.e-4 / 9.e-3, step.step); + EXPECT_SOFT_EQ(1.e-4 / 9.e-3, to_cm(step.step)); // Increase the distance to interaction so range limits the step length phys.interaction_mfp(1); step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(range_action, step.action); - EXPECT_SOFT_EQ(0.48853333333333326, step.step); + EXPECT_SOFT_EQ(0.48853333333333326, to_cm(step.step)); } { PhysicsTrackView phys = this->init_track( @@ -146,7 +148,7 @@ TEST_F(PhysicsStepUtilsTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(range_action, step.action); - EXPECT_SOFT_EQ(0.0016666666666666663, step.step); + EXPECT_SOFT_EQ(0.0016666666666666663, to_cm(step.step)); } { PhysicsTrackView phys = this->init_track(&material, @@ -158,13 +160,13 @@ TEST_F(PhysicsStepUtilsTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(discrete_action, step.action); - EXPECT_SOFT_EQ(1.e-6 / 9.e-1, step.step); + EXPECT_SOFT_EQ(1.e-6 / 9.e-1, to_cm(step.step)); // Increase the distance to interaction so range limits the step length phys.interaction_mfp(1); step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(range_action, step.action); - EXPECT_SOFT_EQ(1.4285714285714282e-5, step.step); + EXPECT_SOFT_EQ(1.4285714285714282e-5, to_cm(step.step)); } { PhysicsTrackView phys = this->init_track(&material, @@ -175,7 +177,7 @@ TEST_F(PhysicsStepUtilsTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(range_action, step.action); - EXPECT_SOFT_EQ(0.014285714285714284, step.step); + EXPECT_SOFT_EQ(0.014285714285714284, to_cm(step.step)); } { PhysicsTrackView phys = this->init_track( @@ -184,13 +186,13 @@ TEST_F(PhysicsStepUtilsTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(discrete_action, step.action); - EXPECT_SOFT_EQ(1.e-4 / 9.e-3, step.step); + EXPECT_SOFT_EQ(1.e-4 / 9.e-3, to_cm(step.step)); // Increase the distance to interaction so range limits the step length phys.interaction_mfp(1); step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(range_action, step.action); - EXPECT_SOFT_EQ(0.48853333333333326, step.step); + EXPECT_SOFT_EQ(0.48853333333333326, to_cm(step.step)); } { // Test absurdly low energy (1 + E = 1) @@ -200,7 +202,7 @@ TEST_F(PhysicsStepUtilsTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(range_action, step.action); - EXPECT_SOFT_EQ(5.2704627669473019e-12, step.step); + EXPECT_SOFT_EQ(5.2704627669473019e-12, to_cm(step.step)); } { // Celerino should have infinite step with no action @@ -210,7 +212,8 @@ TEST_F(PhysicsStepUtilsTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(ActionId{}, step.action); - EXPECT_SOFT_EQ(std::numeric_limits::infinity(), step.step); + EXPECT_SOFT_EQ(std::numeric_limits::infinity(), + to_cm(step.step)); } } @@ -221,6 +224,7 @@ TEST_F(PhysicsStepUtilsTest, calc_mean_energy_loss) ParticleTrackView particle( this->particle()->host_ref(), par_state.ref(), TrackSlotId{0}); + // input: cm; output: MeV auto calc_eloss = [&](PhysicsTrackView& phys, real_type step) -> real_type { // Calculate and store the energy loss range to PhysicsTrackView auto ppid = phys.eloss_ppid(); @@ -229,8 +233,9 @@ TEST_F(PhysicsStepUtilsTest, calc_mean_energy_loss) real_type range = calc_range(particle.energy()); phys.dedx_range(range); - MevEnergy result = calc_mean_energy_loss(particle, phys, step); - return result.value(); + auto result + = calc_mean_energy_loss(particle, phys, step * units::centimeter); + return value_as(result); }; { @@ -245,7 +250,7 @@ TEST_F(PhysicsStepUtilsTest, calc_mean_energy_loss) { PhysicsTrackView phys = this->init_track( &material, MaterialId{0}, &particle, "celeriton", MevEnergy{10}); - real_type const eloss_rate = 0.2 + 0.4; + real_type const eloss_rate = (0.2 + 0.4); // MeV / cm // Tiny step: should still be linear loss (single process) EXPECT_SOFT_EQ(eloss_rate * 1e-6, calc_eloss(phys, 1e-6)); @@ -265,12 +270,12 @@ TEST_F(PhysicsStepUtilsTest, calc_mean_energy_loss) { PhysicsTrackView phys = this->init_track( &material, MaterialId{0}, &particle, "electron", MevEnergy{1e-3}); - real_type const eloss_rate = 0.5; + real_type const eloss_rate = 0.5; // MeV / cm // Low energy particle which loses all its energy over the step will // call inverse lookup. Remaining range will be zero and eloss will be // equal to the pre-step energy. - real_type step = particle.energy().value() / eloss_rate; + real_type step = value_as(particle.energy()) / eloss_rate; EXPECT_SOFT_EQ(1e-3, calc_eloss(phys, step)); } } @@ -295,7 +300,7 @@ TEST_F(PhysicsStepUtilsTest, phys.interaction_mfp(1); StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); - EXPECT_SOFT_EQ(1. / 3.e-4, step.step); + EXPECT_SOFT_EQ(1. / 3.e-4, to_cm(step.step)); // Testing cheat. PhysicsTrackView::PhysicsStateRef state_shortcut(phys_state.ref()); @@ -322,7 +327,7 @@ TEST_F(PhysicsStepUtilsTest, StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); - EXPECT_SOFT_EQ(0.48853333333333326, step.step); + EXPECT_SOFT_EQ(0.48853333333333326, to_cm(step.step)); // Testing cheat. PhysicsTrackView::PhysicsStateRef state_shortcut(phys_state.ref()); @@ -404,7 +409,7 @@ class StepLimiterTest : public PhysicsStepUtilsTest PhysicsOptions opts; opts.min_range = inf; // Use analytic range instead of scaled - opts.fixed_step_limiter = 1e-3; + opts.fixed_step_limiter = 1e-3 * units::centimeter; return opts; } }; @@ -436,7 +441,7 @@ TEST_F(StepLimiterTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(discrete_action, step.action); - EXPECT_SOFT_EQ(1. / 3.e-4, step.step); + EXPECT_SOFT_EQ(1. / 3.e-4, to_cm(step.step)); } { PhysicsTrackView phys = this->init_track( @@ -446,12 +451,12 @@ TEST_F(StepLimiterTest, calc_physics_step_limit) StepLimit step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(range_action, step.action); - EXPECT_SOFT_EQ(0.00016666666666666663, step.step); + EXPECT_SOFT_EQ(0.00016666666666666663, to_cm(step.step)); particle.energy(MevEnergy{1e-1}); step = calc_physics_step_limit(material, particle, phys, pstep); EXPECT_EQ(fixed_step_action, step.action); - EXPECT_SOFT_EQ(0.001, step.step); + EXPECT_SOFT_EQ(0.001, to_cm(step.step)); } } //---------------------------------------------------------------------------// diff --git a/test/celeritas/user/Diagnostic.test.cc b/test/celeritas/user/Diagnostic.test.cc index 817a37949d..c0b4f05e18 100644 --- a/test/celeritas/user/Diagnostic.test.cc +++ b/test/celeritas/user/Diagnostic.test.cc @@ -7,6 +7,7 @@ //---------------------------------------------------------------------------// #include "corecel/cont/Span.hh" #include "corecel/io/StringUtils.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/em/UrbanMscParams.hh" #include "celeritas/ext/GeantPhysicsOptions.hh" #include "celeritas/global/ActionRegistry.hh" @@ -38,7 +39,7 @@ class SimpleComptonDiagnosticTest : public SimpleTestBase, { Primary p; p.energy = MevEnergy{10.0}; - p.position = {-22, 0, 0}; + p.position = from_cm(Real3{-22, 0, 0}); p.direction = {1, 0, 0}; p.time = 0; std::vector result(count, p); @@ -81,7 +82,7 @@ class TestEm3DiagnosticTest : public TestEm3Base, public DiagnosticTestBase { Primary p; p.energy = MevEnergy{10.0}; - p.position = {-22, 0, 0}; + p.position = from_cm(Real3{-22, 0, 0}); p.direction = {1, 0, 0}; p.time = 0; std::vector result(count, p); diff --git a/test/celeritas/user/MctruthTestBase.cc b/test/celeritas/user/MctruthTestBase.cc index a94167be4c..490e18c761 100644 --- a/test/celeritas/user/MctruthTestBase.cc +++ b/test/celeritas/user/MctruthTestBase.cc @@ -11,6 +11,7 @@ #include "corecel/cont/Span.hh" #include "corecel/io/LogContextException.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/global/Stepper.hh" #include "celeritas/user/StepCollector.hh" @@ -91,7 +92,10 @@ auto MctruthTestBase::run(size_type num_tracks, size_type num_steps) result.track.push_back(s.track); result.step.push_back(s.step); result.volume.push_back(s.volume); - result.pos.insert(result.pos.end(), std::begin(s.pos), std::end(s.pos)); + for (auto pos_v : s.pos) + { + result.pos.push_back(to_cm(pos_v)); + } result.dir.insert(result.dir.end(), std::begin(s.dir), std::end(s.dir)); } return result; diff --git a/test/celeritas/user/StepCollector.test.cc b/test/celeritas/user/StepCollector.test.cc index 4476dc8312..906d3efcd2 100644 --- a/test/celeritas/user/StepCollector.test.cc +++ b/test/celeritas/user/StepCollector.test.cc @@ -9,6 +9,7 @@ #include "corecel/cont/Span.hh" #include "corecel/io/LogContextException.hh" +#include "celeritas/UnitUtils.hh" #include "celeritas/em/UrbanMscParams.hh" #include "celeritas/geo/GeoParams.hh" #include "celeritas/global/ActionRegistry.hh" @@ -95,7 +96,7 @@ class TestEm3CollectorTestBase : public TestEm3Base, { Primary p; p.energy = MevEnergy{10.0}; - p.position = {-22, 0, 0}; + p.position = from_cm(Real3{-22, 0, 0}); p.direction = {1, 0, 0}; p.time = 0; std::vector result(count, p);