add function to sample scintillation

src/legendoptics/lar.py

@@ -47,6 +47,9 @@ class ArScintLiftime(NamedTuple):
     singlet: Quantity
     triplet: Quantity
 
+    def as_tuple(self) -> tuple[Quantity, Quantity]:
+        return (self.singlet, self.triplet)
+
 
 def lar_dielectric_constant_bideau_mehu(
     λ: Quantity,
@@ -248,7 +251,9 @@ def lar_lifetimes(
     return ArScintLiftime(singlet=5.95 * u.ns, triplet=triplet)
 
 
-def lar_scintillation_params(flat_top_yield: Quantity = 31250 / u.MeV) -> ScintConfig:
+def lar_scintillation_params(
+    flat_top_yield: Quantity = 31250 / u.MeV,
+) -> ScintConfig:
     r"""Scintillation yield (approx. inverse of the mean energy to produce a UV photon).
 
     This depends on the nature of the impinging particles, the field configuration
@@ -273,9 +278,14 @@ def lar_scintillation_params(flat_top_yield: Quantity = 31250 / u.MeV) -> ScintC
     Excitation ratio:
     For example, for nuclear recoils it should be 0.75
     nominal value for electrons and gammas: 0.23 (WArP data)
+
+    See Also
+    --------
+    .lar_fano_factor
     """
     return ScintConfig(
         flat_top=flat_top_yield,
+        fano_factor=lar_fano_factor(),
         particles=[
             ScintParticle("electron", yield_factor=0.8, exc_ratio=0.23),
             ScintParticle("alpha", yield_factor=0.7, exc_ratio=1),
@@ -411,7 +421,6 @@ def pyg4_lar_attach_scintillation(
     See Also
     --------
     .lar_scintillation_params
-    .lar_fano_factor
     .lar_emission_spectrum
     .lar_lifetimes
     """
@@ -435,9 +444,11 @@ def pyg4_lar_attach_scintillation(
     lar_mat.addConstPropertyPint("SCINTILLATIONTIMECONSTANT1", lifetimes.singlet)
     lar_mat.addConstPropertyPint("SCINTILLATIONTIMECONSTANT2", lifetimes.triplet)
 
+    scint_params = lar_scintillation_params(flat_top_yield)
+
     # the fano factor is the ratio between variance and mean. Geant4 calculates
     # σ = RESOLUTIONSCALE × √mean, so we have to take the root of the fano factor
     # here to keep it consistent.
-    lar_mat.addConstPropertyPint("RESOLUTIONSCALE", np.sqrt(lar_fano_factor()))
+    lar_mat.addConstPropertyPint("RESOLUTIONSCALE", np.sqrt(scint_params.fano_factor))
 
-    pyg4_def_scint_by_particle_type(lar_mat, lar_scintillation_params(flat_top_yield))
+    pyg4_def_scint_by_particle_type(lar_mat, scint_params)

src/legendoptics/pen.py

@@ -196,6 +196,7 @@ def pyg4_pen_attach_scintillation(mat, reg) -> None:
     # setup scintillation response just for electrons.
     scint_config = ScintConfig(
         flat_top=pen_scint_light_yield(),
+        fano_factor=None,
         particles=[
             ScintParticle("electron", yield_factor=1, exc_ratio=None),
         ],

src/legendoptics/pyg4utils.py

@@ -55,6 +55,12 @@ def _def_scint_particle(
 def pyg4_def_scint_by_particle_type(mat, scint_cfg: ScintConfig) -> None:
     """Define a full set of particles for scintillation."""
     for particle in scint_cfg.particles:
+        if not particle.valid_geant_particle():
+            msg = (
+                f"Unknown particle type {particle.name} for ScintillationByParticleType"
+            )
+            raise ValueError(msg)
+
         _def_scint_particle(
             mat,
             particle.name.upper(),

src/legendoptics/utils.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 
-from typing import NamedTuple
+from typing import Literal, NamedTuple, Optional, get_args, get_type_hints
 
 import numpy as np
 import pint
@@ -109,13 +109,89 @@ def __call__(self, pts: Quantity) -> Quantity:
 class ScintParticle(NamedTuple):
     """Configuration for the scintillation yield relative to the flat-top yield."""
 
-    name: str
+    name: Literal["deuteron", "triton", "alpha", "ion", "electron"]
     yield_factor: float
-    exc_ratio: float | None
+    exc_ratio: Optional[float]  # noqa: UP007
+
+    def valid_geant_particle(self) -> bool:
+        return self.name.lower() in get_args(get_type_hints(ScintParticle)["name"])
 
 
 class ScintConfig(NamedTuple):
     """Scintillation yield parameters, depending on the particle types."""
 
     flat_top: Quantity
+    fano_factor: Optional[float]  # noqa: UP007
     particles: list[ScintParticle]
+
+    def get_particle(self, name: str) -> Optional[ScintParticle]:  # noqa: UP007
+        for p in self.particles:
+            if p.name.upper() == name.upper():
+                return p
+        return None
+
+
+def scintillate(
+    scint_config: ScintConfig,
+    time_components: tuple[Quantity, ...],
+    particle: str,
+    edep: Quantity,
+) -> Quantity:
+    """Generates a Poisson/Gauss-distributed number of photons according to the
+    scintillation yield formula, as implemented in Geant4.
+
+    Returns
+    -------
+    array of scintillation time stamps, relative to the point in time of energy
+    deposition.
+    """
+    if not edep.check("[energy]"):
+        msg = "Edep must have dimensionality energy"
+        raise ValueError(msg)
+
+    # validate that we have a valid configuration of
+    part_with_exc_ratio = [1 for p in scint_config.particles if p.exc_ratio is not None]
+    if len(part_with_exc_ratio) not in (0, len(scint_config.particles)):
+        msg = "Not all particles have exc_ratio defined"
+        raise ValueError(msg)
+    has_2_timeconstants = len(part_with_exc_ratio) > 0
+    if len(time_components) != (1 + int(has_2_timeconstants)):
+        msg = "Number of time_components is not as required"
+        raise ValueError(msg)
+
+    # get the particle from its type string.
+    part = scint_config.get_particle(particle)
+    if part is None:
+        msg = "Request for scintillation yield for particle type without correct entry"
+        raise ValueError(msg)
+
+    mean_num_phot = (
+        scint_config.flat_top * part.yield_factor * edep
+    ).to_reduced_units()
+
+    # derive the actual number of photons generated in this step.
+    gen = np.random.default_rng()
+    if mean_num_phot > 10:
+        fano = scint_config.fano_factor if scint_config.fano_factor is not None else 1
+        sigma = np.sqrt(fano * mean_num_phot)
+        num_photons = int(gen.normal(mean_num_phot, sigma) + 0.5)
+    else:
+        num_photons = gen.poisson(mean_num_phot)
+
+    if num_photons <= 0:
+        return np.empty((0,)) * u.ns
+
+    # derive number of photon for all time components.
+    yields = [num_photons]
+    if has_2_timeconstants:
+        yields = [part.exc_ratio, 1 - part.exc_ratio]
+        yields = [int(num_photons * y) for y in yields]
+        yields[-1] = num_photons - sum(yields[0:-1])  # to keep the sum constant.
+
+    # now, calculate the timestamps of each
+    times = []
+    for exc_ratio, scint_t in zip(yields, time_components):
+        num_phot = int(exc_ratio)
+        times.append(-scint_t * np.log(gen.uniform(size=(num_phot,))))
+
+    return np.sort(np.concatenate(times))