From e8ae7358699ed06d5083004c917fd99fdfc5d8dd Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 8 May 2026 17:47:43 +0000 Subject: [PATCH 1/2] Add inversegaussian_nIRF function and test from ac/nIRF_model branch Agent-Logs-Url: https://github.com/aidancrilly/NeSST/sessions/e59eb958-e216-4764-8208-348bbfedcd5d Co-authored-by: aidancrilly <20111148+aidancrilly@users.noreply.github.com> --- src/NeSST/time_of_flight.py | 100 +++++++++++++++++++++++++++++++++++ tests/test_time_of_flight.py | 21 ++++++++ 2 files changed, 121 insertions(+) diff --git a/src/NeSST/time_of_flight.py b/src/NeSST/time_of_flight.py index 55542d2..552db6a 100644 --- a/src/NeSST/time_of_flight.py +++ b/src/NeSST/time_of_flight.py @@ -227,6 +227,106 @@ def base(t_detected, En): return base +def inversegaussian_nIRF( + scint_thickness, ni_scin=8.79e28, CH_ratio=8 / 18, E_lower=0.05e6, E_upper=25.0e6, NE_interp=1000 +): + if scint_thickness != 10e-2 or ni_scin != 8.79e28: + raise NotImplementedError("Current inverse gaussian nIRF fit coefficients for 1e29 1/m^3 and 10cm only!") + + E_range = np.linspace(E_lower, E_upper, NE_interp) + sig_H = mat_dict["H"].sigma_tot(E_range) + sig_C = mat_dict["C12"].sigma_tot(E_range) + + def sig_CH(E): + sig_barns = CH_ratio * np.interp(E, E_range, sig_C) + (1 - CH_ratio) * np.interp(E, E_range, sig_H) + return sig_barns * 1e-28 + + def IGtail_nIRF_bestfit_coeffs(En): + """ + Best fit coefficients as a function of neutron energy + + NB these are specific to a scintillator geometry + For different geometry you must re-perform MCNP/Geant/Scintillator1DMonteCarlo sims + + Analysis performed by A. Crilly, 2025 + """ + E_MeV = En / 1e6 + f = 0.46 * np.ones_like(E_MeV) + A = 0.25 * np.ones_like(E_MeV) + + mu_inverse_ns = [ + 0.45918122858399824, + 0.6557307634639333, + 0.8333427566991481, + 0.8986859864960112, + 1.146202470556479, + 1.2341761264319626, + 1.2686641782526762, + 2.8986506749332044, + 1.4799758893943886, + 2.3823877852359687, + 3.84910816578929, + 3.1966724468523355, + 4.284122083989886, + 3.452849096581845, + 4.838842934652534, + 9.999999999999998, + ] + + lamb_inverse_ns = [ + 0.7552223497006166, + 1.0369636485550817, + 1.1703584280444446, + 1.4085037422394058, + 1.445604993812616, + 1.464511388605816, + 1.6795667261078404, + 1.3844968569428273, + 1.7564016389637123, + 1.6801129259460668, + 1.6345197013154014, + 1.7884346810613823, + 1.7955059618288889, + 1.8762280752322453, + 1.9247350785402442, + 1.8572104736139436, + ] + Egrid = 1.0 + np.arange(len(mu_inverse_ns)) + + mu = np.interp(E_MeV, Egrid, mu_inverse_ns) * 1e9 + lamb = np.interp(E_MeV, Egrid, lamb_inverse_ns) * 1e9 + return f, A, mu, lamb + + def base(t_detected, En): + vn = Ekin_2_beta(En, Mn) * c + t_transit = scint_thickness / vn + + tt_d, tt_a = np.meshgrid(t_detected, t_detected, indexing="ij") + _, tt_t = np.meshgrid(t_detected, t_transit, indexing="ij") + + f, A, mu, lamb = IGtail_nIRF_bestfit_coeffs(En) + + top_hat_mat = np.eye(t_detected.size) + np.heaviside(tt_d - tt_a, 0.0) - np.heaviside(tt_d - (tt_a + tt_t), 1.0) + exp_E_arg = f * vn * ni_scin * sig_CH(En) + main_response = np.exp(-(tt_d - tt_a) * exp_E_arg[None, :]) * top_hat_mat + + t_shift = tt_d - (tt_a + tt_t) + tail_hat = np.heaviside(t_shift, 0.5) + t_shift[t_shift < 0.0] = 0.0 + prefactor = lamb / mu + t_coeff = 2 * mu**2 / lamb + exp_arg = prefactor[None, :] * (1 - np.sqrt(1 + t_coeff[None, :] * t_shift)) + tail_response = tail_hat * A[None, :] * np.exp(exp_arg) + + R = main_response + tail_response + + row_sum = R.sum(axis=1, keepdims=True) + row_sum[row_sum == 0] = 1 # avoid div-by-zero + return R / row_sum + + return base + + def decaying_gaussian_kernel(FWHM, tau, shift_sigma=2.0): """Single exponential tail multiplied by Gaussian.""" sig = FWHM / 2.355 diff --git a/tests/test_time_of_flight.py b/tests/test_time_of_flight.py index ce86beb..471981c 100644 --- a/tests/test_time_of_flight.py +++ b/tests/test_time_of_flight.py @@ -49,3 +49,24 @@ def test_IRF_normalisation(kernel_fn): integral_signal = np.trapezoid(y=signal, x=normt_det) assert np.isclose(integral_signal, 1.0, rtol=1e-3) + + +def test_inversegaussian_nIRF_normalisation(): + """ + For uniform sensitivity, the inversegaussian_nIRF response matrix rows should sum to 1 + (i.e. the response is normalised). + """ + nToF_distance = 20.0 # m + flat_sens = nst.time_of_flight.get_unity_sensitivity() + total_IRF = nst.time_of_flight.inversegaussian_nIRF(scint_thickness=10e-2) + + E_ntof = np.linspace(1.0e6, 10.0e6, 500) + DDmean, _, DDvar = nst.DDprimspecmoments(Tion=5.0e3) + dNdE = nst.QBrysk(E_ntof, DDmean, DDvar) + + test_20m_nToF = nst.time_of_flight.nToF(nToF_distance, flat_sens, total_IRF) + t_det, normt_det, signal = test_20m_nToF.get_signal(E_ntof, dNdE) + + integral_signal = np.trapezoid(y=signal, x=normt_det) + + assert np.isclose(integral_signal, 1.0, rtol=0.15) From f167a485176744f7539c3d3174800129cfd4d237 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 8 May 2026 17:48:51 +0000 Subject: [PATCH 2/2] Fix error message and tighten normalisation test tolerance Agent-Logs-Url: https://github.com/aidancrilly/NeSST/sessions/e59eb958-e216-4764-8208-348bbfedcd5d Co-authored-by: aidancrilly <20111148+aidancrilly@users.noreply.github.com> --- src/NeSST/time_of_flight.py | 2 +- tests/test_time_of_flight.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NeSST/time_of_flight.py b/src/NeSST/time_of_flight.py index 552db6a..4dc2010 100644 --- a/src/NeSST/time_of_flight.py +++ b/src/NeSST/time_of_flight.py @@ -231,7 +231,7 @@ def inversegaussian_nIRF( scint_thickness, ni_scin=8.79e28, CH_ratio=8 / 18, E_lower=0.05e6, E_upper=25.0e6, NE_interp=1000 ): if scint_thickness != 10e-2 or ni_scin != 8.79e28: - raise NotImplementedError("Current inverse gaussian nIRF fit coefficients for 1e29 1/m^3 and 10cm only!") + raise NotImplementedError("Current inverse gaussian nIRF fit coefficients for 8.79e28 1/m^3 and 10cm only!") E_range = np.linspace(E_lower, E_upper, NE_interp) sig_H = mat_dict["H"].sigma_tot(E_range) diff --git a/tests/test_time_of_flight.py b/tests/test_time_of_flight.py index 471981c..062c933 100644 --- a/tests/test_time_of_flight.py +++ b/tests/test_time_of_flight.py @@ -69,4 +69,4 @@ def test_inversegaussian_nIRF_normalisation(): integral_signal = np.trapezoid(y=signal, x=normt_det) - assert np.isclose(integral_signal, 1.0, rtol=0.15) + assert np.isclose(integral_signal, 1.0, rtol=1e-2)