[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