From d320ae2213a47bd3f3c74c8c6e2ebe4985b91f4c Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Wed, 3 Jul 2024 21:46:23 +0200 Subject: [PATCH] add function to sample scintillation --- src/legendoptics/lar.py | 19 ++++++-- src/legendoptics/pen.py | 1 + src/legendoptics/pyg4utils.py | 6 +++ src/legendoptics/utils.py | 86 +++++++++++++++++++++++++++++++++-- 4 files changed, 103 insertions(+), 9 deletions(-) diff --git a/src/legendoptics/lar.py b/src/legendoptics/lar.py index 181b5b1..72e36ed 100644 --- a/src/legendoptics/lar.py +++ b/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) diff --git a/src/legendoptics/pen.py b/src/legendoptics/pen.py index 4ff1df8..98fcc03 100644 --- a/src/legendoptics/pen.py +++ b/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), ], diff --git a/src/legendoptics/pyg4utils.py b/src/legendoptics/pyg4utils.py index 5af681c..e5ff969 100644 --- a/src/legendoptics/pyg4utils.py +++ b/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(), diff --git a/src/legendoptics/utils.py b/src/legendoptics/utils.py index 707b33a..f1a32bc 100644 --- a/src/legendoptics/utils.py +++ b/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 @@ -62,8 +62,8 @@ def __init__( self, idx: Quantity, vals: Quantity, - min_idx: Quantity | None = None, - max_idx: Quantity | None = None, + min_idx: Optional[Quantity] = None, + max_idx: Optional[Quantity] = None, ): # Filter the supplied data points. f = np.full(idx.shape, True) @@ -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] + + 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] particles: list[ScintParticle] + + def get_particle(self, name: str) -> Optional[ScintParticle]: + 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))