detectors/JFCalibration2
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

jungfrau calibration

Browse Source
This commit is contained in:
mazzol_a
2026-07-22 16:46:36 +02:00
parent edc8cfa5fd
commit 62752822b9
19 changed files with 1534 additions and 311 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CMakeLists.txt
-56
View File
@@ -1,56 +0,0 @@
cmake_minimum_required(VERSION 3.15)
project(
jungfraucalibration
DESCRIPTION "helper functions for Jungfrau calibration software"
HOMEPAGE_URL "https://gitea.psi.ch/detectors/JFCalibration2"
LANGUAGES CXX)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
include(FetchContent)
option(JF_FETCH_AARE "Fetch aare library from github" ON)
if(JF_FETCH_AARE)
FetchContent_Declare(aare GIT_REPOSITORY https://github.com/slsdetectorgroup/aare
GIT_TAG dev/jungfraucalibration) #use newest version - change to main once merged
FetchContent_MakeAvailable(aare)
install(
TARGETS aare_core
EXPORT ${TARGETS_EXPORT_NAME}
)
install(TARGETS aare_compiler_flags
EXPORT ${TARGETS_EXPORT_NAME}
) #mmh is this the way to go I think I will define the same compiler options twice now? can directly use aare_compiler_flags instead of compiler_flags
message(STATUS "target: aare")
else()
#set(AARE_INSTALL_PATH "/usr/local/aare" CACHE PATH "Installation directory for AARE")
list(APPEND CMAKE_PREFIX_PATH ${AARE_INSTALL_PATH})
message(STATUS "looking for aare in: ${CMAKE_PREFIX_PATH}")
find_package(aare REQUIRED)
if(TARGET aare_core)
message(STATUS "found aare_core target")
else()
message(FATAL_ERROR "aare_core target was not found!")
endif()
#message(STATUS "found aare: ${AARE_INCLUDE_DIRS}")
endif()
set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/BadChannels.cpp)
add_library(jungfraucalibration STATIC ${SourceFiles} ${PUBLICHEADERS})
target_include_directories(
jungfraucalibration PUBLIC "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>")
target_link_libraries(jungfraucalibration PUBLIC aare_core) # TODO: should it be public?
# TODO: add other compile options?
target_compile_features(jungfraucalibration PRIVATE cxx_std_17)
#add_subdirectory(examples)
include/BadChannels.hpp
-14
View File
@@ -1,14 +0,0 @@
#include "aare/JungfrauDataFile.hpp"
#include "aare/NDArray.hpp"
using namespace aare;
namespace jungfraucalibration
{
NDArray<bool, 2> CreateBadChannelPixelMask(JungfrauDataFile &pedestal_file, const size_t num_pedestals_g0, const size_t num_pedestals_g1, const size_t num_pedestals_g2);
NDArray<bool, 2> CreateBadChannelPixelMask(JungfrauDataFile &pedestals_g0_file, const JungfrauDataFile &pedestals_g1_file, const JungfrauDataFile &pedestal_g2_file);
} // namespace jungfraucalibration
include/JungfrauCalibration.hpp
View File
jfcal/JungfrauCalibration.py
+502
View File
@@ -0,0 +1,502 @@
import math
from pathlib import Path
from pyexpat import model
import time
from aare import JungfrauDataFile, File, Jungfrau, random_pixel
from jfcal.utils import get_first_file
import numpy as np
from dataclasses import dataclass
import logging
from functools import cached_property
from aare import PedestalTrackingPixelHistogram
#import pickle
from aare import GaussianChargeSharingKb
from jfcal.JungfrauFitParameters import FitParams, GAINTYPE
from jfcal.JungfrauCalibrationParameters import JungfrauCalibrationParameters
from jfcal.JungfrauCalibrationResult import JungfrauCalibrationResult
@dataclass
class HistogramParameters:
adu_min : float = 600.5 # minimum ADU value to consider for the histogram
adu_max : float = 1400.5 # maximum ADU value to consider for the histogram
bin_width : int = 4 # width of each bin in the histogram
sigma_multiplier : float = 1.0 # used to determine the threshold for updating the pedestal - residuals less than sigma_multiplier*std are considered - for sigma_multiplier <= 0 no pedestal update
batch_size : int = 100 # number of frames to process in each batch - frames in a batch are loaded sequentially - once a batch was loaded processing of the histogram calculation starts
n_threads : int = 16 # number of threads to use for calculating the histogram
logger = logging.getLogger(__name__)
class JungfrauCalibration:
def __init__(self, calibration_params: JungfrauCalibrationParameters, beam_energy : float = 8.0):
"""
Initialize calibration.
Parameters
----------
calibration_params : JungfrauCalibrationParameters
The calibration parameters.
beam_energy : float, optional
The energy of the beam in keV. Default is 8.0 keV.
"""
self.calibration_params = calibration_params
self._beam_energy : float = beam_energy # energy of the beam in keV
self.bad_channel_mask : np.ndarray = None
self._ph_computer : PedestalTrackingPixelHistogram
# TODO: easier to save but stored twice - e.g. also in ph_computer - pickle might really be the better option to save ph computer - but patching together mean etc. also has computational cost
self.histogram_values : np.ndarray
self.histogram_bin_edges : np.ndarray
self.pedestal_mean : np.ndarray
self.pedestal_std : np.ndarray
self._n_rows : int = Jungfrau.rows
self._n_cols : int = Jungfrau.cols
self.initial_fit_params : FitParams = None
# If one wants to make those parameters configurable can still set fitmodel from outside class
self.fitmodel = GaussianChargeSharingKb(max_calls = 4000, compute_errors = False, tolerance = 1e-2) # K_\alpha, K_\beta, gaussian charge sharing model
self.fit_result : dict = None # dictionary with "par" storing the fitted parameters and "chi2" storing the chi2 values for each pixel
self.calibration_result : JungfrauCalibrationResult = JungfrauCalibrationResult(self.beam_energy)
@property
def calibration_params(self) -> JungfrauCalibrationParameters:
return self._calibration_params
@calibration_params.setter
def calibration_params(self, calibration_params : JungfrauCalibrationParameters):
if not isinstance(calibration_params, JungfrauCalibrationParameters):
raise ValueError("calibration_params must be an instance of JungfrauCalibrationParameters.")
self._calibration_params = calibration_params
@property
def beam_energy(self) -> float:
return self._beam_energy
@property
def ph_computer(self) -> PedestalTrackingPixelHistogram:
return self._ph_computer
@ph_computer.setter
def ph_computer(self, ph_computer: PedestalTrackingPixelHistogram):
if not isinstance(ph_computer, PedestalTrackingPixelHistogram):
raise ValueError("ph_computer must be an instance of PedestalTrackingPixelHistogram.")
self.__dict__.pop("bin_centers", None)
self._ph_computer = ph_computer
def calculate_bad_pixels_mask(self, save_result : bool = False):
"""
Calculate the bad pixels mask for the Jungfrau detector and store it under self.bad_channel_mask. A pixel is considered bad if it does not have the expected gain in one of the pedestal frames.
Params:
save_result : bool, optional
Whether to save the bad pixels mask to a file. Default is False.
"""
if(self.calibration_params.num_pedestals_g0 is not None and self.calibration_params.num_pedestals_g1 is not None and self.calibration_params.num_pedestals_g2 is not None):
pedestal_file = JungfrauDataFile(get_first_file(self.calibration_params.pedestal_file_dir, self.calibration_params.pedestal_file_prefix))
# potential ROI
self._n_rows, self._n_cols = pedestal_file.rows, pedestal_file.cols
t0 = time.perf_counter()
_, g0_pedestal_frames = pedestal_file.read_n(self.calibration_params.num_pedestals_g0) # TODO: option to only read gain? - mmh reading things twice from filesystem also bad - needed for pedestal calculation - but only for G0?
t1 = time.perf_counter()-t0
logger.info(f'G0 pedestal file read took {t1:.2f}s or {self.calibration_params.num_pedestals_g0/t1:.2f} frames/s', flush=True)
_, g1_pedestal_frames = pedestal_file.read_n(self.calibration_params.num_pedestals_g1)
t1 = time.perf_counter()-t0
logger.info(f'G1 pedestal file read took {t1:.2f}s or {self.calibration_params.num_pedestals_g1/t1:.2f} frames/s', flush=True)
_, g2_pedestal_frames = pedestal_file.read_n(self.calibration_params.num_pedestals_g2)
t1 = time.perf_counter()-t0
logger.info(f'G2 pedestal file read took {t1:.2f}s or {self.calibration_params.num_pedestals_g2/t1:.2f} frames/s', flush=True)
else:
pedestal_g0_file = JungfrauDataFile(get_first_file(self.calibration_params.pedestal_file_dir, self.calibration_params.pedestal_g0_file_prefix))
pedestal_g1_file = JungfrauDataFile(get_first_file(self.calibration_params.pedestal_file_dir, self.calibration_params.pedestal_g1_file_prefix))
pedestal_g2_file = JungfrauDataFile(get_first_file(self.calibration_params.pedestal_file_dir, self.calibration_params.pedestal_g2_file_prefix))
if(pedestal_g0_file.rows != pedestal_g1_file.rows or pedestal_g0_file.rows != pedestal_g2_file.rows or pedestal_g0_file.cols != pedestal_g1_file.cols or pedestal_g0_file.cols != pedestal_g1_file.cols or pedestal_g1_file.cols != pedestal_g2_file.cols):
raise ValueError("Pedestal files have different dimensions.")
self._n_rows, self._n_cols = pedestal_g0_file.rows, pedestal_g0_file.cols
t0 = time.perf_counter()
_, g0_pedestal_frames = pedestal_g0_file.read_n(pedestal_g0_file.total_frames())
t1 = time.perf_counter()-t0
logger.info(f'G0 pedestal file read took {t1:.2f}s or {pedestal_g0_file.total_frames()/t1:.2f} frames/s')
t0 = time.perf_counter()
_, g1_pedestal_frames = pedestal_g1_file.read_n(pedestal_g1_file.total_frames())
t1 = time.perf_counter()-t0
logger.info(f'G1 pedestal file read took {t1:.2f}s or {pedestal_g1_file.total_frames()/t1:.2f} frames/s')
t0 = time.perf_counter()
_, g2_pedestal_frames = pedestal_g2_file.read_n(pedestal_g2_file.total_frames())
t1 = time.perf_counter()-t0
logger.info(f'G2 pedestal file read took {t1:.2f}s or {pedestal_g2_file.total_frames()/t1:.2f} frames/s')
# get gain from each pixel - update mask
bad_channel_mask = np.zeros((self._n_rows, self._n_cols), dtype=bool) # bad channels pixel mask
t0 = time.perf_counter()
mask0 = np.any(g0_pedestal_frames >> 14 != 0, axis=0)
mask1 = np.any(g1_pedestal_frames >> 14 != 1, axis=0)
mask2 = np.any(g2_pedestal_frames >> 14 != 3, axis=0)
self.bad_channel_mask = mask0 | mask1 | mask2
t1 = time.perf_counter()-t0
logger.info(f'Bad channel mask calculation took {t1:.2f}s')
if save_result:
self.save_bad_pixel_mask(self.calibration_params.output_dir / self.calibration_params.bad_pixel_mask_output_file)
def save_bad_pixel_mask(self, output_file_path : Path):
"""
Save the bad pixels mask to a file.
Params:
output_file_path:
Path : The path to the file to save the bad pixels mask to.
"""
if self.bad_channel_mask is None:
raise ValueError("Bad channel mask has not been calculated yet. Please run calculate_bad_pixels_mask() first.")
np.save(output_file_path, self.bad_channel_mask)
logger.info(f'Bad channel mask saved to {output_file_path}')
def load_bad_pixel_mask(self, bad_pixel_mask_file_path : Path):
"""
Load the bad pixels mask from a file.
Params:
bad_pixel_mask_file_path:
Path : The path to the file containing the bad pixels mask.
"""
if not bad_pixel_mask_file_path.exists():
raise ValueError(f"Bad pixel mask file {bad_pixel_mask_file_path} does not exist.")
self.bad_channel_mask = np.load(bad_pixel_mask_file_path)
def calculate_pedestal(self, adu_min : float = 600.5, adu_max : float = 1400.5, bin_width : int = 4, sigma_multiplier : float = 1.0, batch_size : int = 100, n_threads : int = 16):
"""
Calculate the pedestal for each pixel and sets up the Histogram class for further histogram calculations.
Params:
adu_min : float, optional
The minimum ADU value to consider for the histogram. Default is 600.5.
adu_max : float, optional
The maximum ADU value to consider for the histogram. Default is 1400.5.
bin_width : int, optional
The width of each bin in the histogram. Default is 4.
sigma_multiplier : float, optional
Used to determine the threshold for updating the pedestal. Residuals less than sigma_multiplier * std are considered. For sigma_multiplier <= 0, there is no pedestal update. Default is 1.0.
batch_size : int, optional
The number of frames to process in each batch. Frames in a batch are loaded sequentially. Once a batch was loaded, processing of the histogram calculation starts. Default is 100.
n_threads : int, optional
The number of threads to use for calculating the histogram. Default is 16.
"""
if(self.calibration_params.num_pedestals_g0 is not None):
pedestal_file_name = get_first_file(self.calibration_params.pedestal_file_dir, self.calibration_params.pedestal_file_prefix)
frames_to_read = self.calibration_params.num_pedestals_g0 # read twice now - bad - rather store - second option to pass view in PedestalTrackingHistogram
else:
pedestal_file_name = get_first_file(self.calibration_params.pedestal_file_dir, self.calibration_params.pedestal_g0_file_prefix)
frames_to_read = -1 # read all frames in file
file = JungfrauDataFile(pedestal_file_name)
if(file.rows != self._n_rows or file.cols != self._n_cols):
raise ValueError(f"Pedestal file {pedestal_file_name} has different dimensions than expected. Expected: ({self._n_rows}, {self._n_cols}), got: ({file.rows}, {file.cols})")
num_bins = int((self._adu_max - self._adu_min) / self._bin_width) # number of bins in histogram
if (self.bad_channel_mask is None):
logger.warning("Bad channel mask is not set. All pixels will be considered as good pixels for pedestal calculation.")
# create histogram computer
self._ph_computer = PedestalTrackingPixelHistogram(rows = self._n_rows, cols = self._n_cols, n_bins = num_bins, xmin = adu_min, xmax = adu_max, n_threads = n_threads, max_pending = batch_size, n_sigma = sigma_multiplier, mask = self.bad_channel_mask)
t0 = time.perf_counter()
self._ph_computer.process_pedestal_file(pedestal_file_name, max_frames=frames_to_read, verbose=True)
t1 = time.perf_counter()-t0
logger.info(f'Pedestal calculation took {t1:.2f}s or {frames_to_read/t1:.2f} frames/s')
# TODO: set pedestal mean and std already here? - consistent with usage but recalculated in calculate_histogram
def calculate_histogram(self, max_frames = None, save_result : bool = False):
"""
Calculate the histogram of each pixel. Update pedestal if residual data-pedestal within threshold std*sigma_multiplier.
Params:
max_frames : int, optional
The maximum number of frames to process. If None, all frames in the raw file are processed. Default is None.
save_result : bool, optional
Whether to save the histogram to a file. Default is False.
"""
# TODO: do we want one large histogram or one for pedestal peak, beta_peak, alpha_peak?
fname = get_first_file(self.calibration_params.raw_file_dir, self.calibration_params.raw_file_prefix_G0)
max_frames = self.max_frames if self.max_frames is not None else File(fname).total_frames
self._ph_computer.fill_from_file(fname, max_frames = max_frames, verbose = True)
self.histogram_values = self._ph_computer.values()
self.histogram_bin_edges = self._ph_computer.bin_edges()
self.pedestal_mean = self._ph_computer.pedestal_mean()
self.pedestal_std = self._ph_computer.pedestal_std()
if save_result:
self.save_histogram(self.calibration_params.output_dir / self.calibration_params.histogram_output_file)
def save_histogram(self, histogram_file_path : Path):
"""
Save the histogram to a file.
Params:
histogram_file_path:
Path : The path to the file to save the histogram to.
"""
if self.histogram_values is None or self.histogram_bin_edges is None or self.pedestal_mean is None or self.pedestal_std is None:
raise ValueError("Histogram has not been calculated yet. Please run calculate_histogram() first.")
# TODO: guess pickle is the best for intermediate saving - but maybe store as numpy as well - need to add support for pybind- a lot to store
#with open(output_dir / output_file_name, 'wb') as f:
# pickle.dump(self.ph_computer, f)
np.savez_compressed(histogram_file_path, values = self.histogram_values, bin_edges = self.histogram_bin_edges, pedestal_mean = self.pedestal_mean, pedestal_std = self.pedestal_std)
logger.info(f'Histogram saved to {histogram_file_path}')
def load_histogram(self, histogram_file_path : Path):
"""
Load the histogram from a file.
Params:
histogram_file_path:
Path : The path to the file containing the histogram.
"""
if not histogram_file_path.exists():
raise ValueError(f"Histogram file {histogram_file_path} does not exist.")
self.histogram_values, self.histogram_bin_edges, self.pedestal_mean, self.pedestal_std = np.load(histogram_file_path).values()
@cached_property
def bin_centers(self) -> np.ndarray:
"""
Get the bin centers of the histogram.
Returns:
np.ndarray: The bin centers of the histogram.
"""
if self.histogram_bin_edges is None:
raise ValueError("Histogram has not been calculated yet. Please run calculate_histogram() first.")
return self.histogram_bin_edges[:-1] + 0.5 * self.ph_computer.bin_width
def estimate_initital_fit_params(self, gain_type : GAINTYPE = GAINTYPE.G0, elastic_scattering : bool = False):
"""
Estimate the initial fit parameters for the fitting.
Params:
gain_type : GAINTYPE, optional
The gain type to use for the initial fit parameters. Default is GAINTYPE.G0
elastic_scattering : bool, optional
Whether to include elastic scattering when estimating fit parameters. Default is False
"""
self.initial_fit_params = FitParams(gain_type = gain_type)
# estimate K_alpha_mean, K_alpha_peak
adc_counts = self.histogram_values # histogram values for each pixel
adc_bin_centers = self.bin_centers
mean_K_alpha_mean : float = 0.0
mean_K_alpha_amplitude : float = 0.0
mean_elastic_scattering_intercept : float = 0.0
mean_elastic_scattering_slope : float = 0.0
num_estimates : int = 3 # number of random pixels to estimate initial fit parameters from
# TODO: is this neccessary is one enough?
for i in range(num_estimates):
is_bad_pixel : bool = True
while(is_bad_pixel):
pixel = random_pixel(0, self._n_rows, 0, self._n_cols)
is_bad_pixel = self.bad_channel_mask[pixel[0], pixel[1]]
# TODO: maybe have sepearte function that only estimates amplitude and mean of K_alpha_peak
elastic_scattering_intercept, elastic_scattering_slope, K_alpha_mean, _, K_alpha_amplitude, _, _, _ = self.fitmodel.estimate_par(adc_bin_centers, adc_counts[pixel[0], pixel[1], :], elastic_scattering)
mean_K_alpha_mean += K_alpha_mean
mean_K_alpha_amplitude += K_alpha_amplitude
mean_elastic_scattering_intercept += elastic_scattering_intercept
mean_elastic_scattering_slope += elastic_scattering_slope
self.initial_fit_params.K_alpha_mean.value = mean_K_alpha_mean / num_estimates
self.initial_fit_params.amplitude_K_alpha.value = mean_K_alpha_amplitude / num_estimates
self.initial_fit_params.elastic_scattering_intercept.value = mean_elastic_scattering_intercept / num_estimates
self.initial_fit_params.elastic_scattering_slope.value = mean_elastic_scattering_slope / num_estimates
def fit_function(self, fit_elastic_scattering : bool = False, initial_fit_params : FitParams = None, save_result : bool = False):
"""
Fit two Gaussian peaks for K_\alpha and K_beta including charge sharing to the histogram of each pixel.
Params:
fit_elastic_scattering : bool, optional
Whether to fit the elastic scattering parameters. Default is False.
initial_fit_params : FitParams, optional
The initial fit parameters to use for the fitting. If None, the initial fit parameters estimated by estimate_initital_fit_params() are used. Default is None.
save_result : bool, optional
Whether to save the fitted parameters to a file. Default is False.
"""
initial_fit_params = initial_fit_params if initial_fit_params is not None else self.initial_fit_params
if initial_fit_params is None:
logger.warning("Initial fit parameters are not set. Estimating initial fit parameters.")
self.estimate_initital_fit_params(elastic_scattering = fit_elastic_scattering)
self.fitmodel.SetParameter("p0", initial_fit_params.elastic_scattering_intercept.value)
self.fitmodel.SetParameter("p1", initial_fit_params.elastic_scattering_slope.value)
self.fitmodel.SetParameter("mu", initial_fit_params.K_alpha_mean.value)
self.fitmodel.SetParameter("sigma", initial_fit_params.sigma.value)
self.fitmodel.SetParameter("N", initial_fit_params.amplitude_K_alpha.value)
self.fitmodel.SetParameter("C", initial_fit_params.ratio_amplitude_charge_sharing.value)
self.fitmodel.SetParameter("kb_mean", initial_fit_params.ratio_mean_K_beta.value)
self.fitmodel.SetParameter("kb_frac", initial_fit_params.ratio_amplitude_K_beta.value)
for param_idx, param in enumerate(initial_fit_params):
if param.lower_bound is not None and param.upper_bound is not None:
self.fitmodel.SetParLimits(param_idx, param.lower_bound, param.upper_bound)
elif param.lower_bound is not None:
self.fitmodel.SetParLimits(param_idx, param.lower_bound, math.inf)
elif param.upper_bound is not None:
self.fitmodel.SetParLimits(param_idx, -math.inf, param.upper_bound)
else:
pass
if(not fit_elastic_scattering):
self.fitmodel.FixParameter("p0", self.initial_fit_params.elastic_scattering_intercept.value)
self.fitmodel.FixParameter("p1", self.initial_fit_params.elastic_scattering_slope.value)
adc_counts = self.histogram_values # histogram values for each pixel
adc_bins = self.bin_centers # histogram bin centers for each pixel
t0 = time.perf_counter()
res = self.fitmodel.fit(adc_bins, adc_counts, np.sqrt(adc_counts), n_threads = self.num_threads)
t = time.perf_counter()-t0
logger.info(f'Fit took {t:.2f}s or {adc_counts.shape[0]*adc_counts.shape[1]/t:.2f} pixels/s')
self.fit_result = res
if save_result:
self.save_fitted_parameters(self.calibration_params.output_dir / self.calibration_params.fit_params_output_file)
def save_fitted_parameters(self, fitted_parameters_file_path : Path):
"""
Save the fitted parameters to a file.
Params:
fitted_parameters_file_path:
Path : The path to the file to save the fitted parameters to.
"""
if self.fit_result is None:
raise ValueError("Fit result has not been calculated yet. Please run fit_function() first.")
np.savez_compressed(fitted_parameters_file_path, par = self.fit_result["par"], chi2 = self.fit_result["chi2"])
logger.info(f'Fitted parameters saved to {fitted_parameters_file_path}')
def load_fitted_parameters(self, fitted_parameters_file_path : Path):
"""
Load the fitted parameters from a file.
Params:
fitted_parameters_file_path:
Path : The path to the file containing the fitted parameters.
"""
if not fitted_parameters_file_path.exists():
raise ValueError(f"Fitted parameters file {fitted_parameters_file_path} does not exist.")
self.fit_result = np.load(fitted_parameters_file_path)
def calculate_G0(self):
"""
Calculate the G0 gain of the Jungfrau detector.
"""
self.calibration_result.G0 = self.fit_results["par"][:,:, 2]/self.beam_energy
def calibrate_G0(self, histogram_params : HistogramParameters = HistogramParameters(), max_frames : int = None, elastic_scattering : bool = False, high_gain0 : bool = False, save_intermediate_results : bool = False):
"""
Calibrate the G0 gain of the Jungfrau detector.
Params:
histogram_params : HistogramParameters, optional
The parameters for the histogram calculation. Default is HistogramParameters().
max_frames : int, optional
The maximum number of frames to take into account for histogram calculation. If None, all frames in the raw file are used.
elastic_scattering : bool, optional
Whether to include elastic scattering when estimating fit parameters. Default is False.
high_gain0 : bool, optional
Whether to use high gain 0 for the calibration. Default is False.
save_intermediate_results : bool, optional
Whether to save intermediate results (bad pixel mask, histogram, fitted parameters) to files.
"""
gain_type = GAINTYPE.HG0 if high_gain0 else GAINTYPE.G0
self.calculate_bad_pixels_mask(save_intermediate_results)
self.calculate_pedestal(adu_min = histogram_params.adu_min, adu_max = histogram_params.adu_max, bin_width = histogram_params.bin_width, sigma_multiplier = histogram_params.sigma_multiplier, batch_size = histogram_params.batch_size, n_threads = histogram_params.n_threads)
self.calculate_histogram(save_intermediate_results)
self.estimate_initital_fit_params(gain_type = gain_type, elastic_scattering = elastic_scattering)
self.fit_function(elastic_scattering = elastic_scattering, save_result = save_intermediate_results)
self.calculate_G0()
def calibrate_HG0(self, histogram_params : HistogramParameters = HistogramParameters(), max_frames : int = None, elastic_scattering : bool = False, save_intermediate_results : bool = False):
"""
Calibrate the HG0 gain of the Jungfrau detector.
Params:
histogram_params : HistogramParameters, optional
The parameters for the histogram calculation. Default is HistogramParameters().
max_frames : int, optional
The maximum number of frames to take into account for histogram calculation. If None, all frames in the raw file are used.
elastic_scattering : bool, optional
Whether to include elastic scattering when estimating fit parameters. Default is False.
save_intermediate_results : bool, optional
Whether to save intermediate results (bad pixel mask, histogram, fitted parameters) to files.
"""
self.calibrate_G0(histogram_params, max_frames, elastic_scattering, high_gain0 = True, save_intermediate_results = save_intermediate_results)
jfcal/JungfrauCalibrationParameters.py
+138
View File
@@ -0,0 +1,138 @@
from pathlib import Path
from dataclasses import dataclass
import json
import logging
logger = logging.getLogger(__name__)
# depercated decorator
def deprecated(message : str):
"""
Decorator to mark functions as deprecated. It will result in a warning being emitted when the function is used.
"""
def decorator(func):
def new_func(*args, **kwargs):
logger.warning(f"Call to deprecated function {func.__name__}. {message}")
return func(*args, **kwargs)
return new_func
return decorator
@dataclass
class JungfrauInputParameters:
"""
A class to hold input parameters for the Jungfrau detector calibration.
"""
_pedestal_file_dir : Path = None
#deprecated - use pedestal_g0_file_prefix instead
_pedestal_file_prefix : str = None
_num_pedestals_g0 : int = 1000
_num_pedestals_g1 : int = 1000
_num_pedestals_g2 : int = 1000
pedestal_g0_file_prefix : str = None
pedestal_g1_file_prefix : str = None
pedestal_g2_file_prefix : str = None
_raw_file_dir : Path = None
raw_file_prefix : str = None
# outputs
_output_dir : Path = Path.cwd() # default to current working directory
bad_pixel_mask_output_file : str = "bad_pixels_mask.npy"
histogram_output_file : str = "histogram.npz"
fit_params_output_file : str = "fit_parameters.npz"
@classmethod
def from_config(cls, config_file : Path):
"""
Create an instance of JungfrauInputParameters from a json configuration file.
"""
cls = cls()
if( not (config_file.exists() and config_file.is_file())):
raise ValueError(f"Configuration file {config_file} does not exist or is not a file.")
with open(config_file, 'r') as f:
config = json.load(f)
for key, value in config.items():
if hasattr(cls, key):
if(value is not None):
setattr(cls, key, value)
else:
logger.warning(f"Unknown configuration parameter {key} in {config_file}.")
# TODO: maybe add HistogramParameters to JungfrauInputParameters and load them here as well. and all other parameters e.g. high_gain0, elastic_scattering etc.
@property
def pedestal_file_dir(self) -> Path:
return self._pedestal_file_dir
@pedestal_file_dir.setter
def pedestal_file_dir(self, filepath : Path):
if not filepath.exists():
raise ValueError(f"Pedestal file directory {filepath} does not exist.")
self._pedestal_file_dir = filepath
@deprecated("Setting pedestal_file_prefix is deprecated use a distinct file for each gain instead.")
@property
def pedestal_file_prefix(self) -> str:
return self._pedestal_file_prefix
@pedestal_file_prefix.setter
def pedestal_file_prefix(self, file_prefix : str):
self._pedestal_file_prefix = file_prefix
@property
def num_pedestals_g0(self) -> int:
return self._num_pedestals_g0
@deprecated("Setting num_pedestals_g0 is deprecated use different files for each gain instead.")
@num_pedestals_g0.setter
def num_pedestals_g0(self, num_pedestals : int):
self._num_pedestals_g0 = num_pedestals
@property
def num_pedestals_g1(self) -> int:
return self._num_pedestals_g1
@deprecated("Setting num_pedestals_g1 is deprecated use different files for each gain instead.")
@num_pedestals_g1.setter
def num_pedestals_g1(self, num_pedestals : int):
self._num_pedestals_g1 = num_pedestals
@property
def num_pedestals_g2(self) -> int:
return self._num_pedestals_g2
@deprecated("Setting num_pedestals_g2 is deprecated use different files for each gain instead.")
@num_pedestals_g2.setter
def num_pedestals_g2(self, num_pedestals : int):
self._num_pedestals_g2 = num_pedestals
@property
def raw_file_dir(self) -> Path:
return self._raw_file_dir
@raw_file_dir.setter
def raw_file_dir(self, filepath : Path):
if not filepath.exists():
raise ValueError(f"Raw file directory {filepath} does not exist.")
self._raw_file_dir = filepath
@property
def output_dir(self) -> Path:
return self._output_dir
@output_dir.setter
def output_dir(self, filepath : Path):
if not filepath.exists():
raise ValueError(f"Output directory {filepath} does not exist.")
self._output_dir = filepath
jfcal/JungfrauCalibrationResult.py
+21
View File
@@ -0,0 +1,21 @@
import numpy as np
class JungfrauCalibrationResult:
def __init__(self, beam_energy : float):
"""
Class to store the results of the Jungfrau calibration.
Parameters
----------
beam_energey : float
The energy of the beam in keV.
"""
self.beam_energy : float = beam_energy # in keV -> should be frozen set by calibration
self.gain0 : np.ndarray
self.gain1 : np.ndarray
self.gain2 : np.ndarray
jfcal/JungfrauFitParameters.py
+69
View File
@@ -0,0 +1,69 @@
from dataclasses import dataclass
from enum import Enum
from typing import ClassVar
class GAINTYPE(Enum):
G0 = 0 # gain 0
HG0 = 1 # high gain 0
#G1 = 2 # gain 1
#G2 = 3 # gain 2
@dataclass
class Parameter:
value : float # the value of the parameter
lower_bound : float = None # a parameter limit of None defaults to -infinity
upper_bound : float = None # a parameter limit of None defaults to infinity
@dataclass
class FitParams:
_CHARGE_SHARING_G0 : ClassVar[float] = 16.0 # sigma gaussian peaks used as initial guess
_CHARGE_SHARING_HG0 : ClassVar[float] = 29.0 # sigma gaussian peaks used as initial guess
_RATIO_CHARGE_SHARING_K_ALPHA_G0 : ClassVar[float] = 0.17 # ratio of charge sharing peak to K_alpha peak used as initial guess
_RATIO_CHARGE_SHARING_K_ALPHA_HG0 : ClassVar[float] = 0.14
_RATIO_MEAN_K_BETA : ClassVar[float] = 1.12 # ratio of K_beta mean to K_alpha mean used as initial guess # 8.04 keV, 8.9 keV - 8.9/8.04 = 1.106
_RATIO_AMPLITUDE_K_BETA_G0 : ClassVar[float] = 0.12 # ratio of K_beta amplitude to K_alpha amplitude used as initial guess
_RATIO_AMPLITUDE_K_BETA_HG0 : ClassVar[float] = 0.14
# module parameters for elastic scattering - moduled as linear function - only important for trailing edge of K_\beta peak
elastic_scattering_intercept : Parameter = Parameter(0.0, None, None)
elastic_scattering_slope : Parameter = Parameter(0.0, None, None)
# module parameters for K_\alpha peak
K_alpha_mean : Parameter = Parameter(None, 850, 1350) # TODO: should there be a rule? for the bounds
amplitude_K_alpha : Parameter = Parameter(None, 0, 500)
# noise
sigma : Parameter | None = None
# module parameters for charge sharing
ratio_amplitude_charge_sharing : Parameter | None = None # amplitude_{charge_sharing}/amplitude_{K_\alpha}
# module parameters for K_\beta peak
ratio_mean_K_beta : Parameter = Parameter(_RATIO_MEAN_K_BETA, 1.02, 1.30) # \mu_{K_\beta}/\mu_{K_\alpha}
ratio_amplitude_K_beta : Parameter | None = None # \amplitude_{K_\beta}/\amplitude_{K_\alpha}
def __post_init__(self, gain_type : GAINTYPE = GAINTYPE.G0):
if gain_type == GAINTYPE.G0:
self.sigma = Parameter(FitParams._CHARGE_SHARING_G0, 5, 50)
self.ratio_charge_sharing = Parameter(FitParams._RATIO_CHARGE_SHARING_K_ALPHA_G0, None, None)
self.ratio_amplitude_K_beta = Parameter(FitParams._RATIO_AMPLITUDE_K_BETA_G0, 0.05, 0.4)
elif gain_type == GAINTYPE.HG0:
self.sigma = Parameter(FitParams._CHARGE_SHARING_HG0, 5, 50)
self.ratio_charge_sharing = Parameter(FitParams._RATIO_CHARGE_SHARING_K_ALPHA_HG0, None, None)
self.ratio_amplitude_K_beta = Parameter(FitParams._RATIO_AMPLITUDE_K_BETA_HG0, 0.05, 0.4)
else:
raise ValueError(f"Unsupported gain type: {gain_type}")
def __iter__(self):
yield self.elastic_scattering_intercept
yield self.elastic_scattering_slope
yield self.K_alpha_mean
yield self.sigma
yield self.amplitude_K_alpha
yield self.ratio_amplitude_charge_sharing
yield self.ratio_mean_K_beta
yield self.ratio_amplitude_K_beta
jfcal/PlotHelpers.py
+139
View File
@@ -0,0 +1,139 @@
from pprint import pp
import matplotlib.pyplot as plt
import numpy as np
from aare import add_colorbar
def plot_histogram(histogram_data : np.ndarray, bin_edges : np.ndarray, Count_range : tuple[int, int] = None, bin_range : tuple[int, int] = None, label : str = None, title : str = None, xlabel : str = None, axis : plt.Axes = None) -> plt.Axes:
"""
Plot the histogram of the pedestal values.
Parameters
----------
histogram_data : np.ndarray
The histogram data to plot.
bin_edges : np.ndarray
The edges of the bins for the histogram.
Count_range : tuple[int, int], optional
The range of counts to display on the y-axis. Default is None.
bin_range : tuple[int, int], optional
The range of bin values to display on the x-axis. If None, it will be set to the range of the bin edges. Default is None.
label : str, optional
The label for the histogram. Default is None.
title : str, optional
The title for the plot. Default is None.
xlabel : str, optional
The label for the x-axis. Default is None.
axis : plt.Axes, optional
The axis to plot on. If None, a new figure and axis will be created.
Returns
-------
plt.Axes
The axis with the histogram plot.
"""
if axis is None:
fig, ax = plt.subplots(figsize = (8,5))
else:
ax = axis
ax.stairs(histogram_data, bin_edges, label = label, zorder = 3)
if bin_range is None:
bin_range = (bin_edges[0]-0.01*bin_edges[0], bin_edges[-1]+0.01*bin_edges[-1])
if Count_range is not None:
ax.set_ylim(*Count_range)
ax.set_xlim(*bin_range)
ax.grid(zorder = 0)
ax.set_title(title)
ax.set_xlabel(xlabel)
ax.set_ylabel('Counts')
if label is not None:
ax.legend()
return ax
def plot_fitted_function(histogram_data : np.ndarray, bin_edges : np.ndarray, function : callable, Count_range : tuple[int, int] = None, bin_range : tuple[int, int] = None, label : str = None, xlabel : str = None, title : str = None, axis : plt.Axes = None) -> plt.Axes:
"""
Plot the histogram of the pedestal values along with the fitted function.
Parameters
----------
histogram_data : np.ndarray
The histogram data to plot.
bin_edges : np.ndarray
The edges of the bins for the histogram.
function : callable
The fitted function to plot.
Count_range : tuple[int, int], optional
The range of counts to display on the y-axis. Default is None.
bin_range : tuple[int, int], optional
The range of bin values to display on the x-axis. If None, it will be set to the range of the bin edges. Default is None.
label : str, optional
The label for the fitted function. Default is None.
xlabel : str, optional
The label for the x-axis. Default is None.
title : str, optional
The title for the plot. Default is None.
axis : plt.Axes, optional
The axis to plot on. If None, a new figure and axis will be created.
Returns
-------
plt.Axes
The axis with the histogram and fitted function plot.
"""
ax = plot_histogram(histogram_data, bin_edges, Count_range = Count_range, bin_range = bin_range, xlabel = xlabel, axis = axis, title = title)
bin_centers = bin_edges[:-1] + np.diff(bin_edges)/2
ax.plot(bin_centers, function(bin_centers), label=label, color='red', zorder=4)
if label is not None:
ax.legend(prop={"family": "monospace", "size": 10}, loc= 'upper left')
return ax
def plot_parameter(fit_parameter : np.ndarray, parameter_name : str, suppress_outliers : bool = True) -> None:
"""
Plot the parameter for all pixels.
Parameters
----------
fit_parameter : np.ndarray
The parameter values to plot.
parameter_name : str
The name of the parameter to plot - used for plot title.
suppress_outliers : bool, optional
Whether to suppress outliers in the plot (True, don't plot outliers). Default is True.
"""
fig, ax = plt.subplots(figsize = (15,5))
im = ax.imshow(fit_parameter)
#suppress outliers for colorbar
if suppress_outliers:
mean = np.mean(fit_parameter)
std = np.std(fit_parameter)
im.set_clim(mean-3*std,mean+3*std)
ax.set_xlabel('Pixel X')
ax.set_ylabel('Pixel Y')
ax.set_title(parameter_name)
add_colorbar(ax, im)
plt.show()
jfcal/Plotter.py
+213
View File
@@ -0,0 +1,213 @@
from jfcal import JungfrauCalibration
from jfcal.PlotHelpers import plot_histogram, plot_fitted_function, plot_parameter
from jfcal.utils import create_histogram_from_data
import matplotlib.pyplot as plt
from enum import Enum
class Parameters(Enum):
ELASTIC_SCATTERING_INTERCEPT = "elastic_scattering_intercept"
ELASTIC_SCATTERING_SLOPE = "elastic_scattering_slope"
K_ALPHA_MEAN = "k_alpha_mean"
SIGMA = "sigma"
K_ALPHA_AMPLITUDE = "k_alpha_amplitude"
RATIO_AMPLITUDE_CHARGE_SHARING = "ratio_amplitude_charge_sharing"
RATIO_MEAN_K_BETA = "ratio_mean_k_beta"
RATIO_AMPLITUDE_K_BETA = "ratio_amplitude_k_beta"
class Plotter:
def __init__(self, jungfrau_calibration : JungfrauCalibration):
self.calibration = jungfrau_calibration
def plot_ADU_histogram(self, pixel : tuple[int, int], Count_range : tuple[int, int] = None, ADU_range : tuple[int, int] = None, label = None, axis : plt.Axes = None, show_plot : bool = True) -> plt.Axes:
"""
Plot the histogram of ADU values for a specific pixel.
Parameters
----------
pixel : tuple[int, int]
The pixel coordinates to plot the histogram for.
Count_range : tuple[int, int], optional
The range of counts to display on the y-axis. Default is None.
ADU_range : tuple[int, int], optional
The range of ADU values to display on the x-axis. If None, it will be set to the range of the bin edges. Default is None.
label : str, optional
The label for the histogram. Default is None.
axis : plt.Axes, optional
The axis to plot on. If None, a new figure and axis will be created.
show_plot : bool, optional
Whether to display the plot. Default is True.
Returns
-------
plt.Axes
The axis with the histogram plot.
"""
ax = plot_histogram(self.calibration.histogram_values[pixel[0], pixel[1], :], self.calibration.histogram_bin_edges, Count_range = Count_range, bin_range = ADU_range, label = label, title = f"ADU Histogram for pixel {pixel}", xlabel = "ADU", axis = axis)
if show_plot:
plt.show()
return ax
def plot_fitted_function(self, pixel : tuple[int, int], Count_range : tuple[int, int] = None, ADU_range : tuple[int, int] = None, axis : plt.Axes = None, show_plot : bool = True) -> plt.Axes:
"""
Plot the histogram of the ADU values along with the fitted function for a specific pixel.
Parameters
----------
pixel : tuple[int, int]
The pixel coordinates to plot the histogram and fitted function for.
Count_range : tuple[int, int], optional
The range of counts to display on the y-axis. Default is None.
ADU_range : tuple[int, int], optional
The range of ADU values to display on the x-axis. If None, it will be set to the range of the bin edges. Default is None.
axis : plt.Axes, optional
The axis to plot on. If None, a new figure and axis will be created.
show_plot : bool, optional
Whether to display the plot. Default is True.
Returns
-------
plt.Axes
The axis with the histogram and fitted function plot.
"""
function_parameters = self.calibration.fit_results["par"][pixel[0],pixel[1],:]
variable_name_width = 8
variable_width = 8
decimal_places = 3
alpha = "\u03B1"
beta = "\u03B2"
sigma = "\u03C3"
mu = "\u03BC"
fit_label = (
f"{f'k{alpha}_{mu}:':<{variable_name_width}}{function_parameters[2]:>{variable_width}.{decimal_places}f}\n"
f"{f'{sigma}:':<{variable_name_width}}{function_parameters[3]:>{variable_width}.{decimal_places}f}\n"
f"{f'k{alpha}_N:':<{variable_name_width}}{function_parameters[4]:>{variable_width}.{decimal_places}f}\n"
f"{f'C:':<{variable_name_width}}{function_parameters[5]:>{variable_width}.{decimal_places}f}\n"
f"{f'k{beta}_m:':<{variable_name_width}}{function_parameters[6]:>{variable_width}.{decimal_places}f}\n"
f"{f'k{beta}_f:':<{variable_name_width}}{function_parameters[7]:>{variable_width}.{decimal_places}f}"
)
func = lambda x : self.calibration.fitmodel(x, function_parameters)
ax = plot_fitted_function(self.calibration.histogram_values[pixel[0], pixel[1], :], self.calibration.histogram_bin_edges, function=func, Count_range = Count_range, bin_range = ADU_range, label = fit_label, xlabel = "ADU", title = f"Fitted Function for pixel {pixel}", axis = axis)
if show_plot:
plt.show()
return ax
def plot_fitted_parameter(self, parameter_name : str, suppress_outliers: bool = True):
"""
Plot the parameter for all pixels.
Parameters
----------
parameter_name : str
The name of the parameter to plot. Must be one of the following:
"elastic_scattering_intercept", "elastic_scattering_slope", "k_alpha_mean", "sigma", "k_alpha_amplitude", "ratio_amplitude_charge_sharing", "ratio_mean_k_beta", "ratio_amplitude_k_beta".
suppress_outliers : bool, optional
Whether to suppress outliers in the plot. Default is True.
"""
match parameter_name:
case Parameters.ELASTIC_SCATTERING_INTERCEPT.value:
plot_parameter(self.calibration.fit_results["par"][:, :, 0], parameter_name="Elastic scattering intercept", suppress_outliers=suppress_outliers)
case Parameters.ELASTIC_SCATTERING_SLOPE.value:
plot_parameter(self.calibration.fit_results["par"][:, :, 1], parameter_name="Elastic scattering slope", suppress_outliers=suppress_outliers)
case Parameters.K_ALPHA_MEAN.value:
plot_parameter(self.calibration.fit_results["par"][:, :, 2], parameter_name="Cu K_alpha mean", suppress_outliers=suppress_outliers)
case Parameters.SIGMA.value:
plot_parameter(self.calibration.fit_results["par"][:, :, 3], parameter_name="Charge sharing sigma", suppress_outliers=suppress_outliers)
case Parameters.K_ALPHA_AMPLITUDE.value:
plot_parameter(self.calibration.fit_results["par"][:, :, 4], parameter_name="Cu K_alpha amplitude", suppress_outliers=suppress_outliers)
case Parameters.RATIO_AMPLITUDE_CHARGE_SHARING.value:
plot_parameter(self.calibration.fit_results["par"][:, :, 5], parameter_name="Charge sharing amplitude ratio", suppress_outliers=suppress_outliers)
case Parameters.RATIO_MEAN_K_BETA.value:
plot_parameter(self.calibration.fit_results["par"][:, :, 6], parameter_name="Cu K_beta mean ratio", suppress_outliers=suppress_outliers)
case Parameters.RATIO_AMPLITUDE_K_BETA.value:
plot_parameter(self.calibration.fit_results["par"][:, :, 7], parameter_name="Cu K_beta amplitude ratio", suppress_outliers=suppress_outliers)
case _:
raise ValueError(f"Unknown parameter name: {parameter_name}. Valid options are: {[param.value for param in Parameters]}")
def plot_parameter_histogram(self, parameter_name : str, bin_range : tuple[float, float] = None, bin_width : float = None, axis : plt.Axes = None, show_plot : bool = True) -> plt.Axes:
"""
Plot the histogram of the fitted parameter for all pixels.
Parameters
----------
parameter_name : str
The name of the parameter to plot. Must be one of the following:
"elastic_scattering_intercept", "elastic_scattering_slope", "k_alpha_mean", "sigma", "k_alpha_amplitude", "ratio_amplitude_charge_sharing", "ratio_mean_k_beta", "ratio_amplitude_k_beta".
bin_range : tuple[float, float], optional
The range of bin values to display on the x-axis. If None, it will be set to the range of the bin edges. Default is None.
bin_width : float, optional
The width of the bins for the histogram. If None, it will be set to the default bin width of bin_range / 200. Default is None.
axis : plt.Axes, optional
The axis to plot on. If None, a new figure and axis will be created.
show_plot : bool, optional
Whether to display the plot. Default is True.
Returns
-------
plt.Axes
The axis with the histogram plot.
"""
match parameter_name:
case Parameters.ELASTIC_SCATTERING_INTERCEPT.value:
data = self.calibration.fit_results["par"][:, :, 0].flatten()
histogram = create_histogram_from_data(data, bin_range = bin_range, bin_width = bin_width)
ax = plot_histogram(histogram, histogram.axes[0].edges[:], xlabel=r"Elastic scattering intercept", axis=axis)
case Parameters.ELASTIC_SCATTERING_SLOPE.value:
data = self.calibration.fit_results["par"][:, :, 1].flatten()
histogram = create_histogram_from_data(data, bin_range = bin_range, bin_width = bin_width)
ax = plot_histogram(histogram, histogram.axes[0].edges[:], xlabel=r"Elastic scattering slope", axis=axis)
case Parameters.K_ALPHA_MEAN.value:
data = self.calibration.fit_results["par"][:, :, 2].flatten()
histogram = create_histogram_from_data(data, bin_range = bin_range, bin_width = bin_width)
ax = plot_histogram(histogram, histogram.axes[0].edges[:], xlabel=r"$k_\alpha$ mean", axis=axis)
case Parameters.SIGMA.value:
data = self.calibration.fit_results["par"][:, :, 3].flatten()
histogram = create_histogram_from_data(data, bin_range = bin_range, bin_width = bin_width)
ax = plot_histogram(histogram, histogram.axes[0].edges[:], xlabel=r"$\sigma$", axis=axis)
case Parameters.K_ALPHA_AMPLITUDE.value:
data = self.calibration.fit_results["par"][:, :, 4].flatten()
histogram = create_histogram_from_data(data, bin_range = bin_range, bin_width = bin_width)
ax = plot_histogram(histogram, histogram.axes[0].edges[:], xlabel=r"$k_\alpha$ amplitude", axis=axis)
case Parameters.RATIO_AMPLITUDE_CHARGE_SHARING.value:
data = self.calibration.fit_results["par"][:, :, 5].flatten()
histogram = create_histogram_from_data(data, bin_range = bin_range, bin_width = bin_width)
ax = plot_histogram(histogram, histogram.axes[0].edges[:], xlabel=r"Ratio amplitude charge sharing", axis=axis)
case Parameters.RATIO_MEAN_K_BETA.value:
data = self.calibration.fit_results["par"][:, :, 6].flatten()
histogram = create_histogram_from_data(data, bin_range = bin_range, bin_width = bin_width)
ax = plot_histogram(histogram, histogram.axes[0].edges[:], xlabel=r"Ratio $k_\beta$ mean", axis=axis)
case Parameters.RATIO_AMPLITUDE_K_BETA.value:
data = self.calibration.fit_results["par"][:, :, 7].flatten()
histogram = create_histogram_from_data(data, bin_range = bin_range, bin_width = bin_width)
ax = plot_histogram(histogram, histogram.axes[0].edges[:], xlabel=r"Ratio $k_\beta$ amplitude", axis=axis)
case _:
raise ValueError(f"Unknown parameter name: {parameter_name}. Valid options are: {[param.value for param in Parameters]}")
if show_plot:
plt.show()
return ax
jfcal/__init__.py
+14 -1
View File
@@ -1 +1,14 @@
from .helpers import *
from .helpers import *
from .JungfrauCalibration import *
from .JungfrauFitParameters import *
from .JungfrauCalibrationParameters import *
from .JungfrauCalibrationResult import *
from .PlotHelpers import *
from .Plotter import *
jfcal/helpers.py
+1 -5
View File
@@ -50,8 +50,4 @@ def get_fname(path, gain, label = 'CuFluo', index = -1):
return file_sets[index] #
def get_fname(path : Path, file_prefix : str):
first_file = min(path.glob(str(Path)+f'{file_prefix}_*', default=None))
if first_file is None:
raise ValueError(f"No files found in {path} with prefix {file_prefix}")
return first_file
jfcal/utils.py
+66
View File
@@ -0,0 +1,66 @@
import glob
from pathlib import Path
import boost_histogram as bh
import numpy as np
def get_first_file(path : Path, file_prefix : str):
"""
Get the first file with lowest index in directory that matches the given prefix.
"""
first_file = min(path.glob(f'{file_prefix}*'), default=None)
if first_file is None:
raise ValueError(f"No files found in {path} with prefix {file_prefix}")
return first_file
def save_fit_parameters(fit_params : dict, output_file : Path):
"""
Save the fit parameters to a text file.
Parameters
----------
fit_params : dict
The fit parameters to save.
output_file : Path
The path to the output file.
"""
with open(output_file, 'w') as f:
for key, value in fit_params.items():
f.write(f"{key}: {value}\n")
def create_histogram_from_data(data : np.ndarray, bin_range : tuple[float, float] = None, bin_width : float = None) -> bh.Histogram:
"""
Create a histogram from the given data.
Parameters
----------
data : np.ndarray
The data to create the histogram from.
bin_range : tuple[float, float], optional
The range of the bins. Default is None, which means the range is determined from the data.
bin_width : float, optional
The width of each bin. Default is None, which means the number of bins is determined automatically.
Returns
-------
bh.Histogram
The created boost histogram.
"""
if bin_range is None:
min = np.min(data)
max = np.max(data)
bin_range = (min - 0.05*(max - min), max + 0.05*(max - min)) # add 5% margin to the range
if bin_width is None:
bins = 200 # 0.5 %
else:
bins = int((bin_range[1] - bin_range[0]) / bin_width)
hist = bh.Histogram(bh.axis.Regular(bins, bin_range[0], bin_range[1]))
hist.fill(data)
return hist
jungfrauinput.json
+16
View File
@@ -0,0 +1,16 @@
{
"pedestal_file_dir": "/mnt/sls_det_storage/jungfrau_calib/data/Module_749_Calib",
"pedestal_g0_file_prefix": null,
"pedestal_g1_file_prefix": null,
"pedestal_g2_file_prefix": null,
"pedestal_file_prefix": "pedeHG0_M749_2026-05-18_000000.dat",
"num_pedestals_g0": 1000,
"num_pedestals_g1": 1000,
"num_pedestals_g2": 1000,
"raw_file_dir": "/mnt/sls_det_storage/jungfrau_calib/data/Module_749_Calib",
"raw_file_prefix": "CuFluoHG0_M749_2026-05-18_",
"output_dir": "/mnt/sls_det_storage/jungfrau_calib/data/Module_749_Calib_output",
"bad_pixel_mask_output_file": null,
"histogram_output_file": null,
"fit_params_output_file": null
}
notebooks/JFCalibration2.ipynb
+355
View File
File diff suppressed because one or more lines are too long
python/jungfrau.config
View File
python/src/JungfrauCalibration.py
-96
View File
@@ -1,96 +0,0 @@
from pathlib import Path
import JungfrauCalibrationParameters
from aare import JungfrauDataFile
from helpers import get_first_file
import numpy as np
from aare.calibration import get_gain
class JungfrauCalibration:
def __init__(self, calibration_params: JungfrauCalibrationParameters):
self.calibration_params = calibration_params
self.bad_channel_mask : np.ndarray
self.histogram : np.ndarray
# TODO: maybe pass num pedestals, pedestal file instead of calibration params?
def calculate_bad_pixels_mask(self) -> np.ndarray:
"""
Calculate the bad pixels mask for the Jungfrau detector.
Returns:
np.ndarray: A boolean array where True indicates a bad pixel.
"""
if(self.calibration_params.num_pedestals_g0 is not None and self.calibration_params.num_pedestals_g1 is not None and self.calibration_params.num_pedestals_g2 is not None):
jungfrau_file = JungfrauDataFile(get_first_file(self.calibration_params.pedestal_file_dir, self.calibration_params.pedestal_file_prefix))
g0_pedestal_frames = jungfrau_file.read_n(self.calibration_params.num_pedestals_g0) # TODO: option to only read gain? - mmh reading things twice from filesystem also bad
g1_pedestal_frames = jungfrau_file.read_n(self.calibration_params.num_pedestals_g1)
g2_pedestal_frames = jungfrau_file.read_n(self.calibration_params.num_pedestals_g2)
# get gain from each pixel - update mask
max_frames = max(self.calibration_params.num_pedestals_g0, self.calibration_params.num_pedestals_g1, self.calibration_params.num_pedestals_g2)
bad_channel_mask = np.zeros((jungfrau_file.rows(), jungfrau_file.cols()), dtype=bool) # bad channels pixel mask
# TODO loop to etensive
return bad_channel_mask
def calculate_histogram(self):
"""
Calculate the histogram of each pixel.
Returns:
np.ndarray: The histogram for each pixel value
"""
def fit_function(self):
"""
Fit a Gaussian to the histogram of each pixel.
Returns:
np.ndarray: The fitted Gaussian parameters for each pixel.
"""
def calibrate_G0(self):
"""
Calibrate the G0 gain of the Jungfrau detector.
Returns:
np.ndarray: The calibrated G0 gain values for each pixel.
"""
self.calculate_bad_pixels_mask()
self.calculate_histogram()
self.fit_function()
## additional plot methods for visualizing the histogram and fitted Gaussian parameters can be added here - depens how fast not neccessary to compute on the fly
def main():
calibration_params = JungfrauCalibrationParameters()
calibration_params.pedestal_file_dir = Path("/mnt/sls_det_storage/jungfrau_calib/data/Module_708_Calib")
calibration_params.pedestal_file_prefix = "pedeG0_M708_2025-12-09_"
calibration_params.num_pedestals_g0 = 1000
calibration_params.num_pedestals_g1 = 1000
calibration_params.num_pedestals_g2 = 1000
python/src/JungfrauCalibrationParameters.py
-77
View File
@@ -1,77 +0,0 @@
from pathlib import Path
class JungfrauCalibrationParameters:
"""
A class to hold calibration parameters for the Jungfrau detector.
"""
@property
def pedestal_file_dir(self) -> Path:
return self._pedestal_file_dir
@pedestal_file_dir.setter
def pedestal_file_dir(self, filepath : Path):
if not filepath.exists():
raise ValueError(f"Pedestal file directory {filepath} does not exist.")
self._pedestal_file_dir = filepath
@property
def pedestal_g0_file_prefix(self) -> str:
return self._pedestal_g0_file_prefix
@pedestal_g0_file_prefix.setter
def pedestal_g0_file_prefix(self, file_prefix : str):
self._pedestal_g0_file_prefix = file_prefix
# TODO: add deprecated decorator
@property
def pedestal_file_prefix(self) -> str:
return self._pedestal_file_prefix
@pedestal_file_prefix.setter
def pedestal_file_prefix(self, file_prefix : str):
self._pedestal_file_prefix = file_prefix
@property
def num_pedestals_g0(self) -> int:
return self._num_pedestals_g0
@num_pedestals_g0.setter
def num_pedestals_g0(self, num_pedestals : int):
self._num_pedestals_g0 = num_pedestals
@property
def num_pedestals_g1(self) -> int:
return self._num_pedestals_g1
@num_pedestals_g1.setter
def num_pedestals_g1(self, num_pedestals : int):
self._num_pedestals_g1 = num_pedestals
@property
def num_pedestals_g2(self) -> int:
return self._num_pedestals_g2
@num_pedestals_g2.setter
def num_pedestals_g2(self, num_pedestals : int):
self._num_pedestals_g2 = num_pedestals
@property
def raw_file_dir(self) -> Path:
return self._raw_file_dir
@raw_file_dir.setter
def raw_file_dir(self, filepath : Path):
if not filepath.exists():
raise ValueError(f"Raw file directory {filepath} does not exist.")
self._raw_file_dir = filepath
@property
def raw_file_prefix(self) -> str:
return self._raw_file_prefix
@raw_file_prefix.setter
def raw_file_prefix(self, file_prefix : str):
self._raw_file_prefix = file_prefix
# TODO add a read_config method
python/src/helpers.py
-11
View File
@@ -1,11 +0,0 @@
import glob
from pathlib import Path
def get_first_file(path : Path, file_prefix : str):
"""
Get the first file with lowest index in directory that matches the given prefix.
"""
first_file = min(path.glob(str(Path)+f'{file_prefix}*', default=None))
if first_file is None:
raise ValueError(f"No files found in {path} with prefix {file_prefix}")
return first_file
src/BadChannels.cpp
-51
View File
@@ -1,51 +0,0 @@
#include "BadChannels.hpp"
#include "aare/calibration.hpp"
using namespace aare;
namespace jungfraucalibration
{
NDArray<bool, 2> CreateBadChannelPixelMask(JungfrauDataFile &pedestal_file, const size_t num_pedestals_g0, const size_t num_pedestals_g1, const size_t num_pedestals_g2)
{
auto pedestals_g0 = pedestal_file.read_n(num_pedestals_g0);
auto pedestals_g1 = pedestal_file.read_n(num_pedestals_g1);
auto pedestals_g2 = pedestal_file.read_n(num_pedestals_g2);
const size_t rows = pedestal_file.rows();
const size_t cols = pedestal_file.cols();
// get gain from each pixel - update mask
NDArray<bool, 2> bad_channel_mask({static_cast<ssize_t>(rows), static_cast<ssize_t>(cols)}, false); // bad channels pixel mask
size_t max_frames = std::max({num_pedestals_g0, num_pedestals_g1, num_pedestals_g2});
// TODO is this more efficient e.g. all three frames fit into cache instead of doing one file at the time?
for (size_t frame_idx = 0; frame_idx < max_frames; ++frame_idx)
{
for (size_t row = 0; row < rows; ++row)
{
for (size_t col = 0; col < cols; ++col)
{
// TODO: nicer to access element from frame directly instead of view()? What is the type?
if (frame_idx < num_pedestals_g0 && get_gain(pedestals_g0[frame_idx].view<uint16_t>()(row, col)) != 0)
{
bad_channel_mask(row, col) = true;
}
if (frame_idx < num_pedestals_g1 && get_gain(pedestals_g1[frame_idx].view<uint16_t>()(row, col)) != 1)
{
bad_channel_mask(row, col) = true;
}
if (frame_idx < num_pedestals_g2 && get_gain(pedestals_g2[frame_idx].view<uint16_t>()(row, col)) != 2)
{
bad_channel_mask(row, col) = true;
}
}
}
}
return bad_channel_mask;
}
} // namespace jungfraucalibration