Geant4-Celeritas offloading template

This is a template for Geant4 applications with Celeritas physics offloading capabilities. It shows how to link Celeritas against Geant4 in the CMakeLists.txt and the Geant4 classes needed to initialize Celeritas, offload events, and recover step information.

Dependencies

• Geant4 v11 or newer

• Celeritas v0.5 or newer with CELERITAS_USE_Geant4=ON

Build and run

$ mkdir build
$ cd build
$ cmake ..
$ make
$ export CELER_DISABLE_PARALLEL=1 # if Celeritas is built with MPI
$ ./run-offload

Example classes

MakeCelerOptions

Build Celeritas integration options before the beginning of the run.

ActionInitialization

Build and BuildForMaster construct the Celeritas integration interface with user-defined options.

RunAction

BeginOfRunAction initializes Celeritas global shared data on master and worker threads, setting up a tracking manager under the hood. EndOfRunAction clears data and finalizes Celeritas data.

SensitiveDetector

ProcessHits: is currently the only Celeritas callback interface to Geant4; at each step, Celeritas sends data back as a G4Step to be processed by Geant4.

StepDiagnostic

This advanced and experimental class (expect it to break at major version changes) demonstrates an efficient GPU-based stepping action that integrates with the main Geant4 code.

Geant4 integration examples

These small examples demonstrate how to offload tracks to Celeritas in a serial or multithreaded environment using:

1. A concrete G4VTrackingManager

2. A concrete G4UserTrackingAction user action class

3. A concrete G4VFastSimulationModel

The Geant4 interface is the only relevant part of Celeritas here. The key components are global celeritas::SetupOptions and celeritas::SharedParams, integrated with celeritas::TrackingManagerIntegration.

CMake infrastructure

project(CeleritasAccelExample VERSION 0.0.1 LANGUAGES CXX)
cmake_policy(VERSION 3.12...3.30)

find_package(Celeritas 0.6 REQUIRED)
find_package(Geant4 REQUIRED)

if(NOT CELERITAS_USE_Geant4)
  message(SEND_ERROR "This example requires Geant4 support "
    "to be enabled in Celeritas")
endif()

add_executable(simple-offload simple-offload.cc)
celeritas_target_link_libraries(simple-offload
  Celeritas::accel
  ${Geant4_LIBRARIES}
)

if(Geant4_VERSION VERSION_LESS 11.1)
  message(WARNING "Fast simulation offload requires Geant4 11.1 or higher")
endif()

add_executable(fastsim-offload fastsim-offload.cc)
celeritas_target_link_libraries(fastsim-offload
  Celeritas::accel
  ${Geant4_LIBRARIES}
)

if(Geant4_VERSION VERSION_LESS 11.0)
  message(WARNING "G4VTrackingManager offload requires Geant4 11.0 or higher")
else()
  add_executable(trackingmanager-offload trackingmanager-offload.cc)
  celeritas_target_link_libraries(trackingmanager-offload
    Celeritas::accel
    ${Geant4_LIBRARIES}
  )
endif()

Main Executables

The executables are less robust (and minimally documented) versions of the Integrated Geant4 application (celer-g4) app. The use of global variables rather than shared pointers is easier to implement but may be more problematic with experiment frameworks or other apps that use a task-based runner.

Offload using a concrete G4VTrackingManager

The tracking manager is the preferred way of offloading tracks to Celeritas. It requires Geant4 11.0 or higher.

#include <algorithm>
#include <iterator>
#include <type_traits>

// Geant4
#include <FTFP_BERT.hh>
#include <G4Box.hh>
#include <G4LogicalVolume.hh>
#include <G4Material.hh>
#include <G4PVPlacement.hh>
#include <G4ParticleGun.hh>
#include <G4ParticleTable.hh>
#include <G4RunManagerFactory.hh>
#include <G4SDManager.hh>
#include <G4SystemOfUnits.hh>
#include <G4ThreeVector.hh>
#include <G4UserEventAction.hh>
#include <G4UserRunAction.hh>
#include <G4UserTrackingAction.hh>
#include <G4VUserActionInitialization.hh>
#include <G4VUserDetectorConstruction.hh>
#include <G4VUserPrimaryGeneratorAction.hh>
#include <G4Version.hh>

// Celeritas
#include <accel/AlongStepFactory.hh>
#include <accel/LocalTransporter.hh>
#include <accel/SetupOptions.hh>
#include <accel/TrackingManagerConstructor.hh>
#include <accel/TrackingManagerIntegration.hh>
#include <corecel/Assert.hh>
#include <corecel/Macros.hh>
#include <corecel/io/Logger.hh>

namespace
{
//---------------------------------------------------------------------------//
class SensitiveDetector final : public G4VSensitiveDetector
{
  public:
    explicit SensitiveDetector(std::string name)
        : G4VSensitiveDetector{std::move(name)}
    {
    }

    double edep() const { return edep_; }

  protected:
    void Initialize(G4HCofThisEvent*) final { edep_ = 0; }
    bool ProcessHits(G4Step* step, G4TouchableHistory*) final
    {
        CELER_ASSERT(step);
        edep_ += step->GetTotalEnergyDeposit();
        return true;
    }

  private:
    double edep_{0};
};

// Simple (not best practice) way of accessing SD
G4ThreadLocal SensitiveDetector const* global_sd{nullptr};

//---------------------------------------------------------------------------//
class DetectorConstruction final : public G4VUserDetectorConstruction
{
  public:
    DetectorConstruction()
        : aluminum_{new G4Material{
              "Aluminium", 13., 26.98 * g / mole, 2.700 * g / cm3}}
    {
    }

    G4VPhysicalVolume* Construct() final
    {
        CELER_LOG_LOCAL(status) << "Setting up detector";
        auto* box = new G4Box("world", 100 * cm, 100 * cm, 100 * cm);
        auto* lv = new G4LogicalVolume(box, aluminum_, "world");
        world_lv_ = lv;
        auto* pv = new G4PVPlacement(
            0, G4ThreeVector{}, lv, "world", nullptr, false, 0);
        return pv;
    }

    void ConstructSDandField() final
    {
        auto* sd_manager = G4SDManager::GetSDMpointer();
        auto detector = std::make_unique<SensitiveDetector>("example-sd");
        world_lv_->SetSensitiveDetector(detector.get());
        global_sd = detector.get();
        sd_manager->AddNewDetector(detector.release());
    }

  private:
    G4Material* aluminum_{nullptr};
    G4LogicalVolume* world_lv_{nullptr};
};

//---------------------------------------------------------------------------//
// Generate 200 MeV pi+
class PrimaryGeneratorAction final : public G4VUserPrimaryGeneratorAction
{
  public:
    PrimaryGeneratorAction()
    {
        auto* g4particle_def
            = G4ParticleTable::GetParticleTable()->FindParticle(211);
        gun_.SetParticleDefinition(g4particle_def);
        gun_.SetParticleEnergy(200 * MeV);
        gun_.SetParticlePosition(G4ThreeVector{0, 0, 0});  // origin
        gun_.SetParticleMomentumDirection(G4ThreeVector{1, 0, 0});  // +x
    }

    void GeneratePrimaries(G4Event* event) final
    {
        CELER_LOG_LOCAL(status) << "Generating primaries";
        gun_.GeneratePrimaryVertex(event);
    }

  private:
    G4ParticleGun gun_;
};

//---------------------------------------------------------------------------//
class RunAction final : public G4UserRunAction
{
  public:
    void BeginOfRunAction(G4Run const* run) final
    {
        celeritas::TrackingManagerIntegration::Instance().BeginOfRunAction(run);
    }
    void EndOfRunAction(G4Run const* run) final
    {
        celeritas::TrackingManagerIntegration::Instance().EndOfRunAction(run);
    }
};

//---------------------------------------------------------------------------//
class EventAction final : public G4UserEventAction
{
  public:
    void BeginOfEventAction(G4Event const*) final {}
    void EndOfEventAction(G4Event const* event) final
    {
        // Log total energy deposition
        if (global_sd)
        {
            CELER_LOG_LOCAL(info)
                << "Total energy deposited for event " << event->GetEventID()
                << ": " << (global_sd->edep() / CLHEP::MeV) << " MeV";
        }
        else
        {
            CELER_LOG_LOCAL(error) << "Global SD was not set";
        }
    }
};

//---------------------------------------------------------------------------//
class ActionInitialization final : public G4VUserActionInitialization
{
  public:
    void BuildForMaster() const final
    {
        celeritas::TrackingManagerIntegration::Instance().BuildForMaster();

        CELER_LOG_LOCAL(status) << "Constructing user actions";

        this->SetUserAction(new RunAction{});
    }
    void Build() const final
    {
        celeritas::TrackingManagerIntegration::Instance().Build();

        CELER_LOG_LOCAL(status) << "Constructing user actions";

        this->SetUserAction(new PrimaryGeneratorAction{});
        this->SetUserAction(new RunAction{});
        this->SetUserAction(new EventAction{});
    }
};

//---------------------------------------------------------------------------//
/*!
 * Construct options for Celeritas.
 */
celeritas::SetupOptions MakeOptions()
{
    celeritas::SetupOptions opts;
    // NOTE: these numbers are appropriate for CPU execution
    opts.max_num_tracks = 2024;
    opts.initializer_capacity = 2024 * 128;
    // Celeritas does not support EmStandard MSC physics above 200 MeV
    opts.ignore_processes = {"CoulombScat"};

    // Use a uniform (zero) magnetic field
    opts.make_along_step = celeritas::UniformAlongStepFactory();

    // Export a GDML file with the problem setup and SDs
    opts.geometry_output_file = "simple-example.gdml";
    // Save diagnostic file to a unique name
    opts.output_file = "trackingmanager-offload.out.json";
    return opts;
}

//---------------------------------------------------------------------------//
}  // namespace

int main()
{
    auto run_manager = std::unique_ptr<G4RunManager>{
        G4RunManagerFactory::CreateRunManager()};

    run_manager->SetUserInitialization(new DetectorConstruction{});

    auto& tmi = celeritas::TrackingManagerIntegration::Instance();

    // Use FTFP_BERT, but use Celeritas tracking for e-/e+/g
    auto* physics_list = new FTFP_BERT{/* verbosity = */ 0};
    physics_list->RegisterPhysics(
        new celeritas::TrackingManagerConstructor(&tmi));
    run_manager->SetUserInitialization(physics_list);
    run_manager->SetUserInitialization(new ActionInitialization());

    tmi.SetOptions(MakeOptions());

    run_manager->Initialize();
    run_manager->BeamOn(2);

    return 0;
}

Offload using a concrete G4UserTrackingAction

#include <algorithm>
#include <iterator>
#include <type_traits>
#include <FTFP_BERT.hh>
#include <G4Box.hh>
#include <G4Electron.hh>
#include <G4Gamma.hh>
#include <G4LogicalVolume.hh>
#include <G4Material.hh>
#include <G4PVPlacement.hh>
#include <G4ParticleDefinition.hh>
#include <G4ParticleGun.hh>
#include <G4ParticleTable.hh>
#include <G4Positron.hh>
#include <G4SystemOfUnits.hh>
#include <G4Threading.hh>
#include <G4ThreeVector.hh>
#include <G4Track.hh>
#include <G4TrackStatus.hh>
#include <G4Types.hh>
#include <G4UserEventAction.hh>
#include <G4UserRunAction.hh>
#include <G4UserTrackingAction.hh>
#include <G4VUserActionInitialization.hh>
#include <G4VUserDetectorConstruction.hh>
#include <G4VUserPrimaryGeneratorAction.hh>
#include <G4Version.hh>
#if G4VERSION_NUMBER >= 1200
#    include <G4RunManagerFactory.hh>
#else
#    include <G4MTRunManager.hh>
#endif

#include <accel/AlongStepFactory.hh>
#include <accel/LocalTransporter.hh>
#include <accel/SetupOptions.hh>
#include <accel/SharedParams.hh>
#include <accel/UserActionIntegration.hh>
#include <corecel/Macros.hh>
#include <corecel/io/Logger.hh>

using celeritas::UserActionIntegration;

namespace
{
//---------------------------------------------------------------------------//
class DetectorConstruction final : public G4VUserDetectorConstruction
{
  public:
    DetectorConstruction()
        : aluminum_{new G4Material{
              "Aluminium", 13., 26.98 * g / mole, 2.700 * g / cm3}}
    {
    }

    G4VPhysicalVolume* Construct() final
    {
        CELER_LOG_LOCAL(status) << "Setting up geometry";
        auto* box = new G4Box("world", 100 * cm, 100 * cm, 100 * cm);
        auto* lv = new G4LogicalVolume(box, aluminum_, "world");
        auto* pv = new G4PVPlacement(
            0, G4ThreeVector{}, lv, "world", nullptr, false, 0);
        return pv;
    }

    void ConstructSDandField() final {}

  private:
    G4Material* aluminum_;
};

//---------------------------------------------------------------------------//
// Generate 200 MeV pi+
class PrimaryGeneratorAction final : public G4VUserPrimaryGeneratorAction
{
  public:
    PrimaryGeneratorAction()
    {
        auto g4particle_def
            = G4ParticleTable::GetParticleTable()->FindParticle(211);
        gun_.SetParticleDefinition(g4particle_def);
        gun_.SetParticleEnergy(200 * MeV);
        gun_.SetParticlePosition(G4ThreeVector{0, 0, 0});  // origin
        gun_.SetParticleMomentumDirection(G4ThreeVector{1, 0, 0});  // +x
    }

    void GeneratePrimaries(G4Event* event) final
    {
        CELER_LOG_LOCAL(status) << "Generating primaries";
        gun_.GeneratePrimaryVertex(event);
    }

  private:
    G4ParticleGun gun_;
};

//---------------------------------------------------------------------------//
class RunAction final : public G4UserRunAction
{
  public:
    void BeginOfRunAction(G4Run const* run) final
    {
        UserActionIntegration::Instance().BeginOfRunAction(run);
    }
    void EndOfRunAction(G4Run const* run) final
    {
        UserActionIntegration::Instance().EndOfRunAction(run);
    }
};

//---------------------------------------------------------------------------//
class EventAction final : public G4UserEventAction
{
  public:
    void BeginOfEventAction(G4Event const* event) final
    {
        UserActionIntegration::Instance().BeginOfEventAction(event);
    }

    void EndOfEventAction(G4Event const* event) final
    {
        UserActionIntegration::Instance().EndOfEventAction(event);
    }
};

//---------------------------------------------------------------------------//
class TrackingAction final : public G4UserTrackingAction
{
    void PreUserTrackingAction(G4Track const* track) final
    {
        UserActionIntegration::Instance().PreUserTrackingAction(
            const_cast<G4Track*>(track));
    }
};

//---------------------------------------------------------------------------//
class ActionInitialization final : public G4VUserActionInitialization
{
  public:
    void BuildForMaster() const final
    {
        UserActionIntegration::Instance().BuildForMaster();

        CELER_LOG_LOCAL(status) << "Constructing user actions";

        this->SetUserAction(new RunAction{});
    }
    void Build() const final
    {
        UserActionIntegration::Instance().Build();

        CELER_LOG_LOCAL(status) << "Constructing user actions";

        this->SetUserAction(new PrimaryGeneratorAction{});
        this->SetUserAction(new RunAction{});
        this->SetUserAction(new EventAction{});
        this->SetUserAction(new TrackingAction{});
    }
};

//---------------------------------------------------------------------------//
/*!
 * Construct options for Celeritas.
 */
celeritas::SetupOptions MakeOptions()
{
    celeritas::SetupOptions opts;

    opts.make_along_step = celeritas::UniformAlongStepFactory();

    // NOTE: since no SD is enabled, we must manually disable Celeritas hit
    // processing
    opts.sd.enabled = false;

    // NOTE: these numbers are appropriate for CPU execution
    opts.max_num_tracks = 2024;
    opts.initializer_capacity = 2024 * 128;
    // Celeritas does not support EmStandard MSC physics above 200 MeV
    opts.ignore_processes = {"CoulombScat"};

    opts.output_file = "simple-offload.out.json";
    return opts;
}

//---------------------------------------------------------------------------//
}  // namespace

int main()
{
    auto run_manager = [] {
#if G4VERSION_NUMBER >= 1200
        return std::unique_ptr<G4RunManager>{
            G4RunManagerFactory::CreateRunManager()};
#else
        return std::make_unique<G4RunManager>();
#endif
    }();

    run_manager->SetUserInitialization(new DetectorConstruction{});
    run_manager->SetUserInitialization(new FTFP_BERT{/* verbosity = */ 0});
    run_manager->SetUserInitialization(new ActionInitialization());

    UserActionIntegration::Instance().SetOptions(MakeOptions());

    run_manager->Initialize();
    run_manager->BeamOn(2);

    return 0;
}

Offload using a concrete G4VFastSimulationModel

#include <algorithm>
#include <iterator>
#include <type_traits>
#include <FTFP_BERT.hh>
#include <G4Box.hh>
#include <G4Electron.hh>
#include <G4FastSimulationPhysics.hh>
#include <G4Gamma.hh>
#include <G4LogicalVolume.hh>
#include <G4Material.hh>
#include <G4PVPlacement.hh>
#include <G4ParticleDefinition.hh>
#include <G4ParticleGun.hh>
#include <G4ParticleTable.hh>
#include <G4Positron.hh>
#include <G4Region.hh>
#include <G4RegionStore.hh>
#include <G4SystemOfUnits.hh>
#include <G4Threading.hh>
#include <G4ThreeVector.hh>
#include <G4Track.hh>
#include <G4TrackStatus.hh>
#include <G4Types.hh>
#include <G4UserEventAction.hh>
#include <G4UserRunAction.hh>
#include <G4UserTrackingAction.hh>
#include <G4VUserActionInitialization.hh>
#include <G4VUserDetectorConstruction.hh>
#include <G4VUserPrimaryGeneratorAction.hh>
#include <G4Version.hh>
#if G4VERSION_NUMBER >= 1200
#    include <G4RunManagerFactory.hh>
#else
#    include <G4MTRunManager.hh>
#endif

#include <accel/AlongStepFactory.hh>
#include <accel/FastSimulationIntegration.hh>
#include <accel/FastSimulationModel.hh>
#include <accel/SetupOptions.hh>
#include <corecel/Macros.hh>
#include <corecel/io/Logger.hh>

using celeritas::FastSimulationIntegration;

namespace
{
//---------------------------------------------------------------------------//
class DetectorConstruction final : public G4VUserDetectorConstruction
{
  public:
    DetectorConstruction()
        : aluminum_{new G4Material{
              "Aluminium", 13., 26.98 * g / mole, 2.700 * g / cm3}}
    {
    }

    G4VPhysicalVolume* Construct() final
    {
        CELER_LOG_LOCAL(status) << "Setting up geometry";
        auto* box = new G4Box("world", 100 * cm, 100 * cm, 100 * cm);
        auto* lv = new G4LogicalVolume(box, aluminum_, "world");
        auto* pv = new G4PVPlacement(
            0, G4ThreeVector{}, lv, "world", nullptr, false, 0);
        return pv;
    }

    void ConstructSDandField() final
    {
        CELER_LOG_LOCAL(status)
            << R"(Creating FastSimulationModel for default region)";
        G4Region* default_region = G4RegionStore::GetInstance()->GetRegion(
            "DefaultRegionForTheWorld");
        // Underlying GVFastSimulationModel constructor handles ownership, so
        // we must ignore the returned pointer...
        new celeritas::FastSimulationModel(default_region);
    }

  private:
    G4Material* aluminum_;
};

//---------------------------------------------------------------------------//
// Generate 200 MeV pi+
class PrimaryGeneratorAction final : public G4VUserPrimaryGeneratorAction
{
  public:
    PrimaryGeneratorAction()
    {
        auto g4particle_def
            = G4ParticleTable::GetParticleTable()->FindParticle(211);
        gun_.SetParticleDefinition(g4particle_def);
        gun_.SetParticleEnergy(200 * MeV);
        gun_.SetParticlePosition(G4ThreeVector{0, 0, 0});  // origin
        gun_.SetParticleMomentumDirection(G4ThreeVector{1, 0, 0});  // +x
    }

    void GeneratePrimaries(G4Event* event) final
    {
        CELER_LOG_LOCAL(status) << "Generating primaries";
        gun_.GeneratePrimaryVertex(event);
    }

  private:
    G4ParticleGun gun_;
};

//---------------------------------------------------------------------------//
class RunAction final : public G4UserRunAction
{
  public:
    void BeginOfRunAction(G4Run const* run) final
    {
        FastSimulationIntegration::Instance().BeginOfRunAction(run);
    }
    void EndOfRunAction(G4Run const* run) final
    {
        FastSimulationIntegration::Instance().EndOfRunAction(run);
    }
};

//---------------------------------------------------------------------------//
class ActionInitialization final : public G4VUserActionInitialization
{
  public:
    void BuildForMaster() const final
    {
        FastSimulationIntegration::Instance().BuildForMaster();

        CELER_LOG_LOCAL(status) << "Constructing user actions";

        this->SetUserAction(new RunAction{});
    }
    void Build() const final
    {
        FastSimulationIntegration::Instance().Build();

        CELER_LOG_LOCAL(status) << "Constructing user actions";

        this->SetUserAction(new PrimaryGeneratorAction{});
        this->SetUserAction(new RunAction{});
    }
};

//---------------------------------------------------------------------------//
/*!
 * Construct options for Celeritas.
 */
celeritas::SetupOptions MakeOptions()
{
    celeritas::SetupOptions opts;
    // NOTE: these numbers are appropriate for CPU execution
    opts.max_num_tracks = 2024;
    opts.initializer_capacity = 2024 * 128;
    // Celeritas does not support EmStandard MSC physics above 200 MeV
    opts.ignore_processes = {"CoulombScat"};

    // NOTE: since no SD is enabled, we must manually disable Celeritas hit
    // processing
    opts.sd.enabled = false;

    // Use a uniform (zero) magnetic field
    opts.make_along_step = celeritas::UniformAlongStepFactory();

    // Export a GDML file with the problem setup and SDs
    opts.geometry_output_file = "simple-example.gdml";
    // Save diagnostic file to a unique name
    opts.output_file = "fastsim-offload.out.json";
    return opts;
}

//---------------------------------------------------------------------------//
}  // namespace

int main()
{
    auto run_manager = [] {
#if G4VERSION_NUMBER >= 1200
        return std::unique_ptr<G4RunManager>{
            G4RunManagerFactory::CreateRunManager()};
#else
        return std::make_unique<G4RunManager>();
#endif
    }();

    run_manager->SetUserInitialization(new DetectorConstruction{});

    // We must add support for fast simulation models to the Physics List
    // NOTE: we have to explicitly name the particles and this should be a
    // superset of what Celeritas can offload
    auto physics_list = new FTFP_BERT{/* verbosity = */ 0};
    auto fast_physics = new G4FastSimulationPhysics();
    fast_physics->ActivateFastSimulation("e-");
    fast_physics->ActivateFastSimulation("e+");
    fast_physics->ActivateFastSimulation("gamma");
    physics_list->RegisterPhysics(fast_physics);
    run_manager->SetUserInitialization(physics_list);

    run_manager->SetUserInitialization(new ActionInitialization());

    FastSimulationIntegration::Instance().SetOptions(MakeOptions());

    run_manager->Initialize();
    run_manager->BeamOn(2);

    return 0;
}