[Gate-users] Scatter dependent on range cut – bug?
David Leibold
D.Leibold at tudelft.nl
Tue May 17 23:42:08 CEST 2022
In case the attachment was automatically removed, here are the files:
MAIN.MAC:
# MATERIAL DATABASE
/gate/geometry/setMaterialDatabase ./GateMaterials.db
# WORLD
/gate/world/geometry/setXLength 800. cm
/gate/world/geometry/setYLength 800. cm
/gate/world/geometry/setZLength 800. cm
/gate/world/setMaterial Vacuum
/gate/world/vis/forceWireframe
#######################################################
# SYSTEM GEOMETRY
#######################################################
# CT SCANNER
/gate/world/daughters/name CTscanner
/gate/world/daughters/insert box
/gate/CTscanner/placement/setTranslation 0.0 0.0 -230.5 mm
/gate/CTscanner/geometry/setXLength 400 mm
/gate/CTscanner/geometry/setYLength 400 mm
/gate/CTscanner/geometry/setZLength 1.5 cm
/gate/CTscanner/setMaterial Vacuum
# CT SCANNER -> MODULE
/gate/CTscanner/daughters/name module
/gate/CTscanner/daughters/insert box
/gate/module/geometry/setXLength 400 mm
/gate/module/geometry/setYLength 400 mm
/gate/module/geometry/setZLength 1.5 cm
/gate/module/setMaterial Vacuum
# CT SCANNER -> MODULE -> CLUSTER
/gate/module/daughters/name cluster
/gate/module/daughters/insert box
/gate/cluster/geometry/setXLength 400 mm
/gate/cluster/geometry/setYLength 400 mm
/gate/cluster/geometry/setZLength 1.5 cm
/gate/cluster/setMaterial Vacuum
# CT SCANNER -> MODULE -> CLUSTER -> PIXEL
/gate/cluster/daughters/name pixel
/gate/cluster/daughters/insert box
/gate/pixel/geometry/setXLength 400 mm
/gate/pixel/geometry/setYLength 400 mm
/gate/pixel/geometry/setZLength 1.5 cm
/gate/pixel/setMaterial CsI
# ATTACHING THE LEVELS TO THE SYSTEM
/gate/systems/CTscanner/module/attach module
/gate/systems/CTscanner/cluster_0/attach cluster
/gate/systems/CTscanner/pixel_0/attach pixel
# Attach a sensitive detector:
/gate/pixel/attachCrystalSD
#######################################################
# PHASE SPACE ACTOR
#######################################################
# First define its volume:
/gate/world/daughters/name PhS_BeforePlane
/gate/world/daughters/insert box
/gate/PhS_BeforePlane/setMaterial Vacuum
/gate/PhS_BeforePlane/geometry/setXLength 400 mm
/gate/PhS_BeforePlane/geometry/setYLength 400 mm
/gate/PhS_BeforePlane/geometry/setZLength 1 nm
/gate/PhS_BeforePlane/placement/setTranslation 0 0 -222 mm
# Add the phase space actor to the volume:
/gate/actor/addActor PhaseSpaceActor PhS_b
/gate/actor/PhS_b/attachTo PhS_BeforePlane
/gate/actor/PhS_b/save ../output/detector.root
#######################################################
# PHANTOM GEOMETRY
#######################################################
/gate/world/daughters/name testVolume
/gate/world/daughters/insert box
/gate/testVolume/geometry/setXLength 200 mm
/gate/testVolume/geometry/setYLength 200 mm
/gate/testVolume/geometry/setZLength 200 mm
/gate/testVolume/setMaterial Water
/gate/testVolume/attachPhantomSD
#######################################################
# PHYSICS
#######################################################
/gate/physics/addPhysicsList emstandard_opt4
#/gate/physics/Gamma/SetCutInRegion pixel 10 mm
#/gate/physics/Electron/SetCutInRegion pixel 10 mm
#/gate/physics/Gamma/SetCutInRegion testVolume 0.001 mm
#/gate/physics/Electron/SetCutInRegion testVolume 0.001 mm
#######################################################
# INITIALISE SIMULATION
#######################################################
/gate/run/initialize
#######################################################
# DETECTOR MODEL
#######################################################
/gate/digitizer/Singles/insert adder
/gate/digitizer/Singles/insert readout
/gate/digitizer/Singles/readout/setDepth 3
#######################################################
# DEFINE SOURCE
#######################################################
/gate/source/addSource xraypoint
/gate/source/xraypoint/gps/particle gamma
/gate/source/xraypoint/gps/position 0. 0. 617 mm
/gate/source/xraypoint/gps/ene/type Mono
/gate/source/xraypoint/gps/ene/mono 120 keV
/gate/source/xraypoint/gps/pos/type Point
/gate/source/xraypoint/gps/ang/mintheta 0 deg
/gate/source/xraypoint/gps/ang/maxtheta 0 deg
/gate/source/xraypoint/gps/ang/type iso
/gate/source/list
/gate/application/setTotalNumberOfPrimaries 100000
#######################################################
# SPECIFY OUTPUT FORMAT
#######################################################
/gate/output/tree/enable
/gate/output/tree/addFileName ../output/0.root
# Hits:
/gate/output/tree/hits/enable
# Singles:
/gate/output/tree/addCollection Singles
#######################################################
# Random engine
#######################################################
/gate/random/setEngineName MersenneTwister
# Use a random seed:
/gate/random/setEngineSeed 987654321
#######################################################
# START ACQUISITION
#######################################################
# In case you want a parametrised macro:
/gate/application/setTimeSlice 1 s
/gate/application/setTimeStart 0 s
/gate/application/setTimeStop 1 s
/gate/application/startDAQ
GATEMATERIALS.DB:
[Elements]
Hydrogen: S= H ; Z= 1. ; A= 1.01 g/mole
Oxygen: S= O ; Z= 8. ; A= 16.00 g/mole
Iodine: S= I ; Z= 53. ; A= 126.90 g/mole
Cesium: S= Cs ; Z= 55. ; A= 132.905 g/mole
[Materials]
Vacuum: d=0.000001 mg/cm3 ; n=1
+el: name=Hydrogen ; n=1
Water: d=1.00 g/cm3; n=2 ; state=liquid
+el: name=Hydrogen ; n=2
+el: name=Oxygen; n=1
CsI: d=4.51 g/cm3 ; n=2 ; state=Solid
+el: name=Cesium ; n=1
+el: name=Iodine ; n=1
EVALUATION.PY
import matplotlib.pyplot as plt
import numpy as np
import pandas
pandas.set_option('display.max_columns', None)
import uproot
import os
##############################################################
# H I T S
##############################################################
# Load the hits registered by the crystalSD detector:
rootfile = uproot.open("./Simulation/output/0.hits.root")
crystalSD = rootfile['tree'].arrays(filter_name=["*"], library='pd')
# For each particle, extract only the first interaction
# with the crystalSD; this is done by checking when
# the eventID changes in the array.
# First, a boolean mask array is created...
temp = (crystalSD['eventID'][1:].to_numpy() - crystalSD['eventID'][:-1].to_numpy())>0
onlyFirstHitBoolean = np.zeros(temp.size + 1, dtype=bool)
onlyFirstHitBoolean[0] = 1
onlyFirstHitBoolean[1:] = temp
# ...and then applied:
crystalSD_onlyFirstHit = crystalSD.loc[onlyFirstHitBoolean,:]
assert crystalSD_onlyFirstHit.shape[0] == len(set(crystalSD['eventID']))
# Determine the number of photons hitting the detector:
numberOfIncidentPhotons = len(set(crystalSD['eventID']))
# Extract the number of events whose flags indicate that the
# photon has undergone Compton or Rayleigh interaction or which
# is a secondary:
numberOfScatterAndSecondaries = np.sum( (crystalSD_onlyFirstHit['nPhantomCompton']>0) | (crystalSD_onlyFirstHit['nPhantomRayleigh']>0) | (crystalSD_onlyFirstHit['parentID']>0))
# Extract the number of events which hit the detector off-centre:
numberOfHitsOffCentre = crystalSD_onlyFirstHit[ (np.abs(crystalSD_onlyFirstHit['posX']) > 1E-4) | (np.abs(crystalSD_onlyFirstHit['posY']) > 1E-4) ].shape[0]
##############################################################
# P H A S E S P A C E A C T O R
##############################################################
# Load the result of the phase space actor in front of the
# detector:
rootfile = uproot.open("./Simulation/output/detector.root")
phs = rootfile['PhaseSpace'].arrays(filter_name="*", library='pd')
# Extract only those photons that travel from the direction of
# the source into the direciton of the detector; this way,
# photons scattered in the detector and passing back through
# the phase space actor in front of it are omitted:
phs_posDir = phs[phs['dZ'] < 0]
# Count the number of events that pass through the phase space
# actor off-centre:
numberPhaseSpaceOPffCentre = np.sum( ((np.abs(phs_posDir['X'])>1E-4) | (np.abs(phs_posDir['Y'])>1E-4)) )
##############################################################
# P R I N T R E S U L T S
##############################################################
print("Number of photons incident on the detector: {:>6}".format(numberOfIncidentPhotons))
print("Number of scattered particles and secondaries as registered by Gate: {:>6}".format(numberOfScatterAndSecondaries))
print("Number of registered hits off-centre: {:>6}".format(numberOfHitsOffCentre))
print("Number of events registered by the phase space actor off-centre: {:>6}".format(numberPhaseSpaceOPffCentre))
##############################################################
# PLOT RESULTS
##############################################################
# The following sections hold the results obtained for
# setting either a range cut for the detector pixels only
# or for the water phantom only, and plots them.
#-----------------------#
# Range cut in detector
#-----------------------#
#Number of photons incident on the detector:
numberOfIncidentPhotonsArray_DetectorRC = [13709, 13709, 13709, 13599, 13578, 13620]
# Number of scattered particles and secondaries as registered by Gate:
numberOfScatterAndSecondariesArray_DetectorRC = [9733, 9733, 9733, 2398, 2217, 2174]
# Number of registered hits off-centre:
numberOfHitsOffCentreArray_DetectorRC = [9730, 9730, 9730, 9697, 9722, 9668]
# Number of events registered by the phase space actor off-centre:
numberPhaseSpaceOPffCentre_DetectorRC = [9788, 9788, 9788, 9750, 9773, 9733]
### Plot
fig, ax = plt.subplots()
x = np.arange(len(numberOfIncidentPhotonsArray_DetectorRC))
ax.plot(x, numberOfIncidentPhotonsArray_DetectorRC, marker='*', label='incident photons')
ax.plot(x, numberOfScatterAndSecondariesArray_DetectorRC, marker='o', label='Gate: scatter & secondaries')
ax.plot(x, numberOfHitsOffCentreArray_DetectorRC, marker='s', label='hits off-centre')
ax.plot(x, numberPhaseSpaceOPffCentre_DetectorRC, marker='d', label='off-centre in phase-space')
ax.set_xlabel("Range cut [mm]")
ax.set_ylabel("Number of events")
ax.legend(loc='upper left', bbox_to_anchor=(1.0, 0.5))
plt.xticks(x, rangeCutArray)
# ax.set_ylim(13000, 14000)
# plt.savefig("temp.png", dpi=300, bbox_inches='tight')
plt.show()
#-----------------------#
# Range cut in phantom
#-----------------------#
# Number of photons incident on the detector:
numberOfIncidentPhotonsArray_PhantomRC = [13709, 13709, 13709, 13550, 13503, 13568]
# Number of scattered particles and secondaries as registered by Gate:
numberOfScatterAndSecondariesArray_PhantomRC = [9733, 9733, 9733, 9689, 9426, 9482]
# Number of registered hits off-centre:
numberOfHitsOffCentreArray_PhantomRC = [9730, 9730, 9730, 9681, 9446, 9523]
# Number of events registered by the phase space actor off-centre:
numberPhaseSpaceOPffCentre_PhantomRC = [9788, 9788, 9788, 9742, 9504, 9595]
rangeCutArray = ['None', '10', '1', '0.1', '0.01', '0.001']
### Plot
fig, ax = plt.subplots()
x = np.arange(len(numberOfIncidentPhotonsArray_PhantomRC))
ax.plot(x, numberOfIncidentPhotonsArray_PhantomRC, marker='*', label='incident photons')
ax.plot(x, numberOfScatterAndSecondariesArray_PhantomRC, marker='o', label='Gate: scatter & secondaries')
ax.plot(x, numberOfHitsOffCentreArray_PhantomRC, marker='s', label='hits off-centre')
ax.plot(x, numberPhaseSpaceOPffCentre_PhantomRC, marker='d', label='off-centre in phase-space')
ax.set_xlabel("Range cut [mm]")
ax.set_ylabel("Number of events")
ax.legend(loc='upper left', bbox_to_anchor=(1.0, 0.5))
plt.xticks(x, rangeCutArray)
# ax.set_ylim(9000, 10000)
# plt.savefig("temp.png", dpi=300, bbox_inches='tight')
plt.show()
Kind regards,
David Leibold
On 17 May2022, at 23:34, David Leibold <D.Leibold at tudelft.nl<mailto:D.Leibold at tudelft.nl>> wrote:
*** E-MAIL SECURITY WARNING *** This e-mail contained one or several attachments which are not allowed as per the company's e-mail policies. Following files were removed: minimum_example.zip If you need these files, please arrange alternate means of transfer with the sender. *** END OF SECURITY WARNING ***
From: David Leibold <D.Leibold at tudelft.nl<mailto:D.Leibold at tudelft.nl>>
Subject: [Gate-users] Scatter dependent on range cut – bug?
Date: 17 May 2022 at 23:34:06 CEST
To: gate-users <gate-users at lists.opengatecollaboration.org<mailto:gate-users at lists.opengatecollaboration.org>>
Dear Gaters,
I am simulating a cone beam CT setup where I output the Hits / Singles registered in the detector, which I then split into primary and scatter contribution based on the nPhantomCompton, nPhantomRayleigh and parentID parameters. These indicate the number of Compton or Rayleigh scatter events in the phantom and whether a particle was created by the source or via creation of secondary particles.
We noticed that the scatter-to-primary ratio in our simulation was far below the values reported in literature (more than a factor 10). I was able to show that for certain range cuts not every scattered photon exhibited a nPhantomCompton or nPhantomRayleigh value that would have indicated that it was indeed scattered. In the following I will try to explain what I mean by that and which behaviour of the simulation I can’t explain.
In the following figure you will see on the left the simulation setup:
* A monoenergetic source emitting a pencil beam of 120 keV photons,
* A water cube serving as a phantom,
* A phase space actor in front of a realistic (two-dimensional) detector,
* And a CTdetector of a realistic material (CsI), to which a crystalSD detector is attached.
With the output of the crystalSD detector (Hits/Singles) and the phase space actor, one can simply plot the coordinates of the intersection point of the incident photons with the detector plane. Any photon that does not go through the centre axis must then either have been scattered in the water phantom or be a secondary. (*)
Additionally, the crystalSD detector outputs the nPhantomCompton, nPhantomRayleigh and parentID parameters, which I used initially to separate primaries from scatter.
<setup.png>
(*) Note: In the Hits dataset, I only evaluate the first entry for a given eventID, so scatter inside the detector does not influence this evaluation (I guess?). In the data of the phase space actor, only photons flying in the direction of the detector are evaluated.
In the following the range cut inside the detector is now varied. I then evaluated:
* How many photons end up in the detector, which is evaluated by counting the number of unique eventIDs in the Hits dataset (label “incident photons”),
* How many photons are flagged by Gate as scatter or as secondary particle, based on the nPhantomCompton, nPhantomRayleigh and parentID parameters (label “Gate: scatter&secondaries”),
* How many photons are outside the centre axis, based on the their intersection with the phase space actor (label “off-centre in phase-space”),
* and how many photons are outside the centre axis, based on their intersection with the realistic detector (i.e., using their coordinates as registered in the Hits output, label “hits off-centre").
The following plot shows the results:
<stp_d.png>
As one can see, below a certain range cut the number of events flagged by Gate as scatter or secondaries (via the nPhantomCompton, nPhantomRayleigh and parentID parameters) drops considerably and deviates from the evaluation based on the number of photons that are off-centre. I have no explanation for this phenomenon.
Now let’s zoom into the the number of incident photons:
<stp_d_ni.png>
Interestingly, the number of photons incident on the detector varies with the range cut, which to me is unexpected. Sure, the number of interactions inside the detector increases with decreasing range cut, but this should have no effect on the number of incident photons.
Here is a zoom into the number of off-axis photons:
<stp_d_oa.png>
I don’t have an explanation why the phase space actor and the crystalSD have a slightly different number of off-axis photons, but this is at the moment not my main concern. Again, one can see that their number changes with decreasing range cut.
Last but not least, I kept the range cut in the detector at the default value and instead changed the range cut inside the phantom. Here is the result:
<stp_p.png>
In this case, the three different ways to extract the number of scatter and primaries agree more or less. Again, one can clearly see a change in the number of total incident photons with smaller range cuts.
So far, I am unable to explain the behaviour shown above, and I am not sure whether I did something wrong in my Gate simulation or whether this is a bug. Please find attached the Gate macros that I used to create this minimum example, and also the Python script I used to evaluate the resulting data. FYI, I use Gate version 9.1, and I also observe this behaviour with Gate version 9.0.
Any help or suggestions would be greatly appreciated. If you think that this is a bug, then I am happy to submit a bug report on GitHub.
Thanks a lot for your time in advance!
Kind regards,
David Leibold
_______________________________________________
Gate-users mailing list
Gate-users at lists.opengatecollaboration.org<mailto:Gate-users at lists.opengatecollaboration.org>
https://urldefense.proofpoint.com/v2/url?u=http-3A__lists.opengatecollaboration.org_mailman_listinfo_gate-2Dusers&d=DwIGaQ&c=XYzUhXBD2cD-CornpT4QE19xOJBbRy-TBPLK0X9U2o8&r=lfXR69GFfS7NT-Wp5HZqbtBDrbYoilDdmNEG7fKW7aM&m=u-88sQlFkb5QZWIFRFDGaf23kouGI_yHp1d_nyjk8ofTcTaowkzO2vXlFxwPJ6y1&s=goJiMnDgWKBlBQrrhCtVaHQH_jWiQXSi-qQI3yGY54M&e=
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.opengatecollaboration.org/pipermail/gate-users/attachments/20220517/92540b1c/attachment-0001.html>
More information about the Gate-users
mailing list