Add exception for primary and basic test

src/RMGParticleFilterOutputScheme.cc

@@ -59,10 +59,19 @@ std::optional<G4ClassificationOfNewTrack> RMGParticleFilterOutputScheme::
   // If a keep volume is specified only keep particle if in that volume.
   if (fKeepVolumes.find(pv_name) != fKeepVolumes.end()) return std::nullopt;
 
+  // If this is the primary particle we can not kill it without crashing the simulation
+  if (aTrack->GetTrackID() == 0) {
+    RMGLog::OutDev(RMGLog::warning,
+        "ParticleFilter: The primary particle was specified to kill. "
+        "This would cause the simulation to crash... so simulating it anyways. ");
+    return std::nullopt;
+  }
+
   // We land here if
   // i) Particle is marked to kill.
   // ii) No Kill volume specified or the particle is in the kill volume.
   // iii) Particle is not in the keep volume.
+  // iiii) Particle is not the primary particle.
   RMGLog::OutDev(RMGLog::debug, "Filtering out particle with PDG code ", pdg,
       " in RMGParticleFilterOutputScheme");
   return fKill;

tests/CMakeLists.txt

@@ -13,6 +13,7 @@ add_subdirectory(distances)
 add_subdirectory(germanium)
 add_subdirectory(internals)
 add_subdirectory(output)
+add_subdirectory(particlefilter)
 add_subdirectory(python)
 add_subdirectory(vertex)
 

tests/particlefilter/CMakeLists.txt

@@ -0,0 +1,14 @@
+# collect auxiliary files
+file(
+  GLOB _aux
+  RELATIVE ${PROJECT_SOURCE_DIR}
+  macros/* gdml/*.gdml gdml/*.xml run-test-*.sh run-test-*.py)
+
+# copy them to the build area
+foreach(_file ${_aux})
+  configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY)
+endforeach()
+
+add_test(NAME particlefilter/test-particlefilter
+         COMMAND ${PYTHONPATH} ./run-test-particlefilter.py $<TARGET_FILE:remage-cli-cpp>)
+set_tests_properties(particlefilter/test-particlefilter PROPERTIES LABELS extra)

tests/particlefilter/gdml/geometry.gdml

@@ -0,0 +1,45 @@
+<?xml version="1.0" ?>
+<gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">
+	<define/>
+	<materials/>
+	<solids>
+		<box name="world" x="2" y="2" z="2" lunit="m"/>
+		<box name="scint" x="0.5" y="1" z="1" lunit="m"/>
+		<box name="det" x="0.1" y="0.5" z="0.5" lunit="m"/>
+	</solids>
+	<structure>
+		<volume name="det">
+			<materialref ref="G4_Ge"/>
+			<solidref ref="det"/>
+		</volume>
+		<volume name="scint1">
+			<materialref ref="G4_lAr"/>
+			<solidref ref="scint"/>
+			<physvol name="det1">
+				<volumeref ref="det"/>
+			</physvol>
+		</volume>
+		<volume name="scint2">
+			<materialref ref="G4_lAr"/>
+			<solidref ref="scint"/>
+			<physvol name="det2">
+				<volumeref ref="det"/>
+			</physvol>
+		</volume>
+		<volume name="world">
+			<materialref ref="G4_Galactic"/>
+			<solidref ref="world"/>
+			<physvol name="scint1">
+				<volumeref ref="scint1"/>
+				<position name="scint1_pos" x="-255.000000000000000" y="0.000000000000000" z="0.000000000000000" unit="mm"/>
+			</physvol>
+			<physvol name="scint2">
+				<volumeref ref="scint2"/>
+				<position name="scint2_pos" x="255.000000000000000" y="0.000000000000000" z="0.000000000000000" unit="mm"/>
+			</physvol>
+		</volume>
+	</structure>
+	<setup name="Default" version="1.0">
+		<world ref="world"/>
+	</setup>
+</gdml>

tests/particlefilter/macros/gammafilter.mac

@@ -0,0 +1,20 @@
+/RMG/Geometry/GDMLDisableOverlapCheck true
+
+/RMG/Output/ActivateOutputScheme Track
+/RMG/Output/ActivateOutputScheme ParticleFilter
+
+/run/initialize
+
+/RMG/Output/ParticleFilter/AddParticle 22
+
+/RMG/Generator/Confine UnConfined
+
+/RMG/Generator/Select GPS
+/gps/pos/type Volume
+/gps/pos/shape Sphere
+/gps/pos/radius 0.5 m
+/gps/particle e-
+/gps/energy       2 MeV
+/gps/ang/type     iso
+
+/run/beamOn 10000

tests/particlefilter/run-test-particlefilter.py

@@ -0,0 +1,24 @@
+#!/bin/env python3
+
+from __future__ import annotations
+
+import subprocess
+import sys
+
+from lgdo import lh5
+
+rmg = sys.argv[1]
+macro = "macros/gammafilter.mac"
+output_lh5 = "gammafilter.lh5"
+
+# run remage, produce lh5 output.
+subprocess.run(
+    [rmg, "-g", "gdml/geometry.gdml", "-o", output_lh5, "-w", "--", macro],
+    check=False,
+)
+
+# check that we get to stable isotopes.
+tracks = lh5.read("stp/tracks", output_lh5).view_as("pd")
+particle_numbers = tracks["particle"].value_counts()
+assert particle_numbers[11] != 0  # primaries.
+assert particle_numbers.get(22, 0) == 0  # gammas