diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 144fd57..f3e419c 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -42,7 +42,8 @@ jobs: t29_optical_digi, t30_dna, t31_vpgTLE-tt, - t32_isotopes] + t32_isotopes, + t33_invert_filter] steps: - name: Checkout github repo diff --git a/t33_invert_filter/data/GateMaterials.db b/t33_invert_filter/data/GateMaterials.db new file mode 100644 index 0000000..5e97f1b --- /dev/null +++ b/t33_invert_filter/data/GateMaterials.db @@ -0,0 +1,27 @@ +[Elements] +Hydrogen: S= H ; Z= 1. ; A= 1.01 g/mole +Carbon: S= C ; Z= 6. ; A= 12.01 g/mole +Nitrogen: S= N ; Z= 7. ; A= 14.01 g/mole +Oxygen: S= O ; Z= 8. ; A= 16.00 g/mole +Argon: S= Ar ; Z= 18. ; A= 39.95 g/mole +Tungsten: S= W ; Z= 74. ; A= 183.84 g/mole +Lead: S= Pb ; Z= 82. ; A= 207.20 g/mole + +[Materials] +Vacuum: d=0.000001 mg/cm3 ; n=1 + +el: name=Hydrogen ; n=1 + +Super_dense: d=19300 g/cm3 ; n=1 ; state=solid + +el: name=Tungsten ; n=1 + +PWO: d=8.28 g/cm3; n=3 ; state=Solid + +el: name=Lead ; n=1 + +el: name=Tungsten ; n=1 + +el: name=Oxygen ; n=4 + +Air: d=1.29 mg/cm3 ; n=4 ; state=gas + +el: name=Nitrogen ; f=0.755268 + +el: name=Oxygen ; f=0.231781 + +el: name=Argon ; f=0.012827 + +el: name=Carbon ; f=0.000124 + diff --git a/t33_invert_filter/data/spectrum.txt b/t33_invert_filter/data/spectrum.txt new file mode 100644 index 0000000..0f8346f --- /dev/null +++ b/t33_invert_filter/data/spectrum.txt @@ -0,0 +1,3 @@ +1 0 +0.1 0.5 +0.2 0.5 \ No newline at end of file diff --git a/t33_invert_filter/mac/actor.mac b/t33_invert_filter/mac/actor.mac new file mode 100644 index 0000000..41e3bad --- /dev/null +++ b/t33_invert_filter/mac/actor.mac @@ -0,0 +1,53 @@ + +/control/doif {filter_id} > -1 /gate/actor/addActor KillActor myKillActor +/control/doif {filter_id} > -1 /gate/actor/myKillActor/attachTo phantom +#/control/doif {filter_id} > -1 /gate/actor/myKillActor/save output/kill.txt + +/control/doif {filter_id} == 0 /gate/actor/myKillActor/addFilter angleFilter +/control/doif {filter_id} == 0 /gate/actor/myKillActor/angleFilter/setAngle 30 +/control/doif {filter_id} == 0 /gate/actor/myKillActor/angleFilter/setDirection 0 0 1 + +/control/doif {filter_id} == 1 /gate/actor/myKillActor/addFilter angleFilter +/control/doif {filter_id} == 1 /gate/actor/myKillActor/angleFilter/setAngle 30 +/control/doif {filter_id} == 1 /gate/actor/myKillActor/angleFilter/setDirection 0 0 1 +/control/doif {filter_id} == 1 /gate/actor/myKillActor/angleFilter/invert + +/control/doif {filter_id} == 2 /gate/actor/myKillActor/addFilter energyFilter +/control/doif {filter_id} == 2 /gate/actor/myKillActor/energyFilter/setEmin 0.195 MeV +/control/doif {filter_id} == 2 /gate/actor/myKillActor/energyFilter/setEmax 0.205 MeV + +/control/doif {filter_id} == 3 /gate/actor/myKillActor/addFilter energyFilter +/control/doif {filter_id} == 3 /gate/actor/myKillActor/energyFilter/setEmin 0.195 MeV +/control/doif {filter_id} == 3 /gate/actor/myKillActor/energyFilter/setEmax 0.205 MeV +/control/doif {filter_id} == 3 /gate/actor/myKillActor/energyFilter/invert + +/control/doif {filter_id} == 4 /gate/actor/myKillActor/addFilter materialFilter +/control/doif {filter_id} == 4 /gate/actor/myKillActor/materialFilter/addMaterial Air +/control/doif {filter_id} == 4 /gate/phantom_part/setMaterial Air + + +/control/doif {filter_id} == 5 /gate/actor/myKillActor/addFilter materialFilter +/control/doif {filter_id} == 5 /gate/actor/myKillActor/materialFilter/addMaterial Air +/control/doif {filter_id} == 5 /gate/actor/myKillActor/materialFilter/invert +/control/doif {filter_id} == 5 /gate/phantom_part/setMaterial Air + +/control/doif {filter_id} == 6 /gate/actor/myKillActor/addFilter volumeFilter +/control/doif {filter_id} == 6 /gate/actor/myKillActor/volumeFilter/addVolume phantom_part + + +/control/doif {filter_id} == 7 /gate/actor/myKillActor/addFilter volumeFilter +/control/doif {filter_id} == 7 /gate/actor/myKillActor/volumeFilter/addVolume phantom_part +/control/doif {filter_id} == 7 /gate/actor/myKillActor/volumeFilter/invert + +# IDfilter is broken, so no test needed +#/control/doif {filter_id} == 8 /gate/actor/myKillActor/addFilter IDFilter +#/control/doif {filter_id} == 9 /gate/actor/myKillActor/addFilter IDFilter +#/control/doif {filter_id} == 9 /gate/actor/myKillActor/IDFilter/invert + +/control/doif {filter_id} == 10 /gate/actor/myKillActor/addFilter particleFilter +/control/doif {filter_id} == 10 /gate/actor/myKillActor/particleFilter/addParticle e- + +/control/doif {filter_id} == 11 /gate/actor/myKillActor/addFilter particleFilter +/control/doif {filter_id} == 11 /gate/actor/myKillActor/particleFilter/addParticle e- +/control/doif {filter_id} == 11 /gate/actor/myKillActor/particleFilter/invert + diff --git a/t33_invert_filter/mac/geometry.mac b/t33_invert_filter/mac/geometry.mac new file mode 100644 index 0000000..647d30f --- /dev/null +++ b/t33_invert_filter/mac/geometry.mac @@ -0,0 +1,53 @@ +/gate/world/geometry/setXLength 20 cm +/gate/world/geometry/setYLength 20 cm +/gate/world/geometry/setZLength 40 cm +/gate/world/vis/setVisible 0 +/gate/world/setMaterial Vacuum + +/gate/world/daughters/name scanner +/gate/world/daughters/systemType scanner +/gate/world/daughters/insert box +/gate/scanner/placement/setTranslation 0 0 15 cm +/gate/scanner/geometry/setXLength 20 cm +/gate/scanner/geometry/setYLength 20 cm +/gate/scanner/geometry/setZLength 10 cm +/gate/scanner/vis/setColor red +/gate/scanner/vis/setVisible 1 + +/gate/scanner/daughters/name crystal +/gate/scanner/daughters/insert box +/gate/crystal/setMaterial Super_dense +/control/doif {filter_id} > 9 /gate/crystal/setMaterial PWO +/gate/crystal/placement/setTranslation 0 0 0 cm +/gate/crystal/geometry/setXLength 20 cm +/gate/crystal/geometry/setYLength 20 cm +/gate/crystal/geometry/setZLength 10 cm + +/gate/world/daughters/name phantom +/gate/world/daughters/insert box +/gate/phantom/placement/setTranslation 0 0 0 cm +/gate/phantom/geometry/setXLength 20 cm +/gate/phantom/geometry/setYLength 20 cm +/gate/phantom/geometry/setZLength 1 cm +/gate/phantom/vis/setColor blue +/gate/phantom/vis/setVisible 1 +/gate/phantom/setMaterial Vacuum + +/gate/phantom/daughters/name phantom_part +/gate/phantom/daughters/insert box +/gate/phantom_part/placement/setTranslation 5 0 0 cm +/gate/phantom_part/geometry/setXLength 10 cm +/gate/phantom_part/geometry/setYLength 20 cm +/gate/phantom_part/geometry/setZLength 1 cm +/gate/phantom_part/vis/setColor red +/gate/phantom_part/vis/setVisible 1 +/gate/phantom_part/setMaterial Vacuum + +/gate/crystal/attachCrystalSD +/gate/systems/scanner/level1/attach crystal +/gate/phantom/attachPhantomSD +/gate/phantom_part/attachPhantomSD + +#/gate/digitizerMgr/crystal/SinglesDigitizer/Singles/insert adder +#/gate/digitizerMgr/crystal/SinglesDigitizer/Singles/insert readout +#/gate/digitizerMgr/crystal/SinglesDigitizer/Singles/readout/setDepth 1 \ No newline at end of file diff --git a/t33_invert_filter/mac/main.mac b/t33_invert_filter/mac/main.mac new file mode 100644 index 0000000..d8cc142 --- /dev/null +++ b/t33_invert_filter/mac/main.mac @@ -0,0 +1,40 @@ +/control/verbose 1 +# basic geometry + +/gate/geometry/setMaterialDatabase data/GateMaterials.db + +/control/execute mac/geometry.mac + +/gate/random/setEngineName MersenneTwister +#fixed seed for reproducibility +#/gate/random/setEngineSeed 15498 +/gate/random/setEngineSeed auto +/gate/random/verbose 0 + +# physics +/gate/physics/addPhysicsList emstandard_opt3 +/gate/physics/setDEDXBinning 500 +/gate/physics/setLambdaBinning 500 +/gate/physics/processList Enabled +/gate/physics/processList Initialized + +#actors + +/control/execute mac/actor.mac + +/gate/run/initialize + +/control/execute mac/source.mac +/control/if {filter_id} > 9 mac/source_electron.mac +#/control/execute mac/source_electron.mac + +# output +/control/execute mac/output.mac + +#/control/execute mac/vis.mac + +/gate/application/setTimeStart 0.0 s +/gate/application/setTimeStop 1 s +/gate/application/setTimeSlice 1 s + +/gate/application/startDAQ \ No newline at end of file diff --git a/t33_invert_filter/mac/output.mac b/t33_invert_filter/mac/output.mac new file mode 100644 index 0000000..248af13 --- /dev/null +++ b/t33_invert_filter/mac/output.mac @@ -0,0 +1,44 @@ + +/gate/output/tree/enable +#/gate/output/tree/addFileName ./output/filter_id_{filter_id}.txt +/gate/output/tree/addFileName ./output/filter_id_{filter_id}.npy +/gate/output/tree/hits/enable + +#PDGEncoding,trackID,parentID,trackLocalTime,time,runID,eventID,sourceID,primaryID,posX,posY,posZ,localPosX,localPosY,localPosZ,momDirX,momDirY,momDirZ,edep,stepLength,trackLength,rotationAngle,axialPos,processName,comptVolName,RayleighVolName,volumeID[0],volumeID[1],volumeID[2],volumeID[3],volumeID[4],volumeID[5],volumeID[6],volumeID[7],volumeID[8],volumeID[9],sourcePosX,sourcePosY,sourcePosZ,nPhantomCompton,nCrystalCompton,nPhantomRayleigh,nCrystalRayleigh,,level1ID,level2ID,level3ID,level4ID,level5ID,photonID,sourceType,decayType,gammaType +#/gate/output/tree/hits/branches/PDGEncoding/disable +#/gate/output/tree/hits/branches/trackID/disable +/gate/output/tree/hits/branches/parentID/disable +/gate/output/tree/hits/branches/trackLocalTime/disable +#/gate/output/tree/hits/branches/time/disable +/gate/output/tree/hits/branches/runID/disable +#/gate/output/tree/hits/branches/eventID/disable +/gate/output/tree/hits/branches/sourceID/disable +/gate/output/tree/hits/branches/primaryID/disable +#/gate/output/tree/hits/branches/posX/disable +#/gate/output/tree/hits/branches/posY/disable +#/gate/output/tree/hits/branches/posZ/disable +/gate/output/tree/hits/branches/localPosX/disable +/gate/output/tree/hits/branches/localPosY/disable +/gate/output/tree/hits/branches/localPosZ/disable +/gate/output/tree/hits/branches/momDirX/disable +/gate/output/tree/hits/branches/momDirY/disable +/gate/output/tree/hits/branches/momDirZ/disable +#/gate/output/tree/hits/branches/edep/disable +/gate/output/tree/hits/branches/stepLength/disable +/gate/output/tree/hits/branches/trackLength/disable +/gate/output/tree/hits/branches/rotationAngle/disable +/gate/output/tree/hits/branches/axialPos/disable +#/gate/output/tree/hits/branches/processName/disable +/gate/output/tree/hits/branches/comptVolName/disable +/gate/output/tree/hits/branches/RayleighVolName/disable +/gate/output/tree/hits/branches/sourcePosX/disable +/gate/output/tree/hits/branches/sourcePosY/disable +/gate/output/tree/hits/branches/sourcePosZ/disable +#/gate/output/tree/hits/branches/nPhantomCompton/disable +#/gate/output/tree/hits/branches/nCrystalCompton/disable +#/gate/output/tree/hits/branches/nPhantomRayleigh/disable +#/gate/output/tree/hits/branches/nCrystalRayleigh/disable +/gate/output/tree/hits/branches/photonID/disable +/gate/output/tree/hits/branches/sourceType/disable +/gate/output/tree/hits/branches/decayType/disable +/gate/output/tree/hits/branches/gammaType/disable diff --git a/t33_invert_filter/mac/source.mac b/t33_invert_filter/mac/source.mac new file mode 100644 index 0000000..795ab59 --- /dev/null +++ b/t33_invert_filter/mac/source.mac @@ -0,0 +1,15 @@ +/gate/source/addSource myGammaSource gps +/gate/source/myGammaSource/gps/pos/type Point +/gate/source/myGammaSource/gps/particle gamma +/gate/source/myGammaSource/gps/ene/type UserSpectrum +/gate/source/myGammaSource/gps/setSpectrumFile data/spectrum.txt +/gate/source/myGammaSource/setActivity 1000. becquerel +/gate/source/myGammaSource/gps/pos/centre 0. 0. -5. cm +/gate/source/myGammaSource/gps/ang/type iso +/gate/source/myGammaSource/gps/ang/maxtheta 180 deg +/gate/source/myGammaSource/gps/ang/mintheta 135 deg +/gate/source/myGammaSource/gps/ang/minphi 0 deg +/gate/source/myGammaSource/gps/ang/maxphi 360 deg +/gate/source/myGammaSource/attachTo world + +/gate/source/myGammaSource/visualize 100 red \ No newline at end of file diff --git a/t33_invert_filter/mac/source_electron.mac b/t33_invert_filter/mac/source_electron.mac new file mode 100644 index 0000000..5872ade --- /dev/null +++ b/t33_invert_filter/mac/source_electron.mac @@ -0,0 +1,14 @@ +/gate/source/addSource mySource gps +/gate/source/mySource/gps/pos/type Point +/gate/source/mySource/gps/particle e- +/gate/source/mySource/gps/ene/type Mono +/gate/source/mySource/gps/ene/mono 100 keV +/gate/source/mySource/setActivity 1000. becquerel +/gate/source/mySource/gps/pos/centre 0. 0. -5. cm +/gate/source/mySource/gps/ang/type iso +/gate/source/mySource/gps/ang/maxtheta 180 deg +/gate/source/mySource/gps/ang/mintheta 135 deg +/gate/source/mySource/gps/ang/minphi 0 deg +/gate/source/mySource/gps/ang/maxphi 360 deg +/gate/source/mySource/attachTo world +/gate/source/mySource/visualize 100 red \ No newline at end of file diff --git a/t33_invert_filter/mac/vis.mac b/t33_invert_filter/mac/vis.mac new file mode 100644 index 0000000..035cb90 --- /dev/null +++ b/t33_invert_filter/mac/vis.mac @@ -0,0 +1,21 @@ +# +# ****** Visualization +# + +#/vis/open OGLSX +/vis/open OGLSQt 1920x1080-0+0 +/vis/viewer/set/viewpointThetaPhi 0 0 deg +/vis/viewer/zoom 1 +/vis/viewer/set/style w +/vis/drawVolume +/vis/viewer/flush +#/tracking/verbose 0 +/tracking/storeTrajectory 1 +/vis/scene/add/hits +/vis/scene/add/trajectories +/vis/scene/endOfEventAction accumulate -1 +/vis/scene/add/axes + +#/vis/viewer/set/viewpointThetaPhi 0 0 +#vis/viewer/panTo 0 430 mm +#/vis/viewer/zoom 10 diff --git a/t33_invert_filter/readme.md b/t33_invert_filter/readme.md new file mode 100644 index 0000000..b521016 --- /dev/null +++ b/t33_invert_filter/readme.md @@ -0,0 +1,16 @@ +This test is used to check the validity of the invert command for following filters for actors: + +| Filter | normal filter_id | Invert filter_id | +| --- | --- | --- | +| AngleFilter | 0 | 1 | +| EnergyFilter | 2 | 3 | +| MaterialFilter | 4 | 5 | +| VolumeFilter | 6 | 7 | +| IDFilter | 8 | 9 | +| ParticleFilter | 10 | 11 | + + + +The test is done by using a simple geometry with a single point source, a volume called phantom, which the killer actor attach to, and a piece of material act as detector. + +Angle and energy filter's test is trivial. Material and volume filter is tested with a daughter volume of phantom. Since the IDFilter is broken right now, nothing will be done for it right now. Particle filter is tested by adding a electron source, and carefully remove all secondary hits. \ No newline at end of file diff --git a/t33_invert_filter/runAnalysis.py b/t33_invert_filter/runAnalysis.py new file mode 100644 index 0000000..a762ec1 --- /dev/null +++ b/t33_invert_filter/runAnalysis.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python3 +# ----------------------------------------------------------------------------- +# Copyright (C): OpenGATE Collaboration +# This software is distributed under the terms +# of the GNU Lesser General Public Licence (LGPL) +# See LICENSE.md for further details +# ----------------------------------------------------------------------------- + +import click +import os +import numpy as np +import gatetools as gt +import matplotlib.pyplot as plt + +CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) + +@click.command(context_settings=CONTEXT_SETTINGS) +@click.argument('output_folders', + nargs=-1, + required=True, + type=click.Path(exists=True, file_okay=True, dir_okay=True)) + +def analyse_command_line(output_folders, **kwargs): + # logger + gt.logging_conf(**kwargs) + # Run the analysis with the command line (click) + # the return code is 0 (fail) or 1 (success) + r = analyse_all_folders(output_folders) + print(f'Last test return is: {r}') + + +def analyse_all_folders(output_folders): + r = True + for folder in output_folders: + folder_basename = os.path.basename(folder) + if folder_basename.split("-")[-1] != "output": + if float(folder_basename.split("-")[-1]) < 9.4: + continue # feature add in 9.3 + r = r and analyse_folder(folder) + return r + +def analyse_folder(folder): + r = True + for filter_id in range(12): + filepath = os.path.join(folder, f"filter_id_{filter_id}.hits.npy") + if os.path.exists(filepath): + data = get_primary(np.load(filepath)) + elif filter_id != 8 and filter_id != 9: + raise(f"File {filepath} does not exist") + + if filter_id == 0 or filter_id == 1: + r = r and analyse_angle(data, filter_id == 1) + + if filter_id == 2 or filter_id == 3: + r = r and analyse_energy(data, filter_id == 3) + + if filter_id == 4 or filter_id == 5: + r = r and analyse_material(data, filter_id == 5) + + if filter_id == 6 or filter_id == 7: + r = r and analyse_volume(data, filter_id == 7) + + if filter_id == 8 or filter_id == 9: + r = r and analyse_ID(data, filter_id == 9) + + if filter_id == 10 or filter_id == 11: + r = r and analyse_particle(data, filter_id == 11) + return r + +def get_primary(data): + is_primary = data["nCrystalRayleigh"]==0 + is_primary = np.logical_and(is_primary, data["nCrystalCompton"]==0) + is_primary = np.logical_and(is_primary, data["nPhantomRayleigh"]==0) + is_primary = np.logical_and(is_primary, data["nPhantomCompton"]==0) + + return np.compress(is_primary,data) + + + +def analyse_angle(data, invert:bool): + distance2 = data["posX"]**2 + data["posY"]**2 + radius2 = 150**2*np.tan(np.pi/6)**2 + if invert: + return (distance2.max() - radius2 < 1) + else: + return (distance2.min() - radius2 > -1) + + +def analyse_energy(data, invert:bool): + if invert: + return data["edep"].min() > 0.195 + else: + return data["edep"].max() < 0.195 + +def plot_scatter(data, invert:bool, name): + plt.figure() + plt.scatter(data["posX"], data["posY"]) + plt.axis('equal') + plt.savefig(f"{name}_{invert}.svg") + +def analyse_material(data, invert:bool): + if invert: + return data["posX"].min() > 0 + else: + return data["posX"].max() < 0 + + +def analyse_volume(data, invert:bool): + if invert: + return data["posX"].min() > 0 + else: + return data["posX"].max() < 0 + +def analyse_ID(data, invert:bool): + return True + +def analyse_particle(data, invert:bool): + + # first remove secondary by bremsstrahlung effect + event_id = data["eventID"] + dupulicate = np.insert((event_id[1:]-event_id[:-1] == 0), 0, False) + data = np.compress(np.logical_not(dupulicate),data) + + if invert: + return data["PDGEncoding"].min()==11 and data["PDGEncoding"].max()==11 + else: + return data["PDGEncoding"].min()==22 and data["PDGEncoding"].max()==22 + + + + + +# ----------------------------------------------------------------------------- +if __name__ == '__main__': + analyse_command_line(["output"]) diff --git a/t33_invert_filter/runTest.sh b/t33_invert_filter/runTest.sh new file mode 100755 index 0000000..48ec6df --- /dev/null +++ b/t33_invert_filter/runTest.sh @@ -0,0 +1,12 @@ +Gate mac/main.mac -a [filter_id,0] +Gate mac/main.mac -a [filter_id,1] +Gate mac/main.mac -a [filter_id,2] +Gate mac/main.mac -a [filter_id,3] +Gate mac/main.mac -a [filter_id,4] +Gate mac/main.mac -a [filter_id,5] +Gate mac/main.mac -a [filter_id,6] +Gate mac/main.mac -a [filter_id,7] +#Gate mac/main.mac -a [filter_id,8] +#Gate mac/main.mac -a [filter_id,9] +Gate mac/main.mac -a [filter_id,10] +Gate mac/main.mac -a [filter_id,11]