Geant4 isotope source energy being single value of 1GeV

[JunwenDiao]

Hi
I was generating decay events using 207Bi and found that the energy depositions of all the events are single valued at 1 GeV：

@Yuan-Ru-Lin mentioned that the issue can be caused by the implementation of function SSDGenerator(source::SolidStateDetectors.IsotopeSource) in the Geant4 extension
the scripts and configs to reproduce the problem can be found in Archive.zip

[oschulz]

@Yuan-Ru-Lin do you have a fix for this?

[Yuan-Ru-Lin]

I came up with a workaround when I simulated for the HADES dataset. It's a workaround because I had to modify the generator of an app after I configure and initialize it, so it doesn't fit to the current API.

Also, the generator that I use is a G4JLGeneralParticleSource, not sure whether it would cause any difference.

function getapp(; nthreads = Threads.nthreads())::G4JLApplication
    app = G4JLApplication(
        # We have to change to the directory where the geometry is, so that the
        # relative paths in the GDML file are correct
        detector=cd(joinpath(@__DIR__, "../geometry/V08682A/ba_HS4/top_0r_78z/geometry")) do
            G4JLDetectorGDML(
                "./main_ba_HS4_top.gdml",
                validate_schema=false
            ) end,
        simdata=MySimData(),
        physics_type=LBE,
        # The generator is just a stub. We will be set it up later
        generator=G4JLGeneralParticleSource(),
        nthreads=nthreads,
        ## Actions below
        beginrunaction_method=beginrun,
        stepaction_method=stepaction,
        endeventaction_method=endevent,
        endrunaction_method=endrun,
        sdetectors=["Detector" => G4JLSensitiveDetector(
            "Detector_SD", MySDData(), processhits_method=_processHits, initialize_method=_initialize,
        )],
    )

    configure(app)
    initialize(app)

    # None of the sources will decay without this. It has to be executed after `initialize`.
    # See also https://geant4-forum.web.cern.ch/t/radioactive-decay-completely-not-working/11798/4
    ui`/process/had/rdm/thresholdForVeryLongDecayTime 1.0e+60 year`

    return app
end

function setupTh228HS2(app::G4JLApplication)
    source = app.generator.data.gps |> GetCurrentSource
    SetParticleDefinition(source, GetIon(90, 228, 0.0))
    SetMonoEnergy(source |> GetEneDist, 0.0MeV)
    ang = source |> GetAngDist
    SetAngDistType(ang, "iso")
    SetMinTheta(ang, 150.0deg)
    SetMaxTheta(ang, 180.0deg)
    pos = source |> GetPosDist
    SetPosDisType(pos, "Volume")
    SetPosDisShape(pos, "Cylinder")
    SetRadius(pos, 0.5)
    SetHalfZ(pos, 0.5)
    SetCentreCoords(pos, G4ThreeVector(0.0, 0.0, -47.0mm))
    ConfineSourceToVolume(pos, "Source_PV")
end

[Yuan-Ru-Lin]

Okay here's a WIP that generates nontrivial energy spectrum.

It doesn't look like a typical 207-Bi spectrum, though.

[oschulz]

Can we ask one of the Geant4 experts? @gipert?

[gipert]

I can have a look if you point me to the relevant part of source code

[Yuan-Ru-Lin]

@gipert Here is the relevant part. Note that the variable gun is of type CxxPtr{G4SingleParticleSource} instead of CxxPtr{G4ParticleGun}. I will rename it once the result of the physics test looks reasonable.

SolidStateDetectors.jl/ext/Geant4/g4jl_application.jl

Lines 100 to 138 in eeb5aa4

	function SSDGenerator(source::SolidStateDetectors.IsotopeSource; kwargs...)
	iszero(source.direction) && throw(ArgumentError("The direction of the source cannot be zero."))
	data = GeneratorData(;kwargs...)
	function _init(data::GeneratorData, ::Any)
	gun = data.gun = move!(G4SingleParticleSource())
	data.position = G4ThreeVector(
	source.position.x * Geant4.SystemOfUnits.meter,
	source.position.y * Geant4.SystemOfUnits.meter,
	source.position.z * Geant4.SystemOfUnits.meter
	)
	SetParticlePosition(gun, data.position)
	data.direction = G4ThreeVector(source.direction.x, source.direction.y, source.direction.z)
	SetParticleMomentumDirection(GetAngDist(gun), data.direction)
	SetMonoEnergy(gun |> GetEneDist, 0.0 * Geant4.SystemOfUnits.MeV)
	end
	function _gen(evt::G4Event, data::GeneratorData)::Nothing
	if isnothing(data.particle) # late initialize (after physics processes)
	data.particle = GetIon(source.Z, source.A, Float64(source.excitEnergy))
	SetParticleDefinition(data.gun, data.particle)
	SetParticleCharge(data.gun, source.ionCharge)
	end
	if !iszero(source.opening_angle)
	d::CartesianVector = normalize(source.direction)
	a::CartesianVector = normalize(d × (abs(d.x) == 1 ? CartesianVector(0,1,0) : CartesianVector(1,0,0)))
	b::CartesianVector = normalize(a × d)
	ϕ = rand()*2π
	θ = acos(1 - (1 - cos(source.opening_angle))*rand())
	v = (cos(θ) * d + sin(θ) * (cos(ϕ) * a + sin(ϕ) * b))
	direction = G4ThreeVector(v.x, v.y, v.z)
	SetParticleMomentumDirection(GetAngDist(data.gun), direction)
	end
	GeneratePrimaryVertex(data.gun, CxxPtr(evt))
	end
	G4JLPrimaryGenerator("SSDGenerator", data; init_method=_init, generate_method=_gen)
	end

[Yuan-Ru-Lin]

Okay the gamma spectrum I showed above does make sense. According to this data sheet, there should be three gamma lines (570 keV, 1064 keV and 1770 keV) and they are all there. I was looking at energy spectrum collected by other type of detectors.

[Yuan-Ru-Lin]

@oschulz I don't think there needs to be additional modification; feel free to merge the PR and close this issue if it looks good to you

[fhagemann]

You also the the X-ray lines at 74keV and higher.
However, I think the spectrum is subject to pile-up, almost every gamma line shows a pileup peak, 74keV right of its original energy. Might it be that we are summing energies that would be more than 10 microseconds apart, that would result in two distinct waveforms to a single energy value. This could also be why you saw the monoenergetic 1GeV before, if it would just sum up all energies given by a decay cascade, no matter how much in time they would be apart.

To my knowledge, we do not extract the timing information from the Geant4 simulation (yet?). So, someone might wanna update the code to write out the correct timing information to cross-check this.

[oschulz]

So, someone might wanna update the code to write out the correct timing information to cross-check this.

@Yuan-Ru-Lin could you take a look at that, before we merge this?

[Yuan-Ru-Lin]

@fhagemann is right. Pileup events¹ do contribute to those peaks to the right of the gamma lines.

Almost all the discontinuities² have negative time differences, which I suppose is because those are two G4Tracks evolving independently of each other.

Anyhow, I suggest we just store the time information and let users decide what to do with it.

Footnotes

1. I uncommented and made some slight modifications to the definition of thit — thit = (step |> GetPostStepPoint |> GetGlobalTime) / Geant4.SystemOfUnits.nanosecond * u"ns". Pileup events are defined to be those whose hits do not occur at exactly the same time. ↩

2. Some statistics — about 18% of the events are pileup events and a handful of them have more than one discontinuity in time.

   julia> sort(countmap(length.(findall.(!=(0.0u"ns"), Δt))))
   OrderedCollections.OrderedDict{Int64, Int64} with 9 entries:
       0 => 40972
       1 => 9008
       2 => 8
       3 => 2
       5 => 2
       6 => 2
       7 => 2
       8 => 2
       9 => 2
   

   ↩

[fhagemann]

Nice work @Yuan-Ru-Lin
I would suggest to merge the PR as is. We could at some point think about also to pass the timing information to SSD to be able to simulate hard pile-up events (and maybe even soft pile-up events), but let's not worry about this for this issue/PR.

[fhagemann]

Closed by #508