Scintillating Detector Example

Taken from the original at settwi/g4-basic-scintillation and adapted to the Geant4.jl interface

This is an example with optical photons and customized physics list, as well as with a couple of sensitive detectors (for the crystal and silicon) and some simple analysis of the results.

using Geant4
using Geant4.PhysicalConstants
using Geant4.SystemOfUnits
using FHist, Printf, Plots
#using CairoMakie, Rotations, IGLWrap_jll  # to force loading G4Vis extension

Detector Geometry

The detector construction and parameters are in separated files. You also visualize the detector at this moment.

include(joinpath(@__DIR__, "Parameters.jl"))
include(joinpath(@__DIR__, "Detector.jl"))

crysdet = ScintDetector()

#world = scintConstruct(crysdet)
#img = draw(world[])
#display("image/png", img)

ScintDetector(50.0, true, true, 20.0, 10.0, 2.0, 1.0, 1.0, 10.0, true, CxxPtr{G4Material}(Ptr{G4Material} @0x0000000000000000), CxxPtr{G4Material}(Ptr{G4Material} @0x0000000000000000), CxxPtr{G4Material}(Ptr{G4Material} @0x0000000000000000), CxxPtr{G4Material}(Ptr{G4Material} @0x0000000000000000), CxxPtr{G4Material}(Ptr{G4Material} @0x0000000000000000), CxxPtr{G4Material}(Ptr{G4Material} @0x0000000000000000), CxxPtr{G4Material}(Ptr{G4Material} @0x0000000000000000))

Physics

We construct a physics list starting from the pre-defined FTFP_BERT list, replacing the EM part with G4EmStandardPhysics_option4, and registering the optical photon physics G4OpticalPhysics. The constructor ScintPhysicsList() will be called by the toolkit at the adequate moment in the initialization of the application.

struct ScintPhysicsList <: G4VUserPhysicsList
    function ScintPhysicsList(verbose)
        pl = FTFP_BERT(verbose)
        ReplacePhysics(pl, move!(G4EmStandardPhysics_option4(verbose)))
        RegisterPhysics(pl, move!(G4OpticalPhysics(verbose)))
        # need to enable scintillation
        optpar = G4OpticalParameters!Instance()
        SetProcessActivation(optpar, "Scintillation", true);
        # I have found Cherenkov radiation to be error-prone
        SetProcessActivation(optpar, "Cerenkov", true);
        return pl
    end 
end

Particle Gun

The primary particle generator is a simple particle gun.

particlegun = G4JLGunGenerator(particle = "gamma", 
                               energy = 30keV, 
                               direction = G4ThreeVector(0, 0, -1), 
                               position = G4ThreeVector(0, 0, 2cm))

G4JLGunGenerator("ParticleGun", Geant4.G4JLParticleGunData(nothing, "gamma", G4ThreeVector(0.0,0.0,-1.0), G4ThreeVector(0.0,0.0,20.0), 0.03), Geant4.var"#init#19"(), Geant4.var"#gen#20"(), G4JLGeneratorAction[])

Simulation Data

It is normally filled by the user actions. We define a set of histograms using the FHist.jl package and set of counters. We need to provide a method add! to be able to reduce the output in case we run in multi-threading mode.

We also provide a method to plot the results.

#---Simulation Data-------------(normally filled by actions)---------------------------------------
const Hist1D64 = Hist1D{Float64, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}
mutable struct ScintSimData <: G4JLSimulationData
    scintPhotonsPerEvent::Int64
    scintPhotonsHisto::Hist1D64
    siHitsHisto::Hist1D64
    crysEdepHisto::Hist1D64
    ScintSimData() = new()
end
#---Addd function-----------------------------------------------------------------------------------
function add!(x::ScintSimData, y::ScintSimData)
    x.scintPhotonsHisto += y.scintPhotonsHisto
    x.siHitsHisto += y.siHitsHisto
    x.crysEdepHisto += y.crysEdepHisto
end
#---Plot Simulation data----------------------------------------------------------------------------
function do_plot(data::ScintSimData)
    (;scintPhotonsHisto, siHitsHisto, crysEdepHisto) = data
    lay = @layout [° °; °]
    plot(layout=lay, show=true, size=(1000,700))
    plot!(subplot=1, scintPhotonsHisto, title="# scintillating photons/event", xlabel="# photons", show=true)
    plot!(subplot=2, siHitsHisto, title="# hits in Silicon/event", xlabel="# hits")
    plot!(subplot=3, crysEdepHisto, title="energy deposited in crystal", xlabel="keV")
end

too many parameters for type Hist1D

Stacktrace:
 [1] top-level scope
   @ In[5]:2

User Actions and Application definition

include(joinpath(@__DIR__, "UserActions.jl"))

#--------------------------------------------------------------------------------------------------
#---Create the Application-------------------------------------------------------------------------
#--------------------------------------------------------------------------------------------------
app = G4JLApplication(; detector = crysdet,                           # detector with parameters
                        simdata = ScintSimData(),                     # simulation data structure
                        generator = particlegun,                      # primary particle generator
                        nthreads = 0,                                 # # of threads (0 = no MT)
                        physics_type = ScintPhysicsList,              # what physics list to instantiate
                        stepaction_method = stepping,                 # step action method
                        begineventaction_method = beginevent,         # begin-event action (initialize per-event data)
                        endeventaction_method = endevent,             # end-event action (fill histograms per event)
                        beginrunaction_method = beginrun,             # begin run action
                        endrunaction_method = endrun,                 # end run action
                        sdetectors = ["si_log" => silicon_SD,
                                      "cebr3_log" => crystal_SD]      # mapping of LVs to SDs (+ means multiple LVs with same name)
                      );

Configure, Initialize

configure(app)
initialize(app)

Run and Plot results

#ui`/tracking/verbose 1`
beamOn(app,1000)

do_plot(app.simdata[1])