pargen: calibration routines fixes (#583)

* split dt cuts out of aoe cuts

* fix name for LQ_Corrected

* binning fix

* lq fail fix

* finite check on weights

* finite check on weights

* lq fix and pc

* bugfixes, add gauss linear guess etc. add pval checks, add csqr output, add cal fit with errors plot

* add no timecorr mode, bugfix for dt corr plot

* add calibration errs

* improve aoe fitting a la energy_cal

src/pygama/math/functions/hpge_peak.py

@@ -210,7 +210,7 @@ def hpge_get_mode(self, pars: np.ndarray, cov: np.ndarray = None) -> tuple:
         # We need to ditch the x_lo and x_hi columns and rows
         if cov is not None:
             cov = np.array(cov)
-            dropped_cov = cov[:, 2:][2:, :]
+            dropped_cov = cov[2:, 2:]
 
             return hpge_peak_mode(
                 pars[mu_idx],
@@ -228,7 +228,7 @@ def hpge_get_mode(self, pars: np.ndarray, cov: np.ndarray = None) -> tuple:
         if cov is None:
             return pars[mu_idx]
         else:
-            return np.sqrt(cov[mu_idx][mu_idx])
+            return pars[mu_idx], np.sqrt(cov[mu_idx][mu_idx])
 
 
 # hpge_peak.get_fwhm = hpge_get_fwhm

src/pygama/pargen/AoE_cal.py

@@ -343,6 +343,10 @@ def unbinned_aoe_fit(
 
     bounds = aoe_peak_bounds(pdf, x0)
 
+    gof_hist, gof_bins, gof_var = pgh.get_hist(
+        aoe[(aoe < fmax) & (aoe > fmin)], bins=100, range=(fmin, fmax)
+    )
+
     # Full fit using Gaussian signal with Gaussian tail background
     c = cost.ExtendedUnbinnedNLL(aoe[(aoe < fmax) & (aoe > fmin)], pdf.pdf_ext)
     m = Minuit(c, *x0)
@@ -354,12 +358,72 @@ def unbinned_aoe_fit(
     m.migrad()
     m.hesse()
 
-    if np.isnan(m.errors).all():
-        try:
-            m.simplex.migrad()
-            m.hesse()
-        except Exception:
-            return return_nans(pdf)
+    valid1 = (
+        m.valid
+        & (~np.isnan(np.array(m.errors)[mask]).any())
+        & (~(np.array(m.errors)[mask] == 0).all())
+    )
+
+    cs = pgf.goodness_of_fit(
+        gof_hist,
+        gof_bins,
+        gof_var,
+        pdf.get_pdf,
+        m.values,
+        method="Pearson",
+        scale_bins=True,
+    )
+    cs = (cs[0], cs[1] + len(np.where(mask)[0]))
+
+    fit1 = (m.values, m.errors, m.covariance, cs, pdf, mask, valid1, m)
+
+    m2 = Minuit(c, *x0)
+    for arg, val in bounds.items():
+        m2.limits[arg] = val
+    fixed, mask = aoe_peak_fixed(pdf)
+    for fix in fixed:
+        m2.fixed[fix] = True
+    m2.simplex().migrad()
+    m2.hesse()
+
+    valid2 = (
+        m2.valid
+        & (~np.isnan(np.array(m2.errors)[mask]).any())
+        & (~(np.array(m2.errors)[mask] == 0).all())
+    )
+    cs2 = pgf.goodness_of_fit(
+        gof_hist,
+        gof_bins,
+        gof_var,
+        pdf.get_pdf,
+        m2.values,
+        method="Pearson",
+        scale_bins=True,
+    )
+    cs2 = (cs2[0], cs2[1] + len(np.where(mask)[0]))
+
+    fit2 = (m2.values, m2.errors, m2.covariance, cs2, pdf, mask, valid2, m2)
+
+    frac_errors1 = np.sum(np.abs(np.array(m.errors)[mask] / np.array(m.values)[mask]))
+    frac_errors2 = np.sum(np.abs(np.array(m2.errors)[mask] / np.array(m2.values)[mask]))
+
+    if valid2 is False:
+        fit = fit1
+    elif valid1 is False:
+        fit = fit2
+    elif cs[0] * 1.05 < cs2[0]:
+        fit = fit1
+
+    elif cs2[0] * 1.05 < cs[0]:
+        fit = fit2
+
+    elif frac_errors1 < frac_errors2:
+        fit = fit1
+
+    elif frac_errors1 > frac_errors2:
+        fit = fit2
+    else:
+        fit = fit1
 
     if display > 1:
         aoe = aoe[(aoe < fmax) & (aoe > fmin)]
@@ -374,9 +438,9 @@ def unbinned_aoe_fit(
         xs = np.linspace(fmin, fmax, 1000)
         counts, bins, bars = plt.hist(aoe, bins=nbins, histtype="step", label="Data")
         dx = np.diff(bins)
-        plt.plot(xs, pdf.get_pdf(xs, *m.values) * dx[0], label="Full fit")
+        plt.plot(xs, pdf.get_pdf(xs, *fit[0]) * dx[0], label="Full fit")
         pdf.components = True
-        sig, bkg = pdf.get_pdf(xs, *m.values)
+        sig, bkg = pdf.get_pdf(xs, *fit[0])
         pdf.components = False
         plt.plot(xs, sig * dx[0], label="Signal")
         plt.plot(xs, bkg * dx[0], label="Background")
@@ -391,18 +455,18 @@ def unbinned_aoe_fit(
 
         plt.figure()
         bin_centers = (bins[1:] + bins[:-1]) / 2
-        res = (pdf.pdf(bin_centers, *m.values) * dx[0]) - counts
+        res = (pdf.pdf(bin_centers, *fit[0]) * dx[0]) - counts
         plt.plot(
             bin_centers,
             [re / count if count != 0 else re for re, count in zip(res, counts)],
             label="Normalised Residuals",
         )
         plt.legend(loc="upper left")
         plt.show()
-        return m.values, m.errors, m.covariance
+        return fit[0], fit[1], fit[2]
 
     else:
-        return m.values, m.errors, m.covariance
+        return fit[0], fit[1], fit[2]
 
 
 def fit_time_means(tstamps, means, sigmas):
@@ -957,6 +1021,15 @@ def time_correction(
                             )
                         }
 
+                    elif mode == "none":
+                        time_dict = {
+                            time: mean
+                            for time, mean in zip(
+                                np.array(self.timecorr_df.index),
+                                np.nanmean(self.timecorr_df["mean"]),
+                            )
+                        }
+
                     else:
                         raise ValueError("unknown mode")
 
@@ -1864,6 +1937,7 @@ def plot_aoe_mean_time(
             aoe_class.timecorr_df["mean"],
             yerr=aoe_class.timecorr_df["mean_err"],
             linestyle=" ",
+            marker="x",
         )
 
         grouped_means = [
@@ -1923,6 +1997,7 @@ def plot_aoe_res_time(
             aoe_class.timecorr_df["res"],
             yerr=aoe_class.timecorr_df["res_err"],
             linestyle=" ",
+            marker="x",
         )
     except Exception:
         pass
@@ -2005,7 +2080,7 @@ def drifttime_corr_plot(
 
         hist, bins, var = pgh.get_hist(
             final_df[aoe_class.dt_param],
-            dx=10,
+            dx=32,
             range=(
                 np.nanmin(final_df[aoe_class.dt_param]),
                 np.nanmax(final_df[aoe_class.dt_param]),