add functions to write g4gps spectra

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)

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)

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)

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")