diff --git a/nain4/src/meson.build b/nain4/src/meson.build index a269507c..0637b71d 100644 --- a/nain4/src/meson.build +++ b/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' diff --git a/nain4/src/n4-will-become-external-lib.cc b/nain4/src/n4-will-become-external-lib.cc new file mode 100644 index 00000000..093e9ef5 --- /dev/null +++ b/nain4/src/n4-will-become-external-lib.cc @@ -0,0 +1,90 @@ + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +std::vector 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 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; + } +} diff --git a/nain4/src/n4-will-become-external-lib.hh b/nain4/src/n4-will-become-external-lib.hh new file mode 100644 index 00000000..3dfbf1b3 --- /dev/null +++ b/nain4/src/n4-will-become-external-lib.hh @@ -0,0 +1,16 @@ +#pragma once + +#include +#include +#include + +struct test_config { + G4VUserPhysicsList* physics; + G4Material* material; + std::string particle_name; + double particle_energy; + std::vector distances; +}; + + +std::vector measure_abslength(test_config const&); diff --git a/templates/basic/test/test-LXe.cc b/templates/basic/test/test-LXe.cc index ea31696e..f08be706 100644 --- a/templates/basic/test/test-LXe.cc +++ b/templates/basic/test/test-LXe.cc @@ -1,6 +1,7 @@ #include "LXe.hh" #include +#include #include #include @@ -28,44 +29,11 @@ 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}); @@ -73,43 +41,15 @@ TEST_CASE("liquid xenon properties", "[.xfail][xenon][properties]") { 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)); } }