diff --git a/docs/Project.toml b/docs/Project.toml index e857dff..f017778 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,6 +5,7 @@ DisplayAs = "0b91fe84-8a4c-11e9-3e1d-67c38462b6d6" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FHist = "68837c9b-b678-4cd5-9925-8a54edc8f695" Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44" +Geant4_jll = "872b6946-528a-5ac7-9145-d37eec569368" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" IGLWrap_jll = "283677c1-8365-580c-84e5-ef4b5d190868" IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" diff --git a/docs/make.jl b/docs/make.jl index cb99e93..eb25f6f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -38,10 +38,10 @@ end basic_mds = process_literate("B1", "B2a", "B2aVis", "B3a") extend_mds = process_literate("GPS", "RE03", "TestEm3", "Solids") -advanced_mds = process_literate("TPCSim", "HBC30", "WaterPhantom", "UserLib") +advanced_mds = process_literate("TPCSim", "HBC30", "WaterPhantom", "UserLib", "JuliaAction") extra_mds = create_extras("B2aDetector.jl", "B2aVisSettings.jl", "B3Detector.jl", "GPSDetector.jl", - "RE03Detector.jl", "TestEm3Detector.jl", "HBC30Detector.jl", "UserLibrary.cpp") - + "RE03Detector.jl", "TestEm3Detector.jl", "HBC30Detector.jl", "UserLibrary.cpp", + "MyCode.jl", "G4example.cpp") examples_mds = [] append!(examples_mds, basic_mds, extend_mds, advanced_mds) diff --git a/docs/src/examples/G4example.cpp b/docs/src/examples/G4example.cpp new file mode 100644 index 0000000..1ec2656 --- /dev/null +++ b/docs/src/examples/G4example.cpp @@ -0,0 +1,170 @@ +#include "G4VUserDetectorConstruction.hh" +#include "G4VUserPrimaryGeneratorAction.hh" +#include "G4UserEventAction.hh" +#include "G4UserRunAction.hh" +#include "G4ParticleGun.hh" +#include "G4RunManager.hh" +#include "G4UImanager.hh" +#include "G4RunManagerFactory.hh" +#include "QBBC.hh" +#include "G4NistManager.hh" +#include "G4Box.hh" +#include "G4LogicalVolume.hh" +#include "G4PVPlacement.hh" +//#include "G4SystemOfUnits.hh" +#include "globals.hh" +#include "G4VUserActionInitialization.hh" +#include +#include + +//---Detector construction class------------------------------------------------------------------- +class DetectorConstruction : public G4VUserDetectorConstruction +{ + public: + DetectorConstruction() = default; + ~DetectorConstruction() override = default; + G4VPhysicalVolume* Construct() override { + auto nist = G4NistManager::Instance(); + auto world_mat = nist->FindOrBuildMaterial("G4_AIR"); + auto core_mat = nist->FindOrBuildMaterial("G4_WATER"); + auto world_size = 1.0*CLHEP::m; + auto solidWorld = new G4Box("World", world_size, world_size, world_size); + auto logicWorld = new G4LogicalVolume(solidWorld, world_mat, "World"); + auto physWorld = new G4PVPlacement(0, G4ThreeVector(), logicWorld, "World", 0, false, 0); + auto core_size = 0.5*CLHEP::m; + auto solidCore = new G4Box("Core", core_size, core_size, core_size); + auto logicCore = new G4LogicalVolume(solidCore, core_mat, "Core"); + new G4PVPlacement(0, G4ThreeVector(), logicCore, "Core", logicWorld, false, 0); + return physWorld; + } +}; + +//---Primary generator action class---------------------------------------------------------------- +class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction +{ +public: + PrimaryGeneratorAction() { fParticleGun = new G4ParticleGun(); } + ~PrimaryGeneratorAction() { delete fParticleGun; } + void GeneratePrimaries(G4Event* anEvent) override { + fPrimaryParticleName = fParticleGun->GetParticleDefinition()->GetParticleName(); + fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); + fParticleGun->SetParticlePosition(G4ThreeVector(0.,0.,-1.*CLHEP::m)); + fParticleGun->GeneratePrimaryVertex(anEvent); + } + G4ParticleGun* GetParticleGun() {return fParticleGun;}; + const G4String& GetParticleName() { return fPrimaryParticleName;} +private: + G4String fPrimaryParticleName; + G4ParticleGun* fParticleGun; +}; + +typedef void (*stepaction_f)(const G4Step*); +class SteppingAction : public G4UserSteppingAction +{ +public: + SteppingAction() { + action_jl = (stepaction_f)(jl_unbox_voidpointer(jl_eval_string("@cfunction(stepping_action, Nothing, (CxxPtr{G4Step},))"))); + std::cout << "=====> " << action_jl << std::endl; + } + ~SteppingAction() {} + void UserSteppingAction(const G4Step* step) override { action_jl(step); } +private: + stepaction_f action_jl; +}; + +typedef void (*eventaction_f)(const G4Event*); +class EventAction : public G4UserEventAction +{ + public: + EventAction() { + beginevent_jl = (eventaction_f)(jl_unbox_voidpointer(jl_eval_string("@cfunction(begin_of_event_action, Nothing, (CxxPtr{G4Event},))"))); + endevent_jl = (eventaction_f)(jl_unbox_voidpointer(jl_eval_string("@cfunction(end_of_event_action, Nothing, (CxxPtr{G4Event},))"))); + } + ~EventAction() override = default; + + void BeginOfEventAction(const G4Event* event) override { beginevent_jl(event); } + void EndOfEventAction(const G4Event* event) override { endevent_jl(event); } + private: + eventaction_f beginevent_jl; + eventaction_f endevent_jl; +}; + +typedef void (*runaction_f)(const G4Run*); +class RunAction : public G4UserRunAction +{ + public: + RunAction() { + beginrun_jl = (runaction_f)(jl_unbox_voidpointer(jl_eval_string("@cfunction(begin_of_run_action, Nothing, (CxxPtr{G4Run},))"))); + endrun_jl = (runaction_f)(jl_unbox_voidpointer(jl_eval_string("@cfunction(end_of_run_action, Nothing, (CxxPtr{G4Run},))"))); + } + ~RunAction() override = default; + + void BeginOfRunAction(const G4Run* run) override { beginrun_jl(run); } + void EndOfRunAction(const G4Run* run) override { endrun_jl(run); } + + private: + runaction_f beginrun_jl; + runaction_f endrun_jl; +}; + +//---Action initialization class------------------------------------------------------------------- +class ActionInitialization : public G4VUserActionInitialization +{ + public: + ActionInitialization() = default; + ~ActionInitialization() override = default; + void BuildForMaster() const override {} + void Build() const override { + SetUserAction(new PrimaryGeneratorAction); + SetUserAction(new RunAction); + SetUserAction(new EventAction); + SetUserAction(new SteppingAction); + } +}; + +JULIA_DEFINE_FAST_TLS // only define this once, in an executable (not in a shared library) if you want fast code. + +//----Main program--------------------------------------------------------------------------------- +int main(int, char**) +{ + //--- Required to setup the Julia context + jl_init(); + /* run Julia commands */ + jl_eval_string("include(\"MyCode.jl\")"); + if (jl_exception_occurred()) { + std::cout << "=====> " << jl_typeof_str(jl_exception_occurred()) << std::endl; + } + + //---Construct the default run manager + auto runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial); + + //---Set mandatory initialization classes + // Detector construction + runManager->SetUserInitialization(new DetectorConstruction()); + + // Physics list + runManager->SetUserInitialization(new QBBC(0)); + + // User action initialization + runManager->SetUserInitialization(new ActionInitialization()); + + // Initialize G4 kernel + runManager->Initialize(); + + // Get the pointer to the User Interface manager + auto UImanager = G4UImanager::GetUIpointer(); + UImanager->ApplyCommand("/control/verbose 1"); + UImanager->ApplyCommand("/run/verbose 1"); + //UImanager->ApplyCommand("/event/verbose 0"); + //UImanager->ApplyCommand("/tracking/verbose 1"); + UImanager->ApplyCommand("/gun/particle e+"); + UImanager->ApplyCommand("/gun/energy 100 MeV"); + UImanager->ApplyCommand("/run/beamOn 100000"); + + // Job termination + delete runManager; + + // strongly recommended: notify Julia that the program is about to terminate. + // this allows Julia time to cleanup pending write requests and run all finalizers + jl_atexit_hook(0); +} diff --git a/docs/src/examples/JuliaAction_lit.jl b/docs/src/examples/JuliaAction_lit.jl new file mode 100644 index 0000000..794141e --- /dev/null +++ b/docs/src/examples/JuliaAction_lit.jl @@ -0,0 +1,94 @@ +# # Calling G4 actions in Julia +# +# This is a very simple example of calling user actions in Julia from a C++ +# Geant4 application. +# We define the user actions in Julia language in the file [`MyCode.jl`](@ref) +# and call them from the C++ application. The name and signatures of the functions +# are important since the C++ will associate them in the corresponding inherited +# classes. +# +# The C++ code is a single file [`G4example.cpp`](@ref) that defines the Geant4 the minimal +# set of classes to run a simulation. +# - The main program is responsible of initializing Julia by calling `julia_init` and +# loading the Julia code executing. +# ```cpp +# jl_init() +# jl_eval_string("include(\"MyCode.jl\")"); +# ``` +# - Each constructor of a user action class needs to initialize a C++ function pointer to the +# corresponding Julia function. This is done in the constructor to avoid any dynamic dispatch +# at runtime. For example, for the `EventAction` class: +# ```cpp +# typedef void (*eventaction_f)(const G4Event*); +# class EventAction : public G4UserEventAction { +# public: +# EventAction() { +# beginevent_jl = (eventaction_f)(jl_unbox_voidpointer(jl_eval_string("@cfunction(begin_of_event_action, Nothing, (CxxPtr{G4Event},))"))); +# endevent_jl = (eventaction_f)(jl_unbox_voidpointer(jl_eval_string("@cfunction(end_of_event_action, Nothing, (CxxPtr{G4Event},))"))); +# } +# ... +# private: +# eventaction_f beginevent_jl; +# eventaction_f endevent_jl; +# }; +# ``` +# - Finally the actions are called in the corresponding Geant4 classes. For example in the `EventAction` class: +# ```cpp +# void EventAction::BeginOfEventAction(const G4Event* event) { +# beginevent_jl(event); +# } +# ... +# ``` + +#md # !!! note "Note that" +#md # You can also download this example as a +#md # [Jupyter notebook](JuliaAction.ipynb) and a plain +#md # [Julia source file](JuliaAction.jl). +#md # +#md # The C++ code is available as a [source file](G4example.cpp) and the +#md # Julia code is available as a [source file](MyCode.jl). +#md # +#md # #### Table of contents +#md # ```@contents +#md # Pages = ["JuliaAction.md"] +#md # Depth = 2:3 +#md # ``` + +# ## Loading the necessary Julia modules +using Geant4_jll # Needed to locate the Geant4 installation directory +#md import DisplayAs: PNG #hide + +# ## Building G4Example Application +# The custom library is defined in the C++ file `G4example.cpp`. It is a single file to +# facilitate the building of the executable. +# +# The attribute `Geant4_jll.artifact_dir` provides the path to the Geant4 installation directory. +# Sources are in the same location as this script. +cd(@__DIR__) +g4prefix = Geant4_jll.artifact_dir +jlprefix = dirname(Sys.BINDIR); + +# We use the executables `geant4-config` and `julia-config.jl` to get the needed +# libraries and compiler/linker flags. + +g4libs = read(`$g4prefix/bin/geant4-config --libs`, String) |> split +filter!(x -> x != "-lG4gdml", g4libs) +jllibs = read(`$jlprefix/share/julia/julia-config.jl --ldlibs`, String) |> split +append!(jllibs, ["-L$jlprefix/lib"]) +cflags = read(`$g4prefix/bin/geant4-config --cflags`, String) |> split +ldflags = ["-Wl,-rpath,$g4prefix/lib", "-Wl,-rpath,$jlprefix/lib"]; +Sys.KERNEL == :Linux && append!(ldflags, ["-Wl,--no-as-needed"]); + +# Run the compilation and link command +Base.run(`c++ -O2 -fPIC $cflags -I$jlprefix/include/julia $ldflags $g4libs $jllibs + -o G4example $(@__DIR__)/G4example.cpp`).exitcode == 0 || error("Compilation failed"); + +# ## Run the application +# We need to set the variable `JULIA_PROJECT` pointing to correctly setup Julia environment. +withenv("JULIA_PROJECT" => "@.") do + Base.run(`./G4example`).exitcode == 0 || error("Execution failed"); +end + +# ## Display the results +println("=====> The file edepHist.png should have been saved") +#md # ![](edepHist.png) diff --git a/docs/src/examples/MyCode.jl b/docs/src/examples/MyCode.jl new file mode 100644 index 0000000..4bafe50 --- /dev/null +++ b/docs/src/examples/MyCode.jl @@ -0,0 +1,39 @@ +#---Load the needed Julia modules------------------------------------------------------------------ + +using Geant4 +using GeometryBasics +using FHist +using Plots + +println("=====> Loading MyCode.jl") + +edepHist = H1D("Event total Edep distribution", 100, 0., 110.) +edep = 0.0 + +function end_of_event_action(event) + push!(edepHist, edep) + return # This is mandatory to force to return nothing +end + +function begin_of_event_action(event) + global edep = 0.0 + return # This is mandatory to force to return nothing +end + +function stepping_action(step) + global edep += GetTotalEnergyDeposit(step) + return # This is mandatory to force to return nothing +end + +function begin_of_run_action(run) + return # This is mandatory to force to return nothing +end + +function end_of_run_action(run) + println("=====> End of run") + h = edepHist + img = plot(h.hist, title=h.title) + savefig(img, "edepHist.png") + println("=====> edepHist.png saved") + return # This is mandatory to force to return nothing +end diff --git a/docs/src/releases.md b/docs/src/releases.md index ee3c80d..bf8a912 100644 --- a/docs/src/releases.md +++ b/docs/src/releases.md @@ -2,6 +2,7 @@ ### 0.2.x - New Examples: - UserLib: how to build and call a user custom library providing additional Geant4 functionally that is not provided by the set of wrapped classes + - JuliaAction: emmbeding Julia in a C++ application. In this example we call user actions implemented in Julia ### 0.2.0 - New Features