diff --git a/app/celer-g4/SensitiveDetector.cc b/app/celer-g4/SensitiveDetector.cc index 9a0aa1a343..75e587264b 100644 --- a/app/celer-g4/SensitiveDetector.cc +++ b/app/celer-g4/SensitiveDetector.cc @@ -7,7 +7,6 @@ #include "SensitiveDetector.hh" #include -#include #include #include #include @@ -17,6 +16,7 @@ #include "corecel/Assert.hh" #include "geocel/g4/Convert.hh" +#include "celeritas/ext/GeantUnits.hh" #include "GlobalSetup.hh" @@ -89,8 +89,8 @@ bool SensitiveDetector::ProcessHits(G4Step* g4step, G4TouchableHistory*) // reloaded hit.volume = log_vol->GetInstanceID(); hit.copy_num = phys_vol->GetCopyNo(); - hit.energy_dep = convert_from_geant(edep, CLHEP::MeV); - hit.time = convert_from_geant(pre_step->GetGlobalTime(), CLHEP::s); + hit.energy_dep = units::ClhepEnergy{edep}.value(); + hit.time = native_value_from(units::ClhepTime{pre_step->GetGlobalTime()}); collection_->insert(new SensitiveHit(hit)); return true; diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index d797729cce..ac51a82f46 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -233,6 +233,9 @@ set(DOXYGEN_CLASS_GRAPH NO) set(DOXYGEN_INCLUDED_BY_GRAPH NO) set(DOXYGEN_INCLUDE_GRAPH NO) set(DOXYGEN_INTERNAL_DOCS NO) +list(APPEND DOXYGEN_EXCLUDE + "${PROJECT_SOURCE_DIR}/src/celeritas/io/" +) set(DOXYGEN_ALIASES "cite{1}=\"\\verbatim embed:rst:inline :cite:p:`\\1` \\endverbatim \"" diff --git a/doc/implementation/units-constants.rst b/doc/implementation/units-constants.rst index 50859a265a..4f2f9beb8a 100644 --- a/doc/implementation/units-constants.rst +++ b/doc/implementation/units-constants.rst @@ -50,7 +50,7 @@ using the Quantity class and helper functions. .. doxygentypedef:: celeritas::RealQuantity .. doxygenfunction:: celeritas::native_value_to -.. doxygenfunction:: celeritas::native_value_from(Quantity quant) noexcept +.. doxygenfunction:: celeritas::native_value_from .. doxygenfunction:: celeritas::value_as .. doxygenfunction:: celeritas::zero_quantity @@ -64,6 +64,11 @@ to use ``sincospi`` under the hood for improved precision. .. doxygentypedef:: celeritas::Turn_t .. doxygentypedef:: celeritas::RealTurn +Arrays of quantities can be constructed with the following helper function and +transformed to/from native values with overloads of the functions above. + +.. doxygenfunction:: celeritas::make_quantity_array + .. _api_units: Units diff --git a/src/accel/CartMapMagneticField.cc b/src/accel/CartMapMagneticField.cc index ad2f5db498..0aaf03e6ff 100644 --- a/src/accel/CartMapMagneticField.cc +++ b/src/accel/CartMapMagneticField.cc @@ -12,11 +12,10 @@ #include "corecel/Types.hh" #include "corecel/cont/Array.hh" -#include "corecel/math/ArrayUtils.hh" +#include "corecel/math/ArrayQuantity.hh" #include "geocel/GeantGeoUtils.hh" -#include "geocel/g4/Convert.hh" #include "celeritas/Types.hh" -#include "celeritas/ext/GeantUnits.hh" +#include "celeritas/UnitTypes.hh" #include "celeritas/field/CartMapField.hh" #include "celeritas/field/CartMapFieldInput.hh" #include "celeritas/field/CartMapFieldParams.hh" @@ -40,16 +39,18 @@ MakeCartMapFieldInput(G4Field const& field, CartMapFieldParams::Input field_input; // Convert from Geant4 units to native units - field_input.x.min = convert_from_geant(params.x.min, clhep_length); - field_input.x.max = convert_from_geant(params.x.max, clhep_length); + // (NOTE: params currently use real_type not double) + using ClhepLength = Quantity; + field_input.x.min = native_value_from(ClhepLength{params.x.min}); + field_input.x.max = native_value_from(ClhepLength{params.x.max}); field_input.x.num = params.x.num; - field_input.y.min = convert_from_geant(params.y.min, clhep_length); - field_input.y.max = convert_from_geant(params.y.max, clhep_length); + field_input.y.min = native_value_from(ClhepLength{params.y.min}); + field_input.y.max = native_value_from(ClhepLength{params.y.max}); field_input.y.num = params.y.num; - field_input.z.min = convert_from_geant(params.z.min, clhep_length); - field_input.z.max = convert_from_geant(params.z.max, clhep_length); + field_input.z.min = native_value_from(ClhepLength{params.z.min}); + field_input.z.max = native_value_from(ClhepLength{params.z.max}); field_input.z.num = params.z.num; // Prepare field data storage @@ -67,7 +68,7 @@ MakeCartMapFieldInput(G4Field const& field, G4double const dy = (params.y.max - params.y.min) / (params.y.num - 1); G4double const dz = (params.z.max - params.z.min) / (params.z.num - 1); - // Position calculator for Cartesian grid + // Position calculator for Cartesian grid (G4 coords) auto position_calculator = [&](size_type ix, size_type iy, size_type iz) { G4double x = params.x.min + ix * dx; G4double y = params.y.min + iy * dy; @@ -79,7 +80,9 @@ MakeCartMapFieldInput(G4Field const& field, auto field_converter = [](Array const& bfield, Array const&, real_type cur_bfield[3]) { - auto bfield_native = convert_from_geant(bfield.data(), clhep_field); + using ClhepField = Quantity; + auto bfield_native + = native_value_from(make_quantity_array(bfield)); std::copy(bfield_native.cbegin(), bfield_native.cend(), cur_bfield); }; diff --git a/src/accel/CylMapMagneticField.cc b/src/accel/CylMapMagneticField.cc index da64af2d93..d3c409f9c2 100644 --- a/src/accel/CylMapMagneticField.cc +++ b/src/accel/CylMapMagneticField.cc @@ -65,10 +65,11 @@ MakeCylMapFieldInput(G4Field const& field, field_input.grid_z.reserve(z_grid.size()); // Convert from geant units to native units + using ClhepLength = Quantity; std::transform(r_grid.cbegin(), r_grid.cend(), std::back_inserter(field_input.grid_r), - [](auto r) { return convert_from_geant(r, clhep_length); }); + [](auto r) { return native_value_from(ClhepLength{r}); }); // Convert phi values to Turn type std::transform(phi_values.cbegin(), phi_values.cend(), @@ -77,7 +78,7 @@ MakeCylMapFieldInput(G4Field const& field, std::transform(z_grid.cbegin(), z_grid.cend(), std::back_inserter(field_input.grid_z), - [](auto z) { return convert_from_geant(z, clhep_length); }); + [](auto z) { return native_value_from(ClhepLength{z}); }); size_type const nr = field_input.grid_r.size(); size_type const nphi = field_input.grid_phi.size(); @@ -104,13 +105,13 @@ MakeCylMapFieldInput(G4Field const& field, auto field_converter = [](Array const& bfield, Array const& pos, real_type cur_bfield[3]) { + using ClhepField = Quantity; + auto bfield_native + = native_value_from(make_quantity_array(bfield)); double const phi = std::atan2(pos[1], pos[0]); EnumArray bfield_cyl; - cartesian_to_cylindrical(bfield, phi, bfield_cyl); - auto bfield_cyl_native - = convert_from_geant(bfield_cyl.data(), clhep_field); - std::copy( - bfield_cyl_native.cbegin(), bfield_cyl_native.cend(), cur_bfield); + cartesian_to_cylindrical(bfield_native, phi, bfield_cyl); + std::copy(bfield_cyl.cbegin(), bfield_cyl.cend(), cur_bfield); }; // Sample field using common utility diff --git a/src/accel/HepMC3PrimaryGenerator.cc b/src/accel/HepMC3PrimaryGenerator.cc index 6756b5c809..761a7fa6cd 100644 --- a/src/accel/HepMC3PrimaryGenerator.cc +++ b/src/accel/HepMC3PrimaryGenerator.cc @@ -62,8 +62,8 @@ class PrimaryInserter G4PrimaryParticle* primary = new G4PrimaryParticle(par.pid()); // Set the primary directlon - auto dir = make_unit_vector(Array{p.x(), p.y(), p.z()}); - primary->SetMomentumDirection(convert_to_geant(dir, 1)); + auto dir = make_unit_vector(Array{p.x(), p.y(), p.z()}); + primary->SetMomentumDirection(to_g4vector(dir)); // Set the kinetic energy primary->SetKineticEnergy(p.e() - p.m()); diff --git a/src/accel/LocalTransporter.cc b/src/accel/LocalTransporter.cc index 85a33b6848..7456e8608f 100644 --- a/src/accel/LocalTransporter.cc +++ b/src/accel/LocalTransporter.cc @@ -284,7 +284,8 @@ void LocalTransporter::Push(G4Track& g4track) << "' particles"); track.position = convert_from_geant(g4track.GetPosition(), clhep_length); - track.direction = convert_from_geant(g4track.GetMomentumDirection(), 1); + track.direction = static_array_cast( + to_array(g4track.GetMomentumDirection())); track.time = convert_from_geant(g4track.GetGlobalTime(), clhep_time); track.weight = g4track.GetWeight(); diff --git a/src/accel/PGPrimaryGeneratorAction.cc b/src/accel/PGPrimaryGeneratorAction.cc index b8bf392d5c..96f406ef89 100644 --- a/src/accel/PGPrimaryGeneratorAction.cc +++ b/src/accel/PGPrimaryGeneratorAction.cc @@ -87,8 +87,10 @@ void PGPrimaryGeneratorAction::GeneratePrimaries(G4Event* event) gun_.SetParticleDefinition( G4ParticleTable::GetParticleTable()->FindParticle( pdg_[p.particle_id.unchecked_get()].get())); - gun_.SetParticlePosition(convert_to_geant(p.position, clhep_length)); - gun_.SetParticleMomentumDirection(convert_to_geant(p.direction, 1)); + gun_.SetParticlePosition( + convert_to_geant(p.position)); + gun_.SetParticleMomentumDirection( + to_g4vector(static_array_cast(p.direction))); gun_.SetParticleEnergy(convert_to_geant(p.energy, CLHEP::MeV)); gun_.GeneratePrimaryVertex(event); diff --git a/src/celeritas/Quantities.hh b/src/celeritas/Quantities.hh index e16bf7ec91..c0ea64d0c3 100644 --- a/src/celeritas/Quantities.hh +++ b/src/celeritas/Quantities.hh @@ -30,6 +30,15 @@ using AmuMass = RealQuantity; //! Special faux quantity for overloading cross section calculation using LogMevEnergy = RealQuantity; +//---------------------------------------------------------------------------// +//!@{ +//! \name Quantities for Geant4 conversion +using ClhepLength = Quantity; +using ClhepField = Quantity; +using ClhepTime = Quantity; +using ClhepEnergy = Quantity; +//!@} + //---------------------------------------------------------------------------// //!@{ //! \name Quantities for manual input and/or test harnesses diff --git a/src/celeritas/ext/GeantUnits.hh b/src/celeritas/ext/GeantUnits.hh index b0d46f282d..d3b01e0071 100644 --- a/src/celeritas/ext/GeantUnits.hh +++ b/src/celeritas/ext/GeantUnits.hh @@ -14,7 +14,7 @@ namespace celeritas { //---------------------------------------------------------------------------// -// CONSTANTS +// DEPRECATED: use ClhepField, ClhepTime //---------------------------------------------------------------------------// //! Value of a unit Celeritas field in the CLHEP unit system inline constexpr double clhep_field{1 / units::ClhepTraits::BField::value()}; diff --git a/src/celeritas/ext/detail/HitProcessor.cc b/src/celeritas/ext/detail/HitProcessor.cc index ba738f54b9..94b342fcb7 100644 --- a/src/celeritas/ext/detail/HitProcessor.cc +++ b/src/celeritas/ext/detail/HitProcessor.cc @@ -185,8 +185,7 @@ HitProcessor::HitProcessor(SPConstVecLV detector_volumes, detectors_[i] = lv->GetSensitiveDetector(); CELER_VALIDATE(detectors_[i], << "no sensitive detector is attached to volume '" - << lv->GetName() << "'@" - << static_cast(lv)); + << StreamableLV{lv}); } CELER_ENSURE(!detectors_.empty()); diff --git a/src/celeritas/ext/detail/NaviTouchableUpdater.cc b/src/celeritas/ext/detail/NaviTouchableUpdater.cc index b1abe63923..4c487b4f39 100644 --- a/src/celeritas/ext/detail/NaviTouchableUpdater.cc +++ b/src/celeritas/ext/detail/NaviTouchableUpdater.cc @@ -94,8 +94,9 @@ bool NaviTouchableUpdater::operator()(Real3 const& pos, CELER_EXPECT(lv); CELER_EXPECT(touchable); - auto g4pos = convert_to_geant(pos, clhep_length); - auto g4dir = convert_to_geant(dir, 1); + auto g4pos + = convert_to_geant(static_array_cast(pos)); + auto g4dir = to_g4vector(static_array_cast(dir)); // Locate pre-step point navi_->LocateGlobalPointAndUpdateTouchable(g4pos, diff --git a/src/celeritas/field/UniformFieldParams.cc b/src/celeritas/field/UniformFieldParams.cc index 64f3573edc..91e9e338bf 100644 --- a/src/celeritas/field/UniformFieldParams.cc +++ b/src/celeritas/field/UniformFieldParams.cc @@ -16,6 +16,7 @@ #include "corecel/cont/VariantUtils.hh" #include "corecel/data/CollectionBuilder.hh" #include "corecel/math/Algorithms.hh" +#include "corecel/math/ArrayQuantity.hh" #include "corecel/math/ArrayUtils.hh" #include "geocel/VolumeCollectionBuilder.hh" #include "geocel/VolumeIdBuilder.hh" @@ -69,10 +70,8 @@ validated_field_data(UniformFieldParams::Input const& inp) validate_input(inp.driver_options); HostVal result; - for (auto i : range(inp.strength.size())) - { - result.field[i] = native_value_from(units::FieldTesla{inp.strength[i]}); - } + result.field = native_value_from( + make_quantity_array(inp.strength)); result.options = inp.driver_options; return result; } diff --git a/src/celeritas/phys/PrimaryGenerator.cc b/src/celeritas/phys/PrimaryGenerator.cc index 1a70eac72d..9030845377 100644 --- a/src/celeritas/phys/PrimaryGenerator.cc +++ b/src/celeritas/phys/PrimaryGenerator.cc @@ -38,8 +38,7 @@ auto make_energy_sampler(inp::EnergyDistribution const& i) return_as(Overload{ [](inp::MonoenergeticDistribution const& me) { CELER_VALIDATE(me.value > 0, - << "invalid primary generator " - "energy " + << R"(invalid primary generator energy )" << me.value); return DeltaDistribution{static_cast(me.value)}; }, @@ -56,16 +55,17 @@ auto make_energy_sampler(inp::EnergyDistribution const& i) auto make_position_sampler(inp::ShapeDistribution const& i) { CELER_ASSUME(!i.valueless_by_exception()); - return std::visit( - return_as(Overload{ - [](inp::PointDistribution const& ps) { - return DeltaDistribution{array_cast(ps.value)}; - }, - [](inp::UniformBoxDistribution const& ubs) { - return UniformBoxDistribution{array_cast(ubs.lower), - array_cast(ubs.upper)}; - }}), - i); + return std::visit(return_as(Overload{ + [](inp::PointDistribution const& ps) { + return DeltaDistribution{ + static_array_cast(ps.value)}; + }, + [](inp::UniformBoxDistribution const& ubs) { + return UniformBoxDistribution{ + static_array_cast(ubs.lower), + static_array_cast(ubs.upper)}; + }}), + i); } //---------------------------------------------------------------------------// @@ -81,10 +81,11 @@ auto make_direction_sampler(inp::AngleDistribution const& i) return IsotropicDistribution{}; }, [](inp::MonodirectionalDistribution const& ma) { - CELER_VALIDATE(is_soft_unit_vector(ma.value), - << "primary generator angle is " - "not a unit vector"); - return DeltaDistribution{array_cast(ma.value)}; + CELER_VALIDATE( + is_soft_unit_vector(ma.value), + << R"(primary generator angle is not a unit vector)"); + return DeltaDistribution{ + static_array_cast(ma.value)}; }}), i); } diff --git a/src/corecel/cont/Array.hh b/src/corecel/cont/Array.hh index ec467aeace..924410cdde 100644 --- a/src/corecel/cont/Array.hh +++ b/src/corecel/cont/Array.hh @@ -179,6 +179,45 @@ CELER_CEF bool operator!=(Array const& lhs, Array const& rhs) return !(lhs == rhs); } +//---------------------------------------------------------------------------// +// FREE FUNCTIONS +//---------------------------------------------------------------------------// +/*! + * Convert a C array from type \c T2 to an Array of \c T1. + * + * The standard library version of this function is available since C++20. + */ +template +CELER_CONSTEXPR_FUNCTION auto to_array(T (&x)[N]) +{ + static_assert(!std::is_array_v, + "to_array elements cannot be multidimensional"); + static_assert(std::is_constructible_v, + "to_array elements must be copy constructible"); + + Array, N> result; + for (size_type i = 0; i != N; ++i) + { + result[i] = x[i]; + } + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Convert an array from type \c T2 to \c T1. + */ +template +CELER_CONSTEXPR_FUNCTION Array static_array_cast(Array const& x) +{ + Array result; + for (size_type i = 0; i != N; ++i) + { + result[i] = static_cast(x[i]); + } + return result; +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/corecel/math/ArrayOperators.hh b/src/corecel/math/ArrayOperators.hh index 8cb5f53480..fd13d268f4 100644 --- a/src/corecel/math/ArrayOperators.hh +++ b/src/corecel/math/ArrayOperators.hh @@ -22,44 +22,44 @@ namespace celeritas { //---------------------------------------------------------------------------// -#define CELER_DEFINE_ARRAY_ASSIGN(TOKEN) \ - template \ - inline CELER_FUNCTION Array& operator TOKEN(Array& x, \ - Array const& y) \ - { \ - for (size_type i = 0; i != N; ++i) \ - { \ - x[i] TOKEN y[i]; \ - } \ - return x; \ - } \ - \ - template> \ - inline CELER_FUNCTION Array& operator TOKEN(Array& x, \ - T2 const& y) \ - { \ - for (size_type i = 0; i != N; ++i) \ - { \ - x[i] TOKEN y; \ - } \ - return x; \ +#define CELER_DEFINE_ARRAY_ASSIGN(TOKEN) \ + template \ + CELER_CONSTEXPR_FUNCTION Array& operator TOKEN( \ + Array& x, Array const& y) \ + { \ + for (size_type i = 0; i != N; ++i) \ + { \ + x[i] TOKEN y[i]; \ + } \ + return x; \ + } \ + \ + template> \ + CELER_CONSTEXPR_FUNCTION Array& operator TOKEN(Array& x, \ + T2 const& y) \ + { \ + for (size_type i = 0; i != N; ++i) \ + { \ + x[i] TOKEN y; \ + } \ + return x; \ } -#define CELER_DEFINE_ARRAY_ARITHM(TOKEN) \ - template \ - inline CELER_FUNCTION Array operator TOKEN(Array const& x, \ - Array const& y) \ - { \ - Array result{x}; \ - return (result TOKEN## = y); \ - } \ - \ - template> \ - inline CELER_FUNCTION Array operator TOKEN(Array const& x, \ - T2 const& y) \ - { \ - Array result{x}; \ - return (result TOKEN## = y); \ +#define CELER_DEFINE_ARRAY_ARITHM(TOKEN) \ + template \ + CELER_CONSTEXPR_FUNCTION Array operator TOKEN(Array const& x, \ + Array const& y) \ + { \ + Array result{x}; \ + return (result TOKEN## = y); \ + } \ + \ + template> \ + CELER_CONSTEXPR_FUNCTION Array operator TOKEN(Array const& x, \ + T2 const& y) \ + { \ + Array result{x}; \ + return (result TOKEN## = y); \ } //---------------------------------------------------------------------------// @@ -82,7 +82,8 @@ CELER_DEFINE_ARRAY_ARITHM(/) //! Left-multiply by scalar template> -inline CELER_FUNCTION Array operator*(T2 const& y, Array const& x) +CELER_CONSTEXPR_FUNCTION Array +operator*(T2 const& y, Array const& x) { return x * y; } @@ -92,7 +93,7 @@ inline CELER_FUNCTION Array operator*(T2 const& y, Array const& x) * Unary negation. */ template -inline CELER_FUNCTION Array operator-(Array const& x) +CELER_CONSTEXPR_FUNCTION Array operator-(Array const& x) { Array result; for (size_type i = 0; i != N; ++i) diff --git a/src/corecel/math/ArrayQuantity.hh b/src/corecel/math/ArrayQuantity.hh new file mode 100644 index 0000000000..848a49cbeb --- /dev/null +++ b/src/corecel/math/ArrayQuantity.hh @@ -0,0 +1,115 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file corecel/math/ArrayQuantity.hh +//! \brief Create and convert arrays of quantities +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/Types.hh" +#include "corecel/cont/Array.hh" + +#include "Quantity.hh" + +namespace celeritas +{ + +//---------------------------------------------------------------------------// +/*! + * Construct an array of quantities from raw values. + * + * This helper function allows concise construction of arrays of quantities: + * \code + * auto distances = make_quantity_array(1.0, 2.5, 3.7); + * \endcode + */ +template +CELER_CONSTEXPR_FUNCTION Array +make_quantity_array(Args const&... args) noexcept +{ + return {Q{args}...}; +} + +//! \cond (CELERITAS_DOC_DEV) +//---------------------------------------------------------------------------// +/*! + * Construct an array of quantities from raw values. + * + * This helper function allows concise construction of arrays of quantities: + * \code + * auto pos = make_quantity_array(hardcoded_pos_cm); + * \endcode + */ +template +CELER_CONSTEXPR_FUNCTION Array +make_quantity_array(Array const& arr) noexcept +{ + Array result; + for (size_type i = 0; i < N; ++i) + { + result[i] = Q{arr[i]}; + } + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Convert an array of quantities to native values. + * + * This applies native_value_from element-wise to each component. + */ +template +CELER_CONSTEXPR_FUNCTION auto +native_value_from(Array, N> const& quant) noexcept +{ + using common_type = typename Quantity::common_type; + Array result; + for (size_type i = 0; i < N; ++i) + { + result[i] = native_value_from(quant[i]); + } + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Convert an array of native values to an array of quantities. + * + * This applies native_value_to element-wise to each component. + */ +template +CELER_CONSTEXPR_FUNCTION auto native_value_to(Array const& value) noexcept +{ + Array result; + for (size_type i = 0; i < N; ++i) + { + result[i] = native_value_to(value[i]); + } + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Convert an array of quantities to native values. + * + * This applies native_value_from element-wise to each component. + */ +template +CELER_CONSTEXPR_FUNCTION auto value_as(Array const& quant) noexcept + -> std::enable_if_t, Array> +{ + Array result; + for (size_type i = 0; i < N; ++i) + { + result[i] = value_as(quant[i]); + } + return result; +} + +//! \endcond +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/corecel/math/ArrayUtils.hh b/src/corecel/math/ArrayUtils.hh index 14f68d79b9..37b01daf22 100644 --- a/src/corecel/math/ArrayUtils.hh +++ b/src/corecel/math/ArrayUtils.hh @@ -85,12 +85,6 @@ template [[nodiscard]] inline CELER_FUNCTION Array rotate(Array const& dir, Array const& rot); -//---------------------------------------------------------------------------// -// Convert an array from type \c T2 to \c T1 -template -[[nodiscard]] inline CELER_FUNCTION Array -array_cast(Array const& x); - //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// @@ -348,20 +342,5 @@ rotate(Array const& dir, Array const& rot) return make_unit_vector(result); } -//---------------------------------------------------------------------------// -/*! - * Convert an array from type \c T2 to \c T1. - */ -template -CELER_FUNCTION Array array_cast(Array const& x) -{ - Array result; - for (size_type i = 0; i != N; ++i) - { - result[i] = static_cast(x[i]); - } - return result; -} - //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/corecel/random/distribution/DistributionInserter.cc b/src/corecel/random/distribution/DistributionInserter.cc index da400b745c..23e17e7578 100644 --- a/src/corecel/random/distribution/DistributionInserter.cc +++ b/src/corecel/random/distribution/DistributionInserter.cc @@ -55,7 +55,7 @@ ThreedDistributionId DistributionInserter::operator()( inp::DeltaDistribution> const& d) { DeltaDistributionRecord> record; - record.value = array_cast(d.value); + record.value = static_array_cast(d.value); auto id = CollectionBuilder{&data_.delta_real3}.push_back(record); return (*this)(ThreedDistributionType::delta, id.get()); } @@ -80,8 +80,8 @@ ThreedDistributionId DistributionInserter::operator()(inp::UniformBoxDistribution const& d) { UniformBoxDistributionRecord record; - record.upper = array_cast(d.upper); - record.lower = array_cast(d.lower); + record.upper = static_array_cast(d.upper); + record.lower = static_array_cast(d.lower); auto id = CollectionBuilder{&data_.uniform_box}.push_back(record); return (*this)(ThreedDistributionType::uniform_box, id.get()); } diff --git a/src/geocel/GeantGeoParams.cc b/src/geocel/GeantGeoParams.cc index 0694c0ee25..861a50553e 100644 --- a/src/geocel/GeantGeoParams.cc +++ b/src/geocel/GeantGeoParams.cc @@ -939,9 +939,14 @@ void GeantGeoParams::build_metadata() "impl volume", make_logical_vol_labels(vi_mapper_, this->lv_offset())}; surfaces_ = make_surface_vec(*this); + using lengthunits::ClhepLength; auto clhep_bbox = this->get_clhep_bbox(); - bbox_ = {convert_from_geant(clhep_bbox.lower().data(), clhep_length), - convert_from_geant(clhep_bbox.upper().data(), clhep_length)}; + auto to_native_real3 = [](Array const& arr) { + return static_array_cast( + native_value_from(make_quantity_array(arr))); + }; + bbox_ = {to_native_real3(clhep_bbox.lower()), + to_native_real3(clhep_bbox.upper())}; CELER_ENSURE(bbox_); CELER_ENSURE(data_); } diff --git a/src/geocel/GeantGeoParams.hh b/src/geocel/GeantGeoParams.hh index 5476c1488d..9d79626ce5 100644 --- a/src/geocel/GeantGeoParams.hh +++ b/src/geocel/GeantGeoParams.hh @@ -157,7 +157,7 @@ class GeantGeoParams final : public GeoParamsInterface, G4VPhysicalVolume* world() { return data_.world; } //!@} - // Get the world extents in Geant4 units + // Get the world extents in Geant4 units and precision BoundingBox get_clhep_bbox() const; //// DATA ACCESS //// diff --git a/src/geocel/detail/LengthQuantities.hh b/src/geocel/detail/LengthQuantities.hh new file mode 100644 index 0000000000..246040eee0 --- /dev/null +++ b/src/geocel/detail/LengthQuantities.hh @@ -0,0 +1,55 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file geocel/detail/LengthQuantities.hh +//! \brief Length unit types and quantity aliases for geometry interop +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/math/ArrayQuantity.hh" // IWYU pragma: export +#include "corecel/math/Quantity.hh" // IWYU pragma: export + +#include "LengthUnits.hh" + +namespace celeritas +{ +namespace lengthunits +{ +//---------------------------------------------------------------------------// +//!@{ +//! \name Length unit types + +//! Centimeter length (CGS/Gaussian units) +struct Centimeter +{ + static CELER_CONSTEXPR_FUNCTION Constant value() + { + return lengthunits::centimeter; + } + static char const* label() { return "cm"; } +}; + +//! Millimeter length (CLHEP units) +struct Millimeter +{ + static CELER_CONSTEXPR_FUNCTION Constant value() + { + return lengthunits::millimeter; + } + static char const* label() { return "mm"; } +}; + +//!@} + +//---------------------------------------------------------------------------// +//!@{ +//! \name Length quantity aliases + +using CmLength = RealQuantity; +using ClhepLength = Quantity; + +//!@} +//---------------------------------------------------------------------------// +} // namespace lengthunits +} // namespace celeritas diff --git a/src/geocel/detail/LengthUnits.hh b/src/geocel/detail/LengthUnits.hh index 2adbdc7f8f..b40f6dbc5b 100644 --- a/src/geocel/detail/LengthUnits.hh +++ b/src/geocel/detail/LengthUnits.hh @@ -9,7 +9,6 @@ #include "corecel/Config.hh" -#include "corecel/Types.hh" #include "corecel/math/Constant.hh" namespace celeritas diff --git a/src/geocel/g4/Convert.hh b/src/geocel/g4/Convert.hh index 88d54db97f..853622fc6c 100644 --- a/src/geocel/g4/Convert.hh +++ b/src/geocel/g4/Convert.hh @@ -6,35 +6,54 @@ //---------------------------------------------------------------------------// #pragma once -#include +#include "corecel/math/ArrayOperators.hh" +#include "corecel/math/ArrayQuantity.hh" +#include "geocel/detail/LengthQuantities.hh" -#include "corecel/Types.hh" -#include "corecel/cont/Array.hh" -#include "corecel/math/Constant.hh" -#include "geocel/Types.hh" -#include "geocel/detail/LengthUnits.hh" +#include "GeantTypes.hh" namespace celeritas { //---------------------------------------------------------------------------// -// CONSTANTS +// FREE FUNCTIONS //---------------------------------------------------------------------------// -//! Value of a unit Celeritas length in the CLHEP unit system -inline constexpr double clhep_length{1 / lengthunits::millimeter}; +//! Convert via a quantity to native Geant4 types/units +template +G4ThreeVector convert_to_geant(Array const& v) +{ + return to_g4vector(value_as(native_value_to(v))); +} //---------------------------------------------------------------------------// -// FREE FUNCTIONS +//! Convert via a quantity to native Geant4 types/units +template +Array convert_from_geant(G4ThreeVector const& v) +{ + return native_value_from(make_quantity_array(to_array(v))); +} + //---------------------------------------------------------------------------// -/*! - * Convert a value from Geant4/CLHEP to Celeritas native units. - */ -template -constexpr inline T convert_from_geant(T const& val, T units) +//! Convert via a quantity to native Geant4 types/units +template +double convert_to_geant(real_type v) { - return val / units; + return value_as(native_value_to(v)); +} + +//---------------------------------------------------------------------------// +//! Convert via a quantity to native Geant4 types/units +template +real_type convert_from_geant(double v) +{ + return native_value_from(Q(v)); } //---------------------------------------------------------------------------// +// DEPRECATED +//---------------------------------------------------------------------------// +//! Value of a unit Celeritas length in the CLHEP unit system +inline constexpr double clhep_length{1 / lengthunits::millimeter}; + /*! * Convert a value from Geant4 with CLHEP units. */ @@ -49,9 +68,7 @@ constexpr inline double convert_from_geant(double val, double units) */ inline Real3 convert_from_geant(G4ThreeVector const& vec, double units) { - return {static_cast(vec[0] / units), - static_cast(vec[1] / units), - static_cast(vec[2] / units)}; + return static_array_cast(to_array(vec) / units); } //---------------------------------------------------------------------------// @@ -60,54 +77,26 @@ inline Real3 convert_from_geant(G4ThreeVector const& vec, double units) */ inline Real3 convert_from_geant(double const vec[3], double units) { - return {static_cast(vec[0] / units), - static_cast(vec[1] / units), - static_cast(vec[2] / units)}; + return static_array_cast(Array{vec[0], vec[1], vec[2]} / units); } //---------------------------------------------------------------------------// /*! * Convert a native Celeritas quantity to a Geant4 value with CLHEP units. */ -template -constexpr inline double convert_to_geant(T const& val, double units) -{ - return val * units; -} - -//---------------------------------------------------------------------------// -/*! - * Convert a native Celeritas quantity to a Geant4 value. - */ constexpr inline double convert_to_geant(real_type val, double units) { - return double{val} * units; + return val * units; } //---------------------------------------------------------------------------// /*! * Convert a native Celeritas 3-vector to a Geant4 equivalent. */ -template -inline G4ThreeVector convert_to_geant(Array const& arr, double units) -{ - return {arr[0] * units, arr[1] * units, arr[2] * units}; -} - -//! \cond (CELERITAS_DOC_DEV) -//---------------------------------------------------------------------------// -/*! - * Set y += a * x . - */ -inline void axpy(double a, G4ThreeVector const& x, G4ThreeVector* y) +inline G4ThreeVector convert_to_geant(Real3 const& arr, double units) { - CELER_EXPECT(y); - for (int i = 0; i < 3; ++i) - { - (*y)[i] = a * x[i] + (*y)[i]; - } + return to_g4vector(static_array_cast(arr) * units); } -//! \endcond //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/geocel/g4/GeantGeoTrackView.hh b/src/geocel/g4/GeantGeoTrackView.hh index aa72f715a4..a29d0f0e2b 100644 --- a/src/geocel/g4/GeantGeoTrackView.hh +++ b/src/geocel/g4/GeantGeoTrackView.hh @@ -141,6 +141,7 @@ class GeantGeoTrackView TrackSlotId parent; //!< Parent track with existing geometry ::celeritas::Real3 const& dir; //!< New direction }; + using ClhepLength = lengthunits::ClhepLength; //// DATA //// @@ -197,9 +198,9 @@ GeantGeoTrackView::GeantGeoTrackView(ParamsRef const& params, , safety_radius_(states.safety_radius[tid]) , touch_handle_(states.nav_state.touch_handle(tid)) , navi_(states.nav_state.navigator(tid)) - , g4pos_(convert_to_geant(pos_, clhep_length)) - , g4dir_(convert_to_geant(dir_, 1)) - , g4safety_(convert_to_geant(safety_radius_, clhep_length)) + , g4pos_(convert_to_geant(pos_)) + , g4dir_(to_g4vector(dir_)) + , g4safety_(convert_to_geant(safety_radius_)) { } @@ -224,8 +225,8 @@ GeantGeoTrackView& GeantGeoTrackView::operator=(Initializer_t const& init) next_step_ = 0; safety_radius_ = -1; // Assume *not* on a boundary - g4pos_ = convert_to_geant(pos_, clhep_length); - g4dir_ = convert_to_geant(dir_, 1); + g4pos_ = convert_to_geant(pos_); + g4dir_ = to_g4vector(dir_); g4safety_ = -1; navi_.LocateGlobalPointAndUpdateTouchable(g4pos_, @@ -377,8 +378,7 @@ auto GeantGeoTrackView::normal() const -> Real3 G4ThreeVector norm = navi_.GetGlobalExitNormal(g4pos_, &valid); CELER_ASSERT(valid); - // TODO: convert_from_geant uses celeritas::real_type, not double - Real3 result{norm[0], norm[1], norm[2]}; + Real3 result = to_array(norm); CELER_ENSURE(is_soft_unit_vector(result)); return result; } @@ -406,20 +406,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_length); + real_type g4step = convert_to_geant(max_step); 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_length); + safety_radius_ = native_value_from(ClhepLength{g4safety_}); CELER_ASSERT(!this->is_on_boundary()); } // Update result Propagation result; - result.distance = convert_from_geant(g4step, clhep_length); + result.distance = native_value_from(ClhepLength{g4step}); if (result.distance <= max_step) { result.boundary = true; @@ -467,9 +467,9 @@ auto GeantGeoTrackView::find_safety(real_type max_step) -> real_type CELER_EXPECT(max_step > 0); if (safety_radius_ < max_step) { - real_type g4step = convert_to_geant(max_step, clhep_length); + real_type g4step = convert_to_geant(max_step); g4safety_ = navi_.ComputeSafety(g4pos_, g4step); - safety_radius_ = max(convert_from_geant(g4safety_, clhep_length), 0.0); + safety_radius_ = max(native_value_from(ClhepLength{g4safety_}), 0.0); } return safety_radius_; @@ -485,7 +485,7 @@ void GeantGeoTrackView::move_to_boundary() // Move next step axpy(next_step_, dir_, &pos_); - axpy(convert_to_geant(next_step_, clhep_length), g4dir_, &g4pos_); + axpy(convert_to_geant(next_step_), g4dir_, &g4pos_); next_step_ = 0; safety_radius_ = 0; g4safety_ = 0; @@ -531,7 +531,7 @@ void GeantGeoTrackView::move_internal(real_type dist) // Move and update next_step axpy(dist, dir_, &pos_); - axpy(convert_to_geant(dist, clhep_length), g4dir_, &g4pos_); + axpy(convert_to_geant(dist), g4dir_, &g4pos_); next_step_ -= dist; navi_.LocateGlobalPointWithinVolume(g4pos_); @@ -549,7 +549,7 @@ void GeantGeoTrackView::move_internal(real_type dist) void GeantGeoTrackView::move_internal(Real3 const& pos) { pos_ = pos; - g4pos_ = convert_to_geant(pos_, clhep_length); + g4pos_ = convert_to_geant(pos_); next_step_ = 0; navi_.LocateGlobalPointWithinVolume(g4pos_); @@ -568,7 +568,7 @@ void GeantGeoTrackView::set_dir(Real3 const& newdir) { CELER_EXPECT(is_soft_unit_vector(newdir)); dir_ = newdir; - g4dir_ = convert_to_geant(newdir, 1); + g4dir_ = to_g4vector(newdir); next_step_ = 0; } diff --git a/src/geocel/g4/GeantTypes.hh b/src/geocel/g4/GeantTypes.hh new file mode 100644 index 0000000000..98ded397fe --- /dev/null +++ b/src/geocel/g4/GeantTypes.hh @@ -0,0 +1,46 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file geocel/g4/GeantTypes.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/Assert.hh" +#include "corecel/Types.hh" +#include "corecel/cont/Array.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +//! Convert a Geant4 vector to a Celeritas array of the same precision +inline Array to_array(G4ThreeVector const& v) +{ + return {v.x(), v.y(), v.z()}; +} + +//---------------------------------------------------------------------------// +//! Convert a native Celeritas 3-vector to a Geant4 equivalent +inline G4ThreeVector to_g4vector(Array const& arr) +{ + return {arr[0], arr[1], arr[2]}; +} + +//! \cond (CELERITAS_DOC_DEV) +//---------------------------------------------------------------------------// +//! Let y <- a * x + y +inline void axpy(double a, G4ThreeVector const& x, G4ThreeVector* y) +{ + CELER_EXPECT(y); + for (int i = 0; i < 3; ++i) + { + (*y)[i] = std::fma(a, x[i], (*y)[i]); + } +} +//! \endcond + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/geocel/vg/VecgeomParams.cc b/src/geocel/vg/VecgeomParams.cc index e436d36004..30a9846971 100644 --- a/src/geocel/vg/VecgeomParams.cc +++ b/src/geocel/vg/VecgeomParams.cc @@ -54,8 +54,8 @@ #include "geocel/detail/MakeLabelVector.hh" #include "VecgeomData.hh" // IWYU pragma: associated +#include "VecgeomTypes.hh" -#include "detail/VecgeomCompatibility.hh" #include "detail/VecgeomSetup.hh" static_assert(std::is_same_v, @@ -585,7 +585,7 @@ VecgeomParams::VecgeomParams(vecgeom::GeoManager const& geo, auto bbox_mgr = ABBoxManager_t::Instance(); VgReal3 lower, upper; bbox_mgr.ComputeABBox(&world, &lower, &upper); - return BBox{detail::to_array(lower), detail::to_array(upper)}; + return BBox{to_array(lower), to_array(upper)}; }(); } diff --git a/src/geocel/vg/VecgeomTrackView.hh b/src/geocel/vg/VecgeomTrackView.hh index 682694767f..e8bb6293d8 100644 --- a/src/geocel/vg/VecgeomTrackView.hh +++ b/src/geocel/vg/VecgeomTrackView.hh @@ -25,8 +25,6 @@ #include "VecgeomData.hh" #include "VecgeomTypes.hh" -#include "detail/VecgeomCompatibility.hh" - #if CELERITAS_VECGEOM_SURFACE # include "detail/SurfNavigator.hh" #elif CELERITAS_VECGEOM_VERSION < 0x020000 @@ -292,7 +290,7 @@ VecgeomTrackView::operator=(Initializer_t const& init) // LocatePointIn sets `vgstate_` constexpr bool contains_point = true; Navigator::LocatePointIn( - world, detail::to_vector(pos_), vgstate_, contains_point); + world, to_vgvector(pos_), vgstate_, contains_point); return *this; } @@ -434,9 +432,8 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step() { // SPECIAL CASE: find distance to interior from outside world volume auto const& world_pv = this->world(); - next_step_ = world_pv.DistanceToIn(detail::to_vector(pos_), - detail::to_vector(dir_), - vecgeom::kInfLength); + next_step_ = world_pv.DistanceToIn( + to_vgvector(pos_), to_vgvector(dir_), vecgeom::kInfLength); next_step_ = max(next_step_, this->extra_push()); vgnext_.Clear(); @@ -470,8 +467,8 @@ CELER_FUNCTION Propagation VecgeomTrackView::find_next_step(real_type max_step) *next_surf_ = null_surface(); } // TODO: vgnext is simply copied and the boundary flag optionally set - next_step_ = Navigator::ComputeStepAndNextVolume(detail::to_vector(pos_), - detail::to_vector(dir_), + next_step_ = Navigator::ComputeStepAndNextVolume(to_vgvector(pos_), + to_vgvector(dir_), max_step, vgstate_, vgnext_ @@ -532,7 +529,7 @@ CELER_FUNCTION real_type VecgeomTrackView::find_safety(real_type max_radius) CELER_EXPECT(!this->is_on_boundary()); CELER_EXPECT(max_radius > 0); - real_type safety = Navigator::ComputeSafety(detail::to_vector(this->pos()), + real_type safety = Navigator::ComputeSafety(to_vgvector(this->pos()), vgstate_ #if VECGEOM_VERSION >= 0x200000 , @@ -580,8 +577,8 @@ CELER_FUNCTION void VecgeomTrackView::cross_boundary() if (!CELERITAS_VECGEOM_SURFACE || !vgstate_.IsOutside()) { // In surf model, relocation does not work from [OUTSIDE] - Navigator::RelocateToNextVolume(detail::to_vector(this->pos_), - detail::to_vector(this->dir_), + Navigator::RelocateToNextVolume(to_vgvector(this->pos_), + to_vgvector(this->dir_), #if CELERITAS_VECGEOM_SURFACE *next_surf_, #endif diff --git a/src/geocel/vg/VecgeomTypes.hh b/src/geocel/vg/VecgeomTypes.hh index 405a618c31..31e3043b55 100644 --- a/src/geocel/vg/VecgeomTypes.hh +++ b/src/geocel/vg/VecgeomTypes.hh @@ -135,5 +135,29 @@ using VgNavStateImpl = VgNavIndex; //! High level VecGeom navigation state using VgNavState = vecgeom::NavigationState; +//---------------------------------------------------------------------------// +// CONVERSION FUNCTIONS +//---------------------------------------------------------------------------// +/*! + * Create a Vector3D from a length-3 array. + */ +template +inline CELER_FUNCTION auto to_vgvector(Array const& a) + -> VgVector3 +{ + return {a[0], a[1], a[2]}; +} + +//---------------------------------------------------------------------------// +/*! + * Create a length-3 array from a VecGeom vector. + */ +template +inline CELER_FUNCTION auto to_array(vecgeom::Vector3D const& vgv) + -> Array +{ + return {vgv[0], vgv[1], vgv[2]}; +} + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/geocel/vg/detail/VecgeomCompatibility.hh b/src/geocel/vg/detail/VecgeomCompatibility.hh deleted file mode 100644 index 77dd04aec3..0000000000 --- a/src/geocel/vg/detail/VecgeomCompatibility.hh +++ /dev/null @@ -1,56 +0,0 @@ -//------------------------------ -*- C++ -*- -------------------------------// -// Copyright Celeritas contributors: see top-level COPYRIGHT file for details -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file geocel/vg/detail/VecgeomCompatibility.hh -//---------------------------------------------------------------------------// -#pragma once - -#include - -#include "corecel/Macros.hh" -#include "corecel/Types.hh" -#include "corecel/cont/Array.hh" -#include "corecel/cont/Span.hh" -#include "geocel/vg/VecgeomTypes.hh" - -namespace celeritas -{ -namespace detail -{ -//---------------------------------------------------------------------------// -/*! - * Create a Vector3D from a length-3 span. - */ -template -inline CELER_FUNCTION auto to_vector(Span s) - -> VgVector3, MemSpace::native> -{ - return {s[0], s[1], s[2]}; -} - -//---------------------------------------------------------------------------// -/*! - * Create a Vector3D from a length-3 array. - */ -template -inline CELER_FUNCTION auto to_vector(Array const& arr) - -> VgVector3 -{ - return to_vector(celeritas::make_span(arr)); -} - -//---------------------------------------------------------------------------// -/*! - * Create a length-3 array from a VecGeom vector. - */ -template -inline CELER_FUNCTION auto to_array(vecgeom::Vector3D const& vgv) - -> Array -{ - return {vgv[0], vgv[1], vgv[2]}; -} - -//---------------------------------------------------------------------------// -} // namespace detail -} // namespace celeritas diff --git a/src/orange/g4org/SolidConverter.cc b/src/orange/g4org/SolidConverter.cc index 84db185dc3..6ffe3ec315 100644 --- a/src/orange/g4org/SolidConverter.cc +++ b/src/orange/g4org/SolidConverter.cc @@ -56,6 +56,7 @@ #include "corecel/math/ArraySoftUnit.hh" #include "corecel/math/SoftEqual.hh" #include "corecel/sys/TypeDemangler.hh" +#include "geocel/g4/GeantTypes.hh" #include "orange/orangeinp/CsgObject.hh" #include "orange/orangeinp/IntersectRegion.hh" #include "orange/orangeinp/PolySolid.hh" @@ -195,7 +196,7 @@ EnclosedPolar enclosed_pol_from(S const& solid) -> std::pair { CELER_EXPECT(axis.z() > 0); - CELER_EXPECT(is_soft_unit_vector(convert_from_geant(axis))); + CELER_EXPECT(is_soft_unit_vector(to_array(axis))); return {native_value_to(std::acos(axis.z())), atan2turn(axis.y(), axis.x())}; @@ -413,8 +414,10 @@ auto SolidConverter::cuttubs(arg_type solid_base) -> result_type real_type const hh = scale_(solid.GetZHalfLength()); // Get bottom and top normal vectors - auto const b_norm = convert_from_geant(solid.GetLowNorm()); - auto const t_norm = convert_from_geant(solid.GetHighNorm()); + auto const b_norm + = static_array_cast(to_array(solid.GetLowNorm())); + auto const t_norm + = static_array_cast(to_array(solid.GetHighNorm())); // Optional inner cylinder std::optional inner; diff --git a/src/orange/g4org/Transformer.hh b/src/orange/g4org/Transformer.hh index f3f506bfb4..3e27e014d6 100644 --- a/src/orange/g4org/Transformer.hh +++ b/src/orange/g4org/Transformer.hh @@ -12,7 +12,9 @@ #include #include "geocel/Types.hh" -#include "orange/transform/VariantTransform.hh" +#include "geocel/g4/Convert.hh" +#include "orange/transform/Transformation.hh" +#include "orange/transform/Translation.hh" #include "Scaler.hh" @@ -52,14 +54,14 @@ class Transformer inline Translation operator()(G4ThreeVector const& t) const; // Convert a pure rotation - inline Transformation operator()(G4RotationMatrix const& rot) const; + inline Transformation operator()(G4RotationMatrix const& g4rm) const; // Convert a translation + rotation inline Transformation - operator()(G4ThreeVector const& t, G4RotationMatrix const& rot) const; + operator()(G4ThreeVector const& t, G4RotationMatrix const& g4rm) const; // Convert a more general transform (includes reflection) - inline Transformation operator()(G4Transform3D const& tran) const; + inline Transformation operator()(G4Transform3D const& g4tr) const; // Convert an affine transform inline Transformation operator()(G4AffineTransform const& at) const; @@ -70,24 +72,6 @@ class Transformer Scaler const& scale_; }; -//---------------------------------------------------------------------------// -// FREE FUNCTIONS -//---------------------------------------------------------------------------// -// Convert a ThreeVector -inline Real3 convert_from_geant(G4ThreeVector const& vec); - -//---------------------------------------------------------------------------// -// Convert three doubles to a Real3 -inline Real3 convert_from_geant(double x, double y, double z); - -//---------------------------------------------------------------------------// -// Convert a rotation matrix -inline SquareMatrixReal3 convert_from_geant(G4RotationMatrix const& rot); - -//---------------------------------------------------------------------------// -// Get the transpose/inverse of a rotation matrix -inline SquareMatrixReal3 transposed_from_geant(G4RotationMatrix const& rot); - //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// @@ -109,25 +93,29 @@ auto Transformer::operator()(G4ThreeVector const& t) const -> Translation /*! * Create a transform from a translation plus rotation. */ -auto Transformer::operator()(G4ThreeVector const& trans, - G4RotationMatrix const& rot) const +auto Transformer::operator()(G4ThreeVector const& g4t, + G4RotationMatrix const& g4rm) const -> Transformation { - return Transformation{convert_from_geant(rot), scale_.to(trans)}; + SquareMatrixReal3 mat{Real3(g4rm.xx(), g4rm.xy(), g4rm.xz()), + Real3(g4rm.yx(), g4rm.yy(), g4rm.yz()), + Real3(g4rm.zx(), g4rm.zy(), g4rm.zz())}; + + return Transformation{mat, scale_.to(g4t)}; } //---------------------------------------------------------------------------// /*! * Convert a more general transform (including possibly reflection). */ -Transformation Transformer::operator()(G4Transform3D const& tran) const +Transformation Transformer::operator()(G4Transform3D const& g4tr) const { - SquareMatrixReal3 rot{convert_from_geant(tran.xx(), tran.xy(), tran.xz()), - convert_from_geant(tran.yx(), tran.yy(), tran.yz()), - convert_from_geant(tran.zx(), tran.zy(), tran.zz())}; + SquareMatrixReal3 mat{Real3(g4tr.xx(), g4tr.xy(), g4tr.xz()), + Real3(g4tr.yx(), g4tr.yy(), g4tr.yz()), + Real3(g4tr.zx(), g4tr.zy(), g4tr.zz())}; - return Transformation{rot, - scale_.to(tran.dx(), tran.dy(), tran.dz())}; + return Transformation{mat, + scale_.to(g4tr.dx(), g4tr.dy(), g4tr.dz())}; } //---------------------------------------------------------------------------// @@ -139,50 +127,13 @@ Transformation Transformer::operator()(G4Transform3D const& tran) const auto Transformer::operator()(G4AffineTransform const& affine) const -> Transformation { - return Transformation{transposed_from_geant(affine.NetRotation()), - scale_.to(affine.NetTranslation())}; -} - -//---------------------------------------------------------------------------// -/*! - * Convert a ThreeVector. - */ -Real3 convert_from_geant(G4ThreeVector const& vec) -{ - return convert_from_geant(vec[0], vec[1], vec[2]); -} + // *Transpose* the rotation matrix + auto const& g4rm = affine.NetRotation(); + SquareMatrixReal3 mat{Real3(g4rm.xx(), g4rm.yx(), g4rm.zx()), + Real3(g4rm.xy(), g4rm.yy(), g4rm.zy()), + Real3(g4rm.xz(), g4rm.yz(), g4rm.zz())}; -//---------------------------------------------------------------------------// -/*! - * Convert three doubles to a Real3. - */ -Real3 convert_from_geant(double x, double y, double z) -{ - return Real3{{static_cast(x), - static_cast(y), - static_cast(z)}}; -} - -//---------------------------------------------------------------------------// -/*! - * Convert a rotation matrix. - */ -SquareMatrixReal3 convert_from_geant(G4RotationMatrix const& rot) -{ - return {convert_from_geant(rot.xx(), rot.xy(), rot.xz()), - convert_from_geant(rot.yx(), rot.yy(), rot.yz()), - convert_from_geant(rot.zx(), rot.zy(), rot.zz())}; -} - -//---------------------------------------------------------------------------// -/*! - * Get a transposed rotation matrix. - */ -SquareMatrixReal3 transposed_from_geant(G4RotationMatrix const& rot) -{ - return {convert_from_geant(rot.xx(), rot.yx(), rot.zx()), - convert_from_geant(rot.xy(), rot.yy(), rot.zy()), - convert_from_geant(rot.xz(), rot.yz(), rot.zz())}; + return Transformation{mat, scale_.to(affine.NetTranslation())}; } //---------------------------------------------------------------------------// diff --git a/test/accel/HepMC3PrimaryGenerator.test.cc b/test/accel/HepMC3PrimaryGenerator.test.cc index 203564b869..ae4face516 100644 --- a/test/accel/HepMC3PrimaryGenerator.test.cc +++ b/test/accel/HepMC3PrimaryGenerator.test.cc @@ -67,7 +67,7 @@ class HepMC3PrimaryGeneratorTest : public SimpleCmsTestBase result.energy.push_back(p->GetKineticEnergy()); result.mass.push_back(p->GetMass()); - auto dir = convert_from_geant(p->GetMomentumDirection(), 1); + auto dir = to_array(p->GetMomentumDirection()); result.dir.insert(result.dir.end(), dir.begin(), dir.end()); } } diff --git a/test/accel/IntegrationTestBase.cc b/test/accel/IntegrationTestBase.cc index cca0ae1a9c..669a0a1440 100644 --- a/test/accel/IntegrationTestBase.cc +++ b/test/accel/IntegrationTestBase.cc @@ -432,8 +432,8 @@ auto LarSphereIntegrationMixin::make_primary_input() const -> PrimaryInput PrimaryInput result; result.pdg = {pdg::electron()}; result.energy = inp::MonoenergeticDistribution{10}; // [MeV] - result.shape - = inp::PointDistribution{array_cast(from_cm({99, 0.1, 0}))}; + result.shape = inp::PointDistribution{ + static_array_cast(from_cm({99, 0.1, 0}))}; result.angle = inp::IsotropicDistribution{}; result.num_events = 4; // Overridden with BeamOn result.primaries_per_event = 10; @@ -503,8 +503,8 @@ auto TestEm3IntegrationMixin::make_primary_input() const -> PrimaryInput PrimaryInput result; result.pdg = {pdg::electron()}; result.energy = inp::MonoenergeticDistribution{100}; // [MeV] - result.shape - = inp::PointDistribution{array_cast(from_cm({-22, 0, 0}))}; + result.shape = inp::PointDistribution{ + static_array_cast(from_cm({-22, 0, 0}))}; result.angle = inp::MonodirectionalDistribution{{1, 0, 0}}; result.num_events = 2; result.primaries_per_event = 1; @@ -552,8 +552,8 @@ auto OpNoviceIntegrationMixin::make_primary_input() const -> PrimaryInput PrimaryInput result; result.pdg = {pdg::positron()}; result.energy = inp::MonoenergeticDistribution{0.5}; // [MeV] - result.shape - = inp::PointDistribution{array_cast(from_cm({0., 0., 0.}))}; + result.shape = inp::PointDistribution{ + static_array_cast(from_cm({0., 0., 0.}))}; result.angle = inp::MonodirectionalDistribution{{1., 0., 0.}}; result.num_events = 12; // Overridden with BeamOn result.primaries_per_event = 10; diff --git a/test/accel/TrackingManagerIntegration.test.cc b/test/accel/TrackingManagerIntegration.test.cc index 7558d5911b..e6d0828c3b 100644 --- a/test/accel/TrackingManagerIntegration.test.cc +++ b/test/accel/TrackingManagerIntegration.test.cc @@ -346,7 +346,7 @@ class LarSphereOptical : public LarSphere auto result = LarSphereIntegrationMixin::make_primary_input(); result.shape = inp::PointDistribution{ - array_cast(from_cm({0.1, 0.1, 0}))}; + static_array_cast(from_cm({0.1, 0.1, 0}))}; result.primaries_per_event = 1; result.energy = inp::MonoenergeticDistribution{2}; // [MeV] return result; @@ -625,8 +625,8 @@ auto OpticalSurfaces::make_primary_input() const -> PrimaryInput { PrimaryInput result; result.pdg = {pdg::positron()}; - result.shape - = inp::PointDistribution{array_cast(from_cm({30, 0, 0}))}; + result.shape = inp::PointDistribution{ + static_array_cast(from_cm({30, 0, 0}))}; result.angle = inp::MonodirectionalDistribution{{-1, 0, 0}}; result.energy = inp::MonoenergeticDistribution{100}; // [MeV] result.primaries_per_event = 1; diff --git a/test/accel/UserActionIntegration.test.cc b/test/accel/UserActionIntegration.test.cc index da230dc1c4..af909a568f 100644 --- a/test/accel/UserActionIntegration.test.cc +++ b/test/accel/UserActionIntegration.test.cc @@ -139,8 +139,8 @@ auto LarSphereOpticalOffload::make_primary_input() const -> PrimaryInput { auto result = LarSphereIntegrationMixin::make_primary_input(); - result.shape - = inp::PointDistribution{array_cast(from_cm({0.1, 0.1, 0}))}; + result.shape = inp::PointDistribution{ + static_array_cast(from_cm({0.1, 0.1, 0}))}; result.primaries_per_event = 1; result.energy = inp::MonoenergeticDistribution{1}; // [MeV] return result; diff --git a/test/celeritas/Units.test.cc b/test/celeritas/Units.test.cc index 9a61e75890..be806580b7 100644 --- a/test/celeritas/Units.test.cc +++ b/test/celeritas/Units.test.cc @@ -110,10 +110,10 @@ TEST(UnitsTest, clhep) } double g4_native = 2.5 * CLHEP::tesla; - auto celer_native = convert_from_geant(g4_native, clhep_field); + auto celer_native = convert_from_geant(g4_native); EXPECT_SOFT_EQ(2.5 * units::tesla, celer_native); EXPECT_SOFT_EQ(2.5, native_value_to(celer_native).value()); - EXPECT_SOFT_EQ(g4_native, convert_to_geant(celer_native, clhep_field)); + EXPECT_SOFT_EQ(g4_native, convert_to_geant(celer_native)); } { double g4_native = 1.5 * CLHEP::s; diff --git a/test/corecel/cont/Array.test.cc b/test/corecel/cont/Array.test.cc index 7e1c670d51..c15048ca1f 100644 --- a/test/corecel/cont/Array.test.cc +++ b/test/corecel/cont/Array.test.cc @@ -120,6 +120,27 @@ TEST(ArrayTest, two_level) EXPECT_VEC_EQ((Int3{7, 8, 9}), x[2]); } +TEST(ArrayTest, to_array) +{ + long vals[] = {1, 2, 3}; + auto arr = to_array(vals); + EXPECT_TRUE((std::is_same_v, decltype(arr)>)); + EXPECT_VEC_EQ((Array{1, 2, 3}), arr); +} + +TEST(ArrayTest, casts) +{ + // Test up- and down-casting + Array dbl = {1.4, 3.1}; + auto flt = static_array_cast(dbl); + EXPECT_TRUE((std::is_same_v, decltype(flt)>)); + EXPECT_VEC_EQ((Array{1.4f, 3.1f}), flt); + + flt[1] = 2.3f; + dbl = static_array_cast(flt); + EXPECT_DOUBLE_EQ(static_cast(2.3f), dbl[1]); +} + TEST(EnumArrayTest, all) { EnumArray x = {1, 3, 2}; diff --git a/test/corecel/math/Quantity.test.cc b/test/corecel/math/Quantity.test.cc index 2660265b79..9de4897b03 100644 --- a/test/corecel/math/Quantity.test.cc +++ b/test/corecel/math/Quantity.test.cc @@ -10,6 +10,7 @@ #include "corecel/Config.hh" +#include "corecel/math/ArrayQuantity.hh" #include "corecel/math/QuantityIO.json.hh" #include "corecel/math/Turn.hh" @@ -229,6 +230,27 @@ TEST(QuantityTest, swappiness) EXPECT_EQ(144, native_value_from(gross)); } +TEST(QuantityTest, array) +{ + using Dozen3 = Array; + using Int3 = Array; + + auto dozens = make_quantity_array(1, 0, 2); + EXPECT_TRUE((std::is_same_v)); + EXPECT_EQ(Dozen{2}, dozens[2]); + auto d = native_value_from(dozens); + EXPECT_TRUE((std::is_same_v)); + EXPECT_EQ(Int3(12, 0, 24), d); + + auto dozens2 = make_quantity_array(Int3{1, 2, 3}); + EXPECT_TRUE((std::is_same_v)); + EXPECT_EQ(Int3(12, 24, 36), native_value_from(dozens2)); + + auto dva = value_as(dozens2); + EXPECT_TRUE((std::is_same_v)); + EXPECT_EQ(Int3(1, 2, 3), dva); +} + TEST(QuantityTest, io) { { diff --git a/test/orange/g4org/Transformer.test.cc b/test/orange/g4org/Transformer.test.cc index 278715c274..3bc2ac1e4a 100644 --- a/test/orange/g4org/Transformer.test.cc +++ b/test/orange/g4org/Transformer.test.cc @@ -31,12 +31,12 @@ constexpr real_type mm{::celeritas::lengthunits::millimeter}; Real3 from_geant(G4ThreeVector const& tv) { - return convert_from_geant(tv); + return static_array_cast(to_array(tv)); } G4ThreeVector to_geant(Real3 const& rv) { - return {rv[0], rv[1], rv[2]}; + return to_g4vector(static_array_cast(rv)); } class TransformerTest : public ::celeritas::test::Test @@ -68,13 +68,10 @@ TEST_F(TransformerTest, affine_transform) // Construct Celeritas matrix auto const mat = make_rotation(rot_axis, rot_turn); - { - auto actual = convert_from_geant(g4mat); - // Check raw matrix conversion - EXPECT_VEC_SOFT_EQ(mat[0], actual[0]); - EXPECT_VEC_SOFT_EQ(mat[1], actual[1]); - EXPECT_VEC_SOFT_EQ(mat[2], actual[2]); - } + // Check raw matrix conversion + EXPECT_VEC_SOFT_EQ(mat[0], Real3(g4mat.xx(), g4mat.xy(), g4mat.xz())); + EXPECT_VEC_SOFT_EQ(mat[1], Real3(g4mat.yx(), g4mat.yy(), g4mat.yz())); + EXPECT_VEC_SOFT_EQ(mat[2], Real3(g4mat.zx(), g4mat.zy(), g4mat.zz())); EXPECT_VEC_SOFT_EQ(gemv(mat, Real3{1, 0, 0}), from_geant(g4mat(G4ThreeVector(1, 0, 0))));