Added UserLib example

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.

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

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();
+  }
+}
+

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
+}

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)

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)
+