diff --git a/docs/validation/CMakeLists.txt b/docs/validation/CMakeLists.txt index 08d1bc04f..17a2071d3 100644 --- a/docs/validation/CMakeLists.txt +++ b/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) diff --git a/docs/validation/gammacorr.md b/docs/validation/gammacorr.md new file mode 100644 index 000000000..6768ef33a --- /dev/null +++ b/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. +``` diff --git a/docs/validation/index.md b/docs/validation/index.md index 68ae5c283..7bbcbec4a 100644 --- a/docs/validation/index.md +++ b/docs/validation/index.md @@ -18,6 +18,7 @@ confinement distances nist observables-ge +gammacorr processes-ib ``` diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 0e8b1011b..9814c1c0a 100644 --- a/tests/CMakeLists.txt +++ b/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) diff --git a/tests/gammacorr/CMakeLists.txt b/tests/gammacorr/CMakeLists.txt new file mode 100644 index 000000000..ca03b9265 --- /dev/null +++ b/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() diff --git a/tests/gammacorr/gdml/geometry.gdml b/tests/gammacorr/gdml/geometry.gdml new file mode 100644 index 000000000..e6b3c8679 --- /dev/null +++ b/tests/gammacorr/gdml/geometry.gdml @@ -0,0 +1,25 @@ + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/gammacorr/gdml/geometry.py b/tests/gammacorr/gdml/geometry.py new file mode 100644 index 000000000..ce9fb4ca7 --- /dev/null +++ b/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") diff --git a/tests/gammacorr/test_angcorr.py b/tests/gammacorr/test_angcorr.py new file mode 100644 index 000000000..5280cbb68 --- /dev/null +++ b/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")