diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index b7bcc43c..abfc5041 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -37,7 +37,7 @@ jobs: run: | mkdir build cd build - cmake -DRMG_BUILD_EXAMPLES=1 .. + cmake -DRMG_BUILD_EXAMPLES=1 -DBUILD_TESTING=ON -DRMG_BUILD_DOCS=OFF .. make make install diff --git a/CMakeLists.txt b/CMakeLists.txt index c58a5f74..99a23183 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -199,6 +199,7 @@ if(RMG_BUILD_EXAMPLES) add_subdirectory(examples) endif() +message(STATUS "Test/validation suite support: BUILD_TESTING=${BUILD_TESTING}") include(CTest) add_subdirectory(tests) diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index ec11670e..79bb4e2d 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -956,6 +956,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 * `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 * `ForceContainmentCheck` – If true (or omitted argument), perform a containment check even after sampling from a natively sampleable object. This is only an extra sanity check that does not alter the behaviour. @@ -978,7 +979,16 @@ Select sampling mode for volume confinement * **Parameter** – `mode` * **Parameter type** – `s` * **Omittable** – `False` - * **Candidates** – `IntersectPhysicalWithGeometrical UnionAll` + * **Candidates** – `IntersectPhysicalWithGeometrical UnionAll SubtractGeometrical` + +### `/RMG/Generator/Confinement/FirstSamplingVolume` + +Select the type of volume which will be sampled first for intersections + +* **Parameter** – `type` + * **Parameter type** – `s` + * **Omittable** – `False` + * **Candidates** – `Physical Geometrical Unset` ### `/RMG/Generator/Confinement/MaxSamplingTrials` @@ -1031,6 +1041,7 @@ Commands for setting geometrical volumes up for primary confinement **Commands:** * `AddSolid` – Add geometrical solid to sample primaries from +* `AddExcludedSolid` – Add geometrical solid to exclude samples from * `CenterPositionX` – Set center position (X coordinate) * `CenterPositionY` – Set center position (Y coordinate) * `CenterPositionZ` – Set center position (Z coordinate) @@ -1039,6 +1050,15 @@ Commands for setting geometrical volumes up for primary confinement Add geometrical solid to sample primaries from +* **Parameter** – `solid` + * **Parameter type** – `s` + * **Omittable** – `False` + * **Candidates** – `Sphere Cylinder Box` + +### `/RMG/Generator/Confinement/Geometrical/AddExcludedSolid` + +Add geometrical solid to exclude samples from + * **Parameter** – `solid` * **Parameter type** – `s` * **Omittable** – `False` diff --git a/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index 4074a004..48d56227 100644 --- a/include/RMGVertexConfinement.hh +++ b/include/RMGVertexConfinement.hh @@ -60,9 +60,15 @@ class RMGVertexConfinement : public RMGVVertexGenerator { enum class SamplingMode { kIntersectPhysicalWithGeometrical, - kUnionAll + kUnionAll, + kSubtractGeometrical, }; + enum class VolumeType { + kPhysical, + kGeometrical, + kUnset, + }; RMGVertexConfinement(); void BeginOfRunAction(const G4Run* run) override; @@ -79,6 +85,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator { void Reset(); inline void SetSamplingMode(SamplingMode mode) { fSamplingMode = mode; } + inline void SetFirstSamplingVolumeType(VolumeType type) { fFirstSamplingVolumeType = type; } inline std::vector& GetGeometricalSolidDataList() { return fGeomVolumeData; @@ -93,6 +100,8 @@ class RMGVertexConfinement : public RMGVVertexGenerator { // NOTE: G4 volume/solid pointers should be fully owned by G4, avoid trying to delete them ~SampleableObject() = default; [[nodiscard]] bool IsInside(const G4ThreeVector& vertex) const; + [[nodiscard]] bool Sample(G4ThreeVector& vertex, int max_attempts, bool sample_on_surface, + bool force_containment_check, long int& n_trials) const; G4VPhysicalVolume* physical_volume = nullptr; G4VSolid* sampling_solid = nullptr; @@ -123,7 +132,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator { for (size_t i = 0; i < other.size(); ++i) this->emplace_back(other.at(i)); } - std::vector data; + std::vector data = {}; double total_volume = 0; double total_surface = 0; }; @@ -144,13 +153,14 @@ class RMGVertexConfinement : public RMGVVertexGenerator { }; void InitializePhysicalVolumes(); - void InitializeGeometricalVolumes(); + void InitializeGeometricalVolumes(bool useExcludedVolumes); bool ActualGenerateVertex(G4ThreeVector& v); std::vector fPhysicalVolumeNameRegexes; std::vector fPhysicalVolumeCopyNrRegexes; std::vector fGeomVolumeData; + std::vector fExcludedGeomVolumeData; static G4Mutex fGeometryMutex; // the final geometry data is shared between all threads and protected by fGeometryMutex. @@ -159,19 +169,27 @@ class RMGVertexConfinement : public RMGVVertexGenerator { // mutate the G4SolidStore temporarily, e.g. G4SubstractionSolid::GetCubicVolume() static SampleableObjectCollection fPhysicalVolumes; static SampleableObjectCollection fGeomVolumeSolids; + static SampleableObjectCollection fExcludedGeomVolumeSolids; + static bool fVolumesInitialized; SamplingMode fSamplingMode = SamplingMode::kUnionAll; + VolumeType fFirstSamplingVolumeType = VolumeType::kUnset; + bool fOnSurface = false; bool fForceContainmentCheck = false; - + bool fLastSolidExcluded = false; // counters used for the current run. long fTrials = 0; std::chrono::nanoseconds fVertexGenerationTime; std::vector> fMessengers; void SetSamplingModeString(std::string mode); + void SetFirstSamplingVolumeTypeString(std::string type); + void AddGeometricalVolumeString(std::string solid); + void AddExcludedGeometricalVolumeString(std::string solid); + GenericGeometricalSolidData& SafeBack( std::optional solid_type = std::nullopt); diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 36248d83..6a7ca91d 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -62,7 +62,7 @@ add_custom_command( cp ${CMAKE_CURRENT_BINARY_DIR}/cpp_config.build.py # now we want to use the cpp_config for the build area ${CMAKE_CURRENT_SOURCE_DIR}/remage/cpp_config.py - COMMAND ${VENV_DIR}/bin/python -m uv -q pip install --reinstall ${_pkg_install} + COMMAND ${VENV_DIR}/bin/python -m uv pip -q install --reinstall ${_pkg_install} DEPENDS python-virtualenv ${PYTHON_SOURCES} COMMENT "Installing remage Python wrapper into the virtual environment") diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index fc860329..91adc083 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -42,7 +42,9 @@ G4Mutex RMGVertexConfinement::fGeometryMutex = G4MUTEX_INITIALIZER; // state shared between threads and protected by fGeometryMutex. See .hh file for details! RMGVertexConfinement::SampleableObjectCollection RMGVertexConfinement::fGeomVolumeSolids = {}; +RMGVertexConfinement::SampleableObjectCollection RMGVertexConfinement::fExcludedGeomVolumeSolids = {}; RMGVertexConfinement::SampleableObjectCollection RMGVertexConfinement::fPhysicalVolumes = {}; + bool RMGVertexConfinement::fVolumesInitialized = false; // This structure must contain at least a non-null pointer, between the first @@ -119,6 +121,54 @@ bool RMGVertexConfinement::SampleableObject::IsInside(const G4ThreeVector& verte return false; } +bool RMGVertexConfinement::SampleableObject::Sample(G4ThreeVector& vertex, int max_attempts, + bool sample_on_surface, bool force_containment_check, long int& n_trials) const { + + if (this->physical_volume) { + RMGLog::OutFormatDev(RMGLog::debug, "Chosen random volume: '{}[{}]'", + this->physical_volume->GetName(), this->physical_volume->GetCopyNo()); + } else { + RMGLog::OutFormatDev(RMGLog::debug, "Chosen random volume: '{}'", + this->sampling_solid->GetName()); + } + RMGLog::OutDev(RMGLog::debug, "Maximum attempts to find a good vertex: ", max_attempts); + + int calls = 0; + + + if (this->containment_check) { + vertex = this->translation + + this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, sample_on_surface); + + while (!this->IsInside(vertex) and calls++ < max_attempts) { + n_trials++; + vertex = this->translation + + this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, sample_on_surface); + RMGLog::OutDev(RMGLog::debug, "Vertex was not inside, new vertex: ", vertex / CLHEP::cm, " cm"); + } + if (calls >= max_attempts) { + RMGLog::Out(RMGLog::error, "Exceeded maximum number of allowed iterations (", max_attempts, + "), check that your volumes are efficiently ", + "sampleable and try, eventually, to increase the threshold through the dedicated ", + "macro command. Returning dummy vertex"); + return false; + } + } else { + vertex = this->translation + + this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, sample_on_surface); + RMGLog::OutDev(RMGLog::debug, "Generated vertex: ", vertex / CLHEP::cm, " cm"); + if (force_containment_check && !this->IsInside(vertex)) { + RMGLog::OutDev(RMGLog::error, + "Generated vertex not inside sampling volumes (forced containment check): ", + vertex / CLHEP::cm, " cm"); + } + } + + RMGLog::OutDev(RMGLog::debug, "Found good vertex ", vertex / CLHEP::cm, " cm", " after ", calls, + " iterations, returning"); + return true; +} + bool RMGVertexConfinement::SampleableObjectCollection::IsInside(const G4ThreeVector& vertex) const { auto navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); @@ -128,6 +178,7 @@ bool RMGVertexConfinement::SampleableObjectCollection::IsInside(const G4ThreeVec return false; } + template void RMGVertexConfinement::SampleableObjectCollection::emplace_back(Args&&... args) { @@ -153,6 +204,7 @@ void RMGVertexConfinement::SampleableObjectCollection::emplace_back(Args&&... ar RMGVertexConfinement::RMGVertexConfinement() : RMGVVertexGenerator("VolumeConfinement") { + RMGLog::OutFormat(RMGLog::detail, "Create RMGVertexConfinment object"); this->DefineCommands(); } @@ -361,42 +413,46 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { fPhysicalVolumes.size(), std::string(G4BestUnit(fPhysicalVolumes.total_volume, "Volume"))); } -void RMGVertexConfinement::InitializeGeometricalVolumes() { +void RMGVertexConfinement::InitializeGeometricalVolumes(bool use_excluded_volumes) { + + // Select the appropriate containers based on the option - // if collections are not empty, assume initialization to be already done and skip - if (!fGeomVolumeSolids.empty() or fGeomVolumeData.empty()) return; + auto& volume_solids = use_excluded_volumes ? fExcludedGeomVolumeSolids : fGeomVolumeSolids; + const auto& volume_data = use_excluded_volumes ? fExcludedGeomVolumeData : fGeomVolumeData; - // no physical volume is specified nor at initialization or later - for (const auto& d : fGeomVolumeData) { + // If collections are not empty, assume initialization to be already done and skip + if (!volume_solids.empty() || volume_data.empty()) return; + + // Initialize volumes based on the data + for (const auto& d : volume_data) { if (d.solid_type == GeometricalSolidType::kSphere) { - fGeomVolumeSolids.emplace_back(nullptr, G4RotationMatrix(), d.volume_center, - new G4Sphere("RMGVertexConfinement::fGeomSamplingShape::Sphere/" + - std::to_string(fGeomVolumeSolids.size() + 1), + volume_solids.emplace_back(nullptr, G4RotationMatrix(), d.volume_center, + new G4Sphere("RMGVertexConfinement::SamplingShape::Sphere/" + + std::to_string(volume_solids.size() + 1), d.sphere_inner_radius, d.sphere_outer_radius, 0, CLHEP::twopi, 0, CLHEP::pi)); } else if (d.solid_type == GeometricalSolidType::kCylinder) { - fGeomVolumeSolids.emplace_back(nullptr, G4RotationMatrix(), d.volume_center, - new G4Tubs("RMGVertexConfinement::fGeomSamplingShape::Cylinder/" + - std::to_string(fGeomVolumeSolids.size() + 1), + volume_solids.emplace_back(nullptr, G4RotationMatrix(), d.volume_center, + new G4Tubs("RMGVertexConfinement::SamplingShape::Cylinder/" + + std::to_string(volume_solids.size() + 1), d.cylinder_inner_radius, d.cylinder_outer_radius, 0.5 * d.cylinder_height, d.cylinder_starting_angle, d.cylinder_spanning_angle)); } else if (d.solid_type == GeometricalSolidType::kBox) { - fGeomVolumeSolids.emplace_back(nullptr, G4RotationMatrix(), d.volume_center, - new G4Box("RMGVertexConfinement::fGeomSamplingShape::Box/" + - std::to_string(fGeomVolumeSolids.size() + 1), + volume_solids.emplace_back(nullptr, G4RotationMatrix(), d.volume_center, + new G4Box("RMGVertexConfinement::SamplingShape::Box/" + + std::to_string(volume_solids.size() + 1), 0.5 * d.box_x_length, 0.5 * d.box_y_length, 0.5 * d.box_z_length)); - } - // else if (...) - else { + } else { RMGLog::OutFormat(RMGLog::error, "Geometrical solid '{}' not known! (Implement me?)", d.solid_type); } - fGeomVolumeSolids.back().containment_check = false; + volume_solids.back().containment_check = false; - RMGLog::Out(RMGLog::detail, "Added geometrical solid of type '", - fGeomVolumeSolids.back().sampling_solid->GetEntityType(), "' with volume ", - G4BestUnit(fGeomVolumeSolids.back().volume, "Volume"), "and surface ", - G4BestUnit(fGeomVolumeSolids.back().surface, "Surface")); + RMGLog::Out(RMGLog::detail, "Added geometrical solid ", + 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")); } RMGLog::Out(RMGLog::detail, "Will sample points in the ", fOnSurface ? "surface" : "bulk", @@ -405,15 +461,22 @@ void RMGVertexConfinement::InitializeGeometricalVolumes() { void RMGVertexConfinement::Reset() { // take lock, just not to race with the fVolumesInitialized flag elsewhere. + + G4AutoLock lock(&fGeometryMutex); fVolumesInitialized = true; fPhysicalVolumes.clear(); fGeomVolumeSolids.clear(); + fExcludedGeomVolumeSolids.clear(); fPhysicalVolumeNameRegexes.clear(); fPhysicalVolumeCopyNrRegexes.clear(); fGeomVolumeData.clear(); + fExcludedGeomVolumeData.clear(); + fSamplingMode = SamplingMode::kUnionAll; + fFirstSamplingVolumeType = VolumeType::kUnset; + fOnSurface = false; fForceContainmentCheck = false; } @@ -437,7 +500,9 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { G4AutoLock lock(&fGeometryMutex); if (!fVolumesInitialized) { this->InitializePhysicalVolumes(); - this->InitializeGeometricalVolumes(); + this->InitializeGeometricalVolumes(true); + this->InitializeGeometricalVolumes(false); + fVolumesInitialized = true; } } @@ -458,49 +523,47 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { "either no physical or no geometrical volumes have been added"); } - // choose a volume component randomly - SampleableObject choice_nonconst; - bool physical_first; - if (fOnSurface) { - physical_first = fGeomVolumeSolids.total_surface > fPhysicalVolumes.total_surface; - choice_nonconst = physical_first ? fPhysicalVolumes.SurfaceWeightedRand() - : fGeomVolumeSolids.SurfaceWeightedRand(); - } else { - physical_first = fGeomVolumeSolids.total_volume > fPhysicalVolumes.total_volume; - choice_nonconst = physical_first ? fPhysicalVolumes.VolumeWeightedRand() - : fGeomVolumeSolids.VolumeWeightedRand(); - } - const auto& choice = choice_nonconst; - - // shoot in the first region - while (calls++ < RMGVVertexGenerator::fMaxAttempts) { - fTrials++; - - if (choice.containment_check) { // this can effectively happen only with physical volumes, at the moment - while (!fPhysicalVolumes.IsInside(vertex) and calls++ < RMGVVertexGenerator::fMaxAttempts) { - fTrials++; - vertex = choice.translation + - choice.rotation * RMGGeneratorUtil::rand(choice.sampling_solid, fOnSurface); - } - if (calls >= RMGVVertexGenerator::fMaxAttempts) { - RMGLog::Out(RMGLog::error, "Exceeded maximum number of allowed iterations (", - RMGVVertexGenerator::fMaxAttempts, "), check that your volumes are efficiently ", - "sampleable and try, eventually, to increase the threshold through the dedicated ", - "macro command. Returning dummy vertex"); - vertex = RMGVVertexGenerator::kDummyPrimaryPosition; - return false; - } + int calls = 0; + + while (calls <= RMGVVertexGenerator::fMaxAttempts) { + calls++; + + // chose a volume + SampleableObject choice_nonconst; + bool physical_first; + + // for surface events the user has to chose which volume type to sample first + if (fOnSurface) { + + if (fFirstSamplingVolumeType == VolumeType::kUnset) + RMGLog::Out(RMGLog::fatal, "For surface sampling vertex confinment mode it is required ", + "to set the type of volume to sample first (geometrical or physical)."); + + + physical_first = fFirstSamplingVolumeType == VolumeType::kPhysical; + choice_nonconst = physical_first ? fPhysicalVolumes.SurfaceWeightedRand() + : fGeomVolumeSolids.SurfaceWeightedRand(); + } else { - vertex = choice.translation + - choice.rotation * RMGGeneratorUtil::rand(choice.sampling_solid, fOnSurface); - if (fForceContainmentCheck) { - auto is_inside = physical_first ? fPhysicalVolumes.IsInside(vertex) - : fGeomVolumeSolids.IsInside(vertex); - if (!is_inside) - RMGLog::OutDev(RMGLog::error, - "Generated vertex not inside sampling volumes (forced containment check): ", - vertex / CLHEP::cm, " cm"); - } + // for volume sampling the user can specify the volume to sample first else the set with smaller total volume is used + physical_first = fFirstSamplingVolumeType == VolumeType::kUnset + ? fGeomVolumeSolids.total_volume > fPhysicalVolumes.total_volume + : fFirstSamplingVolumeType == VolumeType::kPhysical; + + choice_nonconst = physical_first ? fPhysicalVolumes.VolumeWeightedRand() + : fGeomVolumeSolids.VolumeWeightedRand(); + } + + const auto& choice = choice_nonconst; + + // generate a candidate vertex + bool success = + choice.Sample(vertex, fMaxAttempts, fOnSurface, fForceContainmentCheck, fTrials); + + if (!success) { + RMGLog::Out(RMGLog::error, "Sampling unsuccessful return dummy vertex"); + vertex = RMGVVertexGenerator::kDummyPrimaryPosition; + return false; } // is it also in the other volume class (geometrical/physical)? @@ -523,6 +586,104 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { return false; } + case SamplingMode::kSubtractGeometrical: { + // strategy: sample a point in the geometrical user volume or the + // physical user volume depending on the smallest surface/volume. Then + // check if it is inside the other volume + + if (fGeomVolumeSolids.empty() and fPhysicalVolumes.empty()) { + RMGLog::Out(RMGLog::fatal, "'SubtractGeometrical' mode is set but ", + " no physical and no geometrical volumes have been added"); + } + + bool has_physical = not fPhysicalVolumes.empty(); + bool has_geometrical = not fGeomVolumeSolids.empty(); + + if (fFirstSamplingVolumeType == VolumeType::kPhysical and not has_physical) + RMGLog::Out(RMGLog::fatal, "Cannot sample first from physical volumes when none are set."); + + if (fFirstSamplingVolumeType == VolumeType::kGeometrical and not has_geometrical) + RMGLog::Out(RMGLog::fatal, + "Cannot sample first from geometrical volumes when none are set."); + + int calls = 0; + + while (calls <= RMGVVertexGenerator::fMaxAttempts) { + calls++; + + // chose a volume + SampleableObject choice_nonconst; + bool physical_first; + + if (fOnSurface) { + + // check the user has set which volume to samplg first from + if (fFirstSamplingVolumeType == VolumeType::kUnset) + RMGLog::Out(RMGLog::fatal, "For surface sampling vertex confinment mode it is required ", + "to set the type of volume to sample first (geometrical or physical)."); + + physical_first = fFirstSamplingVolumeType == VolumeType::kPhysical; + + + choice_nonconst = physical_first ? fPhysicalVolumes.SurfaceWeightedRand() + : fGeomVolumeSolids.SurfaceWeightedRand(); + } else { + + + // if both physical and geometrical volumes are present and order is not set chose based on the total volume + if (fFirstSamplingVolumeType == VolumeType::kUnset && has_geometrical && has_physical) + physical_first = fGeomVolumeSolids.total_volume > fPhysicalVolumes.total_volume; + else if (has_physical && not has_geometrical) physical_first = true; + else if (has_geometrical && not has_physical) physical_first = false; + else fFirstSamplingVolumeType == VolumeType::kPhysical; + + choice_nonconst = physical_first ? fPhysicalVolumes.VolumeWeightedRand() + : fGeomVolumeSolids.VolumeWeightedRand(); + } + + const auto& choice = choice_nonconst; + + // generate a candidate vertex + bool success = + choice.Sample(vertex, fMaxAttempts, fOnSurface, fForceContainmentCheck, fTrials); + + if (!success) { + RMGLog::Out(RMGLog::error, "Sampling unsuccessful return dummy vertex"); + vertex = RMGVVertexGenerator::kDummyPrimaryPosition; + return false; + } + + // check for intersections + bool accept = false; + + if (physical_first && has_geometrical) { + if (fGeomVolumeSolids.IsInside(vertex)) accept = true; + } else if (not physical_first && has_physical) { + if (fPhysicalVolumes.IsInside(vertex)) accept = true; + } else { + accept = true; + } + RMGLog::Out(RMGLog::debug, accept ? "Chosen vertex passes intersection criteria " + : "Chosen vertex fails intersection criteria. "); + + // now check for subtractions + + if (accept && !fExcludedGeomVolumeSolids.IsInside(vertex)) return true; + + RMGLog::Out(RMGLog::debug, "Chosen vertex fails intersection criteria."); + } + + if (calls >= RMGVVertexGenerator::fMaxAttempts) { + RMGLog::Out(RMGLog::error, "Exceeded maximum number of allowed iterations (", + RMGVVertexGenerator::fMaxAttempts, "), check that your volumes are efficiently ", + "sampleable and try, eventually, to increase the threshold through the dedicated ", + "macro command. Returning dummy vertex"); + } + + // everything has failed so return the dummy vertex + vertex = RMGVVertexGenerator::kDummyPrimaryPosition; + return false; + } case SamplingMode::kUnionAll: { // strategy: just sample uniformly in/on all geometrical and physical volumes @@ -535,68 +696,19 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { "no physical or no geometrical volumes have been added"); } + // chose a volume to sample from const auto choice = fOnSurface ? all_volumes.SurfaceWeightedRand() : all_volumes.VolumeWeightedRand(); - if (choice.physical_volume) { - RMGLog::OutFormatDev(RMGLog::debug, "Chosen random volume: '{}[{}]'", - choice.physical_volume->GetName(), choice.physical_volume->GetCopyNo()); - } else { - RMGLog::OutFormatDev(RMGLog::debug, "Chosen random volume: '{}'", - choice.sampling_solid->GetName()); - } - RMGLog::OutDev(RMGLog::debug, - "Maximum attempts to find a good vertex: ", RMGVVertexGenerator::fMaxAttempts); - - while (calls++ < RMGVVertexGenerator::fMaxAttempts) { - fTrials++; - - - if (choice.containment_check) { - vertex = choice.translation + - choice.rotation * RMGGeneratorUtil::rand(choice.sampling_solid, fOnSurface); - while (!choice.IsInside(vertex) and calls++ < RMGVVertexGenerator::fMaxAttempts) { - fTrials++; - vertex = choice.translation + - choice.rotation * RMGGeneratorUtil::rand(choice.sampling_solid, fOnSurface); - RMGLog::OutDev(RMGLog::debug, "Vertex was not inside, new vertex: ", vertex / CLHEP::cm, - " cm"); - } - if (calls >= RMGVVertexGenerator::fMaxAttempts) { - RMGLog::Out(RMGLog::error, "Exceeded maximum number of allowed iterations (", - RMGVVertexGenerator::fMaxAttempts, "), check that your volumes are efficiently ", - "sampleable and try, eventually, to increase the threshold through the dedicated ", - "macro command. Returning dummy vertex"); - vertex = RMGVVertexGenerator::kDummyPrimaryPosition; - return false; - } - } else { - vertex = choice.translation + - choice.rotation * RMGGeneratorUtil::rand(choice.sampling_solid, fOnSurface); - RMGLog::OutDev(RMGLog::debug, "Generated vertex: ", vertex / CLHEP::cm, " cm"); - if (fForceContainmentCheck && !choice.IsInside(vertex)) { - RMGLog::OutDev(RMGLog::error, - "Generated vertex not inside sampling volumes (forced containment check): ", - vertex / CLHEP::cm, " cm"); - } - } - - RMGLog::OutDev(RMGLog::debug, "Found good vertex ", vertex / CLHEP::cm, " cm", " after ", - calls, " iterations, returning"); - return true; - } + // do the sampling + bool success = choice.Sample(vertex, fMaxAttempts, fOnSurface, fForceContainmentCheck, fTrials); - if (calls >= RMGVVertexGenerator::fMaxAttempts) { - RMGLog::Out(RMGLog::error, "Exceeded maximum number of allowed iterations (", - RMGVVertexGenerator::fMaxAttempts, - "), check that your volumes are efficiently sampleable and ", - "try, eventually, to increase the threshold through the dedicated macro command. " - "Returning dummy vertex"); + if (!success) { + RMGLog::Out(RMGLog::error, "Sampling unsuccessful return dummy vertex"); + vertex = RMGVVertexGenerator::kDummyPrimaryPosition; + return false; } - - // everything has failed so return the dummy vertex - vertex = RMGVVertexGenerator::kDummyPrimaryPosition; - return false; + return true; } default: { @@ -613,6 +725,12 @@ void RMGVertexConfinement::SetSamplingModeString(std::string mode) { } catch (const std::bad_cast&) { return; } } +void RMGVertexConfinement::SetFirstSamplingVolumeTypeString(std::string type) { + try { + this->SetFirstSamplingVolumeType(RMGTools::ToEnum(type, "volume type")); + } catch (const std::bad_cast&) { return; } +} + void RMGVertexConfinement::AddPhysicalVolumeNameRegex(std::string name, std::string copy_nr) { if (copy_nr.empty()) copy_nr = ".*"; // for default arg from messenger. if (!fPhysicalVolumes.empty()) { @@ -627,22 +745,46 @@ void RMGVertexConfinement::AddGeometricalVolumeString(std::string solid) { GenericGeometricalSolidData data; data.solid_type = RMGTools::ToEnum(solid, "solid type"); fGeomVolumeData.push_back(data); + fLastSolidExcluded = false; +} + +void RMGVertexConfinement::AddExcludedGeometricalVolumeString(std::string solid) { + GenericGeometricalSolidData data; + data.solid_type = RMGTools::ToEnum(solid, "solid type"); + fExcludedGeomVolumeData.push_back(data); + fLastSolidExcluded = true; + + if (fSamplingMode != SamplingMode::kSubtractGeometrical) { + RMGLog::Out(RMGLog::fatal, "Cannot set excluded geometrical regions in any other sampling mode", + "' than kSubtractGeometrical"); + } } RMGVertexConfinement::GenericGeometricalSolidData& RMGVertexConfinement::SafeBack( std::optional solid_type) { - if (fGeomVolumeData.empty()) { - RMGLog::Out(RMGLog::fatal, "Must call /RMG/Generator/Confinement/Geometrical/AddSolid", + + + // checks - need to be performed both for GeomVolume and ExcludedGeomVolume + auto& volume_solids = fLastSolidExcluded ? fExcludedGeomVolumeSolids : fGeomVolumeSolids; + auto& volume_data = fLastSolidExcluded ? fExcludedGeomVolumeData : fGeomVolumeData; + auto command_name = fLastSolidExcluded ? "AddExcludedSolid" : "AddSolid"; + + if (volume_data.empty()) { + RMGLog::Out(RMGLog::fatal, "Must call /RMG/Generator/Confinement/Geometrical/", command_name, "' before setting any geometrical parameter value"); } - if (!fGeomVolumeSolids.empty()) { + + if (!volume_solids.empty()) { RMGLog::Out(RMGLog::fatal, "Solids for vertex confinement have already been initialized, no change possible!"); } - if (solid_type.has_value() && fGeomVolumeData.back().solid_type != solid_type) { + + + if (solid_type.has_value() && volume_data.back().solid_type != solid_type) { RMGLog::OutFormat(RMGLog::fatal, "Trying to modify non-{0} as {0}", solid_type.value()); } - return fGeomVolumeData.back(); + + return volume_data.back(); } void RMGVertexConfinement::BeginOfRunAction(const G4Run*) { @@ -663,6 +805,7 @@ void RMGVertexConfinement::EndOfRunAction(const G4Run* run) { void RMGVertexConfinement::DefineCommands() { + fMessengers.push_back(std::make_unique(this, "/RMG/Generator/Confinement/", "Commands for controlling primary confinement")); @@ -687,6 +830,15 @@ void RMGVertexConfinement::DefineCommands() { .SetStates(G4State_PreInit, G4State_Idle) .SetToBeBroadcasted(true); + fMessengers.back() + ->DeclareMethod("FirstSamplingVolume", &RMGVertexConfinement::SetFirstSamplingVolumeTypeString) + .SetGuidance("Select the type of volume which will be sampled first for intersections") + .SetParameterName("type", false) + .SetCandidates(RMGTools::GetCandidates()) + .SetStates(G4State_PreInit, G4State_Idle) + .SetToBeBroadcasted(true); + + fMessengers.back() ->DeclareProperty("MaxSamplingTrials", fMaxAttempts) .SetGuidance("Set maximum number of attempts for sampling primary positions in a volume") @@ -728,6 +880,15 @@ void RMGVertexConfinement::DefineCommands() { .SetStates(G4State_PreInit, G4State_Idle) .SetToBeBroadcasted(true); + fMessengers.back() + ->DeclareMethod("AddExcludedSolid", &RMGVertexConfinement::AddExcludedGeometricalVolumeString) + .SetGuidance("Add geometrical solid to exclude samples from") + .SetParameterName("solid", false) + .SetCandidates(RMGTools::GetCandidates()) + .SetStates(G4State_PreInit, G4State_Idle) + .SetToBeBroadcasted(true); + + // FIXME: see comment in .hh fMessengers.back() ->DeclareMethodWithUnit("CenterPositionX", "cm", &RMGVertexConfinement::SetGeomVolumeCenterX) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 2107bbe2..9acd21e0 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -2,9 +2,10 @@ file( GLOB _aux RELATIVE ${PROJECT_SOURCE_DIR} - macros/*.mac gdml/*.gdml gdml/*.xml) + macros/*.mac macros/*.json gdml/*.gdml gdml/*.xml *.py) # copy them to the build area + foreach(_file ${_aux}) configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) endforeach() @@ -27,3 +28,69 @@ foreach(_mac ${_macros}) set_tests_properties(confinement-vis/${_mac} PROPERTIES SKIP_REGULAR_EXPRESSION "couldn't open display") endforeach() + +# generate the GDML file +add_test(NAME confinment-ge/gen-gdml COMMAND ${PYTHONPATH} make_ge_array_gdml.py) +set_tests_properties(confinment-ge/gen-gdml PROPERTIES LABELS extra FIXTURES_SETUP + confine-gdml-fixture) + +# test on HPGe containment +add_test(NAME confinment-ge/gen-output COMMAND ${REMAGE_PYEXE} -g gdml/ge-array.gdml -w -o + test-confine.lh5 -- macros/test-ge-confine.mac) +set_tests_properties( + confinment-ge/gen-output PROPERTIES LABELS extra FIXTURES_SETUP confine-output-fixture + FIXTURES_REQUIRED confine-gdml-fixture) + +add_test(NAME confinment-ge/relative-fraction COMMAND ${PYTHONPATH} ./test_confinment_ge.py) + +set_tests_properties(confinment-ge/relative-fraction PROPERTIES LABELS extra FIXTURES_REQUIRED + confine-output-fixture) + +# test on lar intersection +add_test(NAME confinment-lar/intersection-gen-output + COMMAND ${REMAGE_PYEXE} -g gdml/ge-array.gdml -w -o test-confine-lar-in.lh5 -- + macros/lar-in.mac) +set_tests_properties( + confinment-lar/intersection-gen-output + PROPERTIES LABELS extra FIXTURES_SETUP confine-lar-in-output-fixture FIXTURES_REQUIRED + confine-gdml-fixture) + +add_test(NAME confinment-lar/intersection COMMAND ${PYTHONPATH} ./test_lar_intersection.py) + +set_tests_properties(confinment-lar/intersection PROPERTIES LABELS extra FIXTURES_REQUIRED + confine-lar-in-output-fixture) + +# generate subtraction output +add_test(NAME confinment-lar/subtraction-gen-output + COMMAND ${REMAGE_PYEXE} -g gdml/ge-array.gdml -w -o test-confine-lar-out.lh5 -- + macros/lar-out.mac) + +set_tests_properties( + confinment-lar/subtraction-gen-output + PROPERTIES LABELS extra FIXTURES_SETUP confine-lar-out-output-fixture FIXTURES_REQUIRED + confine-gdml-fixture) + +# analyse subtraction +add_test(NAME confinment-lar/subtraction COMMAND ${PYTHONPATH} ./test_lar_subtraction.py + test-confine-lar-out.lh5 lar-sub-check.output.pdf) + +set_tests_properties(confinment-lar/subtraction PROPERTIES LABELS extra FIXTURES_REQUIRED + confine-lar-out-output-fixture) + +# generate intersection and subtraction output +add_test(NAME confinment-lar/intersection-and-subtraction-gen-output + COMMAND ${REMAGE_PYEXE} -g gdml/ge-array.gdml -w -o test-confine-lar-in-out.lh5 -- + macros/lar-in-and-out.mac) + +set_tests_properties( + confinment-lar/intersection-and-subtraction-gen-output + PROPERTIES LABELS extra FIXTURES_SETUP confine-lar-in-out-output-fixture FIXTURES_REQUIRED + confine-gdml-fixture) + +# analyse intersection and subtraction +add_test(NAME confinment-lar/intersection-and-subtraction + COMMAND ${PYTHONPATH} ./test_lar_subtraction.py test-confine-lar-in-out.lh5 + lar-int-and-sub-check.output.pdf) + +set_tests_properties(confinment-lar/intersection-and-subtraction + PROPERTIES LABELS extra FIXTURES_REQUIRED confine-lar-ibn-out-output-fixture) diff --git a/tests/confinement/macros/B99000A.json b/tests/confinement/macros/B99000A.json new file mode 100644 index 00000000..d63577d3 --- /dev/null +++ b/tests/confinement/macros/B99000A.json @@ -0,0 +1,36 @@ +{ + "name": "B99000A", + "type": "bege", + "production": { + "manufacturer": "Test", + "enrichment": { + "val": 0.75, + "unc": 0.05 + } + }, + "geometry": { + "height_in_mm": 40.0, + "radius_in_mm": 35.0, + "groove": { + "depth_in_mm": 2.0, + "radius_in_mm": { + "outer": 12.0, + "inner": 10.0 + } + }, + "pp_contact": { + "radius_in_mm": 7.5, + "depth_in_mm": 0 + }, + "taper": { + "top": { + "angle_in_deg": 0.0, + "height_in_mm": 0.0 + }, + "bottom": { + "angle_in_deg": 0.0, + "height_in_mm": 0.0 + } + } + } +} diff --git a/tests/confinement/macros/C99000A.json b/tests/confinement/macros/C99000A.json new file mode 100644 index 00000000..c6b58b33 --- /dev/null +++ b/tests/confinement/macros/C99000A.json @@ -0,0 +1,44 @@ +{ + "name": "C99000A", + "type": "coax", + "production": { + "manufacturer": "Test", + "enrichment": { + "val": 0.75, + "unc": 0.05 + } + }, + "geometry": { + "height_in_mm": 80, + "radius_in_mm": 40, + "borehole": { + "radius_in_mm": 7, + "depth_in_mm": 70 + }, + "groove": { + "depth_in_mm": 2, + "radius_in_mm": { + "outer": 20, + "inner": 17 + } + }, + "pp_contact": { + "radius_in_mm": 17, + "depth_in_mm": 0 + }, + "taper": { + "top": { + "angle_in_deg": 45, + "height_in_mm": 5 + }, + "bottom": { + "angle_in_deg": 45, + "height_in_mm": 2 + }, + "borehole": { + "angle_in_deg": 0, + "height_in_mm": 0 + } + } + } +} diff --git a/tests/confinement/macros/P99000A.json b/tests/confinement/macros/P99000A.json new file mode 100644 index 00000000..ee684f2f --- /dev/null +++ b/tests/confinement/macros/P99000A.json @@ -0,0 +1,29 @@ +{ + "name": "P99000A", + "type": "ppc", + "production": { + "manufacturer": "Test", + "enrichment": { + "val": 0.75, + "unc": 0.05 + } + }, + "geometry": { + "height_in_mm": 45.0, + "radius_in_mm": 35.0, + "pp_contact": { + "radius_in_mm": 2, + "depth_in_mm": 2 + }, + "taper": { + "top": { + "angle_in_deg": 0, + "height_in_mm": 0 + }, + "bottom": { + "angle_in_deg": 30, + "height_in_mm": 10 + } + } + } +} diff --git a/tests/confinement/macros/V99000A.json b/tests/confinement/macros/V99000A.json new file mode 100644 index 00000000..ab1b0906 --- /dev/null +++ b/tests/confinement/macros/V99000A.json @@ -0,0 +1,44 @@ +{ + "name": "V99000A", + "type": "icpc", + "production": { + "manufacturer": "Test", + "enrichment": { + "val": 0.75, + "unc": 0.05 + } + }, + "geometry": { + "height_in_mm": 70, + "radius_in_mm": 35, + "borehole": { + "radius_in_mm": 5, + "depth_in_mm": 55 + }, + "groove": { + "depth_in_mm": 1, + "radius_in_mm": { + "outer": 10, + "inner": 9 + } + }, + "pp_contact": { + "radius_in_mm": 3, + "depth_in_mm": 2 + }, + "taper": { + "top": { + "angle_in_deg": 10, + "height_in_mm": 10 + }, + "bottom": { + "angle_in_deg": 0, + "height_in_mm": 0 + }, + "borehole": { + "angle_in_deg": 0, + "height_in_mm": 0 + } + } + } +} diff --git a/tests/confinement/macros/lar-in-and-out.mac b/tests/confinement/macros/lar-in-and-out.mac new file mode 100644 index 00000000..c5976acf --- /dev/null +++ b/tests/confinement/macros/lar-in-and-out.mac @@ -0,0 +1,25 @@ + +/RMG/Geometry/RegisterDetector Scintillator LAr 0 +/RMG/Output/NtuplePerDetector false +/run/initialize + +/RMG/Generator/Confine Volume + +/RMG/Generator/Confinement/SamplingMode SubtractGeometrical +/RMG/Generator/Confinement/Physical/AddVolume LAr + +/RMG/Generator/Confinement/Geometrical/AddSolid Cylinder +/RMG/Generator/Confinement/Geometrical/CenterPositionX 0 mm +/RMG/Generator/Confinement/Geometrical/CenterPositionY 0 mm +/RMG/Generator/Confinement/Geometrical/CenterPositionZ 50 mm +/RMG/Generator/Confinement/Geometrical/Cylinder/OuterRadius 170 mm +/RMG/Generator/Confinement/Geometrical/Cylinder/Height 400 mm + +# add lar volumes +/control/execute macros/lar-out-coordinates.mac + +/RMG/Generator/Select GPS +/gps/particle e- +/gps/energy 1 eV + +/run/beamOn 500000 diff --git a/tests/confinement/macros/lar-in.mac b/tests/confinement/macros/lar-in.mac new file mode 100644 index 00000000..5a6af878 --- /dev/null +++ b/tests/confinement/macros/lar-in.mac @@ -0,0 +1,18 @@ + +/RMG/Geometry/RegisterDetector Scintillator LAr 0 +/RMG/Output/NtuplePerDetector false +/run/initialize + +/RMG/Generator/Confine Volume + +/RMG/Generator/Confinement/SamplingMode IntersectPhysicalWithGeometrical +/RMG/Generator/Confinement/Physical/AddVolume LAr + +# add lar volumes +/control/execute macros/lar-in-coordinates.mac + +/RMG/Generator/Select GPS +/gps/particle e- +/gps/energy 1 eV + +/run/beamOn 500000 diff --git a/tests/confinement/macros/lar-out.mac b/tests/confinement/macros/lar-out.mac new file mode 100644 index 00000000..d46d77db --- /dev/null +++ b/tests/confinement/macros/lar-out.mac @@ -0,0 +1,18 @@ + +/RMG/Geometry/RegisterDetector Scintillator LAr 0 +/RMG/Output/NtuplePerDetector false +/run/initialize + +/RMG/Generator/Confine Volume + +/RMG/Generator/Confinement/SamplingMode SubtractGeometrical +/RMG/Generator/Confinement/Physical/AddVolume LAr + +# add lar volumes +/control/execute macros/lar-out-coordinates.mac + +/RMG/Generator/Select GPS +/gps/particle e- +/gps/energy 1 eV + +/run/beamOn 500000 diff --git a/tests/confinement/macros/test-ge-confine.mac b/tests/confinement/macros/test-ge-confine.mac new file mode 100644 index 00000000..53b8473b --- /dev/null +++ b/tests/confinement/macros/test-ge-confine.mac @@ -0,0 +1,17 @@ +/control/execute macros/detectors-fake.mac +/RMG/Output/NtuplePerDetector false + +/run/initialize + +/RMG/Generator/Confine Volume +/RMG/Generator/Confinement/Physical/AddVolume B.* +/RMG/Generator/Confinement/Physical/AddVolume C.* +/RMG/Generator/Confinement/Physical/AddVolume P.* +/RMG/Generator/Confinement/Physical/AddVolume V.* + + +/RMG/Generator/Select GPS +/gps/particle e- +/gps/energy 1 eV + +/run/beamOn 500000 diff --git a/tests/confinement/make_ge_array_gdml.py b/tests/confinement/make_ge_array_gdml.py new file mode 100644 index 00000000..067f63a5 --- /dev/null +++ b/tests/confinement/make_ge_array_gdml.py @@ -0,0 +1,136 @@ +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np +import pyg4ometry as pg4 +import pygeomtools as pytools +from legendhpges import make_hpge + +# read the configs +out_gdml = "gdml/ge-array.gdml" +det_macro = "macros/detectors-fake.mac" +config_dict = {} +heights = {} +for det_type in ["B", "P", "V", "C"]: + with Path.open(Path(f"macros/{det_type}99000A.json")) as file: + config_dict[det_type] = json.load(file) + heights[det_type] = config_dict[det_type]["geometry"]["height_in_mm"] + + +def add_hpge(lar, reg, angle, radius, z, idx, dtype): + x = radius * np.sin(np.deg2rad(angle)) + y = radius * np.cos(np.deg2rad(angle)) + logical_detector = make_hpge(config_dict[dtype], name=f"{dtype}{idx}", registry=reg) + logical_detector.pygeom_color_rgba = (0, 1, 1, 0.2) + physical_detector = pg4.geant4.PhysicalVolume( + [0, 0, 0], [x, y, z], logical_detector, f"{dtype}{idx}", lar, reg + ) + + physical_detector.pygeom_active_dector = pytools.RemageDetectorInfo( + "germanium", + idx, + config_dict[dtype], + ) + return idx + 1 + + +# construct geometry +reg = pg4.geant4.Registry() +ws = pg4.geant4.solid.Box("ws", 3000, 3000, 3000, reg, lunit="mm") +wl = pg4.geant4.LogicalVolume(ws, "G4_Galactic", "wl", reg) +wl.pygeom_color_rgba = (0.1, 1, 0.1, 0.5) + +reg.setWorld(wl) + + +# lar +lar_s = pg4.geant4.solid.Tubs( + "LAr_s", 0, 200, 800, 0, 2 * np.pi, registry=reg, lunit="mm" +) +lar_l = pg4.geant4.LogicalVolume(lar_s, "G4_lAr", "LAr_l", registry=reg) +lar_l.pygeom_color_rgba = (1, 0.1, 0, 0.2) +pg4.geant4.PhysicalVolume([0, 0, 0], [0, 0, 0], lar_l, "LAr", wl, registry=reg) + + +# hpge strings +string_radius = 85 +string_angles = [0, 80, 150, 220, 290] +offset = 200 + +n = 0 +det_info = [ + { + "angle": string_angles[0], + "radius": string_radius, + "detectors": ["V", "V", "V", "P"], + "offsets": [5, 10, 3, 5], + }, + { + "angle": string_angles[1], + "radius": string_radius, + "detectors": ["B", "B", "B", "B", "B", "B", "P"], + "offsets": [5, 4, 4, 3, 2, 5, 7], + }, + { + "angle": string_angles[2], + "radius": string_radius, + "detectors": ["V", "V", "B", "B"], + "offsets": [10, 8, 7, 2], + }, + { + "angle": string_angles[3], + "radius": string_radius, + "detectors": ["V", "V", "V"], + "offsets": [4, 3, 10], + }, + { + "angle": string_angles[4], + "radius": string_radius, + "detectors": ["B", "C", "C", "V"], + "offsets": [3, 7, 8, 5], + }, +] + + +n = 0 +lines = [] +for string in det_info: + angle = string["angle"] + radius = string["radius"] + z = offset + for d, o in zip(string["detectors"], string["offsets"]): + z -= heights[d] + n = add_hpge(lar_l, reg, angle, radius, z, n, d) + z -= o + + x = radius * np.sin(np.deg2rad(angle)) + y = radius * np.cos(np.deg2rad(angle)) + + lines.append("/RMG/Generator/Confinement/Geometrical/AddSolid Cylinder ") + lines.append(f"/RMG/Generator/Confinement/Geometrical/CenterPositionX {x} mm") + lines.append(f"/RMG/Generator/Confinement/Geometrical/CenterPositionY {y} mm ") + lines.append("/RMG/Generator/Confinement/Geometrical/CenterPositionZ 50 mm ") + lines.append("/RMG/Generator/Confinement/Geometrical/Cylinder/OuterRadius 44 mm ") + lines.append("/RMG/Generator/Confinement/Geometrical/Cylinder/Height 400 mm \n") + +lines_exclude = [line.replace("AddSolid", "AddExcludedSolid") for line in lines] + +# write the coordinates of the lar volumes +with Path("macros/lar-in-coordinates.mac").open("w") as f: + for line in lines: + f.write(line + "\n") + +with Path("macros/lar-out-coordinates.mac").open("w") as f: + for line in lines_exclude: + f.write(line + "\n") + +pytools.detectors.write_detector_auxvals(reg) +pytools.geometry.check_registry_sanity(reg, reg) + + +w = pg4.gdml.Writer() +w.addDetector(reg) +w.write(out_gdml) +pytools.detectors.generate_detector_macro(reg, det_macro) diff --git a/tests/confinement/test_confinment_ge.py b/tests/confinement/test_confinment_ge.py new file mode 100644 index 00000000..82feedd1 --- /dev/null +++ b/tests/confinement/test_confinment_ge.py @@ -0,0 +1,172 @@ +from __future__ import annotations + +import copy + +import awkward as ak +import legendhpges as hpges +import numpy as np +import pyg4ometry as pg4 +from lgdo import lh5 +from matplotlib import pyplot as plt +from pygeomtools import get_sensvol_metadata +from scipy import stats + +plt.rcParams["lines.linewidth"] = 1 +plt.rcParams["font.size"] = 12 + +gdml = "gdml/ge-array.gdml" +outfile = "test-confine.lh5" + +# get the geometry +reg = pg4.gdml.Reader(gdml).getRegistry() +reg_tmp = pg4.geant4.Registry() +detectors = list(reg.physicalVolumeDict.keys()) + +detectors = [det for det in detectors if det[0] in ["V", "P", "B", "C"]] + +det_map = { + det: { + "uint": int(det[1:]), + "pos": reg.physicalVolumeDict[det].position.eval(), + "hpge": hpges.make_hpge( + get_sensvol_metadata(reg, det), name=det, registry=reg_tmp + ), + } + for idx, det in enumerate(detectors) +} + +# append with the uint and the local positions +vertices = lh5.read_as("stp/vertices", outfile, "ak") + +uint = np.array(np.full_like(vertices.time, -1), dtype=int) +xlocal = np.array(1000 * vertices.xloc) +ylocal = np.array(1000 * vertices.yloc) +zlocal = np.array(1000 * vertices.zloc) + +positions = np.array( + np.transpose( + np.vstack([vertices.xloc * 1000, vertices.yloc * 1000, vertices.zloc * 1000]) + ) +) +for det in det_map: + local_positions = copy.copy(positions) + local_positions -= det_map[det]["pos"] + + is_inside = np.full(len(uint), False) + is_inside[uint == -1] = det_map[det]["hpge"].is_inside(local_positions[uint == -1]) + + uint[is_inside] = det_map[det]["uint"] + xlocal[is_inside] -= det_map[det]["pos"][0] + ylocal[is_inside] -= det_map[det]["pos"][1] + zlocal[is_inside] -= det_map[det]["pos"][2] + +vertices["uid"] = uint +vertices["xlocal"] = xlocal +vertices["ylocal"] = ylocal +vertices["zlocal"] = zlocal + +steps = lh5.read_as("stp/germanium", outfile, "ak") +steps = ak.unflatten(steps, ak.run_lengths(steps.evtid)) +uid = ak.fill_none(ak.firsts(steps.det_uid, axis=-1), -1) +evtid = ak.fill_none(ak.firsts(steps.evtid, axis=-1), -1) +hits = ak.Array({"evtid": evtid, "uid": uid}) + + +def make_plot(vert, hit): + fraction_vert = [] + fraction_vert_err = [] + fraction_hit = [] + expected_fraction = [] + names = [] + bad_uid = [] + + n_tot = len(vert) + + vol_tot = 0 + N = 0 + + for det in det_map: + vol_tot += float(det_map[det]["hpge"].volume.magnitude) + + n_sel_vert = len(vert[vert.uid == det_map[det]["uint"]]) + n_sel_hit = len(hit[hit.uid == det_map[det]["uint"]]) + + names.append(det) + bad_uid.append(det_map[det]["uint"]) + + fraction_vert.append(100 * n_sel_vert / n_tot) + fraction_vert_err.append(100 * np.sqrt(n_sel_vert) / n_tot) + fraction_hit.append(100 * n_sel_hit / n_tot) + expected_fraction.append(float(det_map[det]["hpge"].volume.magnitude)) + N += 1 + + expected_fraction = 100 * np.array(expected_fraction) / vol_tot + fraction_vert = np.array(fraction_vert) + fraction_vert_err = np.array(fraction_vert_err) + + # calculate the p-value + # technically the correct proobability distribution is a multinomial since Ntot is fixed + # I think that a product of Poisson will work fine since the overall poisson term is probably negligible, + # just results in N - 1 not N degrees of freedom. Wilkes theorem should surely hold.. + + test_stat = 2 * np.sum( + n_tot * expected_fraction / 100 + - (fraction_vert * n_tot / 100) + * (1 - np.log(fraction_vert / expected_fraction)) + ) + + # should follow a chi2 distribution with N -1 dof + + p = stats.chi2.sf(test_stat, N - 1) + sigma = stats.norm.ppf(1 - p) + + fig, ax = plt.subplots(2, 1, figsize=(12, 8), sharex=True) + + ax[0].errorbar( + np.arange(len(names)), + fraction_vert, + yerr=fraction_vert_err, + fmt=".", + label="Vertices", + ) + ax[0].errorbar(np.arange(len(names)) - 0.1, fraction_hit, fmt="x", label="Hits") + + ax[0].errorbar(np.arange(len(names)), expected_fraction, fmt=".", label="Expected") + ax[0].set_ylabel("Fraction of vertices [%]") + ax[0].set_xticks(np.arange(len(names)), names, rotation=90, fontsize=10) + ax[0].legend() + ax[0].grid() + + # residual + ax[1].errorbar( + np.arange(len(names)), + 100 * (fraction_vert - expected_fraction) / expected_fraction, + yerr=100 * fraction_vert_err / expected_fraction, + fmt=".", + label="Vertices", + ) + ax[1].errorbar( + np.arange(len(names)) - 0.1, + 100 * (fraction_hit - expected_fraction) / expected_fraction, + fmt="x", + label="Hits", + ) + ax[1].set_xlabel("Channel Index") + ax[1].set_ylabel("Relative difference [%]") + ax[1].set_xticks(np.arange(len(names)), names, rotation=90, fontsize=10) + ax[1].axhline(y=0, color="red") + ax[1].grid() + fig.suptitle(f"confinment check for HPGes, p = {100*p:.1e} %") + caption = "The fraction of the vertices found inside each HPGe. This is compared to the expectation which is that the number " + caption += "should be proportional to the volume of the HPGe. The top panel shows the fraction in each detector " + caption += r"while the lower panel shows the relative difference in % from the expectation." + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + + return p, sigma + + +p, sigma = make_plot(vertices, hits) +plt.savefig("confinement-ge.output.pdf") + +assert sigma < 5 diff --git a/tests/confinement/test_lar_intersection.py b/tests/confinement/test_lar_intersection.py new file mode 100644 index 00000000..8a1cdc92 --- /dev/null +++ b/tests/confinement/test_lar_intersection.py @@ -0,0 +1,241 @@ +from __future__ import annotations + +import legendhpges as hpges +import numpy as np +import pyg4ometry as pg4 +from lgdo import lh5 +from matplotlib import pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages +from pygeomtools.detectors import get_sensvol_metadata +from scipy import stats + +plt.rcParams["lines.linewidth"] = 1 +plt.rcParams["font.size"] = 12 + +# create the registry + +reg_tmp = pg4.geant4.Registry() +reg = pg4.gdml.Reader("gdml/ge-array.gdml").getRegistry() + +cylinder_height = 400 +cylinder_radius = 44 + +# extract the volume of each lar cylinder +vols = {} +detectors = list(reg.physicalVolumeDict.keys()) +detectors = [det for det in detectors if det[0] in ["V", "P", "B", "C"]] + + +for det in detectors: + vols[det[0]] = hpges.make_hpge( + get_sensvol_metadata(reg, f"{det}"), registry=reg_tmp, name=det + ).volume.magnitude + +# information on the detectors +string_radius = 85 +string_angles = [0, 80, 150, 220, 290] +det_info = [ + { + "angle": string_angles[0], + "radius": string_radius, + "detectors": ["V", "V", "V", "P"], + "offsets": [5, 10, 3, 5], + }, + { + "angle": string_angles[1], + "radius": string_radius, + "detectors": ["B", "B", "B", "B", "B", "B", "P"], + "offsets": [5, 4, 4, 3, 2, 5, 7], + }, + { + "angle": string_angles[2], + "radius": string_radius, + "detectors": ["V", "V", "B", "B"], + "offsets": [10, 8, 7, 2], + }, + { + "angle": string_angles[3], + "radius": string_radius, + "detectors": ["V", "V", "V"], + "offsets": [4, 3, 10], + }, + { + "angle": string_angles[4], + "radius": string_radius, + "detectors": ["B", "C", "C", "V"], + "offsets": [3, 7, 8, 5], + }, +] + +# get the coordinates of the lar cylinders +vol_tot = cylinder_height * np.pi * cylinder_radius**2 +x = [] +y = [] +det_vol = [] +for det in det_info: + angle = det["angle"] + radius = det["radius"] + x.append(string_radius * np.sin(np.deg2rad(angle))) + y.append(string_radius * np.cos(np.deg2rad(angle))) + det_vol.append(sum([vols[dtype] for dtype in det["detectors"]])) + +total_inside = 5 * vol_tot - sum(det_vol) + + +# read the output +outfile = "test-confine-lar-in.lh5" + +vertices = lh5.read_as("stp/vertices", outfile, "ak") + +# append which string the vertices are in + +strings = np.array(np.full_like(vertices.evtid, -1)) + +for idx, (xt, yt) in enumerate(zip(x, y)): + is_in = ( + np.sqrt((1000 * vertices.xloc - xt) ** 2 + (1000 * vertices.yloc - yt) ** 2) + < cylinder_radius + ) + strings[is_in] = idx + +vertices["strings"] = strings + +# make some plots of the coordinates + +with PdfPages("lar-in-check.output.pdf") as pdf: + # x-y + fig, ax = plt.subplots(figsize=(8, 6)) + + cmap = plt.cm.BuPu + cmap.set_under("w", 1) + hist = plt.hist2d( + np.array(1000 * vertices.xloc), + np.array(1000 * vertices.yloc), + bins=[100, 100], + range=[[-150, 150], [-150, 150]], + vmin=0.5, + cmap=cmap, + ) + ax.set_xlabel("x position [mm]") + ax.set_ylabel("y position [mm]") + ax.axis("equal") + + # set the caption + caption = "The x-y positions of the vertices for the simulations of the liquid argon inside imaginary mini-shrouds (per string). " + caption += "The vertices should be arranged in cylinders of radius 44 mm around each string. A higher density of points should be found near the edges." + + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + + pdf.savefig() + + # x-z + fig, ax = plt.subplots(1, 5, figsize=(10, 6), sharey=True) + for s in range(5): + min_x = min(1000 * vertices.xloc[vertices.strings == s]) + hist = ax[s].hist2d( + np.array(1000 * vertices.xloc[vertices.strings == s]), + np.array(1000 * vertices.zloc[vertices.strings == s]), + bins=[100, 100], + range=[[-150, 150], [-160, 260]], + vmin=0.5, + cmap=cmap, + ) + ax[s].set_xlabel("x position [mm]") + if s == 0: + ax[s].set_ylabel("z position [mm]") + plt.tight_layout() + + # set the caption + caption = "The x-z positions of the vertices for the simulations of the liquid argon inside imaginary mini-shrouds (per string). " + caption += "The vertices should be arranged in cylinders (center 50mm and height 400mm) and the regions with HPGe should be visible." + + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + + pdf.savefig() + + # now plot the fraction in each detector + fraction_vert = [] + fraction_vert_err = [] + names = [] + expected_fraction = [] + vert = vertices + n_tot = len(vert) + + for s, v in enumerate(det_vol): + n_sel_vert = len(vert[vert.strings == s]) + fraction_vert.append(100 * n_sel_vert / n_tot) + fraction_vert_err.append(100 * np.sqrt(n_sel_vert) / n_tot) + expected_fraction.append(100 * (vol_tot - v) / total_inside) + + fraction_vert = np.array(fraction_vert) + fraction_vert_err = np.array(fraction_vert_err) + expected_fraction = np.array(expected_fraction) + + # get the p-value + + test_stat = 2 * np.sum( + n_tot * expected_fraction / 100 + - (fraction_vert * n_tot / 100) + * (1 - np.log(fraction_vert / expected_fraction)) + ) + + # should follow a chi2 distribution with N -1 dof + + N = len(det_vol) + p = stats.chi2.sf(test_stat, N - 1) + sigma = stats.norm.ppf(1 - p) + + # make the plot + fig, ax = plt.subplots(2, 1, figsize=(12, 6), sharex=True) + + ax[0].errorbar( + np.arange(5), + fraction_vert, + yerr=fraction_vert_err, + fmt=".", + label="Vertices", + ) + ax[0].errorbar( + np.arange(5), + expected_fraction, + fmt="x", + label="Expected", + ) + + ax[0].set_ylabel("Fraction of vertices [%]") + ax[0].set_xticks( + np.arange(5), [f"string {i}" for i in range(5)], rotation=90, fontsize=13 + ) + ax[0].legend() + ax[0].grid() + + # make the residual + ax[1].errorbar( + np.arange(5), + 100 * (fraction_vert - expected_fraction) / expected_fraction, + yerr=100 * fraction_vert_err / expected_fraction, + fmt=".", + ) + + ax[1].set_ylabel("Relative Difference [%]") + + ax[0].set_xticks( + np.arange(5), [f"string {i}" for i in range(5)], rotation=90, fontsize=13 + ) + ax[1].axhline(y=0, color="red") + ax[1].grid() + fig.suptitle(f"Intersection check ({sigma:.1f} $\sigma$)") + + caption = "The fraction of the vertices found inside each cylinder. This is compared to the expectation which is that the number " + caption += "should be proportional to the volume of the cylinder minus the total HPGe volume. The top panel shows the fraction in each string " + caption += ( + r"while the lower panel shows the relative difference in % from the expectation" + ) + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + + pdf.savefig() + +assert sigma < 5 diff --git a/tests/confinement/test_lar_subtraction.py b/tests/confinement/test_lar_subtraction.py new file mode 100644 index 00000000..f0a3c04a --- /dev/null +++ b/tests/confinement/test_lar_subtraction.py @@ -0,0 +1,79 @@ +from __future__ import annotations + +import argparse +import sys + +import numpy as np +from lgdo import lh5 +from matplotlib import pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages + +plt.rcParams["lines.linewidth"] = 1 +plt.rcParams["font.size"] = 12 + +parser = argparse.ArgumentParser(description="Check on LAr subtraction") +parser.add_argument("outfile", type=str, help="File name to process") +parser.add_argument("plot", type=str, help="Where to save the plot") + +args = parser.parse_args() + +outfile = args.outfile + +vertices = lh5.read_as("stp/vertices", outfile, "ak") + +with PdfPages(args.plot) as pdf: + # x-y + fig, ax = plt.subplots(figsize=(8, 6)) + + cmap = plt.cm.BuPu + cmap.set_under("w", 1) + hist = plt.hist2d( + np.array(1000 * vertices.xloc), + np.array(1000 * vertices.yloc), + bins=[100, 100], + range=[[-300, 300], [-300, 300]], + vmin=0.5, + cmap=cmap, + ) + ax.set_xlabel("x position [mm]") + ax.set_ylabel("y position [mm]") + ax.axis("equal") + + # set the caption + caption = "The x-y positions of the vertices for the simulations of the liquid argon outside imaginary mini-shrouds (per string). " + caption += "The vertices should be excluded from cylinders of radius 44 mm around each string. A higher density of points should be found near the edges." + + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + + pdf.savefig() + + # x-z + fig, ax = plt.subplots(1, 1, figsize=(10, 6), sharey=True) + + min_x = min(1000 * vertices.xloc) + hist = ax.hist2d( + np.array(1000 * vertices.xloc), + np.array(1000 * vertices.zloc), + bins=[100, 100], + range=[[-250, 250], [-500, 500]], + vmin=0.5, + cmap=cmap, + ) + ax.set_xlabel("x position [mm]") + + plt.tight_layout() + + # set the caption + caption = "The x-z positions of the vertices for the simulations of the liquid argon inside imaginary mini-shrouds (per string). " + caption += ( + "The vertices should be excluded from cylinders (center 50mm and height 400mm)" + ) + + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + + pdf.savefig() + + +sys.exit(0)