diff --git a/doc/implementation/geant4-interface/low-level.rst b/doc/implementation/geant4-interface/low-level.rst index 6106a78f8e..3bfc6694a3 100644 --- a/doc/implementation/geant4-interface/low-level.rst +++ b/doc/implementation/geant4-interface/low-level.rst @@ -36,3 +36,13 @@ Utility interfaces .. doxygenclass:: celeritas::ScopedGeantLogger .. doxygenclass:: celeritas::ScopedGeantExceptionHandler + +.. _g4_views: + +Views +^^^^^ + +.. doxygenclass:: celeritas::GeantParticleView +.. doxygenclass:: celeritas::GeantStepPointView +.. doxygenclass:: celeritas::GeantTrackView +.. doxygenclass:: celeritas::GeantStepView diff --git a/src/accel/LocalTransporter.cc b/src/accel/LocalTransporter.cc index 8248fc4f49..e202646f40 100644 --- a/src/accel/LocalTransporter.cc +++ b/src/accel/LocalTransporter.cc @@ -25,6 +25,7 @@ #include "corecel/cont/Span.hh" #include "corecel/io/BuildOutput.hh" #include "corecel/io/Logger.hh" +#include "corecel/math/QuantityIO.hh" #include "corecel/sys/Device.hh" #include "corecel/sys/ScopedProfiling.hh" #include "corecel/sys/ScopedSignalHandler.hh" @@ -35,6 +36,7 @@ #include "celeritas/Quantities.hh" #include "celeritas/Types.hh" #include "celeritas/ext/GeantSd.hh" +#include "celeritas/ext/GeantTrackView.hh" #include "celeritas/ext/GeantUnits.hh" #include "celeritas/ext/detail/HitProcessor.hh" #include "celeritas/global/ActionSequence.hh" @@ -245,56 +247,49 @@ void LocalTransporter::Push(G4Track& g4track) ScopedProfiling profile_this{"push"}; - if (Real3 pos = convert_from_geant(g4track.GetPosition(), 1); - !is_inside(bbox_, pos)) + GeantTrackView gtv{g4track}; + + if (!is_inside(bbox_, gtv.pos())) { // Primary may have been created by a particle generator outside the // geometry - double energy - = convert_from_geant(g4track.GetKineticEnergy(), CLHEP::MeV); CELER_LOG_LOCAL(error) - << "Discarding track outside world bounds: " << energy - << " MeV from " << g4track.GetDefinition()->GetParticleName() - << " at " << pos << " along " - << convert_from_geant(g4track.GetMomentumDirection(), 1); + << "Discarding track outside world bounds: " << gtv.energy() + << " from " << gtv.particle().name() << " at " << gtv.pos() + << " along " << gtv.dir(); - buffer_accum_.lost_energy += energy; + buffer_accum_.lost_energy += gtv.energy().value(); ++buffer_accum_.lost_primaries; return; } - Primary track; - - PDGNumber const pdg{g4track.GetDefinition()->GetPDGEncoding()}; - track.particle_id = particles_->find(pdg); + Primary offloaded; // Generate Celeritas-specific PrimaryID if (hit_processor_) { - track.primary_id + offloaded.primary_id = hit_processor_->track_processor().register_primary(g4track); } - track.energy = units::MevEnergy( - convert_from_geant(g4track.GetKineticEnergy(), CLHEP::MeV)); + offloaded.energy = gtv.energy(); + offloaded.particle_id = particles_->find(gtv.particle().pdg()); + offloaded.position = gtv.pos(); + offloaded.direction = gtv.dir(); + offloaded.time = gtv.time(); + offloaded.weight = gtv.weight(); - CELER_VALIDATE(track.particle_id, - << "cannot offload '" - << g4track.GetDefinition()->GetParticleName() + CELER_VALIDATE(offloaded.particle_id, + << "cannot offload '" << gtv.particle().name() << "' particles"); - 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_time); - track.weight = g4track.GetWeight(); - /*! * \todo Eliminate event ID from primary. */ - track.event_id = EventId{0}; + offloaded.event_id = EventId{0}; - buffer_.push_back(track); - buffer_accum_.energy += track.energy.value(); + buffer_.push_back(offloaded); + buffer_accum_.energy += offloaded.energy.value(); if (buffer_.size() >= auto_flush_) { this->Flush(); diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index 555e15cc39..09677cd0ad 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -224,6 +224,8 @@ if(CELERITAS_USE_Geant4) ext/GeantSd.cc ext/GeantSdOutput.cc ext/GeantSetup.cc + ext/GeantStepPointView.cc + ext/GeantStepView.cc ext/detail/EmStandardPhysics.cc ext/detail/GeantMicroXsCalculator.cc ext/detail/GeantModelImporter.cc diff --git a/src/celeritas/ext/GeantParticleView.hh b/src/celeritas/ext/GeantParticleView.hh index bc46588a20..f6c2e015f8 100644 --- a/src/celeritas/ext/GeantParticleView.hh +++ b/src/celeritas/ext/GeantParticleView.hh @@ -6,13 +6,17 @@ //---------------------------------------------------------------------------// #pragma once +#include #include #include #include "corecel/math/Quantity.hh" +#include "geocel/g4/Convert.hh" #include "celeritas/UnitTypes.hh" #include "celeritas/phys/PDGNumber.hh" +#include "GeantUnits.hh" + namespace celeritas { //---------------------------------------------------------------------------// @@ -72,11 +76,9 @@ auto GeantParticleView::decay_constant() const -> real_type return 0; } - // CLHEP time unit system - using Time = Quantity; - - // Decay constant is 1/lifetime - return 1 / native_value_from(Time{pd_.GetPDGLifeTime()}); + // Decay constant is 1/lifetime (lifetime is in CLHEP time units) + real_type lifetime = convert_from_geant(pd_.GetPDGLifeTime(), clhep_time); + return 1 / lifetime; } //---------------------------------------------------------------------------// diff --git a/src/celeritas/ext/GeantStepPointView.cc b/src/celeritas/ext/GeantStepPointView.cc new file mode 100644 index 0000000000..879056f758 --- /dev/null +++ b/src/celeritas/ext/GeantStepPointView.cc @@ -0,0 +1,86 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/ext/GeantStepPointView.cc +//---------------------------------------------------------------------------// +#include "GeantStepPointView.hh" + +#include + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Update attributes from logical volume. + */ +void GeantStepPointView::update_from_volume(G4LogicalVolume const& lv) +{ + sp_.SetMaterial(lv.GetMaterial()); + sp_.SetMaterialCutsCouple(lv.GetMaterialCutsCouple()); + sp_.SetSensitiveDetector(lv.GetSensitiveDetector()); +} + +//---------------------------------------------------------------------------// +/*! + * Update attributes from the touchable's logical volume if possible. + * + * If the step point has an associated touchable, and that touchable is inside + * the geometry, it updates. Otherwise, it clears the corresponding attributes. + */ +void GeantStepPointView::update_from_volume() +{ + G4LogicalVolume const* lv = nullptr; + if (auto* touch = sp_.GetTouchable()) + { + // The physical volume could be null if post-step is outside + if (auto* pv = touch->GetVolume()) + { + lv = pv->GetLogicalVolume(); + } + } + if (lv) + { + this->update_from_volume(*lv); + } + else + { + sp_.SetMaterial(nullptr); + sp_.SetMaterialCutsCouple(nullptr); + sp_.SetSensitiveDetector(nullptr); + } +} + +//---------------------------------------------------------------------------// +/*! + * Update mass and charge from particle definition. + */ +void GeantStepPointView::update_from_particle(GeantParticleView const& particle) +{ + sp_.SetMass(particle.mass().value() * CLHEP::MeV); + sp_.SetCharge(particle.charge().value() * CLHEP::eplus); +} + +//---------------------------------------------------------------------------// +/*! + * Clear unsupported attributes to invalid sentinel values. + * + * This sets attributes that Celeritas does not currently track to sentinel + * values to indicate they are unavailable to sensitive detectors. + */ +void GeantStepPointView::clear_unsupported() +{ + // Time since track was created + sp_.SetLocalTime(std::numeric_limits::infinity()); + // Time in rest frame since track was created + sp_.SetProperTime(std::numeric_limits::infinity()); + // Speed (TODO: use ParticleView) + sp_.SetVelocity(std::numeric_limits::infinity()); + // Safety distance + sp_.SetSafety(std::numeric_limits::infinity()); + // Polarization (default to zero) + sp_.SetPolarization(G4ThreeVector()); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/ext/GeantStepPointView.hh b/src/celeritas/ext/GeantStepPointView.hh new file mode 100644 index 0000000000..f6d54f9593 --- /dev/null +++ b/src/celeritas/ext/GeantStepPointView.hh @@ -0,0 +1,179 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/ext/GeantStepPointView.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/Assert.hh" +#include "corecel/Types.hh" +#include "corecel/math/Quantity.hh" +#include "geocel/g4/Convert.hh" +#include "celeritas/UnitTypes.hh" + +#include "GeantParticleView.hh" +#include "GeantUnits.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Access and modify step point data from Geant4 with Celeritas units. + * + * This provides a uniform interface to G4StepPoint data using Celeritas types + * and units. Geant4 data are all in double precision. + */ +class GeantStepPointView +{ + public: + //!@{ + //! \name Type aliases + using Energy = Quantity; + using real_type = double; + //!@} + + public: + // Construct from G4StepPoint + explicit GeantStepPointView(G4StepPoint& step_point) : sp_(step_point) {} + + //!@{ + //! \name Accessors + + // Position in native Celeritas length units + inline Real3 pos() const; + + // Momentum direction (unit vector) + inline Real3 dir() const; + + // Kinetic energy [MeV] + inline Energy energy() const; + + // Global time in native Celeritas time units + inline real_type time() const; + + //! Statistical weight + real_type weight() const { return sp_.GetWeight(); } + + //!@} + //!@{ + //! \name Mutators + + // Set position in native Celeritas length units + inline void pos(Real3 const& position); + + // Set momentum direction (unit vector) + inline void dir(Real3 const& direction); + + // Set kinetic energy [MeV] + inline void energy(Energy kinetic_energy); + + // Set global time in native Celeritas time units + inline void time(real_type global_time); + + // Set statistical weight + void weight(real_type w) { sp_.SetWeight(w); } + + // Update attributes from logical volume + void update_from_volume(G4LogicalVolume const& lv); + + // Update attributes from touchable's logical volume + void update_from_volume(); + + // Update mass and charge from particle definition + void update_from_particle(GeantParticleView const& particle); + + // Clear unsupported attributes to invalid sentinel values + void clear_unsupported(); + + //!@} + + //! Access underlying G4 object + G4StepPoint& step_point() { return sp_; } + + private: + G4StepPoint& sp_; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Get position in native Celeritas length units. + */ +Real3 GeantStepPointView::pos() const +{ + return convert_from_geant(sp_.GetPosition(), clhep_length); +} + +//---------------------------------------------------------------------------// +/*! + * Get momentum direction. + */ +Real3 GeantStepPointView::dir() const +{ + return convert_from_geant(sp_.GetMomentumDirection(), 1); +} + +//---------------------------------------------------------------------------// +/*! + * Get kinetic energy in MeV. + */ +GeantStepPointView::Energy GeantStepPointView::energy() const +{ + return Energy{convert_from_geant(sp_.GetKineticEnergy(), CLHEP::MeV)}; +} + +//---------------------------------------------------------------------------// +/*! + * Get global time in native Celeritas time units. + */ +real_type GeantStepPointView::time() const +{ + return convert_from_geant(sp_.GetGlobalTime(), clhep_time); +} + +//---------------------------------------------------------------------------// +/*! + * Set position in native Celeritas length units. + */ +void GeantStepPointView::pos(Real3 const& position) +{ + sp_.SetPosition(convert_to_geant(position, clhep_length)); +} + +//---------------------------------------------------------------------------// +/*! + * Set momentum direction. + */ +void GeantStepPointView::dir(Real3 const& direction) +{ + sp_.SetMomentumDirection(convert_to_geant(direction, 1)); +} + +//---------------------------------------------------------------------------// +/*! + * Set kinetic energy in MeV. + */ +void GeantStepPointView::energy(Energy kinetic_energy) +{ + CELER_EXPECT(kinetic_energy >= zero_quantity()); + sp_.SetKineticEnergy(convert_to_geant(kinetic_energy.value(), CLHEP::MeV)); + // TODO: update speed based on mass, KE +} + +//---------------------------------------------------------------------------// +/*! + * Set global time in native Celeritas time units. + */ +void GeantStepPointView::time(real_type global_time) +{ + CELER_EXPECT(global_time >= 0); + sp_.SetGlobalTime(convert_to_geant(global_time, clhep_time)); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/ext/GeantStepView.cc b/src/celeritas/ext/GeantStepView.cc new file mode 100644 index 0000000000..e1889ce32e --- /dev/null +++ b/src/celeritas/ext/GeantStepView.cc @@ -0,0 +1,80 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/ext/GeantStepView.cc +//---------------------------------------------------------------------------// +#include "GeantStepView.hh" + +#include "GeantTrackView.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Update track from step data. + * + * Copies step length and step point data to the track. This is similar to + * \c G4Step::UpdateTrack but applies only to attributes we know about and + * safely handles null pointers. + */ +void GeantStepView::update_track() +{ + CELER_EXPECT(s_.GetTrack()); + + GeantTrackView track{*s_.GetTrack()}; + GeantParticleView particle_view = track.particle(); + + // Update pre-step point if present + if (G4StepPoint* pre_step = s_.GetPreStepPoint()) + { + GeantStepPointView{*pre_step}.update_from_particle(particle_view); + track.track().SetTouchableHandle(pre_step->GetTouchableHandle()); + } + + // Update post-step point and track from post-step if present + if (G4StepPoint* post_step = s_.GetPostStepPoint()) + { + GeantStepPointView post_view{*post_step}; + post_view.update_from_particle(particle_view); + + // Copy post-step state to track + track.time(post_view.time()); + track.pos(post_view.pos()); + track.energy(post_view.energy()); + track.dir(post_view.dir()); + track.weight(post_view.weight()); + + track.track().SetNextTouchableHandle(post_step->GetTouchableHandle()); + track.track().SetVelocity(post_step->GetVelocity()); + } +} + +//---------------------------------------------------------------------------// +/*! + * Delete a step point. + * + * This sets the step point to null if supported by the Geant4 version, + * otherwise does nothing (no reset available before v11.0.3). + */ +void GeantStepView::delete_step_point(StepPoint sp) +{ + CELER_EXPECT(sp != StepPoint::size_); + +#if G4VERSION_NUMBER >= 1103 + if (sp == StepPoint::pre) + { + s_.ResetPreStepPoint(nullptr); + } + else + { + s_.ResetPostStepPoint(nullptr); + } +#else + // No "reset" available before v11.0.3 + CELER_DISCARD(sp); +#endif +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/ext/GeantStepView.hh b/src/celeritas/ext/GeantStepView.hh new file mode 100644 index 0000000000..bda885e3a1 --- /dev/null +++ b/src/celeritas/ext/GeantStepView.hh @@ -0,0 +1,170 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/ext/GeantStepView.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/Assert.hh" +#include "corecel/math/Quantity.hh" +#include "geocel/g4/Convert.hh" +#include "celeritas/Types.hh" +#include "celeritas/UnitTypes.hh" + +#include "GeantStepPointView.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Access and modify step data from Geant4 with Celeritas units. + * + * This provides a uniform interface to G4Step data using Celeritas types and + * units. Geant4 data are all in double precision. + */ +class GeantStepView +{ + public: + //!@{ + //! \name Type aliases + using Energy = Quantity; + using real_type = double; + //!@} + + public: + // Construct from G4Step + explicit GeantStepView(G4Step& step) : s_(step) {} + + //!@{ + //! \name Accessors + + // Total energy deposited during step [MeV] + inline Energy energy_deposition() const; + + // Step length in native Celeritas length units + inline real_type step_length() const; + + // Pre-step point accessor + inline GeantStepPointView pre_step() const; + + // Post-step point accessor + inline GeantStepPointView post_step() const; + + // Whether the step point is valid (not null) + inline bool has_step_point(StepPoint sp) const; + + // Step point accessor by enum + inline GeantStepPointView step_point(StepPoint sp) const; + + //!@} + //!@{ + //! \name Mutators + + // Set total energy deposited during step [MeV] + inline void energy_deposition(Energy edep); + + // Set step length in native Celeritas length units + inline void step_length(real_type length); + + // Update track from step data + void update_track(); + + // Delete a step point (set to null or clear based on Geant4 version) + void delete_step_point(StepPoint sp); + + //!@} + + private: + G4Step& s_; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Get total energy deposited during step in MeV. + */ +auto GeantStepView::energy_deposition() const -> Energy +{ + return Energy{convert_from_geant(s_.GetTotalEnergyDeposit(), CLHEP::MeV)}; +} + +//---------------------------------------------------------------------------// +/*! + * Get step length in native Celeritas length units. + */ +real_type GeantStepView::step_length() const +{ + return convert_from_geant(s_.GetStepLength(), clhep_length); +} + +//---------------------------------------------------------------------------// +/*! + * Set total energy deposited during step in MeV. + */ +void GeantStepView::energy_deposition(Energy edep) +{ + s_.SetTotalEnergyDeposit(convert_to_geant(edep.value(), CLHEP::MeV)); +} + +//---------------------------------------------------------------------------// +/*! + * Set step length in native Celeritas length units. + */ +void GeantStepView::step_length(real_type length) +{ + s_.SetStepLength(convert_to_geant(length, clhep_length)); + if (s_.GetTrack()) + { + // Set on track as well + s_.GetTrack()->SetStepLength(s_.GetStepLength()); + } +} + +//---------------------------------------------------------------------------// +/*! + * Get pre-step point. + */ +GeantStepPointView GeantStepView::pre_step() const +{ + CELER_EXPECT(this->has_step_point(StepPoint::pre)); + return GeantStepPointView{*s_.GetPreStepPoint()}; +} + +//---------------------------------------------------------------------------// +/*! + * Get post-step point. + */ +GeantStepPointView GeantStepView::post_step() const +{ + CELER_EXPECT(this->has_step_point(StepPoint::post)); + return GeantStepPointView{*s_.GetPostStepPoint()}; +} + +//---------------------------------------------------------------------------// +/*! + * Whether the step point is valid (not null). + */ +bool GeantStepView::has_step_point(StepPoint sp) const +{ + CELER_EXPECT(sp != StepPoint::size_); + return (sp == StepPoint::pre ? s_.GetPreStepPoint() : s_.GetPostStepPoint()) + != nullptr; +} + +//---------------------------------------------------------------------------// +/*! + * Get step point by enum. + */ +GeantStepPointView GeantStepView::step_point(StepPoint sp) const +{ + CELER_EXPECT(sp != StepPoint::size_); + return sp == StepPoint::pre ? this->pre_step() : this->post_step(); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/ext/GeantTrackView.hh b/src/celeritas/ext/GeantTrackView.hh new file mode 100644 index 0000000000..18d1f37590 --- /dev/null +++ b/src/celeritas/ext/GeantTrackView.hh @@ -0,0 +1,231 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/ext/GeantTrackView.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/Types.hh" +#include "corecel/math/Quantity.hh" +#include "geocel/g4/Convert.hh" +#include "celeritas/UnitTypes.hh" + +#include "GeantParticleView.hh" +#include "GeantUnits.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +template +class GeantTrackView +{ + static_assert(W != Ownership::value, "GeantTrackView cannot own data"); +}; + +//---------------------------------------------------------------------------// +/*! + * Access track data from Geant4 with Celeritas units. + * + * This provides a uniform interface to G4Track data using Celeritas types and + * units. Geant4 data are all in double precision. + * + * The const_reference version provides read-only access, while the reference + * version adds setters. + */ +template<> +class GeantTrackView +{ + public: + //!@{ + //! \name Type aliases + using Energy = Quantity; + using real_type = double; + //!@} + + public: + // Construct from G4Track + explicit GeantTrackView(G4Track const& track) : t_(track) {} + + // Get particle definition view + inline GeantParticleView particle() const; + + // Position in native Celeritas length units + inline Real3 pos() const; + + // Momentum direction (unit vector) + inline Real3 dir() const; + + // Kinetic energy [MeV] + inline Energy energy() const; + + // Global time in native Celeritas time units + inline real_type time() const; + + //! Statistical weight + real_type weight() const { return t_.GetWeight(); } + + //! Access the G4 track directly + G4Track const& track() const { return t_; } + + //! Access the G4 track directly (const) + G4Track const& ctrack() const { return t_; } + + private: + G4Track const& t_; +}; + +//---------------------------------------------------------------------------// +/*! + * Mutable track view with setters. + */ +template<> +class GeantTrackView + : public GeantTrackView +{ + using Base = GeantTrackView; + + public: + //!@{ + //! \name Type aliases + using Energy = typename Base::Energy; + using real_type = typename Base::real_type; + //!@} + + public: + // Construct from mutable G4Track + explicit GeantTrackView(G4Track& track) : Base(track) {} + + // Bring base class accessors into scope + using Base::dir; + using Base::energy; + using Base::pos; + using Base::time; + using Base::weight; + + // Mutators + inline void pos(Real3 const& position); + inline void dir(Real3 const& direction); + inline void energy(Energy e); + inline void time(real_type t); + + //! Set statistical weight + void weight(real_type w) { this->track().SetWeight(w); } + + using Base::ctrack; + using Base::track; + //! Access mutable track reference (safe: constructed from non-const) + G4Track& track() { return const_cast(this->ctrack()); } +}; + +//---------------------------------------------------------------------------// +// TYPE ALIASES +//---------------------------------------------------------------------------// + +using GeantTrackViewConst = GeantTrackView; +using GeantTrackViewMutable = GeantTrackView; + +//---------------------------------------------------------------------------// +// DEDUCTION GUIDES +//---------------------------------------------------------------------------// + +// Deduce const_reference from const G4Track& +GeantTrackView(G4Track const&) -> GeantTrackView; + +// Deduce reference from mutable G4Track& +GeantTrackView(G4Track&) -> GeantTrackView; + +// Doxygen fails to deduce correct templated class +#if !defined(__DOXYGEN__) || __DOXYGEN__ > 0x011600 +//---------------------------------------------------------------------------// +// INLINE ACCESSOR DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Get particle definition view. + */ +GeantParticleView GeantTrackView::particle() const +{ + CELER_EXPECT(t_.GetDefinition()); + return GeantParticleView{*t_.GetDefinition()}; +} + +//---------------------------------------------------------------------------// +/*! + * Get position in native Celeritas length units. + */ +Real3 GeantTrackView::pos() const +{ + return convert_from_geant(t_.GetPosition(), clhep_length); +} + +//---------------------------------------------------------------------------// +/*! + * Get momentum direction. + */ +Real3 GeantTrackView::dir() const +{ + return convert_from_geant(t_.GetMomentumDirection(), 1); +} + +//---------------------------------------------------------------------------// +/*! + * Get kinetic energy in MeV. + */ +auto GeantTrackView::energy() const -> Energy +{ + return Energy{convert_from_geant(t_.GetKineticEnergy(), CLHEP::MeV)}; +} + +//---------------------------------------------------------------------------// +/*! + * Get global time in native Celeritas time units. + */ +real_type GeantTrackView::time() const +{ + return convert_from_geant(t_.GetGlobalTime(), clhep_time); +} + +//---------------------------------------------------------------------------// +// INLINE MUTATOR DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Set position in native Celeritas length units. + */ +void GeantTrackView::pos(Real3 const& position) +{ + this->track().SetPosition(convert_to_geant(position, clhep_length)); +} + +//---------------------------------------------------------------------------// +/*! + * Set momentum direction. + */ +void GeantTrackView::dir(Real3 const& direction) +{ + this->track().SetMomentumDirection(convert_to_geant(direction, 1)); +} + +//---------------------------------------------------------------------------// +/*! + * Set kinetic energy in MeV. + */ +void GeantTrackView::energy(Energy e) +{ + this->track().SetKineticEnergy(convert_to_geant(e.value(), CLHEP::MeV)); +} + +//---------------------------------------------------------------------------// +/*! + * Set global time in native Celeritas time units. + */ +void GeantTrackView::time(real_type t) +{ + this->track().SetGlobalTime(convert_to_geant(t, clhep_time)); +} +#endif + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/ext/detail/HitProcessor.cc b/src/celeritas/ext/detail/HitProcessor.cc index 63f75ae4bb..23c48ef7d2 100644 --- a/src/celeritas/ext/detail/HitProcessor.cc +++ b/src/celeritas/ext/detail/HitProcessor.cc @@ -25,16 +25,17 @@ #include "corecel/cont/EnumArray.hh" #include "corecel/cont/Range.hh" #include "corecel/io/Logger.hh" +#include "corecel/math/QuantityIO.hh" #include "corecel/sys/ScopedProfiling.hh" #include "corecel/sys/TraceCounter.hh" #include "geocel/GeantGeoParams.hh" -#include "geocel/g4/Convert.hh" #include "celeritas/Types.hh" #include "celeritas/user/DetectorSteps.hh" #include "celeritas/user/StepData.hh" #include "LevelTouchableUpdater.hh" -#include "../GeantUnits.hh" +#include "../GeantStepPointView.hh" +#include "../GeantStepView.hh" namespace celeritas { @@ -106,72 +107,40 @@ HitProcessor::HitProcessor(SPConstVecLV detector_volumes, CELER_LOG(debug) << "Setting up thread-local hit processor for " << detector_volumes_->size() << " sensitive detectors"; -#if G4VERSION_NUMBER >= 1103 -# define HP_CLEAR_STEP_POINT(CMD) step_->CMD(nullptr) -#else -# define HP_CLEAR_STEP_POINT(CMD) /* no "reset" before v11.0.3 */ -#endif - -#define HP_SETUP_POINT(LOWER, TITLE) \ - do \ - { \ - if (!selection.points[StepPoint::LOWER]) \ - { \ - HP_CLEAR_STEP_POINT(Reset##TITLE##StepPoint); \ - } \ - else \ - { \ - auto* sp = step_->Get##TITLE##StepPoint(); \ - sp->SetStepStatus(fUserDefinedLimit); \ - step_points_[StepPoint::LOWER] = sp; \ - } \ - } while (0) - - HP_SETUP_POINT(pre, Pre); - HP_SETUP_POINT(post, Post); -#undef HP_SETUP_POINT -#undef HP_CLEAR_STEP_POINT - - for (auto p : range(StepPoint::size_)) + GeantStepView step_view{*step_}; + for (auto sp : range(StepPoint::size_)) { - if (locate_touchable[p]) + if (!selection.points[sp]) { - // Create touchable handle for this step point - touch_handle_[p] = new G4TouchableHistory; - CELER_ASSERT(step_points_[p]); - step_points_[p]->SetTouchableHandle(touch_handle_[p]); + step_view.delete_step_point(sp); + CELER_ASSERT(!locate_touchable[sp]); } - if (locate_touchable[p] && !update_touchable_) + else { - CELER_EXPECT(selection.points[p].volume_instance_ids); - // FIXME: pass geant geo into this constructor - auto ggeo = ::celeritas::global_geant_geo().lock(); - CELER_ASSERT(ggeo); - update_touchable_ - = std::make_unique(std::move(ggeo)); + auto point_view = step_view.step_point(sp); + point_view.clear_unsupported(); + step_points_[sp] = &point_view.step_point(); + if (locate_touchable[sp]) + { + // Create touchable handle for this step point + touch_handle_[sp] = new G4TouchableHistory; + step_points_[sp]->SetTouchableHandle(touch_handle_[sp]); + if (!update_touchable_) + { + CELER_EXPECT(selection.points[sp].volume_instance_ids); + // FIXME: pass geant geo into this constructor + auto ggeo = ::celeritas::global_geant_geo().lock(); + CELER_ASSERT(ggeo); + update_touchable_ = std::make_unique( + std::move(ggeo)); + } + } } } // Set invalid values for unsupported SD attributes step_->SetNonIonizingEnergyDeposit( -std::numeric_limits::infinity()); - for (G4StepPoint* p : step_points_) - { - if (!p) - { - continue; - } - // Time since track was created - p->SetLocalTime(std::numeric_limits::infinity()); - // Time in rest frame since track was created - p->SetProperTime(std::numeric_limits::infinity()); - // Speed (TODO: use ParticleView) - p->SetVelocity(std::numeric_limits::infinity()); - // Safety distance - p->SetSafety(std::numeric_limits::infinity()); - // Polarization (default to zero) - p->SetPolarization(G4ThreeVector()); - } // Convert logical volumes (global) to sensitive detectors (thread local) detectors_.resize(detector_volumes_->size()); @@ -243,19 +212,17 @@ void HitProcessor::operator()(DetectorStepOutput const& out, size_type i) const { CELER_EXPECT(!out.detector.empty()); CELER_EXPECT(i < out.size()); -#define HP_SET(SETTER, OUT, UNITS) \ - do \ - { \ - if (!OUT.empty()) \ - { \ - SETTER(convert_to_geant(OUT[i], UNITS)); \ - } \ - } while (0) - G4LogicalVolume const* lv = this->detector_volume(out.detector[i]); - - HP_SET(step_->SetTotalEnergyDeposit, out.energy_deposition, CLHEP::MeV); - HP_SET(step_->SetStepLength, out.step_length, clhep_length); + GeantStepView step_view{*step_}; + if (!out.energy_deposition.empty()) + { + step_view.energy_deposition( + GeantStepView::Energy{out.energy_deposition[i]}); + } + if (!out.step_length.empty()) + { + step_view.step_length(out.step_length[i]); + } for (auto sp : range(StepPoint::size_)) { @@ -274,55 +241,55 @@ void HitProcessor::operator()(DetectorStepOutput const& out, size_type i) const if (CELER_UNLIKELY(!success)) { // Inconsistent touchable: skip this energy deposition - CELER_LOG_LOCAL(error) - << "Omitting energy deposition of " - << step_->GetTotalEnergyDeposit() / CLHEP::MeV << " [MeV]"; + CELER_LOG_LOCAL(error) << "Omitting energy deposition of " + << step_view.energy_deposition(); return; } } - HP_SET(g4sp->SetGlobalTime, out.points[sp].time, clhep_time); - HP_SET(g4sp->SetPosition, out.points[sp].pos, clhep_length); - HP_SET(g4sp->SetKineticEnergy, out.points[sp].energy, CLHEP::MeV); - HP_SET(g4sp->SetMomentumDirection, out.points[sp].dir, 1); + auto const& celer_sp = out.points[sp]; + GeantStepPointView sp_view{*g4sp}; + +#define HP_SET_STEP_POINT_ATTR(ATTR) \ + if (!celer_sp.ATTR.empty()) \ + { \ + sp_view.ATTR(celer_sp.ATTR[i]); \ + } + HP_SET_STEP_POINT_ATTR(time); + HP_SET_STEP_POINT_ATTR(pos); + HP_SET_STEP_POINT_ATTR(energy); + HP_SET_STEP_POINT_ATTR(dir); if (!out.weight.empty()) { - g4sp->SetWeight(out.weight[i]); + // Celeritas weight does not currently change across a step + sp_view.weight(out.weight[i]); } - G4LogicalVolume const* point_lv = [&]() -> G4LogicalVolume const* { - if (sp == StepPoint::pre) - return lv; - - // NOTE: post-step volume is only fetched if we're locating the - // touchable - if (auto* touch = g4sp->GetTouchable()) - { - // The physical volume could be null if post-step is outside - if (auto* pv = touch->GetVolume()) - { - return pv->GetLogicalVolume(); - } - } - return nullptr; - }(); +#undef HP_SET_STEP_POINT_ATTR - if (point_lv) + // Copy attributes from logical volume + if (sp == StepPoint::pre) { - // Copy attributes from logical volume - g4sp->SetMaterial(point_lv->GetMaterial()); - g4sp->SetMaterialCutsCouple(point_lv->GetMaterialCutsCouple()); - g4sp->SetSensitiveDetector(point_lv->GetSensitiveDetector()); + G4LogicalVolume const* lv = this->detector_volume(out.detector[i]); + CELER_ASSERT(lv); + // Use lv already known from the in-volume detector + sp_view.update_from_volume(*lv); + } + else + { + // Look up LV from the touchable + sp_view.update_from_volume(); } } -#undef HP_SET if (!out.particle.empty()) { G4Track& g4track = track_processor_.restore_track( out.particle[i], !out.primary_id.empty() ? out.primary_id[i] : PrimaryId{}); - this->update_track(g4track); + CELER_ASSERT(&g4track == step_->GetTrack()); + // Copy step information to the corresponding track + GeantStepView{*step_}.update_track(); } if (step_post_status_) @@ -337,51 +304,6 @@ void HitProcessor::operator()(DetectorStepOutput const& out, size_type i) const this->detector(out.detector[i])->Hit(step_); } -//---------------------------------------------------------------------------// -/*! - * Recreate the track from the particle ID and saved post-step data. - * - * This is a bit like \c G4Step::UpdateTrack . - */ -void HitProcessor::update_track(G4Track& track) const -{ - // Copy data from step to track - track.SetStepLength(step_->GetStepLength()); - - G4ParticleDefinition const& pd = *track.GetParticleDefinition(); - - for (G4StepPoint* p : step_points_) - { - if (!p) - { - continue; - } - - // Copy data from track to step points - p->SetMass(pd.GetPDGMass()); - p->SetCharge(pd.GetPDGCharge()); - } - - if (G4StepPoint* pre_step = step_points_[StepPoint::pre]) - { - // Copy data from post-step to track - track.SetTouchableHandle(pre_step->GetTouchableHandle()); - } - - if (G4StepPoint* post_step = step_points_[StepPoint::post]) - { - // Copy data from post-step to track - track.SetGlobalTime(post_step->GetGlobalTime()); - track.SetPosition(post_step->GetPosition()); - track.SetKineticEnergy(post_step->GetKineticEnergy()); - track.SetMomentumDirection(post_step->GetMomentumDirection()); - track.SetWeight(post_step->GetWeight()); - - track.SetNextTouchableHandle(post_step->GetTouchableHandle()); - track.SetVelocity(post_step->GetVelocity()); - } -} - //---------------------------------------------------------------------------// } // namespace detail } // namespace celeritas diff --git a/src/celeritas/ext/detail/HitProcessor.hh b/src/celeritas/ext/detail/HitProcessor.hh index e7959d56c7..683f7684c0 100644 --- a/src/celeritas/ext/detail/HitProcessor.hh +++ b/src/celeritas/ext/detail/HitProcessor.hh @@ -138,8 +138,6 @@ class HitProcessor //! Accumulated number of hits size_type num_hits_{0}; - - void update_track(G4Track&) const; }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/optical/SimTrackView.hh b/src/celeritas/optical/SimTrackView.hh index 77cc5451eb..2e7ac08e70 100644 --- a/src/celeritas/optical/SimTrackView.hh +++ b/src/celeritas/optical/SimTrackView.hh @@ -19,6 +19,11 @@ namespace optical //---------------------------------------------------------------------------// /*! * Simulation properties for a single track. + * + * Celeritas currently stores only the global time in the lab frame: the time + * since the creation of a track should be handled separately, and the + * "proper time" (local time in the inertial frame of the track since its + * creation) is not currently tracked. */ class SimTrackView {