Skip to content

Add helpers for array type, precision, and quantity conversion #2223

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sethrj wants to merge 22 commits into
base: develop
Choose a base branch
from
Open

Add helpers for array type, precision, and quantity conversion #2223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
a155cee
Revert "REVERTME future: array quantity"
sethrj Jan 30, 2026
95df594
Rename static_array_cast and move to Array
sethrj Jan 30, 2026
0a87eb9
Add to_array for C array-to-Array conversion
sethrj Jan 28, 2026
d38c08e
Use to_vgvector for vecgeom
sethrj Jan 28, 2026
a34fb8b
Use constexpr for array operators
sethrj Jan 28, 2026
3541a25
Fix g4app SD units
sethrj Jan 28, 2026
05f5d60
Add value_as for array quanity
sethrj Jan 31, 2026
7ea63aa
Refactor geant conversion to array types vs precision
sethrj Jan 31, 2026
4009f79
Update ORANGE g4 conversion
sethrj Jan 31, 2026
472a1d0
Merge remote-tracking branch 'upstream/develop' into quantity-array
sethrj Feb 1, 2026
8c89d6c
Upgrade some uses of convert_to_geant
sethrj Feb 1, 2026
262eee9
use LV utility
sethrj Feb 1, 2026
f3aebe6
Restore clhep_length for compatibility
sethrj Feb 1, 2026
b2337f9
fixup! Upgrade some uses of convert_to_geant
sethrj Feb 1, 2026
559f024
Add quantities and fix single precisoin build
sethrj Feb 1, 2026
e3b47c5
Update casts
sethrj Feb 1, 2026
bc87762
Remove unnecessary template arg
sethrj Feb 1, 2026
3f4a067
formatting
sethrj Feb 1, 2026
7291409
Merge remote-tracking branch 'upstream/develop' into quantity-array
sethrj Feb 1, 2026
90f3db3
Fix docs
sethrj Feb 2, 2026
ac22663
Restore deprecated function signature
sethrj Feb 2, 2026
45a4fec
Address feedback
sethrj Feb 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions app/celer-g4/SensitiveDetector.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@
#include "SensitiveDetector.hh"

#include <memory>
#include <CLHEP/Units/SystemOfUnits.h>
#include <G4HCofThisEvent.hh>
#include <G4SDManager.hh>
#include <G4Step.hh>
Expand All @@ -17,6 +16,7 @@

#include "corecel/Assert.hh"
#include "geocel/g4/Convert.hh"
#include "celeritas/ext/GeantUnits.hh"

#include "GlobalSetup.hh"

Expand Down Expand Up @@ -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()});
Comment on lines +92 to +93
Copy link
Contributor

@amandalund amandalund Feb 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe the former should also be native_value_from (just to make it clearer to developers that it's converting units) even though we know the energy units are both MeV.


collection_->insert(new SensitiveHit(hit));
return true;
Expand Down
3 changes: 3 additions & 0 deletions doc/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 \""
Expand Down
7 changes: 6 additions & 1 deletion doc/implementation/units-constants.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<UnitT, ValueT> quant) noexcept
.. doxygenfunction:: celeritas::native_value_from
.. doxygenfunction:: celeritas::value_as

.. doxygenfunction:: celeritas::zero_quantity
Expand All @@ -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
Expand Down
25 changes: 14 additions & 11 deletions src/accel/CartMapMagneticField.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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<units::Millimeter, real_type>;
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
Expand All @@ -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;
Expand All @@ -79,7 +80,9 @@ MakeCartMapFieldInput(G4Field const& field,
auto field_converter = [](Array<G4double, 3> const& bfield,
Array<G4double, 4> const&,
real_type cur_bfield[3]) {
auto bfield_native = convert_from_geant(bfield.data(), clhep_field);
using ClhepField = Quantity<units::ClhepUnitBField, double>;
Copy link
Contributor

@amandalund amandalund Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use ClhepField from Quantities.hh?

Copy link
Member Author

@sethrj sethrj Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not the same: real_type vs Geant4 value (double ).

Copy link
Contributor

@amandalund amandalund Feb 5, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The ones in Quantities are double too, were they meant to be real_type?

auto bfield_native
= native_value_from(make_quantity_array<ClhepField>(bfield));
std::copy(bfield_native.cbegin(), bfield_native.cend(), cur_bfield);
};

Expand Down
15 changes: 8 additions & 7 deletions src/accel/CylMapMagneticField.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<units::Millimeter, double>;
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(),
Expand All @@ -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();
Expand All @@ -104,13 +105,13 @@ MakeCylMapFieldInput(G4Field const& field,
auto field_converter = [](Array<G4double, 3> const& bfield,
Array<G4double, 4> const& pos,
real_type cur_bfield[3]) {
using ClhepField = Quantity<units::ClhepUnitBField, G4double>;
auto bfield_native
= native_value_from(make_quantity_array<ClhepField>(bfield));
double const phi = std::atan2(pos[1], pos[0]);
EnumArray<CylAxis, G4double> 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
Expand Down
4 changes: 2 additions & 2 deletions src/accel/HepMC3PrimaryGenerator.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,8 +62,8 @@ class PrimaryInserter
G4PrimaryParticle* primary = new G4PrimaryParticle(par.pid());

// Set the primary directlon
auto dir = make_unit_vector(Array<double, 3>{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());
Expand Down
3 changes: 2 additions & 1 deletion src/accel/LocalTransporter.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<real_type>(
to_array(g4track.GetMomentumDirection()));
track.time = convert_from_geant(g4track.GetGlobalTime(), clhep_time);
track.weight = g4track.GetWeight();

Expand Down
6 changes: 4 additions & 2 deletions src/accel/PGPrimaryGeneratorAction.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<lengthunits::ClhepLength>(p.position));
gun_.SetParticleMomentumDirection(
to_g4vector(static_array_cast<double>(p.direction)));
gun_.SetParticleEnergy(convert_to_geant(p.energy, CLHEP::MeV));
gun_.GeneratePrimaryVertex(event);

Expand Down
9 changes: 9 additions & 0 deletions src/celeritas/Quantities.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,15 @@ using AmuMass = RealQuantity<Amu>;
//! Special faux quantity for overloading cross section calculation
using LogMevEnergy = RealQuantity<LogMev>;

//---------------------------------------------------------------------------//
//!@{
//! \name Quantities for Geant4 conversion
using ClhepLength = Quantity<Millimeter, double>;
using ClhepField = Quantity<ClhepUnitBField, double>;
using ClhepTime = Quantity<Nanosecond, double>;
using ClhepEnergy = Quantity<Mev, double>;
//!@}

//---------------------------------------------------------------------------//
//!@{
//! \name Quantities for manual input and/or test harnesses
Expand Down
2 changes: 1 addition & 1 deletion src/celeritas/ext/GeantUnits.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()};
Expand Down
3 changes: 1 addition & 2 deletions src/celeritas/ext/detail/HitProcessor.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<void const*>(lv));
<< StreamableLV{lv});
}

CELER_ENSURE(!detectors_.empty());
Expand Down
5 changes: 3 additions & 2 deletions src/celeritas/ext/detail/NaviTouchableUpdater.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<units::ClhepLength>(static_array_cast<double>(pos));
auto g4dir = to_g4vector(static_array_cast<double>(dir));

// Locate pre-step point
navi_->LocateGlobalPointAndUpdateTouchable(g4pos,
Expand Down
7 changes: 3 additions & 4 deletions src/celeritas/field/UniformFieldParams.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -69,10 +70,8 @@ validated_field_data(UniformFieldParams::Input const& inp)
validate_input(inp.driver_options);

HostVal<UniformFieldParamsData> 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<units::FieldTesla>(inp.strength));
result.options = inp.driver_options;
return result;
}
Expand Down
33 changes: 17 additions & 16 deletions src/celeritas/phys/PrimaryGenerator.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,7 @@ auto make_energy_sampler(inp::EnergyDistribution const& i)
return_as<PrimaryGenerator::EnergySampler>(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<real_type>(me.value)};
},
Expand All @@ -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<PrimaryGenerator::PositionSampler>(Overload{
[](inp::PointDistribution const& ps) {
return DeltaDistribution{array_cast<real_type>(ps.value)};
},
[](inp::UniformBoxDistribution const& ubs) {
return UniformBoxDistribution{array_cast<real_type>(ubs.lower),
array_cast<real_type>(ubs.upper)};
}}),
i);
return std::visit(return_as<PrimaryGenerator::PositionSampler>(Overload{
[](inp::PointDistribution const& ps) {
return DeltaDistribution{
static_array_cast<real_type>(ps.value)};
},
[](inp::UniformBoxDistribution const& ubs) {
return UniformBoxDistribution{
static_array_cast<real_type>(ubs.lower),
static_array_cast<real_type>(ubs.upper)};
}}),
i);
}

//---------------------------------------------------------------------------//
Expand All @@ -81,10 +81,11 @@ auto make_direction_sampler(inp::AngleDistribution const& i)
return IsotropicDistribution<real_type>{};
},
[](inp::MonodirectionalDistribution const& ma) {
CELER_VALIDATE(is_soft_unit_vector(ma.value),
<< "primary generator angle is "
"not a unit vector");
return DeltaDistribution{array_cast<real_type>(ma.value)};
CELER_VALIDATE(
is_soft_unit_vector(ma.value),
<< R"(primary generator angle is not a unit vector)");
return DeltaDistribution{
static_array_cast<real_type>(ma.value)};
}}),
i);
}
Expand Down
39 changes: 39 additions & 0 deletions src/corecel/cont/Array.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -179,6 +179,45 @@ CELER_CEF bool operator!=(Array<T, N> const& lhs, Array<T, N> 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<class T, size_type N>
CELER_CONSTEXPR_FUNCTION auto to_array(T (&x)[N])
{
static_assert(!std::is_array_v<T>,
"to_array elements cannot be multidimensional");
static_assert(std::is_constructible_v<T, T&>,
"to_array elements must be copy constructible");

Array<std::remove_cv_t<T>, 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<class T1, class T2, size_type N>
CELER_CONSTEXPR_FUNCTION Array<T1, N> static_array_cast(Array<T2, N> const& x)
{
Array<T1, N> result;
for (size_type i = 0; i != N; ++i)
{
result[i] = static_cast<T1>(x[i]);
}
return result;
}

//---------------------------------------------------------------------------//
} // namespace celeritas

Expand Down
Loading
Loading