Add measure_abslength generic test utility

nain4/src/meson.build

@@ -46,6 +46,7 @@ nain4_includes = [ 'n4-all.hh'
                  , 'n4-exceptions.hh'
                  , 'n4-geometry-iterators.hh'
                  , 'n4-geometry.hh'
+                 , 'n4-will-become-external-lib.hh'
                  , 'n4-inspect.hh'
                  , 'n4-main.hh'
                  , 'n4-mandatory.hh'
@@ -68,6 +69,7 @@ nain4_includes = [ 'n4-all.hh'
 nain4_sources = [ 'n4-boolean-shape.cc'
                 , 'n4-constants.cc'
                 , 'n4-geometry-iterators.cc'
+                , 'n4-will-become-external-lib.cc'
                 , 'n4-mandatory.cc'
                 , 'n4-material.cc'
                 , 'n4-place.cc'

nain4/src/n4-will-become-external-lib.cc

@@ -0,0 +1,90 @@
+
+#include <n4-will-become-external-lib.hh>
+#include <n4-inspect.hh>
+#include <n4-material.hh>
+#include <n4-sequences.hh>
+#include <n4-shape.hh>
+#include <n4-stream.hh>
+
+#include <G4PrimaryVertex.hh>
+#include <G4RandomDirection.hh>
+#include <G4VUserPhysicsList.hh>
+
+
+std::vector<double> measure_abslength(test_config const& config) {
+  auto air    = n4::material("G4_AIR");
+  auto sphere = [material=config.material, air] (auto radius) {
+    return [material, air, radius] () {
+      auto lab = n4::box   ("LAB"   ).half_cube(1.1*radius).place(air     )        .now();
+                 n4::sphere("Sphere").r        (    radius).place(material).in(lab).now();
+      return lab;
+    };
+  };
+
+  // --- Generator -----
+  auto shoot_gamma = [&config](G4Event* event) {
+     auto vertex = new G4PrimaryVertex{};
+     auto p      = config.particle_energy * G4RandomDirection();
+     auto particle = n4::find_particle(config.particle_name);
+     vertex -> SetPrimary(new G4PrimaryParticle(particle, p.x(), p.y(), p.z()));
+     event -> AddPrimaryVertex(vertex);
+  };
+
+  // --- Count unscathed gammas (in stepping action) -----
+  size_t unscathed = 0;
+  auto count_unscathed = [&unscathed, initial_energy=config.particle_energy](auto step) {
+    auto energy = step -> GetTrack() -> GetTotalEnergy();
+    if (energy > 0.999999 * initial_energy) {
+      auto name = step -> GetPreStepPoint() -> GetTouchable() -> GetVolume() -> GetName();
+      if (name == "LAB") { unscathed++; }
+    }
+  };
+
+  // --- Eliminate secondaries (in stacking action)  -----
+  auto kill_secondaries = [](auto track) {
+    auto kill = track -> GetParentID() > 0;
+    return kill > 0 ? G4ClassificationOfNewTrack::fKill : G4ClassificationOfNewTrack::fUrgent;
+  };
+
+  auto create_actions = [&] {
+    return (new n4::actions{shoot_gamma})
+         -> set((new n4::stacking_action) -> classify(kill_secondaries))
+         -> set (new n4::stepping_action{count_unscathed});
+  };
+
+  // ----- Initialize and run Geant4 -------------------------------------------
+  {
+    n4::silence _{G4cout};
+
+    auto rm = n4::run_manager::create()
+      .fake_ui ()
+      .physics (config.physics)
+      .geometry(sphere(1))
+      .actions(create_actions)
+      .run()
+      ;
+
+    auto events = 10000;
+    auto radii  = config.distances;
+    std::vector<double> result;
+    result.reserve(radii.size());
+
+    // --- Infer attenuation length by gathering statistics for given radius -------------
+    auto estimate_att_length = [&unscathed, rm, &sphere, events, &result] (auto radius) {
+      unscathed = 0;
+
+      rm -> replace_geometry(sphere(radius)).run(events);
+
+      auto ratio = unscathed / (1.0 * events);
+      auto expected_attenuation_length = 3.74 * cm;
+      auto attenuation_length = - radius / log(ratio);
+      result.push_back(attenuation_length);
+   };
+
+    // --- Check attenuation length across range of radii --------------------------------
+    for (auto radius : radii) {
+      estimate_att_length(radius);
+    }
+    return result;
+  }
+}

nain4/src/n4-will-become-external-lib.hh

@@ -0,0 +1,16 @@
+#pragma once
+
+#include <G4Material.hh>
+#include <G4ParticleDefinition.hh>
+#include <G4VUserPhysicsList.hh>
+
+struct test_config {
+  G4VUserPhysicsList* physics;
+  G4Material*         material;
+  std::string         particle_name;
+  double              particle_energy;
+  std::vector<double> distances;
+};
+
+
+std::vector<double> measure_abslength(test_config const&);

templates/basic/test/test-LXe.cc

@@ -1,6 +1,7 @@
 #include "LXe.hh"
 
 #include <n4-all.hh>
+#include <n4-will-become-external-lib.hh>
 
 #include <G4OpticalPhysics.hh>
 #include <G4PVPlacement.hh>
@@ -28,88 +29,27 @@ TEST_CASE("dummy test2", "[dummy]") {
 
 TEST_CASE("liquid xenon properties", "[.xfail][xenon][properties]") {
   // --- Geometry -----
-  auto air = n4::material("G4_AIR");
   auto LXe = LXe_with_properties();
   REQUIRE_THAT( LXe -> GetDensity() / (g / cm3)
               , WithinRel(2.98, 1e-6)
               );
 
-  auto xe_sphere = [&air, &LXe] (auto radius) {
-    return [&air, &LXe, radius] () {
-      auto lab = n4::box   ("LAB"   ).half_cube(1.1*radius).place(air)        .now();
-                 n4::sphere("Sphere").r        (    radius).place(LXe).in(lab).now();
-      return lab;
-    };
-  };
-
-  // --- Generator -----
-  auto shoot_gamma = [](G4Event* event) {
-     auto vertex = new G4PrimaryVertex{};
-     auto p      = 511*keV * G4RandomDirection();
-     vertex -> SetPrimary(new G4PrimaryParticle(n4::find_particle("gamma"), p.x(), p.y(), p.z()));
-     event -> AddPrimaryVertex(vertex);
-  };
-
-  // --- Count unscathed gammas (in stepping action) -----
-  size_t unscathed = 0;
-  auto count_unscathed = [&unscathed](auto step) {
-    auto energy = step -> GetTrack() -> GetTotalEnergy();
-    if (WithinRel(energy, 1e-6).match(511*keV)) { // ignore post-Compton gammas
-      auto name = step -> GetPreStepPoint() -> GetTouchable() -> GetVolume() -> GetName();
-      if (name == "LAB") { unscathed++; }
-    }
-  };
-
-  // --- Eliminate secondaries (in stacking action)  -----
-  auto kill_secondaries = [](auto track) {
-    auto kill = track -> GetParentID() > 0;
-    return kill > 0 ? G4ClassificationOfNewTrack::fKill : G4ClassificationOfNewTrack::fUrgent;
-  };
-
   auto our_optical_physics = [&] {
     auto phys_list = new FTFP_BERT{0};
     phys_list -> ReplacePhysics(new G4EmStandardPhysics_option4{0});
     phys_list -> ReplacePhysics(new G4OpticalPhysics{0});
     return phys_list;
   };
 
-  auto create_actions = [&] {
-    return (new n4::actions{shoot_gamma})
-         -> set((new n4::stacking_action) -> classify(kill_secondaries))
-         -> set (new n4::stepping_action{count_unscathed});
-  };
-
-  // ----- Initialize and run Geant4 -------------------------------------------
-  {
-    n4::silence _{G4cout};
-
-    auto rm = n4::run_manager::create()
-      .fake_ui ()
-      .physics (our_optical_physics)
-      .geometry(xe_sphere(1))
-      .actions(create_actions)
-      .run()
-      ;
-
-    auto events = 10000;
-    auto radii  = n4::scale_by(cm, {1, 2, 3, 4, 5, 6, 7, 8});
-
-    // --- Infer attenuation length by gathering statistics for given radius -------------
-    auto check_att_length = [&unscathed, rm, &xe_sphere, events] (auto radius) {
-      unscathed = 0;
-
-      rm -> replace_geometry(xe_sphere(radius)).run(events);
 
-      auto ratio = unscathed / (1.0 * events);
-      auto expected_attenuation_length = 3.74 * cm;
-      auto attenuation_length = - radius / log(ratio);
-      REQUIRE_THAT( attenuation_length
-                  , WithinRel(expected_attenuation_length, 0.05));
-   };
+  auto abs_lengths = measure_abslength(test_config{ .physics         = our_optical_physics()
+                                                  , .material        = LXe
+                                                  , .particle_name   = "gamma"
+                                                  , .particle_energy = 511 * keV}
+                                                  , .distances       = n4::scale_by(cm , {1, 2, 3, 4, 5, 6, 7, 8});
 
-    // --- Check attenuation length across range of radii --------------------------------
-    for (auto radius : radii) {
-      check_att_length(radius);
-    }
+  auto expected_abs_length = 3.74 * cm;
+  for (auto abs_length : abs_lengths) {
+    CHECK_THAT(abs_length, WithinRel(expected_abs_length, 0.05));
   }
 }