Added JuliaAction example

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"

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)
 

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 <julia.h>
+#include <iostream>
+
+//---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);
+}

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)

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  

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