diff --git a/docs/manual/confinement.md b/docs/manual/confinement.md index 69a6b8577..8398f074d 100644 --- a/docs/manual/confinement.md +++ b/docs/manual/confinement.md @@ -168,6 +168,19 @@ All sampling modes described above are available, with few notes/limitations: with a small but usually acceptable performance cost, however very large values can significantly slow down the simulation. +## Weighting modes + +For the sampling in volumes, _remage_ by default weights the probability to find +a sampled vertex in one of the volumes by the ratio of their volume to the total +sampling volume. Other weighting modes can be selected by the user with a macro +command: + +- by the mass of the volumes can with the command + . +- by the mass of a given isotope (specified by the pair of proton and neutron + number) in the materials with the command + . + ## Vertices from external files For more complicated vertex confinement _remage_ supports the possibility to @@ -202,6 +215,14 @@ capabilities described above: /RMG/Generator/Confinement/Physical/AddVolume [BCPV]\w+ ``` +This can now be used as an example to weigh the volumes not by their volume: It +could be useful to weigh the different volumes by the mass fraction of the +isotope $^{76}$Ge, i.e. for simulating double beta decays with bxdecay0: + +```geant4 +/RMG/Generator/Confinement/SampleWeightByMassIsotope 32 76 +``` + ### Active LAr region of LEGEND-200 This is a typical example of the sampling from a volume intersection. This will diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index 598dce5b5..c52c95bd9 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -567,6 +567,8 @@ 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 +* `SampleWeightByMassIsotope` – Weigh the different volumes by mass of the given isotope (specified by proton and neutron numbers) * `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 @@ -591,6 +593,28 @@ This is disabled by default * **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/SampleWeightByMassIsotope` + +Weigh the different volumes by mass of the given isotope (specified by proton and neutron numbers) + +* **Parameter** – `Z` + * **Parameter type** – `i` + * **Omittable** – `False` +* **Parameter** – `N` + * **Parameter type** – `i` + * **Omittable** – `False` +* **Allowed states** – `PreInit Idle` + ### `/RMG/Generator/Confinement/SamplingMode` Select sampling mode for volume confinement diff --git a/include/RMGPrimaryTransformer.hh b/include/RMGPrimaryTransformer.hh index 95a77a021..e4e6548a6 100644 --- a/include/RMGPrimaryTransformer.hh +++ b/include/RMGPrimaryTransformer.hh @@ -21,7 +21,7 @@ class RMGPrimaryTransformer : public G4PrimaryTransformer { public: - RMGPrimaryTransformer() : G4PrimaryTransformer() { + RMGPrimaryTransformer() { // this suppresses the ZeroPolarization warning emitted for optical photons. nWarn = 11; }; diff --git a/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index bda130ce8..105f13f9e 100644 --- a/include/RMGVertexConfinement.hh +++ b/include/RMGVertexConfinement.hh @@ -113,6 +113,12 @@ class RMGVertexConfinement : public RMGVVertexGenerator { void SetSamplingMode(SamplingMode mode) { fSamplingMode = mode; } void SetFirstSamplingVolumeType(VolumeType type) { fFirstSamplingVolumeType = type; } + void SetWeightByMass(bool mode) { fWeightByMass = mode; } + void SetWeightByMassIsotope(int z, int n) { + fWeightByMass = true; + fWeightByMassIsotopeZ = z; + fWeightByMassIsotopeN = n; + } std::vector& GetGeometricalSolidDataList() { return fGeomVolumeData; @@ -238,12 +244,15 @@ class RMGVertexConfinement : public RMGVVertexGenerator { */ void GetDirection(G4ThreeVector& dir, G4ThreeVector& pos) const; + void RecalcMass(int z, int n); + G4VPhysicalVolume* physical_volume = nullptr; G4VSolid* sampling_solid = nullptr; G4RotationMatrix rotation; G4ThreeVector translation; double volume = -1; + double mass = -1; double surface = -1; bool surface_sample = false; @@ -266,8 +275,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 @@ -279,10 +289,16 @@ class RMGVertexConfinement : public RMGVVertexGenerator { void clear() { data.clear(); } void insert(SampleableObjectCollection& other) { for (size_t i = 0; i < other.size(); ++i) this->emplace_back(other.at(i)); + this->total_volume += other.total_volume; + this->total_mass += other.total_mass; + this->total_surface += other.total_surface; } + void recalc_total(bool weigh_by_mass, int mass_isotope_z, int mass_istotope_n); + std::vector data; double total_volume = 0; + double total_mass = 0; double total_surface = 0; }; @@ -311,6 +327,9 @@ class RMGVertexConfinement : public RMGVVertexGenerator { SamplingMode fSamplingMode = SamplingMode::kUnionAll; VolumeType fFirstSamplingVolumeType = VolumeType::kUnset; + bool fWeightByMass = false; + int fWeightByMassIsotopeZ = 0; + int fWeightByMassIsotopeN = 0; bool fOnSurface = false; bool fForceContainmentCheck = false; diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index 6b0e1be4b..cbf42f35d 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -77,12 +77,36 @@ RMGVertexConfinement::SampleableObject::SampleableObject( } } this->volume = cubic_volume; + if (physvol) this->RecalcMass(0, 0); this->surface = solid->GetSurfaceArea(); } +void RMGVertexConfinement::SampleableObject::RecalcMass(int z, int n) { + if (!this->physical_volume) return; + + auto mat = this->physical_volume->GetLogicalVolume()->GetMaterial(); + double mfrac = 1.; + if (z > 0) { + mfrac = 0; + for (size_t i = 0; i < mat->GetNumberOfElements(); i++) { + auto elm = mat->GetElement(static_cast(i)); + for (size_t j = 0; j < elm->GetNumberOfIsotopes(); j++) { + auto iso = elm->GetIsotope(static_cast(j)); + if (iso->GetZ() == z && iso->GetN() == n) { + mfrac += mat->GetFractionVector()[i] * elm->GetRelativeAbundanceVector()[j]; + } + } + } + } + this->mass = this->volume * mfrac * mat->GetDensity(); +} + const RMGVertexConfinement::SampleableObject& RMGVertexConfinement::SampleableObjectCollection::SurfaceWeightedRand() const { - if (data.empty()) RMGLog::OutDev(RMGLog::fatal, "Cannot sample from an empty collection"); + if (data.empty()) [[unlikely]] + RMGLog::OutDev(RMGLog::fatal, "Cannot sample from an empty collection"); + if (this->total_surface == 0) [[unlikely]] + RMGLog::OutDev(RMGLog::fatal, "Cannot sample from a collection with no total weight"); auto choice = this->total_surface * G4UniformRand(); double w = 0; @@ -101,16 +125,25 @@ 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"); + const auto total_weight = weight_by_mass ? this->total_mass : this->total_volume; - auto choice = this->total_volume * G4UniformRand(); + if (total_weight == 0) [[unlikely]] + RMGLog::OutDev(RMGLog::fatal, "Cannot sample from a collection with no total weight"); + + 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 ", @@ -391,27 +424,50 @@ template void RMGVertexConfinement::SampleableObjectCollection::emplace_back(Args&&... args) { this->data.emplace_back(std::forward(args)...); +} - const auto& _v = this->data.back().volume; - const auto& _s = this->data.back().surface; +void RMGVertexConfinement::SampleableObjectCollection::recalc_total( + bool weigh_by_mass, + int mass_isotope_z, + int mass_istotope_n +) { - if (_v > 0) this->total_volume += _v; - else - 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." - ); + this->total_volume = 0; + this->total_mass = 0; + this->total_surface = 0; - if (_s > 0) this->total_surface += _s; - else - 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." - ); + for (auto& v : this->data) { + this->total_volume += v.volume; + if (v.volume <= 0) { + 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." + ); + } + + v.RecalcMass(mass_isotope_z, mass_istotope_n); + this->total_mass += v.mass; + if (v.mass <= 0 && weigh_by_mass) { + 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." + ); + } + + this->total_surface += v.surface; + if (v.surface <= 0) { + 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." + ); + } + } } /* ========================================================================================== */ @@ -441,7 +497,7 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { RMGLog::OutFormat( RMGLog::detail, " · '{}[{}]', volume = {}", - volume->GetName().c_str(), + volume->GetName(), volume->GetCopyNo(), std::string(G4BestUnit(fPhysicalVolumes.data.back().volume, "Volume")) ); @@ -615,12 +671,27 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { // invalidates old iterators). for (const auto& s : new_obj_from_inspection) { fPhysicalVolumes.emplace_back(s); } - RMGLog::OutFormat( - RMGLog::detail, - "Sampling from {} physical volumes, volume = {}", - fPhysicalVolumes.size(), - std::string(G4BestUnit(fPhysicalVolumes.total_volume, "Volume")) - ); + // calculate the total volume/surface/mass. + fPhysicalVolumes.recalc_total(fWeightByMass, fWeightByMassIsotopeZ, fWeightByMassIsotopeN); + + if (!fWeightByMass) { + RMGLog::OutFormat( + RMGLog::detail, + "Sampling from {} physical volumes, volume = {}", + fPhysicalVolumes.size(), + std::string(G4BestUnit(fPhysicalVolumes.total_volume, "Volume")) + ); + } else { + RMGLog::OutFormat( + RMGLog::detail, + "Sampling from {} physical volumes, mass = {}", + fPhysicalVolumes.size(), + std::string(G4BestUnit(fPhysicalVolumes.total_mass, "Mass")) + ); + } + if (fWeightByMass && fOnSurface) { + RMGLog::Out(RMGLog::fatal, "cannot sample from surface while weighting by mass"); + } } void RMGVertexConfinement::InitializeGeometricalVolumes(bool use_excluded_volumes) { @@ -689,25 +760,28 @@ void RMGVertexConfinement::InitializeGeometricalVolumes(bool use_excluded_volume volume_solids.back().native_sample = true; volume_solids.back().surface_sample = fOnSurface; - RMGLog::Out( + RMGLog::OutFormat( RMGLog::detail, - "Added geometrical solid ", + "Added geometrical solid of type '{}' with volume = {} and surface = {}", use_excluded_volumes ? "(excluded) " : " ", - "of type '", volume_solids.back().sampling_solid->GetEntityType(), - "' with volume ", - G4BestUnit(volume_solids.back().volume, "Volume"), - "and surface ", - G4BestUnit(volume_solids.back().surface, "Surface") + std::string(G4BestUnit(volume_solids.back().volume, "Volume")), + std::string(G4BestUnit(volume_solids.back().surface, "Surface")) ); } + // calculate the total volume/surface/mass. + volume_solids.recalc_total(fWeightByMass, fWeightByMassIsotopeZ, fWeightByMassIsotopeN); + RMGLog::Out( RMGLog::detail, "Will sample points in the ", fOnSurface ? "surface" : "bulk", " of geometrical volumes" ); + if (fWeightByMass && fOnSurface) { + RMGLog::Out(RMGLog::fatal, "cannot sample from surface while weighting by mass"); + } } void RMGVertexConfinement::Reset() { @@ -804,8 +878,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; @@ -893,8 +967,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; @@ -963,7 +1037,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); @@ -1094,7 +1168,6 @@ void RMGVertexConfinement::EndOfRunAction(const G4Run* run) { void RMGVertexConfinement::DefineCommands() { - fMessengers.push_back( std::make_unique( this, @@ -1118,6 +1191,24 @@ 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("SampleWeightByMassIsotope", &RMGVertexConfinement::SetWeightByMassIsotope) + .SetGuidance("Weigh the different volumes by mass of the given isotope (specified by proton and neutron numbers)") + .SetParameterName(0, "Z", false, false) + .SetParameterName(1, "N", false, false) + .SetStates(G4State_PreInit, G4State_Idle) + .SetToBeBroadcasted(true); + fMessengers.back() ->DeclareMethod("SamplingMode", &RMGVertexConfinement::SetSamplingModeString) .SetGuidance("Select sampling mode for volume confinement") diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 923345485..d7f4e1000 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -238,3 +238,7 @@ foreach(det ${_solids}) MPLCONFIGDIR=${CMAKE_SOURCE_DIR}/tests) endforeach() + +add_test(NAME confinement/weighting COMMAND ${PYTHONPATH} test_weighting.py) +set_tests_properties(confinement/weighting PROPERTIES LABELS extra ENVIRONMENT + MPLCONFIGDIR=${CMAKE_SOURCE_DIR}/tests) diff --git a/tests/confinement/gdml/geometry_weight.gdml b/tests/confinement/gdml/geometry_weight.gdml new file mode 100644 index 000000000..c4624ecd8 --- /dev/null +++ b/tests/confinement/gdml/geometry_weight.gdml @@ -0,0 +1,59 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/confinement/gdml/geometry_weight.py b/tests/confinement/gdml/geometry_weight.py new file mode 100644 index 000000000..2563b42c0 --- /dev/null +++ b/tests/confinement/gdml/geometry_weight.py @@ -0,0 +1,49 @@ +from __future__ import annotations + +from math import pi + +from pyg4ometry import gdml, geant4 + +reg = geant4.Registry() + +# dummy isotopes/elements with numbers that are easy to calculate. +iso1 = geant4.Isotope("isotope1", 32, 76, 76, reg) +iso2 = geant4.Isotope("isotope2", 32, 75, 75, reg) + +el1 = geant4.ElementIsotopeMixture("element1", "X", 2, reg) +el1.add_isotope(iso1, 0.5) +el1.add_isotope(iso2, 0.5) + +el2 = geant4.ElementIsotopeMixture("element2", "X", 2, reg) +el2.add_isotope(iso1, 0.25) +el2.add_isotope(iso2, 0.75) + +mat1 = geant4.MaterialCompound("mat1", 10, 1, reg) +mat1.add_element_massfraction(el1, 1) + +mat2 = geant4.MaterialCompound("mat2", 5, 1, reg) +mat2.add_element_massfraction(el2, 1) + +# world volume (air) +world_side = 10 # m +world_s = geant4.solid.Box("World", world_side, world_side, world_side, reg, lunit="m") +world_l = geant4.LogicalVolume(world_s, "G4_Pb", "World", reg) +reg.setWorld(world_l) + +# 2 boxes with one smaller +box_s = geant4.solid.Box("Box", 1, 1, 1, reg, lunit="m") +box1_l = geant4.LogicalVolume(box_s, mat1, "Box1", reg) +geant4.PhysicalVolume( + [pi / 4, -pi / 4, pi / 4], [2, 0, 0, "m"], box1_l, "Box1", world_l, reg +) + +box2_l = geant4.LogicalVolume(box_s, mat2, "Box2", reg) +geant4.PhysicalVolume( + [pi / 4, -pi / 4, pi / 4], [-2, 0, 0, "m"], box2_l, "Box2", world_l, reg +) + +world_l.checkOverlaps(recursive=True) + +w = gdml.Writer() +w.addDetector(reg) +w.write("geometry_weight.gdml") diff --git a/tests/confinement/test_weighting.py b/tests/confinement/test_weighting.py new file mode 100644 index 000000000..f422210b2 --- /dev/null +++ b/tests/confinement/test_weighting.py @@ -0,0 +1,75 @@ +from __future__ import annotations + +import numpy as np +from lgdo import lh5 +from remage import remage_run + + +def run_sim(run_name: str, cmd: str, gdml: str): + macro = f""" + /RMG/Geometry/RegisterDetector Scintillator Box1 1 + /RMG/Geometry/RegisterDetector Scintillator Box2 2 + + /RMG/Output/NtupleUseVolumeName true + + /run/initialize + + /RMG/Generator/Confine Volume + /RMG/Generator/Confinement/ForceContainmentCheck true + /RMG/Generator/Confinement/Physical/AddVolume Box1 + /RMG/Generator/Confinement/Physical/AddVolume Box2 + + {cmd} + + /RMG/Generator/Select GPS + /gps/particle e- + /gps/energy 1 eV + + /run/beamOn 20000""" + + remage_run( + macro.split("\n"), + gdml_files=gdml, + output=f"{run_name}.lh5", + overwrite_output=True, + ) + + +# Define the different detector registration commands to test +gdml_default = "gdml/geometry_weight.gdml" +runs = [ + ( + "test-weight-volume", + "/RMG/Generator/Confinement/SampleWeightByMass false", + ), + ( + "test-weight-mass", + "/RMG/Generator/Confinement/SampleWeightByMass true", + ), + ( + "test-weight-isotope", + "/RMG/Generator/Confinement/SampleWeightByMassIsotope 32 76", + ), +] +file_vols = "test-weight-volume.lh5" +file_mass = "test-weight-mass.lh5" +file_isot = "test-weight-isotope.lh5" + +# Run each simulation +for run_name, register_command in runs: + run_sim(run_name, register_command, gdml_default) + + +def get_ratio(file): + box1 = lh5.read_as("/stp/Box1", file, "ak") + box2 = lh5.read_as("/stp/Box2", file, "ak") + print(file, len(box1) / len(box2)) + return len(box1) / len(box2) + + +exp_vol = 1 +assert np.abs(get_ratio(file_vols) - exp_vol) / exp_vol < 0.1 +exp_mass = 10 / 5 +assert np.abs(get_ratio(file_mass) - exp_mass) / exp_mass < 0.1 +exp_isot = (10 * 0.5) / (5 * 0.25) +assert np.abs(get_ratio(file_isot) - exp_isot) / exp_isot < 0.1