Add inversegaussian_nIRF from ac/nIRF_model

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 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)
+    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

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=1e-2)