From 8dd3883e9adeea362aa5b706a5fed24bb05c7ccf Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Thu, 9 Jan 2025 17:01:26 +0100 Subject: [PATCH 01/26] implement generic surface sampling --- include/RMGVertexConfinement.hh | 19 +- src/RMGVertexConfinement.cc | 223 +++++++++++++++--- tests/confinement/CMakeLists.txt | 5 +- .../test-surface-sampler-methods.cpp | 157 ++++++++++++ 4 files changed, 368 insertions(+), 36 deletions(-) create mode 100644 tests/confinement/test-surface-sampler-methods.cpp diff --git a/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index 48d5622..b6dfd29 100644 --- a/include/RMGVertexConfinement.hh +++ b/include/RMGVertexConfinement.hh @@ -69,6 +69,8 @@ class RMGVertexConfinement : public RMGVVertexGenerator { kGeometrical, kUnset, }; + + RMGVertexConfinement(); void BeginOfRunAction(const G4Run* run) override; @@ -96,12 +98,20 @@ class RMGVertexConfinement : public RMGVVertexGenerator { SampleableObject() = default; SampleableObject(const SampleableObject&) = default; SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r, G4ThreeVector t, G4VSolid* s, - bool cc = true); + bool ns = false, bool ss = false); // 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, + [[nodiscard]] bool Sample(G4ThreeVector& vertex, int max_attempts, bool force_containment_check, long int& n_trials) const; + [[nodiscard]] bool GenerateSurfacePoint(G4ThreeVector& vertex, int max_samples, + int n_max) const; + + // methods for the generic surface sampling + std::vector GetIntersections(const G4ThreeVector start, + const G4ThreeVector dir) const; + void GetDirection(G4ThreeVector& dir, G4ThreeVector& pos) const; + G4VPhysicalVolume* physical_volume = nullptr; G4VSolid* sampling_solid = nullptr; @@ -109,7 +119,9 @@ class RMGVertexConfinement : public RMGVVertexGenerator { G4ThreeVector translation; double volume = -1; double surface = -1; - bool containment_check = true; + bool surface_sample = false; + bool native_sample = false; + int max_num_intersections = -1; }; struct SampleableObjectCollection { @@ -179,6 +191,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator { bool fOnSurface = false; bool fForceContainmentCheck = false; bool fLastSolidExcluded = false; + int fMaxNumIntersections = -1; // counters used for the current run. long fTrials = 0; std::chrono::nanoseconds fVertexGenerationTime; diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index 91adc08..3d3fb8d 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -52,8 +52,8 @@ bool RMGVertexConfinement::fVolumesInitialized = false; // - physical volumes get always a bounding box assigned, but at later time // - purely geometrical volumes only have the sampling_solid member defined RMGVertexConfinement::SampleableObject::SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r, - G4ThreeVector t, G4VSolid* s, bool cc) - : rotation(r), translation(t), containment_check(cc) { + G4ThreeVector t, G4VSolid* s, bool ns, bool ss) + : rotation(r), translation(t), native_sample(ns), surface_sample(ss) { // at least one volume must be specified if (!v and !s) RMGLog::Out(RMGLog::error, "Invalid pointers given to constructor"); @@ -121,8 +121,125 @@ bool RMGVertexConfinement::SampleableObject::IsInside(const G4ThreeVector& verte return false; } + + +// get the intersections between the solid and an arbitrary line, defined by start and dir +std::vector RMGVertexConfinement::SampleableObject::GetIntersections(G4ThreeVector start, + G4ThreeVector dir) const { + + + // check if the physical volume exists + if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) + RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject where the physical volume is not set", + "this probably means you are trying to generically sample a geometrical volume which " + "instead", + "should be natively sampled"); + + auto solid = this->physical_volume->GetLogicalVolume()->GetSolid(); + + + G4double dist = 0; + std::vector intersections = {}; + G4ThreeVector int_point; + + int_point = start; + int counter = 0; + + dist = 0; + + while (dist < kInfinity) { + + dist = (counter % 2) == 0 ? solid->DistanceToIn(int_point, dir) + : solid->DistanceToOut(int_point, dir); + + int_point = int_point + dist * dir; + + if (dist < kInfinity) intersections.push_back(int_point); + counter++; + } + + if (intersections.size() % 2) throw; + + return intersections; +} + + +// generate a random position and direction for the surface sampling +void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, + G4ThreeVector& pos) const { + + + if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) + RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject where the physical volume is not set", + "this probably means you are trying to generically sample a geometrical volume which " + "instead", + "should be natively sampled"); + + // get the bounding radius + G4double bounding_radius = + physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); + + + // start on z-axis, pointing down. + pos = G4ThreeVector(0.0, 0.0, bounding_radius); + dir = G4ThreeVector(0.0, 0.0, -1.0); + + // push in rho direction by some impact parameter + G4double diskPhi = 2.0 * CLHEP::pi * G4UniformRand(); + G4double diskR = sqrt(G4UniformRand()) * bounding_radius; + + pos += G4ThreeVector(cos(diskPhi) * diskR, sin(diskPhi) * diskR, 0); + + // now rotate pos and dir by some random direction + G4double theta = acos(2.0 * G4UniformRand() - 1.0); + G4double phi = 2.0 * CLHEP::pi * G4UniformRand(); + + pos.rotateY(theta); + pos.rotateZ(phi); + dir.rotateY(theta); + dir.rotateZ(phi); +} + + +// generate with the generic sampler the actual surface point +bool RMGVertexConfinement::SampleableObject::GenerateSurfacePoint(G4ThreeVector& vertex, + int max_attempts, int n_max) const { + + int calls = 0; + G4ThreeVector dir; + G4ThreeVector pos; + + while (calls < max_attempts) { + + // chose a random line + this->GetDirection(dir, pos); + + // find the intersections + std::vector intersections = this->GetIntersections(pos, dir); + + if (intersections.size() == 0) continue; + + // pick one weighting by n_max and return it + int random_int = static_cast(n_max * G4UniformRand()); + + if (random_int <= intersections.size() - 1) { + vertex = intersections[random_int]; + return true; + } + calls++; + } + + RMGLog::Out(RMGLog::error, "Exceeded maximum number of allowed iterations (", max_attempts, + "), check that your surfaces are efficiently ", + "sampleable and try, eventually, to increase the threshold through the dedicated ", + "macro command. Returning dummy vertex"); + + return false; +} + + bool RMGVertexConfinement::SampleableObject::Sample(G4ThreeVector& vertex, int max_attempts, - bool sample_on_surface, bool force_containment_check, long int& n_trials) const { + bool force_containment_check, long int& n_trials) const { if (this->physical_volume) { RMGLog::OutFormatDev(RMGLog::debug, "Chosen random volume: '{}[{}]'", @@ -135,15 +252,33 @@ bool RMGVertexConfinement::SampleableObject::Sample(G4ThreeVector& vertex, int m int calls = 0; + // possible sampling strategies + // 1) native sampling + // 2) volume sampling with containment check + // 3) general surface sampling + - if (this->containment_check) { + if (this->native_sample) { vertex = this->translation + - this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, sample_on_surface); + this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, this->surface_sample); + 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"); + } + } else if (this->surface_sample) { + bool success = this->GenerateSurfacePoint(vertex, this->max_num_intersections, max_attempts); + + if (not success) return false; + + } else { + vertex = this->translation + this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, false); while (!this->IsInside(vertex) and calls++ < max_attempts) { n_trials++; - vertex = this->translation + - this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, sample_on_surface); + vertex = + this->translation + this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, false); RMGLog::OutDev(RMGLog::debug, "Vertex was not inside, new vertex: ", vertex / CLHEP::cm, " cm"); } if (calls >= max_attempts) { @@ -153,23 +288,14 @@ bool RMGVertexConfinement::SampleableObject::Sample(G4ThreeVector& vertex, int m "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(); for (const auto& o : data) { @@ -282,16 +408,39 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { // if there are no daughters one can avoid doing containment checks if (log_vol->GetNoDaughters() > 0) { RMGLog::OutDev(RMGLog::debug, "Has daughters, containment check needed"); - el.containment_check = true; + el.native_sample = false; + if (fOnSurface) + RMGLog::OutDev(RMGLog::error, + "Surface sampling is not implemented for volumes with daughters"); + } else { RMGLog::OutDev(RMGLog::debug, "Has no daughters, no containment check needed"); - el.containment_check = false; + el.native_sample = true; } + el.surface_sample = fOnSurface; } // if it's not sampleable, cannot perform native surface sampling else if (fOnSurface) { - RMGLog::Out(RMGLog::fatal, "Physical volume '", el.physical_volume->GetName(), - "' does not satisfy requirements for native surface sampling"); + + if (log_vol->GetNoDaughters() > 0) + RMGLog::OutDev(RMGLog::error, + "Surface sampling is not implemented for volumes with daughters"); + + RMGLog::Out(RMGLog::debug, "Physical volume '", el.physical_volume->GetName(), + "' does not satisfy requirements for native surface sampling so resort to general " + "surface sampler"); + + el.native_sample = false; + el.surface_sample = true; + el.max_num_intersections = fMaxNumIntersections; + + if (fMaxNumIntersections < 2) + RMGLog::Out(RMGLog::fatal, + " for generic surface sampling MaxNumIntersections, the maximum number of lines a ", + "line can intersect with the surface must be set with " + "/RMG/Generator/Confinement/MaxNumberOfIntersections", + "Note: this can be an overestimate."); + } // if we have a subtraction solid and the first one is supported for // sampling, use it but check for containment @@ -300,14 +449,15 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { RMGLog::OutDev(RMGLog::debug, "Is a subtraction solid, sampling from constituent solid with containment check"); el.sampling_solid = solid->GetConstituentSolid(0); - el.containment_check = true; + el.native_sample = false; + el.surface_sample = false; } // use bounding solid for all other cases else { RMGLog::OutDev(RMGLog::debug, "Is not sampleable natively, need a bounding box with ", "containment check"); - el.containment_check = true; - + el.native_sample = false; + el.surface_sample = false; // to get a guaranteed bounding solid we rely on G4VSolid::BoundingLimits() // the function, implemented for each G4 solid, calculates the dimensions // of a bounding box. NOTE: it's possible to obtain a radius through @@ -331,7 +481,7 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { el.sampling_solid = new G4Box(el.physical_volume->GetName() + "/RMGVertexConfinement::fBoundingBox", bb_x, bb_y, bb_z); - } // sampling_solid and containment_check must hold a valid value at this point + } // sampling_solid and native_sample and surface_sample must hold a valid value at this point // determine solid transformation w.r.t. world volume reference @@ -401,7 +551,7 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { for (size_t i = 1; i < trees.size(); ++i) { const auto& v = trees.at(i); new_obj_from_inspection.emplace_back(el.physical_volume, v.vol_global_rotation, - v.vol_global_translation, el.sampling_solid, el.containment_check); + v.vol_global_translation, el.sampling_solid, el.native_sample, el.surface_sample); } } @@ -446,7 +596,8 @@ void RMGVertexConfinement::InitializeGeometricalVolumes(bool use_excluded_volume d.solid_type); } - volume_solids.back().containment_check = false; + volume_solids.back().native_sample = true; + volume_solids.back().surface_sample = fOnSurface; RMGLog::Out(RMGLog::detail, "Added geometrical solid ", use_excluded_volumes ? "(excluded) " : " ", "of type '", @@ -557,8 +708,7 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { const auto& choice = choice_nonconst; // generate a candidate vertex - bool success = - choice.Sample(vertex, fMaxAttempts, fOnSurface, fForceContainmentCheck, fTrials); + bool success = choice.Sample(vertex, fMaxAttempts, fForceContainmentCheck, fTrials); if (!success) { RMGLog::Out(RMGLog::error, "Sampling unsuccessful return dummy vertex"); @@ -644,8 +794,7 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { const auto& choice = choice_nonconst; // generate a candidate vertex - bool success = - choice.Sample(vertex, fMaxAttempts, fOnSurface, fForceContainmentCheck, fTrials); + bool success = choice.Sample(vertex, fMaxAttempts, fForceContainmentCheck, fTrials); if (!success) { RMGLog::Out(RMGLog::error, "Sampling unsuccessful return dummy vertex"); @@ -701,7 +850,7 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { fOnSurface ? all_volumes.SurfaceWeightedRand() : all_volumes.VolumeWeightedRand(); // do the sampling - bool success = choice.Sample(vertex, fMaxAttempts, fOnSurface, fForceContainmentCheck, fTrials); + bool success = choice.Sample(vertex, fMaxAttempts, fForceContainmentCheck, fTrials); if (!success) { RMGLog::Out(RMGLog::error, "Sampling unsuccessful return dummy vertex"); @@ -847,6 +996,16 @@ void RMGVertexConfinement::DefineCommands() { .SetStates(G4State_PreInit, G4State_Idle) .SetToBeBroadcasted(true); + + fMessengers.back() + ->DeclareProperty("MaxNumberOfIntersections", fMaxNumIntersections) + .SetGuidance("Set maximum number of intersections of a line with the surface. Note: can be " + "set to an overestimate. ") + .SetParameterName("N", false) + .SetRange("N > 1") + .SetStates(G4State_PreInit, G4State_Idle) + .SetToBeBroadcasted(true); + fMessengers.back() ->DeclareProperty("ForceContainmentCheck", fForceContainmentCheck) .SetGuidance("If true (or omitted argument), perform a containment check even after sampling " diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 9acd21e..d1293dc 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -5,7 +5,6 @@ file( 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() @@ -94,3 +93,7 @@ add_test(NAME confinment-lar/intersection-and-subtraction set_tests_properties(confinment-lar/intersection-and-subtraction PROPERTIES LABELS extra FIXTURES_REQUIRED confine-lar-ibn-out-output-fixture) + +# test generic surface sampling methods +add_executable(test-surface-sampler-methods EXCLUDE_FROM_ALL test-surface-sampler-methods.cpp) +target_link_libraries(test-surface-sampler-methods PUBLIC remage Geant4) diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp new file mode 100644 index 0000000..2471684 --- /dev/null +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -0,0 +1,157 @@ +#include +#include +#include + +#include "G4Material.hh" +#include "G4NistManager.hh" +#include "G4PVPlacement.hh" +#include "G4SubtractionSolid.hh" +#include "G4SystemOfUnits.hh" +#include "G4ThreeVector.hh" +#include "G4Tubs.hh" +#include "G4VSolid.hh" +#include "Randomize.hh" + +#include "RMGVertexConfinement.hh" + +G4VPhysicalVolume* get_phy_volume(const G4VSolid* solid, std::string Material) { + G4NistManager* nist = G4NistManager::Instance(); + G4Material* Material = nist->FindOrBuildMaterial(mat); + G4LogicalVolume* logicalVolume = new G4LogicalVolume(solid, Material, "log_vol"); + + G4ThreeVector position = G4ThreeVector(0, 0, 0); // Origin + G4RotationMatrix* rotation = nullptr; + G4PVPlacement* physicalVolume = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0), logicalVolume, + "phy_vol", nullptr, false, 0); + + return physicalVolume; +} + + +int main(int argc, char* argv[]) { + // Check if exactly one argument is provided + std::string test_type = argv[1]; + + if (test_type == "test-intersections-basic") { + + // make a basic geometry - first just a G4Tubs + G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); + auto tubby_phy = get_phy_volume(tubby, "G4_SKIN_ICRP"); + RMGVertexConfinement::SampleableObject obj = RMGVertexConfinement::SampleableObject(tubby_phy, + G4RotationMatrix(), G4ThreeVector(), nullptr); + + + // shoot along x + std::vector ints = + obj.GetIntersections(G4ThreeVector(60 * mm, 0, 0), G4ThreeVector(-1, 0, 0)); + + if (ints.size() != 2) { + std::cout << "The number of intersections should be 2" << std::endl; + return 1; + } + // shoot along z + ints = obj.GetIntersections(G4ThreeVector(0 * mm, 0, 60 * mm), G4ThreeVector(0, 0, -1)); + + if (ints.size() != 2) { + std::cout << "The number of intersections should be 2" << std::endl; + return 1; + } + + // miss the object + ints = obj.GetIntersections(G4ThreeVector(60 * mm, 0, 0), G4ThreeVector(0, -1, 0)); + + if (ints.size() != 0) { + std::cout << "The number of intersections should be 0" << std::endl; + return 1; + } + + // now generate a bunch randomly + int i = 0; + while (i < 10000) { + G4ThreeVector dir; + G4ThreeVector pos; + obj.GetDirection(dir, pos); + + int ints = obj.GetIntersections(pos, dir).size(); + + if (ints != 0 and ints != 2) { + std::cout << "The number of intersections can only be 0 or 2 not " << ints << std::endl; + return 1; + } + i++; + } + + // tests passed + return 0; + + + } else if (test_type == "test-intersections-subtraction") + + G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); + G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); + G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); + auto subby_phy = get_phy_volume(subby, "G4_SKIN_ICRP"); + RMGVertexConfinement::SampleableObject obj = RMGVertexConfinement::SampleableObject(subby_phy, + G4RotationMatrix(), G4ThreeVector(), nullptr); + + + // shoot along x + + std::vector ints = + obj.GetIntersections(G4ThreeVector(60 * mm, 0, 0), G4ThreeVector(-1, 0, 0)); + + if (ints.size() != 4) { + std::cout << "The number of intersections should be 4" << std::endl; + return 1; + } + + // shoot along z + + std::vector ints = + obj.GetIntersections(G4ThreeVector(15 * mm, 0, 60 * mm), G4ThreeVector(0, 0, -1)); + + if (ints.size() != 2) { + std::cout << "The number of intersections should be 2" << std::endl; + return 1; + } + + int i = 0; + while (i < 10000) { + G4ThreeVector dir; + G4ThreeVector pos; + obj.GetDirection(dir, pos); + + int ints = obj.GetIntersections(pos, dir).size(); + + if (ints != 0 and ints != 2 and ints != 4) { + std::cout << "The number of intersections can only be 0, 2 or 4 not " << ints << std::endl; + return 1; + } + i++; + } + + else if (test_type == "test-containment") { + G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); + G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); + G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); + auto subby_phy = get_phy_volume(subby, "G4_SKIN_ICRP"); + RMGVertexConfinement::SampleableObject obj = RMGVertexConfinement::SampleableObject(subby_phy, + G4RotationMatrix(), G4ThreeVector(), nullptr); + + int i = 0; + while (i < 10000) { + G4ThreeVector pos; + bool success = obj.GenerateSurfacePoint(&pos, 4, 200); + + if (not obj.IsInside(&pos)) { + + std::cout << "the sampled position is not inside the solid" << std::endl; + return 1; + } + i++; + } + } + else { return 0; } + + return 0; +} From 8b0fa2e13506fd8c3bc70e74d6b0dca530f31e04 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Sat, 11 Jan 2025 22:25:13 +0100 Subject: [PATCH 02/26] [tests] add some basic tests for surface generator --- src/RMGVertexConfinement.cc | 10 +- tests/confinement/CMakeLists.txt | 52 ++++++++- tests/confinement/macros/vis-surface.mac | 60 +++++++++++ tests/confinement/make_basic_solids.py | 100 ++++++++++++++++++ .../test-surface-sampler-methods.cpp | 80 +++++++------- 5 files changed, 260 insertions(+), 42 deletions(-) create mode 100644 tests/confinement/macros/vis-surface.mac create mode 100644 tests/confinement/make_basic_solids.py diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index 3d3fb8d..c85a5b9 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -268,10 +268,16 @@ bool RMGVertexConfinement::SampleableObject::Sample(G4ThreeVector& vertex, int m vertex / CLHEP::cm, " cm"); } } else if (this->surface_sample) { - bool success = this->GenerateSurfacePoint(vertex, this->max_num_intersections, max_attempts); - + bool success = this->GenerateSurfacePoint(vertex, max_attempts, this->max_num_intersections); + vertex = this->translation + this->rotation * vertex; if (not success) return false; + if (force_containment_check && !this->IsInside(vertex)) { + RMGLog::OutDev(RMGLog::error, + "Generated vertex not inside sampling volumes (forced containment check): ", + vertex / CLHEP::cm, " cm"); + } + } else { vertex = this->translation + this->rotation * RMGGeneratorUtil::rand(this->sampling_solid, false); diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index d1293dc..cb9d2f1 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -96,4 +96,54 @@ set_tests_properties(confinment-lar/intersection-and-subtraction # test generic surface sampling methods add_executable(test-surface-sampler-methods EXCLUDE_FROM_ALL test-surface-sampler-methods.cpp) -target_link_libraries(test-surface-sampler-methods PUBLIC remage Geant4) +target_link_libraries(test-surface-sampler-methods PUBLIC remage BxDecay0::BxDecay0) + +add_test(NAME confinement/test-surface-sampler-methods-build + COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target + test-surface-sampler-methods) +set_tests_properties(confinement/test-surface-sampler-methods-build + PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) + +# add the tests +set(_surface_tests test-intersections-basic test-intersections-subtraction test-containment) + +foreach(_test ${_surface_tests}) + + add_test(NAME confinment-surface/${_test} COMMAND ./test-surface-sampler-methods ${_test}) + + set_tests_properties(confinment-surface/${_test} PROPERTIES LABELS extra FIXTURES_REQUIRED + test-surface-sampler-methods-fixture) + +endforeach() + +# now run some simulations and test the outputs + +# generate the GDML file +add_test(NAME confinment-surface/gen-gdml COMMAND ${PYTHONPATH} make_basic_solids.py) +set_tests_properties(confinment-surface/gen-gdml PROPERTIES LABELS extra FIXTURES_SETUP + confine-surface-gdml-fixture) + +# visualise + +set(_solids tubby box con uni sub trd) + +foreach(det ${_solids}) + + add_test( + NAME confinement-surface/gen-macro-${det} + COMMAND /bin/sh -c + "sed 's/\${det}/${det}/g' macros/vis-surface.mac > macros/vis-surface-${det}.mac") + + set_tests_properties( + confinement-surface/gen-macro-${det} + PROPERTIES LABELS extra FIXTURES_SETUP confine-surface-macro-${det}-fixture FIXTURES_REQUIRED + confine-surface-gdml-fixture) + + add_test(NAME confinment-surface/${det}-vis + COMMAND ${REMAGE_PYEXE} macros/vis-surface-${det}.mac macros/_vis-export.mac -g + gdml/simple-solids.gdml -w -l detail) + + set_tests_properties( + confinment-surface/${det}-vis PROPERTIES LABELS extra FIXTURES_REQUIRED + confine-surface-macro-${det}-fixture) +endforeach() diff --git a/tests/confinement/macros/vis-surface.mac b/tests/confinement/macros/vis-surface.mac new file mode 100644 index 0000000..4c6f942 --- /dev/null +++ b/tests/confinement/macros/vis-surface.mac @@ -0,0 +1,60 @@ +/RMG/Output/NtuplePerDetector false +/RMG/Geometry/RegisterDetector Germanium ${det} 001 + +/run/initialize + +/vis/open OGL +/vis/viewer/reset +/vis/viewer/set/globalLineWidthScale 1.5 +/vis/viewer/set/style wireframe + +# draw the detector +/vis/drawVolume + + +# set colors + +/vis/geometry/set/colour tubby 0 1 0.2 1 .6 +/vis/geometry/set/colour box 0 0.7 1 0.5 .6 +/vis/geometry/set/colour trd 0 0.2 0.2 0.5 .6 +/vis/geometry/set/colour con 0 0.0 0.2 0.8 .6 +/vis/geometry/set/colour sub 0 1 0.2 0.5 .6 +/vis/geometry/set/colour uni 0 1 0.2 0.1 .6 + +/vis/scene/add/text 10 0 30 cm 30 0 0 "This shows the primary position of events generated on the surface of ${det}." + +# set the viewpoint +/vis/scene/add/axes 0 0 0 0.1 m +/vis/viewer/set/viewpointVector 0.9 0.9 1 +/vis/viewer/set/upVector 0 0 1 + +# add trajectories and control the coloring +/vis/scene/add/trajectories smooth +/vis/scene/endOfEventAction accumulate +/vis/modeling/trajectories/create/drawByCharge +/vis/modeling/trajectories/drawByCharge-0/default/setDrawStepPts true +/vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 10 + +# set the background to white +/vis/viewer/set/background white +/vis/geometry/set/forceSolid +/vis/geometry/set/forceWireframe + + +# now generate the primaries +/RMG/Generator/Select GPS +/RMG/Generator/Confine Volume +/RMG/Generator/Confinement/SampleOnSurface true +/RMG/Generator/Confinement/ForceContainmentCheck true + + +/RMG/Generator/Confinement/Physical/AddVolume ${det} +/RMG/Generator/Confinement/MaxNumberOfIntersections 6 + +/gps/particle e- +/gps/ang/type iso + +/gps/energy 1 eV +/run/beamOn 1000 + +/control/alias export-fn vis-surface-${det}.output.pdf diff --git a/tests/confinement/make_basic_solids.py b/tests/confinement/make_basic_solids.py new file mode 100644 index 0000000..177ec2b --- /dev/null +++ b/tests/confinement/make_basic_solids.py @@ -0,0 +1,100 @@ +from __future__ import annotations + +import numpy as np +import pint +import pyg4ometry as pg4 + +u = pint.UnitRegistry() + +# read the configs + +out_gdml = "gdml/simple-solids.gdml" + + +reg = pg4.geant4.Registry() +ws = pg4.geant4.solid.Box("ws", 1000, 1000, 300, reg, lunit="mm") +wl = pg4.geant4.LogicalVolume(ws, "G4_lAr", "wl", reg) +wl.pygeom_color_rgba = (0.1, 1, 0.1, 0.5) +reg.setWorld(wl) + + +# tubby +tubby = pg4.geant4.solid.Tubs( + "tubby", + pRMin=0, + pRMax=60, + pDz=60, + pSPhi=0, + pDPhi=2 * np.pi, + registry=reg, + lunit="mm", +) +tubby_l = pg4.geant4.LogicalVolume(tubby, "G4_ADIPOSE_TISSUE_ICRP", "tubby", reg) +tubby_l.pygeom_color_rgba = (1, 0.2, 1, 1) +tubby_p = pg4.geant4.PhysicalVolume([0, 0, 0], [-200, 0, 0], tubby_l, "tubby", wl, reg) + +# box +box = pg4.geant4.solid.Box("box", pX=50, pY=50, pZ=100, registry=reg, lunit="mm") +box_l = pg4.geant4.LogicalVolume(box, "G4_CHLOROFORM", "box", reg) +box_l.pygeom_color_rgba = (0.7, 1, 0.5, 1) +box_p = pg4.geant4.PhysicalVolume([0, 0, 0], [-50, 0, 0], box_l, "box", wl, reg) + + +# trd +trd = pg4.geant4.solid.Trd( + "trd", pDx1=50, pDx2=10, pDy1=50, pDy2=10, pDz=100, lunit="mm", registry=reg +) +trd_l = pg4.geant4.LogicalVolume(trd, "G4_TESTIS_ICRP", "trd", reg) +trd_l.pygeom_color_rgba = (0.2, 0.2, 0.5, 1) +trd_p = pg4.geant4.PhysicalVolume([0, 0, 0], [100, 0, 0], trd_l, "trd", wl, reg) + + +# cone +con = pg4.geant4.solid.Cons( + "con", + pRmin2=0, + pRmax2=0, + pRmin1=0, + pRmax1=60, + pDz=100, + pSPhi=0, + pDPhi=2 * np.pi, + registry=reg, + lunit="mm", +) +con_l = pg4.geant4.LogicalVolume(con, "G4_BONE_CORTICAL_ICRP", "con", reg) +con_l.pygeom_color_rgba = (0.0, 0.2, 0.8, 1) +con_p = pg4.geant4.PhysicalVolume([0, 0, 0], [-200, 150, 0], con_l, "con", wl, reg) + +# subtraction +small_tubby = pg4.geant4.solid.Tubs( + "small_tubby", + pRMin=0, + pRMax=20, + pDz=60, + pSPhi=0, + pDPhi=2 * np.pi, + lunit="mm", + registry=reg, +) +sub = pg4.geant4.solid.Subtraction( + "sub", tubby, small_tubby, tra2=[[0, 0, 0], [0, 0, 0]], registry=reg +) +sub_l = pg4.geant4.LogicalVolume(sub, "G4_UREA", "sub", reg) +sub_l.pygeom_color_rgba = (1, 0.2, 0.5, 1) +sub_p = pg4.geant4.PhysicalVolume([0, 0, 0], [-50, 150, 0], sub_l, "sub", wl, reg) + +small_box = pg4.geant4.solid.Box( + "small_box", pX=20, pY=20, pZ=50, registry=reg, lunit="mm" +) +uni = pg4.geant4.solid.Union( + "uni", box, small_box, tra2=[[0, 0, 0], [0, 0, 60]], registry=reg +) +uni_l = pg4.geant4.LogicalVolume(uni, "G4_BLOOD_ICRP", "uni", reg) +uni_l.pygeom_color_rgba = (1, 0.2, 0.1, 1) +uni_p = pg4.geant4.PhysicalVolume([0, 0, 0], [100, 150, 0], uni_l, "uni", wl, reg) + + +w = pg4.gdml.Writer() +w.addDetector(reg) +w.write(out_gdml) diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index 2471684..2b45ce7 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -2,6 +2,7 @@ #include #include +#include "G4LogicalVolume.hh" #include "G4Material.hh" #include "G4NistManager.hh" #include "G4PVPlacement.hh" @@ -14,9 +15,9 @@ #include "RMGVertexConfinement.hh" -G4VPhysicalVolume* get_phy_volume(const G4VSolid* solid, std::string Material) { +G4VPhysicalVolume* get_phy_volume(G4VSolid* solid, std::string Material_string) { G4NistManager* nist = G4NistManager::Instance(); - G4Material* Material = nist->FindOrBuildMaterial(mat); + G4Material* Material = nist->FindOrBuildMaterial(Material_string); G4LogicalVolume* logicalVolume = new G4LogicalVolume(solid, Material, "log_vol"); G4ThreeVector position = G4ThreeVector(0, 0, 0); // Origin @@ -85,52 +86,51 @@ int main(int argc, char* argv[]) { return 0; - } else if (test_type == "test-intersections-subtraction") - + } else if (test_type == "test-intersections-subtraction") { G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); - G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); - G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); - auto subby_phy = get_phy_volume(subby, "G4_SKIN_ICRP"); - RMGVertexConfinement::SampleableObject obj = RMGVertexConfinement::SampleableObject(subby_phy, - G4RotationMatrix(), G4ThreeVector(), nullptr); + G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); + G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); + auto subby_phy = get_phy_volume(subby, "G4_SKIN_ICRP"); + RMGVertexConfinement::SampleableObject obj = RMGVertexConfinement::SampleableObject(subby_phy, + G4RotationMatrix(), G4ThreeVector(), nullptr); - // shoot along x + // shoot along x - std::vector ints = - obj.GetIntersections(G4ThreeVector(60 * mm, 0, 0), G4ThreeVector(-1, 0, 0)); + std::vector ints = + obj.GetIntersections(G4ThreeVector(60 * mm, 0, 0), G4ThreeVector(-1, 0, 0)); - if (ints.size() != 4) { - std::cout << "The number of intersections should be 4" << std::endl; - return 1; - } + if (ints.size() != 4) { + std::cout << "The number of intersections should be 4" << std::endl; + return 1; + } - // shoot along z + // shoot along z - std::vector ints = - obj.GetIntersections(G4ThreeVector(15 * mm, 0, 60 * mm), G4ThreeVector(0, 0, -1)); + ints = obj.GetIntersections(G4ThreeVector(15 * mm, 0, 60 * mm), G4ThreeVector(0, 0, -1)); - if (ints.size() != 2) { - std::cout << "The number of intersections should be 2" << std::endl; - return 1; - } + if (ints.size() != 2) { + std::cout << "The number of intersections should be 2" << std::endl; + return 1; + } - int i = 0; - while (i < 10000) { - G4ThreeVector dir; - G4ThreeVector pos; - obj.GetDirection(dir, pos); + int i = 0; + while (i < 10000) { + G4ThreeVector dir; + G4ThreeVector pos; + obj.GetDirection(dir, pos); - int ints = obj.GetIntersections(pos, dir).size(); + int num_ints = obj.GetIntersections(pos, dir).size(); - if (ints != 0 and ints != 2 and ints != 4) { - std::cout << "The number of intersections can only be 0, 2 or 4 not " << ints << std::endl; - return 1; + if (num_ints != 0 and num_ints != 2 and num_ints != 4) { + std::cout << "The number of intersections can only be 0, 2 or 4 not " << num_ints + << std::endl; + return 1; + } + i++; } - i++; - } - else if (test_type == "test-containment") { + } else if (test_type == "test-containment") { G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); @@ -141,17 +141,19 @@ int main(int argc, char* argv[]) { int i = 0; while (i < 10000) { G4ThreeVector pos; - bool success = obj.GenerateSurfacePoint(&pos, 4, 200); + bool success = obj.GenerateSurfacePoint(pos, 200, 4); - if (not obj.IsInside(&pos)) { + if (subby->Inside(pos) != EInside::kSurface) { - std::cout << "the sampled position is not inside the solid" << std::endl; + std::string side = (subby->Inside(pos) == EInside::kInside) ? "Inside" : "Outside"; + std::cout << "the sampled position is not inside the solid it is " << side << std::endl; return 1; } i++; } + } else { + return 0; } - else { return 0; } return 0; } From 207204692685927963c014a0f92638c3c4dd33ff Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Mon, 13 Jan 2025 01:19:09 +0100 Subject: [PATCH 03/26] [tests] add a more sophisticated test - but ... it fails for the unions --- tests/confinement/CMakeLists.txt | 49 +++- tests/confinement/macros/gen-surface.mac | 21 ++ tests/confinement/test_basic_surface.py | 351 +++++++++++++++++++++++ 3 files changed, 416 insertions(+), 5 deletions(-) create mode 100644 tests/confinement/macros/gen-surface.mac create mode 100644 tests/confinement/test_basic_surface.py diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index cb9d2f1..82c79ba 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -101,6 +101,7 @@ target_link_libraries(test-surface-sampler-methods PUBLIC remage BxDecay0::BxDec add_test(NAME confinement/test-surface-sampler-methods-build COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target test-surface-sampler-methods) + set_tests_properties(confinement/test-surface-sampler-methods-build PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) @@ -130,14 +131,14 @@ set(_solids tubby box con uni sub trd) foreach(det ${_solids}) add_test( - NAME confinement-surface/gen-macro-${det} + NAME confinement-surface/gen-vis-macro-${det} COMMAND /bin/sh -c "sed 's/\${det}/${det}/g' macros/vis-surface.mac > macros/vis-surface-${det}.mac") set_tests_properties( - confinement-surface/gen-macro-${det} - PROPERTIES LABELS extra FIXTURES_SETUP confine-surface-macro-${det}-fixture FIXTURES_REQUIRED - confine-surface-gdml-fixture) + confinement-surface/gen-vis-macro-${det} + PROPERTIES LABELS extra FIXTURES_SETUP confine-surface-vis-macro-${det}-fixture + FIXTURES_REQUIRED confine-surface-gdml-fixture) add_test(NAME confinment-surface/${det}-vis COMMAND ${REMAGE_PYEXE} macros/vis-surface-${det}.mac macros/_vis-export.mac -g @@ -145,5 +146,43 @@ foreach(det ${_solids}) set_tests_properties( confinment-surface/${det}-vis PROPERTIES LABELS extra FIXTURES_REQUIRED - confine-surface-macro-${det}-fixture) + confine-surface-vis-macro-${det}-fixture) endforeach() + +# run simulations and check + +foreach(det ${_solids}) + + add_test( + NAME confinement-surface/gen-macro-${det} + COMMAND /bin/sh -c + "sed 's/\${det}/${det}/g' macros/gen-surface.mac > macros/gen-surface-${det}.mac") + + set_tests_properties( + confinement-surface/gen-macro-${det} + PROPERTIES LABELS extra FIXTURES_SETUP confine-surface-macro-${det}-fixture FIXTURES_REQUIRED + confine-surface-gdml-fixture) + + add_test(NAME confinment-surface/${det}-gen + COMMAND ${REMAGE_PYEXE} macros/gen-surface-${det}.mac -g gdml/simple-solids.gdml -w -l + detail -o test-confine-surface-${det}-out.lh5) + + set_tests_properties( + confinment-surface/${det}-gen + PROPERTIES LABELS extra FIXTURES_SETUP confine-surface-${det} FIXTURES_REQUIRED + confine-surface-macro-${det}-fixture) +endforeach() + +# required fixtures +set(required_fixtures "") +foreach(det ${_solids}) + list(APPEND required_fixtures "confine-surface-${det}") +endforeach() + +# Join the fixtures into a semicolon-separated string +string(JOIN ";" required_fixtures_str ${required_fixtures}) + +add_test(NAME confinment-surface/relative-fractions COMMAND ${PYTHONPATH} ./test_basic_surface.py) + +set_tests_properties(confinment-surface/relative-fractions + PROPERTIES LABELS extra FIXTURES_REQUIRED required_fixtures) diff --git a/tests/confinement/macros/gen-surface.mac b/tests/confinement/macros/gen-surface.mac new file mode 100644 index 0000000..d376dd5 --- /dev/null +++ b/tests/confinement/macros/gen-surface.mac @@ -0,0 +1,21 @@ +/RMG/Output/NtuplePerDetector false +/RMG/Geometry/RegisterDetector Germanium ${det} 001 + +/run/initialize + + +# now generate the primaries +/RMG/Generator/Select GPS +/RMG/Generator/Confine Volume +/RMG/Generator/Confinement/SampleOnSurface true +/RMG/Generator/Confinement/ForceContainmentCheck true + + +/RMG/Generator/Confinement/Physical/AddVolume ${det} +/RMG/Generator/Confinement/MaxNumberOfIntersections 6 + +/gps/particle e- +/gps/ang/type iso + +/gps/energy 1 eV +/run/beamOn 500000 diff --git a/tests/confinement/test_basic_surface.py b/tests/confinement/test_basic_surface.py new file mode 100644 index 0000000..00febce --- /dev/null +++ b/tests/confinement/test_basic_surface.py @@ -0,0 +1,351 @@ +from __future__ import annotations + +import copy + +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 scipy import stats + +plt.rcParams["lines.linewidth"] = 1 +plt.rcParams["font.size"] = 12 + + +gdml = "gdml/simple-solids.gdml" + +reg = pg4.gdml.Reader(gdml).getRegistry() + + +def add_local_pos(vertices, pos): + xlocal = np.array(1000 * vertices.xloc) + ylocal = np.array(1000 * vertices.yloc) + zlocal = np.array(1000 * vertices.zloc) + + vertices["xlocal"] = xlocal - pos[0] + vertices["ylocal"] = ylocal - pos[1] + vertices["zlocal"] = zlocal - pos[2] + return vertices + + +tol = 1e-6 +select_sides = { + "tubby": { + "func": [ + lambda x, y, z: (abs(z - 30) < tol) & (np.sqrt(x**2 + y**2) < 60), # bottom + lambda x, y, z: (abs(np.sqrt(x**2 + y**2) - 60) < tol) + & (abs(z) < 30 - tol), # side + lambda x, y, z: (abs(z + 30) < tol) & (np.sqrt(x**2 + y**2) < 60), # top + ], + "area": [ + np.pi * 60**2, + 2 * np.pi * 60 * 60, + np.pi * 60**2, + ], + "order": [0, 1, 2], + "nice_name": "G4Tubs", + }, + "sub": { + "func": [ + lambda x, y, z: (abs(z + 30) < tol) + & (np.sqrt(x**2 + y**2) < 60) + & (np.sqrt(x**2 + y**2) > 20), # bottom + lambda x, y, z: (abs(z - 30) < tol) + & (np.sqrt(x**2 + y**2) < 60) + & (np.sqrt(x**2 + y**2) > 20), # top + lambda x, y, z: (abs(np.sqrt(x**2 + y**2) - 60) < tol) + & (abs(z) < 30 - tol), # outside + lambda x, y, z: (abs(np.sqrt(x**2 + y**2) - 20) < tol) + & (abs(z) < 30 - tol), # inside + ], + "area": [ + np.pi * (60**2 - 20**2), # bottom + np.pi * (60**2 - 20**2), # top + 2 * np.pi * 60 * 60, # outside + 2 * np.pi * 60 * 20, + ], # inside + "order": [1, 2, 3, 0], + "nice_name": "subtraction of two G4Tubs", + }, + "con": { + "func": [ + lambda x, y, z: (abs(z + 50) < tol) & (np.sqrt(x**2 + y**2) < 60), # bottom + lambda x, y, z: (abs(z + 100 * (np.sqrt(x**2 + y**2) / 60) - 50) < tol) + & (abs(z) < 50 - tol), # side, + ], + "area": [ + np.pi * 60**2, # bottom + np.pi * 60 * np.sqrt(60**2 + 100**2), # side + ], + "order": [1, 0], + "nice_name": "G4Cone", + }, + "box": { + "func": [ + lambda x, y, z: (abs(z - 50) < tol) & (abs(x) < 25) & (abs(y) < 25), # top + lambda x, y, z: (abs(x - 25) < tol) + & (abs(y) < 25) + & (abs(z) < 50 - tol), # sides + lambda x, y, z: (abs(y - 25) < tol) & (abs(x) < 25) & (abs(z) < 50 - tol), + lambda x, y, z: (abs(x + 25) < tol) & (abs(y) < 25) & (abs(z) < 50 - tol), + lambda x, y, z: (abs(y + 25) < tol) & (abs(x) < 25) & (abs(z) < 50 - tol), + lambda x, y, z: (abs(z + 50) < tol) & (abs(x) < 25) & (abs(y) < 25), + ], # bottom + "area": [50**2, 100 * 50, 100 * 50, 100 * 50, 100 * 50, 50**2], + "order": [0, 1, 2, 3, 4, 5], + "nice_name": "G4Box", + }, + "uni": { + "func": [ + lambda x, y, z: (abs(z + 50) < tol) + & (abs(x) < 25) + & (abs(y) < 25), # bottom + lambda x, y, z: (abs(x - 25) < tol) + & (abs(y) < 25) + & (abs(z) < 50 - tol), # sides + lambda x, y, z: (abs(y - 25) < tol) & (abs(x) < 25) & (abs(z) < 50 - tol), + lambda x, y, z: (abs(x + 25) < tol) & (abs(y) < 25) & (abs(z) < 50 - tol), + lambda x, y, z: (abs(y + 25) < tol) & (abs(x) < 25) & (abs(z) < 50 - tol), + lambda x, y, z: (abs(z - 50) < tol) + & (abs(x) < 25) + & (abs(y) < 25) + & ((abs(x) > 10) | (abs(y) > 10)), # top + lambda x, y, z: (abs(x - 10) < tol) + & (abs(y) < 10) + & (abs(z - 60) < 25), # small sides + lambda x, y, z: (abs(y - 10) < tol) & (abs(x) < 10) & (abs(z - 60) < 25), + lambda x, y, z: (abs(x + 10) < tol) & (abs(y) < 10) & (abs(z - 60) < 25), + lambda x, y, z: (abs(y + 10) < tol) & (abs(x) < 10) & (abs(z - 60) < 25), + lambda x, y, z: (abs(z - 85) < tol) & (abs(x) < 10) & (abs(y) < 10), + ], + "area": [ + 50 * 50.0, + 100 * 50.0, + 100 * 50.0, + 100 * 50.0, + 100 * 50.0, + 50 * 50.0 - 10 * 10.0, + 20 * 50.0, + 20 * 50.0, + 20 * 50.0, + 20 * 50.0, + 20 * 20.0, + ], + "order": [5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 0], + "nice_name": "Union of two G4Boxs", + }, + "trd": { + "func": [ + lambda x, y, z: (abs(z + 50) < tol) & (abs(x) < 25) & (abs(y) < 25), # base + lambda x, y, z: (abs(z - 50) < tol) & (abs(x) < 5) & (abs(y) < 5), # top + lambda x, y, z: (y > 5) + & (y < 25) + & (y > x) + & (y > -x) + & (abs(z) < 50 - tol), # sides + lambda x, y, z: (-y > 5) + & (-y < 25) + & (-y > x) + & (-y > -x) + & (abs(z) < 50 - tol), + lambda x, y, z: (x > 5) + & (x < 25) + & (x > y) + & (x > -y) + & (abs(z) < 50 - tol), + lambda x, y, z: (-x > 5) + & (-x < 25) + & (-x > y) + & (-x > -y) + & (abs(z) < 50 - tol), + ], + "area": [ + 50 * 50.0, + 10 * 10.0, + np.sqrt(20 * 20 + 100 * 100) * (10 + 50) / 2.0, + np.sqrt(20 * 20 + 100 * 100) * (10 + 50) / 2.0, + np.sqrt(20 * 20 + 100 * 100) * (10 + 50) / 2.0, + np.sqrt(20 * 20 + 100 * 100) * (10 + 50) / 2.0, + ], + "order": [1, 2, 4, 3, 5, 0], + "nice_name": "G4Trapezoid", + }, +} + +with PdfPages("simple-solids-surface.output.pdf") as pdf: + for dtype in ["tubby", "con", "sub", "uni", "box", "trd"]: + # get positions + pos = reg.physicalVolumeDict[dtype].position.eval() + + # read vertices and hits + outfile = f"test-confine-surface-{dtype}-out.lh5" + vertices = lh5.read_as("stp/vertices", outfile, "ak") + hits = lh5.read_as("stp/germanium", outfile, "ak") + + hits = add_local_pos(hits, pos) + vertices = add_local_pos(vertices, pos) + indices = np.array(np.full_like(vertices.time, -1)) + + # search for vertices being close to the sides + funcs = select_sides[dtype]["func"] + nice_name = select_sides[dtype]["nice_name"] + + for idx, fun in enumerate(funcs): + x = vertices.xlocal + y = vertices.ylocal + z = vertices.zlocal + is_close = fun( + np.array(vertices.xlocal), + np.array(vertices.ylocal), + np.array(vertices.zlocal), + ) + indices[is_close] = idx + + if len(indices[indices == -1]) > 0: + print(f"{dtype} has primaries not close to any side") + + vertices["idx"] = indices + vertices["rlocal"] = np.sqrt(vertices.xlocal**2 + vertices.ylocal**2) + n = 1000 + vert_small = vertices[0:n] + + # 3D plot + fig = plt.figure(constrained_layout=True, figsize=(8, 6)) + ax = fig.add_subplot(111, projection="3d", computed_zorder=False) + + order = copy.copy(select_sides[dtype]["order"]) + order.reverse() + + for idx in range(len(funcs)): + sel = vert_small.idx == order[idx] + scatter = ax.scatter( + vert_small[sel].xlocal, + vert_small[sel].ylocal, + vert_small[sel].zlocal, + zorder=idx, + s=5, + ) + + sel = vert_small.idx == -1 + scatter = ax.scatter( + vert_small[sel].xlocal, vert_small[sel].ylocal, vert_small[sel].zlocal, s=5 + ) + ax.view_init(elev=20, azim=30) + + ax.set_xlabel("x position [mm]") + ax.set_ylabel("y position [mm]") + ax.set_zlabel("z position [mm]") + + caption = f"The position of primary vertices for {nice_name} in 3D space. \n" + caption += "Primaries are grouped by the surface with each shown in a different color. " + + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + pdf.savefig() + + # 2D plots + fig, ax = plt.subplots(1, 3, figsize=(8, 6)) + for idx, var in enumerate( + [("xlocal", "ylocal"), ("rlocal", "zlocal"), ("xlocal", "zlocal")] + ): + for idx2 in range(len(funcs)): + sel = vert_small.idx == order[idx2] + ax[idx].scatter(vert_small[sel][var[0]], vert_small[sel][var[1]], s=2) + + sel = vert_small.idx == -1 + ax[idx].scatter(vert_small[sel][var[0]], vert_small[sel][var[1]], s=2) + + ax[idx].set_xlabel(f" {var[0]} [mm]") + ax[idx].set_ylabel(f" {var[1]} [mm]") + ax[idx].axis("equal") + plt.tight_layout() + + caption = f"The position of primary vertices for {nice_name} in 2D space for different projections. \n" + caption += "Primaries are grouped by the surface with each shown in a different color. " + + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + pdf.savefig() + + # statistical test + fraction_vert = [] + fraction_vert_err = [] + + names = [] + expected_fraction = [] + vert = vertices + n_tot = len(vertices) + tot = sum(select_sides[dtype]["area"]) + + for idx, a in enumerate(select_sides[dtype]["area"]): + n_sel_vert = len(vertices[vertices.idx == idx]) + 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 * a / tot) + + 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(select_sides[dtype]["area"]) + p = stats.chi2.sf(test_stat, N - 1) + sigma = stats.norm.ppf(1 - p) + + # make the plot + fig, ax = plt.subplots(2, 1, figsize=(8, 6), sharex=True) + + ax[0].errorbar( + np.arange(N), + fraction_vert, + yerr=fraction_vert_err, + fmt=".", + label="Vertices", + ) + ax[0].errorbar( + np.arange(N), + expected_fraction, + fmt="x", + label="Expected", + ) + + ax[0].set_ylabel("Fraction of vertices [%]") + ax[0].set_xticks( + np.arange(N), [f"face {i}" for i in range(N)], rotation=90, fontsize=13 + ) + ax[0].legend() + ax[0].grid() + + # make the residual + ax[1].errorbar( + np.arange(N), + 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(N), [f"face {i}" for i in range(N)], rotation=90, fontsize=13 + ) + ax[1].axhline(y=0, color="red") + ax[1].grid() + fig.suptitle(f"Surface uniformity test for {nice_name} ({sigma:.1f} $\\sigma$)") + + caption = "The fraction of the vertices found on each face of the shapes." + caption += "This should be proportional to the surface area. The top panel shows the fraction in each face " + 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() From 7d9151e1edf252d040900c77d6f1f48d515307d9 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Mon, 13 Jan 2025 12:26:32 +0100 Subject: [PATCH 04/26] [docs] add new macro commands --- docs/rmg-commands.md | 10 ++++++++++ tests/confinement/macros/vis-surface.mac | 2 +- 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index 79bb4e2..a9db9dd 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -958,6 +958,7 @@ Commands for controlling primary confinement * `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 +* `MaxNumberOfIntersections` – Set maximum number of intersections of a line with the surface. Note: can be set to an overestimate. * `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. ### `/RMG/Generator/Confinement/Reset` @@ -999,6 +1000,15 @@ Set maximum number of attempts for sampling primary positions in a volume * **Parameter type** – `i` * **Omittable** – `False` +### `/RMG/Generator/Confinement/MaxNumberOfIntersections` + +Set maximum number of intersections of a line with the surface. Note: can be set to an overestimate. + +* **Range of parameters** – `N > 1` +* **Parameter** – `N` + * **Parameter type** – `i` + * **Omittable** – `False` + ### `/RMG/Generator/Confinement/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. diff --git a/tests/confinement/macros/vis-surface.mac b/tests/confinement/macros/vis-surface.mac index 4c6f942..7d03be5 100644 --- a/tests/confinement/macros/vis-surface.mac +++ b/tests/confinement/macros/vis-surface.mac @@ -36,7 +36,7 @@ /vis/modeling/trajectories/drawByCharge-0/default/setStepPtsSize 10 # set the background to white -/vis/viewer/set/background white +#/vis/viewer/set/background white /vis/geometry/set/forceSolid /vis/geometry/set/forceWireframe From de50cb9e7d9185def42886cf1c8860b3674dd70b Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Mon, 13 Jan 2025 12:34:31 +0100 Subject: [PATCH 05/26] [tests] try to fix the tests in CI --- tests/confinement/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 82c79ba..bb772c4 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -95,6 +95,8 @@ set_tests_properties(confinment-lar/intersection-and-subtraction PROPERTIES LABELS extra FIXTURES_REQUIRED confine-lar-ibn-out-output-fixture) # test generic surface sampling methods +find_package(BxDecay0 ${RMG_BXDECAY0_MINIMUM_VERSION} REQUIRED COMPONENTS Geant4) + add_executable(test-surface-sampler-methods EXCLUDE_FROM_ALL test-surface-sampler-methods.cpp) target_link_libraries(test-surface-sampler-methods PUBLIC remage BxDecay0::BxDecay0) From 81c5d9bf1a9773a300266ccdf927a3726782457a Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Mon, 13 Jan 2025 12:39:50 +0100 Subject: [PATCH 06/26] some formatting fixes --- src/RMGVertexConfinement.cc | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index c85a5b9..b4a217a 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -158,7 +158,9 @@ std::vector RMGVertexConfinement::SampleableObject::GetIntersecti counter++; } - if (intersections.size() % 2) throw; + if (intersections.size() % 2) + RMGLog::OutDev(RMGLog::fatal, "Found ", intersections.size(), " intersections. However,", + "the number should always be divisible by two."); return intersections; } @@ -170,10 +172,9 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) - RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject where the physical volume is not set", - "this probably means you are trying to generically sample a geometrical volume which " - "instead", - "should be natively sampled"); + RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject ", + "where the physical volume is not set this probably means you are trying to generically", + "sample a geometrical volume which instead should be natively sampled"); // get the bounding radius G4double bounding_radius = @@ -185,10 +186,10 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, dir = G4ThreeVector(0.0, 0.0, -1.0); // push in rho direction by some impact parameter - G4double diskPhi = 2.0 * CLHEP::pi * G4UniformRand(); - G4double diskR = sqrt(G4UniformRand()) * bounding_radius; + G4double disk_phi = 2.0 * CLHEP::pi * G4UniformRand(); + G4double disk_r = sqrt(G4UniformRand()) * bounding_radius; - pos += G4ThreeVector(cos(diskPhi) * diskR, sin(diskPhi) * diskR, 0); + pos += G4ThreeVector(cos(disk_phi) * disk_r, sin(disk_phi) * disk_r, 0); // now rotate pos and dir by some random direction G4double theta = acos(2.0 * G4UniformRand() - 1.0); From 443dbb38eda615767c328d8bba7662291444f83a Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 00:33:37 +0100 Subject: [PATCH 07/26] [fix] adjust sampling sphere to have the correct center --- src/RMGVertexConfinement.cc | 18 +- tests/confinement/CMakeLists.txt | 17 +- tests/confinement/macros/gen-surface.mac | 2 +- tests/confinement/macros/vis-surface.mac | 1 - tests/confinement/make_basic_solids.py | 2 +- .../test-surface-sampler-methods.cpp | 294 ++++++++++++-- tests/confinement/test_basic_surface.py | 360 +++++++++--------- 7 files changed, 473 insertions(+), 221 deletions(-) diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index b4a217a..e39a604 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -130,10 +130,9 @@ std::vector RMGVertexConfinement::SampleableObject::GetIntersecti // check if the physical volume exists if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) - RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject where the physical volume is not set", - "this probably means you are trying to generically sample a geometrical volume which " - "instead", - "should be natively sampled"); + RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject ", + "where the physical volume is not set this probably means you are trying to generically ", + "sample a geometrical volume which instead should be natively sampled"); auto solid = this->physical_volume->GetLogicalVolume()->GetSolid(); @@ -172,14 +171,16 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) - RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject ", - "where the physical volume is not set this probably means you are trying to generically", + RMGLog::OutDev(RMGLog::fatal, "Cannot generate directions for a SampleableObject ", + "where the physical volume is not set this probably means you are trying to generically ", "sample a geometrical volume which instead should be natively sampled"); // get the bounding radius G4double bounding_radius = physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); + G4ThreeVector barycenter = + physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); // start on z-axis, pointing down. pos = G4ThreeVector(0.0, 0.0, bounding_radius); @@ -190,7 +191,7 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, G4double disk_r = sqrt(G4UniformRand()) * bounding_radius; pos += G4ThreeVector(cos(disk_phi) * disk_r, sin(disk_phi) * disk_r, 0); - + // now rotate pos and dir by some random direction G4double theta = acos(2.0 * G4UniformRand() - 1.0); G4double phi = 2.0 * CLHEP::pi * G4UniformRand(); @@ -199,6 +200,9 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, pos.rotateZ(phi); dir.rotateY(theta); dir.rotateZ(phi); + std::cout<<"barycenter "< #include #include +#include #include "G4LogicalVolume.hh" #include "G4Material.hh" #include "G4NistManager.hh" #include "G4PVPlacement.hh" #include "G4SubtractionSolid.hh" +#include "G4UnionSolid.hh" +#include "G4SubtractionSolid.hh" + + +#include "G4VUserDetectorConstruction.hh" #include "G4SystemOfUnits.hh" #include "G4ThreeVector.hh" #include "G4Tubs.hh" +#include "G4Box.hh" +#include "G4RunManager.hh" +#include "G4UImanager.hh" +#include "G4VisExecutive.hh" +#include "G4UIExecutive.hh" + +#include "G4NistManager.hh" +#include "G4Orb.hh" +#include "G4LogicalVolume.hh" +#include "G4PVPlacement.hh" +#include "G4Box.hh" +#include "G4ThreeVector.hh" +#include "G4Material.hh" +#include "G4Sphere.hh" +#include "G4SystemOfUnits.hh" + +#include "G4VisAttributes.hh" #include "G4VSolid.hh" #include "Randomize.hh" - +#include "QBBC.hh" #include "RMGVertexConfinement.hh" + + + +// class for some basic vis +class MyDetectorConstruction : public G4VUserDetectorConstruction { +public: + + G4double fRadius; + std::vector fPoints; + G4LogicalVolume * fLog; + G4ThreeVector fCenter; + + MyDetectorConstruction(G4ThreeVector center,G4double radius,std::vector points,G4LogicalVolume *solid_log){ + fCenter = center; + fRadius = radius; + fPoints= points; + fLog = solid_log; + }; + + G4VPhysicalVolume* Construct() override { + // Define materials + G4Material* vacuum = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"); + + // World volume + G4Box* worldBox = new G4Box("World", 200*mm,200*mm,200*mm); + G4LogicalVolume* worldLog = new G4LogicalVolume(worldBox, vacuum, "World"); + + // turn off vis + auto worldVisAttr = new G4VisAttributes(); + worldVisAttr->SetVisibility(true); + worldLog->SetVisAttributes(worldVisAttr); + + G4VPhysicalVolume* worldPhys = new G4PVPlacement(nullptr, {}, worldLog, "World", nullptr, false, 0); + + // Sampling sphere + G4Orb* orb = new G4Orb("MyOrb", fRadius); + G4LogicalVolume* orbLog = new G4LogicalVolume(orb, vacuum, "MyOrb"); + + + // place some solids + new G4PVPlacement(nullptr,fCenter, orbLog, "MyOrb", worldLog, false, 0); + new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0), fLog, "sampled_solid", worldLog, false, 0); + + // Visualization attributes for the orb + G4VisAttributes* orbVisAttr = new G4VisAttributes(G4Colour(1.0, 0.0, 0.0)); // Red + orbVisAttr->SetVisibility(true); + + orbLog->SetVisAttributes(orbVisAttr); + G4VisAttributes* VisAttr = new G4VisAttributes(G4Colour(0.0, 1.0, 0.0)); // green + VisAttr->SetVisibility(true); + + fLog->SetVisAttributes(VisAttr); + + // Create some G4ThreeVectors as points + + for (const auto& point :fPoints) { + G4Sphere* pointSphere = new G4Sphere("Point", 0, 1*mm, 0, 2 * M_PI, 0, M_PI); + G4LogicalVolume* pointLog = new G4LogicalVolume(pointSphere, vacuum, "Point"); + new G4PVPlacement(nullptr, point, pointLog, "Point", worldLog, false, 0); + + // Visualization attributes for points + + G4VisAttributes* pointVisAttr = new G4VisAttributes(G4Colour(0,0,1)); // Blue + pointVisAttr->SetVisibility(true); + + pointLog->SetVisAttributes(pointVisAttr); + + + } + + return worldPhys; + } +}; + + G4VPhysicalVolume* get_phy_volume(G4VSolid* solid, std::string Material_string) { G4NistManager* nist = G4NistManager::Instance(); G4Material* Material = nist->FindOrBuildMaterial(Material_string); @@ -33,14 +131,41 @@ int main(int argc, char* argv[]) { // Check if exactly one argument is provided std::string test_type = argv[1]; - if (test_type == "test-intersections-basic") { + // define the solids we need + std::map sampleables; - // make a basic geometry - first just a G4Tubs - G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); - auto tubby_phy = get_phy_volume(tubby, "G4_SKIN_ICRP"); - RMGVertexConfinement::SampleableObject obj = RMGVertexConfinement::SampleableObject(tubby_phy, + + // tubs + G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); + auto tubby_phy = get_phy_volume(tubby, "G4_SKIN_ICRP"); + RMGVertexConfinement::SampleableObject obj_tub = RMGVertexConfinement::SampleableObject(tubby_phy, + G4RotationMatrix(), G4ThreeVector(), nullptr); + + sampleables["tubs"]= obj_tub; + + G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); + G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); + auto subby_phy = get_phy_volume(subby, "G4_SKIN_ICRP"); + RMGVertexConfinement::SampleableObject obj_sub = RMGVertexConfinement::SampleableObject(subby_phy, + G4RotationMatrix(), G4ThreeVector(), nullptr); + + sampleables["sub"]=obj_sub; + + // box + G4Box* box = new G4Box("box", 25 * mm, 25 * mm, 50 * mm); + G4Box* small_box = new G4Box("small_box", 10 * mm, 10 * mm, 50 * mm); + + G4UnionSolid* solid = new G4UnionSolid("solid",box,small_box, nullptr, G4ThreeVector(0,0,60*mm)); + auto union_phy = get_phy_volume(solid, "G4_SKIN_ICRP"); + RMGVertexConfinement::SampleableObject obj_box = RMGVertexConfinement::SampleableObject(union_phy, G4RotationMatrix(), G4ThreeVector(), nullptr); + sampleables["uni"]=obj_box; + + if (test_type == "test-intersections-basic") { + + // make a basic geometry - first just a G4Tubs + auto obj = sampleables["tubs"]; // shoot along x std::vector ints = @@ -87,14 +212,8 @@ int main(int argc, char* argv[]) { } else if (test_type == "test-intersections-subtraction") { - G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); - G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); - G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); - auto subby_phy = get_phy_volume(subby, "G4_SKIN_ICRP"); - RMGVertexConfinement::SampleableObject obj = RMGVertexConfinement::SampleableObject(subby_phy, - G4RotationMatrix(), G4ThreeVector(), nullptr); - - + + auto obj = sampleables["sub"]; // shoot along x std::vector ints = @@ -130,28 +249,155 @@ int main(int argc, char* argv[]) { i++; } - } else if (test_type == "test-containment") { - G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); - G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); - G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); - auto subby_phy = get_phy_volume(subby, "G4_SKIN_ICRP"); - RMGVertexConfinement::SampleableObject obj = RMGVertexConfinement::SampleableObject(subby_phy, - G4RotationMatrix(), G4ThreeVector(), nullptr); + } + + else if (test_type == "test-intersections-union") { + + + auto obj = sampleables["uni"]; + + // shoot along x + + std::vector ints = + obj.GetIntersections(G4ThreeVector(90 * mm, 0, 0), G4ThreeVector(-1, 0, 0)); + if (ints.size() != 2) { + std::cout << "The number of intersections should be 2" << std::endl; + return 1; + } + + // shoot along z + ints = obj.GetIntersections(G4ThreeVector(15 * mm, 0, 95 * mm), G4ThreeVector(0, 0, -1)); + + if (ints.size() != 2) { + std::cout << "The number of intersections should be 2" << std::endl; + return 1; + } + + int i = 0; + while (i < 10000) { + G4ThreeVector dir; + G4ThreeVector pos; + obj.GetDirection(dir, pos); + + int num_ints = obj.GetIntersections(pos, dir).size(); + + if (num_ints != 0 and num_ints != 2 and num_ints!=4) { + std::cout << "The number of intersections can only be 0, 2 or 4 not " << num_ints + << std::endl; + return 1; + } + + i++; + } + + } + else if (test_type == "test-containment-union") { + auto obj = sampleables["uni"]; + int i = 0; while (i < 10000) { G4ThreeVector pos; bool success = obj.GenerateSurfacePoint(pos, 200, 4); - if (subby->Inside(pos) != EInside::kSurface) { + if (obj.physical_volume->GetLogicalVolume()->GetSolid()->Inside(pos) != EInside::kSurface) { - std::string side = (subby->Inside(pos) == EInside::kInside) ? "Inside" : "Outside"; + std::string side = (obj.physical_volume->GetLogicalVolume()->GetSolid()->Inside(pos) == EInside::kInside) ? "Inside" : "Outside"; + std::cout << "the sampled position is not inside the solid it is " << side << " "<GetLogicalVolume()->GetSolid()->Inside(pos) != EInside::kSurface) { + + std::string side = (obj.physical_volume->GetLogicalVolume()->GetSolid()->Inside(pos) == EInside::kInside) ? "Inside" : "Outside"; std::cout << "the sampled position is not inside the solid it is " << side << std::endl; return 1; } i++; } - } else { + } + else if (test_type == "test-points-union") { + + // get some points to plot + auto obj = sampleables["uni"]; + G4double radius = obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); + G4ThreeVector center = obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); + + G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); + + std::vector points; + int i = 0; + while (i < 1000) { + G4ThreeVector pos; + G4ThreeVector dir; + obj.GetDirection(dir,pos); + + points.push_back(pos); + if ((pos-center).mag()SetUserInitialization(new MyDetectorConstruction(center,radius,points,log)); + + auto physicsList = new QBBC; + physicsList->SetVerboseLevel(1); + runManager->SetUserInitialization(physicsList); + + // Initialize visualization + G4VisManager* visManager = new G4VisExecutive(); + visManager->Initialize(); + + // User interface + G4UIExecutive* ui = new G4UIExecutive(argc, argv); + + // Set up visualization commands + G4UImanager* UImanager = G4UImanager::GetUIpointer(); + UImanager->ApplyCommand("/run/initialize "); + + UImanager->ApplyCommand("/vis/open OGL "); + UImanager->ApplyCommand("/vis/verbose warnings"); + + UImanager->ApplyCommand("/vis/drawVolume"); + UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); + UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); + UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); + + UImanager->ApplyCommand("/vis/scene/list "); + UImanager->ApplyCommand("/vis/scene/add/axes 0 0 0 0.1 m"); + UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); + UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); + UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); + + //UImanager->ApplyCommand("/vis/viewer/set/background white"); + UImanager->ApplyCommand("/vis/ogl/export test-points-union.output.pdf"); + + + delete ui; + delete visManager; + delete runManager; + + return 0; + + + + } + else { return 0; } diff --git a/tests/confinement/test_basic_surface.py b/tests/confinement/test_basic_surface.py index 00febce..3b545d0 100644 --- a/tests/confinement/test_basic_surface.py +++ b/tests/confinement/test_basic_surface.py @@ -3,6 +3,8 @@ import copy import numpy as np +import argparse + import pyg4ometry as pg4 from lgdo import lh5 from matplotlib import pyplot as plt @@ -13,6 +15,12 @@ plt.rcParams["font.size"] = 12 +parser = argparse.ArgumentParser(description="Check on surface generation") +parser.add_argument("det", type=str, help="detector type") + +args = parser.parse_args() + + gdml = "gdml/simple-solids.gdml" reg = pg4.gdml.Reader(gdml).getRegistry() @@ -113,11 +121,11 @@ def add_local_pos(vertices, pos): & ((abs(x) > 10) | (abs(y) > 10)), # top lambda x, y, z: (abs(x - 10) < tol) & (abs(y) < 10) - & (abs(z - 60) < 25), # small sides - lambda x, y, z: (abs(y - 10) < tol) & (abs(x) < 10) & (abs(z - 60) < 25), - lambda x, y, z: (abs(x + 10) < tol) & (abs(y) < 10) & (abs(z - 60) < 25), - lambda x, y, z: (abs(y + 10) < tol) & (abs(x) < 10) & (abs(z - 60) < 25), - lambda x, y, z: (abs(z - 85) < tol) & (abs(x) < 10) & (abs(y) < 10), + & (abs(z - 75) < 25), # small sides + lambda x, y, z: (abs(y - 10) < tol) & (abs(x) < 10) & (abs(z - 75) < 25), + lambda x, y, z: (abs(x + 10) < tol) & (abs(y) < 10) & (abs(z - 75) < 25), + lambda x, y, z: (abs(y + 10) < tol) & (abs(x) < 10) & (abs(z - 75) < 25), + lambda x, y, z: (abs(z - 100) < tol) & (abs(x) < 10) & (abs(y) < 10), ], "area": [ 50 * 50.0, @@ -125,7 +133,7 @@ def add_local_pos(vertices, pos): 100 * 50.0, 100 * 50.0, 100 * 50.0, - 50 * 50.0 - 10 * 10.0, + 50 * 50.0 - 20 * 20.0, 20 * 50.0, 20 * 50.0, 20 * 50.0, @@ -172,180 +180,182 @@ def add_local_pos(vertices, pos): "nice_name": "G4Trapezoid", }, } - -with PdfPages("simple-solids-surface.output.pdf") as pdf: - for dtype in ["tubby", "con", "sub", "uni", "box", "trd"]: - # get positions - pos = reg.physicalVolumeDict[dtype].position.eval() - - # read vertices and hits - outfile = f"test-confine-surface-{dtype}-out.lh5" - vertices = lh5.read_as("stp/vertices", outfile, "ak") - hits = lh5.read_as("stp/germanium", outfile, "ak") - - hits = add_local_pos(hits, pos) - vertices = add_local_pos(vertices, pos) - indices = np.array(np.full_like(vertices.time, -1)) - - # search for vertices being close to the sides - funcs = select_sides[dtype]["func"] - nice_name = select_sides[dtype]["nice_name"] - - for idx, fun in enumerate(funcs): - x = vertices.xlocal - y = vertices.ylocal - z = vertices.zlocal - is_close = fun( - np.array(vertices.xlocal), - np.array(vertices.ylocal), - np.array(vertices.zlocal), - ) - indices[is_close] = idx - - if len(indices[indices == -1]) > 0: - print(f"{dtype} has primaries not close to any side") - - vertices["idx"] = indices - vertices["rlocal"] = np.sqrt(vertices.xlocal**2 + vertices.ylocal**2) - n = 1000 - vert_small = vertices[0:n] - - # 3D plot - fig = plt.figure(constrained_layout=True, figsize=(8, 6)) - ax = fig.add_subplot(111, projection="3d", computed_zorder=False) - - order = copy.copy(select_sides[dtype]["order"]) - order.reverse() - - for idx in range(len(funcs)): - sel = vert_small.idx == order[idx] - scatter = ax.scatter( - vert_small[sel].xlocal, - vert_small[sel].ylocal, - vert_small[sel].zlocal, - zorder=idx, - s=5, - ) - - sel = vert_small.idx == -1 - scatter = ax.scatter( - vert_small[sel].xlocal, vert_small[sel].ylocal, vert_small[sel].zlocal, s=5 - ) - ax.view_init(elev=20, azim=30) - - ax.set_xlabel("x position [mm]") - ax.set_ylabel("y position [mm]") - ax.set_zlabel("z position [mm]") - - caption = f"The position of primary vertices for {nice_name} in 3D space. \n" - caption += "Primaries are grouped by the surface with each shown in a different color. " - - plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) - plt.tight_layout(rect=[0, 0.12, 1, 1]) - pdf.savefig() - - # 2D plots - fig, ax = plt.subplots(1, 3, figsize=(8, 6)) - for idx, var in enumerate( - [("xlocal", "ylocal"), ("rlocal", "zlocal"), ("xlocal", "zlocal")] - ): - for idx2 in range(len(funcs)): - sel = vert_small.idx == order[idx2] - ax[idx].scatter(vert_small[sel][var[0]], vert_small[sel][var[1]], s=2) - - sel = vert_small.idx == -1 - ax[idx].scatter(vert_small[sel][var[0]], vert_small[sel][var[1]], s=2) - - ax[idx].set_xlabel(f" {var[0]} [mm]") - ax[idx].set_ylabel(f" {var[1]} [mm]") - ax[idx].axis("equal") - plt.tight_layout() - - caption = f"The position of primary vertices for {nice_name} in 2D space for different projections. \n" - caption += "Primaries are grouped by the surface with each shown in a different color. " - - plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) - plt.tight_layout(rect=[0, 0.12, 1, 1]) - pdf.savefig() - - # statistical test - fraction_vert = [] - fraction_vert_err = [] - - names = [] - expected_fraction = [] - vert = vertices - n_tot = len(vertices) - tot = sum(select_sides[dtype]["area"]) - - for idx, a in enumerate(select_sides[dtype]["area"]): - n_sel_vert = len(vertices[vertices.idx == idx]) - 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 * a / tot) - - 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)) +dtype = args.det + + +with PdfPages(f"simple-solids-surface-{dtype}.output.pdf") as pdf: + + # get positions + pos = reg.physicalVolumeDict[dtype].position.eval() + + # read vertices and hits + outfile = f"test-confine-surface-{dtype}-out.lh5" + vertices = lh5.read_as("stp/vertices", outfile, "ak") + hits = lh5.read_as("stp/germanium", outfile, "ak") + + hits = add_local_pos(hits, pos) + vertices = add_local_pos(vertices, pos) + indices = np.array(np.full_like(vertices.time, -1)) + + # search for vertices being close to the sides + funcs = select_sides[dtype]["func"] + nice_name = select_sides[dtype]["nice_name"] + + for idx, fun in enumerate(funcs): + x = vertices.xlocal + y = vertices.ylocal + z = vertices.zlocal + is_close = fun( + np.array(vertices.xlocal), + np.array(vertices.ylocal), + np.array(vertices.zlocal), ) + indices[is_close] = idx - # should follow a chi2 distribution with N -1 dof + if len(indices[indices == -1]) > 0: + print(f"{dtype} has primaries not close to any side") + assert False + vertices["idx"] = indices + vertices["rlocal"] = np.sqrt(vertices.xlocal**2 + vertices.ylocal**2) + n = 5000 + vert_small = vertices[0:n] - N = len(select_sides[dtype]["area"]) - p = stats.chi2.sf(test_stat, N - 1) - sigma = stats.norm.ppf(1 - p) + # 3D plot + fig = plt.figure(constrained_layout=True, figsize=(8, 6)) + ax = fig.add_subplot(111, projection="3d", computed_zorder=False) - # make the plot - fig, ax = plt.subplots(2, 1, figsize=(8, 6), sharex=True) + order = copy.copy(select_sides[dtype]["order"]) + order.reverse() - ax[0].errorbar( - np.arange(N), - fraction_vert, - yerr=fraction_vert_err, - fmt=".", - label="Vertices", - ) - ax[0].errorbar( - np.arange(N), - expected_fraction, - fmt="x", - label="Expected", - ) - - ax[0].set_ylabel("Fraction of vertices [%]") - ax[0].set_xticks( - np.arange(N), [f"face {i}" for i in range(N)], rotation=90, fontsize=13 - ) - ax[0].legend() - ax[0].grid() - - # make the residual - ax[1].errorbar( - np.arange(N), - 100 * (fraction_vert - expected_fraction) / expected_fraction, - yerr=100 * fraction_vert_err / expected_fraction, - fmt=".", + for idx in range(len(funcs)): + sel = vert_small.idx == order[idx] + scatter = ax.scatter( + vert_small[sel].xlocal, + vert_small[sel].ylocal, + vert_small[sel].zlocal, + zorder=idx, + s=5, ) - ax[1].set_ylabel("Relative Difference [%]") + sel = vert_small.idx == -1 + scatter = ax.scatter( + vert_small[sel].xlocal, vert_small[sel].ylocal, vert_small[sel].zlocal, s=5 + ) + ax.view_init(elev=20, azim=30) + + ax.set_xlabel("x position [mm]") + ax.set_ylabel("y position [mm]") + ax.set_zlabel("z position [mm]") + + caption = f"The position of primary vertices for {nice_name} in 3D space. \n" + caption += "Primaries are grouped by the surface with each shown in a different color. " + + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + pdf.savefig() + + # 2D plots + fig, ax = plt.subplots(1, 3, figsize=(8, 6)) + for idx, var in enumerate( + [("xlocal", "ylocal"), ("rlocal", "zlocal"), ("xlocal", "zlocal")] + ): + for idx2 in range(len(funcs)): + sel = vert_small.idx == order[idx2] + ax[idx].scatter(vert_small[sel][var[0]], vert_small[sel][var[1]], s=2) - ax[0].set_xticks( - np.arange(N), [f"face {i}" for i in range(N)], rotation=90, fontsize=13 - ) - ax[1].axhline(y=0, color="red") - ax[1].grid() - fig.suptitle(f"Surface uniformity test for {nice_name} ({sigma:.1f} $\\sigma$)") - - caption = "The fraction of the vertices found on each face of the shapes." - caption += "This should be proportional to the surface area. The top panel shows the fraction in each face " - 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() + sel = vert_small.idx == -1 + ax[idx].scatter(vert_small[sel][var[0]], vert_small[sel][var[1]], s=2) + + ax[idx].set_xlabel(f" {var[0]} [mm]") + ax[idx].set_ylabel(f" {var[1]} [mm]") + ax[idx].axis("equal") + plt.tight_layout() + + caption = f"The position of primary vertices for {nice_name} in 2D space for different projections. \n" + caption += "Primaries are grouped by the surface with each shown in a different color. " + + plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) + plt.tight_layout(rect=[0, 0.12, 1, 1]) + pdf.savefig() + + # statistical test + fraction_vert = [] + fraction_vert_err = [] + + names = [] + expected_fraction = [] + vert = vertices + n_tot = len(vertices) + tot = sum(select_sides[dtype]["area"]) + + for idx, a in enumerate(select_sides[dtype]["area"]): + n_sel_vert = len(vertices[vertices.idx == idx]) + 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 * a / tot) + + 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(select_sides[dtype]["area"]) + p = stats.chi2.sf(test_stat, N - 1) + sigma = stats.norm.ppf(1 - p) + + # make the plot + fig, ax = plt.subplots(2, 1, figsize=(8, 6), sharex=True) + + ax[0].errorbar( + np.arange(N), + fraction_vert, + yerr=fraction_vert_err, + fmt=".", + label="Vertices", + ) + ax[0].errorbar( + np.arange(N), + expected_fraction, + fmt="x", + label="Expected", + ) + + ax[0].set_ylabel("Fraction of vertices [%]") + ax[0].set_xticks( + np.arange(N), [f"face {i}" for i in range(N)], rotation=90, fontsize=13 + ) + ax[0].legend() + ax[0].grid() + + # make the residual + ax[1].errorbar( + np.arange(N), + 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(N), [f"face {i}" for i in range(N)], rotation=90, fontsize=13 + ) + ax[1].axhline(y=0, color="red") + ax[1].grid() + fig.suptitle(f"Surface uniformity test for {nice_name} ({sigma:.1f} $\\sigma$)") + + caption = "The fraction of the vertices found on each face of the shapes." + caption += "This should be proportional to the surface area. The top panel shows the fraction in each face " + 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() From 3478a24b33f93735b060956c0ce26319b973a2eb Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 00:45:18 +0100 Subject: [PATCH 08/26] formatting --- src/RMGVertexConfinement.cc | 7 +- tests/confinement/CMakeLists.txt | 9 +- .../test-surface-sampler-methods.cpp | 198 +++++++++--------- tests/confinement/test_basic_surface.py | 25 ++- 4 files changed, 119 insertions(+), 120 deletions(-) diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index e39a604..7c7e90c 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -179,7 +179,7 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, G4double bounding_radius = physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); - G4ThreeVector barycenter = + G4ThreeVector barycenter = physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); // start on z-axis, pointing down. @@ -191,7 +191,7 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, G4double disk_r = sqrt(G4UniformRand()) * bounding_radius; pos += G4ThreeVector(cos(disk_phi) * disk_r, sin(disk_phi) * disk_r, 0); - + // now rotate pos and dir by some random direction G4double theta = acos(2.0 * G4UniformRand() - 1.0); G4double phi = 2.0 * CLHEP::pi * G4UniformRand(); @@ -200,9 +200,8 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, pos.rotateZ(phi); dir.rotateY(theta); dir.rotateZ(phi); - std::cout<<"barycenter "< +#include #include #include -#include +#include "G4Box.hh" #include "G4LogicalVolume.hh" #include "G4Material.hh" #include "G4NistManager.hh" +#include "G4Orb.hh" #include "G4PVPlacement.hh" +#include "G4RunManager.hh" +#include "G4Sphere.hh" #include "G4SubtractionSolid.hh" -#include "G4UnionSolid.hh" -#include "G4SubtractionSolid.hh" - - -#include "G4VUserDetectorConstruction.hh" #include "G4SystemOfUnits.hh" #include "G4ThreeVector.hh" #include "G4Tubs.hh" -#include "G4Box.hh" -#include "G4RunManager.hh" -#include "G4UImanager.hh" -#include "G4VisExecutive.hh" #include "G4UIExecutive.hh" - -#include "G4NistManager.hh" -#include "G4Orb.hh" -#include "G4LogicalVolume.hh" -#include "G4PVPlacement.hh" -#include "G4Box.hh" -#include "G4ThreeVector.hh" -#include "G4Material.hh" -#include "G4Sphere.hh" -#include "G4SystemOfUnits.hh" - -#include "G4VisAttributes.hh" +#include "G4UImanager.hh" +#include "G4UnionSolid.hh" #include "G4VSolid.hh" +#include "G4VUserDetectorConstruction.hh" +#include "G4VisAttributes.hh" +#include "G4VisExecutive.hh" #include "Randomize.hh" -#include "QBBC.hh" -#include "RMGVertexConfinement.hh" +#include "RMGVertexConfinement.hh" +#include "QBBC.hh" // class for some basic vis class MyDetectorConstruction : public G4VUserDetectorConstruction { -public: + public: G4double fRadius; std::vector fPoints; - G4LogicalVolume * fLog; - G4ThreeVector fCenter; + G4LogicalVolume* fLog; + G4ThreeVector fCenter; - MyDetectorConstruction(G4ThreeVector center,G4double radius,std::vector points,G4LogicalVolume *solid_log){ + MyDetectorConstruction(G4ThreeVector center, G4double radius, std::vector points, + G4LogicalVolume* solid_log) { fCenter = center; fRadius = radius; - fPoints= points; + fPoints = points; fLog = solid_log; }; G4VPhysicalVolume* Construct() override { - // Define materials - G4Material* vacuum = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"); + // Define materials + G4Material* vacuum = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"); + + // World volume + G4Box* worldBox = new G4Box("World", 200 * mm, 200 * mm, 200 * mm); + G4LogicalVolume* worldLog = new G4LogicalVolume(worldBox, vacuum, "World"); - // World volume - G4Box* worldBox = new G4Box("World", 200*mm,200*mm,200*mm); - G4LogicalVolume* worldLog = new G4LogicalVolume(worldBox, vacuum, "World"); - - // turn off vis - auto worldVisAttr = new G4VisAttributes(); - worldVisAttr->SetVisibility(true); - worldLog->SetVisAttributes(worldVisAttr); - - G4VPhysicalVolume* worldPhys = new G4PVPlacement(nullptr, {}, worldLog, "World", nullptr, false, 0); - - // Sampling sphere - G4Orb* orb = new G4Orb("MyOrb", fRadius); - G4LogicalVolume* orbLog = new G4LogicalVolume(orb, vacuum, "MyOrb"); + // turn off vis + auto worldVisAttr = new G4VisAttributes(); + worldVisAttr->SetVisibility(true); + worldLog->SetVisAttributes(worldVisAttr); + G4VPhysicalVolume* worldPhys = + new G4PVPlacement(nullptr, {}, worldLog, "World", nullptr, false, 0); - // place some solids - new G4PVPlacement(nullptr,fCenter, orbLog, "MyOrb", worldLog, false, 0); - new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0), fLog, "sampled_solid", worldLog, false, 0); + // Sampling sphere + G4Orb* orb = new G4Orb("MyOrb", fRadius); + G4LogicalVolume* orbLog = new G4LogicalVolume(orb, vacuum, "MyOrb"); - // Visualization attributes for the orb - G4VisAttributes* orbVisAttr = new G4VisAttributes(G4Colour(1.0, 0.0, 0.0)); // Red - orbVisAttr->SetVisibility(true); - orbLog->SetVisAttributes(orbVisAttr); - G4VisAttributes* VisAttr = new G4VisAttributes(G4Colour(0.0, 1.0, 0.0)); // green - VisAttr->SetVisibility(true); + // place some solids + new G4PVPlacement(nullptr, fCenter, orbLog, "MyOrb", worldLog, false, 0); + new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0), fLog, "sampled_solid", worldLog, false, 0); - fLog->SetVisAttributes(VisAttr); + // Visualization attributes for the orb + G4VisAttributes* orbVisAttr = new G4VisAttributes(G4Colour(1.0, 0.0, 0.0)); // Red + orbVisAttr->SetVisibility(true); - // Create some G4ThreeVectors as points + orbLog->SetVisAttributes(orbVisAttr); + G4VisAttributes* VisAttr = new G4VisAttributes(G4Colour(0.0, 1.0, 0.0)); // green + VisAttr->SetVisibility(true); - for (const auto& point :fPoints) { - G4Sphere* pointSphere = new G4Sphere("Point", 0, 1*mm, 0, 2 * M_PI, 0, M_PI); - G4LogicalVolume* pointLog = new G4LogicalVolume(pointSphere, vacuum, "Point"); - new G4PVPlacement(nullptr, point, pointLog, "Point", worldLog, false, 0); + fLog->SetVisAttributes(VisAttr); - // Visualization attributes for points + // Create some G4ThreeVectors as points - G4VisAttributes* pointVisAttr = new G4VisAttributes(G4Colour(0,0,1)); // Blue - pointVisAttr->SetVisibility(true); + for (const auto& point : fPoints) { + G4Sphere* pointSphere = new G4Sphere("Point", 0, 1 * mm, 0, 2 * M_PI, 0, M_PI); + G4LogicalVolume* pointLog = new G4LogicalVolume(pointSphere, vacuum, "Point"); + new G4PVPlacement(nullptr, point, pointLog, "Point", worldLog, false, 0); - pointLog->SetVisAttributes(pointVisAttr); + // Visualization attributes for points + G4VisAttributes* pointVisAttr = new G4VisAttributes(G4Colour(0, 0, 1)); // Blue + pointVisAttr->SetVisibility(true); - } + pointLog->SetVisAttributes(pointVisAttr); + } - return worldPhys; + return worldPhys; } }; @@ -132,35 +120,36 @@ int main(int argc, char* argv[]) { std::string test_type = argv[1]; // define the solids we need - std::map sampleables; + std::map sampleables; // tubs G4Tubs* tubby = new G4Tubs("tubby", 0 * mm, 50 * mm, 50 * mm, 0, 360 * deg); auto tubby_phy = get_phy_volume(tubby, "G4_SKIN_ICRP"); RMGVertexConfinement::SampleableObject obj_tub = RMGVertexConfinement::SampleableObject(tubby_phy, - G4RotationMatrix(), G4ThreeVector(), nullptr); + G4RotationMatrix(), G4ThreeVector(), nullptr); - sampleables["tubs"]= obj_tub; + sampleables["tubs"] = obj_tub; G4Tubs* small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); G4SubtractionSolid* subby = new G4SubtractionSolid("subby", tubby, small_tubby); auto subby_phy = get_phy_volume(subby, "G4_SKIN_ICRP"); RMGVertexConfinement::SampleableObject obj_sub = RMGVertexConfinement::SampleableObject(subby_phy, - G4RotationMatrix(), G4ThreeVector(), nullptr); + G4RotationMatrix(), G4ThreeVector(), nullptr); - sampleables["sub"]=obj_sub; + sampleables["sub"] = obj_sub; // box G4Box* box = new G4Box("box", 25 * mm, 25 * mm, 50 * mm); G4Box* small_box = new G4Box("small_box", 10 * mm, 10 * mm, 50 * mm); - G4UnionSolid* solid = new G4UnionSolid("solid",box,small_box, nullptr, G4ThreeVector(0,0,60*mm)); + G4UnionSolid* solid = + new G4UnionSolid("solid", box, small_box, nullptr, G4ThreeVector(0, 0, 60 * mm)); auto union_phy = get_phy_volume(solid, "G4_SKIN_ICRP"); RMGVertexConfinement::SampleableObject obj_box = RMGVertexConfinement::SampleableObject(union_phy, - G4RotationMatrix(), G4ThreeVector(), nullptr); + G4RotationMatrix(), G4ThreeVector(), nullptr); - sampleables["uni"]=obj_box; + sampleables["uni"] = obj_box; if (test_type == "test-intersections-basic") { @@ -212,7 +201,7 @@ int main(int argc, char* argv[]) { } else if (test_type == "test-intersections-subtraction") { - + auto obj = sampleables["sub"]; // shoot along x @@ -249,8 +238,8 @@ int main(int argc, char* argv[]) { i++; } - } - + } + else if (test_type == "test-intersections-union") { @@ -282,7 +271,7 @@ int main(int argc, char* argv[]) { int num_ints = obj.GetIntersections(pos, dir).size(); - if (num_ints != 0 and num_ints != 2 and num_ints!=4) { + if (num_ints != 0 and num_ints != 2 and num_ints != 4) { std::cout << "The number of intersections can only be 0, 2 or 4 not " << num_ints << std::endl; return 1; @@ -291,10 +280,9 @@ int main(int argc, char* argv[]) { i++; } - } - else if (test_type == "test-containment-union") { + } else if (test_type == "test-containment-union") { auto obj = sampleables["uni"]; - + int i = 0; while (i < 10000) { G4ThreeVector pos; @@ -302,14 +290,17 @@ int main(int argc, char* argv[]) { if (obj.physical_volume->GetLogicalVolume()->GetSolid()->Inside(pos) != EInside::kSurface) { - std::string side = (obj.physical_volume->GetLogicalVolume()->GetSolid()->Inside(pos) == EInside::kInside) ? "Inside" : "Outside"; - std::cout << "the sampled position is not inside the solid it is " << side << " "<GetLogicalVolume()->GetSolid()->Inside(pos) == EInside::kInside) + ? "Inside" + : "Outside"; + std::cout << "the sampled position is not inside the solid it is " << side << " " << pos + << std::endl; return 1; } i++; } - } - else if (test_type == "test-containment-subtraction") { + } else if (test_type == "test-containment-subtraction") { auto obj = sampleables["sub"]; int i = 0; @@ -319,19 +310,23 @@ int main(int argc, char* argv[]) { if (obj.physical_volume->GetLogicalVolume()->GetSolid()->Inside(pos) != EInside::kSurface) { - std::string side = (obj.physical_volume->GetLogicalVolume()->GetSolid()->Inside(pos) == EInside::kInside) ? "Inside" : "Outside"; + std::string side = + (obj.physical_volume->GetLogicalVolume()->GetSolid()->Inside(pos) == EInside::kInside) + ? "Inside" + : "Outside"; std::cout << "the sampled position is not inside the solid it is " << side << std::endl; return 1; } i++; } - } - else if (test_type == "test-points-union") { + } else if (test_type == "test-points-union") { // get some points to plot auto obj = sampleables["uni"]; - G4double radius = obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); - G4ThreeVector center = obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); + G4double radius = + obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); + G4ThreeVector center = + obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); @@ -340,12 +335,11 @@ int main(int argc, char* argv[]) { while (i < 1000) { G4ThreeVector pos; G4ThreeVector dir; - obj.GetDirection(dir,pos); + obj.GetDirection(dir, pos); points.push_back(pos); - if ((pos-center).mag()SetUserInitialization(new MyDetectorConstruction(center,radius,points,log)); - + runManager->SetUserInitialization(new MyDetectorConstruction(center, radius, points, log)); + auto physicsList = new QBBC; physicsList->SetVerboseLevel(1); runManager->SetUserInitialization(physicsList); @@ -369,7 +363,7 @@ int main(int argc, char* argv[]) { // Set up visualization commands G4UImanager* UImanager = G4UImanager::GetUIpointer(); UImanager->ApplyCommand("/run/initialize "); - + UImanager->ApplyCommand("/vis/open OGL "); UImanager->ApplyCommand("/vis/verbose warnings"); @@ -377,17 +371,17 @@ int main(int argc, char* argv[]) { UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); - + UImanager->ApplyCommand("/vis/scene/list "); UImanager->ApplyCommand("/vis/scene/add/axes 0 0 0 0.1 m"); UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); - //UImanager->ApplyCommand("/vis/viewer/set/background white"); + // UImanager->ApplyCommand("/vis/viewer/set/background white"); UImanager->ApplyCommand("/vis/ogl/export test-points-union.output.pdf"); - + delete ui; delete visManager; delete runManager; @@ -395,9 +389,7 @@ int main(int argc, char* argv[]) { return 0; - - } - else { + } else { return 0; } diff --git a/tests/confinement/test_basic_surface.py b/tests/confinement/test_basic_surface.py index 3b545d0..f80f1c8 100644 --- a/tests/confinement/test_basic_surface.py +++ b/tests/confinement/test_basic_surface.py @@ -1,10 +1,9 @@ from __future__ import annotations +import argparse import copy import numpy as np -import argparse - import pyg4ometry as pg4 from lgdo import lh5 from matplotlib import pyplot as plt @@ -180,11 +179,10 @@ def add_local_pos(vertices, pos): "nice_name": "G4Trapezoid", }, } -dtype = args.det +dtype = args.det with PdfPages(f"simple-solids-surface-{dtype}.output.pdf") as pdf: - # get positions pos = reg.physicalVolumeDict[dtype].position.eval() @@ -213,8 +211,9 @@ def add_local_pos(vertices, pos): indices[is_close] = idx if len(indices[indices == -1]) > 0: - print(f"{dtype} has primaries not close to any side") - assert False + msg = f"{dtype} has primaries not close to any side" + raise ValueError(msg) + vertices["idx"] = indices vertices["rlocal"] = np.sqrt(vertices.xlocal**2 + vertices.ylocal**2) n = 5000 @@ -248,7 +247,9 @@ def add_local_pos(vertices, pos): ax.set_zlabel("z position [mm]") caption = f"The position of primary vertices for {nice_name} in 3D space. \n" - caption += "Primaries are grouped by the surface with each shown in a different color. " + caption += ( + "Primaries are grouped by the surface with each shown in a different color. " + ) plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) plt.tight_layout(rect=[0, 0.12, 1, 1]) @@ -272,7 +273,9 @@ def add_local_pos(vertices, pos): plt.tight_layout() caption = f"The position of primary vertices for {nice_name} in 2D space for different projections. \n" - caption += "Primaries are grouped by the surface with each shown in a different color. " + caption += ( + "Primaries are grouped by the surface with each shown in a different color. " + ) plt.figtext(0.1, 0.06, caption, wrap=True, ha="left", fontsize=11) plt.tight_layout(rect=[0, 0.12, 1, 1]) @@ -355,7 +358,11 @@ def add_local_pos(vertices, pos): caption = "The fraction of the vertices found on each face of the shapes." caption += "This should be proportional to the surface area. The top panel shows the fraction in each face " - caption += r"while the lower panel shows the relative difference in % from the expectation" + 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 From 12bb478cef310c13ecedbda2cfad0ffb0f80b5af Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 01:18:53 +0100 Subject: [PATCH 09/26] [tests] fix the test of the bounding sphere --- tests/confinement/CMakeLists.txt | 2 +- .../test-surface-sampler-methods.cpp | 146 ++++++++++-------- 2 files changed, 83 insertions(+), 65 deletions(-) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 5fa9d3a..817146f 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -108,7 +108,7 @@ set_tests_properties(confinement/test-surface-sampler-methods-build PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) # add the tests -set(_surface_tests test-points-union test-intersections-basic test-intersections-subtraction +set(_surface_tests test-points-subtraction test-points-basic test-intersections-basic test-intersections-subtraction test-intersections-union test-containment-subtraction test-containment-union) foreach(_test ${_surface_tests}) diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index f92ffcb..4acce34 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -28,7 +28,6 @@ #include "QBBC.hh" - // class for some basic vis class MyDetectorConstruction : public G4VUserDetectorConstruction { public: @@ -56,7 +55,7 @@ class MyDetectorConstruction : public G4VUserDetectorConstruction { // turn off vis auto worldVisAttr = new G4VisAttributes(); - worldVisAttr->SetVisibility(true); + worldVisAttr->SetVisibility(false); worldLog->SetVisAttributes(worldVisAttr); G4VPhysicalVolume* worldPhys = @@ -101,6 +100,73 @@ class MyDetectorConstruction : public G4VUserDetectorConstruction { }; +int RunVis(RMGVertexConfinement::SampleableObject obj,std::string name) +{ + G4double radius = + obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); + G4ThreeVector center = + obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); + + G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); + + std::vector points; + int i = 0; + while (i < 100) { + G4ThreeVector pos; + G4ThreeVector dir; + obj.GetDirection(dir, pos); + + points.push_back(pos); + if ((pos - center).mag() < radius) { + std::cout << "the initial point must be less than the bounding radius" << std::endl; + return 1; + } + i++; + } + G4RunManager* runManager = new G4RunManager(); + + // Set mandatory initialization classes + runManager->SetUserInitialization(new MyDetectorConstruction(center, radius, points, log)); + + auto physicsList = new QBBC; + physicsList->SetVerboseLevel(1); + runManager->SetUserInitialization(physicsList); + runManager->Initialize(); + // Initialize visualization + G4VisManager* visManager = new G4VisExecutive(); + visManager->Initialize(); + + // User interface + + // Set up visualization commands + G4UImanager* UImanager = G4UImanager::GetUIpointer(); + UImanager->ApplyCommand("/run/initialize "); + + UImanager->ApplyCommand("/vis/open OGL "); + UImanager->ApplyCommand("/vis/verbose warnings"); + + UImanager->ApplyCommand("/vis/drawVolume"); + UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); + UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); + UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); + + UImanager->ApplyCommand("/vis/scene/list "); + UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); + UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); + UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); + + UImanager->ApplyCommand("/vis/ogl/export test-points-"+name+".output.pdf"); + + + delete visManager; + delete runManager; + + + return 0; + + +} + G4VPhysicalVolume* get_phy_volume(G4VSolid* solid, std::string Material_string) { G4NistManager* nist = G4NistManager::Instance(); G4Material* Material = nist->FindOrBuildMaterial(Material_string); @@ -323,73 +389,25 @@ int main(int argc, char* argv[]) { // get some points to plot auto obj = sampleables["uni"]; - G4double radius = - obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); - G4ThreeVector center = - obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); - - G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); - - std::vector points; - int i = 0; - while (i < 1000) { - G4ThreeVector pos; - G4ThreeVector dir; - obj.GetDirection(dir, pos); - - points.push_back(pos); - if ((pos - center).mag() < radius) { - std::cout << "the initial point must be less than the bounding radius" << std::endl; - return 1; - } - i++; - } - G4RunManager* runManager = new G4RunManager(); - - // Set mandatory initialization classes - runManager->SetUserInitialization(new MyDetectorConstruction(center, radius, points, log)); - - auto physicsList = new QBBC; - physicsList->SetVerboseLevel(1); - runManager->SetUserInitialization(physicsList); - - // Initialize visualization - G4VisManager* visManager = new G4VisExecutive(); - visManager->Initialize(); - - // User interface - G4UIExecutive* ui = new G4UIExecutive(argc, argv); - - // Set up visualization commands - G4UImanager* UImanager = G4UImanager::GetUIpointer(); - UImanager->ApplyCommand("/run/initialize "); - - UImanager->ApplyCommand("/vis/open OGL "); - UImanager->ApplyCommand("/vis/verbose warnings"); - - UImanager->ApplyCommand("/vis/drawVolume"); - UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); - UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); - UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); - - UImanager->ApplyCommand("/vis/scene/list "); - UImanager->ApplyCommand("/vis/scene/add/axes 0 0 0 0.1 m"); - UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); - UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); - UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); - - // UImanager->ApplyCommand("/vis/viewer/set/background white"); - UImanager->ApplyCommand("/vis/ogl/export test-points-union.output.pdf"); + RunVis(obj,"union"); + } + else if (test_type == "test-points-subtraction") { - delete ui; - delete visManager; - delete runManager; + // get some points to plot + auto obj = sampleables["sub"]; + return RunVis(obj,"subtraction"); - return 0; + } + else if (test_type == "test-points-basic") { + // get some points to plot + auto obj = sampleables["tubs"]; + return RunVis(obj,"simple"); - } else { + } + + else { return 0; } From f459e8de3d35fd7946490b5481bcc069f48050ec Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 01:29:04 +0100 Subject: [PATCH 10/26] fix linking --- tests/confinement/CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 817146f..b991d16 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -100,6 +100,10 @@ find_package(BxDecay0 ${RMG_BXDECAY0_MINIMUM_VERSION} REQUIRED COMPONENTS Geant4 add_executable(test-surface-sampler-methods EXCLUDE_FROM_ALL test-surface-sampler-methods.cpp) target_link_libraries(test-surface-sampler-methods PUBLIC remage BxDecay0::BxDecay0) +if(BxDecay0_FOUND) + target_link_libraries(test-surface-sampler-methods PRIVATE BxDecay0::BxDecay0) +endif() + add_test(NAME confinement/test-surface-sampler-methods-build COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target test-surface-sampler-methods) From 577a3911849a3b08dc334ca6e1aea49962890ef2 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 14 Jan 2025 00:29:24 +0000 Subject: [PATCH 11/26] style: pre-commit fixes --- tests/confinement/CMakeLists.txt | 10 +- .../test-surface-sampler-methods.cpp | 119 +++++++++--------- 2 files changed, 65 insertions(+), 64 deletions(-) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index b991d16..2aa357f 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -112,8 +112,14 @@ set_tests_properties(confinement/test-surface-sampler-methods-build PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) # add the tests -set(_surface_tests test-points-subtraction test-points-basic test-intersections-basic test-intersections-subtraction - test-intersections-union test-containment-subtraction test-containment-union) +set(_surface_tests + test-points-subtraction + test-points-basic + test-intersections-basic + test-intersections-subtraction + test-intersections-union + test-containment-subtraction + test-containment-union) foreach(_test ${_surface_tests}) diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index 4acce34..ee28a84 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -100,71 +100,68 @@ class MyDetectorConstruction : public G4VUserDetectorConstruction { }; -int RunVis(RMGVertexConfinement::SampleableObject obj,std::string name) -{ - G4double radius = - obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); - G4ThreeVector center = - obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); - - G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); - - std::vector points; - int i = 0; - while (i < 100) { - G4ThreeVector pos; - G4ThreeVector dir; - obj.GetDirection(dir, pos); - - points.push_back(pos); - if ((pos - center).mag() < radius) { - std::cout << "the initial point must be less than the bounding radius" << std::endl; - return 1; - } - i++; +int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { + G4double radius = + obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); + G4ThreeVector center = + obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); + + G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); + + std::vector points; + int i = 0; + while (i < 100) { + G4ThreeVector pos; + G4ThreeVector dir; + obj.GetDirection(dir, pos); + + points.push_back(pos); + if ((pos - center).mag() < radius) { + std::cout << "the initial point must be less than the bounding radius" << std::endl; + return 1; } - G4RunManager* runManager = new G4RunManager(); - - // Set mandatory initialization classes - runManager->SetUserInitialization(new MyDetectorConstruction(center, radius, points, log)); - - auto physicsList = new QBBC; - physicsList->SetVerboseLevel(1); - runManager->SetUserInitialization(physicsList); - runManager->Initialize(); - // Initialize visualization - G4VisManager* visManager = new G4VisExecutive(); - visManager->Initialize(); + i++; + } + G4RunManager* runManager = new G4RunManager(); - // User interface + // Set mandatory initialization classes + runManager->SetUserInitialization(new MyDetectorConstruction(center, radius, points, log)); - // Set up visualization commands - G4UImanager* UImanager = G4UImanager::GetUIpointer(); - UImanager->ApplyCommand("/run/initialize "); + auto physicsList = new QBBC; + physicsList->SetVerboseLevel(1); + runManager->SetUserInitialization(physicsList); + runManager->Initialize(); + // Initialize visualization + G4VisManager* visManager = new G4VisExecutive(); + visManager->Initialize(); - UImanager->ApplyCommand("/vis/open OGL "); - UImanager->ApplyCommand("/vis/verbose warnings"); + // User interface - UImanager->ApplyCommand("/vis/drawVolume"); - UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); - UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); - UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); + // Set up visualization commands + G4UImanager* UImanager = G4UImanager::GetUIpointer(); + UImanager->ApplyCommand("/run/initialize "); - UImanager->ApplyCommand("/vis/scene/list "); - UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); - UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); - UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); + UImanager->ApplyCommand("/vis/open OGL "); + UImanager->ApplyCommand("/vis/verbose warnings"); - UImanager->ApplyCommand("/vis/ogl/export test-points-"+name+".output.pdf"); + UImanager->ApplyCommand("/vis/drawVolume"); + UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); + UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); + UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); + UImanager->ApplyCommand("/vis/scene/list "); + UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); + UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); + UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); - delete visManager; - delete runManager; + UImanager->ApplyCommand("/vis/ogl/export test-points-" + name + ".output.pdf"); - return 0; + delete visManager; + delete runManager; + return 0; } G4VPhysicalVolume* get_phy_volume(G4VSolid* solid, std::string Material_string) { @@ -389,25 +386,23 @@ int main(int argc, char* argv[]) { // get some points to plot auto obj = sampleables["uni"]; - RunVis(obj,"union"); + RunVis(obj, "union"); - } - else if (test_type == "test-points-subtraction") { + } else if (test_type == "test-points-subtraction") { // get some points to plot auto obj = sampleables["sub"]; - return RunVis(obj,"subtraction"); + return RunVis(obj, "subtraction"); - } - else if (test_type == "test-points-basic") { + } else if (test_type == "test-points-basic") { // get some points to plot auto obj = sampleables["tubs"]; - return RunVis(obj,"simple"); + return RunVis(obj, "simple"); + + } - } - - else { + else { return 0; } From 3e18805a10b88a0ea1db24fa2af47a01eb571a28 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 11:28:29 +0100 Subject: [PATCH 12/26] [tests] a bit of reformatting --- src/RMGVertexConfinement.cc | 7 +- tests/confinement/CMakeLists.txt | 50 ++++---- .../test-surface-sampler-methods.cpp | 119 +++++++++--------- 3 files changed, 86 insertions(+), 90 deletions(-) diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index 7c7e90c..5734cd5 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -126,8 +126,7 @@ bool RMGVertexConfinement::SampleableObject::IsInside(const G4ThreeVector& verte // get the intersections between the solid and an arbitrary line, defined by start and dir std::vector RMGVertexConfinement::SampleableObject::GetIntersections(G4ThreeVector start, G4ThreeVector dir) const { - - + // check if the physical volume exists if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject ", @@ -136,14 +135,12 @@ std::vector RMGVertexConfinement::SampleableObject::GetIntersecti auto solid = this->physical_volume->GetLogicalVolume()->GetSolid(); - G4double dist = 0; std::vector intersections = {}; G4ThreeVector int_point; int_point = start; int counter = 0; - dist = 0; while (dist < kInfinity) { @@ -169,7 +166,6 @@ std::vector RMGVertexConfinement::SampleableObject::GetIntersecti void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, G4ThreeVector& pos) const { - if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) RMGLog::OutDev(RMGLog::fatal, "Cannot generate directions for a SampleableObject ", "where the physical volume is not set this probably means you are trying to generically ", @@ -300,7 +296,6 @@ bool RMGVertexConfinement::SampleableObject::Sample(G4ThreeVector& vertex, int m } } - RMGLog::OutDev(RMGLog::debug, "Found good vertex ", vertex / CLHEP::cm, " cm", " after ", calls, " iterations, returning"); return true; diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index b991d16..f687ca9 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -96,33 +96,39 @@ set_tests_properties(confinment-lar/intersection-and-subtraction # test generic surface sampling methods find_package(BxDecay0 ${RMG_BXDECAY0_MINIMUM_VERSION} REQUIRED COMPONENTS Geant4) - add_executable(test-surface-sampler-methods EXCLUDE_FROM_ALL test-surface-sampler-methods.cpp) -target_link_libraries(test-surface-sampler-methods PUBLIC remage BxDecay0::BxDecay0) + +target_link_libraries(test-surface-sampler-methods PUBLIC remage) if(BxDecay0_FOUND) target_link_libraries(test-surface-sampler-methods PRIVATE BxDecay0::BxDecay0) -endif() - -add_test(NAME confinement/test-surface-sampler-methods-build - COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target - test-surface-sampler-methods) - -set_tests_properties(confinement/test-surface-sampler-methods-build - PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) - -# add the tests -set(_surface_tests test-points-subtraction test-points-basic test-intersections-basic test-intersections-subtraction - test-intersections-union test-containment-subtraction test-containment-union) -foreach(_test ${_surface_tests}) - - add_test(NAME confinment-surface/${_test} COMMAND ./test-surface-sampler-methods ${_test}) - - set_tests_properties(confinment-surface/${_test} PROPERTIES LABELS extra FIXTURES_REQUIRED - test-surface-sampler-methods-fixture) - -endforeach() + add_test(NAME confinement/test-surface-sampler-methods-build + COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target + test-surface-sampler-methods) + + set_tests_properties(confinement/test-surface-sampler-methods-build + PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) + + # add the tests + set(_surface_tests + test-points-subtraction + test-points-basic + test-intersections-basic + test-intersections-subtraction + test-intersections-union + test-containment-subtraction + test-containment-union) + + foreach(_test ${_surface_tests}) + + add_test(NAME confinment-surface/${_test} COMMAND ./test-surface-sampler-methods ${_test}) + set_tests_properties(confinment-surface/${_test} PROPERTIES LABELS extra FIXTURES_REQUIRED + test-surface-sampler-methods-fixture) + endforeach() +else () + message(WARNING "BxDecay0 not found - skipping tests depending on it.") +endif () # now run some simulations and test the outputs diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index 4acce34..f98421f 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -100,71 +100,68 @@ class MyDetectorConstruction : public G4VUserDetectorConstruction { }; -int RunVis(RMGVertexConfinement::SampleableObject obj,std::string name) -{ - G4double radius = - obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); - G4ThreeVector center = - obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); - - G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); - - std::vector points; - int i = 0; - while (i < 100) { - G4ThreeVector pos; - G4ThreeVector dir; - obj.GetDirection(dir, pos); - - points.push_back(pos); - if ((pos - center).mag() < radius) { - std::cout << "the initial point must be less than the bounding radius" << std::endl; - return 1; - } - i++; +int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { + G4double radius = + obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); + G4ThreeVector center = + obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); + + G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); + + std::vector points; + int i = 0; + while (i < 100) { + G4ThreeVector pos; + G4ThreeVector dir; + obj.GetDirection(dir, pos); + + points.push_back(pos); + if ((pos - center).mag() < radius) { + std::cout << "the initial point must be less than the bounding radius" << std::endl; + return 1; } - G4RunManager* runManager = new G4RunManager(); - - // Set mandatory initialization classes - runManager->SetUserInitialization(new MyDetectorConstruction(center, radius, points, log)); - - auto physicsList = new QBBC; - physicsList->SetVerboseLevel(1); - runManager->SetUserInitialization(physicsList); - runManager->Initialize(); - // Initialize visualization - G4VisManager* visManager = new G4VisExecutive(); - visManager->Initialize(); + i++; + } + G4RunManager* runManager = new G4RunManager(); - // User interface + // Set mandatory initialization classes + runManager->SetUserInitialization(new MyDetectorConstruction(center, radius, points, log)); - // Set up visualization commands - G4UImanager* UImanager = G4UImanager::GetUIpointer(); - UImanager->ApplyCommand("/run/initialize "); + auto physicsList = new QBBC; + physicsList->SetVerboseLevel(1); + runManager->SetUserInitialization(physicsList); + runManager->Initialize(); + // Initialize visualization + G4VisManager* visManager = new G4VisExecutive(); + visManager->Initialize(); - UImanager->ApplyCommand("/vis/open OGL "); - UImanager->ApplyCommand("/vis/verbose warnings"); + // User interface - UImanager->ApplyCommand("/vis/drawVolume"); - UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); - UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); - UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); + // Set up visualization commands + G4UImanager* UImanager = G4UImanager::GetUIpointer(); + UImanager->ApplyCommand("/run/initialize "); - UImanager->ApplyCommand("/vis/scene/list "); - UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); - UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); - UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); + UImanager->ApplyCommand("/vis/open OGL "); + UImanager->ApplyCommand("/vis/verbose warnings"); - UImanager->ApplyCommand("/vis/ogl/export test-points-"+name+".output.pdf"); + UImanager->ApplyCommand("/vis/drawVolume"); + UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); + UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); + UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); + UImanager->ApplyCommand("/vis/scene/list "); + UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); + UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); + UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); - delete visManager; - delete runManager; + UImanager->ApplyCommand("/vis/ogl/export surface-sample-bounding-box-" + name + ".output.pdf"); - return 0; + delete visManager; + delete runManager; + return 0; } G4VPhysicalVolume* get_phy_volume(G4VSolid* solid, std::string Material_string) { @@ -389,25 +386,23 @@ int main(int argc, char* argv[]) { // get some points to plot auto obj = sampleables["uni"]; - RunVis(obj,"union"); + RunVis(obj, "union"); - } - else if (test_type == "test-points-subtraction") { + } else if (test_type == "test-points-subtraction") { // get some points to plot auto obj = sampleables["sub"]; - return RunVis(obj,"subtraction"); + return RunVis(obj, "subtraction"); - } - else if (test_type == "test-points-basic") { + } else if (test_type == "test-points-basic") { // get some points to plot auto obj = sampleables["tubs"]; - return RunVis(obj,"simple"); + return RunVis(obj, "simple"); + + } - } - - else { + else { return 0; } From 257ac05c7e1e0eeb9c97fb56948ceb7ffdce696b Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 11:48:54 +0100 Subject: [PATCH 13/26] add some more comments --- src/RMGVertexConfinement.cc | 15 +++++++++++---- tests/confinement/CMakeLists.txt | 3 +-- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index 5734cd5..91f5829 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -143,6 +143,11 @@ std::vector RMGVertexConfinement::SampleableObject::GetIntersecti int counter = 0; dist = 0; + // Find the number of intersections between the shape and a line. + // Uses the DistanceToIn and DistanceToOut methods to find the next intersection, + // continually call these (depending on if we are inside or outside) until the distance + // becomes kInfinity. + while (dist < kInfinity) { dist = (counter % 2) == 0 ? solid->DistanceToIn(int_point, dir) @@ -171,7 +176,7 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, "where the physical volume is not set this probably means you are trying to generically ", "sample a geometrical volume which instead should be natively sampled"); - // get the bounding radius + // Get the center and radius of a bounding sphere around the shape G4double bounding_radius = physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); @@ -185,7 +190,6 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, // push in rho direction by some impact parameter G4double disk_phi = 2.0 * CLHEP::pi * G4UniformRand(); G4double disk_r = sqrt(G4UniformRand()) * bounding_radius; - pos += G4ThreeVector(cos(disk_phi) * disk_r, sin(disk_phi) * disk_r, 0); // now rotate pos and dir by some random direction @@ -197,11 +201,11 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, dir.rotateY(theta); dir.rotateZ(phi); + // shift by the barycenter of the bounding sphere. pos += barycenter; } -// generate with the generic sampler the actual surface point bool RMGVertexConfinement::SampleableObject::GenerateSurfacePoint(G4ThreeVector& vertex, int max_attempts, int n_max) const { @@ -219,7 +223,10 @@ bool RMGVertexConfinement::SampleableObject::GenerateSurfacePoint(G4ThreeVector& if (intersections.size() == 0) continue; - // pick one weighting by n_max and return it + // The surface sampling algorithm returns N intersections. + // We have to select one, to keep independence of the sampled points + // and weight by the number of intersections. + int random_int = static_cast(n_max * G4UniformRand()); if (random_int <= intersections.size() - 1) { diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index f687ca9..1fdee8d 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -114,6 +114,7 @@ if(BxDecay0_FOUND) set(_surface_tests test-points-subtraction test-points-basic + test-points-union test-intersections-basic test-intersections-subtraction test-intersections-union @@ -138,7 +139,6 @@ set_tests_properties(confinment-surface/gen-gdml PROPERTIES LABELS extra FIXTURE confine-surface-gdml-fixture) # visualise - set(_solids tubby box con uni sub trd) foreach(det ${_solids}) @@ -163,7 +163,6 @@ foreach(det ${_solids}) endforeach() # run simulations and check - foreach(det ${_solids}) add_test( From b61349496fc6a5b5b8f2c61275dc4c0e4eae8dca Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 12:32:50 +0100 Subject: [PATCH 14/26] [docs] adding docstrings for some functions --- src/RMGVertexConfinement.cc | 4 ++-- tests/confinement/CMakeLists.txt | 17 +++++++++-------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index 91f5829..3f39382 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -126,7 +126,7 @@ bool RMGVertexConfinement::SampleableObject::IsInside(const G4ThreeVector& verte // get the intersections between the solid and an arbitrary line, defined by start and dir std::vector RMGVertexConfinement::SampleableObject::GetIntersections(G4ThreeVector start, G4ThreeVector dir) const { - + // check if the physical volume exists if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject ", @@ -226,7 +226,7 @@ bool RMGVertexConfinement::SampleableObject::GenerateSurfacePoint(G4ThreeVector& // The surface sampling algorithm returns N intersections. // We have to select one, to keep independence of the sampled points // and weight by the number of intersections. - + int random_int = static_cast(n_max * G4UniformRand()); if (random_int <= intersections.size() - 1) { diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 1fdee8d..ed3fd45 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -104,11 +104,11 @@ if(BxDecay0_FOUND) target_link_libraries(test-surface-sampler-methods PRIVATE BxDecay0::BxDecay0) add_test(NAME confinement/test-surface-sampler-methods-build - COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target - test-surface-sampler-methods) + COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target + test-surface-sampler-methods) set_tests_properties(confinement/test-surface-sampler-methods-build - PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) + PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) # add the tests set(_surface_tests @@ -124,12 +124,13 @@ if(BxDecay0_FOUND) foreach(_test ${_surface_tests}) add_test(NAME confinment-surface/${_test} COMMAND ./test-surface-sampler-methods ${_test}) - set_tests_properties(confinment-surface/${_test} PROPERTIES LABELS extra FIXTURES_REQUIRED - test-surface-sampler-methods-fixture) + set_tests_properties( + confinment-surface/${_test} PROPERTIES LABELS extra FIXTURES_REQUIRED + test-surface-sampler-methods-fixture) endforeach() -else () - message(WARNING "BxDecay0 not found - skipping tests depending on it.") -endif () +else() + message(WARNING "BxDecay0 not found - skipping tests depending on it.") +endif() # now run some simulations and test the outputs From 3eed479a7dc2a70f5a5f2008556f193804a349a7 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 12:43:06 +0100 Subject: [PATCH 15/26] [docs] add some docstrings for generic surface sampling objects --- include/RMGVertexConfinement.hh | 58 +++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index b6dfd29..0c1f847 100644 --- a/include/RMGVertexConfinement.hh +++ b/include/RMGVertexConfinement.hh @@ -102,14 +102,72 @@ 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; + + /** + * @brief Generate a sample from the solid. + * + * @details Depending on if the solid is a basic one either sample natively, + * or using rejection sampling. Either samples the volume or the surface depending + * on the @c surface_sample member. + * + * @param vertex The sampled vertex. + * @param max_attempts The maximum number of candidate vertices for rejection sampling. + * @param force_containment_check Whether to force a check on where the point is + * inside the solid. + * @param n_trials The total number of trials performed. + * + */ [[nodiscard]] bool Sample(G4ThreeVector& vertex, int max_attempts, bool force_containment_check, long int& n_trials) const; + + + /** + * @brief Generate a point on the surface of the solid. + * + * @details This follows the algorithm from https://arxiv.org/abs/0802.2960. + * - Produce a direction vector corresponding to a uniform flux in a bounding sphere. + * - Find the intersections of this line with the solid. + * - Pick one intersection, or repeat. + * + * @param vertex The sampled vertex, + * @param max_samples The maximum number of attempts to find a valid vertex. + * @param n_max The maximum number of intersections possible for the solid, + * can be an overestimate. + * + */ [[nodiscard]] bool GenerateSurfacePoint(G4ThreeVector& vertex, int max_samples, int n_max) const; // methods for the generic surface sampling + /** + * @brief Get the number of intersections between the solid and the line starting at @c + * start with direction @c dir. + * + * @details This is used in the generic surface sampling algorithm. This function makes use + * of the methods @c GetDistanceToIn(p,v) and @c GetDistanceToOut(p,v) of @c G4VSolid . + * It continually looks for the distance to the next boundary (along the line) + * until this becomes zero indicating there are no more intersections. + * + * @param start The starting vector of the line, note this should be outside the solid. + * @param dir The direction vector. + * + * @returns A vector of the points of intersection. + */ std::vector GetIntersections(const G4ThreeVector start, const G4ThreeVector dir) const; + + + /** + * @brief Get a position and direction for the generic surface sampling algorithm. + * + * @details This generates a point on a bounding sphere, then shifts by some impact + * parameter, following the algorithm from https://arxiv.org/abs/0802.2960. This produces a + * uniform and isotropic flux inside the bounding sphere. + * + * @param dir The direction vector for the point. + * @param pos The initial position for the point. + * + */ void GetDirection(G4ThreeVector& dir, G4ThreeVector& pos) const; From 8cd69658f0299d26e509957a9c68306456fb629b Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 12:56:35 +0100 Subject: [PATCH 16/26] fix cmake --- tests/confinement/CMakeLists.txt | 47 +++++++++----------------------- 1 file changed, 13 insertions(+), 34 deletions(-) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 815b7dd..941a6dd 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -102,19 +102,22 @@ target_link_libraries(test-surface-sampler-methods PUBLIC remage) if(BxDecay0_FOUND) target_link_libraries(test-surface-sampler-methods PRIVATE BxDecay0::BxDecay0) -endif() -add_test(NAME confinement/test-surface-sampler-methods-build - COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target - test-surface-sampler-methods) + add_test(NAME confinement/test-surface-sampler-methods-build + COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target + test-surface-sampler-methods) -set_tests_properties(confinement/test-surface-sampler-methods-build - PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) + set_tests_properties(confinement/test-surface-sampler-methods-build + PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) +else() + message(WARNING "BxDecay0 not found - skipping tests depending on it.") +endif() # add the tests set(_surface_tests test-points-subtraction test-points-basic + test-points-union test-intersections-basic test-intersections-subtraction test-intersections-union @@ -123,34 +126,10 @@ set(_surface_tests foreach(_test ${_surface_tests}) - add_test(NAME confinement/test-surface-sampler-methods-build - COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target - test-surface-sampler-methods) - - set_tests_properties(confinement/test-surface-sampler-methods-build - PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) - - # add the tests - set(_surface_tests - test-points-subtraction - test-points-basic - test-points-union - test-intersections-basic - test-intersections-subtraction - test-intersections-union - test-containment-subtraction - test-containment-union) - - foreach(_test ${_surface_tests}) - - add_test(NAME confinment-surface/${_test} COMMAND ./test-surface-sampler-methods ${_test}) - set_tests_properties( - confinment-surface/${_test} PROPERTIES LABELS extra FIXTURES_REQUIRED - test-surface-sampler-methods-fixture) - endforeach() -else() - message(WARNING "BxDecay0 not found - skipping tests depending on it.") -endif() + add_test(NAME confinment-surface/${_test} COMMAND ./test-surface-sampler-methods ${_test}) + set_tests_properties(confinment-surface/${_test} PROPERTIES LABELS extra FIXTURES_REQUIRED + test-surface-sampler-methods-fixture) +endforeach() # now run some simulations and test the outputs From 913bba733337052e5206ed398acb9f76eefa52c5 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 12:59:55 +0100 Subject: [PATCH 17/26] more fixes to cmake --- tests/confinement/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 941a6dd..413d7f3 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -95,7 +95,8 @@ set_tests_properties(confinment-lar/intersection-and-subtraction PROPERTIES LABELS extra FIXTURES_REQUIRED confine-lar-ibn-out-output-fixture) # test generic surface sampling methods -find_package(BxDecay0 ${RMG_BXDECAY0_MINIMUM_VERSION} REQUIRED COMPONENTS Geant4) +find_package(BxDecay0 ${RMG_BXDECAY0_MINIMUM_VERSION} QUIET COMPONENTS Geant4) + add_executable(test-surface-sampler-methods EXCLUDE_FROM_ALL test-surface-sampler-methods.cpp) target_link_libraries(test-surface-sampler-methods PUBLIC remage) From ce8b2d8967cd651c22a092671bf9a031d2f3e967 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 13:36:40 +0100 Subject: [PATCH 18/26] Update tests/confinement/CMakeLists.txt Co-authored-by: Manuel Huber --- tests/confinement/CMakeLists.txt | 20 +++++++------------- 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 413d7f3..7305661 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -95,24 +95,18 @@ set_tests_properties(confinment-lar/intersection-and-subtraction PROPERTIES LABELS extra FIXTURES_REQUIRED confine-lar-ibn-out-output-fixture) # test generic surface sampling methods -find_package(BxDecay0 ${RMG_BXDECAY0_MINIMUM_VERSION} QUIET COMPONENTS Geant4) - add_executable(test-surface-sampler-methods EXCLUDE_FROM_ALL test-surface-sampler-methods.cpp) - target_link_libraries(test-surface-sampler-methods PUBLIC remage) - if(BxDecay0_FOUND) - target_link_libraries(test-surface-sampler-methods PRIVATE BxDecay0::BxDecay0) + target_link_libraries(test-surface-sampler-methods PRIVATE BxDecay0::BxDecay0_Geant4) +endif() - add_test(NAME confinement/test-surface-sampler-methods-build - COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target - test-surface-sampler-methods) +add_test(NAME confinement/test-surface-sampler-methods-build + COMMAND "${CMAKE_COMMAND}" --build "${CMAKE_BINARY_DIR}" --config "$" --target + test-surface-sampler-methods) - set_tests_properties(confinement/test-surface-sampler-methods-build - PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) -else() - message(WARNING "BxDecay0 not found - skipping tests depending on it.") -endif() +set_tests_properties(confinement/test-surface-sampler-methods-build + PROPERTIES LABELS extra FIXTURES_SETUP test-surface-sampler-methods-fixture) # add the tests set(_surface_tests From a795ac32f410dfc78d932d042a9b3e02916fa1f0 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 15:26:15 +0100 Subject: [PATCH 19/26] [docs] adding more documentation --- include/RMGVertexConfinement.hh | 63 +++++++++++++++++++++++++++++++-- 1 file changed, 60 insertions(+), 3 deletions(-) diff --git a/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index 0c1f847..1b4df80 100644 --- a/include/RMGVertexConfinement.hh +++ b/include/RMGVertexConfinement.hh @@ -33,16 +33,21 @@ class G4VPhysicalVolume; class G4VSolid; + +/** @brief Class for generating vertices in physical or geometrical volumes. + */ class RMGVertexConfinement : public RMGVVertexGenerator { public: + /** @brief Different types of geometrical (user) defined solids. */ enum class GeometricalSolidType { kSphere, kCylinder, kBox, }; + /** @brief Information about the geometrical (user) defined solids. */ struct GenericGeometricalSolidData { GeometricalSolidType solid_type; G4ThreeVector volume_center = G4ThreeVector(0, 0, 0); @@ -58,24 +63,37 @@ class RMGVertexConfinement : public RMGVVertexGenerator { double box_z_length = -1; }; + /** + * @brief Strategy for sampling physical and geometrical volumes. + * + * @details Can be either: + * - @c kIntersectPhysicalWithGeometrical : In which case vertices are generated in the + * intersection of the set of physical and geometrical volumes. + * - @c kUnionAll Generate in the union of all volumes, weighted by surface area / volume. + * - @c kSubtractGeometrical : Similar to @c kIntersectPhysicalWithGeometrical but specified + * regions can also be excluded. + */ enum class SamplingMode { kIntersectPhysicalWithGeometrical, kUnionAll, kSubtractGeometrical, }; + /** @brief Types of volume to sample, either physical (a volume in the geometry), geometrical + * (defined by the user) or unset. */ enum class VolumeType { kPhysical, kGeometrical, kUnset, }; - RMGVertexConfinement(); void BeginOfRunAction(const G4Run* run) override; void EndOfRunAction(const G4Run* run) override; + /** @brief Generate the actual vertex, according to the sampling mode (see \ref + * RMGVertexConfinement::SamplingMode). */ bool GenerateVertex(G4ThreeVector& v) override; // to be used in the messenger class @@ -93,14 +111,38 @@ class RMGVertexConfinement : public RMGVVertexGenerator { return fGeomVolumeData; } + /** + * An object which we can generate position samples in. Based on either a + * @c G4VPhysicalVolume or geometrical volume defined by a @c G4VSolid . The + * sampling can be performed either on the surface or in the volume of the solid. + */ struct SampleableObject { + SampleableObject() = default; SampleableObject(const SampleableObject&) = default; + + /** + * @brief SampleableObject constructor. + * + * @param v The physical volume. + * @param r A rotation matrix for the sampling solid. + * @param t A translation vector for the sampling solid. + * @param s A solid for geometrical volume sampling or for generating candidate points for + * rejection sampling. + * @param ns A flag of whether the solid is natively sampeable. + * @param ss A flag of whether the solid should be sampled on the surface. + */ SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r, G4ThreeVector t, G4VSolid* s, bool ns = false, bool ss = false); // NOTE: G4 volume/solid pointers should be fully owned by G4, avoid trying to delete them ~SampleableObject() = default; + + /** + * @brief Check if the vertex is inside the solid. + * @param vertex The sampled vertex. + * @returns Boolean flag of whether the vertexx is inside the solid. + */ [[nodiscard]] bool IsInside(const G4ThreeVector& vertex) const; /** @@ -109,6 +151,10 @@ class RMGVertexConfinement : public RMGVVertexGenerator { * @details Depending on if the solid is a basic one either sample natively, * or using rejection sampling. Either samples the volume or the surface depending * on the @c surface_sample member. + * - For surface sampling mode the solid is either natively sampled (if this is + * implemented), or is sampled with \ref SampleableObject::GenerateSurfacePoint + * - For volume sampling, if the solid is not natively sampleable, points are generated in a + * bounding box and then rejection sampling is used using \ref SampleableObject::IsInside. * * @param vertex The sampled vertex. * @param max_attempts The maximum number of candidate vertices for rejection sampling. @@ -120,7 +166,6 @@ class RMGVertexConfinement : public RMGVVertexGenerator { [[nodiscard]] bool Sample(G4ThreeVector& vertex, int max_attempts, bool force_containment_check, long int& n_trials) const; - /** * @brief Generate a point on the surface of the solid. * @@ -170,24 +215,36 @@ class RMGVertexConfinement : public RMGVVertexGenerator { */ void GetDirection(G4ThreeVector& dir, G4ThreeVector& pos) const; - G4VPhysicalVolume* physical_volume = nullptr; G4VSolid* sampling_solid = nullptr; G4RotationMatrix rotation; G4ThreeVector translation; + double volume = -1; double surface = -1; + bool surface_sample = false; bool native_sample = false; int max_num_intersections = -1; }; + /** A collection of @c SampleableObjects . Can be used + * to sample from by selecting a volume weighted by surface area + * or volume. + */ struct SampleableObjectCollection { SampleableObjectCollection() = default; inline ~SampleableObjectCollection() { data.clear(); } + /** @brief Select a @c SampleableObject from the collection, weighted by surface area. + * @returns a reference to the chosen @c SampleableObject . + */ [[nodiscard]] const SampleableObject& SurfaceWeightedRand() const; + + /** @brief Select a @c SampleableObject from the collection, weighted by volume. + * @returns a reference to the chosen @c SampleableObject . + */ [[nodiscard]] const SampleableObject& VolumeWeightedRand() const; [[nodiscard]] bool IsInside(const G4ThreeVector& vertex) const; From 5678c73780406b0432543523736e5cd269e26210 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Tue, 14 Jan 2025 18:22:56 +0100 Subject: [PATCH 20/26] [docs] small changes --- docs/developer.md | 8 ++++++-- tests/confinement/CMakeLists.txt | 7 ++++--- tests/confinement/test-surface-sampler-methods.cpp | 3 ++- tests/confinement/test_basic_surface.py | 2 +- tests/confinement/test_confinment_ge.py | 2 +- 5 files changed, 14 insertions(+), 8 deletions(-) diff --git a/docs/developer.md b/docs/developer.md index 0be3fd0..a82aca0 100644 --- a/docs/developer.md +++ b/docs/developer.md @@ -62,9 +62,13 @@ $ cmake -DCMAKE_INSTALL_PREFIX= .. $ make install ``` -```{tip} +:::{warning} +If you want to run the _remage_ tests the cmake flag `-DBUILD_TESTING=ON` is required. +::: + +:::{note} A list of available Make targets can be printed by running `make help`. -``` +::: ## Code style diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 7305661..b83832e 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -70,8 +70,9 @@ set_tests_properties( 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) +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) @@ -92,7 +93,7 @@ add_test(NAME confinment-lar/intersection-and-subtraction 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) + PROPERTIES LABELS extra FIXTURES_REQUIRED confine-lar-in-out-output-fixture) # test generic surface sampling methods add_executable(test-surface-sampler-methods EXCLUDE_FROM_ALL test-surface-sampler-methods.cpp) diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index f98421f..ef757ff 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -154,7 +154,8 @@ int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); - UImanager->ApplyCommand("/vis/ogl/export surface-sample-bounding-box-" + name + ".output.pdf"); + UImanager->ApplyCommand( + "/vis/ogl/export surface-sample-bounding-box-" + name + ".output.pdf"); delete visManager; diff --git a/tests/confinement/test_basic_surface.py b/tests/confinement/test_basic_surface.py index f80f1c8..259aab6 100644 --- a/tests/confinement/test_basic_surface.py +++ b/tests/confinement/test_basic_surface.py @@ -182,7 +182,7 @@ def add_local_pos(vertices, pos): dtype = args.det -with PdfPages(f"simple-solids-surface-{dtype}.output.pdf") as pdf: +with PdfPages(f"confinement.simple-solids-surface-{dtype}.output.pdf") as pdf: # get positions pos = reg.physicalVolumeDict[dtype].position.eval() diff --git a/tests/confinement/test_confinment_ge.py b/tests/confinement/test_confinment_ge.py index 82feedd..e01f4f2 100644 --- a/tests/confinement/test_confinment_ge.py +++ b/tests/confinement/test_confinment_ge.py @@ -167,6 +167,6 @@ def make_plot(vert, hit): p, sigma = make_plot(vertices, hits) -plt.savefig("confinement-ge.output.pdf") +plt.savefig("relative-ge.output.pdf") assert sigma < 5 From 2751f066d881fa08bc4e88e648be854ec237e19c Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Wed, 15 Jan 2025 01:31:23 +0100 Subject: [PATCH 21/26] [tests] first attempt to organise all the test outputs into a LaTeX document --- Testing/Temporary/LastTest.log | 3 + tests/CMakeLists.txt | 1 + tests/basics/CMakeLists.txt | 2 +- tests/basics/template-basic.tex | 28 +++ tests/confinement/CMakeLists.txt | 7 +- tests/confinement/macros/_vis.mac | 2 +- tests/confinement/template-confinement.tex | 171 ++++++++++++++++++ .../test-surface-sampler-methods.cpp | 3 +- tests/germanium/CMakeLists.txt | 2 +- tests/tex/CMakeLists.txt | 27 +++ tests/tex/validation-report.tex | 27 +++ 11 files changed, 264 insertions(+), 9 deletions(-) create mode 100644 Testing/Temporary/LastTest.log create mode 100644 tests/basics/template-basic.tex create mode 100644 tests/confinement/template-confinement.tex create mode 100644 tests/tex/CMakeLists.txt create mode 100644 tests/tex/validation-report.tex diff --git a/Testing/Temporary/LastTest.log b/Testing/Temporary/LastTest.log new file mode 100644 index 0000000..851b12a --- /dev/null +++ b/Testing/Temporary/LastTest.log @@ -0,0 +1,3 @@ +Start testing: Jan 14 00:54 CET +---------------------------------------------------------- +End testing: Jan 14 00:54 CET diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4b1ce4f..90ba5ec 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -14,3 +14,4 @@ add_subdirectory(internals) add_subdirectory(output) add_subdirectory(python) add_subdirectory(vertex) +add_subdirectory(tex) diff --git a/tests/basics/CMakeLists.txt b/tests/basics/CMakeLists.txt index f30b7f8..b701cff 100644 --- a/tests/basics/CMakeLists.txt +++ b/tests/basics/CMakeLists.txt @@ -2,7 +2,7 @@ file( GLOB _aux RELATIVE ${PROJECT_SOURCE_DIR} - macros/*.mac gdml/*.gdml gdml/*.xml) + macros/*.mac gdml/*.gdml gdml/*.xml *.tex) # copy them to the build area foreach(_file ${_aux}) diff --git a/tests/basics/template-basic.tex b/tests/basics/template-basic.tex new file mode 100644 index 0000000..89e2cb2 --- /dev/null +++ b/tests/basics/template-basic.tex @@ -0,0 +1,28 @@ +The first set of checks shown below test that remage is running and has reasonable seeming results. +\IfFileExists{vis-2nbb.output_0000.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{vis-2nbb.output_0000.pdf} + \caption{Geant4 visualisation of a simulation of $2\nu\beta\beta$ decay in a HPGe detector for a typical screeing \ + setup. You should see the trajectories focused around the HPGe detector and mainly electrons (red), with a \ + few gamma (green) being produced. \ + } + \end{figure} +}{ + \noindent \fbox{\textbf{The test of visualising $2\nu\beta\beta$ either failed or was not run. + }} +} + + +\IfFileExists{vis-co60.output_0000.pdf}{ + \begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{vis-co60.output_0000.pdf} + \caption{Geant4 visualisation of a simulation of $^{60}Co$ decay in a source close to a HPGe + detector for a typical screening station setup. You should see $\gamma$ particles (in green) throughout + the experimental setup starting from a source above the detector.} + \end{figure} +}{ + \noindent \fbox{\textbf{The test of visualising $^{60}Co$ either failed or was not run.}} +} diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index b83832e..cf4f530 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -2,7 +2,7 @@ file( GLOB _aux RELATIVE ${PROJECT_SOURCE_DIR} - macros/*.mac macros/*.json gdml/*.gdml gdml/*.xml *.py) + macros/*.mac macros/*.json gdml/*.gdml gdml/*.xml *.py *.tex) # copy them to the build area foreach(_file ${_aux}) @@ -70,9 +70,8 @@ set_tests_properties( 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) +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) diff --git a/tests/confinement/macros/_vis.mac b/tests/confinement/macros/_vis.mac index 894be17..0e6910d 100644 --- a/tests/confinement/macros/_vis.mac +++ b/tests/confinement/macros/_vis.mac @@ -9,7 +9,7 @@ /vis/viewer/set/lineSegmentsPerCircle 10 # /vis/viewer/set/defaultColour black # /vis/viewer/set/background white -# /vis/scene/add/axes 0 0 0 3 m +/vis/scene/add/axes 0 0 0 3 m /vis/viewer/set/upVector 0 0 -1 /vis/viewer/set/viewpointVector -1 1 1 diff --git a/tests/confinement/template-confinement.tex b/tests/confinement/template-confinement.tex new file mode 100644 index 0000000..5472a92 --- /dev/null +++ b/tests/confinement/template-confinement.tex @@ -0,0 +1,171 @@ + +These tests check the generation of primary vertices by remage. +\subsection{Native sampling} +The first checks plot the positions for natively sampleable objects. + + +\IfFileExists{native-surface.output_0000.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{native-surface.output_0000.pdf} + \caption{Geant4 visualisation of a simulation of the primary positions for the surface of + some natively sampled shapes. Primaries should be present on the G4Box (bottom row second from right), + G4Orb (full sphere - above the box) and G4Tubs (left). The primaries should be on the detector surface, although this + is hard to appreciate with Geant4 visualisation. + } + \end{figure} +}{ + \noindent \fbox{\textbf{The test of generating natively sampled surface primaries either failed or did not run. + }} +} + + +\IfFileExists{native-surface.output_0000.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{native-volume.output_0000.pdf} + \caption{Geant4 visualisation of a simulation of the primary positions for the volume of + some natively sampled shapes. Primaries should be present on the G4Box (bottom row second from right), + G4Orb (full sphere - above the box), G4Tubs and G4Sphere (sector of a sphere - bottom center). The primaries should be on the detector surface, although this + be in the detector volume although this is hard to appreciate with Geant4 visualisation. + } + \end{figure} +}{ + \noindent \fbox{\textbf{The test of generating natively sampled volume primaries either failed or did not run. + }} +} + +Finally, we check the generation of points on complex objects. + + +\IfFileExists{complex-volume.output_0000.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{complex-volume.output_0000.pdf} + \caption{Geant4 visualisation of a simulation of the primary positions for the volume of + some complex shapes. A polycone (top center), a box with a hole (center), the union of an orb and box (upper right) + and an orb (right). + } + \end{figure} +}{ + \noindent \fbox{\textbf{The test of generating complex sampled volume primaries either failed or did not run. + }} +} + + + +\subsection{Geometrical volumes and intersections} +The next tests check the generation of points in geometrical (user defined) volumes and intersections / unions of +geometrical and physical volumes. + +\IfFileExists{geometrical-output_0000.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{geometrical-output_0000.pdf} + \caption{Geant4 visualisation of a simulation of the primary positions for some geometrical (user defined ) volumes. + The primaries should be located in a sphere of center $(0,2,2)$ m and radius $0.7$ m (center overlapping box), a partially filled cylinder with center + $(-2,-2,-2)$ m and outer radius $0.7$ m and height 1 m (far right) and a box of center $(0,2,-2)$ m and sides of length 0.7, 0.5 and 1.2 m (top). + } + \end{figure} +}{ + \noindent \fbox{\textbf{The test of generating geometrically volume primaries either failed or did not run. + }} +} + +\IfFileExists{geometrical-and-physical-output_0000.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{geometrical-and-physical-output_0000.pdf} + \caption{Geant4 visualisation of a simulation of the primary positions for the Union of physical volumes defined + by the union of the box and orb (upper left) and the partially filled sphere (lower center) and a geometric volume defined by + sphere of center $(0,2,2)$ m and radius $0.7$ m near the center overlapping the box. + The primaries should be in these three locations.} + \end{figure} +}{ + \noindent \fbox{\textbf{The test of generating geometric or physical volume primaries either failed or did not run. + }} +} +\IfFileExists{geometrical-and-physical-output_0000.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{geometrical-and-physical-output_0000.pdf} + \caption{Geant4 visualisation of a simulation of the primary positions for the intersection of + physical volumes defined by a user defined sphere of center $(1.5,2,2)$ m and radius $0.4$ m and the orb in the lower right. + The primaries should be located in a subset of this orb.} + \end{figure} +}{ + \noindent \fbox{\textbf{The test of generating geometric and physical volume primaries either failed or did not run. + }} +} + +\subsection{Union and intersection tests} +The next test simulates $2\nu\beta\beta$ decay in a hypothetical HPGe array and makes a +statistical comparison of the number of primaries in each volume. +\IfFileExists{relative-ge.output.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.9\textwidth]{relative-ge.output.pdf} + + \end{figure} +}{ + \noindent \fbox{\textbf{Relative HPGe sampling test either failed or did not run + }} +} + +To check the intersection we generate events inside a cylinder around the HPGe detectors and again check the +positions and relative fractions. + +\IfFileExists{lar-in-check.output.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.9\textwidth,page=1]{lar-in-check.output.pdf} + \includegraphics[width=0.9\textwidth,page=2]{lar-in-check.output.pdf} + \includegraphics[width=0.9\textwidth,page=3]{lar-in-check.output.pdf} + + \end{figure} +}{ + \noindent \fbox{\textbf{Intersection of physical and geometrical volume test (LAr cylinder) either failed or did not run + }} +} + +Finally check subtraction by generating events excluding these regions. + +\IfFileExists{lar-sub-check.output.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.9\textwidth,page=1]{lar-sub-check.output.pdf} + \includegraphics[width=0.9\textwidth,page=2]{lar-sub-check.output.pdf} + + \end{figure} +}{ + \noindent \fbox{\textbf{Subtraction of physical and geometrical volume test (LAr cylinder) either failed or did not run + }} + +} +Finally, both an intersection and a subtraction. + +\IfFileExists{lar-int-and-sub-check.output.pdf}{ + + \begin{figure}[h!] + \centering + \includegraphics[width=0.9\textwidth,page=1]{lar-int-and-sub-check.output.pdf} + \includegraphics[width=0.9\textwidth,page=2]{lar-int-and-sub-check.output.pdf} + + \end{figure} +}{ + \noindent \fbox{\textbf{Subtraction of physical and geometrical volume test (LAr cylinder) either failed or did not run + }} + +} + +\subsection{Generic surface sampling} +The final part of this report describes checks on the generic surface sampler algorithm. diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index ef757ff..f98421f 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -154,8 +154,7 @@ int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); - UImanager->ApplyCommand( - "/vis/ogl/export surface-sample-bounding-box-" + name + ".output.pdf"); + UImanager->ApplyCommand("/vis/ogl/export surface-sample-bounding-box-" + name + ".output.pdf"); delete visManager; diff --git a/tests/germanium/CMakeLists.txt b/tests/germanium/CMakeLists.txt index 9fde7ca..6018bf2 100644 --- a/tests/germanium/CMakeLists.txt +++ b/tests/germanium/CMakeLists.txt @@ -1,7 +1,7 @@ file( GLOB _aux RELATIVE ${PROJECT_SOURCE_DIR} - macros/*.mac gdml/*.gdml *.py) + macros/*.mac gdml/*.gdml *.py *.tex) # copy them to the build area foreach(_file ${_aux}) diff --git a/tests/tex/CMakeLists.txt b/tests/tex/CMakeLists.txt new file mode 100644 index 0000000..31afa1d --- /dev/null +++ b/tests/tex/CMakeLists.txt @@ -0,0 +1,27 @@ +file( + GLOB _aux + RELATIVE ${PROJECT_SOURCE_DIR} + *.tex) + +# copy them to the build area +foreach(_file ${_aux}) + configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) +endforeach() + +# copy all files +add_test( + NAME tex/copy-files + COMMAND + /bin/sh -c + "find ../ -type f \\( -name '*.output*' -o -name '*.tex' \\) ! -path '../tex/*' -exec cp {} . \; " +) +set_tests_properties(tex/copy-files PROPERTIES LABELS tex FIXTURES_SETUP copy-files-fixture) + +find_program(PDFLATEX_COMPILER pdflatex) +if(PDFLATEX_COMPILER) + add_test(NAME tex/compile COMMAND ${PDFLATEX_COMPILER} validation-report.tex) + set_tests_properties(tex/copy-files PROPERTIES LABELS tex FIXTURES_REQUIRED copy-files-fxiture) + +else() + message(STATUS "No LaTeX compiler found, report cannot be compiled.") +endif() diff --git a/tests/tex/validation-report.tex b/tests/tex/validation-report.tex new file mode 100644 index 0000000..ecb2787 --- /dev/null +++ b/tests/tex/validation-report.tex @@ -0,0 +1,27 @@ +\documentclass[a4paper,12pt]{article} + +% Packages for enhanced functionality +\usepackage{amsmath} % For math equations +\usepackage{amsfonts} % For special fonts +\usepackage{graphicx} % For including graphics +\usepackage{hyperref} % For hyperlinks +\usepackage{geometry} % To adjust margins + +% Adjust margins +\geometry{left=1in, right=1in, top=1in, bottom=0.5in} + +\title{{\Large \bf remage validation suite.}} +\author{{\normalsize Citation: \url{https://doi.org/10.5281/zenodo.11115662}}} +\date{{\it\normalsize \today}} % Date is automatically set to the current date + +\begin{document} +\maketitle + +\section{Basics} +\input{template-basic} + +\section{Confinement} +\input{template-confinement} + + +\end{document} From 7a2a1d7dc8538e7483ea03ef45c4dfeff7ca3945 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Wed, 15 Jan 2025 01:36:09 +0100 Subject: [PATCH 22/26] update github actions to save the .tex --- .github/workflows/main.yml | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index abfc504..b5d78a6 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -58,8 +58,10 @@ jobs: uses: actions/upload-artifact@v4 with: name: test-output-${{ matrix.container_version }} - # artifacts must have a ".output*" extension - path: build/tests/**/*.output*.* + # artifacts must have a ".output*" extension or .tex + path: | + build/tests/tex/*.output*.* + build/tests/tex/*.tex - name: Run minimal test suite if: ${{ matrix.container_version == 'slim' }} From 458c62160d62e7b8c92d0f38c0a84bb5bb71f680 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Thu, 16 Jan 2025 11:37:23 +0100 Subject: [PATCH 23/26] [tests] update plot formats --- tests/basics/template-basic.tex | 34 +-- tests/confinement/CMakeLists.txt | 2 +- tests/confinement/macros/_vis-export.mac | 3 +- tests/confinement/macros/_vis.mac | 2 +- tests/confinement/macros/complex-surface.mac | 19 ++ tests/confinement/macros/complex-volume.mac | 2 +- .../macros/geometrical-and-physical.mac | 2 +- .../macros/geometrical-or-physical.mac | 2 +- tests/confinement/macros/geometrical.mac | 2 +- tests/confinement/macros/native-volume.mac | 2 +- tests/confinement/macros/vis-surface.mac | 2 +- tests/confinement/template-confinement.tex | 277 +++++++++--------- .../test-surface-sampler-methods.cpp | 2 +- tests/tex/validation-report.tex | 3 +- 14 files changed, 187 insertions(+), 167 deletions(-) create mode 100644 tests/confinement/macros/complex-surface.mac diff --git a/tests/basics/template-basic.tex b/tests/basics/template-basic.tex index 89e2cb2..00fb37c 100644 --- a/tests/basics/template-basic.tex +++ b/tests/basics/template-basic.tex @@ -1,28 +1,12 @@ The first set of checks shown below test that remage is running and has reasonable seeming results. -\IfFileExists{vis-2nbb.output_0000.pdf}{ - \begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{vis-2nbb.output_0000.pdf} - \caption{Geant4 visualisation of a simulation of $2\nu\beta\beta$ decay in a HPGe detector for a typical screeing \ - setup. You should see the trajectories focused around the HPGe detector and mainly electrons (red), with a \ - few gamma (green) being produced. \ - } - \end{figure} -}{ - \noindent \fbox{\textbf{The test of visualising $2\nu\beta\beta$ either failed or was not run. - }} -} +\begin{figure}[h!] + \centering + \includegraphics[width=0.45\textwidth]{vis-2nbb.output_0000.pdf} + \includegraphics[width=0.45\textwidth]{vis-co60.output_0000.pdf} - -\IfFileExists{vis-co60.output_0000.pdf}{ - \begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{vis-co60.output_0000.pdf} - \caption{Geant4 visualisation of a simulation of $^{60}Co$ decay in a source close to a HPGe - detector for a typical screening station setup. You should see $\gamma$ particles (in green) throughout - the experimental setup starting from a source above the detector.} - \end{figure} -}{ - \noindent \fbox{\textbf{The test of visualising $^{60}Co$ either failed or was not run.}} -} + \caption{Geant4 visualisation of a simulation of $2\nu\beta\beta$ decay in a HPGe detector (left) and $^{60}$Co in a source above the detector (right( for a typical screeing \ + setup. For the left figure you should see the trajectories focused around the HPGe detector and mainly electrons (red), with a \ + few gamma (green) being produced. For the right you should see $\gamma$ particles (in green) throughout the experimental setup starting from a source above the detector + } +\end{figure} diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index cf4f530..663f546 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -9,7 +9,7 @@ foreach(_file ${_aux}) configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) endforeach() -set(_macros complex-volume.mac geometrical.mac native-surface.mac geometrical-and-physical.mac +set(_macros complex-volume.mac complex-surface.mac geometrical.mac native-surface.mac geometrical-and-physical.mac geometrical-or-physical.mac native-volume.mac) foreach(_mac ${_macros}) diff --git a/tests/confinement/macros/_vis-export.mac b/tests/confinement/macros/_vis-export.mac index 7b58471..5f5962b 100644 --- a/tests/confinement/macros/_vis-export.mac +++ b/tests/confinement/macros/_vis-export.mac @@ -1 +1,2 @@ -/vis/ogl/export {export-fn} +/vis/ogl/set/printSize 300 300 +/vis/ogl/export {export-fn}.pdf diff --git a/tests/confinement/macros/_vis.mac b/tests/confinement/macros/_vis.mac index 0e6910d..a0e3be6 100644 --- a/tests/confinement/macros/_vis.mac +++ b/tests/confinement/macros/_vis.mac @@ -2,7 +2,7 @@ /RMG/Manager/Logging/LogLevel summary -/vis/open OGL +/vis/open OGL 600x600-0+0 /vis/drawVolume /vis/sceneHandler/attach diff --git a/tests/confinement/macros/complex-surface.mac b/tests/confinement/macros/complex-surface.mac new file mode 100644 index 0000000..aa4d68c --- /dev/null +++ b/tests/confinement/macros/complex-surface.mac @@ -0,0 +1,19 @@ +/control/execute macros/_init.mac + +/RMG/Generator/Confine Volume +/RMG/Generator/Confinement/ForceContainmentCheck true +/RMG/Generator/Confinement/SampleOnSurface true +/RMG/Generator/Confinement/MaxNumberOfIntersections 6 + +/RMG/Generator/Confinement/Physical/AddVolume BoxWithHole +/RMG/Generator/Confinement/Physical/AddVolume BoxAndOrb +#/RMG/Generator/Confinement/Physical/AddVolume MotherOrb +/RMG/Generator/Confinement/Physical/AddVolume Polycone + +/RMG/Generator/Select GPS +/gps/particle e- +/gps/energy 1 eV + +/run/beamOn 5000 + +/control/alias export-fn complex-surface.output diff --git a/tests/confinement/macros/complex-volume.mac b/tests/confinement/macros/complex-volume.mac index cd63da1..964d192 100644 --- a/tests/confinement/macros/complex-volume.mac +++ b/tests/confinement/macros/complex-volume.mac @@ -15,4 +15,4 @@ /run/beamOn 5000 -/control/alias export-fn complex-volume.output.pdf +/control/alias export-fn complex-volume.output diff --git a/tests/confinement/macros/geometrical-and-physical.mac b/tests/confinement/macros/geometrical-and-physical.mac index c19dd3c..06abf86 100644 --- a/tests/confinement/macros/geometrical-and-physical.mac +++ b/tests/confinement/macros/geometrical-and-physical.mac @@ -20,4 +20,4 @@ /run/beamOn 2000 -/control/alias export-fn geometrical-and-physical.output.pdf +/control/alias export-fn geometrical-and-physical.output.jpg diff --git a/tests/confinement/macros/geometrical-or-physical.mac b/tests/confinement/macros/geometrical-or-physical.mac index 8b6340f..60932ef 100644 --- a/tests/confinement/macros/geometrical-or-physical.mac +++ b/tests/confinement/macros/geometrical-or-physical.mac @@ -21,4 +21,4 @@ /run/beamOn 2000 -/control/alias export-fn geometrical-or-physical.output.pdf +/control/alias export-fn geometrical-or-physical.output.jpg diff --git a/tests/confinement/macros/geometrical.mac b/tests/confinement/macros/geometrical.mac index e14805b..f504d5b 100644 --- a/tests/confinement/macros/geometrical.mac +++ b/tests/confinement/macros/geometrical.mac @@ -36,4 +36,4 @@ /run/beamOn 5000 -/control/alias export-fn geometrical.output.pdf +/control/alias export-fn geometrical.output.jpg diff --git a/tests/confinement/macros/native-volume.mac b/tests/confinement/macros/native-volume.mac index b5804af..5c0dcac 100644 --- a/tests/confinement/macros/native-volume.mac +++ b/tests/confinement/macros/native-volume.mac @@ -13,4 +13,4 @@ /run/beamOn 2000 -/control/alias export-fn native-volume.output.pdf +/control/alias export-fn native-volume.output diff --git a/tests/confinement/macros/vis-surface.mac b/tests/confinement/macros/vis-surface.mac index 84dc15b..0a88f41 100644 --- a/tests/confinement/macros/vis-surface.mac +++ b/tests/confinement/macros/vis-surface.mac @@ -56,4 +56,4 @@ /gps/energy 1 eV /run/beamOn 1000 -/control/alias export-fn vis-surface-${det}.output.pdf +/control/alias export-fn vis-surface-${det}.output.png diff --git a/tests/confinement/template-confinement.tex b/tests/confinement/template-confinement.tex index 5472a92..c57873e 100644 --- a/tests/confinement/template-confinement.tex +++ b/tests/confinement/template-confinement.tex @@ -4,56 +4,30 @@ \subsection{Native sampling} The first checks plot the positions for natively sampleable objects. -\IfFileExists{native-surface.output_0000.pdf}{ - - \begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{native-surface.output_0000.pdf} - \caption{Geant4 visualisation of a simulation of the primary positions for the surface of - some natively sampled shapes. Primaries should be present on the G4Box (bottom row second from right), - G4Orb (full sphere - above the box) and G4Tubs (left). The primaries should be on the detector surface, although this - is hard to appreciate with Geant4 visualisation. - } - \end{figure} -}{ - \noindent \fbox{\textbf{The test of generating natively sampled surface primaries either failed or did not run. - }} -} - - -\IfFileExists{native-surface.output_0000.pdf}{ - - \begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{native-volume.output_0000.pdf} - \caption{Geant4 visualisation of a simulation of the primary positions for the volume of - some natively sampled shapes. Primaries should be present on the G4Box (bottom row second from right), - G4Orb (full sphere - above the box), G4Tubs and G4Sphere (sector of a sphere - bottom center). The primaries should be on the detector surface, although this - be in the detector volume although this is hard to appreciate with Geant4 visualisation. - } - \end{figure} -}{ - \noindent \fbox{\textbf{The test of generating natively sampled volume primaries either failed or did not run. - }} -} -Finally, we check the generation of points on complex objects. +\begin{figure}[h!] + \centering + \includegraphics[width=0.45\textwidth]{native-surface.output_0000.png} + \includegraphics[width=0.45\textwidth]{native-volume.output_0000.png} + \caption{Geant4 visualisation of a simulation of the primary positions for the surface (left) and volume (right) of + some natively sampled shapes. Primaries should be present on the G4Box (bottom row second from right), + G4Orb (full sphere - above the box) and G4Tubs (left). The primaries should be on the detector surface or volume, although this + is hard to appreciate with Geant4 visualisation. + } +\end{figure} -\IfFileExists{complex-volume.output_0000.pdf}{ +Finally, we check the generation of points on complex objects. +\begin{figure}[h!] + \centering + \includegraphics[width=0.45\textwidth]{complex-surface.output_0000.png} + \includegraphics[width=0.45\textwidth]{complex-volume.output_0000.png} - \begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{complex-volume.output_0000.pdf} - \caption{Geant4 visualisation of a simulation of the primary positions for the volume of - some complex shapes. A polycone (top center), a box with a hole (center), the union of an orb and box (upper right) - and an orb (right). - } - \end{figure} -}{ - \noindent \fbox{\textbf{The test of generating complex sampled volume primaries either failed or did not run. - }} -} + \caption{Geant4 visualisation of a simulation of the primary positions for the surface (left) or volume (right) of + some complex shapes. A polycone (top center), a box with a hole (center), the union of an orb and box (upper right) + and an orb (right). + } +\end{figure} @@ -61,111 +35,152 @@ \subsection{Geometrical volumes and intersections} The next tests check the generation of points in geometrical (user defined) volumes and intersections / unions of geometrical and physical volumes. -\IfFileExists{geometrical-output_0000.pdf}{ +\begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{geometrical.output_0000.png} + \caption{Geant4 visualisation of a simulation of the primary positions for some geometrical (user defined ) volumes. + The primaries should be located in a sphere of center $(0,2,2)$ m and radius $0.7$ m (center overlapping box), a partially filled cylinder with center + $(-2,-2,-2)$ m and outer radius $0.7$ m and height 1 m (far right) and a box of center $(0,2,-2)$ m and sides of length 0.7, 0.5 and 1.2 m (top). + } +\end{figure} + +\begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{geometrical-or-physical.output_0000.png} + \caption{Geant4 visualisation of a simulation of the primary positions for the Union of physical volumes defined + by the union of the box and orb (upper left) and the partially filled sphere (lower center) and a geometric volume defined by + sphere of center $(0,2,2)$ m and radius $0.7$ m near the center overlapping the box. + The primaries should be in these three locations.} +\end{figure} + + +\begin{figure}[h!] + \centering + \includegraphics[width=0.6\textwidth]{geometrical-and-physical.output_0000.png} + \caption{Geant4 visualisation of a simulation of the primary positions for the intersection of + physical volumes defined by a user defined sphere of center $(1.5,2,2)$ m and radius $0.4$ m and the orb in the lower right. + The primaries should be located in a subset of this orb.} +\end{figure} - \begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{geometrical-output_0000.pdf} - \caption{Geant4 visualisation of a simulation of the primary positions for some geometrical (user defined ) volumes. - The primaries should be located in a sphere of center $(0,2,2)$ m and radius $0.7$ m (center overlapping box), a partially filled cylinder with center - $(-2,-2,-2)$ m and outer radius $0.7$ m and height 1 m (far right) and a box of center $(0,2,-2)$ m and sides of length 0.7, 0.5 and 1.2 m (top). - } - \end{figure} -}{ - \noindent \fbox{\textbf{The test of generating geometrically volume primaries either failed or did not run. - }} -} - -\IfFileExists{geometrical-and-physical-output_0000.pdf}{ - - \begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{geometrical-and-physical-output_0000.pdf} - \caption{Geant4 visualisation of a simulation of the primary positions for the Union of physical volumes defined - by the union of the box and orb (upper left) and the partially filled sphere (lower center) and a geometric volume defined by - sphere of center $(0,2,2)$ m and radius $0.7$ m near the center overlapping the box. - The primaries should be in these three locations.} - \end{figure} -}{ - \noindent \fbox{\textbf{The test of generating geometric or physical volume primaries either failed or did not run. - }} -} -\IfFileExists{geometrical-and-physical-output_0000.pdf}{ - - \begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{geometrical-and-physical-output_0000.pdf} - \caption{Geant4 visualisation of a simulation of the primary positions for the intersection of - physical volumes defined by a user defined sphere of center $(1.5,2,2)$ m and radius $0.4$ m and the orb in the lower right. - The primaries should be located in a subset of this orb.} - \end{figure} -}{ - \noindent \fbox{\textbf{The test of generating geometric and physical volume primaries either failed or did not run. - }} -} \subsection{Union and intersection tests} The next test simulates $2\nu\beta\beta$ decay in a hypothetical HPGe array and makes a statistical comparison of the number of primaries in each volume. -\IfFileExists{relative-ge.output.pdf}{ - \begin{figure}[h!] - \centering - \includegraphics[width=0.9\textwidth]{relative-ge.output.pdf} +\begin{figure}[h!] + \centering + \includegraphics[width=0.9\textwidth]{relative-ge.output.pdf} + \caption{Statistical check on the number of primaries in each HPGe detector.} +\end{figure} - \end{figure} -}{ - \noindent \fbox{\textbf{Relative HPGe sampling test either failed or did not run - }} -} To check the intersection we generate events inside a cylinder around the HPGe detectors and again check the positions and relative fractions. -\IfFileExists{lar-in-check.output.pdf}{ +\begin{figure}[h!] + \centering + \includegraphics[width=0.8\textwidth,page=1]{lar-in-check.output.pdf} + \includegraphics[width=0.8\textwidth,page=2]{lar-in-check.output.pdf} + \includegraphics[width=0.8\textwidth,page=3]{lar-in-check.output.pdf} + \caption{Check on intersections by generating events in a cylinder around each HPGe string.} +\end{figure} - \begin{figure}[h!] - \centering - \includegraphics[width=0.9\textwidth,page=1]{lar-in-check.output.pdf} - \includegraphics[width=0.9\textwidth,page=2]{lar-in-check.output.pdf} - \includegraphics[width=0.9\textwidth,page=3]{lar-in-check.output.pdf} +Finally check subtraction by generating events excluding these regions. +\begin{figure}[h!] + \centering + \includegraphics[width=0.45\textwidth,page=1]{lar-sub-check.output.pdf} + \includegraphics[width=0.45\textwidth,page=1]{lar-int-and-sub-check.output.pdf} + \includegraphics[width=0.45\textwidth,page=2]{lar-sub-check.output.pdf} + \includegraphics[width=0.45\textwidth,page=2]{lar-int-and-sub-check.output.pdf} - \end{figure} -}{ - \noindent \fbox{\textbf{Intersection of physical and geometrical volume test (LAr cylinder) either failed or did not run - }} -} + \caption{Check the vertex positions for a physical volume (the LAr) and subtracted geometrical volumes (cylinders) round each string. + Left figures are just this subtraction while on the right we also intersect with a geometrically (user defined) cylinder the height of the + subtracted cylinders. + } -Finally check subtraction by generating events excluding these regions. +\end{figure} + +\subsection{Generic surface sampling} +The final part of this report describes checks on the generic surface sampler algorithm. +First check the bounding spheres for some solids. -\IfFileExists{lar-sub-check.output.pdf}{ +\begin{figure}[h!] + \centering + \includegraphics[width=0.3\textwidth]{surface-sample-bounding-box-simple.output_0000.png} + \includegraphics[width=0.3\textwidth]{surface-sample-bounding-box-subtraction.output_0000.png} + \includegraphics[width=0.3\textwidth]{surface-sample-bounding-box-union.output_0000.png} - \begin{figure}[h!] - \centering - \includegraphics[width=0.9\textwidth,page=1]{lar-sub-check.output.pdf} - \includegraphics[width=0.9\textwidth,page=2]{lar-sub-check.output.pdf} + \caption{Check the bounding sphere (red) for generating surface events in different volumes, left a cylinder, center the subtraction of two cylinders and left the union of two boxes. We also show the initial points of the lines (used to find intersections). + The bounding sphere should fully contain the solid and the points should be outside (hard to appreciate). + } +\end{figure} - \end{figure} -}{ - \noindent \fbox{\textbf{Subtraction of physical and geometrical volume test (LAr cylinder) either failed or did not run - }} -} -Finally, both an intersection and a subtraction. +\end{multicols} -\IfFileExists{lar-int-and-sub-check.output.pdf}{ +Next, we generate surface events in various shapes and check the location of the primaries. - \begin{figure}[h!] - \centering - \includegraphics[width=0.9\textwidth,page=1]{lar-int-and-sub-check.output.pdf} - \includegraphics[width=0.9\textwidth,page=2]{lar-int-and-sub-check.output.pdf} +% check all the files exist - \end{figure} -}{ - \noindent \fbox{\textbf{Subtraction of physical and geometrical volume test (LAr cylinder) either failed or did not run - }} +\begin{figure}[h!] + \centering + \includegraphics[width=0.45\textwidth]{vis-surface-trd.output_0000.png} + \includegraphics[width=0.45\textwidth]{vis-surface-box.output_0000.png} + \includegraphics[width=0.45\textwidth]{vis-surface-tubby.output_0000.png} + \includegraphics[width=0.45\textwidth]{vis-surface-uni.output_0000.png} + \includegraphics[width=0.45\textwidth]{vis-surface-sub.output_0000.png} + \includegraphics[width=0.45\textwidth]{vis-surface-con.output_0000.png} + \caption{Check the location of primary vertices for points generated on the surface of various shapes. + The primaries should be generated in one solid at a time with the same layout as the figures are arranged. } +\end{figure} -} -\subsection{Generic surface sampling} -The final part of this report describes checks on the generic surface sampler algorithm. +We make a more detailed test in python by extracting the coordinates of generated primaries and checking that they are close to a surface and then finding which. + +% check all the files exist + +\begin{figure}[h!] + \centering + \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-trd.output.pdf} + \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-box.output.pdf} + \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-tubby.output.pdf} + \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-uni.output.pdf} + \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-sub.output.pdf} + \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-con.output.pdf} + \caption{Check of the location of primary vertices for points generated on the surfaces of various shapes. The position of the vertices is plotted in python colouring by face each corresponds to. + You should be able to see that the primaries are generated on the surface and the coloruing should be correct, although the perspective of the 3D plots makes this hard to appreciate. + } +\end{figure} + +% check all the files exist + +\begin{figure}[h!] + \centering + \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-trd.output.pdf} + \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-box.output.pdf} + \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-tubby.output.pdf} + \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-uni.output.pdf} + \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-sub.output.pdf} + \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-con.output.pdf} + + \caption{Check the location of primary vertices for points generated on the surfaces of various shapes, showing 2D projections onto the (x,y), (r,z) and (x,z) planes. The position of the vertices is plotted in python colouring by face each corresponds to. + You should be able to see that the primaries are generated on the surface and the coloring should be correct. + } +\end{figure} + + +Finally, we count the number of events on each face and make a statical test that this is compatible with the surface area. + +\begin{figure}[h!] + \centering + \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-trd.output.pdf} + \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-box.output.pdf} + \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-tubby.output.pdf} + \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-uni.output.pdf} + \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-sub.output.pdf} + \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-con.output.pdf} + + \caption{Fraction of the primary vertices located on each face of the various shapes. This is compared to the expectation + that the fraction is proportional to the surface area. The lower plot shows a residual. These points should be compatible with the expectation with only some statistical scatter. + } +\end{figure} diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index f98421f..b765c0d 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -154,7 +154,7 @@ int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); - UImanager->ApplyCommand("/vis/ogl/export surface-sample-bounding-box-" + name + ".output.pdf"); + UImanager->ApplyCommand("/vis/ogl/export surface-sample-bounding-box-" + name + ".output.png"); delete visManager; diff --git a/tests/tex/validation-report.tex b/tests/tex/validation-report.tex index ecb2787..fbd0c2a 100644 --- a/tests/tex/validation-report.tex +++ b/tests/tex/validation-report.tex @@ -1,4 +1,4 @@ -\documentclass[a4paper,12pt]{article} +\documentclass[a4paper,11pt]{article} % Packages for enhanced functionality \usepackage{amsmath} % For math equations @@ -6,6 +6,7 @@ \usepackage{graphicx} % For including graphics \usepackage{hyperref} % For hyperlinks \usepackage{geometry} % To adjust margins +\usepackage{multicol} % To adjust margins % Adjust margins \geometry{left=1in, right=1in, top=1in, bottom=0.5in} From d8677d10fa1b277e4bd3c4d31b3716e2fed52f7e Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Thu, 16 Jan 2025 12:00:10 +0100 Subject: [PATCH 24/26] cleanup --- Testing/Temporary/LastTest.log | 3 -- include/RMGVertexConfinement.hh | 29 +++++++++------ src/RMGVertexConfinement.cc | 43 +++++++++++----------- tests/CMakeLists.txt | 1 - tests/confinement/CMakeLists.txt | 10 ++++- tests/confinement/template-confinement.tex | 2 +- tests/tex/CMakeLists.txt | 27 -------------- tests/tex/validation-report.tex | 28 -------------- 8 files changed, 48 insertions(+), 95 deletions(-) delete mode 100644 Testing/Temporary/LastTest.log delete mode 100644 tests/tex/CMakeLists.txt delete mode 100644 tests/tex/validation-report.tex diff --git a/Testing/Temporary/LastTest.log b/Testing/Temporary/LastTest.log deleted file mode 100644 index 851b12a..0000000 --- a/Testing/Temporary/LastTest.log +++ /dev/null @@ -1,3 +0,0 @@ -Start testing: Jan 14 00:54 CET ----------------------------------------------------------- -End testing: Jan 14 00:54 CET diff --git a/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index 1b4df80..50b0c2d 100644 --- a/include/RMGVertexConfinement.hh +++ b/include/RMGVertexConfinement.hh @@ -115,26 +115,31 @@ class RMGVertexConfinement : public RMGVVertexGenerator { * An object which we can generate position samples in. Based on either a * @c G4VPhysicalVolume or geometrical volume defined by a @c G4VSolid . The * sampling can be performed either on the surface or in the volume of the solid. + * + * This structure must contain at least a non-null pointer, between the @c physical_volume and + * @c sampling_solid arguments. The idea is that: + * - physical volumes get always a bounding box assigned, but at later time + * - purely geometrical volumes only have the sampling_solid member defined */ struct SampleableObject { - SampleableObject() = default; SampleableObject(const SampleableObject&) = default; /** * @brief SampleableObject constructor. * - * @param v The physical volume. - * @param r A rotation matrix for the sampling solid. - * @param t A translation vector for the sampling solid. - * @param s A solid for geometrical volume sampling or for generating candidate points for - * rejection sampling. - * @param ns A flag of whether the solid is natively sampeable. - * @param ss A flag of whether the solid should be sampled on the surface. + * @param physical_volume The physical volume. + * @param rotation A rotation matrix for the sampling solid. + * @param translation A translation vector for the sampling solid. + * @param sampling_solid A solid for geometrical volume sampling or for generating candidate + * points for rejection sampling. + * @param is_native_sampleable A flag of whether the solid is natively sampeable. + * @param surface_sample A flag of whether the solid should be sampled on the surface. */ - SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r, G4ThreeVector t, G4VSolid* s, - bool ns = false, bool ss = false); + SampleableObject(G4VPhysicalVolume* physical_volume, G4RotationMatrix rotation, + G4ThreeVector translation, G4VSolid* sampling_solid, bool is_native_sampleable = false, + bool surface_sample = false); // NOTE: G4 volume/solid pointers should be fully owned by G4, avoid trying to delete them ~SampleableObject() = default; @@ -201,7 +206,6 @@ class RMGVertexConfinement : public RMGVVertexGenerator { std::vector GetIntersections(const G4ThreeVector start, const G4ThreeVector dir) const; - /** * @brief Get a position and direction for the generic surface sampling algorithm. * @@ -306,7 +310,8 @@ class RMGVertexConfinement : public RMGVVertexGenerator { bool fOnSurface = false; bool fForceContainmentCheck = false; bool fLastSolidExcluded = false; - int fMaxNumIntersections = -1; + int fSurfaceSampleMaxIntersections = -1; + // counters used for the current run. long fTrials = 0; std::chrono::nanoseconds fVertexGenerationTime; diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index 3f39382..92b3522 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -47,19 +47,18 @@ RMGVertexConfinement::SampleableObjectCollection RMGVertexConfinement::fPhysical bool RMGVertexConfinement::fVolumesInitialized = false; -// This structure must contain at least a non-null pointer, between the first -// and the last argument. The idea is that: -// - physical volumes get always a bounding box assigned, but at later time -// - purely geometrical volumes only have the sampling_solid member defined -RMGVertexConfinement::SampleableObject::SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r, - G4ThreeVector t, G4VSolid* s, bool ns, bool ss) - : rotation(r), translation(t), native_sample(ns), surface_sample(ss) { +RMGVertexConfinement::SampleableObject::SampleableObject(G4VPhysicalVolume* physical_volume, + G4RotationMatrix rotation, G4ThreeVector translation, G4VSolid* sampling_solid, + bool is_native_sampleable, bool surface_sample) + : rotation(rotation), translation(translation), native_sample(is_native_sampleable), + surface_sample(surface_sample) { // at least one volume must be specified - if (!v and !s) RMGLog::Out(RMGLog::error, "Invalid pointers given to constructor"); + if (!physical_volume and !sampling_solid) + RMGLog::Out(RMGLog::error, "Invalid pointers given to constructor"); - this->physical_volume = v; - this->sampling_solid = s; + this->physical_volume = physical_volume; + this->sampling_solid = sampling_solid; // should use the physical volume properties, if available const auto& solid = @@ -130,7 +129,8 @@ std::vector RMGVertexConfinement::SampleableObject::GetIntersecti // check if the physical volume exists if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) RMGLog::OutDev(RMGLog::fatal, "Cannot find number of intersections for a SampleableObject ", - "where the physical volume is not set this probably means you are trying to generically ", + "where the physical volume is not set this probably means you are trying to " + "generically ", "sample a geometrical volume which instead should be natively sampled"); auto solid = this->physical_volume->GetLogicalVolume()->GetSolid(); @@ -173,7 +173,8 @@ void RMGVertexConfinement::SampleableObject::GetDirection(G4ThreeVector& dir, if ((not this->physical_volume) or (not this->physical_volume->GetLogicalVolume()->GetSolid())) RMGLog::OutDev(RMGLog::fatal, "Cannot generate directions for a SampleableObject ", - "where the physical volume is not set this probably means you are trying to generically ", + "where the physical volume is not set this probably means you are trying to " + "generically ", "sample a geometrical volume which instead should be natively sampled"); // Get the center and radius of a bounding sphere around the shape @@ -444,13 +445,12 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { el.native_sample = false; el.surface_sample = true; - el.max_num_intersections = fMaxNumIntersections; + el.max_num_intersections = fSurfaceSampleMaxIntersections; - if (fMaxNumIntersections < 2) - RMGLog::Out(RMGLog::fatal, - " for generic surface sampling MaxNumIntersections, the maximum number of lines a ", + if (fSurfaceSampleMaxIntersections < 2) + RMGLog::Out(RMGLog::fatal, " for generic surface sampling SurfaceSampleMaxIntersections, the maximum number of lines a ", "line can intersect with the surface must be set with " - "/RMG/Generator/Confinement/MaxNumberOfIntersections", + "/RMG/Generator/Confinement/SurfaceSampleMaxIntersections", "Note: this can be an overestimate."); } @@ -484,8 +484,8 @@ void RMGVertexConfinement::InitializePhysicalVolumes() { RMGLog::OutDev(RMGLog::debug, "Bounding box coordinates: min = ", lim_min, ", max = ", lim_max); - // the origin of the local coordinates of the non-sampleable solid is not necessarily at its - // barycenter. However, the coordinate origin of a G4Box is always its barycenter. + // the origin of the local coordinates of the non-sampleable solid is not necessarily at + // its barycenter. However, the coordinate origin of a G4Box is always its barycenter. double bb_x = std::max(std::abs(lim_max.getX()), std::abs(lim_min.getX())); double bb_y = std::max(std::abs(lim_max.getY()), std::abs(lim_min.getY())); double bb_z = std::max(std::abs(lim_max.getZ()), std::abs(lim_min.getZ())); @@ -708,7 +708,8 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) { : fGeomVolumeSolids.SurfaceWeightedRand(); } else { - // for volume sampling the user can specify the volume to sample first else the set with smaller total volume is used + // 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; @@ -1010,7 +1011,7 @@ void RMGVertexConfinement::DefineCommands() { fMessengers.back() - ->DeclareProperty("MaxNumberOfIntersections", fMaxNumIntersections) + ->DeclareProperty("SurfaceSampleMaxIntersections", fSurfaceSampleMaxIntersections) .SetGuidance("Set maximum number of intersections of a line with the surface. Note: can be " "set to an overestimate. ") .SetParameterName("N", false) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 90ba5ec..4b1ce4f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -14,4 +14,3 @@ add_subdirectory(internals) add_subdirectory(output) add_subdirectory(python) add_subdirectory(vertex) -add_subdirectory(tex) diff --git a/tests/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 663f546..854ca9b 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -9,8 +9,14 @@ foreach(_file ${_aux}) configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) endforeach() -set(_macros complex-volume.mac complex-surface.mac geometrical.mac native-surface.mac geometrical-and-physical.mac - geometrical-or-physical.mac native-volume.mac) +set(_macros + complex-volume.mac + complex-surface.mac + geometrical.mac + native-surface.mac + geometrical-and-physical.mac + geometrical-or-physical.mac + native-volume.mac) foreach(_mac ${_macros}) add_test(NAME confinement/${_mac} COMMAND ${REMAGE_PYEXE} -g gdml/geometry.gdml -o test-out.root diff --git a/tests/confinement/template-confinement.tex b/tests/confinement/template-confinement.tex index c57873e..76ecc6b 100644 --- a/tests/confinement/template-confinement.tex +++ b/tests/confinement/template-confinement.tex @@ -94,7 +94,7 @@ \subsection{Union and intersection tests} \includegraphics[width=0.45\textwidth,page=2]{lar-int-and-sub-check.output.pdf} \caption{Check the vertex positions for a physical volume (the LAr) and subtracted geometrical volumes (cylinders) round each string. - Left figures are just this subtraction while on the right we also intersect with a geometrically (user defined) cylinder the height of the + Left figures are just this subtraction while on the right we also intersect with a geometrically (user defined) cylinder the height of the subtracted cylinders. } diff --git a/tests/tex/CMakeLists.txt b/tests/tex/CMakeLists.txt deleted file mode 100644 index 31afa1d..0000000 --- a/tests/tex/CMakeLists.txt +++ /dev/null @@ -1,27 +0,0 @@ -file( - GLOB _aux - RELATIVE ${PROJECT_SOURCE_DIR} - *.tex) - -# copy them to the build area -foreach(_file ${_aux}) - configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) -endforeach() - -# copy all files -add_test( - NAME tex/copy-files - COMMAND - /bin/sh -c - "find ../ -type f \\( -name '*.output*' -o -name '*.tex' \\) ! -path '../tex/*' -exec cp {} . \; " -) -set_tests_properties(tex/copy-files PROPERTIES LABELS tex FIXTURES_SETUP copy-files-fixture) - -find_program(PDFLATEX_COMPILER pdflatex) -if(PDFLATEX_COMPILER) - add_test(NAME tex/compile COMMAND ${PDFLATEX_COMPILER} validation-report.tex) - set_tests_properties(tex/copy-files PROPERTIES LABELS tex FIXTURES_REQUIRED copy-files-fxiture) - -else() - message(STATUS "No LaTeX compiler found, report cannot be compiled.") -endif() diff --git a/tests/tex/validation-report.tex b/tests/tex/validation-report.tex deleted file mode 100644 index fbd0c2a..0000000 --- a/tests/tex/validation-report.tex +++ /dev/null @@ -1,28 +0,0 @@ -\documentclass[a4paper,11pt]{article} - -% Packages for enhanced functionality -\usepackage{amsmath} % For math equations -\usepackage{amsfonts} % For special fonts -\usepackage{graphicx} % For including graphics -\usepackage{hyperref} % For hyperlinks -\usepackage{geometry} % To adjust margins -\usepackage{multicol} % To adjust margins - -% Adjust margins -\geometry{left=1in, right=1in, top=1in, bottom=0.5in} - -\title{{\Large \bf remage validation suite.}} -\author{{\normalsize Citation: \url{https://doi.org/10.5281/zenodo.11115662}}} -\date{{\it\normalsize \today}} % Date is automatically set to the current date - -\begin{document} -\maketitle - -\section{Basics} -\input{template-basic} - -\section{Confinement} -\input{template-confinement} - - -\end{document} From 13d2dc510915d57e2dd3b4ad172fb145d6bf1d51 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Thu, 16 Jan 2025 12:43:26 +0100 Subject: [PATCH 25/26] cleanup and removing the tex --- .github/workflows/main.yml | 6 +- tests/basics/template-basic.tex | 12 -- tests/confinement/macros/_vis-export.mac | 2 +- tests/confinement/macros/_vis.mac | 2 +- tests/confinement/macros/complex-surface.mac | 2 +- tests/confinement/macros/gen-surface.mac | 2 +- .../macros/geometrical-and-physical.mac | 2 +- .../macros/geometrical-or-physical.mac | 2 +- tests/confinement/macros/geometrical.mac | 2 +- tests/confinement/macros/native-surface.mac | 2 +- tests/confinement/macros/vis-surface.mac | 4 +- tests/confinement/template-confinement.tex | 186 ------------------ .../test-surface-sampler-methods.cpp | 7 +- 13 files changed, 14 insertions(+), 217 deletions(-) delete mode 100644 tests/basics/template-basic.tex delete mode 100644 tests/confinement/template-confinement.tex diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index b5d78a6..abfc504 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -58,10 +58,8 @@ jobs: uses: actions/upload-artifact@v4 with: name: test-output-${{ matrix.container_version }} - # artifacts must have a ".output*" extension or .tex - path: | - build/tests/tex/*.output*.* - build/tests/tex/*.tex + # artifacts must have a ".output*" extension + path: build/tests/**/*.output*.* - name: Run minimal test suite if: ${{ matrix.container_version == 'slim' }} diff --git a/tests/basics/template-basic.tex b/tests/basics/template-basic.tex deleted file mode 100644 index 00fb37c..0000000 --- a/tests/basics/template-basic.tex +++ /dev/null @@ -1,12 +0,0 @@ -The first set of checks shown below test that remage is running and has reasonable seeming results. - -\begin{figure}[h!] - \centering - \includegraphics[width=0.45\textwidth]{vis-2nbb.output_0000.pdf} - \includegraphics[width=0.45\textwidth]{vis-co60.output_0000.pdf} - - \caption{Geant4 visualisation of a simulation of $2\nu\beta\beta$ decay in a HPGe detector (left) and $^{60}$Co in a source above the detector (right( for a typical screeing \ - setup. For the left figure you should see the trajectories focused around the HPGe detector and mainly electrons (red), with a \ - few gamma (green) being produced. For the right you should see $\gamma$ particles (in green) throughout the experimental setup starting from a source above the detector - } -\end{figure} diff --git a/tests/confinement/macros/_vis-export.mac b/tests/confinement/macros/_vis-export.mac index 5f5962b..07bd5aa 100644 --- a/tests/confinement/macros/_vis-export.mac +++ b/tests/confinement/macros/_vis-export.mac @@ -1,2 +1,2 @@ -/vis/ogl/set/printSize 300 300 +/vis/ogl/set/printMode pixmap /vis/ogl/export {export-fn}.pdf diff --git a/tests/confinement/macros/_vis.mac b/tests/confinement/macros/_vis.mac index a0e3be6..0e6910d 100644 --- a/tests/confinement/macros/_vis.mac +++ b/tests/confinement/macros/_vis.mac @@ -2,7 +2,7 @@ /RMG/Manager/Logging/LogLevel summary -/vis/open OGL 600x600-0+0 +/vis/open OGL /vis/drawVolume /vis/sceneHandler/attach diff --git a/tests/confinement/macros/complex-surface.mac b/tests/confinement/macros/complex-surface.mac index aa4d68c..9bf6b8d 100644 --- a/tests/confinement/macros/complex-surface.mac +++ b/tests/confinement/macros/complex-surface.mac @@ -3,7 +3,7 @@ /RMG/Generator/Confine Volume /RMG/Generator/Confinement/ForceContainmentCheck true /RMG/Generator/Confinement/SampleOnSurface true -/RMG/Generator/Confinement/MaxNumberOfIntersections 6 +/RMG/Generator/Confinement/SurfaceSampleMaxIntersections 6 /RMG/Generator/Confinement/Physical/AddVolume BoxWithHole /RMG/Generator/Confinement/Physical/AddVolume BoxAndOrb diff --git a/tests/confinement/macros/gen-surface.mac b/tests/confinement/macros/gen-surface.mac index 5ad999f..1a97fdf 100644 --- a/tests/confinement/macros/gen-surface.mac +++ b/tests/confinement/macros/gen-surface.mac @@ -12,7 +12,7 @@ /RMG/Generator/Confinement/Physical/AddVolume ${det} -/RMG/Generator/Confinement/MaxNumberOfIntersections 6 +/RMG/Generator/Confinement/SurfaceSampleMaxIntersections 6 /gps/particle e- /gps/ang/type iso diff --git a/tests/confinement/macros/geometrical-and-physical.mac b/tests/confinement/macros/geometrical-and-physical.mac index 06abf86..d5403ec 100644 --- a/tests/confinement/macros/geometrical-and-physical.mac +++ b/tests/confinement/macros/geometrical-and-physical.mac @@ -20,4 +20,4 @@ /run/beamOn 2000 -/control/alias export-fn geometrical-and-physical.output.jpg +/control/alias export-fn geometrical-and-physical.output diff --git a/tests/confinement/macros/geometrical-or-physical.mac b/tests/confinement/macros/geometrical-or-physical.mac index 60932ef..78d7417 100644 --- a/tests/confinement/macros/geometrical-or-physical.mac +++ b/tests/confinement/macros/geometrical-or-physical.mac @@ -21,4 +21,4 @@ /run/beamOn 2000 -/control/alias export-fn geometrical-or-physical.output.jpg +/control/alias export-fn geometrical-or-physical.output diff --git a/tests/confinement/macros/geometrical.mac b/tests/confinement/macros/geometrical.mac index f504d5b..46f6daa 100644 --- a/tests/confinement/macros/geometrical.mac +++ b/tests/confinement/macros/geometrical.mac @@ -36,4 +36,4 @@ /run/beamOn 5000 -/control/alias export-fn geometrical.output.jpg +/control/alias export-fn geometrical.output diff --git a/tests/confinement/macros/native-surface.mac b/tests/confinement/macros/native-surface.mac index 979d570..aad2e7e 100644 --- a/tests/confinement/macros/native-surface.mac +++ b/tests/confinement/macros/native-surface.mac @@ -13,4 +13,4 @@ /run/beamOn 2000 -/control/alias export-fn native-surface.output.pdf +/control/alias export-fn native-surface.output diff --git a/tests/confinement/macros/vis-surface.mac b/tests/confinement/macros/vis-surface.mac index 0a88f41..9ad8d69 100644 --- a/tests/confinement/macros/vis-surface.mac +++ b/tests/confinement/macros/vis-surface.mac @@ -48,7 +48,7 @@ /RMG/Generator/Confinement/Physical/AddVolume ${det} -/RMG/Generator/Confinement/MaxNumberOfIntersections 6 +/RMG/Generator/Confinement/SurfaceSampleMaxIntersections 6 /gps/particle e- /gps/ang/type iso @@ -56,4 +56,4 @@ /gps/energy 1 eV /run/beamOn 1000 -/control/alias export-fn vis-surface-${det}.output.png +/control/alias export-fn vis-surface-${det}.output diff --git a/tests/confinement/template-confinement.tex b/tests/confinement/template-confinement.tex deleted file mode 100644 index 76ecc6b..0000000 --- a/tests/confinement/template-confinement.tex +++ /dev/null @@ -1,186 +0,0 @@ - -These tests check the generation of primary vertices by remage. -\subsection{Native sampling} -The first checks plot the positions for natively sampleable objects. - - - -\begin{figure}[h!] - \centering - \includegraphics[width=0.45\textwidth]{native-surface.output_0000.png} - \includegraphics[width=0.45\textwidth]{native-volume.output_0000.png} - \caption{Geant4 visualisation of a simulation of the primary positions for the surface (left) and volume (right) of - some natively sampled shapes. Primaries should be present on the G4Box (bottom row second from right), - G4Orb (full sphere - above the box) and G4Tubs (left). The primaries should be on the detector surface or volume, although this - is hard to appreciate with Geant4 visualisation. - } -\end{figure} - - -Finally, we check the generation of points on complex objects. -\begin{figure}[h!] - \centering - \includegraphics[width=0.45\textwidth]{complex-surface.output_0000.png} - \includegraphics[width=0.45\textwidth]{complex-volume.output_0000.png} - - \caption{Geant4 visualisation of a simulation of the primary positions for the surface (left) or volume (right) of - some complex shapes. A polycone (top center), a box with a hole (center), the union of an orb and box (upper right) - and an orb (right). - } -\end{figure} - - - -\subsection{Geometrical volumes and intersections} -The next tests check the generation of points in geometrical (user defined) volumes and intersections / unions of -geometrical and physical volumes. - -\begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{geometrical.output_0000.png} - \caption{Geant4 visualisation of a simulation of the primary positions for some geometrical (user defined ) volumes. - The primaries should be located in a sphere of center $(0,2,2)$ m and radius $0.7$ m (center overlapping box), a partially filled cylinder with center - $(-2,-2,-2)$ m and outer radius $0.7$ m and height 1 m (far right) and a box of center $(0,2,-2)$ m and sides of length 0.7, 0.5 and 1.2 m (top). - } -\end{figure} - -\begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{geometrical-or-physical.output_0000.png} - \caption{Geant4 visualisation of a simulation of the primary positions for the Union of physical volumes defined - by the union of the box and orb (upper left) and the partially filled sphere (lower center) and a geometric volume defined by - sphere of center $(0,2,2)$ m and radius $0.7$ m near the center overlapping the box. - The primaries should be in these three locations.} -\end{figure} - - -\begin{figure}[h!] - \centering - \includegraphics[width=0.6\textwidth]{geometrical-and-physical.output_0000.png} - \caption{Geant4 visualisation of a simulation of the primary positions for the intersection of - physical volumes defined by a user defined sphere of center $(1.5,2,2)$ m and radius $0.4$ m and the orb in the lower right. - The primaries should be located in a subset of this orb.} -\end{figure} - - -\subsection{Union and intersection tests} -The next test simulates $2\nu\beta\beta$ decay in a hypothetical HPGe array and makes a -statistical comparison of the number of primaries in each volume. - -\begin{figure}[h!] - \centering - \includegraphics[width=0.9\textwidth]{relative-ge.output.pdf} - \caption{Statistical check on the number of primaries in each HPGe detector.} -\end{figure} - - -To check the intersection we generate events inside a cylinder around the HPGe detectors and again check the -positions and relative fractions. - -\begin{figure}[h!] - \centering - \includegraphics[width=0.8\textwidth,page=1]{lar-in-check.output.pdf} - \includegraphics[width=0.8\textwidth,page=2]{lar-in-check.output.pdf} - \includegraphics[width=0.8\textwidth,page=3]{lar-in-check.output.pdf} - \caption{Check on intersections by generating events in a cylinder around each HPGe string.} -\end{figure} - -Finally check subtraction by generating events excluding these regions. -\begin{figure}[h!] - \centering - \includegraphics[width=0.45\textwidth,page=1]{lar-sub-check.output.pdf} - \includegraphics[width=0.45\textwidth,page=1]{lar-int-and-sub-check.output.pdf} - \includegraphics[width=0.45\textwidth,page=2]{lar-sub-check.output.pdf} - \includegraphics[width=0.45\textwidth,page=2]{lar-int-and-sub-check.output.pdf} - - \caption{Check the vertex positions for a physical volume (the LAr) and subtracted geometrical volumes (cylinders) round each string. - Left figures are just this subtraction while on the right we also intersect with a geometrically (user defined) cylinder the height of the - subtracted cylinders. - } - -\end{figure} - -\subsection{Generic surface sampling} -The final part of this report describes checks on the generic surface sampler algorithm. -First check the bounding spheres for some solids. - -\begin{figure}[h!] - \centering - \includegraphics[width=0.3\textwidth]{surface-sample-bounding-box-simple.output_0000.png} - \includegraphics[width=0.3\textwidth]{surface-sample-bounding-box-subtraction.output_0000.png} - \includegraphics[width=0.3\textwidth]{surface-sample-bounding-box-union.output_0000.png} - - \caption{Check the bounding sphere (red) for generating surface events in different volumes, left a cylinder, center the subtraction of two cylinders and left the union of two boxes. We also show the initial points of the lines (used to find intersections). - The bounding sphere should fully contain the solid and the points should be outside (hard to appreciate). - } -\end{figure} - - -\end{multicols} - -Next, we generate surface events in various shapes and check the location of the primaries. - -% check all the files exist - -\begin{figure}[h!] - \centering - \includegraphics[width=0.45\textwidth]{vis-surface-trd.output_0000.png} - \includegraphics[width=0.45\textwidth]{vis-surface-box.output_0000.png} - \includegraphics[width=0.45\textwidth]{vis-surface-tubby.output_0000.png} - \includegraphics[width=0.45\textwidth]{vis-surface-uni.output_0000.png} - \includegraphics[width=0.45\textwidth]{vis-surface-sub.output_0000.png} - \includegraphics[width=0.45\textwidth]{vis-surface-con.output_0000.png} - \caption{Check the location of primary vertices for points generated on the surface of various shapes. - The primaries should be generated in one solid at a time with the same layout as the figures are arranged. } -\end{figure} - - -We make a more detailed test in python by extracting the coordinates of generated primaries and checking that they are close to a surface and then finding which. - -% check all the files exist - -\begin{figure}[h!] - \centering - \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-trd.output.pdf} - \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-box.output.pdf} - \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-tubby.output.pdf} - \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-uni.output.pdf} - \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-sub.output.pdf} - \includegraphics[width=0.45\textwidth,page=1]{confinement.simple-solids-surface-con.output.pdf} - \caption{Check of the location of primary vertices for points generated on the surfaces of various shapes. The position of the vertices is plotted in python colouring by face each corresponds to. - You should be able to see that the primaries are generated on the surface and the coloruing should be correct, although the perspective of the 3D plots makes this hard to appreciate. - } -\end{figure} - -% check all the files exist - -\begin{figure}[h!] - \centering - \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-trd.output.pdf} - \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-box.output.pdf} - \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-tubby.output.pdf} - \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-uni.output.pdf} - \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-sub.output.pdf} - \includegraphics[width=0.45\textwidth,page=2]{confinement.simple-solids-surface-con.output.pdf} - - \caption{Check the location of primary vertices for points generated on the surfaces of various shapes, showing 2D projections onto the (x,y), (r,z) and (x,z) planes. The position of the vertices is plotted in python colouring by face each corresponds to. - You should be able to see that the primaries are generated on the surface and the coloring should be correct. - } -\end{figure} - - -Finally, we count the number of events on each face and make a statical test that this is compatible with the surface area. - -\begin{figure}[h!] - \centering - \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-trd.output.pdf} - \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-box.output.pdf} - \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-tubby.output.pdf} - \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-uni.output.pdf} - \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-sub.output.pdf} - \includegraphics[width=0.45\textwidth,page=3]{confinement.simple-solids-surface-con.output.pdf} - - \caption{Fraction of the primary vertices located on each face of the various shapes. This is compared to the expectation - that the fraction is proportional to the surface area. The lower plot shows a residual. These points should be compatible with the expectation with only some statistical scatter. - } -\end{figure} diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index b765c0d..cf32bf8 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -140,21 +140,18 @@ int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { // Set up visualization commands G4UImanager* UImanager = G4UImanager::GetUIpointer(); UImanager->ApplyCommand("/run/initialize "); - UImanager->ApplyCommand("/vis/open OGL "); UImanager->ApplyCommand("/vis/verbose warnings"); - UImanager->ApplyCommand("/vis/drawVolume"); UImanager->ApplyCommand("/vis/geometry/set/forceWireframe"); UImanager->ApplyCommand("/vis/viewer/set/viewpointVector 0.7 0.9 0.7"); UImanager->ApplyCommand("/vis/viewer/zoom 1.5"); - UImanager->ApplyCommand("/vis/scene/list "); UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate"); UImanager->ApplyCommand("/vis/viewer/set/globalLineWidthScale 1.5"); UImanager->ApplyCommand("/vis/viewer/set/upVector 0 0 1"); - - UImanager->ApplyCommand("/vis/ogl/export surface-sample-bounding-box-" + name + ".output.png"); + UImanager->ApplyCommand("/vis/ogl/set/printMode pixmap"); + UImanager->ApplyCommand("/vis/ogl/export surface-sample-bounding-box-" + name + ".output.pdf"); delete visManager; From 5dfe3078a559356dd62b6b384b2d02d061adff90 Mon Sep 17 00:00:00 2001 From: Toby Dixon Date: Thu, 16 Jan 2025 12:56:46 +0100 Subject: [PATCH 26/26] [docs] update the macro commands file --- docs/rmg-commands.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index a9db9dd..9b24129 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -958,7 +958,7 @@ Commands for controlling primary confinement * `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 -* `MaxNumberOfIntersections` – Set maximum number of intersections of a line with the surface. Note: can be set to an overestimate. +* `SurfaceSampleMaxIntersections` – Set maximum number of intersections of a line with the surface. Note: can be set to an overestimate. * `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. ### `/RMG/Generator/Confinement/Reset` @@ -1000,7 +1000,7 @@ Set maximum number of attempts for sampling primary positions in a volume * **Parameter type** – `i` * **Omittable** – `False` -### `/RMG/Generator/Confinement/MaxNumberOfIntersections` +### `/RMG/Generator/Confinement/SurfaceSampleMaxIntersections` Set maximum number of intersections of a line with the surface. Note: can be set to an overestimate.