From 4b8781ed6c50cfb2686fa9fea1ff2ee3f2a5b2e8 Mon Sep 17 00:00:00 2001
From: ahmedk <karim.ahmed@xfel.eu>
Date: Thu, 5 Sep 2024 13:50:39 +0200
Subject: [PATCH] use interpolate and safe division
---
...ungfrau_gain_Spectra_Fit_Summary_NBC.ipynb | 7 +++---
src/cal_tools/jungfrau/jungfrau_ff.py | 23 +++++++++++++------
2 files changed, 20 insertions(+), 10 deletions(-)
diff --git a/notebooks/Jungfrau/Jungfrau_gain_Spectra_Fit_Summary_NBC.ipynb b/notebooks/Jungfrau/Jungfrau_gain_Spectra_Fit_Summary_NBC.ipynb
index 3df1df4fe..06be513fe 100644
--- a/notebooks/Jungfrau/Jungfrau_gain_Spectra_Fit_Summary_NBC.ipynb
+++ b/notebooks/Jungfrau/Jungfrau_gain_Spectra_Fit_Summary_NBC.ipynb
@@ -186,7 +186,7 @@
"metadata": {},
"outputs": [],
"source": [
- "display(Markdown(f\"## Display fitting results using {fit_func} model\"))\n",
+ "display(Markdown(f\"## Display fitting results using {fit_func}\"))\n",
"\n",
"def plot_stacked_heatmap(geom, stacked_constants, title):\n",
" _, ax = plt.subplots(figsize=(8, 6))\n",
@@ -196,6 +196,7 @@
" ax=ax,\n",
" vmin=vmin,\n",
" vmax=vmax,\n",
+ " interpolation=\"None\",\n",
" colorbar={'shrink': 1, 'pad': 0.01, 'label': title},\n",
" )\n",
" ax.set_title(title, size=12)\n",
@@ -226,7 +227,8 @@
" data[..., i],\n",
" cmap='viridis',\n",
" vmin=vmin,\n",
- " vmax=vmax\n",
+ " vmax=vmax,\n",
+ " interpolation=\"None\",\n",
" )\n",
" ax.set_title(f'Cell {i}', size=10)\n",
" plt.colorbar(im, ax=ax)\n",
@@ -323,7 +325,6 @@
" av_modules.append(da)\n",
" except Exception as e:\n",
" warning(f\"Error while loading Fitting file {file_path}: {e}\")\n",
- " # Help to avoid plotting missing module.\n",
"\n",
" # Plot stacked heatmap\n",
" plot_stacked_heatmap(geom, stacked_constants, title)\n",
diff --git a/src/cal_tools/jungfrau/jungfrau_ff.py b/src/cal_tools/jungfrau/jungfrau_ff.py
index b451d5ad3..8e4361964 100644
--- a/src/cal_tools/jungfrau/jungfrau_ff.py
+++ b/src/cal_tools/jungfrau/jungfrau_ff.py
@@ -227,7 +227,7 @@ def rebin_histo(hist, bin_centers, rebin_factor):
def fit_histogram(
bin_centers,
- _fit_func,
+ fit_func,
n_sigma,
rebin,
ratio,
@@ -238,8 +238,8 @@ def fit_histogram(
Wrap around function for fitting of histogram
Args:
- x (list, float): - bin centers along x
- _fit_func (string): which function to use for fit
+ bin_centers (list, float): - bin centers along x
+ fit_func (string): which function to use for fit
- CHARGE_SHARING: single peak with charge sharing tail
- CHARGE_SHARING_2: sum of two CHARGE_SHARING
- GAUSS: gaussian function
@@ -261,7 +261,7 @@ def fit_histogram(
GAUSS=fit_gauss,
)
initial_sigma = 15.
- fit_func = _funcs[_fit_func]
+ fit_func = _funcs[fit_func]
n_cells, n_rows, n_cols = histo.shape[1:]
_init_array = np.zeros((n_cells, n_rows, n_cols), dtype=np.float32)
@@ -458,9 +458,12 @@ def fit_charge_sharing(
yerr_fit = yerr[i1:i2]
def cost_function(ped, q, sigma, alpha, norm):
+ # Add a small constant to yerr_fit to avoid division by zero
+ epsilon = 1e-10 # Small constant to prevent division by zero
+ safe_yerr = yerr_fit + epsilon
return np.sum(
((y_fit - charge_sharing(
- x_fit, ped, q, sigma, alpha, norm)) / yerr_fit)**2)
+ x_fit, ped, q, sigma, alpha, norm)) / safe_yerr)**2)
m = Minuit(
cost_function,
@@ -523,11 +526,14 @@ def fit_double_charge_sharing(x, y, yerr, initial_sigma, n_sigma, ratio):
def cost_function(
ped, q_a, sigma_a, alpha, norm_a, q_b, sigma_b, norm_b
):
+ # Add a small constant to yerr_fit to avoid division by zero
+ epsilon = 1e-10 # Small constant to prevent division by zero
+ safe_yerr = yerr_fit + epsilon
return np.sum((
(y_fit - double_peak(
x_fit, ped, q_a, sigma_a, alpha,
norm_a, q_b, sigma_b, norm_b
- )) / yerr_fit)**2)
+ )) / safe_yerr)**2)
limit_alpha = (0.01, 1.)
limit_norm_a = (np.sqrt(0.1 * norm), 2. * norm)
@@ -604,9 +610,12 @@ def fit_gauss(x, y, yerr, initial_sigma, n_sigma, ratio):
yerr_fit = yerr[i1:i2]
def cost_function(amp, mean, sigma):
+ # Add a small constant to yerr_fit to avoid division by zero
+ epsilon = 1e-10 # Small constant to prevent division by zero
+ safe_yerr = yerr_fit + epsilon
return np.sum(((
y_fit - gauss(
- x_fit, amp, mean, sigma)) / yerr_fit)**2)
+ x_fit, amp, mean, sigma)) / safe_yerr)**2)
m = Minuit(
cost_function,
--
GitLab