diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index aa5ca19b9..49a9cdeda 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -271,6 +271,7 @@ Commands for controlling physics processes **Commands:** * `DaughterNucleusMaxLifetime` – Determines which unstable daughter nuclei will be killed, if they are at rest, depending on their lifetime. +* `SkipTracking` – Immediately discard tracks after primary particle generation. This feature is meant for debugging primary generation. * `ResetInitialDecayTime` – If the initial step is a radioactive decay, reset the global time of all its secondary tracks to 0. * `LargeGlobalTimeUncertaintyWarning` – Warn if the global times of tracks get too large to provide the requested time uncertainty. @@ -293,6 +294,16 @@ Set to -1 to disable this feature. * **Candidates** – `s ms us ns ps min h d y second millisecond microsecond nanosecond picosecond minute hour day year` * **Allowed states** – `Idle` +### `/RMG/Processes/Stepping/SkipTracking` + +Immediately discard tracks after primary particle generation. This feature is meant for debugging primary generation. + +* **Parameter** – `boolean` + * **Parameter type** – `b` + * **Omittable** – `True` + * **Default value** – `true` +* **Allowed states** – `Idle` + ### `/RMG/Processes/Stepping/ResetInitialDecayTime` If the initial step is a radioactive decay, reset the global time of all its secondary tracks to 0. diff --git a/docs/validation/CMakeLists.txt b/docs/validation/CMakeLists.txt index 95280175a..393d5b1e5 100644 --- a/docs/validation/CMakeLists.txt +++ b/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 diff --git a/docs/validation/bxdecay0.md b/docs/validation/bxdecay0.md new file mode 100644 index 000000000..34d978d30 --- /dev/null +++ b/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. diff --git a/docs/validation/index.md b/docs/validation/index.md index c15ce78cf..f3f66fdae 100644 --- a/docs/validation/index.md +++ b/docs/validation/index.md @@ -13,6 +13,7 @@ versioninfo :caption: Sections :maxdepth: 2 +bxdecay0 confinement distances nist diff --git a/include/RMGSteppingAction.hh b/include/RMGSteppingAction.hh index eae17a99e..63b18c044 100644 --- a/include/RMGSteppingAction.hh +++ b/include/RMGSteppingAction.hh @@ -41,6 +41,7 @@ class RMGSteppingAction : public G4UserSteppingAction { private: + bool fSkipTracking = false; double fDaughterKillLifetime = -1; std::unique_ptr fMessenger; diff --git a/src/RMGSteppingAction.cc b/src/RMGSteppingAction.cc index a6985b359..9a6f7e364 100644 --- a/src/RMGSteppingAction.cc +++ b/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); } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index ef50faeda..a5330a15f 100644 --- a/tests/CMakeLists.txt +++ b/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() diff --git a/tests/bxdecay0/CMakeLists.txt b/tests/bxdecay0/CMakeLists.txt new file mode 100644 index 000000000..c024116fa --- /dev/null +++ b/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) diff --git a/tests/bxdecay0/gdml/geometry.gdml b/tests/bxdecay0/gdml/geometry.gdml new file mode 100644 index 000000000..996890426 --- /dev/null +++ b/tests/bxdecay0/gdml/geometry.gdml @@ -0,0 +1,31 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/bxdecay0/macros/template.mac b/tests/bxdecay0/macros/template.mac new file mode 100644 index 000000000..46363168b --- /dev/null +++ b/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 diff --git a/tests/bxdecay0/macros/template2.mac b/tests/bxdecay0/macros/template2.mac new file mode 100644 index 000000000..79338f39a --- /dev/null +++ b/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 diff --git a/tests/bxdecay0/plot_decay_modes.py b/tests/bxdecay0/plot_decay_modes.py new file mode 100644 index 000000000..570b10f06 --- /dev/null +++ b/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() diff --git a/tests/bxdecay0/run_sims.py b/tests/bxdecay0/run_sims.py new file mode 100644 index 000000000..66ad89fbf --- /dev/null +++ b/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, + }, + )