diff --git a/README.md b/README.md index e717513..9a6fb52 100644 --- a/README.md +++ b/README.md @@ -43,4 +43,7 @@ Example of a bubble chamber in which we display the particle tracks for an event Example from the original at https://github.com/settwi/g4-basic-scintillation and adapted to the Geant4.jl interface. Introduces optical photons and a custom physics list. It produces histograms as a result. ### AlephTPC -Example that shows how easy is to integrate several packages in Julia. The example uses the package `PYTHIA8` to generate LEP collisions of e+ e- which are then input as primary particles to a `Geant4` simulation. Only the event display is exercised of the simulation. +Example that shows how easy is to integrate several packages in Julia. The example uses the package `PYTHIA8` to generate LEP collisions of e+ e- which are then input as primary particles to a `Geant4` simulation. Only the event display is exercised of the simulation. + +### UserLib +Example on how to build and call a user custom library providing additional Geant4 functionally that is not provided by the set of wrapped classes. diff --git a/advanced/CustomLib/RoundCube.jl b/advanced/CustomLib/RoundCube.jl new file mode 100644 index 0000000..51ab839 --- /dev/null +++ b/advanced/CustomLib/RoundCube.jl @@ -0,0 +1,46 @@ +#---RoundCube definition--------------------------------------------------------------------------- +# A RoundCube is a cube with rounded edges and vertices. The cube is defined by the side length a +# and the radius of the rounded edges and vertices r. +function createRoundCube(a, r) + ca = a/4 + cube = G4Box("cube",ca,ca,ca) + cyl = G4Tubs("cyl",0, r, ca,0, 2π) + orb = G4Orb("orb", r) + + edge = G4Box("edge",r, r, ca) + vert = G4Box("vert",r, r, r) + + acyl1 = G4SubtractionSolid("ucyl", CxxPtr(edge), CxxPtr(cyl), G4Transform3D(G4RotationMatrix(), G4ThreeVector(-r,-r,0))) + acyl2 = G4SubtractionSolid("ucyl", CxxPtr(edge), CxxPtr(cyl), G4Transform3D(G4RotationMatrix(), G4ThreeVector(-r,r,0))) + acyl3 = G4SubtractionSolid("ucyl", CxxPtr(edge), CxxPtr(cyl), G4Transform3D(G4RotationMatrix(), G4ThreeVector(r,r,0))) + aorb = G4SubtractionSolid("uorb", CxxPtr(vert), CxxPtr(orb), G4Transform3D(G4RotationMatrix(), G4ThreeVector(-r,-r, -r))) + + t1 = G4Transform3D(G4RotationMatrix(), G4ThreeVector(ca,ca,0)) + t2 = G4Transform3D(G4RotationMatrix(0,π/2,0), G4ThreeVector(ca,0,ca)) + t3 = G4Transform3D(G4RotationMatrix(0,π/2,π/2), G4ThreeVector(0,ca,ca)) + t4 = G4Transform3D(G4RotationMatrix(), G4ThreeVector(ca,ca,ca)) + + s1 = G4SubtractionSolid("s1", CxxPtr(cube), CxxPtr(acyl1), t1) + s2 = G4SubtractionSolid("s2", CxxPtr(s1), CxxPtr(acyl2), t2) + s3 = G4SubtractionSolid("s3", CxxPtr(s2), CxxPtr(acyl3), t3) + s4 = G4SubtractionSolid("s4", CxxPtr(s3), CxxPtr(aorb), t4) + + t5 = G4Transform3D(G4RotationMatrix(0,0,0), G4ThreeVector(ca,ca,ca)) + t6 = G4Transform3D(G4RotationMatrix(0,π/2,0), G4ThreeVector(ca,ca,-ca)) + t7 = G4Transform3D(G4RotationMatrix(0,-π/2,0), G4ThreeVector(ca,-ca,ca)) + t8 = G4Transform3D(G4RotationMatrix(0,π/2,π/2), G4ThreeVector(ca,-ca,-ca)) + + u1 = G4UnionSolid("u1", CxxPtr(cube), CxxPtr(s4), t5) + u2 = G4UnionSolid("u2", CxxPtr(u1), CxxPtr(s4), t6) + u3 = G4UnionSolid("u3", CxxPtr(u2), CxxPtr(s4), t7) + u4 = G4UnionSolid("u4", CxxPtr(u3), CxxPtr(s4), t8) + + t9 = G4Transform3D(G4RotationMatrix(0,0,π), G4ThreeVector()) + u5 = G4UnionSolid("u5", CxxPtr(u4), CxxPtr(u4), t9) + + # remove finalizers + for a in (cube, cyl, orb, edge, vert, acyl1, acyl2, acyl3, aorb, s1, s2, s3, s4, u1, u2, u3, u4) + a.cpp_object = C_NULL + end + return u5 +end diff --git a/advanced/CustomLib/UserLib.cpp b/advanced/CustomLib/UserLib.cpp new file mode 100644 index 0000000..672c8da --- /dev/null +++ b/advanced/CustomLib/UserLib.cpp @@ -0,0 +1,87 @@ +// Example of a custom solid class that will be called by Julia +#include "G4Box.hh" +#include "G4Tubs.hh" +#include "G4Orb.hh" +#include "G4SubtractionSolid.hh" +#include "G4UnionSolid.hh" + +class RoundCube : public G4VSolid +{ +public: + RoundCube(double a, double r); + virtual ~RoundCube() { delete solid; } + + // Required methods for a custom solid + virtual EInside Inside(const G4ThreeVector& p) const override {return solid->Inside(p);} + virtual G4ThreeVector SurfaceNormal(const G4ThreeVector& p) const override {return solid->SurfaceNormal(p);} + virtual G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v) const override {return solid->DistanceToIn(p, v);} + virtual G4double DistanceToIn(const G4ThreeVector& p) const override {return solid->DistanceToIn(p);} + virtual G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v, + const G4bool calcNorm=false, G4bool *validNorm=0, + G4ThreeVector *n=0) const override {return solid->DistanceToOut(p, v, calcNorm, validNorm, n);} + virtual G4double DistanceToOut(const G4ThreeVector& p) const override {return solid->DistanceToOut(p);} + virtual G4GeometryType GetEntityType() const override {return "RoundCube";} + virtual std::ostream& StreamInfo(std::ostream& os) const override {os << "RoundCube with a = " << a << " r = " << r; return os;} + + virtual G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits& pVoxelLimits, const G4AffineTransform& pTransform, G4double& pMin, G4double& pMax) const override {return solid->CalculateExtent(pAxis, pVoxelLimits, pTransform, pMin, pMax);} + virtual void DescribeYourselfTo (G4VGraphicsScene& scene) const override {return solid->DescribeYourselfTo(scene);} + virtual void BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const override {return solid->BoundingLimits(pMin, pMax);} +private: + G4VSolid* solid; + double a, r; +}; + +RoundCube::RoundCube(double a, double r) : a(a), r(r), G4VSolid("RoundCube") { + G4double ca = a / 4; + G4VSolid* cube = new G4Box("cube", ca, ca, ca); + G4VSolid* cyl = new G4Tubs("cyl", 0, r, ca, 0, 2 * M_PI); + G4VSolid* orb = new G4Orb("orb", r); + + G4VSolid* edge = new G4Box("edge", r, r, ca); + G4VSolid* vert = new G4Box("vert", r, r, r); + G4VSolid* frame = new G4Box("frame", ca, ca, ca); + + G4VSolid* acyl1 = new G4SubtractionSolid("ucyl", edge, cyl, G4Transform3D(G4RotationMatrix(), G4ThreeVector(-r, -r, 0))); + G4VSolid* acyl2 = new G4SubtractionSolid("ucyl", edge, cyl, G4Transform3D(G4RotationMatrix(), G4ThreeVector(-r, r, 0))); + G4VSolid* acyl3 = new G4SubtractionSolid("ucyl", edge, cyl, G4Transform3D(G4RotationMatrix(), G4ThreeVector(r, r, 0))); + G4VSolid* aorb = new G4SubtractionSolid("uorb", vert, orb, G4Transform3D(G4RotationMatrix(), G4ThreeVector(-r, -r, -r))); + + G4Transform3D t1(G4RotationMatrix(), G4ThreeVector(ca, ca, 0)); + G4Transform3D t2(G4RotationMatrix(0, M_PI / 2, 0), G4ThreeVector(ca, 0, ca)); + G4Transform3D t3(G4RotationMatrix(0, M_PI / 2, M_PI / 2), G4ThreeVector(0, ca, ca)); + G4Transform3D t4(G4RotationMatrix(), G4ThreeVector(ca, ca, ca)); + + G4VSolid* s1 = new G4SubtractionSolid("s1", cube, acyl1, t1); + G4VSolid* s2 = new G4SubtractionSolid("s2", s1, acyl2, t2); + G4VSolid* s3 = new G4SubtractionSolid("s3", s2, acyl3, t3); + G4VSolid* s4 = new G4SubtractionSolid("s4", s3, aorb, t4); + + G4Transform3D t5(G4RotationMatrix(0, 0, 0), G4ThreeVector(ca, ca, ca)); + G4Transform3D t6(G4RotationMatrix(0, M_PI / 2, 0), G4ThreeVector(ca, ca, -ca)); + G4Transform3D t7(G4RotationMatrix(0, -M_PI / 2, 0), G4ThreeVector(ca, -ca, ca)); + G4Transform3D t8(G4RotationMatrix(0, M_PI / 2, M_PI / 2), G4ThreeVector(ca, -ca, -ca)); + + G4VSolid* u1 = new G4UnionSolid("u1", frame, s4, t5); + G4VSolid* u2 = new G4UnionSolid("u2", u1, s4, t6); + G4VSolid* u3 = new G4UnionSolid("u3", u2, s4, t7); + G4VSolid* u4 = new G4UnionSolid("u4", u3, s4, t8); + + solid = new G4UnionSolid("u5", u4, u4, G4Transform3D(G4RotationMatrix(0, 0, M_PI), G4ThreeVector())); +} + +extern "C" { + G4VSolid* createRoundCube(double a, double r) { + return new RoundCube(a, r); + } + void deleteRoundCube(G4VSolid* solid) { + delete solid; + } + const char* infoRoundCube(G4VSolid* solid) { + static std::string str; + std::ostringstream os; + solid->StreamInfo(os); + str = os.str(); + return str.c_str(); + } +} + diff --git a/advanced/CustomLib/UserLib.ipynb b/advanced/CustomLib/UserLib.ipynb new file mode 100644 index 0000000..584e77c --- /dev/null +++ b/advanced/CustomLib/UserLib.ipynb @@ -0,0 +1,160 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "# Calling C++ code from Julia\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", + "which is a cube with rounded edges and vertices. The cube is defined by the side length `a` and the\n", + "radius of the rounded edges and vertices `r`." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "## Loading the necessary Julia modules" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "using Geant4\n", + "using Geant4.SystemOfUnits\n", + "using Libdl\n", + "using CairoMakie, Rotations, LinearAlgebra, IGLWrap_jll # to force loading G4Vis extension" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "## Building the custom library\n", + "The custom library is defined in the C++ file `UserLib.cpp`.\n", + "The library provides a function to create a custom solid `RoundCube` and some additional functions\n", + "to interact with the solid.\n", + "\n", + "The attribute `Geant4_jll.artifact_dir` provides the path to the Geant4 installation directory.\n", + "We use only a sub-set of libraries needed to link shared library." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "prefix = Geant4.Geant4_jll.artifact_dir\n", + "dlext = Libdl.dlext\n", + "# Compilation of the custom library\n", + "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", + " -o UserLib.$dlext $(@__DIR__)/UserLib.cpp`).exitcode == 0 || error(\"Compilation failed\")" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "## Define Julia functions to interact with the custom library\n", + "The `@call` macro provides a very convenient way to call C functions (or extern \"C\").\n", + "We define the following functions to interact with the custom library in a more Julia-friendly way:" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "const lib = \"./UserLib.$(Libdl.dlext)\"\n", + "createRoundCube(a,r) = @ccall lib.createRoundCube(a::Float64, r::Float64)::CxxPtr{G4VSolid}\n", + "deleteRoundCube(s::CxxPtr{G4VSolid}) = @ccall lib.deleteRoundCube(s::CxxPtr{G4VSolid})::Cvoid\n", + "infoRoundCube(s::CxxPtr{G4VSolid}) = (@ccall lib.infoRoundCube(s::CxxPtr{G4VSolid})::Cstring) |> unsafe_string" + ], + "metadata": {}, + "execution_count": null + }, + { + "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" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "rcube = createRoundCube(10cm, 1cm)\n", + "img = drawDistanceToOut(rcube[], 100000)\n", + "display(\"image/png\", img)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Get the information about the `RoundCube`" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "info = infoRoundCube(rcube)\n", + "println(info)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Delete the `RoundCube`" + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "deleteRoundCube(rcube)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "---\n", + "\n", + "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" + ], + "metadata": {} + } + ], + "nbformat_minor": 3, + "metadata": { + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.1" + }, + "kernelspec": { + "name": "julia-1.11", + "display_name": "Julia 1.11.1", + "language": "julia" + } + }, + "nbformat": 4 +} diff --git a/advanced/CustomLib/UserLib.jl b/advanced/CustomLib/UserLib.jl new file mode 100644 index 0000000..1620348 --- /dev/null +++ b/advanced/CustomLib/UserLib.jl @@ -0,0 +1,26 @@ +using Geant4 +using Geant4.SystemOfUnits +using Libdl +using CairoMakie, Rotations, LinearAlgebra, IGLWrap_jll # to force loading G4Vis extension + +prefix = Geant4.Geant4_jll.artifact_dir +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 + -o UserLib.$dlext $(@__DIR__)/UserLib.cpp`).exitcode == 0 || error("Compilation failed") + +const lib = "./UserLib.$(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 + +rcube = createRoundCube(10cm, 1cm) +img = drawDistanceToOut(rcube[], 100000) +display(img) + +info = infoRoundCube(rcube) +println(info) + +deleteRoundCube(rcube) diff --git a/advanced/CustomLib/UserLib_lit.jl b/advanced/CustomLib/UserLib_lit.jl new file mode 100644 index 0000000..eb982be --- /dev/null +++ b/advanced/CustomLib/UserLib_lit.jl @@ -0,0 +1,62 @@ +# # Calling C++ code from Julia +# +# 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`, +# which is a cube with rounded edges and vertices. The cube is defined by the side length `a` and the +# radius of the rounded edges and vertices `r`. + +#md # !!! note "Note that" +#md # You can also download this example as a +#md # [Jupyter notebook](UserLib.ipynb) and a plain +#md # [Julia source file](UserLib.jl). +# +#md # #### Table of contents +#md # ```@contents +#md # Pages = ["UserLib.md"] +#md # Depth = 2:3 +#md # ``` + +# ## Loading the necessary Julia modules +using Geant4 +using Geant4.SystemOfUnits +using Libdl +using CairoMakie, Rotations, LinearAlgebra, IGLWrap_jll # to force loading G4Vis extension +#md import DisplayAs: PNG #hide + +# ## Building the custom library +# The custom library is defined in the C++ file `UserLib.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 +## 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") + +# ## 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)" +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) +#jl display(img) +#nb display("image/png", img) +#md PNG(img) +# Get the information about the `RoundCube` +info = infoRoundCube(rcube) +println(info) +# Delete the `RoundCube` +deleteRoundCube(rcube) +