diff --git a/docs/developer.md b/docs/developer.md index 0be3fd07..a82aca04 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/docs/rmg-commands.md b/docs/rmg-commands.md index 79bb4e2d..9b241296 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 +* `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` @@ -999,6 +1000,15 @@ Set maximum number of attempts for sampling primary positions in a volume * **Parameter type** – `i` * **Omittable** – `False` +### `/RMG/Generator/Confinement/SurfaceSampleMaxIntersections` + +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/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index 48d56227..50b0c2d1 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,22 +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 @@ -91,33 +111,144 @@ 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. + * + * 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; - SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r, G4ThreeVector t, G4VSolid* s, - bool cc = true); + + /** + * @brief SampleableObject constructor. + * + * @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* 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; + + /** + * @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; - [[nodiscard]] bool Sample(G4ThreeVector& vertex, int max_attempts, bool sample_on_surface, + + /** + * @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. + * - 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. + * @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; + G4VPhysicalVolume* physical_volume = nullptr; G4VSolid* sampling_solid = nullptr; G4RotationMatrix rotation; 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; }; + /** 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; @@ -179,6 +310,8 @@ class RMGVertexConfinement : public RMGVVertexGenerator { bool fOnSurface = false; bool fForceContainmentCheck = false; bool fLastSolidExcluded = false; + 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 91adc083..92b3522a 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 cc) - : rotation(r), translation(t), containment_check(cc) { +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 = @@ -121,8 +120,134 @@ 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; + + // 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) + : solid->DistanceToOut(int_point, dir); + + int_point = int_point + dist * dir; + + if (dist < kInfinity) intersections.push_back(int_point); + counter++; + } + + if (intersections.size() % 2) + RMGLog::OutDev(RMGLog::fatal, "Found ", intersections.size(), " intersections. However,", + "the number should always be divisible by two."); + + 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 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 center and radius of a bounding sphere around the shape + 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); + dir = G4ThreeVector(0.0, 0.0, -1.0); + + // 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 + 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); + + // shift by the barycenter of the bounding sphere. + pos += barycenter; +} + + +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; + + // 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) { + 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 +260,39 @@ 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, 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); 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,15 +302,6 @@ 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, @@ -169,7 +309,6 @@ bool RMGVertexConfinement::SampleableObject::Sample(G4ThreeVector& vertex, int m return true; } - bool RMGVertexConfinement::SampleableObjectCollection::IsInside(const G4ThreeVector& vertex) const { auto navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); for (const auto& o : data) { @@ -282,16 +421,38 @@ 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 = fSurfaceSampleMaxIntersections; + + 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/SurfaceSampleMaxIntersections", + "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 +461,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 @@ -322,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())); @@ -331,7 +493,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 +563,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 +608,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 '", @@ -545,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; @@ -557,8 +721,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 +807,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 +863,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 +1009,16 @@ void RMGVertexConfinement::DefineCommands() { .SetStates(G4State_PreInit, G4State_Idle) .SetToBeBroadcasted(true); + + fMessengers.back() + ->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) + .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/basics/CMakeLists.txt b/tests/basics/CMakeLists.txt index f30b7f8d..b701cffe 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/confinement/CMakeLists.txt b/tests/confinement/CMakeLists.txt index 9acd21e0..854ca9bc 100644 --- a/tests/confinement/CMakeLists.txt +++ b/tests/confinement/CMakeLists.txt @@ -2,16 +2,21 @@ 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}) 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 - 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 @@ -93,4 +98,97 @@ 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) +target_link_libraries(test-surface-sampler-methods PUBLIC remage) +if(BxDecay0_FOUND) + 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) + +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() + +# 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-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-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 + gdml/simple-solids.gdml -w -l detail) + + set_tests_properties( + confinment-surface/${det}-vis PROPERTIES LABELS extra FIXTURES_REQUIRED + 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) + + add_test(NAME confinment-surface/relative-fractions-${det} + COMMAND ${PYTHONPATH} ./test_basic_surface.py ${det}) + + set_tests_properties(confinment-surface/relative-fractions-${det} + PROPERTIES LABELS extra FIXTURES_REQUIRED confine-surface-${det}) + +endforeach() diff --git a/tests/confinement/macros/_vis-export.mac b/tests/confinement/macros/_vis-export.mac index 7b58471e..07bd5aa1 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/printMode pixmap +/vis/ogl/export {export-fn}.pdf diff --git a/tests/confinement/macros/_vis.mac b/tests/confinement/macros/_vis.mac index 894be177..0e6910db 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/macros/complex-surface.mac b/tests/confinement/macros/complex-surface.mac new file mode 100644 index 00000000..9bf6b8d5 --- /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/SurfaceSampleMaxIntersections 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 cd63da15..964d1926 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/gen-surface.mac b/tests/confinement/macros/gen-surface.mac new file mode 100644 index 00000000..1a97fdfc --- /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/SurfaceSampleMaxIntersections 6 + +/gps/particle e- +/gps/ang/type iso + +/gps/energy 1 eV +/run/beamOn 50000 diff --git a/tests/confinement/macros/geometrical-and-physical.mac b/tests/confinement/macros/geometrical-and-physical.mac index c19dd3c3..d5403ecc 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 diff --git a/tests/confinement/macros/geometrical-or-physical.mac b/tests/confinement/macros/geometrical-or-physical.mac index 8b6340f7..78d74179 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 diff --git a/tests/confinement/macros/geometrical.mac b/tests/confinement/macros/geometrical.mac index e14805b2..46f6daa0 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 diff --git a/tests/confinement/macros/native-surface.mac b/tests/confinement/macros/native-surface.mac index 979d570a..aad2e7e9 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/native-volume.mac b/tests/confinement/macros/native-volume.mac index b5804af0..5c0dcac3 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 new file mode 100644 index 00000000..9ad8d692 --- /dev/null +++ b/tests/confinement/macros/vis-surface.mac @@ -0,0 +1,59 @@ +/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/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/SurfaceSampleMaxIntersections 6 + +/gps/particle e- +/gps/ang/type iso + +/gps/energy 1 eV +/run/beamOn 1000 + +/control/alias export-fn vis-surface-${det}.output diff --git a/tests/confinement/make_basic_solids.py b/tests/confinement/make_basic_solids.py new file mode 100644 index 00000000..775b7186 --- /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, 75]], 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 new file mode 100644 index 00000000..cf32bf88 --- /dev/null +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -0,0 +1,407 @@ +#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 "G4SystemOfUnits.hh" +#include "G4ThreeVector.hh" +#include "G4Tubs.hh" +#include "G4UIExecutive.hh" +#include "G4UImanager.hh" +#include "G4UnionSolid.hh" +#include "G4VSolid.hh" +#include "G4VUserDetectorConstruction.hh" +#include "G4VisAttributes.hh" +#include "G4VisExecutive.hh" +#include "Randomize.hh" + +#include "RMGVertexConfinement.hh" + +#include "QBBC.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(false); + 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; + } +}; + + +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/set/printMode pixmap"); + UImanager->ApplyCommand("/vis/ogl/export surface-sample-bounding-box-" + 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); + 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]; + + // define the solids we need + 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); + + 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 = + 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") { + + auto obj = sampleables["sub"]; + // 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 + + 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 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-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 (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 << " " << pos + << std::endl; + return 1; + } + i++; + } + } else if (test_type == "test-containment-subtraction") { + auto obj = sampleables["sub"]; + + int i = 0; + while (i < 10000) { + G4ThreeVector pos; + bool success = obj.GenerateSurfacePoint(pos, 200, 4); + + 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 << std::endl; + return 1; + } + i++; + } + } else if (test_type == "test-points-union") { + + // get some points to plot + auto obj = sampleables["uni"]; + RunVis(obj, "union"); + + } else if (test_type == "test-points-subtraction") { + + // get some points to plot + auto obj = sampleables["sub"]; + return RunVis(obj, "subtraction"); + + } else if (test_type == "test-points-basic") { + + // get some points to plot + auto obj = sampleables["tubs"]; + return RunVis(obj, "simple"); + + } + + else { + return 0; + } + + return 0; +} diff --git a/tests/confinement/test_basic_surface.py b/tests/confinement/test_basic_surface.py new file mode 100644 index 00000000..259aab62 --- /dev/null +++ b/tests/confinement/test_basic_surface.py @@ -0,0 +1,368 @@ +from __future__ import annotations + +import argparse +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 + + +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() + + +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 - 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, + 100 * 50.0, + 100 * 50.0, + 100 * 50.0, + 100 * 50.0, + 50 * 50.0 - 20 * 20.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", + }, +} +dtype = args.det + + +with PdfPages(f"confinement.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 + + if len(indices[indices == -1]) > 0: + 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 + 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() + + assert sigma < 5 diff --git a/tests/confinement/test_confinment_ge.py b/tests/confinement/test_confinment_ge.py index 82feedd1..e01f4f2b 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 diff --git a/tests/germanium/CMakeLists.txt b/tests/germanium/CMakeLists.txt index 9fde7cab..6018bf27 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})