From 5671ce71359ccd5f06de391fe35eb70de7f11d46 Mon Sep 17 00:00:00 2001 From: Pere Mato Date: Wed, 17 Apr 2024 17:30:38 +0200 Subject: [PATCH] Completed the set of notebooks --- Project.toml | 1 + tutorial/docs/02-wrapped-classes.ipynb | 4 +- tutorial/docs/04-physics-list.ipynb | 36 +- tutorial/docs/05-primary-particles.ipynb | 363 ++++++++++++++- tutorial/docs/06-field.ipynb | 143 +++++- tutorial/docs/07-applications.ipynb | 338 +++++++++++++- tutorial/docs/08-sensitive-detectors.ipynb | 227 ++++++++- tutorial/docs/09-scoring-meshes.ipynb | 195 +++++++- tutorial/docs/10-histograms.ipynb | 18 +- tutorial/docs/11-event-display.ipynb | 138 +++++- tutorial/docs/VisSettings.jl | 9 + tutorial/docs/_toc.yml | 2 + .../docs/examples/TestEm3/DetectorTestEm3.jl | 260 +++++++++++ tutorial/docs/examples/TestEm3/TestEm3.ipynb | 439 ++++++++++++++++++ 14 files changed, 2087 insertions(+), 86 deletions(-) create mode 100644 tutorial/docs/VisSettings.jl create mode 100644 tutorial/docs/examples/TestEm3/DetectorTestEm3.jl create mode 100644 tutorial/docs/examples/TestEm3/TestEm3.ipynb diff --git a/Project.toml b/Project.toml index 7cabd19..c053e4b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,5 +1,6 @@ [deps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" FHist = "68837c9b-b678-4cd5-9925-8a54edc8f695" Geant4 = "559df036-b7a0-42fd-85df-7d5dd9d70f44" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" diff --git a/tutorial/docs/02-wrapped-classes.ipynb b/tutorial/docs/02-wrapped-classes.ipynb index a8c69e2..3cef929 100644 --- a/tutorial/docs/02-wrapped-classes.ipynb +++ b/tutorial/docs/02-wrapped-classes.ipynb @@ -79,7 +79,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The type hierarchy is a follows:\n", + "The type hierarchy is as follows:\n", "```\n", " +-------+\n", " | Any |\n", @@ -422,7 +422,7 @@ "box1 = box2 = zeros(1000)\n", "GC.gc()\n", "@show GetCubicVolume(union) # G4BooleanSolid cashes the volume\n", - "@show DistanceToIn(union, G4ThreeVector(10,10,10)); # this will probabluy crash the program" + "@show DistanceToIn(union, G4ThreeVector(10,10,10)); # now this will not crash the program" ] } ], diff --git a/tutorial/docs/04-physics-list.ipynb b/tutorial/docs/04-physics-list.ipynb index 22e2987..3d8c5e0 100644 --- a/tutorial/docs/04-physics-list.ipynb +++ b/tutorial/docs/04-physics-list.ipynb @@ -51,34 +51,9 @@ "metadata": {}, "outputs": [], "source": [ - "#---The next will be hidden in new versions-------------------------\n", - "# create a empty world filled with vacumm\n", - "struct World <: G4JLDetector end\n", - "using Geant4.PhysicalConstants: universe_mean_density\n", - "using Geant4.SystemOfUnits: g, mole, kelvin, pascal\n", - "function bigbang(::World)\n", - " vacuum = G4Material(\"Vacuum\", z=1., a=1.01g/mole, density=universe_mean_density, state=kStateGas, \n", - " temperature=2.73*kelvin, pressure=3.e-18*pascal)\n", - " G4PVPlacement(nothing, G4ThreeVector(),\n", - " G4LogicalVolume(G4Box(\"world\", 1000,1000,1000), move!(vacuum), \"World\"),\n", - " \"World\", nothing, false, 0, false)\n", - "end \n", - "Geant4.getConstructor(::World) = bigbang\n", - "#---end-------------------------------------------------------------" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "\n", "app = G4JLApplication(\n", - " detector = World(), # should be the default\n", - " generator = G4JLGunGenerator(), # should be the default\n", " physics_type = FTFP_BERT\n", - ");" + ")" ] }, { @@ -120,10 +95,8 @@ "metadata": {}, "outputs": [], "source": [ - "ui`/run/verbose 1`\n", - "app = nothing\n", - "GC.gc()\n", - "G4RunManager!GetRunManager()" + "ui`/control/verbose 1`\n", + "ui`/run/verbose 1`" ] }, { @@ -149,9 +122,8 @@ "outputs": [], "source": [ "app = G4JLApplication(\n", - " detector = World(), # should be the default\n", - " generator = G4JLGunGenerator(), # should be the default\n", " physics_type = MyPhysicsList,\n", + " verbose=1,\n", ")\n", "configure(app)\n", "initialize(app)" diff --git a/tutorial/docs/05-primary-particles.ipynb b/tutorial/docs/05-primary-particles.ipynb index 1f8dbfa..f221fa5 100644 --- a/tutorial/docs/05-primary-particles.ipynb +++ b/tutorial/docs/05-primary-particles.ipynb @@ -4,24 +4,371 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Defining Primary Particles" + "# Defining Primary Particles\n", + "\n", + "A Geant4 application requires a source of primary particles. We have several possibilities to provide a primary particle generator.\n", + "\n", + "- [**Particle Gun**](#gun) - This is the simplest primary particle generator. It is a single type of particle with a fix energy and direction \n", + "- [**General Particle Source**](#gps) - This is a generalization of a particle gun with many more parameters to configure its behavior. See the [G4GeneralParticleSource documentation](https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/GettingStarted/generalParticleSource.html).\n", + "- **Full customization** - The user can write code to decide the primary particles to create, as well as, their energies and directions for each event. \n", + "\n", + "We are going to exercise the three cases in this notebook" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, + "outputs": [], + "source": [ + "using Geant4\n", + "using Geant4.SystemOfUnits" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Do some analysis of energy, positions and directions\n", + "using DataFrames, Plots \n", + "\n", + "function generate(generator, N)\n", + " df = DataFrame(energy = Float32[], \n", + " pos_x = Float32[], pos_y = Float32[], pos_z = Float32[],\n", + " dir_x = Float32[], dir_y = Float32[], dir_z = Float32[])\n", + " for i in 1:N\n", + " evt = G4Event()\n", + " generator.gen_method(evt, generator.data)\n", + " pos = evt |> GetPrimaryVertex |> GetPosition\n", + " dir = evt |> GetPrimaryVertex |> GetPrimary |> GetMomentumDirection\n", + " ene = evt |> GetPrimaryVertex |> GetPrimary |> GetKineticEnergy\n", + " push!(df, (ene, x(pos), y(pos), z(pos), x(dir), y(dir), z(dir)))\n", + " end\n", + " return df\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "## Particle Gun\n", + "We have introduced a new type `G4JLGunGenerator` to facilitate the task of defining a " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "particlegun = G4JLGunGenerator(particle = \"pi+\", \n", + " energy = 330MeV, \n", + " direction = G4ThreeVector(0, 0, 1), \n", + " position = G4ThreeVector(0, 0, 0));" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The actual `G4ParticleGun` object (the `data.gun` field) has not been created yet. It will be done by the `G4RunManager` during initialization. If we do initialize by hand now we will get a error since particles have not been created yet." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Lets build an application with just the generator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "app = G4JLApplication(\n", + " generator = particlegun,\n", + ")\n", + "configure(app)\n", + "initialize(app)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Run one event and see the printout of the initial particle" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ui`/tracking/verbose 1`\n", + "beamOn(app,1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can easily change the parameters of the particle gun with the following methods" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "SetParticleByName(particlegun, \"e-\")\n", + "SetParticleEnergy(particlegun, 10GeV)\n", + "SetParticleMomentumDirection(particlegun, G4ThreeVector(1,0,0))\n", + "\n", + "beamOn(app,1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "## General Particle Source\n", + "You can create a `G4JLGeneralParticleSource` with the same parameters (and same behavior) as the Particle Gun." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "gps = G4JLGeneralParticleSource(particle = \"proton\", \n", + " energy = 10MeV, \n", + " direction = G4ThreeVector(1,0,0), \n", + " position = G4ThreeVector(0,0,0));" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "app = G4JLApplication(\n", + " generator = gps,\n", + ")\n", + "configure(app)\n", + "initialize(app)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df = generate(gps, 1000)\n", + "histogram(df.energy)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Note that** The exception is due to the fact that Geant4 does not like to re-initialize. Please ignore the exception and look at the created particle. Now we can change the parameters in two different ways:\n", + "- using the UI commands\n", + "- instantiating a new `G4JLGeneralParticleSource` with the needed parameters " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "reinitialize(app.generator; particle=\"geantino\",\n", + " energy=2MeV,\n", + " pos=(type=\"Point\", centre=G4ThreeVector(1cm,2cm,1cm)),\n", + " ang=(type=\"iso\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], - "source": [] + "source": [ + "df = generate(gps, 1000)\n", + "histogram(df.energy)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "reinitialize(gps; particle=\"geantino\",\n", + " ene=(type=\"Lin\", min=2MeV, max=10MeV, gradient=1, intercept=1),\n", + " pos=(type=\"Plane\", shape=\"Square\", centre=G4ThreeVector(1cm,2cm,1cm), halfx=2cm, halfy=2cm),\n", + " ang=(type=\"cos\",mintheta=10deg, maxtheta=80deg))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df = generate(gps, 1000)\n", + "histogram(df.energy)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scatter(df.pos_x, df.pos_y, df.pos_z)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "## Custom Primary Generator\n", + "The user can provide fully custom primary particle generator. This is done by defining custom structure to configure the generator and two functions to initialize and generate the primary particles that will be called by the `G4RunManager` for each event.\n", + "\n", + "The user configuration data structure should inherit from `G4JLGeneratorData` abstract type. Here is an example of a generator of a for a rectangle shape origined a mono-energy particles. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# define the data structure with the generator parameters\n", + "mutable struct PlaneSourceData <: G4JLGeneratorData\n", + " particleName::String # particle type\n", + " particlePtr::CxxPtr{G4ParticleDefinition} # keep the pointer to definition for performance\n", + " energy::Float64 # kinetic energy\n", + " halfx::Float64 # rectangle dimensions\n", + " halfy::Float64\n", + " position::G4ThreeVector # rectangle origin\n", + " direction::G4ThreeVector # particle direction\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We defile now a function `PlaneSource` that will create an instance of `G4JLPrimaryGenerator` with the configuration data (using some default values) and the two functions: `init` and `generate`. The `init` converts the particle name to a definition, which is then used for the `generate` function at each event. The `generate` generate random points en X and Y and creates a `G4PrimaryParticle` and `G4PrimaryVertex` which are added to the event given as argument. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# define the constructor with the default parameters\n", + "function PlaneSource(;particle=\"gamma\", energy=0.07MeV, halfx=7cm, halfy=7cm, \n", + " position=G4ThreeVector(0,0,0), \n", + " direction=G4ThreeVector(0,0,1))\n", + " data = PlaneSourceData(particle, CxxPtr{G4ParticleDefinition}(C_NULL), energy, halfx, halfy, position, direction)\n", + "\n", + " function init(data:: PlaneSourceData, app::G4JLApplication)\n", + " data.particlePtr = FindParticle(data.particleName)\n", + " end\n", + "\n", + " function generate( evt::G4Event, data:: PlaneSourceData)::Nothing\n", + " mass = data.particlePtr |> GetPDGMass\n", + " momentum = √((mass + data.energy)^2 - mass^2)\n", + " pvec = momentum * data.direction\n", + " pos = data.position + G4ThreeVector( data.halfx * (rand() - 0.5), data.halfy * (rand() - 0.5), 0)\n", + " primary = G4PrimaryParticle(data.particlePtr, pvec |> x, pvec |> y, pvec |> z )\n", + " vertex = G4PrimaryVertex(pos, 0ns)\n", + " SetPrimary(vertex, move!(primary)) # note that we give up ownership of the objects just created\n", + " AddPrimaryVertex(evt, move!(vertex)) # note that we give up ownership of the objects just created\n", + " end\n", + " G4JLPrimaryGenerator(\"PlaneSource\", data; init_method=init, generate_method=generate)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "planesource = PlaneSource(energy=10MeV, halfx=10cm, halfy=10cm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "app = G4JLApplication(\n", + " generator = planesource,\n", + ")\n", + "configure(app)\n", + "initialize(app)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df = generate(planesource, 1000)\n", + "histogram(df.energy)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scatter(df.pos_x, df.pos_y)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scatter(df.pos_x, df.pos_y, df.pos_z)" + ] } ], "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, "language_info": { - "name": "python" + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" } }, "nbformat": 4, diff --git a/tutorial/docs/06-field.ipynb b/tutorial/docs/06-field.ipynb index 653ed28..bcb281f 100644 --- a/tutorial/docs/06-field.ipynb +++ b/tutorial/docs/06-field.ipynb @@ -4,24 +4,151 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Defining Magnetic Field" + "# Defining Magnetic Field\n", + "\n", + "The user can define either an uniform magnetic field or a custom one." ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, + "outputs": [], + "source": [ + "using Geant4\n", + "using Geant4.SystemOfUnits\n", + "using Geant4.PhysicalConstants" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "## Uniform magnetic field\n", + "The function`G4JLUniformMagField` creates a uniform magnetic field proving a **B** vector. For example: " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bfield = G4JLUniformMagField(G4ThreeVector(0,0, 1.5tesla))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The G4 toolkit will call the function `bfield.getfield(field::G4ThreeVector, pos::G4ThreeVector, data::G4JLFieldData)` each time that is needed. The method `getfield(bfield::G4JLMagneticField, pos::G4ThreeVector)` does it as well." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Will be moved to Geant4.jl\n", + "function getfield(bfield::G4JLMagneticField, pos::G4ThreeVector)\n", + " B = G4ThreeVector()\n", + " bfield.getfield(B, pos, bfield.data)\n", + " return B\n", + "end " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@show getfield(bfield, G4ThreeVector(0,0,0))/tesla" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Custom magnetic Field\n", + " To define a custom field we need:\n", + "- define first a user structure for the parameters that inherits from the abstract type `G4JLFieldData`\n", + "- then, define a function with the signature `(result::G4ThreeVector, position::G4ThreeVector, params::G4JLFieldData)::Nothing`\n", + "- and finally, with all this, instantiate the magnetic field calling the function \n", + " ```\n", + " G4JLMagneticField(, ; getfield_method=)\n", + " ```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], - "source": [] + "source": [ + "using Geant4.SystemOfUnits: ampere\n", + "using Parameters\n", + "@with_kw mutable struct WireFieldData <: G4JLFieldData\n", + " I::Float64 = 1ampere\n", + " wiredir::G4ThreeVector = G4ThreeVector(0,0,1)\n", + "end\n", + "\n", + "function getfield!(field::G4ThreeVector, pos::G4ThreeVector, data::WireFieldData)::Nothing\n", + " r = cross(data.wiredir, pos)\n", + " B = (mu0 * data.I)/(2π*mag2(r)) * r \n", + " assign(field, B)\n", + " return\n", + "end\n", + "\n", + "bfield = G4JLMagneticField(\"WireField\", WireFieldData(); getfield_method=getfield!);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@show getfield(bfield, G4ThreeVector(1,0,0))/tesla\n", + "@show getfield(bfield, G4ThreeVector(2,0,0))/tesla\n", + "@show getfield(bfield, G4ThreeVector(3,0,0))/tesla\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Plots\n", + "\n", + "# Range of distances from the wire\n", + "r_values = range(0.01m, stop=1m, length=100) # Distances from 0.01m to 1m\n", + "\n", + "# Calculate magnetic field strengths corresponding to each distance\n", + "B_values = [ mag(getfield(bfield, G4ThreeVector(r,0,0)))/tesla for r in r_values]\n", + "\n", + "# Plot\n", + "plot(r_values, B_values, xlabel=\"Distance (mm)\", ylabel=\"Magnetic Field (T)\", label=\"\", legend=:bottomright)\n", + "title!(\"Magnetic Field vs. Distance from Wire\")" + ] } ], "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, "language_info": { - "name": "python" + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" } }, "nbformat": 4, diff --git a/tutorial/docs/07-applications.ipynb b/tutorial/docs/07-applications.ipynb index e23005a..8a69ad7 100644 --- a/tutorial/docs/07-applications.ipynb +++ b/tutorial/docs/07-applications.ipynb @@ -4,24 +4,346 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Building Applications" + "# Building Applications\n", + "What is still missing to build complete applications is to add user actions or sensitive detectors or scorers. In this tutorial we have a look at the user actions.\n", + "\n", + "## User Actions\n", + "User actions are native Julia functions that will be called by the Geant4 toolkit at the defined moment during the simulation. They\n", + "\n", + "The following are the possible user actions:\n", + "- **stepping action**. Called on each simulation step. The signature is `(::G4Step, ::G4JLApplication)::Nothing`. Consult the [G4Step](https://geant4.kek.jp/Reference/11.2.0/classG4Step.html) reference manual to see what you can get from it. \n", + "- **pre-tracking action**. Called at the creation of a new participle being tracked. The signature is `(::G4Track, ::G4JLApplication)::Nothing`. Consult the [G4Step](https://geant4.kek.jp/Reference/11.2.0/classG4Track.html) reference manual to see what you can get from it.\n", + "- **post-tracking action**. Called at the end of the particle being tracked. The signature is `(::G4Track, ::G4JLApplication)::Nothing`. Consult the [G4Track](https://geant4.kek.jp/Reference/11.2.0/classG4Track.html) reference manual to see what you can get from it.\n", + "- **begin-event action**. Called at the beginning of each event. The signature is `(::G4Event, ::G4JLApplication)::Nothing`. Consult the [G4Event](https://geant4.kek.jp/Reference/11.2.0/classG4Event.html) reference manual to see what you can get from it.\n", + "- **end-event action**. Called at the end of each event. The signature is `(::G4Event, ::G4JLApplication)::Nothing`. Consult the [G4Event](https://geant4.kek.jp/Reference/11.2.0/classG4Event.html) reference manual to see what you can get from it.\n", + "- **begin-run action**. Called at the beginning of a run. The signature is `(::G4Run, ::G4JLApplication)::Nothing`. Consult the [G4Run](https://geant4.kek.jp/Reference/11.2.0/classG4Run.html) reference manual to see what you can get from it.\n", + "- **end-run action**. Called at the end of a run. The signature is `(::G4Run, ::G4JLApplication)::Nothing`. Consult the [G4Run](https://geant4.kek.jp/Reference/11.2.0/classG4Run.html) reference manual to see what you can get from it.\n", + "- **stack action**. Called when a particle is going to be put back on the particle stack. The signature is `(::G4Track, ::G4JLApplication)::G4ClassificationOfNewTrack`. Consult the [G4Track](https://geant4.kek.jp/Reference/11.2.0/classG4Track.html) reference manual to see what you can get from it.\n", + "\n", + "All user actions receive in addition to the standard G4 arguments a reference to the G4JLApplication object from which the user can get additional information on the application. These are the available field:\n", + "- **detector** - user defined detector parameters object\n", + "- **simdata** - vector of the user defined simulation data objects (one per worker thread plus the accumulated one). The function `getSIMdata(app)` return the simulation data corresponding to the `threadid` of the worker calling the user action.\n", + "- **generator** - primary particle generator\n", + "- **field** - magnetic field\n", + "- **nthreads** - number of worker threads\n", + "- **verbose** - verbosity level\n", + "\n", + "**Ploase note that:** user actions need to be coded as thread-safe (global data should not be modified without any protection). However, this notebook not use multi-threading in order to be able to re-configure the `G4JLApplication` with different actions. \n", + "\n", + "\n" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "vscode": { - "languageId": "julia" + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "using Geant4" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's begin with something very simple. To define a begin event action that prints a message (notice that we use `G4JL_println` instead of the native `println` to ensure thread safety.)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "beginevent (generic function with 1 method)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "function beginevent(evt::G4Event, app::G4JLApplication)\n", + " G4JL_println(\"started event $(evt |> GetEventID)\")\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "**************************************************************\n", + " Geant4 version Name: geant4-11-02 [MT] (8-December-2023)\n", + " Copyright : Geant4 Collaboration\n", + " References : NIM A 506 (2003), 250-303\n", + " : IEEE-TNS 53 (2006), 270-278\n", + " : NIM A 835 (2016), 186-225\n", + " WWW : http://geant4.org/\n", + "**************************************************************\n", + "\n" + ] } - }, + ], + "source": [ + "app = G4JLApplication(begineventaction_method=beginevent)\n", + "configure(app)\n", + "initialize(app)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "started event 0\n", + "started event 1\n", + "started event 2\n", + "started event 3\n", + "started event 4\n" + ] + } + ], + "source": [ + "\n", + "beamOn(app,5)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Please notice that the output message would be prefixed by the working thread ID and the order of execution of the events is not sequential in case we enable multi-threading. " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "endevent (generic function with 1 method)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "function endevent(evt::G4Event, app::G4JLApplication)\n", + " G4JL_println(\"end event $(evt |> GetEventID)\")\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, "outputs": [], - "source": [] + "source": [ + "app = G4JLApplication(begineventaction_method = beginevent,\n", + " endeventaction_method = endevent)\n", + "configure(app)\n", + "initialize(app)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "started event 0\n", + "end event 0\n", + "started event 1\n", + "end event 1\n", + "started event 2\n", + "end event 2\n", + "started event 3\n", + "end event 3\n", + "started event 4\n", + "end event 4\n" + ] + } + ], + "source": [ + "beamOn(app,5)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "stepaction (generic function with 1 method)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "function stepaction(step::G4Step, app::G4JLApplication)\n", + " G4JL_println(\"step with length $(step |> GetStepLength)\")\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "started event 0\n", + "step with length 1000.0\n", + "end event 0\n", + "started event 1\n", + "step with length 1000.0\n", + "end event 1\n", + "started event 2\n", + "step with length 1000.0\n", + "end event 2\n", + "started event 3\n", + "step with length 1000.0\n", + "end event 3\n", + "started event 4\n", + "step with length 1000.0\n", + "end event 4\n" + ] + } + ], + "source": [ + "app = G4JLApplication(begineventaction_method = beginevent,\n", + " endeventaction_method = endevent,\n", + " stepaction_method = stepaction)\n", + "configure(app)\n", + "initialize(app)\n", + "beamOn(app,5)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining Simulation Data\n", + "\n", + "The user can define a custom data structure to collect simulation data during the execution of the user actions. Here is an example:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "started event 0\n", + "end event 0\n", + "started event 1\n", + "end event 1\n", + "started event 2\n", + "end event 2\n", + "started event 3\n", + "end event 3\n", + "started event 4\n", + "end event 4\n" + ] + } + ], + "source": [ + "using DataFrames\n", + "mutable struct MySimData <: G4JLSimulationData\n", + " steps::DataFrame\n", + " MySimData() = new(DataFrame( X = Float64[], Y = Float64[], Z = Float64[],\n", + " KinE = Float64[], dE = Float64[], StepLeng = Float64[]))\n", + "end\n", + "function stepaction(step::G4Step, app::G4JLApplication)\n", + " steps = getSIMdata(app).steps # Get it own copy of the simdata (avoid data racing in MT mode)\n", + " point = step |> GetPreStepPoint |> GetPosition\n", + " push!(steps, (point |> x, point |> y, point |> z,\n", + " step |> GetTrack |> GetKineticEnergy, step |> GetTotalEnergyDeposit, step |> GetStepLength ))\n", + " return\n", + "end\n", + "\n", + "app = G4JLApplication(simdata = MySimData(),\n", + " begineventaction_method = beginevent,\n", + " endeventaction_method = endevent,\n", + " stepaction_method = stepaction)\n", + "configure(app)\n", + "initialize(app)\n", + "beamOn(app,5)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To get the results we can access `app.simdata[1]`" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "MySimData(\u001b[1m5×6 DataFrame\u001b[0m\n", + "\u001b[1m Row \u001b[0m│\u001b[1m X \u001b[0m\u001b[1m Y \u001b[0m\u001b[1m Z \u001b[0m\u001b[1m KinE \u001b[0m\u001b[1m dE \u001b[0m\u001b[1m StepLeng \u001b[0m\n", + " │\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\u001b[90m Float64 \u001b[0m\n", + "─────┼───────────────────────────────────────────────────────────\n", + " 1 │ 0.0 0.0 0.0 10.0 3.03438e-23 1000.0\n", + " 2 │ 0.0 0.0 0.0 10.0 3.03438e-23 1000.0\n", + " 3 │ 0.0 0.0 0.0 10.0 3.03438e-23 1000.0\n", + " 4 │ 0.0 0.0 0.0 10.0 3.03438e-23 1000.0\n", + " 5 │ 0.0 0.0 0.0 10.0 3.03438e-23 1000.0)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "app.simdata[1]" + ] } ], "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, "language_info": { - "name": "python" + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" } }, "nbformat": 4, diff --git a/tutorial/docs/08-sensitive-detectors.ipynb b/tutorial/docs/08-sensitive-detectors.ipynb index 0a933f0..a2ca8fa 100644 --- a/tutorial/docs/08-sensitive-detectors.ipynb +++ b/tutorial/docs/08-sensitive-detectors.ipynb @@ -4,24 +4,237 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Sensitive Detectors" + "# Sensitive Detectors\n", + "The user can define sensitive detectors by defining a custom data structure and 3 callback functions, which will initialize, fill and dispose the defined data structure. Later, the instantiated sensitive detector would be associated to one or more logical volumes of the detector setup.\n", + "\n", + "Instantiating a `G4JLSensitiveDetector` will require the following arguments:\n", + "- **name**. A string to identify the sensitive detector. No default.\n", + "- **sd data**. A instance of a user defined `G4JLSDData` mutable data structure that will passed to each callback invocation.\n", + "- **initialize method**. User method that is called at the beginning of each event. The signature is `(::B2aSDData, ::G4HCofThisEvent)::Nothing`.\n", + "- **endOfEvent method**. User method that is called at the end of each event. The signature is `(::B2aSDData, ::G4HCofThisEvent)::Nothing`.\n", + "- **processHits method**. User method that is called at simulation step that ends at the associated logical volume. The signature is `(::B2aSDData, ::G4Step, ::G4TouchableHistory)::Bool`. Consult the [G4Step](https://geant4.kek.jp/Reference/v11.1.1/classG4Step.html) reference manual to see what you can get from the G4Step. It returns true if a true hit is generated.\n" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, + "outputs": [], + "source": [ + "using Geant4\n", + "using Geant4.SystemOfUnits" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "## Sensitive Detector Data\n", + "Let's define something very simple, for example collecting generated `Hits` (i.e deposit energies) on a sensitive detector. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Hit structure------------------------------------------------------------------\n", + "@enum HitType Si ScintCryst\n", + "struct Hit\n", + " arrivalTime::Float64\n", + " depositedEnergy::Float64\n", + " hittype::HitType\n", + " position::G4ThreeVector\n", + " function Hit(time, pos, edep, typ) # constructor\n", + " hit = new(time, edep, typ, G4ThreeVector())\n", + " assign(hit.position, pos) # this is needed to fill a G4ThreeVector (in C++ = operator)\n", + " return hit\n", + " end\n", + "end\n", + "\n", + "#--SD data structure---------------------------------------------------------------\n", + "struct CrystalSDData <: G4JLSDData\n", + " hitcollection::Vector{Hit} # define a hit collection\n", + " CrystalSDData() = new(Hit[]) # default constructor\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Sensitive Detector Functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Initialize method------------------------------------------------------------------------------\n", + "function crystal_initialize(::G4HCofThisEvent, data::CrystalSDData)::Nothing\n", + " empty!(data.hitcollection) # empty the hit collection at every event\n", + " return\n", + "end\n", + "#---Process Hit method-----------------------------------------------------------------------------\n", + "function crystal_processHits(step::G4Step, ::G4TouchableHistory, data::CrystalSDData)::Bool\n", + " part = step |> GetTrack |> GetParticleDefinition\n", + " edep = step |> GetTotalEnergyDeposit\n", + " edep < 0. && return false\n", + " pos = step |> GetPostStepPoint |> GetPosition\n", + " push!(data.hitcollection, Hit(0., pos, edep, ScintCryst)) # fill the collection with a new Hit\n", + " return true\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## SD Instance\n", + "And finally construct the SD instance" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Create SD instance-----------------------------------------------------------------------------\n", + "crystal_sd = G4JLSensitiveDetector(\"Crystal_SD\", CrystalSDData(); # name an associated data are mandatory\n", + " processhits_method=crystal_processHits, # process hits method (also mandatory)\n", + " initialize_method=crystal_initialize) # intialize method" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To test this created SD we will create a very simple detector geometry" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "struct SimpleDetector <: G4JLDetector ; end\n", + "\n", + "#---Materials----------------------------------------------------------------------------------\n", + "nist = G4NistManager!Instance()\n", + "m_air = FindOrBuildMaterial(nist, \"G4_AIR\")\n", + "m_bgo = FindOrBuildMaterial(nist, \"G4_BGO\")\n", + "\n", + "#---Volumes------------------------------------------------------------------------------------\n", + "worldS = G4Box(\"world\", 1m, 1m, 1m)\n", + "worldLV = G4LogicalVolume(worldS, m_air, \"World\")\n", + "worldPV = G4PVPlacement(nothing, G4ThreeVector(), worldLV, \"World\", nothing, false, 0, false)\n", + "\n", + "crystalS = G4Box(\"world\", 5cm, 5cm, 20cm)\n", + "crystalLV = G4LogicalVolume(crystalS, m_bgo, \"Crystal\")\n", + "crystalPV = G4PVPlacement(nothing, G4ThreeVector(), crystalLV, \"Crystal\", worldLV, false, 0, false)\n", + "\n", + "#---define getConstructor\n", + "Geant4.getConstructor(::SimpleDetector)::Function = (::SimpleDetector) -> worldPV" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Define Application\n", + "We will use the end-of-event user action to see how many hits we have produced. We could also at this moment to store the hits of some statistical data for each event. Be aware that the function is called from different threads. It is the responsibility of the developer to ensure that the function es thread safe. The function `getSDdata` ensures that you get a copy of the SD data that corresponds to the running thread. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---End Event Action-------------------------------------------------------------------------------\n", + "function endeventaction(evt::G4Event, app::G4JLApplication)\n", + " hits = getSDdata(app, \"Crystal_SD\").hitcollection\n", + " eventID = evt |> GetEventID\n", + " if eventID < 10 || eventID % 100 == 0\n", + " G4JL_println(\"Event: $eventID with $(length(hits)) hits stored in this event\")\n", + " end\n", + " return\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We use a simple particle gun" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "particlegun = G4JLGunGenerator(particle = \"e-\", \n", + " energy = 3GeV, \n", + " direction = G4ThreeVector(0,0,1), \n", + " position = G4ThreeVector(0,0,-1m))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Create the Application-------------------------------------------------------------------------\n", + "app = G4JLApplication( detector = SimpleDetector(), # detector with parameters\n", + " generator = particlegun, # primary particle generator\n", + " nthreads = 4, # number of threads (MT)\n", + " physics_type = FTFP_BERT, # what physics list to instantiate\n", + " endeventaction_method = endeventaction, # end event action\n", + " sdetectors = [\"Crystal\" => crystal_sd] # mapping of LVs to SDs (+ means multiple LVs with same name)\n", + " )\n", + "configure(app)\n", + "initialize(app)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "beamOn(app,10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], "source": [] } ], "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, "language_info": { - "name": "python" + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" } }, "nbformat": 4, diff --git a/tutorial/docs/09-scoring-meshes.ipynb b/tutorial/docs/09-scoring-meshes.ipynb index ab52473..a85aba6 100644 --- a/tutorial/docs/09-scoring-meshes.ipynb +++ b/tutorial/docs/09-scoring-meshes.ipynb @@ -4,24 +4,203 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Scoring Meshes" + "# Scoring Meshes\n", + "\n", + "The user can specify scoring meshes to obtain quantities on the defined grid. In Geant4 this is achieved using a set of UI commands. In this Julia interface this functionality has been encapsulated in a number of data structures. The function to create a scoring mesh is [`G4JLScoringMesh`](@ref) and receive as arguments the the type and dimensions of the mesh, the position, the rotation, the number of bins in each dimension, and the quantities to accumulate with eventually some filter conditions.\n" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, + "outputs": [], + "source": [ + "using Geant4\n", + "using Geant4.SystemOfUnits\n", + "using CairoMakie" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's define a simple detector (from the previous notebook on sensitive detectors)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "struct SimpleDetector <: G4JLDetector ; end\n", + "\n", + "#---Materials----------------------------------------------------------------------------------\n", + "nist = G4NistManager!Instance()\n", + "m_air = FindOrBuildMaterial(nist, \"G4_AIR\")\n", + "m_bgo = FindOrBuildMaterial(nist, \"G4_BGO\")\n", + "\n", + "#---Volumes------------------------------------------------------------------------------------\n", + "worldS = G4Box(\"world\", 1m, 1m, 1m)\n", + "worldLV = G4LogicalVolume(worldS, m_air, \"World\")\n", + "worldPV = G4PVPlacement(nothing, G4ThreeVector(), worldLV, \"World\", nothing, false, 0, false)\n", + "\n", + "crystalS = G4Box(\"world\", 5cm, 5cm, 10cm)\n", + "crystalLV = G4LogicalVolume(crystalS, m_bgo, \"Crystal\")\n", + "crystalPV = G4PVPlacement(nothing, G4ThreeVector(), crystalLV, \"Crystal\", worldLV, false, 0, false)\n", + "\n", + "#---define getConstructor\n", + "Geant4.getConstructor(::SimpleDetector)::Function = (::SimpleDetector) -> worldPV" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We create a scoring mesh of the same size as the crystal and with bins 20, 20 and 40 in z direction to collect energy deposit, number of steps produced by gammas, electrons and positrons." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sc1 = G4JLScoringMesh(\"boxMesh_1\",\n", + " BoxMesh(5cm,5cm,20cm),\n", + " bins = (20, 20, 40),\n", + " quantities = [ energyDeposit(\"eDep\")\n", + " nOfStep(\"nOfStepGamma\", filters=[ParticleFilter(\"gammafilter\", \"gamma\")])\n", + " nOfStep(\"nOfStepEMinus\", filters=[ParticleFilter(\"eMinusFilter\", \"e-\")])\n", + " nOfStep(\"nOfStepEPlus\", filters=[ParticleFilter(\"ePlusFilter\", \"e+\")])\n", + " ]\n", + " );" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's initialize now the application" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "particlegun = G4JLGunGenerator(particle = \"e-\", \n", + " energy = 3GeV, \n", + " direction = G4ThreeVector(0,0,1), \n", + " position = G4ThreeVector(0,0,-1m))\n", + "\n", + "app = G4JLApplication(;detector = SimpleDetector(), # detector with parameters\n", + " generator = particlegun, # primary particle generator\n", + " nthreads = 4, # number of threads (MT)\n", + " physics_type = FTFP_BERT, # what physics list to instantiate\n", + " scorers = [sc1] # list of scorers \n", + " );" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "configure(app)\n", + "initialize(app)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "run for 1000 events" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "beamOn(app,1000)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now the scoring mesh variable `sc1` holds the collected information" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], - "source": [] + "source": [ + "sum, sum2, entries = sc1.eDep;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Each of these variables (`sum`, `sum2` and `entries`) is a `20×20×40 Array` corresponding to the deposit energy. \n", + "Let's see a X-Y plane of it in the middle of the Z axis:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sum[:,:,20]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A better way to to plot a heatmap with Makie." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "img = heatmap(sum[:,:,20], color=:thermal, title=\"Edep (XY)\")\n", + "display(\"image/png\", img)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "img = heatmap(sum[:,10,:]', color=:thermal, title=\"Edep (XZ)\")\n", + "display(\"image/png\", img)" + ] } ], "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, "language_info": { - "name": "python" + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" } }, "nbformat": 4, diff --git a/tutorial/docs/10-histograms.ipynb b/tutorial/docs/10-histograms.ipynb index 77f0cd8..7551670 100644 --- a/tutorial/docs/10-histograms.ipynb +++ b/tutorial/docs/10-histograms.ipynb @@ -8,20 +8,28 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": { "vscode": { "languageId": "julia" } }, - "outputs": [], - "source": [] + "source": [ + "This notebook will be done later with the new version of the FHist package" + ] } ], "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, "language_info": { - "name": "python" + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" } }, "nbformat": 4, diff --git a/tutorial/docs/11-event-display.ipynb b/tutorial/docs/11-event-display.ipynb index b4ea167..6c39ff3 100644 --- a/tutorial/docs/11-event-display.ipynb +++ b/tutorial/docs/11-event-display.ipynb @@ -4,24 +4,146 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Event Display" + "# Event Display\n", + "\n", + "For visualization applications, the user can create an instance of `G4JLEventDisplay([settings file])` and give it to the constructor of `G4JLApplication` in the `evtdisplay` attribute.\n", + "\n", + "The constructor accepts a visualization settings file that will overwrite the default settings in the file [Geant4.jl/ext/G4Vis/settings.jl](https://github.com/JuliaHEP/Geant4.jl/blob/master/ext/G4Vis/settings.jl). The format of the settings is Julia `NamedTuple`. Here is an example:\n", + "```julia\n", + "(\n", + " display = (\n", + " show_axis = false,\n", + " ),\n", + " trajectories = (\n", + " color = :yellow,\n", + " ),\n", + ")\n", + "```\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's define a simple detector and particle gun as generator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Geant4\n", + "using Geant4.SystemOfUnits\n", + "using CairoMakie, Rotations, IGLWrap_jll # to force loading G4Vis extension" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "struct SimpleDetector <: G4JLDetector ; end\n", + "\n", + "#---Materials----------------------------------------------------------------------------------\n", + "nist = G4NistManager!Instance()\n", + "m_air = FindOrBuildMaterial(nist, \"G4_AIR\")\n", + "m_bgo = FindOrBuildMaterial(nist, \"G4_BGO\")\n", + "\n", + "#---Volumes------------------------------------------------------------------------------------\n", + "worldS = G4Box(\"world\", 10cm, 10cm, 20cm)\n", + "worldLV = G4LogicalVolume(worldS, m_air, \"World\")\n", + "worldPV = G4PVPlacement(nothing, G4ThreeVector(), worldLV, \"World\", nothing, false, 0, false)\n", + "\n", + "crystalS = G4Box(\"world\", 5cm, 5cm, 10cm)\n", + "crystalLV = G4LogicalVolume(crystalS, m_bgo, \"Crystal\")\n", + "crystalPV = G4PVPlacement(nothing, G4ThreeVector(), crystalLV, \"Crystal\", worldLV, false, 0, false)\n", + "\n", + "#---define getConstructor\n", + "Geant4.getConstructor(::SimpleDetector)::Function = (::SimpleDetector) -> worldPV" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "particlegun = G4JLGunGenerator(particle = \"e-\", \n", + " energy = 1GeV, \n", + " direction = G4ThreeVector(0,0,1), \n", + " position = G4ThreeVector(0,0,-10cm));" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now define an instance of `G4JLEventDisplay`. The visualization settings file is in same directory as this notebook. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Event Display----------------------------------------------------------------------------------\n", + "evtdisplay = G4JLEventDisplay(joinpath(@__DIR__, \"VisSettings.jl\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Create the Application-------------------------------------------------------------------------\n", + "app = G4JLApplication(detector = SimpleDetector(), # simple detector\n", + " generator = particlegun, # primary generator to instantiate\n", + " physics_type = FTFP_BERT, # what physics list to instantiate\n", + " evtdisplay = evtdisplay # event display\n", + " )\n", + "\n", + "configure(app)\n", + "initialize(app)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Event Display is anow ready display events-----------------------------------------------------\n", + "beamOn(app, 1)\n", + "display(\"image/png\", evtdisplay.figure)" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, "outputs": [], - "source": [] + "source": [ + "beamOn(app, 1)\n", + "display(\"image/png\", evtdisplay.figure)" + ] } ], "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, "language_info": { - "name": "python" + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" } }, "nbformat": 4, diff --git a/tutorial/docs/VisSettings.jl b/tutorial/docs/VisSettings.jl new file mode 100644 index 0000000..31b75d2 --- /dev/null +++ b/tutorial/docs/VisSettings.jl @@ -0,0 +1,9 @@ +( + display = ( + backgroundcolor = :black, # Display background color + show_axis = false, + ), + trajectories = ( + color = :yellow, + ), +) \ No newline at end of file diff --git a/tutorial/docs/_toc.yml b/tutorial/docs/_toc.yml index 8b2d0de..791f6fd 100644 --- a/tutorial/docs/_toc.yml +++ b/tutorial/docs/_toc.yml @@ -22,5 +22,7 @@ parts: - file: examples/HBC30/HBC30 - file: examples/WaterPhantom/WaterPhantom - file: examples/Scintillation/Scintillation + - file: examples/TestEm3/TestEm3 + \ No newline at end of file diff --git a/tutorial/docs/examples/TestEm3/DetectorTestEm3.jl b/tutorial/docs/examples/TestEm3/DetectorTestEm3.jl new file mode 100644 index 0000000..85ede11 --- /dev/null +++ b/tutorial/docs/examples/TestEm3/DetectorTestEm3.jl @@ -0,0 +1,260 @@ +using Geant4 +using Geant4.SystemOfUnits + +#-------------------------------------------------------------------------------------------------- +# TestEm3 Detector definition +#-------------------------------------------------------------------------------------------------- + +mutable struct TestEm3Detector <: G4JLDetector + # main input parameters + const fNbOfAbsor::Int64 + const fNbOfLayers::Int64 + const fCalorSizeYZ::Float64 + const fAbsorThickness::Vector{Float64} + + # mutable (computed) detector data + fLayerThickness::Float64 + fCalorThickness::Float64 + fWorldSizeYZ::Float64 + fWorldSizeX::Float64 + + fAbsorMaterial::Vector{CxxPtr{G4Material}} + fWorldMaterial::CxxPtr{G4Material} + + fSolidWorld::G4Box + fLogicWorld::G4LogicalVolume + fPhysiWorld::CxxPtr{G4VPhysicalVolume} + + fSolidCalor::G4Box + fLogicCalor::G4LogicalVolume + fPhysiCalor::CxxPtr{G4VPhysicalVolume} + + fSolidLayer::G4Box + fLogicLayer::G4LogicalVolume + fPhysiLayer::CxxPtr{G4VPhysicalVolume} + + fSolidAbsor::Vector{G4Box} + fLogicAbsor::Vector{G4LogicalVolume} + fPhysiAbsor::Vector{CxxPtr{G4VPhysicalVolume}} + + function TestEm3Detector(;nbOfLayers = 50, + calorSizeYZ = 40cm, + absorThickness = [2.3mm, 5.7mm], + absorMaterial = ["G4_Pb", "G4_lAr"]) + self = new(length(absorThickness), nbOfLayers, calorSizeYZ, absorThickness) + + # Compute derived parameters of the calorimeter + self.fLayerThickness = 0. + for i in 1:self.fNbOfAbsor + self.fLayerThickness += self.fAbsorThickness[i] + end + self.fCalorThickness = self.fNbOfLayers * self.fLayerThickness + self.fWorldSizeX = 1.2 * self.fCalorThickness + self.fWorldSizeYZ = 1.2 * self.fCalorSizeYZ + + DefineMaterials!(self) + + nist = G4NistManager!Instance() + self.fWorldMaterial = FindOrBuildMaterial(nist,"Galactic") + self.fAbsorMaterial = [FindOrBuildMaterial(nist, mat) for mat in absorMaterial] + self + end + +end + +using Geant4.SystemOfUnits:cm3, mole, kelvin, atmosphere, eV, pascal, bar +const perCent=0.01 + +function DefineMaterials!(det::TestEm3Detector) + # to avoid warnings of redefinitions + G4Material!GetMaterial("liquidH2", false) != C_NULL && return + + nist = G4NistManager!Instance() + H = FindOrBuildElement(nist,1) + C = FindOrBuildElement(nist,6) + N = FindOrBuildElement(nist,7) + O = FindOrBuildElement(nist,8) + Si = FindOrBuildElement(nist,14) + Ge = FindOrBuildElement(nist,32) + Sb = FindOrBuildElement(nist,51) + I = FindOrBuildElement(nist,53) + Cs = FindOrBuildElement(nist,55) + Pb = FindOrBuildElement(nist,82) + Bi = FindOrBuildElement(nist,83) + + U5 = G4Isotope("U235", z=92, n=235, a=235.01*g/mole) + U8 = G4Isotope("U238", z=92, n=238, a=238.03*g/mole) + + U = G4Element("enriched Uranium", "U", 2) + AddIsotope(U, move!(U5), 0.90) + AddIsotope(U, move!(U8), 0.10) + + G4Material("liquidH2", z=1., a= 1.008*g/mole, density= 70.8*mg/cm3) + G4Material("Aluminium", z=13., a= 26.98*g/mole, density= 2.700*g/cm3) + G4Material("Titanium", z=22., a= 47.867*g/mole, density= 4.54*g/cm3) + G4Material("Iron", z=26., a= 55.85*g/mole, density= 7.870*g/cm3) + G4Material("Copper", z=29., a= 63.55*g/mole, density= 8.960*g/cm3) + G4Material("Tungsten", z=74., a= 183.85*g/mole, density= 19.30*g/cm3) + G4Material("Gold", z=79., a= 196.97*g/mole, density= 19.32*g/cm3) + G4Material("Uranium", z=92., a= 238.03*g/mole, density= 18.95*g/cm3) + + # define a material from elements. case 1: chemical molecule + + H₂O = G4Material("Water", density=1.000*g/cm3, ncomponents=2) + + AddElement(H₂O, H, natoms=2) + AddElement(H₂O, O, natoms=1) + #SetMeanExcitationEnergy(GetIonisation(H₂O), 78.0*eV) + SetChemicalFormula(H₂O, "H_2O") + + CH = G4Material("Polystyrene", density= 1.032*g/cm3, ncomponents=2) + AddElement(CH, C, natoms=1) + AddElement(CH, H, natoms=1) + + Sci = G4Material("Scintillator", density= 1.032*g/cm3, ncomponents=2) + AddElement(Sci, C, natoms=9) + AddElement(Sci, H, natoms=10) + #SetBirksConstant(GetIonisation(Sci), 0.126*mm/MeV) + + Lct = G4Material("Lucite", density= 1.185*g/cm3, ncomponents=3) + AddElement(Lct, C, fractionmass=.5997) + AddElement(Lct, H, fractionmass=.0807) + AddElement(Lct, O, fractionmass=.3196) + + Sili = G4Material("Silicon", density= 2.330*g/cm3, ncomponents=1) + AddElement(Sili, Si, natoms=1) + + SiO₂ = G4Material("quartz", density= 2.200*g/cm3, ncomponents=2) + AddElement(SiO₂, Si, natoms=1) + AddElement(SiO₂, O , natoms=2) + + G10 = G4Material("NemaG10", density= 1.700*g/cm3, ncomponents=4) + AddElement(G10, Si, natoms=1) + AddElement(G10, O , natoms=2) + AddElement(G10, C , natoms=3) + AddElement(G10, H , natoms=3) + + CsI = G4Material("CsI", density= 4.534*g/cm3, ncomponents=2) + AddElement(CsI, Cs, natoms=1) + AddElement(CsI, I , natoms=1) + #SetMeanExcitationEnergy(GetIonisation(CsI), 553.1*eV) + + BGO = G4Material("BGO", density= 7.10*g/cm3, ncomponents=3) + AddElement(BGO, O , natoms=12) + AddElement(BGO, Ge, natoms= 3) + AddElement(BGO, Bi, natoms= 4) + + SiNx= G4Material("SiNx", density= 3.1 *g/cm3, ncomponents=3) + AddElement(SiNx, Si, natoms=300) + AddElement(SiNx, N, natoms=310) + AddElement(SiNx, H, natoms=6) + + # define gaseous materials using G4 NIST database + Air = FindOrBuildMaterial(nist, "G4_AIR") + ConstructNewGasMaterial(nist, "Air20","G4_AIR", 293kelvin, 1atmosphere) + + lAr = FindOrBuildMaterial(nist, "G4_lAr") + lArEm3 = G4Material("liquidArgon", density= 1.390g/cm3, ncomponents=1) + AddMaterial(lArEm3, lAr, fractionmass=1.0) + + # define a material from elements and others materials (mixture of mixtures) + Lead = G4Material("Lead", density=11.35*g/cm3, ncomponents=1) + AddElement(Lead, Pb, fractionmass=1.0) + + LeadSb = G4Material("LeadSb", density=11.35*g/cm3, ncomponents=2) + AddElement(LeadSb, Sb, fractionmass=4perCent) + AddElement(LeadSb, Pb, fractionmass=96perCent) + + Aerog = G4Material("Aerogel", density= 0.200*g/cm3, ncomponents=3) + AddMaterial(Aerog, SiO₂, fractionmass=62.5*perCent) + AddMaterial(Aerog, H₂O , fractionmass=37.4*perCent) + AddElement(Aerog, C , fractionmass= 0.1*perCent) + + # examples of gas in non STP conditions + CO₂ = G4Material("CarbonicGas", density= 27*mg/cm3, ncomponents=2, state=kStateGas, temperature= 325*kelvin, pressure= 50*atmosphere) + AddElement(CO₂, C, natoms=1) + AddElement(CO₂, O, natoms=2) + + steam = G4Material("WaterSteam", density= 1.0*mg/cm3, ncomponents=1, state=kStateGas, temperature= 273*kelvin, pressure= 1*atmosphere) + AddMaterial(steam, H₂O, fractionmass=1.) + + G4Material("ArgonGas", z=18., a=39.948*g/mole, density= 1.782*mg/cm3, state=kStateGas, temperature=273.15*kelvin, pressure=1*atmosphere) + + # examples of vacuum + universe_mean_density = 1.e-25*g/cm3 + move!(G4Material("Galactic", z=1., a=1.008*g/mole, density=universe_mean_density, state=kStateGas, temperature=2.73*kelvin, pressure=3.e-18*pascal)) + + beam = G4Material("Beam", density=1.e-5*g/cm3, ncomponents=1, state=kStateGas, temperature=273.15*kelvin, pressure=2e-2*bar) + AddMaterial(beam, Air, fractionmass=1.) + +end + +function TestEm3Construct(det::TestEm3Detector)::CxxPtr{G4VPhysicalVolume} + (; fNbOfAbsor, fNbOfLayers, fCalorSizeYZ, fAbsorThickness, fLayerThickness, fCalorThickness, fWorldSizeYZ, fWorldSizeX, fWorldMaterial, fAbsorMaterial) = det + + #isdefined(det,:fSolidWorld) && return det.fPhysiWorld + println("Building Geometry now!!!") + + #---World-------------------------------------------------------------------------------------- + det.fSolidWorld = G4Box("World", fWorldSizeX/2,fWorldSizeYZ/2,fWorldSizeYZ/2) + det.fLogicWorld = G4LogicalVolume(CxxPtr(det.fSolidWorld), fWorldMaterial, "World") + det.fPhysiWorld = G4PVPlacement(nothing, # no rotation + G4ThreeVector(), # at (0,0,0) + det.fLogicWorld, # its fLogical volume + "World", # its name + nothing, # its mother volume + false, # no boolean operation + 0) #copy number + + #---Calorimeter-------------------------------------------------------------------------------- + det.fSolidCalor = G4Box("Calorimeter", fCalorThickness/2,fCalorSizeYZ/2,fCalorSizeYZ/2) + det.fLogicCalor = G4LogicalVolume(CxxPtr(det.fSolidCalor), fWorldMaterial, "Calorimeter") + det.fPhysiCalor = G4PVPlacement(nothing, # no rotation + G4ThreeVector(), # at (0,0,0) + det.fLogicCalor, # its fLogical volume + "Calorimeter", # its name + det.fLogicWorld, # its mother volume + false, # no boolean operation + 0) # copy number + + #---Layers------------------------------------------------------------------------------------- + det.fSolidLayer = G4Box("Layer", fLayerThickness/2,fCalorSizeYZ/2,fCalorSizeYZ/2) + det.fLogicLayer = G4LogicalVolume(CxxPtr(det.fSolidLayer), fWorldMaterial, "Layer") + if fNbOfLayers > 1 + det.fPhysiLayer = G4PVReplica("Layer", + det.fLogicLayer, + det.fLogicCalor, + kXAxis, + fNbOfLayers, + fLayerThickness) + else + det.fPhysiLayer = G4PVPlacement(nothing, + G4ThreeVector(), + det.fLogicLayer, + "Layer", + det.fLogicCalor, + false, + 0) + end + + #---Absorbers---------------------------------------------------------------------------------- + det.fSolidAbsor = [G4Box("Absorber", fAbsorThickness[k]/2,fCalorSizeYZ/2,fCalorSizeYZ/2) for k in 1:fNbOfAbsor] + det.fLogicAbsor = [G4LogicalVolume(CxxPtr(det.fSolidAbsor[k]), fAbsorMaterial[k], GetName(fAbsorMaterial[k])) for k in 1:fNbOfAbsor] + + xcenter = zeros(fNbOfAbsor) + xfront = -0.5*fLayerThickness + for k in 1:fNbOfAbsor + xcenter[k] = xfront + 0.5*fAbsorThickness[k] + xfront += fAbsorThickness[k] + end + det.fPhysiAbsor = [G4PVPlacement(nothing, + G4ThreeVector(xcenter[k],0.,0.), + det.fLogicAbsor[k], + String(GetName(fAbsorMaterial[k])), + det.fLogicLayer, + false, + k) for k in 1:fNbOfAbsor] + return det.fPhysiWorld +end + +Geant4.getConstructor(::TestEm3Detector)::Function = TestEm3Construct diff --git a/tutorial/docs/examples/TestEm3/TestEm3.ipynb b/tutorial/docs/examples/TestEm3/TestEm3.ipynb new file mode 100644 index 0000000..aadbd66 --- /dev/null +++ b/tutorial/docs/examples/TestEm3/TestEm3.ipynb @@ -0,0 +1,439 @@ +{ + "cells": [ + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# TestEM3 Example\n", + "\n", + " - How to collect energy deposition in a sampling calorimeter.\n", + " - How to survey energy flow.\n", + " - How to print stopping power.\n", + " \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Geant4\n", + "using Geant4.SystemOfUnits\n", + "using FHist, Printf, Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "## The Geometry\n", + " \n", + " The calorimeter is a box made of a given number of layers.\n", + " A layer consists of a sequence of various absorbers (maximum MaxAbsor=9).\n", + " The layer is replicated.\n", + " \n", + " Parameters defining the calorimeter :\n", + " - the number of layers,\n", + " - the number of absorbers within a layer,\t\t\n", + " - the material of the absorbers,\n", + " - the thickness of the absorbers,\n", + " - the transverse size of the calorimeter (the input face is a square). \n", + " \n", + " In addition a transverse uniform magnetic field can be applied.\n", + " \n", + " The default geometry is constructed in DetectorConstruction class, but all\n", + " of the above parameters can be modified interactively via the commands \n", + " defined in the DetectorMessenger class.\n", + " \n", + "```\n", + " |<----layer 0---------->|<----layer 1---------->|<----layer 2---------->|\n", + " | | | | |\n", + " ==========================================================================\n", + " || | || | || | ||\n", + " || | || | || | ||\n", + " || abs 1 | abs 2 || abs 1 | abs 2 || abs 1 | abs 2 ||\n", + " || | || | || | ||\n", + " || | || | || | ||\n", + " beam || | || | || | ||\n", + "======> || | || | || | ||\n", + " || | || | || | ||\n", + " || | || | || | ||\n", + " || | || | || | ||\n", + " || | || | || | ||\n", + " || cell 1 | cell 2 || cell 3 | cell 4 || cell 5 | cell 6 ||\n", + " ==========================================================================\n", + " ^ ^ ^ ^ ^ ^ ^\n", + " pln1 pln2 pln3 pln4 pln5 pln6 pln7\n", + " ```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Define Detector Parameters struct--------------------------------------------------------------\n", + "include(joinpath(@__DIR__, \"DetectorTestEm3.jl\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulation Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Define Simulation Data struct------------------------------------------------------------------\n", + "mutable struct TestEm3SimData <: G4JLSimulationData\n", + " #---Run data-----------------------------------------------------------------------------------\n", + " fParticle::CxxPtr{G4ParticleDefinition}\n", + " fEkin::Float64\n", + " \n", + " fChargedStep::Int32\n", + " fNeutralStep::Int32\n", + "\n", + " fN_gamma::Int32\n", + " fN_elec::Int32\n", + " fN_pos::Int32\n", + "\n", + " fEnergyDeposit::Vector{Float64} # Energy deposit per event\n", + " fTrackLengthCh::Vector{Float64} # Track length per event\n", + "\n", + " fEdepEventHistos::Vector{Hist1D}\n", + " fTrackLengthChHistos::Vector{Hist1D}\n", + " fEdepHistos::Vector{Hist1D}\n", + " fAbsorLabel::Vector{String}\n", + "\n", + " TestEm3SimData() = new()\n", + "end\n", + "#---add function-----------------------------------------------------------------------------------\n", + "function add!(x::TestEm3SimData, y::TestEm3SimData)\n", + " x.fChargedStep += y.fChargedStep\n", + " x.fNeutralStep += y.fNeutralStep\n", + " x.fN_gamma += y.fN_gamma\n", + " x.fN_elec += y.fN_elec\n", + " x.fN_pos += y.fN_pos\n", + " x.fEdepEventHistos += y.fEdepEventHistos\n", + " x.fTrackLengthChHistos += y.fTrackLengthChHistos\n", + " x.fEdepHistos += y.fEdepHistos\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Plot the Sumulation data-----------------------------------------------------------------------\n", + "function do_plot(data::TestEm3SimData)\n", + " (;fEdepHistos, fEdepEventHistos, fTrackLengthChHistos, fAbsorLabel) = data\n", + " lay = @layout [°; ° °]\n", + " p = plot(layout=lay, show=true, size=(1000,800))\n", + " for (h, l) in zip(fEdepHistos, fAbsorLabel)\n", + " plot!(subplot=1, h, title=\"Energy Deposition\", xlabel=\"layer #\", label=l, show=true)\n", + " end\n", + " for (h, l) in zip(fEdepEventHistos, fAbsorLabel)\n", + " plot!(subplot=2, h, title=\"Energy/event Distribution\", label=l, xlabel=\"MeV\")\n", + " end\n", + " for (h, l) in zip(fTrackLengthChHistos, fAbsorLabel)\n", + " plot!(subplot=3, h, title=\"Track Lengh Distribution\", label=l, xlabel=\"mm\")\n", + " end\n", + " display(\"image/png\", p)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Particle Gun initialization--------------------------------------------------------------------\n", + "particlegun = G4JLGunGenerator(particle = \"e-\", \n", + " energy = 1GeV, \n", + " direction = G4ThreeVector(1,0,0), \n", + " position = G4ThreeVector(0,0,0)) # temporary potition, will update once the detector is constructed" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## User Actions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Step action------------------------------------------------------------------------------------\n", + "function stepaction(step::G4Step, app::G4JLApplication)::Nothing\n", + " data = getSIMdata(app)\n", + " detector = app.detector\n", + " prepoint = GetPreStepPoint(step)\n", + " endPoint = GetPostStepPoint(step)\n", + " track = GetTrack(step)\n", + " \n", + " # Return if step in not in the world volume\n", + " prepoint |> GetPhysicalVolume |> GetLogicalVolume |> GetMaterial == detector.fWorldMaterial && return nothing\n", + " \n", + " particle = GetDefinition(track)\n", + " charge = GetPDGCharge(particle)\n", + " stepl = 0.\n", + " if charge != 0.\n", + " stepl = GetStepLength(step)\n", + " data.fChargedStep += 1\n", + " else\n", + " data.fNeutralStep += 1\n", + " end\n", + " edep = GetTotalEnergyDeposit(step) * GetWeight(track)\n", + " absorNum = GetCopyNumber(GetTouchable(prepoint), 0)\n", + " layerNum = GetCopyNumber(GetTouchable(prepoint), 1) + 1 # replicas copynumber starts at 0\n", + "\n", + " data.fEnergyDeposit[absorNum] += edep\n", + " data.fTrackLengthCh[absorNum] += stepl \t\n", + "\n", + " push!(data.fEdepHistos[absorNum], layerNum, edep)\n", + " nothing\n", + "end\n", + "\n", + "#---Tracking pre-action----------------------------------------------------------------------------\n", + "let G4Gamma, G4Electron, G4Positron, first=true\n", + "global function pretrackaction(track::G4Track, app::G4JLApplication)::Nothing\n", + " if first\n", + " G4Gamma = FindParticle(\"gamma\")\n", + " G4Electron = FindParticle(\"e-\")\n", + " G4Positron = FindParticle(\"e+\")\n", + " first = false\n", + " end\n", + " data = getSIMdata(app)\n", + " d = GetDefinition(track)\n", + " if d === G4Gamma \n", + " data.fN_gamma += 1\n", + " elseif d === G4Electron\n", + " data.fN_elec +=1\n", + " elseif d === G4Positron\n", + " data.fN_pos += 1\n", + " end\n", + " nothing\n", + "end\n", + "end\n", + "#---Tracking post-action---------------------------------------------------------------------------\n", + "function posttrackaction(track::G4Track, ::G4JLApplication)::Nothing\n", + " return\n", + "end\n", + "#---Begin Run Action-------------------------------------------------------------------------------\n", + "function beginrun(run::G4Run, app::G4JLApplication)::Nothing\n", + " data = getSIMdata(app)\n", + " (; fNbOfAbsor, fNbOfLayers, fAbsorMaterial, fAbsorThickness) = app.detector\n", + " gun = app.generator.data.gun\n", + " data.fParticle = GetParticleDefinition(gun)\n", + " data.fEkin = GetParticleEnergy(gun)\n", + " data.fN_gamma = data.fN_elec = data.fN_pos = 0\n", + " data.fChargedStep = data.fNeutralStep = 0\n", + " # init arrays\n", + " data.fEnergyDeposit = zeros(fNbOfAbsor)\n", + " data.fTrackLengthCh = zeros(fNbOfAbsor)\n", + " data.fEdepHistos = [Hist1D(Float64; bins=0.:1.0:fNbOfLayers) for i in 1:fNbOfAbsor]\n", + " data.fEdepEventHistos = [Hist1D(;bins=0:10:1000) for i in 1:fNbOfAbsor]\n", + " data.fTrackLengthChHistos = [Hist1D(;bins=0:20:2000) for i in 1:fNbOfAbsor]\n", + " data.fAbsorLabel = [\"$(fAbsorThickness[i])mm of $(fAbsorMaterial[i] |> GetName |> String)\" for i in 1:fNbOfAbsor]\n", + " return\n", + "end\n", + "#---End Run Action---------------------------------------------------------------------------------\n", + "function endrun(run::G4Run, app::G4JLApplication)::Nothing\n", + " #---end run action is called for each workwer thread and the master onc\n", + " if G4Threading!G4GetThreadId() == -1 \n", + " data = app.simdata[1]\n", + " #---This is the master thread, so we need to add all the simuation results-----------------\n", + " for d in app.simdata[2:end]\n", + " add!(data, d)\n", + " end\n", + " nEvt = GetNumberOfEvent(run)\n", + " norm = nEvt > 0 ? 1/nEvt : 1.\n", + "\n", + " @printf \"------------------------------------------------------------\\n\"\n", + " @printf \" Beam particle %s E = %.2f GeV\\n\" String(GetParticleName(data.fParticle)) data.fEkin/GeV \n", + " @printf \" Mean number of gamma %.2f\\n\" data.fN_gamma*norm \n", + " @printf \" Mean number of e- %.2f\\n\" data.fN_elec*norm \n", + " @printf \" Mean number of e+ %.2f\\n\" data.fN_pos*norm \n", + " @printf \" Mean number of charged steps %f\\n\" data.fChargedStep*norm \n", + " @printf \" Mean number of neutral steps %f\\n\" data.fNeutralStep*norm \n", + " @printf \"------------------------------------------------------------\\n\"\n", + " else\n", + " G4JL_println(\"end-run for worker $(G4Threading!G4GetThreadId())\") \n", + " end\n", + "end\n", + "#---Begin Event Action-----------------------------------------------------------------------------\n", + "function beginevent(evt::G4Event, app::G4JLApplication)\n", + " data = getSIMdata(app)\n", + " fill!(data.fEnergyDeposit, 0.0)\n", + " fill!(data.fTrackLengthCh, 0.0)\n", + " return\n", + "end\n", + "#---End Event Action-------------------------------------------------------------------------------\n", + "function endevent(evt::G4Event, app::G4JLApplication)\n", + " data = getSIMdata(app)\n", + " (; fNbOfAbsor, fNbOfLayers) = app.detector\n", + " for i in 1:fNbOfAbsor\n", + " push!(data.fEdepEventHistos[i], data.fEnergyDeposit[i])\n", + " push!(data.fTrackLengthChHistos[i], data.fTrackLengthCh[i])\n", + " end\n", + " return\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Create the Application-------------------------------------------------------------------------\n", + "app = G4JLApplication(detector = TestEm3Detector(), # detector with parameters\n", + " simdata = TestEm3SimData(), # simulation data structure\n", + " generator = particlegun, # primary particle generator\n", + " nthreads = 8, # number of threads (MT)\n", + " physics_type = FTFP_BERT, # what physics list to instantiate\n", + " #----Actions--------------------------------\n", + " stepaction_method = stepaction, # step action method\n", + " pretrackaction_method = pretrackaction, # pre-tracking action\n", + " posttrackaction_method = posttrackaction, # post-tracking action\n", + " beginrunaction_method = beginrun, # begin-run action (initialize counters and histograms)\n", + " endrunaction_method = endrun, # end-run action (print summary) \n", + " begineventaction_method = beginevent, # begin-event action (initialize per-event data)\n", + " endeventaction_method = endevent # end-event action (fill histogram per event data)\n", + " );\n", + " \n", + "\n", + "#ui\"/tracking/verbose 1\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#---Configure, Initialize and Run------------------------------------------------------------------ \n", + "configure(app)\n", + "initialize(app)\n", + "SetParticlePosition(particlegun, G4ThreeVector(-app.detector.fWorldSizeX/2,0,0)) # Only now is known the size of the 'world'" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Start the initial run" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "beamOn(app, 1000)\n", + "do_plot(app.simdata[1])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Change the particle gun energy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "SetParticleEnergy(particlegun, 100MeV)\n", + "beamOn(app, 100)\n", + "do_plot(app.simdata[1])" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Change the geometry and re-start the run" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "reinitialize(app, TestEm3Detector(absorThickness = [2.3mm, 5.7mm, 1mm],\n", + " absorMaterial = [\"G4_Pb\", \"G4_lAr\", \"G4_Al\"]))\n", + "beamOn(app, 100)\n", + "do_plot(app.simdata[1])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@time beamOn(app, 10000)" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Change the definition of the action and re-start" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function posttrackaction(track::G4Track, ::G4JLApplication)::Nothing\n", + " G4JL_println(\"Track ID: $(GetTrackID(track)) ended\")\n", + "end\n", + "beamOn(app, 3)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.10.0", + "language": "julia", + "name": "julia-1.10" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.10.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}