Sensitive Detectors

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.

Instantiating a G4JLSensitiveDetector will require the following arguments:

• name. A string to identify the sensitive detector. No default.

• sd data. A instance of a user defined G4JLSDData mutable data structure that will passed to each callback invocation.

• initialize method. User method that is called at the beginning of each event. The signature is (::B2aSDData, ::G4HCofThisEvent)::Nothing.

• endOfEvent method. User method that is called at the end of each event. The signature is (::B2aSDData, ::G4HCofThisEvent)::Nothing.

• 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 reference manual to see what you can get from the G4Step. It returns true if a true hit is generated.

using Geant4
using Geant4.SystemOfUnits

Sensitive Detector Data

Let’s define something very simple, for example collecting generated Hits (i.e deposit energies) on a sensitive detector.

#---Hit structure------------------------------------------------------------------
@enum HitType Si ScintCryst
struct Hit
    arrivalTime::Float64
    depositedEnergy::Float64
    hittype::HitType
    position::G4ThreeVector
    function Hit(time, pos, edep, typ)  # constructor
        hit = new(time, edep, typ, G4ThreeVector())
        assign(hit.position, pos)       # this is needed to fill a G4ThreeVector (in C++ = operator)
        return hit
    end
end

#--SD data structure---------------------------------------------------------------
struct CrystalSDData <: G4JLSDData
    hitcollection::Vector{Hit}      # define a hit collection
    CrystalSDData() = new(Hit[])    # default constructor
end

Sensitive Detector Functions

#---Initialize method------------------------------------------------------------------------------
function crystal_initialize(::G4HCofThisEvent, data::CrystalSDData)::Nothing
    empty!(data.hitcollection)                                # empty the hit collection at every event
    return
end
#---Process Hit method-----------------------------------------------------------------------------
function crystal_processHits(step::G4Step, ::G4TouchableHistory, data::CrystalSDData)::Bool
    part = step |> GetTrack |> GetParticleDefinition
    edep = step |> GetTotalEnergyDeposit
    edep <  0. && return false
    pos = step |> GetPostStepPoint |> GetPosition
    push!(data.hitcollection, Hit(0., pos, edep, ScintCryst))  # fill the collection with a new Hit
    return true
end

crystal_processHits (generic function with 1 method)

SD Instance

And finally construct the SD instance

#---Create SD instance-----------------------------------------------------------------------------
crystal_sd = G4JLSensitiveDetector("Crystal_SD", CrystalSDData();           # name an associated data are mandatory
             processhits_method=crystal_processHits, # process hits method (also mandatory)
             initialize_method=crystal_initialize)   # intialize method

Geant4.G4JLProtoSD{CrystalSDData}("Crystal_SD", CrystalSDData(Hit[]), Main.crystal_processHits, Main.crystal_initialize, nothing)

To test this created SD we will create a very simple detector geometry

struct SimpleDetector <: G4JLDetector ; end

#---Materials----------------------------------------------------------------------------------
nist = G4NistManager!Instance()
m_air = FindOrBuildMaterial(nist, "G4_AIR")
m_bgo = FindOrBuildMaterial(nist, "G4_BGO")

#---Volumes------------------------------------------------------------------------------------
worldS = G4Box("world", 1m, 1m, 1m)
worldLV = G4LogicalVolume(worldS, m_air, "World")
worldPV = G4PVPlacement(nothing, G4ThreeVector(), worldLV, "World", nothing, false, 0, false)

crystalS = G4Box("world", 5cm, 5cm, 20cm)
crystalLV = G4LogicalVolume(crystalS, m_bgo, "Crystal")
crystalPV = G4PVPlacement(nothing, G4ThreeVector(), crystalLV, "Crystal", worldLV, false, 0, false)

#---define getConstructor
Geant4.getConstructor(::SimpleDetector)::Function = (::SimpleDetector) -> worldPV

Define Application

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.

#---End Event Action-------------------------------------------------------------------------------
function endeventaction(evt::G4Event, app::G4JLApplication)
    hits = getSDdata(app, "Crystal_SD").hitcollection
    eventID = evt |> GetEventID
    if eventID < 10 || eventID % 100 == 0
      G4JL_println("Event: $eventID with $(length(hits)) hits stored in this event")
    end
    return
end

endeventaction (generic function with 1 method)

We use a simple particle gun

particlegun = G4JLGunGenerator(particle = "e-", 
                               energy = 3GeV, 
                               direction = G4ThreeVector(0,0,1), 
                               position = G4ThreeVector(0,0,-1m))

G4JLGunGenerator("ParticleGun", Geant4.G4JLParticleGunData(nothing, "e-", G4ThreeVector(0.0,0.0,1.0), G4ThreeVector(0.0,0.0,-1000.0), 3000.0), Geant4.var"#init#19"(), Geant4.var"#gen#20"(), G4JLGeneratorAction[])

#---Create the Application-------------------------------------------------------------------------
app = G4JLApplication( detector = SimpleDetector(),                  # detector with parameters
                       generator = particlegun,                      # primary particle generator
                       nthreads = 4,                                 # number of threads (MT)
                       physics_type = FTFP_BERT,                     # what physics list to instantiate
                       endeventaction_method = endeventaction,       # end event action
                       sdetectors = ["Crystal" => crystal_sd]        # mapping of LVs to SDs (+ means multiple LVs with same name)
                      )
configure(app)
initialize(app)

beamOn(app,10)