diff --git a/include/RMGVertexConfinement.hh b/include/RMGVertexConfinement.hh index 50b0c2d..f45c6b6 100644 --- a/include/RMGVertexConfinement.hh +++ b/include/RMGVertexConfinement.hh @@ -49,7 +49,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator { /** @brief Information about the geometrical (user) defined solids. */ struct GenericGeometricalSolidData { - GeometricalSolidType solid_type; + GeometricalSolidType solid_type = GeometricalSolidType::kBox; G4ThreeVector volume_center = G4ThreeVector(0, 0, 0); double sphere_inner_radius = 0; double sphere_outer_radius = -1; @@ -180,12 +180,12 @@ class RMGVertexConfinement : public RMGVVertexGenerator { * - Pick one intersection, or repeat. * * @param vertex The sampled vertex, - * @param max_samples The maximum number of attempts to find a valid vertex. + * @param max_attempts 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, + [[nodiscard]] bool GenerateSurfacePoint(G4ThreeVector& vertex, int max_attempts, int n_max) const; // methods for the generic surface sampling @@ -203,8 +203,8 @@ class RMGVertexConfinement : public RMGVVertexGenerator { * * @returns A vector of the points of intersection. */ - std::vector GetIntersections(const G4ThreeVector start, - const G4ThreeVector dir) const; + [[nodiscard]] std::vector GetIntersections(G4ThreeVector start, + G4ThreeVector dir) const; /** * @brief Get a position and direction for the generic surface sampling algorithm. @@ -284,7 +284,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator { }; void InitializePhysicalVolumes(); - void InitializeGeometricalVolumes(bool useExcludedVolumes); + void InitializeGeometricalVolumes(bool use_excluded_volumes); bool ActualGenerateVertex(G4ThreeVector& v); std::vector fPhysicalVolumeNameRegexes; diff --git a/src/RMGVertexConfinement.cc b/src/RMGVertexConfinement.cc index 92b3522..e784341 100644 --- a/src/RMGVertexConfinement.cc +++ b/src/RMGVertexConfinement.cc @@ -222,7 +222,7 @@ bool RMGVertexConfinement::SampleableObject::GenerateSurfacePoint(G4ThreeVector& // find the intersections std::vector intersections = this->GetIntersections(pos, dir); - if (intersections.size() == 0) continue; + if (intersections.empty()) continue; // The surface sampling algorithm returns N intersections. // We have to select one, to keep independence of the sampled points @@ -310,11 +310,8 @@ bool RMGVertexConfinement::SampleableObject::Sample(G4ThreeVector& vertex, int m } bool RMGVertexConfinement::SampleableObjectCollection::IsInside(const G4ThreeVector& vertex) const { - auto navigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); - for (const auto& o : data) { - if (o.IsInside(vertex)) return true; - } - return false; + return std::any_of(data.begin(), data.end(), + [&vertex](const auto& o) { return o.IsInside(vertex); }); } diff --git a/tests/confinement/test-surface-sampler-methods.cpp b/tests/confinement/test-surface-sampler-methods.cpp index cf32bf8..0fae96c 100644 --- a/tests/confinement/test-surface-sampler-methods.cpp +++ b/tests/confinement/test-surface-sampler-methods.cpp @@ -5,7 +5,6 @@ #include "G4Box.hh" #include "G4LogicalVolume.hh" -#include "G4Material.hh" #include "G4NistManager.hh" #include "G4Orb.hh" #include "G4PVPlacement.hh" @@ -15,7 +14,6 @@ #include "G4SystemOfUnits.hh" #include "G4ThreeVector.hh" #include "G4Tubs.hh" -#include "G4UIExecutive.hh" #include "G4UImanager.hh" #include "G4UnionSolid.hh" #include "G4VSolid.hh" @@ -47,51 +45,45 @@ class MyDetectorConstruction : public G4VUserDetectorConstruction { G4VPhysicalVolume* Construct() override { // Define materials - G4Material* vacuum = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR"); + auto 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"); + auto worldBox = new G4Box("World", 200 * mm, 200 * mm, 200 * mm); + auto 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); + auto worldPhys = new G4PVPlacement(nullptr, {}, worldLog, "World", nullptr, false, 0); // Sampling sphere - G4Orb* orb = new G4Orb("MyOrb", fRadius); - G4LogicalVolume* orbLog = new G4LogicalVolume(orb, vacuum, "MyOrb"); - + auto orb = new G4Orb("MyOrb", fRadius); + auto 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 + auto 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); + auto 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"); + auto pointSphere = new G4Sphere("Point", 0, 1 * mm, 0, 2 * M_PI, 0, M_PI); + auto 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 + auto pointVisAttr = new G4VisAttributes(G4Colour(0, 0, 1)); // Blue pointVisAttr->SetVisibility(true); - pointLog->SetVisAttributes(pointVisAttr); } @@ -101,16 +93,13 @@ 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(); + double radius = obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentRadius(); + auto center = obj.physical_volume->GetLogicalVolume()->GetSolid()->GetExtent().GetExtentCenter(); - G4LogicalVolume* log = obj.physical_volume->GetLogicalVolume(); + auto log = obj.physical_volume->GetLogicalVolume(); std::vector points; - int i = 0; - while (i < 100) { + for (int i = 0; i < 100; i++) { G4ThreeVector pos; G4ThreeVector dir; obj.GetDirection(dir, pos); @@ -120,9 +109,8 @@ int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { std::cout << "the initial point must be less than the bounding radius" << std::endl; return 1; } - i++; } - G4RunManager* runManager = new G4RunManager(); + auto runManager = new G4RunManager(); // Set mandatory initialization classes runManager->SetUserInitialization(new MyDetectorConstruction(center, radius, points, log)); @@ -132,7 +120,7 @@ int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { runManager->SetUserInitialization(physicsList); runManager->Initialize(); // Initialize visualization - G4VisManager* visManager = new G4VisExecutive(); + auto visManager = new G4VisExecutive(); visManager->Initialize(); // User interface @@ -153,23 +141,17 @@ int RunVis(RMGVertexConfinement::SampleableObject obj, std::string name) { 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); +G4VPhysicalVolume* get_phy_volume(G4VSolid* solid, std::string material_string) { + auto material = G4NistManager::Instance()->FindOrBuildMaterial(material_string); + auto logicalVolume = new G4LogicalVolume(solid, material, "log_vol"); + auto physicalVolume = new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 0), logicalVolume, "phy_vol", + nullptr, false, 0); return physicalVolume; } @@ -182,39 +164,35 @@ int main(int argc, char* argv[]) { // 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 = 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); + auto 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 small_tubby = new G4Tubs("small_tubby", 0 * mm, 10 * mm, 50 * mm, 0, 360 * deg); + auto 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); + auto 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); + auto box = new G4Box("box", 25 * mm, 25 * mm, 50 * mm); + auto 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 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); + auto 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"]; + auto obj = sampleables["tubs"]; // make a basic geometry - first just a G4Tubs // shoot along x std::vector ints = @@ -224,6 +202,7 @@ int main(int argc, char* argv[]) { 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)); @@ -235,36 +214,32 @@ int main(int argc, char* argv[]) { // miss the object ints = obj.GetIntersections(G4ThreeVector(60 * mm, 0, 0), G4ThreeVector(0, -1, 0)); - if (ints.size() != 0) { + if (!ints.empty()) { std::cout << "The number of intersections should be 0" << std::endl; return 1; } // now generate a bunch randomly - int i = 0; - while (i < 10000) { + for (int i = 0; i < 10000; i++) { G4ThreeVector dir; G4ThreeVector pos; obj.GetDirection(dir, pos); - int ints = obj.GetIntersections(pos, dir).size(); + size_t 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 + // shoot along x std::vector ints = obj.GetIntersections(G4ThreeVector(60 * mm, 0, 0), G4ThreeVector(-1, 0, 0)); @@ -274,7 +249,6 @@ int main(int argc, char* argv[]) { } // shoot along z - ints = obj.GetIntersections(G4ThreeVector(15 * mm, 0, 60 * mm), G4ThreeVector(0, 0, -1)); if (ints.size() != 2) { @@ -282,31 +256,25 @@ int main(int argc, char* argv[]) { return 1; } - int i = 0; - while (i < 10000) { + for (int i = 0; i < 10000; i++) { G4ThreeVector dir; G4ThreeVector pos; obj.GetDirection(dir, pos); - int num_ints = obj.GetIntersections(pos, dir).size(); + size_t 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)); @@ -323,28 +291,24 @@ int main(int argc, char* argv[]) { return 1; } - int i = 0; - while (i < 10000) { + for (int i = 0; i < 10000; i++) { G4ThreeVector dir; G4ThreeVector pos; obj.GetDirection(dir, pos); - int num_ints = obj.GetIntersections(pos, dir).size(); + size_t 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) { + for (int i = 0; i < 10000; i++) { G4ThreeVector pos; bool success = obj.GenerateSurfacePoint(pos, 200, 4); @@ -358,13 +322,12 @@ int main(int argc, char* argv[]) { << std::endl; return 1; } - i++; } + } else if (test_type == "test-containment-subtraction") { auto obj = sampleables["sub"]; - int i = 0; - while (i < 10000) { + for (int i = 0; i < 10000; i++) { G4ThreeVector pos; bool success = obj.GenerateSurfacePoint(pos, 200, 4); @@ -377,8 +340,8 @@ int main(int argc, char* argv[]) { 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 @@ -396,11 +359,6 @@ int main(int argc, char* argv[]) { // get some points to plot auto obj = sampleables["tubs"]; return RunVis(obj, "simple"); - - } - - else { - return 0; } return 0;