From f4dc7294c1d456f2d9677bfd9d6932c68772361c Mon Sep 17 00:00:00 2001 From: ggmarshall <72088559+ggmarshall@users.noreply.github.com> Date: Fri, 10 May 2024 15:34:29 +0100 Subject: [PATCH] Pargen bug fixes (#582) --- src/pygama/pargen/energy_cal.py | 19 +++++++++++++------ src/pygama/pargen/lq_cal.py | 3 +++ 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/pygama/pargen/energy_cal.py b/src/pygama/pargen/energy_cal.py index 3959dc781..bdac157ca 100644 --- a/src/pygama/pargen/energy_cal.py +++ b/src/pygama/pargen/energy_cal.py @@ -477,8 +477,8 @@ def hpge_cal_energy_peak_tops( (Polynomial(self.pars) - i).roots() for i in (peaks_kev[0] * 0.9, peaks_kev[-1] * 1.1) ) - euc_min = euc_min[0] - euc_max = euc_max[0] + euc_min = np.nanmin(euc_min) + euc_max = np.nanmax(euc_max) if euc_min < 0: euc_min = 0 @@ -831,8 +831,8 @@ def hpge_fit_energy_peaks( (Polynomial(self.pars) - i).roots() for i in (peaks_kev[0] * 0.9, peaks_kev[-1] * 1.1) ) - euc_min = euc_min[0] - euc_max = euc_max[0] + euc_min = np.nanmin(euc_min) + euc_max = np.nanmax(euc_max) if euc_min < 0: euc_min = 0 if euc_max > np.nanmax(e_uncal) * 1.1: @@ -2484,8 +2484,13 @@ def hpge_fit_energy_cal_func( for n in range(len(energy_scale_pars) - 1): d_mu_d_es += energy_scale_pars[n + 1] * mus ** (n + 1) e_weights = np.sqrt(d_mu_d_es * mu_vars) + mask = np.isfinite(e_weights) poly_pars = ( - Polynomial.fit(mus, energies_kev, deg=deg, w=1 / e_weights).convert().coef + Polynomial.fit( + mus[mask], energies_kev[mask], deg=deg, w=1 / e_weights[mask] + ) + .convert() + .coef ) if fixed is not None: for idx, val in fixed.items(): @@ -2493,7 +2498,9 @@ def hpge_fit_energy_cal_func( pass else: poly_pars[idx] = val - c = cost.LeastSquares(mus, energies_kev, e_weights, poly_wrapper) + c = cost.LeastSquares( + mus[mask], energies_kev[mask], e_weights[mask], poly_wrapper + ) m = Minuit(c, *poly_pars) if fixed is not None: for idx in list(fixed): diff --git a/src/pygama/pargen/lq_cal.py b/src/pygama/pargen/lq_cal.py index bf6f0daea..669330b2b 100644 --- a/src/pygama/pargen/lq_cal.py +++ b/src/pygama/pargen/lq_cal.py @@ -513,6 +513,9 @@ def get_cut_lq_dep(self, df: pd.DataFrame(), lq_param: str, cal_energy_param: st raise (e) log.error("LQ cut determination failed") self.cut_val = np.nan + c = cost.UnbinnedNLL(np.array([0]), gaussian.pdf) + m = Minuit(c, np.full(2, np.nan)) + self.cut_fit_pars = pars = m.values self.update_cal_dicts( {