<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 15 (filtered medium)">
<style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-family:"Calibri",sans-serif;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="#0563C1" vlink="#954F72" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal">Hello all,<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I think I have found a bug in the new tree output functions, but I haven’t confirmed it.<o:p></o:p></p>
<p class="MsoNormal">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.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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…)<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I’ve made a simple macro that recreates the problem:<o:p></o:p></p>
<p class="MsoNormal" style="text-indent:.5in">Two cubes located at X = -10 and X = +10. Both are defined as separate “scanner” type systems.<o:p></o:p></p>
<p class="MsoNormal" style="text-indent:.5in">A point source located at X = 0 emitting individual 511 keV gammas (not back2back).<o:p></o:p></p>
<p class="MsoNormal" style="text-indent:.5in">Output is set to “tree” with both hits and singles enabled.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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).<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">I have looked through code in GateToTree.cc and GateOutputMgr.cc but haven’t seen where this might be happening.
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Can anyone help me with this? I have attached my Gate macro and my python script below.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Thanks!<o:p></o:p></p>
<p class="MsoNormal">Josh<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">/vis/disable<o:p></o:p></p>
<p class="MsoNormal">/gate/physics/Gamma/SetCutInRegion world 1 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/physics/Electron/SetCutInRegion world 1 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/physics/Positron/SetCutInRegion world 1 mm<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># WORLD<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/geometry/setMaterialDatabase Materials.txt<o:p></o:p></p>
<p class="MsoNormal">/gate/world/geometry/setXLength 500 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/world/geometry/setYLength 500 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/world/geometry/setZLength 500 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/world/setMaterial Air<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># DET1<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/world/daughters/name DET1<o:p></o:p></p>
<p class="MsoNormal">/gate/world/daughters/systemType scanner<o:p></o:p></p>
<p class="MsoNormal">/gate/world/daughters/insert box<o:p></o:p></p>
<p class="MsoNormal">/gate/DET1/geometry/setXLength 5 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/DET1/geometry/setYLength 5 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/DET1/geometry/setZLength 5 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/DET1/placement/setTranslation -10 0 0 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/DET1/setMaterial BGO<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># DET2<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/world/daughters/name DET2<o:p></o:p></p>
<p class="MsoNormal">/gate/world/daughters/systemType scanner<o:p></o:p></p>
<p class="MsoNormal">/gate/world/daughters/insert box<o:p></o:p></p>
<p class="MsoNormal">/gate/DET2/geometry/setXLength 5 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/DET2/geometry/setYLength 5 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/DET2/geometry/setZLength 5 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/DET2/placement/setTranslation 10 0 0 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/DET2/setMaterial BGO<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># Sensitive Detectors<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/DET1/attachCrystalSD<o:p></o:p></p>
<p class="MsoNormal">/gate/DET2/attachCrystalSD<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># PHYSICS<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/physics/addPhysicsList emstandard<o:p></o:p></p>
<p class="MsoNormal">/gate/physics/addProcess RadioactiveDecay<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># INITIALISATION<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/random/setEngineName MersenneTwister<o:p></o:p></p>
<p class="MsoNormal">/gate/random/setEngineSeed 0<o:p></o:p></p>
<p class="MsoNormal">/gate/run/initialize<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># SOURCE<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/source/addSource src_SRC1<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/gps/particle gamma<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/gps/ene/type Mono<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/gps/ene/mono 511 keV<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/setForcedUnstableFlag True<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/setForcedHalfLife 6586.2 s<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/setActivity 100000 Bq<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/gps/pos/type Point<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/gps/pos/centre 0.000 0.000 0.000 mm<o:p></o:p></p>
<p class="MsoNormal">/gate/source/src_SRC1/gps/ang/type iso<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># DIGITIZER<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/digitizer/Singles/insert adder<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># OUTPUT<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/output/tree/enable<o:p></o:p></p>
<p class="MsoNormal">/gate/output/tree/addFileName multisystem.npy<o:p></o:p></p>
<p class="MsoNormal">/gate/output/tree/hits/enable<o:p></o:p></p>
<p class="MsoNormal">/gate/output/tree/addCollection Singles<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal"># EXPERIMENT<o:p></o:p></p>
<p class="MsoNormal">#############################################<o:p></o:p></p>
<p class="MsoNormal">/gate/application/setTimeSlice 1.0000 s<o:p></o:p></p>
<p class="MsoNormal">/gate/application/setTimeStart 0.0000 s<o:p></o:p></p>
<p class="MsoNormal">/gate/application/setTimeStop 1.0000 s<o:p></o:p></p>
<p class="MsoNormal">/gate/application/startDAQ<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">import numpy as np<o:p></o:p></p>
<p class="MsoNormal">print("Checking SINGLES...")<o:p></o:p></p>
<p class="MsoNormal">singles = np.load("multisystem.Singles.npy")<o:p></o:p></p>
<p class="MsoNormal">for i in range(singles.size):<o:p></o:p></p>
<p class="MsoNormal"> row = singles[i]<o:p></o:p></p>
<p class="MsoNormal"> X = row["globalPosX"]<o:p></o:p></p>
<p class="MsoNormal"> systemID = row["systemID"]<o:p></o:p></p>
<p class="MsoNormal"> if(((systemID == 0) and (X > 0)) or ((systemID == 1) and (X < 0))):<o:p></o:p></p>
<p class="MsoNormal"> print(" Problem on row " + str(i) + ": globalPosX= " + str(X) + ", systemID= " + str(systemID))<o:p></o:p></p>
<p class="MsoNormal">print("Exporting SINGLES...")<o:p></o:p></p>
<p class="MsoNormal">f=open("singles.txt",'a')<o:p></o:p></p>
<p class="MsoNormal">for n in singles.dtype.names:<o:p></o:p></p>
<p class="MsoNormal"> f.write(str(n) + '\t')<o:p></o:p></p>
<p class="MsoNormal">f.write('\n')<o:p></o:p></p>
<p class="MsoNormal">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')<o:p></o:p></p>
<p class="MsoNormal">f.close()<o:p></o:p></p>
<p class="MsoNormal">print("Checking HITS...")<o:p></o:p></p>
<p class="MsoNormal">hits = np.load("multisystem.hits.npy")<o:p></o:p></p>
<p class="MsoNormal">for i in range(hits.size):<o:p></o:p></p>
<p class="MsoNormal"> row = hits[i]<o:p></o:p></p>
<p class="MsoNormal"> X = row["posX"]<o:p></o:p></p>
<p class="MsoNormal"> systemID = row["systemID"]<o:p></o:p></p>
<p class="MsoNormal"> if(((systemID == 0) and (X > 0)) or ((systemID == 1) and (X < 0))):<o:p></o:p></p>
<p class="MsoNormal"> print(" Problem on row " + str(i) + ": posX= " + str(X) + ", systemID= " + str(systemID))<o:p></o:p></p>
<p class="MsoNormal">print("Exporting HITS...")<o:p></o:p></p>
<p class="MsoNormal">f=open("hits.txt",'a')<o:p></o:p></p>
<p class="MsoNormal">for n in hits.dtype.names:<o:p></o:p></p>
<p class="MsoNormal"> f.write(str(n) + '\t')<o:p></o:p></p>
<p class="MsoNormal">f.write('\n')<o:p></o:p></p>
<p class="MsoNormal">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')<o:p></o:p></p>
<p class="MsoNormal">f.close()<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
</body>
</html>