Add unit tests and validation plots for bxdecay

docs/validation/CMakeLists.txt

@@ -30,6 +30,8 @@ endforeach()
 # TODO: add new images here
 # cmake-format: off
 set(TESTS_IMAGES
+    bxdecay0/double-beta-e-combined-primary-energy.output.png
+    bxdecay0/double-beta-e-primary-opening-angle.output.png
     confinement/native-surface.output.jpeg
     confinement/native-volume.output.jpeg
     confinement/complex-surface.output.jpeg

docs/validation/bxdecay0.md

@@ -0,0 +1,37 @@
+# Double-beta decay physics
+
+In this section we take a look at the primary spectrum produced by the
+_bxdecay0_ package. The main documentation for this extension is available
+[here](https://github.com/BxCppDev/bxdecay0). It aims to reproduce reasonable
+primary particles to be used as generator for double beta decay physics.
+
+## Primary electron properties
+
+The main thing we can compare here are the properties of the primary electrons
+produced. The _bxdecay0_ package actually does not even produce the neutrinos
+emitted in e.g. $2\nu\beta\beta$.
+
+```{subfigure} AB
+:subcaptions: above
+
+:::{image} ./_img/bxdecay0/double-beta-e-combined-primary-energy.output.png
+:width: 95%
+:alt: The combined energy of the two primary electrons.
+:::
+
+:::{image} ./_img/bxdecay0/double-beta-e-primary-opening-angle.output.png
+:width: 95%
+:alt: The opening angle between the two primary electrons.
+:::
+
+Different decay modes for the Ge76 decay.
+```
+
+Here $2\nu\beta\beta$ is the two neutrino double beta decay as predicted by the
+Standard Model. $0\nu\beta\beta$ is the neutrinoless double beta decay channel
+(assuming light majorana neutrino exchange) used for e.g. sensitivity estimates
+by experiments like LEGEND. More exotic physics like the Majoron-emitting
+neutrinoless double beta decay with SI=1 $0\nu\beta\beta\_M1$ shows a
+distinctive continuous electron spectrum. Other exotic physics like the
+$0\nu\beta\beta\_\lambda0$ is a decay mode driven by right-handed currents
+according to the lambda mechanism and shows a distinctive angular correlation.

docs/validation/index.md

@@ -13,6 +13,7 @@ versioninfo
 :caption: Sections
 :maxdepth: 2
 
+bxdecay0
 confinement
 distances
 nist

include/RMGSteppingAction.hh

@@ -41,6 +41,7 @@ class RMGSteppingAction : public G4UserSteppingAction {
 
   private:
 
+    bool fSkipTracking = false;
     double fDaughterKillLifetime = -1;
 
     std::unique_ptr<G4GenericMessenger> fMessenger;

src/RMGSteppingAction.cc

@@ -27,6 +27,11 @@ RMGSteppingAction::RMGSteppingAction(RMGEventAction*) { this->DefineCommands();
 
 void RMGSteppingAction::UserSteppingAction(const G4Step* step) {
 
+  if (fSkipTracking) {
+    step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
+    return;
+  }
+
   // Kill _daughter_ nuclei with a lifetime longer than a user-defined threshold. This applies to
   // the defined half-life of the particle, and not the sampled time to the decay of the secondary
   // nucleus.
@@ -92,6 +97,12 @@ void RMGSteppingAction::DefineCommands() {
       .SetParameterName("max_lifetime", false)
       .SetDefaultValue("-1")
       .SetStates(G4State_Idle);
+
+  fMessenger->DeclareProperty("SkipTracking", fSkipTracking)
+      .SetGuidance("Immediately discard tracks after primary particle generation. This feature is meant for debugging primary generation.")
+      .SetParameterName("boolean", true)
+      .SetDefaultValue("true")
+      .SetStates(G4State_Idle);
 }
 
 

tests/CMakeLists.txt

@@ -24,3 +24,7 @@ add_subdirectory(vertex)
 if(RMG_BUILD_EXAMPLES)
   add_subdirectory(examples)
 endif()
+
+if(BxDecay0_FOUND)
+  add_subdirectory(bxdecay0)
+endif()

tests/bxdecay0/CMakeLists.txt

@@ -0,0 +1,19 @@
+# collect auxiliary files
+file(
+  GLOB _aux
+  RELATIVE ${PROJECT_SOURCE_DIR}
+  macros/* gdml/*.gdml gdml/*.xml run-test-*.sh *.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 bxdecay0/gen-output COMMAND ${PYTHONPATH} ./run_sims.py ${REMAGE_PYEXE})
+set_tests_properties(bxdecay0/gen-output PROPERTIES LABELS extra FIXTURES_SETUP
+                                                    bxdecay0-output-fixture)
+
+add_test(NAME bxdecay0/plot-output COMMAND ${PYTHONPATH} plot_decay_modes.py)
+set_tests_properties(
+  bxdecay0/plot-output PROPERTIES LABELS extra FIXTURES_REQUIRED bxdecay0-output-fixture
+                                  ENVIRONMENT MPLCONFIGDIR=${CMAKE_SOURCE_DIR}/tests)

tests/bxdecay0/gdml/geometry.gdml

@@ -0,0 +1,31 @@
+<?xml version="1.0" ?>
+<gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">
+	<define/>
+	<materials/>
+	<solids>
+		<box name="world" x="0.1" y="0.1" z="0.1" lunit="m"/>
+		<orb name="basic_orb" r="0.04" lunit="m"/>
+	</solids>
+	<structure>
+		<volume name="germanium_det1">
+			<materialref ref="G4_Ge"/>
+			<solidref ref="basic_orb"/>
+		</volume>
+		<volume name="world">
+			<materialref ref="G4_Galactic"/>
+			<solidref ref="world"/>
+			<physvol name="germanium_det1">
+				<volumeref ref="germanium_det1"/>
+			</physvol>
+		</volume>
+	</structure>
+	<userinfo>
+		<auxiliary auxtype="RMG_detector_meta" auxvalue=""/>
+		<auxiliary auxtype="RMG_detector" auxvalue="germanium">
+			<auxiliary auxtype="germanium_det1" auxvalue="101"/>
+		</auxiliary>
+	</userinfo>
+	<setup name="Default" version="1.0">
+		<world ref="world"/>
+	</setup>
+</gdml>

tests/bxdecay0/macros/template.mac

@@ -0,0 +1,15 @@
+/RMG/Geometry/GDMLDisableOverlapCheck true
+/RMG/Geometry/RegisterDetector Germanium germanium_det1 001
+
+/run/initialize
+
+/RMG/Processes/Stepping/SkipTracking
+/RMG/Output/Vertex/StorePrimaryParticleInformation
+
+/RMG/Generator/Confine Volume
+/RMG/Generator/Confinement/Physical/AddVolume germanium_det1
+
+/RMG/Generator/Select BxDecay0
+/RMG/Generator/BxDecay0/DoubleBetaDecay Ge76 {MODE}
+
+/run/beamOn 250000

tests/bxdecay0/macros/template2.mac

@@ -0,0 +1,16 @@
+/RMG/Geometry/GDMLDisableOverlapCheck true
+/RMG/Geometry/RegisterDetector Germanium germanium_det1 001
+
+
+/run/initialize
+
+/RMG/Processes/Stepping/SkipTracking
+/RMG/Output/Vertex/StorePrimaryParticleInformation
+
+/RMG/Generator/Select BxDecay0
+/RMG/Generator/BxDecay0/DoubleBetaDecay Ge76 {MODE}
+
+/RMG/Generator/Confine Volume
+/RMG/Generator/Confinement/Physical/AddVolume germanium_det1
+
+/run/beamOn 250000

tests/bxdecay0/plot_decay_modes.py

@@ -0,0 +1,109 @@
+from __future__ import annotations
+
+import awkward as ak
+import hist
+import matplotlib.pyplot as plt
+import numpy as np
+from lgdo import lh5
+from numpy.linalg import norm
+
+
+# In units of keV
+def get_summed_primary_ekin(filename):
+    data = lh5.read_as("/particles", filename, "ak")
+
+    evt = ak.unflatten(data, ak.run_lengths(data.evtid))
+    # Technically we would have to filter on electrons, but the neutrinos are not even created
+    return ak.sum(evt.ekin * 1000, axis=-1)  # convert from MeV into KeV
+
+
+def plot_energy():
+    fig, ax = plt.subplots()
+
+    ekin_0vbb = get_summed_primary_ekin("0vbb.lh5")[
+        :4000
+    ]  # limit to 4000 events for good looking plots
+    ekin_2vbb = get_summed_primary_ekin("2vbb.lh5")
+    ekin_0vbb_M1 = get_summed_primary_ekin("0vbb_M1.lh5")
+    ekin_0vbb_lambda_0 = get_summed_primary_ekin("0vbb_lambda_0.lh5")[
+        :2500
+    ]  # limit to 2500 events for good looking plots
+
+    modes = [
+        r"$0\nu\beta\beta$",
+        r"$2\nu\beta\beta$",
+        r"$0\nu\beta\beta\_M1$",
+        r"$0\nu\beta\beta\_\lambda0$",
+    ]
+    ekins = [ekin_0vbb, ekin_2vbb, ekin_0vbb_M1, ekin_0vbb_lambda_0]
+
+    for decay_mode, energies in zip(modes, ekins):  # in keV
+        h = hist.new.Reg(80, 0, 2200, name="Summed electron energy [keV]").Double()
+        h.fill(energies)
+        h.plot(ax=ax, yerr=False, label=f"{decay_mode}")
+
+    ax.set_xlabel("Combined primary electron energy [keV]")
+    ax.set_ylabel("Density")
+    ax.legend(loc="upper right")
+    ax.grid()
+
+    fig.savefig("double-beta-e-combined-primary-energy.output.png")
+
+
+def get_primary_electron_angle(filename):
+    data = lh5.read_as("/particles", filename, "ak")
+
+    evt = ak.unflatten(data, ak.run_lengths(data.evtid))
+    px = evt.px
+    py = evt.py
+    pz = evt.pz
+
+    # Build 3-vectors for each particle
+    p1 = np.stack(
+        [ak.to_numpy(px[:, 0]), ak.to_numpy(py[:, 0]), ak.to_numpy(pz[:, 0])], axis=1
+    )
+    p2 = np.stack(
+        [ak.to_numpy(px[:, 1]), ak.to_numpy(py[:, 1]), ak.to_numpy(pz[:, 1])], axis=1
+    )
+    # Dot product and norms
+    dot = np.einsum("ij,ij->i", p1, p2)
+    norm1 = norm(p1, axis=1)
+    norm2 = norm(p2, axis=1)
+    # Angle in radians
+    return np.arccos(np.clip(dot / (norm1 * norm2), -1.0, 1.0))
+
+
+def plot_angles():
+    fig, ax = plt.subplots()
+
+    angles_0vbb = get_primary_electron_angle("0vbb.lh5")
+    angles_2vbb = get_primary_electron_angle("2vbb.lh5")
+    angles_0vbb_M1 = get_primary_electron_angle("0vbb_M1.lh5")
+    angles_0vbb_lambda_0 = get_primary_electron_angle("0vbb_lambda_0.lh5")
+
+    modes = [
+        r"$0\nu\beta\beta$",
+        r"$2\nu\beta\beta$",
+        r"$0\nu\beta\beta\_M1$",
+        r"$0\nu\beta\beta\_\lambda0$",
+    ]
+    angles = [angles_0vbb, angles_2vbb, angles_0vbb_M1, angles_0vbb_lambda_0]
+
+    for decay_mode, angle in zip(modes, angles):
+        # convert to degree
+        h = hist.new.Reg(
+            80, 0, 182, name="Angle between primary electrons [°]"
+        ).Double()
+        h.fill(angle * 180 / np.pi)
+        h.plot(ax=ax, yerr=False, label=f"{decay_mode}")
+
+    ax.set_xlabel("Angle between primary electrons [°]")
+    ax.set_ylabel("Density")
+    ax.legend()
+    ax.grid()
+
+    fig.savefig("double-beta-e-primary-opening-angle.output.png")
+
+
+plot_angles()
+plot_energy()

tests/bxdecay0/run_sims.py

@@ -0,0 +1,28 @@
+from __future__ import annotations
+
+from remage import remage_run
+
+runs = ["2vbb", "0vbb_M1", "0vbb_lambda_0"]
+
+# Run the template 2 once to test the different order of vertex gen macros
+remage_run(
+    macros="macros/template2.mac",
+    gdml_files="gdml/geometry.gdml",
+    output="0vbb.lh5",
+    overwrite_output=True,
+    macro_substitutions={
+        "MODE": "0vbb",
+    },
+)
+
+# Run the other three
+for run_name in runs:
+    remage_run(
+        macros="macros/template.mac",
+        gdml_files="gdml/geometry.gdml",
+        output=f"{run_name}.lh5",
+        overwrite_output=True,
+        macro_substitutions={
+            "MODE": run_name,
+        },
+    )