gamma angular correlation validation

docs/validation/CMakeLists.txt

@@ -111,6 +111,8 @@ set(TESTS_IMAGES
     observables/beta-observables.surf-active-energy.spec.output.png
     observables/beta-observables.surf-max-z.eff.output.png
     observables/beta-observables.surf-max-z.spec.output.png
+    gammacorr/gamma-angular-distribution-27-60.output.png
+    gammacorr/gamma-angular-distribution-81-208.output.png
     processes/innerbrem-gamma-vertex-energy.output.png
     processes/innerbrem-signal.output.png
     processes/innerbrem-difference.output.png)

docs/validation/gammacorr.md

@@ -0,0 +1,23 @@
+# Gamma correlations
+
+This section of the validation report shows the effect of angular correlations
+after emission of multiple gammas after a decay. As an example: The background
+rate of Co-60 depends on both gammas hitting one detector, so the angle between
+the gammas matter.
+
+See [presentation](https://indico.cern.ch/event/372884/contributions/1792361/)
+for more background on the calculation.
+
+```{subfigure} AB
+:subcaptions: above
+
+:::{image} ./_img/gammacorr/gamma-angular-distribution-27-60.output.png
+:width: 300px
+:::
+
+:::{image} ./_img/gammacorr/gamma-angular-distribution-81-208.output.png
+:width: 300px
+:::
+
+Plots of the angle between two emitted gammas after a Co-60 and Tl-208 decay.
+```

docs/validation/index.md

@@ -18,6 +18,7 @@ confinement
 distances
 nist
 observables-ge
+gammacorr
 processes-ib
 ```
 

tests/CMakeLists.txt

@@ -11,6 +11,7 @@ add_subdirectory(basics)
 add_subdirectory(confinement)
 add_subdirectory(detectors)
 add_subdirectory(distances)
+add_subdirectory(gammacorr)
 add_subdirectory(hades)
 add_subdirectory(internals)
 add_subdirectory(nist)

tests/gammacorr/CMakeLists.txt

@@ -0,0 +1,19 @@
+file(
+  GLOB _file_list
+  RELATIVE ${PROJECT_SOURCE_DIR}
+  macros/*.mac gdml/*.gdml *.py)
+
+# copy them to the build area
+foreach(_file ${_file_list})
+  configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY)
+endforeach()
+
+# run python tests using pytest
+add_test(NAME gammacorr/pytest-all COMMAND ${PYTHONPATH} -m pytest -s -vvv .)
+
+set_tests_properties(gammacorr/pytest-all PROPERTIES LABELS "mt;extra" ENVIRONMENT
+                                                     MPLCONFIGDIR=${CMAKE_SOURCE_DIR}/tests)
+
+if(Geant4_VERSION VERSION_LESS "11.3.0")
+  set_tests_properties(gammacorr/pytest-all PROPERTIES DISABLED True)
+endif()

tests/gammacorr/gdml/geometry.gdml

@@ -0,0 +1,25 @@
+<?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>
+		<orb name="world" r="20" lunit="cm"/>
+		<orb name="vacuum" r="15" lunit="cm"/>
+	</solids>
+	<structure>
+		<volume name="vacuum">
+			<materialref ref="G4_Galactic"/>
+			<solidref ref="vacuum"/>
+		</volume>
+		<volume name="world">
+			<materialref ref="G4_Galactic"/>
+			<solidref ref="world"/>
+			<physvol name="vacuum">
+				<volumeref ref="vacuum"/>
+			</physvol>
+		</volume>
+	</structure>
+	<setup name="Default" version="1.0">
+		<world ref="world"/>
+	</setup>
+</gdml>

tests/gammacorr/gdml/geometry.py

@@ -0,0 +1,17 @@
+from __future__ import annotations
+
+import pyg4ometry as pg4
+
+reg = pg4.geant4.Registry()
+
+world_s = pg4.geant4.solid.Orb("world", 20, registry=reg, lunit="cm")
+world_l = pg4.geant4.LogicalVolume(world_s, "G4_Galactic", "world", registry=reg)
+reg.setWorld(world_l)
+
+ge_s = pg4.geant4.solid.Orb("vacuum", 15, registry=reg, lunit="cm")
+ge_l = pg4.geant4.LogicalVolume(ge_s, "G4_Galactic", "vacuum", registry=reg)
+pg4.geant4.PhysicalVolume([0, 0, 0], [0, 0, 0], ge_l, "vacuum", world_l, registry=reg)
+
+w = pg4.gdml.Writer()
+w.addDetector(reg)
+w.write("geometry.gdml")

tests/gammacorr/test_angcorr.py

@@ -0,0 +1,108 @@
+from __future__ import annotations
+
+import awkward as ak
+import hist
+import matplotlib.pyplot as plt
+import numpy as np
+from lgdo import lh5
+from remage import remage_run
+
+macro = """
+/RMG/Output/ActivateOutputScheme Track
+/RMG/Output/ActivateOutputScheme ParticleFilter
+
+/RMG/Processes/EnableGammaAngularCorrelation {angcorr}
+
+/run/initialize
+
+/RMG/Output/ParticleFilter/AddParticle 11
+/RMG/Output/ParticleFilter/AddKillVolume vacuum
+
+/RMG/Generator/Confine UnConfined
+
+/RMG/Generator/Select GPS
+/gps/position     0 0 0
+/gps/particle     ion
+/gps/ion          {Z} {A} 0 {level}
+/gps/energy       0 keV
+/gps/ang/type     iso
+
+/run/beamOn 100000
+"""
+
+
+def simulate(Z, A, level, angcorr):
+    output = f"output-{Z}-{A}-{level}-angcorr-{angcorr}.lh5"
+
+    remage_run(
+        macro.split("\n"),
+        macro_substitutions={"Z": Z, "A": A, "level": level, "angcorr": angcorr},
+        gdml_files="gdml/geometry.gdml",
+        output=output,
+        overwrite_output=True,
+        log_level="summary",
+    )
+
+    return output
+
+
+def test_plot_gammacorr():
+    level = 0
+
+    items = [
+        (27, 60, "$^{60}$Co: 1.17 MeV vs. 1.33 MeV"),
+        (81, 208, "$^{208}$Tl: 0.58 MeV vs. 2.6 MeV"),
+    ]
+
+    for Z, A, title in items:
+        fig, ax = plt.subplots()
+
+        for angcorr in (True, False):
+            remage_output = simulate(Z, A, level, angcorr)
+            # read in track data
+            tracks = lh5.read_as("tracks", remage_output, library="ak")
+
+            # read in dictionary with process ids, to filter later
+            processes = lh5.read("processes", remage_output)
+
+            # group-by event id
+            trk = ak.unflatten(tracks, ak.run_lengths(tracks.evtid))
+
+            # select gamma particles created by a radioactive decay
+            decay_procid = processes[
+                (
+                    "RadioactiveDecay"
+                    if "RadioactiveDecay" in processes
+                    else "Radioactivation"
+                )
+            ].value
+            gtrk = trk[(trk.particle == 22) & (trk.procid == decay_procid)]
+
+            # keep only events with two gammas
+            gtrk = gtrk[ak.num(gtrk) == 2]
+
+            # gtrk has shape (n_events, 2) and fields .px, .py, .pz
+            # Split into the two momenta per event:
+            p1 = gtrk[:, 0]
+            p2 = gtrk[:, 1]
+
+            # Dot product
+            dot = p1.px * p2.px + p1.py * p2.py + p1.pz * p2.pz
+
+            # Norms
+            norm1 = np.sqrt(p1.px**2 + p1.py**2 + p1.pz**2)
+            norm2 = np.sqrt(p2.px**2 + p2.py**2 + p2.pz**2)
+
+            # Cosine of angle
+            cos_theta = dot / (norm1 * norm2)
+
+            h = hist.new.Reg(50, -1, 1).Double().fill(cos_theta)
+
+            h.plot(
+                ax=ax, yerr=False, label=rf"$\gamma$ angular correlations = {angcorr}"
+            )
+
+        ax.set_title(title)
+        ax.set_xlabel(r"$cos(\theta)$")
+        ax.legend()
+        fig.savefig(f"gamma-angular-distribution-{Z}-{A}.output.png")