Skip to content

Commit

Permalink
Added JuliaAction example
Browse files Browse the repository at this point in the history
  • Loading branch information
peremato committed Dec 3, 2024
1 parent 61572b6 commit 90584bb
Show file tree
Hide file tree
Showing 10 changed files with 572 additions and 21 deletions.
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ DisplayAs = "0b91fe84-8a4c-11e9-3e1d-67c38462b6d6"
FHist = "68837c9b-b678-4cd5-9925-8a54edc8f695"
GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
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"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand All @@ -22,6 +23,7 @@ Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
[compat]
FHist = "0.11"
Geant4 = "0.2"
Geant4_jll = "11.2"
Literate = "2.20"

[extras]
Expand Down
30 changes: 24 additions & 6 deletions advanced/CustomLib/UserLib.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
{
"cell_type": "markdown",
"source": [
"# Calling C++ code from Julia\n",
"# Calling Custom C++ library\n",
"\n",
"An example of calling user libraries that can provide additional Geant4 functionally that is not\n",
"provided by the set of wrapped classes. In this example we define a custom solid called `RoundCube`,\n",
Expand Down Expand Up @@ -49,8 +49,24 @@
"cell_type": "code",
"source": [
"prefix = Geant4.Geant4_jll.artifact_dir\n",
"dlext = Libdl.dlext\n",
"# Compilation of the custom library\n",
"dlext = Libdl.dlext;\n",
"# Compilation of the custom library"
],
"metadata": {},
"execution_count": null
},
{
"cell_type": "markdown",
"source": [
"The custom library is defined in the C++ file `UserLib.cpp`. Please note that the\n",
"callable functions are defined with the `extern \"C\"` attribute to avoid name mangling."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"Base.run(`c++ -O2 -shared -fPIC -std=c++17 -I$prefix/include/Geant4\n",
" -Wl,-rpath,$prefix/lib -L$prefix/lib\n",
" -lG4geometry -lG4materials -lG4global -lG4clhep\n",
Expand Down Expand Up @@ -84,16 +100,18 @@
"cell_type": "markdown",
"source": [
"## Testing the custom library\n",
"We create a `RoundCube` with side length `100` and radius `10` and draw the distance to the outside"
"We create a `RoundCube` with side length `100` and radius `10` and draw the distance to the outside of the solid\n",
"from a number of randomly distributed points in a random directions. This should get a nice image of the\n",
"surface `RoundCube`. It is exercising `Inside` and `DistanceToOut` methods of the custom solid."
],
"metadata": {}
},
{
"outputs": [],
"cell_type": "code",
"source": [
"rcube = createRoundCube(10cm, 1cm)\n",
"img = drawDistanceToOut(rcube[], 100000)\n",
"rcube = createRoundCube(10cm, 1cm) # returns a CxxPtr{G4VSolid}\n",
"img = drawDistanceToOut(rcube[], 100000) # implemented in G4Vis ext. It expects a G4VSolid.\n",
"display(\"image/png\", img)"
],
"metadata": {},
Expand Down
7 changes: 4 additions & 3 deletions advanced/CustomLib/UserLib.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,9 @@ using Libdl
using CairoMakie, Rotations, LinearAlgebra, IGLWrap_jll # to force loading G4Vis extension

prefix = Geant4.Geant4_jll.artifact_dir
dlext = Libdl.dlext
dlext = Libdl.dlext;
# Compilation of the custom library

Base.run(`c++ -O2 -shared -fPIC -std=c++17 -I$prefix/include/Geant4
-Wl,-rpath,$prefix/lib -L$prefix/lib
-lG4geometry -lG4materials -lG4global -lG4clhep
Expand All @@ -16,8 +17,8 @@ createRoundCube(a,r) = @ccall lib.createRoundCube(a::Float64, r::Float64)::CxxPt
deleteRoundCube(s::CxxPtr{G4VSolid}) = @ccall lib.deleteRoundCube(s::CxxPtr{G4VSolid})::Cvoid
infoRoundCube(s::CxxPtr{G4VSolid}) = (@ccall lib.infoRoundCube(s::CxxPtr{G4VSolid})::Cstring) |> unsafe_string

rcube = createRoundCube(10cm, 1cm)
img = drawDistanceToOut(rcube[], 100000)
rcube = createRoundCube(10cm, 1cm) # returns a CxxPtr{G4VSolid}
img = drawDistanceToOut(rcube[], 100000) # implemented in G4Vis ext. It expects a G4VSolid.
display(img)

info = infoRoundCube(rcube)
Expand Down
32 changes: 20 additions & 12 deletions advanced/CustomLib/UserLib_lit.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# # Calling C++ code from Julia
# # Calling Custom C++ library
#
# An example of calling user libraries that can provide additional Geant4 functionally that is not
# provided by the set of wrapped classes. In this example we define a custom solid called `RoundCube`,
Expand All @@ -24,36 +24,44 @@ using CairoMakie, Rotations, LinearAlgebra, IGLWrap_jll # to force loading G4Vi
#md import DisplayAs: PNG #hide

# ## Building the custom library
# The custom library is defined in the C++ file `UserLib.cpp`.
# The custom library is defined in the C++ file `UserLibrary.cpp`.
# The library provides a function to create a custom solid `RoundCube` and some additional functions
# to interact with the solid.
#
# The attribute `Geant4_jll.artifact_dir` provides the path to the Geant4 installation directory.
# We use only a sub-set of libraries needed to link shared library.

prefix = Geant4.Geant4_jll.artifact_dir
dlext = Libdl.dlext
dlext = Libdl.dlext;
if Sys.KERNEL == :Linux
ldflags = "-Wl,-rpath,$prefix/lib -Wl,--no-as-needed"
else
ldflags = "-Wl,-rpath,$prefix/lib -Wl"
end
## Compilation of the custom library
Base.run(`c++ -O2 -shared -fPIC -std=c++17 -I$prefix/include/Geant4
-Wl,-rpath,$prefix/lib -L$prefix/lib
-lG4geometry -lG4materials -lG4global -lG4clhep
-o UserLib.$dlext $(@__DIR__)/UserLib.cpp`).exitcode == 0 || error("Compilation failed")
# The custom library is defined in the C++ file [`UserLibrary.cpp`](@ref). Please note that the
# callable functions are defined with the `extern "C"` attribute to avoid name mangling.
Base.run(`c++ -O2 -shared -fPIC -std=c++17 -I$prefix/include/Geant4 $ldflags
-L$prefix/lib -lG4geometry -lG4materials -lG4global -lG4clhep
-o UserLibrary.$dlext $(@__DIR__)/UserLibrary.cpp`).exitcode == 0 || error("Compilation failed")

# ## Define Julia functions to interact with the custom library
# The `@call` macro provides a very convenient way to call C functions (or extern "C").
# We define the following functions to interact with the custom library in a more Julia-friendly way:
const lib = "./UserLib.$(Libdl.dlext)"
const lib = "./UserLibrary.$(Libdl.dlext)"
createRoundCube(a,r) = @ccall lib.createRoundCube(a::Float64, r::Float64)::CxxPtr{G4VSolid}
deleteRoundCube(s::CxxPtr{G4VSolid}) = @ccall lib.deleteRoundCube(s::CxxPtr{G4VSolid})::Cvoid
infoRoundCube(s::CxxPtr{G4VSolid}) = (@ccall lib.infoRoundCube(s::CxxPtr{G4VSolid})::Cstring) |> unsafe_string

# ## Testing the custom library
# We create a `RoundCube` with side length `100` and radius `10` and draw the distance to the outside
rcube = createRoundCube(10cm, 1cm)
img = drawDistanceToOut(rcube[], 100000)
# We create a `RoundCube` with side length `100` and radius `10` and draw the distance to the outside of the solid
# from a number of randomly distributed points in a random directions. This should get a nice image of the
# surface `RoundCube`. It is exercising `Inside` and `DistanceToOut` methods of the custom solid.
rcube = createRoundCube(10cm, 1cm) # returns a CxxPtr{G4VSolid}
img = drawDistanceToOut(rcube[], 100000) # implemented in G4Vis ext. It expects a G4VSolid.
#jl display(img)
#nb display("image/png", img)
#md PNG(img)
#md PNG(img) #hide
# Get the information about the `RoundCube`
info = infoRoundCube(rcube)
println(info)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include "G4SubtractionSolid.hh"
#include "G4UnionSolid.hh"

// The custom solid class needs to inherit from G4VSolid and implement the required methods
class RoundCube : public G4VSolid
{
public:
Expand All @@ -31,6 +32,7 @@ class RoundCube : public G4VSolid
double a, r;
};

// Constructor of the custom solid (RoundCube) implemented as a union and subtractions of several solids
RoundCube::RoundCube(double a, double r) : a(a), r(r), G4VSolid("RoundCube") {
G4double ca = a / 4;
G4VSolid* cube = new G4Box("cube", ca, ca, ca);
Expand Down Expand Up @@ -69,6 +71,7 @@ RoundCube::RoundCube(double a, double r) : a(a), r(r), G4VSolid("RoundCube") {
solid = new G4UnionSolid("u5", u4, u4, G4Transform3D(G4RotationMatrix(0, 0, M_PI), G4ThreeVector()));
}

// Define the callable functions for the custom solid
extern "C" {
G4VSolid* createRoundCube(double a, double r) {
return new RoundCube(a, r);
Expand Down
170 changes: 170 additions & 0 deletions advanced/EmbedJulia/G4example.cpp
Original file line number Diff line number Diff line change
@@ -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);
}
Loading

0 comments on commit 90584bb

Please sign in to comment.