Reset run manager/kernel in single process

For the purposes of creating a collection of thorough, robust and conveniently-usable automated tests, I am looking for the ability to reset the run manager completely.

  • IIUC, constructing an instance of G4RunManager more than once in the same process (even if the previous one has been deleted), is doomed to failure (guaranteed to crash).

    • Is my understanding correct?

    • Is there something about the structure and design of Geant4 that makes this a fundamentally insurmountable problem, even if there were the will in the Geant4 core team, to make it possible?

  • Is there some way to reset an existing run manager to a tabula rasa state?

  • Is there some way to completely replace the geometry/physics/primary generator
    in an existing run manager?

  • Can you offer any wisdom on the topic of writing extensive unit and integration tests for Geant4 applications and their components?

Hi @jacek

Any progress on this?

Thanks.

Digging through the Geant4 source code led me to the conclusion that G4RunManager and G4RunManagerKernel are fundamentally hostile to this idea.

I’m using Catch2 which integrates with cmake's ctest. This gives me two, far from ideal but tolerable options:

  1. ctest runs different Catch2 TEST_CASEs in separate processes. So two different test cases can set up their own run manager, and not interfere with each other at all. Therefore simply typing ctest works—it runs all the tests without any danger of crashing because of multiple run manager instantiations—but it:

    • suppresses the more fine-grained information that Catch2 displays in its output,
    • blocks access to Catch2s test selection mechanism

    so …

  2. I base my test execution on this kind of loop

    make -j && ./my-test-executable --list-test-names-only | while read testname; do ./my-t 
    test-executable "$testname"; done
    

    where my-test-executable is the Catch2 test suite executable produced by make.

    In other words, I use a shell loop to launch each test case in a separate process.

    With a justfile to provide a convenient interface to some variations on this theme (such as ones which cater for Catch2s test selection features), this has been enough to relegate this issue from the Top 10 Banes of my Life, quite some time ago and I haven’t given it much thought since then, so I don’t remember details of the hoops I had to jump through to make it work, and I might be forgetting something important.

Once you get over the hurdles of the tests crashing, in order to actually write some (and keep writing more) tests with their own independent run managers, you have to do battle with the mind-numbing verbosity of the G4 mandatory classes, geometry definition, etc. If you’re describing a huge detector, the overhead of this noise might be bearable if you have a high tolerance for this sort of noise. But when trying to write simple, self-contained tests, this verbosity accounts for around 90% of the code, making a huge barrier to writing tests. Without my Geant4 noise-reducing utilities (nain4 :slight_smile:) I doubt I’d be able to muster the energy to write these tests, even though I believe them to be fundamentally important.

But that's a different story ...

Here’s an example of such a test. It’s a full Geant4 app, inside a Catch2 test case, containing definitions of

  • dectector description (xe_sphere)
  • stacking action (kill_secondaries)
  • stepping action (count_unscathed)

defined in around 60 lines of code

#include <nain4.hh>
#include <catch2/catch.hpp>
// ...

TEST_CASE("liquid xenon properties", "[xenon][properties]") {

  auto LXe = LXe_with_properties();
  CHECK(LXe -> GetDensity() / (g / cm3) == Approx(2.98));

  // --- Geometry: sphere of given radius inside an air envelope
  auto xe_sphere = [&LXe] (auto radius) {
    return [&LXe, radius] () {
      auto air = n4::material("G4_AIR");
      auto l2 = 1.1 * radius;
      auto sphere = n4::volume<G4Orb>("Sphere", LXe, radius);
      auto lab    = n4::volume<G4Box>("LAB"   , air, l2, l2, l2);
      n4::place(sphere).in(lab).now();
      return n4::place(lab).now();
    };
  };

  // --- Generator -----
  auto two_gammas_at_origin = [](auto event){generate_back_to_back_511_keV_gammas(event, {}, 0);};

  // --- Count unscathed gammas (in stepping action) -----
  size_t unscathed = 0;
  auto count_unscathed = [&unscathed](auto step) {
    auto energy = step -> GetTrack() -> GetTotalEnergy();
    if (energy == Approx(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;
  };

  // ----- Initialize and run Geant4 -------------------------------------------
  {
    n4::silence _{G4cout};
    auto run_manager = G4RunManager::GetRunManager();
    n4::use_our_optical_physics(run_manager);
    run_manager -> SetUserInitialization((new n4::actions{two_gammas_at_origin})
         -> set((new n4::stacking_action) -> classify(kill_secondaries))
         -> set (new n4::stepping_action{count_unscathed}));

    // --- Infer attenuation length by gathering statistics for given radius -------------
    auto check_attlength = [&unscathed, run_manager](auto build, auto radius, auto events) {
      unscathed = 0;
      n4::clear_geometry();
      run_manager -> SetUserInitialization(new n4::geometry{build(radius)});
      run_manager -> Initialize();
      run_manager -> BeamOn(events);
      auto gammas_sent = 2.0 * events;
      auto ratio = unscathed / gammas_sent;
      auto expected_attenuation_length = 3.74 * cm;
      auto attenuation_length = - radius / log(ratio);
      CHECK(attenuation_length == Approx(expected_attenuation_length).epsilon(0.05));
    };
    // --- Check attenuation length across range of radii --------------------------------
    size_t events = 10000;
    for (auto r: {1,2,3,4,5,6,7,8}) { check_attlength(xe_sphere, r*cm, events); }
  }
}