diff --git a/notebooks/Jungfrau/Jungfrau_Create_Fit_Spectra_Histos_NBC.ipynb b/notebooks/Jungfrau/Jungfrau_Create_Fit_Spectra_Histos_NBC.ipynb
index 11765a8d7102e7cbd2c795fbbd31f3cb99312dce..f59eceda1c14960f0243ee7c18417431ec8ba09b 100644
--- a/notebooks/Jungfrau/Jungfrau_Create_Fit_Spectra_Histos_NBC.ipynb
+++ b/notebooks/Jungfrau/Jungfrau_Create_Fit_Spectra_Histos_NBC.ipynb
@@ -408,7 +408,7 @@
" # Allocate shared arrays for corrected data. Used in `correct_train()`\n",
" data_corr = context.alloc(shape=adc.shape, dtype=np.float32)\n",
" context.map(correct_train, adc)\n",
- " step_timer.done_step(\"correct_train\")\n",
+ " step_timer.done_step(\"correcting trains per chunk\", print_step=False)\n",
"\n",
" step_timer.start()\n",
" chunks = jungfrau_ff.chunk_multi(data_corr, block_size)\n",
@@ -425,7 +425,7 @@
" icol = i%n_blocks_col * block_size[1]\n",
" h_spectra[da][..., irow:irow+block_size[0], icol:icol+block_size[1]] += h_chk\n",
"\n",
- " step_timer.done_step(\"Histogram created\")"
+ " step_timer.done_step(\"Histogram created per chunk\", print_step=False)"
]
},
{
diff --git a/src/cal_tools/jungfrau/jungfrau_ff.py b/src/cal_tools/jungfrau/jungfrau_ff.py
index f1c421efe2e55f2baf43572d6b3c968e7fceefc5..b451d5ad3d949aa402ec4d459d059ed8f916bb62 100644
--- a/src/cal_tools/jungfrau/jungfrau_ff.py
+++ b/src/cal_tools/jungfrau/jungfrau_ff.py
@@ -3,25 +3,28 @@ from scipy.special import erf
from iminuit import Minuit
-def histogram_errors(x):
+def histogram_errors(bin_counts):
"""
Calculate the errors for histogram bins.
+ This function computes the error for each bin in a histogram,
+ treating positive and negative bin counts differently.
+
Args:
x (np.ndarray): Input array of histogram bin counts.
Returns:
np.ndarray: Array of calculated errors.
"""
- out = np.ones_like(x)
+ out = np.ones_like(bin_counts)
# Handle positive values
- positive_mask = x > 0
- out[positive_mask] = np.sqrt(x[positive_mask])
+ positive_mask = bin_counts > 0
+ out[positive_mask] = np.sqrt(bin_counts[positive_mask])
# Handle negative values
- negative_mask = x < 0
- out[negative_mask] = np.sqrt(-x[negative_mask])
+ negative_mask = bin_counts < 0
+ out[negative_mask] = np.sqrt(-bin_counts[negative_mask])
return out
@@ -73,37 +76,41 @@ def chunk_multi(data, block_size):
return chunks
-def fill_histogram(imgs, h_bins):
+def fill_histogram(image_data, histogram_bins):
"""
- Fills an histogram with shape
- (n_bins-1, memory_cells, n_rows, n_cols)
- from input images.
+ Fill a histogram from input images.
+ This function creates a histogram with shape
+ (n_bins-1, memory_cells, n_rows, n_cols) from the input image data.
Args:
- h_bins (list, float): the binning of the x-axis
- imgs (np.ndarray): image data to histogram
+ histogram_bins (list, float): the binning of the x-axis
+ image_data (np.ndarray): image data to histogram
(trains, memory_cells, row, col).
Returns: histogram bin counts, bin edges.
+
+ Raises
+ TypeError: If imgs is not a numpy ndarray.
+ ValueError: If imgs is not a 4D array.
"""
- if not isinstance(imgs, np.ndarray):
+ if not isinstance(image_data, np.ndarray):
raise TypeError("Expected imgs numpy ndarray type.")
- if imgs.ndim < 4:
+ if image_data.ndim < 4:
raise ValueError("Expected 4D imgs array.")
- n_cells, n_rows, n_cols = imgs.shape[1:]
+ n_cells, n_rows, n_cols = image_data.shape[1:]
h_chk = np.zeros(
- (len(h_bins)-1, n_cells, n_rows, n_cols),
+ (len(histogram_bins)-1, n_cells, n_rows, n_cols),
dtype=np.int32)
e_chk = None
for cell in range(n_cells):
for row in range(n_rows):
for col in range(n_cols):
- this_pix = imgs[:, cell, row, col]
- _h, _e = np.histogram(this_pix, bins=h_bins)
+ this_pix = image_data[:, cell, row, col]
+ _h, _e = np.histogram(this_pix, bins=histogram_bins)
h_chk[..., cell, row, col] += _h
if e_chk is None:
e_chk = np.array(_e)
@@ -237,7 +244,7 @@ def fit_histogram(
- CHARGE_SHARING_2: sum of two CHARGE_SHARING
- GAUSS: gaussian function
n_sigma (int): to calculate threshold of the peak finder as
- thr = n_sigma * sigma0
+ thr = n_sigma * initial_sigma
ratio (float): ratio parameter of the peak finder
histo (ndarray): histogram bin counts with shape
(n_bins, memory_cells, row, col)
@@ -253,12 +260,12 @@ def fit_histogram(
CHARGE_SHARING_2=fit_double_charge_sharing,
GAUSS=fit_gauss,
)
- sigma0 = 15.
+ initial_sigma = 15.
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)
- g0_chk = _init_array.copy()
+ peak_positions = _init_array.copy()
sigma_chk = _init_array.copy()
chi2ndf_chk = _init_array.copy()
alpha_chk = _init_array.copy()
@@ -271,19 +278,19 @@ def fit_histogram(
yerr = histogram_errors(y)
if noise[0, cell, row, col] > 0.:
- sigma0 = noise[0, cell, row, col]
+ initial_sigma = noise[0, cell, row, col]
q, sigma, chi2ndf, alpha = fit_func(
- x_rebinned, y, yerr, sigma0, n_sigma, ratio)
+ x_rebinned, y, yerr, initial_sigma, n_sigma, ratio)
- g0_chk[cell, row, col] = q
+ peak_positions[cell, row, col] = q
sigma_chk[cell, row, col] = sigma
chi2ndf_chk[cell, row, col] = chi2ndf
alpha_chk[cell, row, col] = alpha
- return g0_chk, sigma_chk, chi2ndf_chk, alpha_chk
+ return peak_positions, sigma_chk, chi2ndf_chk, alpha_chk
-def erfbox(x, ped, q0, s0):
+def erfbox(x, ped, q, sigma):
"""
Compute the error function box used in charge sharing calculations.
@@ -291,34 +298,34 @@ def erfbox(x, ped, q0, s0):
erfbox is mostly flat between ped and q0 and quickly goes to zero
outside the range [ped, q0]. The function integral is normalized to 1.
- The slope of the box sides is given by erf(x)/s0; if s0 == 0 then the
+ The slope of the box sides is given by erf(x)/sigma; if sigma == 0 then the
function goes abruptly to 0 at the edges.
Args:
x (array-like): Input values
ped (float): Pedestal value
q0 (float): Charge value
- s0 (float): Width parameter
+ sigma (float): Width parameter
Returns:
array-like: Computed erfbox values
"""
- if ped == q0:
+ if ped == q:
return np.zeros_like(x)
- if ped > q0:
- q0, ped = ped, q0
-
- slope = lambda x: erf(x / abs(s0)) if s0 != 0 else 2 * np.heaviside(x, 1) # noqa
- return (slope(x - ped) - slope(x - q0)) * 0.5 / (q0 - ped)
-
+ if ped > q:
+ q, ped = ped, q
-def double_peak(x, ped, q_a, sigma_a, alpha, norm_a, q_b, sigma_b, norm_b):
- k_a = charge_sharing(x, ped, q_a, sigma_a, alpha, norm_a)
- k_b = charge_sharing(x, ped, q_b, sigma_b, alpha, norm_b)
- return k_a + k_b
+ slope = lambda x: erf(x / abs(sigma)) if sigma != 0 else 2 * np.heaviside(x, 1) # noqa
+ return (slope(x - ped) - slope(x - q)) * 0.5 / (q - ped)
def gauss(x, amp, mean, sigma):
+ """
+ Compute a Gaussian function.
+
+ This function calculates the values of a Gaussian distribution
+ for given input parameters.
+ """
z = (x - mean) / sigma
return amp * np.exp(-0.5 * z**2)
@@ -328,11 +335,7 @@ def charge_sharing(x, ped, q, sigma, alpha, norm):
Model the detector response to a single photon with charge sharing.
This function combines a Gaussian peak
- (representing full charge collection)
- with a charge sharing tail.
-
- More details on JOURNAL OF SYNCHROTRON RADIATION 23, 1462-1473 (2016)
- - f(x) = norm * (gaussian_peak + charge_sharing_tail) / ((1 + α)²)
+ (representing full charge collection) with a charge sharing tail.
Args:
x (array-like): Input charge values
@@ -344,22 +347,42 @@ def charge_sharing(x, ped, q, sigma, alpha, norm):
Returns:
array-like: Modeled detector response
- """
- # A small threshold to avoid division by zero
- small_value_threshold = 1e-10
- log_term = (x - ped) / (q - ped)
- # Avoid log of zero or negative numbers
- log_term = np.maximum(log_term, small_value_threshold)
+ Notes: More details on JOURNAL OF SYNCHROTRON RADIATION 23, 1462-1473 (2016)
+ - f(x) = norm * (gaussian_peak + charge_sharing_tail) / ((1 + α)²)
+ """
A = ((1. - alpha)**2) / (sigma * np.sqrt(2. * np.pi))
B = 4. * alpha * (1. - alpha)
C = 4. * alpha**2
- charge_sh_tail = (B - C * np.log1p(log_term)) * erfbox(x, ped, q, sigma)
+ charge_sh_tail = (
+ B - C * np.log1p((x - ped) / (q - ped))) * erfbox(x, ped, q, sigma)
return norm * (charge_sh_tail + gauss(x, A, q, sigma)) / ((1. + alpha)**2)
+def double_peak(x, ped, q_a, sigma_a, alpha, norm_a, q_b, sigma_b, norm_b):
+ """
+ Model the detector response for two overlapping peaks with charge sharing.
+
+ This function combines two charge sharing models to represent two
+ overlapping peaks in the detector response.
+ """
+ k_a = charge_sharing(x, ped, q_a, sigma_a, alpha, norm_a)
+ k_b = charge_sharing(x, ped, q_b, sigma_b, alpha, norm_b)
+ return k_a + k_b
+
+
def set_fit_range(x, x1, x2):
+ """
+ Set the range of x values to be used in fitting.
+
+ This function determines the indices of x that fall within the
+ specified range [x1, x2].
+
+ Notes: If x2 <= x1, the function ensures at least two points are included
+ in the range. The returned x values are cast to np.float64 for
+ higher precision in fitting.
+ """
i1 = 0
i2 = len(x)
i = 0
@@ -389,7 +412,7 @@ def set_fit_range(x, x1, x2):
def fit_charge_sharing(
- x, y, yerr, sigma0, n_sigma, ratio, alpha0=0.3
+ x, y, yerr, initial_sigma, n_sigma, ratio, alpha0=0.3
):
"""
Perform charge sharing fit using iminuit optimization.
@@ -401,12 +424,13 @@ def fit_charge_sharing(
x (array-like): Input charge values
y (array-like): Observed counts or intensities
yerr (array-like): Errors on y values
- sigma0 (float): Initial guess for peak width
+ initial_sigma (float): Initial guess for peak width
n_sigma (float): Number of sigmas for threshold calculation
(not used in this implementation)
ratio (float): Ratio parameter (not used in this implementation)
alpha0: Initial guess for alpha
- Returns:AC
+
+ Returns:
tuple: (q, sigma, chi2ndf, alpha)
q (float): Fitted peak position
sigma (float): Fitted peak width
@@ -422,15 +446,14 @@ def fit_charge_sharing(
norm = np.sum(y) * (x[1] - x[0])
# Use _peak_position function to find initial q0
- _peaks, _ = _peak_position(x, y, thr=n_sigma*sigma0, ratio=ratio)
+ _peaks, _ = _peak_position(x, y, thr=n_sigma*initial_sigma, ratio=ratio)
if len(_peaks) > 0:
q0 = np.min(x[_peaks])
else:
return -1, -1, -1, -1
- x_fit, i1, i2 = set_fit_range(x, 0.5 * q0, q0 + 2. * sigma0)
+ x_fit, i1, i2 = set_fit_range(x, 0.5 * q0, q0 + 2. * initial_sigma)
- x_fit = x_fit
y_fit = y[i1:i2]
yerr_fit = yerr[i1:i2]
@@ -441,12 +464,12 @@ def fit_charge_sharing(
m = Minuit(
cost_function,
- ped=0., q=q0, sigma=sigma0, alpha=alpha0, norm=norm,
- error_ped=0.1, error_q=0.1*q0, error_sigma=0.1*sigma0,
+ ped=0., q=q0, sigma=initial_sigma, alpha=alpha0, norm=norm,
+ error_ped=0.1, error_q=0.1*q0, error_sigma=0.1*initial_sigma,
error_alpha=0.1, error_norm=0.1*norm,
fix_ped=True,
- limit_q=(q0 - sigma0, q0 + sigma0),
- limit_sigma=(0.1 * sigma0, 2. * sigma0),
+ limit_q=(q0 - initial_sigma, q0 + initial_sigma),
+ limit_sigma=(0.1 * initial_sigma, 2. * initial_sigma),
limit_alpha=(0.01 * alpha0, 1.),
limit_norm=(np.sqrt(0.1 * norm), 3. * norm),
errordef=1, # for least squares cost function
@@ -454,19 +477,18 @@ def fit_charge_sharing(
)
m.migrad()
- m.hesse()
if m.get_fmin().is_valid:
q = m.values['q']
sigma = m.values['sigma']
alpha = m.values['alpha']
- chi2ndf = m.fval / (len(x) - len(m.values))
+ chi2ndf = m.fval / (len(x_fit) - len(m.values))
else:
q, sigma, chi2ndf, alpha = -1000, -1000, -1000, -1000
return q, sigma, chi2ndf, alpha
-def fit_double_charge_sharing(x, y, yerr, sigma0, n_sigma, ratio):
+def fit_double_charge_sharing(x, y, yerr, initial_sigma, n_sigma, ratio):
"""
Fits histogram with sum of two single photon functions with
charge sharing tail.
@@ -475,9 +497,9 @@ def fit_double_charge_sharing(x, y, yerr, sigma0, n_sigma, ratio):
x (array, float): x values
y (array, float): y values
yerr (array, float): errors of the y values
- sigma0 (float): rough estimate of peak variance
+ initial_sigma (float): rough estimate of peak variance
n_sigma (int): to calculate threshold of the peak finder as
- thr = n_sigma * sigma0
+ thr = n_sigma * initial_sigma
ratio (float): ratio parameter of the peak finder
Returns:
@@ -487,14 +509,14 @@ def fit_double_charge_sharing(x, y, yerr, sigma0, n_sigma, ratio):
alpha0 = 0.3
norm = np.sum(y) * (x[1] - x[0])
- _peaks, _ = _peak_position(x, y, thr=n_sigma*sigma0, ratio=ratio)
+ _peaks, _ = _peak_position(x, y, thr=n_sigma*initial_sigma, ratio=ratio)
if len(_peaks) > 0:
q0 = np.min(x[_peaks])
else:
return -1, -1, -1, -1
q_b = 1.1 * q0
- x_fit, i1, i2 = set_fit_range(x, 0.5 * q0, q_b + sigma0)
+ x_fit, i1, i2 = set_fit_range(x, 0.5 * q0, q_b + initial_sigma)
y_fit = y[i1:i2]
yerr_fit = yerr[i1:i2]
@@ -509,25 +531,25 @@ def fit_double_charge_sharing(x, y, yerr, sigma0, n_sigma, ratio):
limit_alpha = (0.01, 1.)
limit_norm_a = (np.sqrt(0.1 * norm), 2. * norm)
- limit_q_b = (q_b - sigma0, q_b + sigma0)
- limit_sigma = (0.1 * sigma0, 2. * sigma0)
+ limit_q_b = (q_b - initial_sigma, q_b + initial_sigma)
+ limit_sigma = (0.1 * initial_sigma, 2. * initial_sigma)
limit_norm_b = (np.sqrt(0.1 * 0.2 * norm), 2 * 0.2 * norm)
m = Minuit(
cost_function,
ped=0.,
q_a=q0,
- sigma_a=sigma0,
+ sigma_a=initial_sigma,
alpha=alpha0,
norm_a=norm,
q_b=q_b,
- sigma_b=sigma0,
+ sigma_b=initial_sigma,
norm_b=norm,
- error_ped=0.1, error_q_a=0.1*q0, error_sigma_a=0.1*sigma0,
+ error_ped=0.1, error_q_a=0.1*q0, error_sigma_a=0.1*initial_sigma,
error_alpha=0.1, error_norm_a=0.1*norm,
- error_q_b=0.1*q_b, error_sigma_b=0.1*sigma0,
+ error_q_b=0.1*q_b, error_sigma_b=0.1*initial_sigma,
error_norm_b=0.1*norm,
- limit_q_a=(q0 - sigma0, q0 + sigma0),
+ limit_q_a=(q0 - initial_sigma, q0 + initial_sigma),
limit_sigma_a=limit_sigma,
limit_alpha=limit_alpha,
limit_norm_a=limit_norm_a,
@@ -540,19 +562,18 @@ def fit_double_charge_sharing(x, y, yerr, sigma0, n_sigma, ratio):
)
m.migrad()
- m.hesse()
if m.get_fmin().is_valid:
q = m.values['q_a']
sigma = m.values['sigma_a']
alpha = m.values['alpha']
- chi2ndf = m.fval / (len(x) - len(m.values))
+ chi2ndf = m.fval / (len(x_fit) - len(m.values))
else:
q, sigma, chi2ndf, alpha = -1000, -1000, -1000, -1000
return q, sigma, chi2ndf, alpha
-def fit_gauss(x, y, yerr, sigma0, n_sigma, ratio):
+def fit_gauss(x, y, yerr, initial_sigma, n_sigma, ratio):
"""
Fits histogram with a gaussian function
@@ -560,9 +581,9 @@ def fit_gauss(x, y, yerr, sigma0, n_sigma, ratio):
x (array, float): x values
y (array, float): y values
yerr (array, float): errors of the y values
- sigma0 (float): rough estimate of peak variance
+ initial_sigma (float): rough estimate of peak variance
n_sigma (int): to calculate threshold of the peak finder as
- thr = n_sigma * sigma0
+ thr = n_sigma * initial_sigma
ratio (float): ratio parameter of the peak finder
Returns:
@@ -570,15 +591,15 @@ def fit_gauss(x, y, yerr, sigma0, n_sigma, ratio):
(last one alway == 0)
all of them are kept for compatibility with the wrap around function
"""
- norm = np.sum(y) * (x[1] - x[0])/np.sqrt(2.*np.pi*sigma0**2)
+ norm = np.sum(y) * (x[1] - x[0])/np.sqrt(2.*np.pi*initial_sigma**2)
- _peaks, _ = _peak_position(x, y, thr=n_sigma*sigma0, ratio=ratio)
+ _peaks, _ = _peak_position(x, y, thr=n_sigma*initial_sigma, ratio=ratio)
if len(_peaks) > 0:
q0 = np.min(x[_peaks])
else:
return -1, -1, -1, -1
- x_fit, i1, i2 = set_fit_range(x, q0 - sigma0, q0 + 2.*sigma0)
+ x_fit, i1, i2 = set_fit_range(x, q0 - initial_sigma, q0 + 2.*initial_sigma)
y_fit = y[i1:i2]
yerr_fit = yerr[i1:i2]
@@ -589,23 +610,22 @@ def fit_gauss(x, y, yerr, sigma0, n_sigma, ratio):
m = Minuit(
cost_function,
- amp=norm, mean=q0, sigma=sigma0,
- error_amp=0.1*norm, error_mean=0.1*sigma0, error_sigma=0.1*sigma0,
+ amp=norm, mean=q0, sigma=initial_sigma,
+ error_amp=0.1*norm, error_mean=0.1*initial_sigma, error_sigma=0.1*initial_sigma,
limit_amp=(0.01*norm, 2.*norm),
- limit_mean=(q0-sigma0, q0+sigma0),
- limit_sigma=(0.1*sigma0, 3.*sigma0),
+ limit_mean=(q0-initial_sigma, q0+initial_sigma),
+ limit_sigma=(0.1*initial_sigma, 3.*initial_sigma),
errordef=1, # for least squares cost function
print_level=0,
)
m.migrad()
- m.hesse()
if m.get_fmin().is_valid:
q = m.values['mean']
sigma = m.values['sigma']
alpha = 0.
- chi2ndf = m.fval / (len(x) - len(m.values))
+ chi2ndf = m.fval / (len(x_fit) - len(m.values))
else:
q, sigma, chi2ndf, alpha = -1000, -1000, -1000, -1000
return q, sigma, chi2ndf, alpha
diff --git a/src/cal_tools/step_timing.py b/src/cal_tools/step_timing.py
index a68587dcbac336e9a85f7c1fafa5d44d9c823c11..1ce001c8c6cccc5965b320e11238b9362eebf494 100644
--- a/src/cal_tools/step_timing.py
+++ b/src/cal_tools/step_timing.py
@@ -17,12 +17,13 @@ class StepTimer:
self.t0 = t
self.t_latest = t
- def done_step(self, step_name):
+ def done_step(self, step_name, print_step=True):
"""Record a step timing, since the last call to done_step() or start()
"""
t = perf_counter()
self.steps[step_name].append(t - self.t_latest)
- print(f"{step_name}: {t - self.t_latest:.01f} s")
+ if print_step:
+ print(f"{step_name}: {t - self.t_latest:.01f} s")
self.t_latest = t
def timespan(self):