From 758b39582b2be350d3f1ee26be3082da5ff2afdb Mon Sep 17 00:00:00 2001
From: Manuel Huber <ManuelHu@users.noreply.github.com>
Date: Mon, 18 Sep 2023 23:40:33 +0200
Subject: [PATCH] lar: emission interpolation in defining function (#36)

* lar: move emission spectrum interpolation to property function

* pyg4units: update supported length units for Geant4 11.1.0

* fix tests calling lar_emission_spectrum

* style: pre-commit fixes

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
---
 docs/source/exts/optics_plot_extension.py |  5 ++++-
 src/legendoptics/lar.py                   | 27 ++++++++++++++---------
 src/legendoptics/pyg4utils.py             |  3 ++-
 tests/test_import.py                      | 17 +++++++-------
 tests/test_plot.py                        |  5 +++--
 5 files changed, 34 insertions(+), 23 deletions(-)

diff --git a/docs/source/exts/optics_plot_extension.py b/docs/source/exts/optics_plot_extension.py
index 00841bf..5710d04 100644
--- a/docs/source/exts/optics_plot_extension.py
+++ b/docs/source/exts/optics_plot_extension.py
@@ -48,7 +48,10 @@ def do_plot(
 
     if "call_x" in options:
         # special case for LAr properties
-        x = np.linspace(112 * u.nm, 650 * u.nm, num=200)
+        lim = [112 * u.nm, 650 * u.nm]
+        if "xlim" in options:
+            lim = [xl * u.nm for xl in options["xlim"]]
+        x = np.linspace(*lim, num=200)
         ys = obj(x)
         # wrap the result in a tuple, if needed
         ys = ys if isinstance(ys, tuple) else (ys,)
diff --git a/src/legendoptics/lar.py b/src/legendoptics/lar.py
index 4c95466..d62ba4c 100644
--- a/src/legendoptics/lar.py
+++ b/src/legendoptics/lar.py
@@ -132,12 +132,21 @@ def lar_refractive_index(λ: Quantity, method: str = "cern2020") -> Quantity:
     return np.sqrt(lar_dielectric_constant(λ, method))
 
 
-def lar_emission_spectrum() -> tuple[Quantity, Quantity]:
+def lar_emission_spectrum(λ: Quantity) -> Quantity:
     """Return the LAr emission spectrum, adapted from [Heindl2010]_.
 
-    .. optics-plot::
+    .. optics-plot:: {'call_x': True, 'xlim': [116, 141]}
     """
-    return readdatafile("lar_emission_heindl2010.dat")
+    heindl = readdatafile("lar_emission_heindl2010.dat")
+
+    # sample the measured emission spectrum and avoid the fluctuations below 115 nm.
+    scint_em = InterpolatingGraph(
+        *heindl,
+        min_idx=115 * u.nm,
+        max_idx=150 * u.nm,
+    )(λ)
+
+    return scint_em
 
 
 def lar_fano_factor() -> float:
@@ -256,11 +265,11 @@ def lar_scintillation_params(flat_top_yield: Quantity = 31250 / u.MeV) -> ScintC
     relative to the one of flat top particles is:
 
     .. math::
-        Y_\texrm{e} &= 0.8 Y
+        Y_\textrm{e} &= 0.8 Y
 
-        Y_\texrm{alpha} &= 0.7 Y
+        Y_\textrm{alpha} &= 0.7 Y
 
-        Y_\texrm{recoils} &= 0.2\textrm{--}0.4
+        Y_\textrm{recoils} &= 0.2\textrm{--}0.4
 
     Excitation ratio:
     For example, for nuclear recoils it should be 0.75
@@ -412,11 +421,7 @@ def pyg4_lar_attach_scintillation(
     λ_peak = pyg4_sample_λ(116 * u.nm, 141 * u.nm)
 
     # sample the measured emission spectrum.
-    scint_em = InterpolatingGraph(
-        *lar_emission_spectrum(),
-        min_idx=115 * u.nm,
-        max_idx=150 * u.nm,
-    )(λ_peak)
+    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
diff --git a/src/legendoptics/pyg4utils.py b/src/legendoptics/pyg4utils.py
index 14a74c6..248449a 100644
--- a/src/legendoptics/pyg4utils.py
+++ b/src/legendoptics/pyg4utils.py
@@ -96,7 +96,8 @@ def _val_pint_to_gdml(v):
         v = v.m_as(base_unit)
         return unit, v
 
-    length_u = ["m", "cm", "mm", "um"]
+    # Only as of Geant4 11.1.0, `um` and `nm` are supported.
+    length_u = ["km", "m", "cm", "mm", "um", "nm"]
     dimless_props = [
         "RINDEX",
         "WLSCOMPONENT",
diff --git a/tests/test_import.py b/tests/test_import.py
index c573dd0..c9282d4 100644
--- a/tests/test_import.py
+++ b/tests/test_import.py
@@ -17,14 +17,15 @@ def test_import():
     import legendoptics.tpb
 
     legendoptics.lar.lar_fano_factor()
-    legendoptics.lar.lar_emission_spectrum()
-    wvl = np.arange(110, 400, 20) * u.nm
-    legendoptics.lar.lar_dielectric_constant_bideau_mehu(wvl)
-    legendoptics.lar.lar_dielectric_constant_cern2020(wvl)
-    legendoptics.lar.lar_dielectric_constant(wvl)
-    legendoptics.lar.lar_refractive_index(wvl)
-    legendoptics.lar.lar_rayleigh(wvl, 87 * u.K)
-    legendoptics.lar.lar_abs_length(wvl)
+    λ_peak = np.arange(116, 141, 20) * u.nm
+    legendoptics.lar.lar_emission_spectrum(λ_peak)
+    λ = np.arange(110, 400, 20) * u.nm
+    legendoptics.lar.lar_dielectric_constant_bideau_mehu(λ)
+    legendoptics.lar.lar_dielectric_constant_cern2020(λ)
+    legendoptics.lar.lar_dielectric_constant(λ)
+    legendoptics.lar.lar_refractive_index(λ)
+    legendoptics.lar.lar_rayleigh(λ, 87 * u.K)
+    legendoptics.lar.lar_abs_length(λ)
     legendoptics.lar.lar_peak_attenuation_length()
     legendoptics.lar.lar_lifetimes()
     legendoptics.lar.lar_scintillation_params()
diff --git a/tests/test_plot.py b/tests/test_plot.py
index 28b86fb..4fddf93 100644
--- a/tests/test_plot.py
+++ b/tests/test_plot.py
@@ -2,7 +2,8 @@
 import numpy as np
 import pint
 
-from legendoptics.lar import lar_emission_spectrum, lar_refractive_index
+from legendoptics.copper import copper_reflectivity
+from legendoptics.lar import lar_refractive_index
 from legendoptics.plot import plot_continuous_prop, plot_discrete_prop
 
 u = pint.get_application_registry()
@@ -15,4 +16,4 @@ def test_plot_continuous_prop():
 
 def test_plot_discrete_prop():
     fig, ax = plt.subplots()
-    plot_discrete_prop(ax, lar_emission_spectrum)
+    plot_discrete_prop(ax, copper_reflectivity)