diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index 1dd42c1f..9efac083 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -496,6 +496,7 @@ Commands for controlling primary confinement * `Reset` – Reset all parameters of vertex confinement, so that it can be reconfigured. * `SampleOnSurface` – If true (or omitted argument), sample on the surface of solids +* `SampleWeightByMass` – If true (or omitted argument), weigh the different volumes by mass and not by volume * `SamplingMode` – Select sampling mode for volume confinement * `FirstSamplingVolume` – Select the type of volume which will be sampled first for intersections * `MaxSamplingTrials` – Set maximum number of attempts for sampling primary positions in a volume @@ -518,6 +519,16 @@ If true (or omitted argument), sample on the surface of solids * **Default value** – `true` * **Allowed states** – `PreInit Idle` +### `/RMG/Generator/Confinement/SampleWeightByMass` + +If true (or omitted argument), weigh the different volumes by mass and not by volume + +* **Parameter** – `boolean` + * **Parameter type** – `b` + * **Omittable** – `True` + * **Default value** – `true` +* **Allowed states** – `PreInit Idle` + ### `/RMG/Generator/Confinement/SamplingMode` Select sampling mode for volume confinement diff --git a/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index bda130ce..083ec958 100644 --- a/include/RMGVertexConfinement.hh +++ b/include/RMGVertexConfinement.hh @@ -113,6 +113,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator { void SetSamplingMode(SamplingMode mode) { fSamplingMode = mode; } void SetFirstSamplingVolumeType(VolumeType type) { fFirstSamplingVolumeType = type; } + void SetWeightByMass(bool mode) { fWeightByMass = mode; } std::vector& GetGeometricalSolidDataList() { return fGeomVolumeData; @@ -244,6 +245,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator { G4ThreeVector translation; double volume = -1; + double mass = -1; double surface = -1; bool surface_sample = false; @@ -266,8 +268,9 @@ class RMGVertexConfinement : public RMGVVertexGenerator { /** @brief Select a @c SampleableObject from the collection, weighted by volume. * @returns a reference to the chosen @c SampleableObject . + * @param weight_by_mass A flag of whether the volume weighting should be done by mass and not by volume. */ - [[nodiscard]] const SampleableObject& VolumeWeightedRand() const; + [[nodiscard]] const SampleableObject& VolumeWeightedRand(bool weight_by_mass) const; [[nodiscard]] bool IsInside(const G4ThreeVector& vertex) const; // emulate @c std::vector @@ -283,7 +286,11 @@ class RMGVertexConfinement : public RMGVVertexGenerator { std::vector data; double total_volume = 0; + bool total_volume_all = true; + double total_mass = 0; + bool total_mass_all = true; double total_surface = 0; + bool total_surface_all = true; }; private: @@ -311,6 +318,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator { SamplingMode fSamplingMode = SamplingMode::kUnionAll; VolumeType fFirstSamplingVolumeType = VolumeType::kUnset; + bool fWeightByMass = false; bool fOnSurface = false; bool fForceContainmentCheck = false; diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index e69623a8..6f84f7cc 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -77,6 +77,9 @@ RMGVertexConfinement::SampleableObject::SampleableObject( } } this->volume = cubic_volume; + if (physvol) { + this->mass = cubic_volume * physvol->GetLogicalVolume()->GetMaterial()->GetDensity(); + } this->surface = solid->GetSurfaceArea(); } @@ -101,16 +104,31 @@ const RMGVertexConfinement::SampleableObject& RMGVertexConfinement::SampleableOb return this->data.back(); } -const RMGVertexConfinement::SampleableObject& RMGVertexConfinement::SampleableObjectCollection::VolumeWeightedRand() const { +const RMGVertexConfinement::SampleableObject& RMGVertexConfinement::SampleableObjectCollection::VolumeWeightedRand( + bool weight_by_mass +) const { + + if (data.empty()) [[unlikely]] + RMGLog::OutDev(RMGLog::fatal, "Cannot sample from an empty collection"); - if (data.empty()) RMGLog::OutDev(RMGLog::fatal, "Cannot sample from an empty collection"); + if (weight_by_mass && !this->total_mass_all) [[unlikely]] { + RMGLog::Out( + RMGLog::fatal, + "One of the sampled solids has no mass attribute, ", + "will not add it to the total mass of the collection. ", + "this will affect sampling from multiple solids." + ); + } + + const auto total_weight = weight_by_mass ? this->total_mass : this->total_volume; - auto choice = this->total_volume * G4UniformRand(); + auto choice = total_weight * G4UniformRand(); double w = 0; for (const auto& o : this->data) { - if (choice > w and choice <= w + o.volume) return o; - w += o.volume; - if (w >= this->total_volume) { + const auto weight = weight_by_mass ? o.mass : o.volume; + if (choice > w and choice <= w + weight) return o; + w += weight; + if (w >= total_weight) [[unlikely]] { RMGLog::Out( RMGLog::error, "Sampling from collection of sampleables failed unexpectedly ", @@ -389,25 +407,32 @@ void RMGVertexConfinement::SampleableObjectCollection::emplace_back(Args&&... ar this->data.emplace_back(std::forward(args)...); const auto& _v = this->data.back().volume; + const auto& _m = this->data.back().mass; const auto& _s = this->data.back().surface; if (_v > 0) this->total_volume += _v; - else + else { + this->total_volume_all = false; RMGLog::Out( RMGLog::error, "One of the sampled solids has no volume attribute, ", "will not add it to the total volume of the collection. ", "this will affect sampling from multiple solids." ); + } + if (_m > 0) this->total_mass += _m; + else this->total_mass_all = false; if (_s > 0) this->total_surface += _s; - else + else { + this->total_surface_all = false; RMGLog::Out( RMGLog::error, "One of the sampled solids has no surface attribute, ", "will not add it to the total surface of the collection. ", "this will affect sampling from multiple solids." ); + } } /* ========================================================================================== */ @@ -796,8 +821,8 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { ? fGeomVolumeSolids.total_volume > fPhysicalVolumes.total_volume : fFirstSamplingVolumeType == VolumeType::kPhysical; - choice_nonconst = physical_first ? fPhysicalVolumes.VolumeWeightedRand() - : fGeomVolumeSolids.VolumeWeightedRand(); + choice_nonconst = physical_first ? fPhysicalVolumes.VolumeWeightedRand(fWeightByMass) + : fGeomVolumeSolids.VolumeWeightedRand(fWeightByMass); } const auto& choice = choice_nonconst; @@ -885,8 +910,8 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { else if (has_geometrical && not has_physical) physical_first = false; else physical_first = (fFirstSamplingVolumeType == VolumeType::kPhysical); - choice_nonconst = physical_first ? fPhysicalVolumes.VolumeWeightedRand() - : fGeomVolumeSolids.VolumeWeightedRand(); + choice_nonconst = physical_first ? fPhysicalVolumes.VolumeWeightedRand(fWeightByMass) + : fGeomVolumeSolids.VolumeWeightedRand(fWeightByMass); } const auto& choice = choice_nonconst; @@ -955,7 +980,7 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { // chose a volume to sample from const auto choice = fOnSurface ? all_volumes.SurfaceWeightedRand() - : all_volumes.VolumeWeightedRand(); + : all_volumes.VolumeWeightedRand(fWeightByMass); // do the sampling bool success = choice.Sample(vertex, fMaxAttempts, fForceContainmentCheck, fTrials); @@ -1109,6 +1134,16 @@ void RMGVertexConfinement::DefineCommands() { .SetStates(G4State_PreInit, G4State_Idle) .SetToBeBroadcasted(true); + fMessengers.back() + ->DeclareMethod("SampleWeightByMass", &RMGVertexConfinement::SetWeightByMass) + .SetGuidance( + "If true (or omitted argument), weigh the different volumes by mass and not by volume" + ) + .SetParameterName("boolean", true) + .SetDefaultValue("true") + .SetStates(G4State_PreInit, G4State_Idle) + .SetToBeBroadcasted(true); + fMessengers.back() ->DeclareMethod("SamplingMode", &RMGVertexConfinement::SetSamplingModeString) .SetGuidance("Select sampling mode for volume confinement")