[Gate-users] Possible Bug: Gate Tree Output and Multiple Systems
Josh Knowland
jknowland at lucernodynamics.com
Thu May 5 00:06:02 CEST 2022
Hello all,
I think I have found a bug in the new tree output functions, but I haven't confirmed it.
I'm using GATE 9.1 with GEANT4 10.7.0. I looked at the release notes for Gate 9.2, but didn't see anything that might fix my problem.
I'm trying to run a multisystem simulation using the new tree output with numpy. The numpy output includes a column that indicates which system each hit/single is associated with - very helpful! However, a (small) percentage of my singles are being assigned to the wrong system based on their recorded coordinates. (or, maybe the coordinates are wrong and the systemID is right...)
I've made a simple macro that recreates the problem:
Two cubes located at X = -10 and X = +10. Both are defined as separate "scanner" type systems.
A point source located at X = 0 emitting individual 511 keV gammas (not back2back).
Output is set to "tree" with both hits and singles enabled.
Based on this macro, all hits/singles for System_0 should have X values less than 0 and all hits/singles for System_1 should have X values greater than 0. I've written a python script to check the numpy arrays for events that don't meet this requirement and also write the numpy arrays to text files for inspection.
Interestingly, I do not see any hits that are assigned to the wrong system, but I do see singles that are incorrectly assigned (going both ways - X values less than 0 being assigned to System_1 and vice versa).
I have looked through code in GateToTree.cc and GateOutputMgr.cc but haven't seen where this might be happening.
Can anyone help me with this? I have attached my Gate macro and my python script below.
Thanks!
Josh
/vis/disable
/gate/physics/Gamma/SetCutInRegion world 1 mm
/gate/physics/Electron/SetCutInRegion world 1 mm
/gate/physics/Positron/SetCutInRegion world 1 mm
#############################################
# WORLD
#############################################
/gate/geometry/setMaterialDatabase Materials.txt
/gate/world/geometry/setXLength 500 mm
/gate/world/geometry/setYLength 500 mm
/gate/world/geometry/setZLength 500 mm
/gate/world/setMaterial Air
#############################################
# DET1
#############################################
/gate/world/daughters/name DET1
/gate/world/daughters/systemType scanner
/gate/world/daughters/insert box
/gate/DET1/geometry/setXLength 5 mm
/gate/DET1/geometry/setYLength 5 mm
/gate/DET1/geometry/setZLength 5 mm
/gate/DET1/placement/setTranslation -10 0 0 mm
/gate/DET1/setMaterial BGO
#############################################
# DET2
#############################################
/gate/world/daughters/name DET2
/gate/world/daughters/systemType scanner
/gate/world/daughters/insert box
/gate/DET2/geometry/setXLength 5 mm
/gate/DET2/geometry/setYLength 5 mm
/gate/DET2/geometry/setZLength 5 mm
/gate/DET2/placement/setTranslation 10 0 0 mm
/gate/DET2/setMaterial BGO
#############################################
# Sensitive Detectors
#############################################
/gate/DET1/attachCrystalSD
/gate/DET2/attachCrystalSD
#############################################
# PHYSICS
#############################################
/gate/physics/addPhysicsList emstandard
/gate/physics/addProcess RadioactiveDecay
#############################################
# INITIALISATION
#############################################
/gate/random/setEngineName MersenneTwister
/gate/random/setEngineSeed 0
/gate/run/initialize
#############################################
# SOURCE
#############################################
/gate/source/addSource src_SRC1
/gate/source/src_SRC1/gps/particle gamma
/gate/source/src_SRC1/gps/ene/type Mono
/gate/source/src_SRC1/gps/ene/mono 511 keV
/gate/source/src_SRC1/setForcedUnstableFlag True
/gate/source/src_SRC1/setForcedHalfLife 6586.2 s
/gate/source/src_SRC1/setActivity 100000 Bq
/gate/source/src_SRC1/gps/pos/type Point
/gate/source/src_SRC1/gps/pos/centre 0.000 0.000 0.000 mm
/gate/source/src_SRC1/gps/ang/type iso
#############################################
# DIGITIZER
#############################################
/gate/digitizer/Singles/insert adder
#############################################
# OUTPUT
#############################################
/gate/output/tree/enable
/gate/output/tree/addFileName multisystem.npy
/gate/output/tree/hits/enable
/gate/output/tree/addCollection Singles
#############################################
# EXPERIMENT
#############################################
/gate/application/setTimeSlice 1.0000 s
/gate/application/setTimeStart 0.0000 s
/gate/application/setTimeStop 1.0000 s
/gate/application/startDAQ
import numpy as np
print("Checking SINGLES...")
singles = np.load("multisystem.Singles.npy")
for i in range(singles.size):
row = singles[i]
X = row["globalPosX"]
systemID = row["systemID"]
if(((systemID == 0) and (X > 0)) or ((systemID == 1) and (X < 0))):
print(" Problem on row " + str(i) + ": globalPosX= " + str(X) + ", systemID= " + str(systemID))
print("Exporting SINGLES...")
f=open("singles.txt",'a')
for n in singles.dtype.names:
f.write(str(n) + '\t')
f.write('\n')
np.savetxt(f,singles,fmt='%i %i %i %f %f %f %f %f %f %i %i %i %i %i %i %i %i %i %i %f %f %i %i %i %i %s %s %f %f %i')
f.close()
print("Checking HITS...")
hits = np.load("multisystem.hits.npy")
for i in range(hits.size):
row = hits[i]
X = row["posX"]
systemID = row["systemID"]
if(((systemID == 0) and (X > 0)) or ((systemID == 1) and (X < 0))):
print(" Problem on row " + str(i) + ": posX= " + str(X) + ", systemID= " + str(systemID))
print("Exporting HITS...")
f=open("hits.txt",'a')
for n in hits.dtype.names:
f.write(str(n) + '\t')
f.write('\n')
np.savetxt(f,hits,fmt='%i %i %i %f %f %i %i %i %i %f %f %f %f %f %f %f %f %f %f %f %f %f %f %s %s %s %i %i %i %i %i %i %i %i %i %i %f %f %f %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i')
f.close()
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.opengatecollaboration.org/pipermail/gate-users/attachments/20220504/a66d3d01/attachment-0001.html>
More information about the Gate-users
mailing list