diff --git a/src/legendoptics/utils.py b/src/legendoptics/utils.py index ed1afa3..3d2774e 100644 --- a/src/legendoptics/utils.py +++ b/src/legendoptics/utils.py @@ -169,6 +169,8 @@ def scintillate( scint_config.flat_top * part.yield_factor * edep ).to_reduced_units() + time_components = [t.to(u.microsecond).m for t in time_components] + # derive the actual number of photons generated in this step. gen = np.random.default_rng() if mean_num_phot > 10: @@ -181,17 +183,18 @@ def scintillate( if num_photons <= 0: return np.empty((0,)) * u.ns - # derive number of photon for all time components. + # derive number of photons 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,)))) + # now, calculate the timestamps of each generated photon. + times = np.log(gen.uniform(size=num_photons)) + start = 0 + for num_phot, scint_t in zip(yields, time_components): + times[start : start + num_phot] *= -scint_t + start += num_phot - return np.sort(np.concatenate(times)) + return u.Quantity(times, u.microsecond) diff --git a/tests/test_scintillate.py b/tests/test_scintillate.py new file mode 100644 index 0000000..26c4e79 --- /dev/null +++ b/tests/test_scintillate.py @@ -0,0 +1,24 @@ +"""Test the attaching of all properties to Geant4 materials.""" + +from __future__ import annotations + +import pint + +from legendoptics import lar, utils + +u = pint.get_application_registry() + + +def test_scintillate(): + utils.scintillate( + lar.lar_scintillation_params(), + lar.lar_lifetimes().as_tuple(), + "electron", + 10 * u.keV, + ) + utils.scintillate( + lar.lar_scintillation_params(), + lar.lar_lifetimes().as_tuple(), + "ion", + 10 * u.keV, + )