From f22278a9e7e7ff74bfd0d15863027b2fc94da49e Mon Sep 17 00:00:00 2001
From: Manuel Huber <info@manuelhu.de>
Date: Mon, 9 Sep 2024 19:08:54 +0200
Subject: [PATCH] add functions to write g4gps spectra

---
 src/legendoptics/lar.py   | 22 ++++++++++++++++++++++
 src/legendoptics/pen.py   | 21 +++++++++++++++++++++
 src/legendoptics/utils.py | 26 ++++++++++++++++++++++++++
 tests/test_g4gps.py       | 19 +++++++++++++++++++
 4 files changed, 88 insertions(+)
 create mode 100644 tests/test_g4gps.py

diff --git a/src/legendoptics/lar.py b/src/legendoptics/lar.py
index 30776f3..ff24e52 100644
--- a/src/legendoptics/lar.py
+++ b/src/legendoptics/lar.py
@@ -38,6 +38,7 @@
     InterpolatingGraph,
     ScintConfig,
     ScintParticle,
+    g4gps_write_emission_spectrum,
     readdatafile,
 )
 
@@ -455,3 +456,24 @@ def pyg4_lar_attach_scintillation(
     lar_mat.addConstPropertyPint("RESOLUTIONSCALE", np.sqrt(lar_fano_factor()))
 
     pyg4_def_scint_by_particle_type(lar_mat, lar_scintillation_params(flat_top_yield))
+
+
+def g4gps_lar_emissions_spectrum(filename: str) -> None:
+    """Write a LAr emission energy spectrum for G4GeneralParticleSource.
+
+    See Also
+    --------
+    .lar_emission_spectrum
+    utils.g4gps_write_emission_spectrum
+    """
+    from legendoptics.pyg4utils import pyg4_sample_λ
+
+    λ_peak = pyg4_sample_λ(116 * u.nm, 141 * u.nm)
+
+    # sample the measured emission spectrum.
+    scint_em = lar_emission_spectrum(λ_peak)
+    # make sure that the scintillation spectrum is zero at the boundaries.
+    scint_em[0] = 0
+    scint_em[-1] = 0
+
+    g4gps_write_emission_spectrum(filename, λ_peak, scint_em)
diff --git a/src/legendoptics/pen.py b/src/legendoptics/pen.py
index 4ff1df8..c98bb0b 100644
--- a/src/legendoptics/pen.py
+++ b/src/legendoptics/pen.py
@@ -23,6 +23,7 @@
     InterpolatingGraph,
     ScintConfig,
     ScintParticle,
+    g4gps_write_emission_spectrum,
     readdatafile,
 )
 
@@ -201,3 +202,23 @@ def pyg4_pen_attach_scintillation(mat, reg) -> None:
         ],
     )
     pyg4_def_scint_by_particle_type(mat, scint_config)
+
+
+def g4gps_pen_emissions_spectrum(filename: str) -> None:
+    """Write a PEN emission energy spectrum for G4GeneralParticleSource.
+
+    See Also
+    --------
+    .pen_wls_emission
+    utils.g4gps_write_emission_spectrum
+    """
+    from legendoptics.pyg4utils import pyg4_sample_λ
+
+    # sample the measured emission spectrum.
+    λ_scint = pyg4_sample_λ(350 * u.nm, 650 * u.nm, 200)
+    scint_em = InterpolatingGraph(*pen_wls_emission(), min_idx=350 * u.nm)(λ_scint)
+    # make sure that the scintillation spectrum is zero at the boundaries.
+    scint_em[0] = 0
+    scint_em[-1] = 0
+
+    g4gps_write_emission_spectrum(filename, λ_scint, scint_em)
diff --git a/src/legendoptics/utils.py b/src/legendoptics/utils.py
index 707b33a..d10f8be 100644
--- a/src/legendoptics/utils.py
+++ b/src/legendoptics/utils.py
@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import logging
 from typing import NamedTuple
 
 import numpy as np
@@ -8,6 +9,7 @@
 from importlib_resources import files
 from pint import Quantity
 
+log = logging.getLogger(__name__)
 u = pint.get_application_registry()
 
 
@@ -119,3 +121,27 @@ class ScintConfig(NamedTuple):
 
     flat_top: Quantity
     particles: list[ScintParticle]
+
+
+def g4gps_write_emission_spectrum(
+    filename: str, λ_peak: Quantity, scint_em: Quantity
+) -> None:
+    """Write a energy spectrum for use with G4GeneralParticleSource
+
+    It can be used like this in a Geant4 macro:
+
+    .. code ::
+
+        /gps/ene/type     Arb
+        /gps/ene/diffspec true
+        /gps/hist/type    arb
+        /gps/hist/file    <filename>
+        /gps/hist/inter   Lin
+    """
+    with u.context("sp"):
+        pointwise = np.array([λ_peak.to("MeV").m, scint_em.m]).T
+
+    if pointwise.shape[0] > 1024:
+        log.warning("G4GeneralParticleSource spectrum can only have 1024 bins.")
+
+    np.savetxt(filename, pointwise)
diff --git a/tests/test_g4gps.py b/tests/test_g4gps.py
new file mode 100644
index 0000000..7788764
--- /dev/null
+++ b/tests/test_g4gps.py
@@ -0,0 +1,19 @@
+"""Test the exporting of G4GPS spectra."""
+
+from __future__ import annotations
+
+import pint
+
+u = pint.get_application_registry()
+
+
+def test_g4gps_lar(tmp_path) -> None:
+    import legendoptics.lar
+
+    legendoptics.lar.g4gps_lar_emissions_spectrum(tmp_path / "lar.csv")
+
+
+def test_g4gps_pen(tmp_path) -> None:
+    import legendoptics.pen
+
+    legendoptics.pen.g4gps_pen_emissions_spectrum(tmp_path / "pen.csv")