garnet.reduction package#

Submodules#

garnet.reduction.background module#

class Normalization(plan)[source]#

Bases: object

Methods

combine

combine_parallel

get_file

normalize

normalize_parallel

symmetry_name

validate_params

combine(files)[source]#
static combine_parallel(plan, files)[source]#
get_file(file, ws='')[source]#
normalize()[source]#
static normalize_parallel(plan, runs, proc)[source]#
symmetry_name(file)[source]#
validate_params()[source]#

garnet.reduction.crystallography module#

garnet.reduction.data module#

class BaseDataModel(instrument_config)[source]#

Bases: object

Methods

add_UBW(ws, ub_file, projections[, run_number])

Attach sample UB and projection matrix to workspace

bin_Q_sample_to_hkl(md, peaks, hkl_histo[, ...])

Bin Q-sample data into a rough HKL histogram for a diagnostic check of whether the data is well aligned with the UB from a peaks workspace.

bin_events_in_Q(event_name, centers, ...)

Histogram the raw events about each of several Q-sample centers into a cube, straight from the events rather than by converting to MD: only the pixels whose lines Q = k u cross each cube, and only the momenta that can reach it, are read.

bin_in_Q(md, extents, bins, projections)

Histogram data into Q-space.

calculate_binning_from_bins(xmin, xmax, bins)

Determine the binning from the number of bins.

calculate_binning_from_step(xmin, xmax, step)

Determine the binning from step size.

combine_Q_sample(combine, merge)

Merge Q-sample workspaces into one.

combine_files(IPTS, process_map)

Combine runs together based on crystal orientation.

combine_histogram_files(files, merge[, n_proc])

Merge histogram files into a single accumulated workspace.

combine_histograms(ws, merge)

Add two histogram workspaces together.

combine_splits(md, log_name, log_vals, ...)

Remove all accumulation workspaces.

delete_workspace(ws)

Delete workspace.

divide_histograms(ws, num, den)

Divide two histogram workspaces.

extract_axis_info(ws)

Obtain the axis information from a histogram.

extract_bin_info(ws)

Obtain the bin information from a histogram.

extract_counts(ws)

Obtain the bin counts from a histogram.

file_names(IPTS, runs)

Complete file paths.

get_file_name_list(IPTS, runs)

Complete list of file paths.

get_hkl_limits(ws[, d_min])

The minimum and maximum Q-values.

get_min_max_values()

The minimum and maximum Q-values.

get_resolution_in_Q(lamda, two_theta)

Obtain the wavelength and detector-dependent Q-resolution

load_clear_UB(filename, ws[, run_number])

Load UB from file and replace.

load_histograms(filename, ws)

Load histograms file.

log_split_info(ws, log_name, log_limits, ...)

Generate split information for filtering events by log values.

merge_Q_sample(filenames, filename, merge)

Merge Q-sample files into one.

save_UB(filename, ws)

Save UB to file.

save_histograms(filename, ws[, sample_logs])

Save histograms file.

set_goniometer(ws)

Set the goniomter motor angles

subtract_histograms(ws, ws1, ws2)

Difference between two histograms.

update_raw_path(plan)

Set additional parameters for data file.

check_volume_preservation

clear_norm

get_counting_rate

slice_extents

slice_roi

update_logs_for_time

update_wavelength

workspace_exists

add_UBW(ws, ub_file, projections, run_number=None)[source]#

Attach sample UB and projection matrix to workspace

Parameters:
wsstr

Name of histogram to be added.

filenamestr

Name of UB file with extension .mat.

projectionslist

Axis projections.

bin_Q_sample_to_hkl(md, peaks, hkl_histo, bins=201, d_min=None)[source]#

Bin Q-sample data into a rough HKL histogram for a diagnostic check of whether the data is well aligned with the UB from a peaks workspace.

Parameters:
mdstr

Name of Q-sample MD workspace.

peaksstr

Name of peaks workspace providing the UB matrix.

hkl_histostr

Name of output HKL histogram workspace.

binsint, optional

Number of bins per dimension. The default is 201.

d_minfloat, optional

Minimum d-spacing to bound hkl by instead of the instrument Q-coverage. The default is None.

bin_events_in_Q(event_name, centers, half_widths, R, bins)[source]#

Histogram the raw events about each of several Q-sample centers into a cube, straight from the events rather than by converting to MD: only the pixels whose lines Q = k u cross each cube, and only the momenta that can reach it, are read.

Parameters:
event_namestr

Name of the events, in momentum, with “detectors” from preprocess_detectors.

centerslist

Q-sample center of each cube.

half_widthslist

Half the side of each cube.

Rarray

Goniometer matrix.

binsint

Number of bins along each side.

Returns:
boxeslist

For each cube, the counts and the Q-sample coordinates of the bin centers (each indexed as [i, j, k] along the three axes).

bin_in_Q(md, extents, bins, projections)[source]#

Histogram data into Q-space.

Parameters:
mdstr

3D Q-space data.

extentslist

Min/max pairs for each dimension.

binslist

Number of bins for each dimension.

projectionslist

Direction of projection for each dimension.

calculate_binning_from_bins(xmin, xmax, bins)[source]#

Determine the binning from the number of bins.

Parameters:
xminfloat

Minimum bin center.

xmaxfloat

Maximum bin center.

binsint

Number of bins.

Returns:
min_edgefloat

Minimum bin edge.

max_edgefloat

Maximum bin edge.

stepfloat

Bin step.

calculate_binning_from_step(xmin, xmax, step)[source]#

Determine the binning from step size.

Parameters:
xminfloat

Minimum bin center.

xmaxfloat

Maximum bin center.

stepfloat

Bin step.

Returns:
min_edgefloat

Minimum bin edge.

max_edgefloat

Maximum bin edge.

binsint

Number of bins.

check_volume_preservation(ws)[source]#
clear_norm(md)[source]#
combine_Q_sample(combine, merge)[source]#

Merge Q-sample workspaces into one.

Parameters:
combinelist

Name of Q-sample workspaces to be accumulated.

mergestr

Name of combined Q-sample workspace.

combine_files(IPTS, process_map)[source]#

Combine runs together based on crystal orientation.

combine_histogram_files(files, merge, n_proc=None)[source]#

Merge histogram files into a single accumulated workspace.

The first file is loaded as a full MDHistoWorkspace to provide the bin geometry, sample logs, and experiment info for the merge. The signal and error-squared arrays of the remaining files are read directly with h5py instead of being loaded through Mantid, since only those two arrays are needed for the sum, and are accumulated onto the first workspace’s arrays. Each remaining file’s own per-index “log_<n>” sample logs (see combine_splits) are read back the same way and reapplied, since they otherwise only survive from the first file.

Parameters:
fileslist

Histogram filenames to combine. The first is used as the base.

mergestr

Name of the combined workspace to create.

n_procint, optional

Number of worker processes used to read the remaining files. The default is None (serial reads).

combine_histograms(ws, merge)[source]#

Add two histogram workspaces together.

Parameters:
wsstr

Name of histogram to be added.

mergestr

Name of histogram to be accumulated.

combine_splits(md, log_name, log_vals, log_units, index, runs)[source]#

Remove all accumulation workspaces.

Parameters:
mdstr

Base name of histogram.

delete_workspace(ws)[source]#

Delete workspace.

Parameters:
wsstr

Workspace to remove.

divide_histograms(ws, num, den)[source]#

Divide two histogram workspaces.

Parameters:
wsstr

Name of resulting histogram.

numstr

Name of numerator histogram.

denstr

Name of denominator histogram.

extract_axis_info(ws)[source]#

Obtain the axis information from a histogram.

Parameters:
wsstr

Name of histogram.

Returns:
UB3x3-matrix

UB-matrix.

W3x3-matrix

Projection matrix.

titleslist

Axis names and units.

x0, x1, …array

Bin center coordinates.

extract_bin_info(ws)[source]#

Obtain the bin information from a histogram.

Parameters:
wsstr

Name of histogram.

Returns:
signalarray

Data signal.

errorarray

Data uncertanies.

x0, x1, …array

Dense bin center coordinates.

extract_counts(ws)[source]#

Obtain the bin counts from a histogram.

Parameters:
wsstr

Name of histogram.

Returns:
Countsarray

Number of events.

file_names(IPTS, runs)[source]#

Complete file paths.

Parameters:
IPTSint

Proposal number.

runslist, int

List of run number(s).

Returns:
filenamesstr

Comma separated filepaths.

get_counting_rate()[source]#
get_file_name_list(IPTS, runs)[source]#

Complete list of file paths.

Parameters:
IPTSint

Proposal number.

runslist, int

List of run number(s).

Returns:
filenameslist

List of filepaths.

get_hkl_limits(ws, d_min=None)[source]#

The minimum and maximum Q-values.

Parameters:
wsstr

Name of workspace with an oriented lattice.

d_minfloat, optional

Minimum d-spacing to bound hkl by instead of the instrument Q-coverage. The default is None.

Returns:
hkl_min_vals: list

Minumum hkl.

hkl_max_vals: list

Maximum hkl.

get_min_max_values()[source]#

The minimum and maximum Q-values.

Returns:
Q_min_vals: list

Minumum Q.

Q_max_vals: list

Maximum Q.

get_resolution_in_Q(lamda, two_theta)[source]#

Obtain the wavelength and detector-dependent Q-resolution

Parameters:
lamdafloat

Incident wavelength.

two_thetafloat

Detector scattering angle in degrees.

Returns:
dQarray

Q-resolution.

load_clear_UB(filename, ws, run_number=None)[source]#

Load UB from file and replace.

Parameters:
filenamestr

Name of UB file with extension .mat.

wsstr, optional

Name of data.

run_numberstr, optional

Run number to replace starred expression in filename.

load_histograms(filename, ws)[source]#

Load histograms file.

Parameters:
filenamestr

Name of peaks file with extension .nxs.

wsstr

Name of histogram to be added.

log_split_info(ws, log_name, log_limits, log_bins, log_ramp='Both')[source]#

Generate split information for filtering events by log values.

Parameters:
wsstr

Workspace with log time and values.

log_namestr

Name of the log.

log_limitslist, float

Min/max values of log (bin center).

log_binsint

Number of equally space bins.

log_rampstr

Filter log by ramp direction. The default is Both.

merge_Q_sample(filenames, filename, merge)[source]#

Merge Q-sample files into one.

Parameters:
filenameslist

Name of Q-sample filenames to be combined.

filename: str

Name of Q-sample filename to be saved.

mergestr

Name of Q-sample workspace to be accumulated.

save_UB(filename, ws)[source]#

Save UB to file.

Parameters:
filenamestr

Name of UB file with extension .mat.

wsstr, optional

Name of data.

save_histograms(filename, ws, sample_logs=False)[source]#

Save histograms file.

Parameters:
filenamestr

Name of peaks file with extension .nxs.

wsstr

Name of histogram to be added.

set_goniometer(ws)[source]#

Set the goniomter motor angles

Parameters:
wsstr, optional

Name of raw data.

slice_extents(UB, hkl)[source]#
slice_roi(md, UB, hkl)[source]#
subtract_histograms(ws, ws1, ws2)[source]#

Difference between two histograms.

Parameters:
wsstr

Name of resulting histogram.

ws1str

Name of first histogram.

ws2str

Name of second histogram.

update_logs_for_time(ws)[source]#
update_raw_path(plan)[source]#

Set additional parameters for data file.

Parameters:
plandict

Reduction plan.

update_wavelength(wl)[source]#
workspace_exists(ws)[source]#
DataModel(instrument_config)[source]#
class LaueData(instrument_config)[source]#

Bases: BaseDataModel

Methods

add_UBW(ws, ub_file, projections[, run_number])

Attach sample UB and projection matrix to workspace

apply_calibration(event_name, ...[, ...])

Apply detector calibration.

apply_mask(event_name, detector_mask[, save])

Apply detector mask.

bin_Q_sample_to_hkl(md, peaks, hkl_histo[, ...])

Bin Q-sample data into a rough HKL histogram for a diagnostic check of whether the data is well aligned with the UB from a peaks workspace.

bin_events_in_Q(event_name, centers, ...)

Histogram the raw events about each of several Q-sample centers into a cube, straight from the events rather than by converting to MD: only the pixels whose lines Q = k u cross each cube, and only the momenta that can reach it, are read.

bin_in_Q(md, extents, bins, projections)

Histogram data into Q-space.

calculate_binning_from_bins(xmin, xmax, bins)

Determine the binning from the number of bins.

calculate_binning_from_step(xmin, xmax, step)

Determine the binning from step size.

calculate_maximum_Q([scale])

Update maximum Q.

combine_Q_sample(combine, merge)

Merge Q-sample workspaces into one.

combine_files(IPTS, process_map)

Combine runs together based on crystal orientation.

combine_histogram_files(files, merge[, n_proc])

Merge histogram files into a single accumulated workspace.

combine_histograms(ws, merge)

Add two histogram workspaces together.

combine_splits(md, log_name, log_vals, ...)

Remove all accumulation workspaces.

convert_to_Q_lab(event_name, md_name[, ...])

Convert raw data to Q-lab.

convert_to_Q_sample(event_name, md_name[, ...])

Convert raw data to Q-sample.

convert_to_hkl(event_name, md_name[, ...])

Convert raw data to hkl.

crop_for_normalization(event_name[, tolerance])

Convert units to momentum and crop to momentum band.

delete_workspace(ws)

Delete workspace.

divide_histograms(ws, num, den)

Divide two histogram workspaces.

extract_axis_info(ws)

Obtain the axis information from a histogram.

extract_bin_info(ws)

Obtain the bin information from a histogram.

extract_counts(ws)

Obtain the bin counts from a histogram.

file_names(IPTS, runs)

Complete file paths.

filter_events(ws, run, runs)

Split workspaces according to log.

get_counting_rate(event_name)

Counting rate.

get_file_name_list(IPTS, runs)

Complete list of file paths.

get_hkl_limits(ws[, d_min])

The minimum and maximum Q-values.

get_min_max_values()

The minimum and maximum Q-values.

get_resolution_in_Q(lamda, two_theta)

Obtain the wavelength and detector-dependent Q-resolution

load_background(filename, event_name[, ...])

Load a background file and scale to data.

load_clear_UB(filename, ws[, run_number])

Load UB from file and replace.

load_data(event_name, IPTS, runs[, time_cut])

Load raw data into time-of-flight vs counts.

load_generate_normalization(vanadium_file, ...)

Load a vanadium file and generate normalization data.

load_histograms(filename, ws)

Load histograms file.

log_split_info(ws, log_name, log_limits, ...)

Generate split information for filtering events by log values.

mask_to_bank(event_name, bank_name)

Mask to bank

merge_Q_sample(filenames, filename, merge)

Merge Q-sample files into one.

normalize_to_hkl(md, projections, extents, bins)

Normalize to binned hkl.

preprocess_detector_banks(bank)

Generate detector coordinates for each bank.

preprocess_detectors([ws])

Generate detector coordinates.

save_UB(filename, ws)

Save UB to file.

save_histograms(filename, ws[, sample_logs])

Save histograms file.

set_goniometer(ws)

Set the goniomter motor angles

subtract_histograms(ws, ws1, ws2)

Difference between two histograms.

update_raw_path(plan)

Set additional parameters for data file.

absorption_correction

check_volume_preservation

clear_norm

convert_to_TOF

get_monitor_ratio

grouping_list

slice_extents

slice_roi

subtract_background

update_logs_for_time

update_wavelength

workspace_exists

absorption_correction(event_name, shapestl, angles, material)[source]#
apply_calibration(event_name, detector_calibration, tube_calibration=None, goniometer_calibration=None)[source]#

Apply detector calibration.

Parameters:
event_namestr

Name of raw event_name data.

detector_calibrationstr

Detector calibration as either .xml or .DetCal.

tube_calibrationstr, optional

CORELLI-only tube calibration. The default is None.

goniometer_calibrationstr, optional

Goniometer calibration as .xml. The default is None.

apply_mask(event_name, detector_mask, save=False)[source]#

Apply detector mask.

Parameters:
event_namestr

Name of raw event_name data.

detector_maskstr

Detector mask as .xml.

calculate_maximum_Q(scale=1)[source]#

Update maximum Q.

convert_to_Q_lab(event_name, md_name, lorentz_corr=False)[source]#

Convert raw data to Q-lab.

Parameters:
event_namestr

Name of raw event name data.

md_namestr

Name of Q-sample workspace.

lorentz_corrbool, optional

Apply Lorentz correction. The default is False.

convert_to_Q_sample(event_name, md_name, lorentz_corr=False, preproc_dets=None)[source]#

Convert raw data to Q-sample.

Parameters:
event_namestr

Name of raw event name data.

md_namestr

Name of Q-sample workspace.

lorentz_corrbool, optional

Apply Lorentz correction. The default is False.

preproc_dets: str, optional

Preprocess detector information workspace.

convert_to_TOF(event_name)[source]#
convert_to_hkl(event_name, md_name, lorentz_corr=False, preproc_dets=None)[source]#

Convert raw data to hkl.

Parameters:
event_namestr

Name of raw event name data.

md_namestr

Name of Q-sample workspace.

lorentz_corrbool, optional

Apply Lorentz correction. The default is False.

preproc_dets: str, optional

Preprocess detector information workspace.

crop_for_normalization(event_name, tolerance=0.01)[source]#

Convert units to momentum and crop to momentum band.

event_namestr

Name of raw event data.

factorfloat

Compression tolerance.

filter_events(ws, run, runs)[source]#

Split workspaces according to log. If no splitting, index according to run number.

Parameters:
wsstr

Workspace to split.

runint

Run number.

runslist

All run numbers

Returns:
indiceslist, str

Split index.

workspaceslist, float

Split workspace.

get_counting_rate(event_name)[source]#

Counting rate.

Parameters:
event_namestr

Name of raw event_name data.

Returns:
cntrtfloat

Counting rate in events per proton charge.

get_monitor_ratio()[source]#
grouping_list(ws, cols, rows, lite_c, lite_r, group_c, group_r)[source]#
load_background(filename, event_name, detector_calibration=None, tube_calibration=None, save=False, remove=False)[source]#

Load a background file and scale to data.

Parameters:
filenamestr

Background file.

event_namestr

Name of raw event data.

removebool

Delete raw background workspace.

load_data(event_name, IPTS, runs, time_cut=None)[source]#

Load raw data into time-of-flight vs counts.

Parameters:
event_namestr

Name of raw event_name data.

IPTSint

Proposal number.

runslist, int

List of run number(s).

time_cut: float, optional

Time cut off for faster loading.

load_generate_normalization(vanadium_file, flux_file, save=False)[source]#

Load a vanadium file and generate normalization data.

Parameters:
vanadium_filestr

Solid angle file.

flux_filestr

Flux file.

mask_to_bank(event_name, bank_name)[source]#

Mask to bank

Parameters:
event_namestr

Name of raw event_name data.

bank_namestr

Bank name.

normalize_to_hkl(md, projections, extents, bins, symmetry=None, zero=False)[source]#

Normalize to binned hkl.

Parameters:
mdstr

3D Q-sample data.

projectionslist

Projection axis vectors.

extentslist

Min/max pairs defining the bin center limits.

binslist

Number of bins.

symmetrystr, optional

Laue point group. The default is None.

zerobool, optional

Discard the normalized result, keeping data/norm at zero. Used when the projection axes could not be derived (e.g. a MillerIndex that is not elastically reachable). The default is False.

preprocess_detector_banks(bank)[source]#

Generate detector coordinates for each bank.

Parameters:
bankstr

Bank-only event data.

preprocess_detectors(ws=None)[source]#

Generate detector coordinates.

Parameters:
event_namestr

Workspace with instrument data.

subtract_background(event_name)[source]#
class MonochromaticData(instrument_config)[source]#

Bases: BaseDataModel

Methods

add_UBW(ws, ub_file, projections[, run_number])

Attach sample UB and projection matrix to workspace

bin_Q_sample_to_hkl(md, peaks, hkl_histo[, ...])

Bin Q-sample data into a rough HKL histogram for a diagnostic check of whether the data is well aligned with the UB from a peaks workspace.

bin_events_in_Q(event_name, centers, ...)

Histogram the raw events about each of several Q-sample centers into a cube, straight from the events rather than by converting to MD: only the pixels whose lines Q = k u cross each cube, and only the momenta that can reach it, are read.

bin_in_Q(md, extents, bins, projections)

Histogram data into Q-space.

calculate_binning_from_bins(xmin, xmax, bins)

Determine the binning from the number of bins.

calculate_binning_from_step(xmin, xmax, step)

Determine the binning from step size.

combine_Q_sample(combine, merge)

Merge Q-sample workspaces into one.

combine_files(IPTS, process_map)

Combine runs together based on crystal orientation.

combine_histogram_files(files, merge[, n_proc])

Merge histogram files into a single accumulated workspace.

combine_histograms(ws, merge)

Add two histogram workspaces together.

combine_splits(md, log_name, log_vals, ...)

Remove all accumulation workspaces.

convert_to_Q_sample(histo_name, md_name[, ...])

Convert raw data to Q-sample.

delete_workspace(ws)

Delete workspace.

divide_histograms(ws, num, den)

Divide two histogram workspaces.

extract_axis_info(ws)

Obtain the axis information from a histogram.

extract_bin_info(ws)

Obtain the bin information from a histogram.

extract_counts(ws)

Obtain the bin counts from a histogram.

file_names(IPTS, runs)

Complete file paths.

get_file_name_list(IPTS, runs)

Complete list of file paths.

get_hkl_limits(ws[, d_min])

The minimum and maximum Q-values.

get_min_max_values()

The minimum and maximum Q-values.

get_resolution_in_Q(lamda, two_theta)

Obtain the wavelength and detector-dependent Q-resolution

load_background(filename[, histo_name])

Load a background file and scale to data.

load_clear_UB(filename, ws[, run_number])

Load UB from file and replace.

load_data(histo_name, IPTS, runs)

Load raw data into detector counts vs rotation index.

load_generate_normalization(filename[, ...])

Load a vanadium file and generate normalization data.

load_histograms(filename, ws)

Load histograms file.

log_split_info(ws, log_name, log_limits, ...)

Generate split information for filtering events by log values.

merge_Q_sample(filenames, filename, merge)

Merge Q-sample files into one.

normalize_to_hkl(ws, projections, extents, bins)

Normalizae to binned hkl.

save_UB(filename, ws)

Save UB to file.

save_histograms(filename, ws[, sample_logs])

Save histograms file.

set_goniometer(ws)

Set the goniomter motor angles

subtract_histograms(ws, ws1, ws2)

Difference between two histograms.

update_raw_path(plan)

Set additional parameters for data file.

check_volume_preservation

clear_norm

get_counting_rate

slice_extents

slice_roi

update_logs_for_time

update_wavelength

workspace_exists

convert_to_Q_sample(histo_name, md_name, lorentz_corr=False)[source]#

Convert raw data to Q-sample.

Parameters:
histo_namestr

Name of raw histogram data.

md_namestr

Name of Q-sample workspace.

lorentz_corrbool, optional

Apply Lorentz correction. The default is False.

load_background(filename, histo_name=None)[source]#

Load a background file and scale to data.

Parameters:
filenamestr

Background file.

histo_namestr

Name of raw histogram data.

load_data(histo_name, IPTS, runs)[source]#

Load raw data into detector counts vs rotation index.

Parameters:
histo_namestr

Name of raw histogram data.

IPTSint

Proposal number.

runslist, int

List of run number(s).

groupingstr, optional

Options for grouping pixels.

load_generate_normalization(filename, histo_name=None)[source]#

Load a vanadium file and generate normalization data. Provided a histogram workspace name, generate corresponding shape.

Parameters:
filenamestr

Vanadium file.

histo_namestr, optional

Name of raw histogram data.

normalize_to_hkl(ws, projections, extents, bins, symmetry=None)[source]#

Normalizae to binned hkl.

Parameters:
wsstr

3D detector counts vs rotation index data.

projectionslist

Projection axis vectors.

extentslist

Min/max pairs defining the bin center limits.

binslist

Number of bins.

symmetrystr, optional

Laue point group. The default is None.

garnet.reduction.ellipsoid module#

class PeakEllipsoid[source]#

Bases: object

Methods

S_matrix(r0, r1, r2, u0, u1, u2)

Ellipsoid covariance-like matrix S = U diag(r0^2, r1^2, r2^2) U^T.

U_matrix(u0, u1, u2)

Ellipsoid orientation rotation matrix from an axis-angle rotation vector.

background_profile(x0, x1, x2, d, n, b, m)

Projected 1D background intensity and Poisson uncertainty for each axis.

calculate_intensity(A, H, r0, r1, r2, u0, u1, u2)

Total fitted intensity for a mode: A * gaussian_integral(inv_S, mode).

centroid_inverse_covariance(c0, c1, c2, r0, ...)

Package the peak center and inverse ellipsoid covariance from scalar parameters.

coerce_weight(target[, weight])

Normalize a weight argument into an array matching the target's shape.

coerce_weights(targets[, weights])

Apply coerce_weight element-wise across parallel lists of targets and weights.

collect_mode_fit_metrics(x0, x1, x2, c, ...)

Fitted profile, chi-squared, and box-sum intensity for each requested mode.

copy_combine()

Snapshot self.params into self.combine_params for estimate_envelope to resume from.

counts_to_intensity_uncertainty(d, n)

Convert raw counts and normalization into normalized intensity and its uncertainty.

data_norm(d, n, v[, rel_err])

Normalize raw counts by monitor counts and estimate their uncertainty.

det_S(r0, r1, r2, u0, u1, u2)

Determinant of the ellipsoid covariance-like matrix S, i.e. (r0*r1*r2)**2.

ellipsoid_covariance(inv_S[, mode, perc])

Rescale (a projection of) inv_S to an arbitrary confidence-contour percentile.

estimate_center_weighted(d, n, kernel[, ...])

Estimate a 1D peak-center index via matched-filter correlation with a kernel.

estimate_envelope(x0, x1, x2, d_int, n_int, wgt)

Build all mode projections, seed initial amplitudes, and iteratively fit the ellipsoid.

estimate_intensity(x0, x1, x2, c, inv_S, ...)

Background-subtracted box-sum intensity and uncertainty within the 1-sigma peak region (I_1d/I_2d/I_3d).

estimate_peak_strength()

Fitted 1D-mode amplitude and background values as float arrays.

extract_amplitude_background()

Collect the fitted 1D-mode amplitude and background parameters.

extract_intensity(d, n, pk, bkg[, vol_frac, ...])

Normalization-aware background-subtracted box-sum intensity (the I_ell estimator).

extract_raw_intensity(counts, pk, bkg)

Simple background-subtracted box-sum intensity from raw counts, with Poisson error.

extract_result(args_1d, args_2d, args_3d, xmod)

Finalize the fit: compute per-mode metrics and populate the post-fit result attributes.

fit(x0_prof, x1_proj, x2_proj, d, n, dx, ...)

Top-level fit entry point: crop data to the valid detector region and fit the ellipsoid.

fitted_profile(x0, x1, x2, d, n, c, S[, p, ...])

Iteratively matched-filter fit a 1D Gaussian profile along axis 0 (the I_prof estimator).

gaussian(x0, x1, x2, c, inv_S[, mode])

Unnormalized (peak-value-1) Gaussian profile for a given mode.

gaussian_integral(inv_S[, mode])

Normalization integral sqrt((2*pi)**k * det) of the unit-peak Gaussian for a given mode.

gaussian_integral_jac_S(inv_S, d_inv_S[, mode])

Derivative of gaussian_integral with respect to a shape parameter group.

gaussian_jac_S(x0, x1, x2, c, inv_S, d_inv_S)

Derivative of the unnormalized gaussian profile with respect to a shape parameter group.

gaussian_jac_c(x0, x1, x2, c, inv_S[, mode])

Derivative of the unnormalized gaussian profile with respect to the center c0, c1, c2.

integrate(x0, x1, x2, d, n, c, S[, b, m])

Final post-fit intensity integration combining box-sum, ellipsoid, and profile estimators.

inv_S_deriv_domega(r0, r1, r2, domega[, delta])

Finite-difference derivatives of inv_S with respect to domega_x, domega_y, domega_z.

inv_S_deriv_r(r0, r1, r2, u0, u1, u2)

Analytic derivatives of inv_S with respect to each radius r0, r1, r2.

inv_S_deriv_size_shape(r0, r1, r2, u0, u1, u2)

Analytic derivatives of inv_S with respect to log_size, shape_1, shape_2.

inv_S_matrix(r0, r1, r2, u0, u1, u2)

Inverse ellipsoid covariance-like matrix inv_S = U diag(1/r0^2, 1/r1^2, 1/r2^2) U^T.

jacobian(params, args_1d, args_2d, args_3d)

Full stacked Jacobian matching the residual ordering of residual.

jacobian_1d(params, x0, x1, x2, ds, ns[, ...])

Column-stacked Poisson-deviance Jacobians across all three 1D projections.

jacobian_2d(params, x0, x1, x2, ds, ns[, ...])

Column-stacked Poisson-deviance Jacobians across all three 2D projections.

jacobian_3d(params, x0, x1, x2, d, n[, w, ...])

Poisson-deviance Jacobian for the full 3D volume.

jacobian_mode_poisson(params, x0, x1, x2, d, ...)

Analytic Jacobian of residual_mode_poisson with respect to the varying fit parameters.

matched_filter(d, n, kernel, mask[, bkg_meas])

Matched-filter Poisson maximum-likelihood fit of an intensity and background level.

missing_support_fraction(x0, x1, x2, c, inv_S, n)

Fraction of the peak's total Gaussian mass not captured by valid voxels.

mode_model_counts(params, x0, x1, x2, d, n, ...)

Poisson model counts mu = n*(A*gaussian + B) [+ bkg] for one mode.

normalize(x0, x1, x2, d, n, w[, mode])

Project raw data onto a mode and normalize it into a weighted intensity/uncertainty.

peak_roi(x0, x1, x2, c, S, scale[, p])

Build peak/background/combined masks and a p-confidence Gaussian kernel around an ellipsoid.

poisson_deviance_fit(x0, x1, x2, c, inv_S, ...)

Reduced Poisson deviance goodness-of-fit statistic within 1-sigma of the peak center.

poisson_deviance_residual_factor(d, mu[, w, eps])

Signed Poisson deviance residual r = sign(mu-d)*sqrt(2*(mu-d+d*log(d/mu))) and dr/dmu.

prior_jacobian(params)

Analytic Jacobian of prior_residual with respect to the varying fit parameters.

prior_residual(params)

Whitened regularization residuals pulling the fit toward the resolution-model prior.

profile_fit(d, n, kernel, mask[, bkg_meas, ...])

Properly weighted Poisson profile fit of intensity and background from a fixed kernel.

profile_project(x0, x1, x2, d, n, w[, mode])

Project raw counts, normalization, variance proxy, and weight onto a given mode.

project_background(bkg, mode)

Sum a 3D background-count array onto a given mode's projection.

quick_gaussian(x0, x1, x2, d, n[, mode, b, m])

Quick amplitude/background estimate for a mode, seeding self.params.

residual(params, args_1d, args_2d, args_3d)

Full stacked residual vector combining all modes and the prior.

residual_1d(params, x0, x1, x2, ds, ns[, ...])

Concatenated Poisson-deviance residuals across all three 1D projections.

residual_2d(params, x0, x1, x2, ds, ns[, ...])

Concatenated Poisson-deviance residuals across all three 2D projections.

residual_3d(params, x0, x1, x2, d, n[, w, ...])

Poisson-deviance residual vector for the full 3D volume.

residual_mode_poisson(params, x0, x1, x2, d, ...)

Flattened Poisson-deviance residual vector for one mode, for use in least-squares fitting.

shape_from_params(params)

Reconstruct absolute radii and orientation rotation vector from the size/shape/orientation parameters.

subtract_profile(d, n[, perc])

Subtract a per-axis-0-slice background percentile from the raw counts.

sweep(args_1d, args_2d, args_3d[, protocol, ...])

Run one lmfit/scipy.optimize.least_squares minimization pass, accepted only if SNR improves.

uniform_mode_weights(ys, es[, val])

Build uniform per-mode weight arrays scaled by the total valid point count.

update_constraints(x0, x1, x2, dx)

(Re)initialize self.params bounds and starting values from the data extent.

update_estimate(shape[, sigma_c])

Seed self.params and the resolution-model prior from an initial ellipsoid estimate.

vech6(M)

Half-vectorize a symmetric 3x3 matrix into its 6 unique entries.

voxel_volume(x0, x1, x2)

Volume of a single voxel, the product of the per-axis voxel spacings.

voxels(x0, x1, x2)

Voxel spacing along each axis, from adjacent-element differences of the meshgrid.

S_matrix(r0, r1, r2, u0, u1, u2)[source]#

Ellipsoid covariance-like matrix S = U diag(r0^2, r1^2, r2^2) U^T.

Parameters:
r0, r1, r2float

Principal radii along the rotated axes.

u0, u1, u2float

Orientation rotation vector components.

Returns:
Sndarray of shape (3, 3)

Ellipsoid covariance-like matrix.

U_matrix(u0, u1, u2)[source]#

Ellipsoid orientation rotation matrix from an axis-angle rotation vector.

Parameters:
u0, u1, u2float

Axis-angle rotation vector components.

Returns:
Undarray of shape (3, 3)

Proper orthogonal rotation matrix.

background_profile(x0, x1, x2, d, n, b, m)[source]#

Projected 1D background intensity and Poisson uncertainty for each axis.

Crops inputs to the fit region (self._fit_crop) so coordinates match self.best_prof, rescales the background from monitor to signal units by the median ratio m/n, then projects onto each 1D axis normalized by projected normalization counts.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates (uncropped).

dndarray

Raw observed event counts, 3D array (uncropped).

nndarray

Normalization (monitor/solid-angle) counts, same shape as d (uncropped).

bndarray or None

Measured background counts, same shape as d (uncropped).

mndarray or None

Measured background monitor counts, same shape as d (uncropped).

Returns:
resultlist of tuple or None

3 tuples (x, y_bkg, e_bkg), one per axis (1d_0, 1d_1, 1d_2). None if b or m is None, or the median monitor scale factor is not finite and positive.

calculate_intensity(A, H, r0, r1, r2, u0, u1, u2, mode='3d')[source]#

Total fitted intensity for a mode: A * gaussian_integral(inv_S, mode).

Parameters:
Afloat

Fitted Gaussian amplitude for this mode.

Hfloat

Unused (reserved).

r0, r1, r2float

Principal radii along the rotated axes.

u0, u1, u2float

Orientation rotation vector components.

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

Returns:
intensityfloat

Total integrated intensity A * gaussian_integral(…).

centroid_inverse_covariance(c0, c1, c2, r0, r1, r2, u0, u1, u2)[source]#

Package the peak center and inverse ellipsoid covariance from scalar parameters.

Parameters:
c0, c1, c2float

Peak center coordinates.

r0, r1, r2float

Principal radii along the rotated axes.

u0, u1, u2float

Orientation rotation vector components.

Returns:
cndarray of shape (3,)

Peak center [c0, c1, c2].

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

coerce_weight(target, weight=None)[source]#

Normalize a weight argument into an array matching the target’s shape.

Parameters:
targetarray_like

Array whose shape the returned weight array should match.

weightNone, scalar, or array_like, optional

None gives all-ones; a scalar is broadcast; array-like is coerced to a float ndarray as-is.

Returns:
weightndarray

Float array of weights matching target’s shape.

coerce_weights(targets, weights=None)[source]#

Apply coerce_weight element-wise across parallel lists of targets and weights.

Parameters:
targetssequence of array_like

Arrays whose shapes the returned weights should match.

weightsNone or sequence, optional

If None, each target gets an all-ones weight; otherwise a same-length sequence, each element passed to coerce_weight.

Returns:
weightslist of ndarray

One coerced weight array per target.

collect_mode_fit_metrics(x0, x1, x2, c, inv_S, mode_data, bkg_offset=None)[source]#

Fitted profile, chi-squared, and box-sum intensity for each requested mode.

For each mode, evaluates A*gaussian + B (plus any bkg_offset), masks invalid points, and computes the reduced Poisson deviance (poisson_deviance_fit) and box-sum intensity/uncertainty (estimate_intensity).

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

mode_datadict

Mode string -> (y, e) observed normalized intensity and uncertainty.

bkg_offsetdict, optional

Mode string -> additional per-voxel background-rate offset; modes absent get no offset.

Returns:
metricsdict

Mode string -> [I0, s0, chi2, fit_triplet], where I0/s0 are the box-sum intensity/uncertainty, chi2 the reduced Poisson deviance, and fit_triplet = (y_fit, y, e).

copy_combine()[source]#

Snapshot self.params into self.combine_params for estimate_envelope to resume from.

counts_to_intensity_uncertainty(d, n)[source]#

Convert raw counts and normalization into normalized intensity and its uncertainty.

Uncertainty uses the sqrt(d + 1) Poisson approximation (avoids zero uncertainty at d == 0).

Parameters:
dndarray

Raw observed event counts.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

Returns:
yndarray

Normalized intensity d / n, NaN where n <= 0.

endarray

Uncertainty sqrt(clip(d, 0, None) + 1) / n, NaN where n <= 0.

data_norm(d, n, v, rel_err=30)[source]#

Normalize raw counts by monitor counts and estimate their uncertainty.

Voxels with non-positive or non-finite normalization are masked to NaN in-place in d, n, and v before normalizing.

Parameters:
dndarray

Raw observed event counts. Modified in place.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d. Modified in place.

vndarray

Variance proxy (e.g. clipped counts), same shape as d. Modified in place.

rel_errfloat, optional

Percentile of v used as a variance floor. Default 30.

Returns:
y_intndarray

Normalized intensity d / n.

e_intndarray

Estimated uncertainty sqrt(v + percentile(v, rel_err)) / n.

det_S(r0, r1, r2, u0, u1, u2)[source]#

Determinant of the ellipsoid covariance-like matrix S, i.e. (r0*r1*r2)**2.

Rotation-invariant, so computed directly from the radii.

Parameters:
r0, r1, r2float

Principal radii along the rotated axes.

u0, u1, u2float

Unused (det(S) does not depend on orientation).

Returns:
detfloat

Determinant of S.

ellipsoid_covariance(inv_S, mode='3d', perc=99.7)[source]#

Rescale (a projection of) inv_S to an arbitrary confidence-contour percentile.

inv_S is on the 99.7%-contour convention; this converts it (or the marginal/slice selected by mode) to the inverse-variance for percentile perc of the corresponding chi-squared distribution.

Parameters:
inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

percfloat, optional

Target confidence-contour percentile (0-100). Default 99.7.

Returns:
inv_varndarray or float

Rescaled inverse-variance: (3, 3) or (2, 2) matrix for “3d”/ “2d_*” modes, scalar for “1d_*” modes.

estimate_center_weighted(d, n, kernel, frac=0.9, min_weight=1e-12)[source]#

Estimate a 1D peak-center index via matched-filter correlation with a kernel.

Cross-correlates the background-subtracted counts with kernel (normalized by sqrt(correlate(n^2, kernel^2))) and returns the index of maximum correlation, falling back to the array midpoint if that maximum falls in the outer eighths (an unreliable edge estimate).

Parameters:
dndarray

1D raw observed event counts.

nndarray

1D normalization (monitor/solid-angle) counts, same shape as d.

kernelndarray

1D correlation kernel (e.g. a Gaussian template).

fracfloat, optional

Unused, reserved for future use. Default 0.9.

min_weightfloat, optional

Minimum correlated normalization-squared value to accept a correlation value. Default 1e-12.

Returns:
indexint

Estimated integer index of the peak center along d.

estimate_envelope(x0, x1, x2, d_int, n_int, wgt, report_fit=False, b=None, m=None)[source]#

Build all mode projections, seed initial amplitudes, and iteratively fit the ellipsoid.

Restores self.combine_params if available, projects the data onto all seven modes, quick-estimates each mode’s amplitude/background (quick_gaussian), converts a measured background into per-mode background-count arrays, then runs a multiresolution schedule of sweep calls: geometry parameters are unfrozen in stages (amplitude/background only, then + center, + size, + shape, + orientation). Since the 1D/2D modes are projected along the ellipsoid’s own principal axes, they are clean (separable) marginals of the 3D Gaussian rather than a correlated view of it, so the 1D/2D auxiliary residual weights (alpha_1d/alpha_2d, relative to self.mode_weights_1d/ self.mode_weights_2d) stay at full strength throughout, instead of being wound down to avoid biasing the 3D fit.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

d_intndarray

Raw observed event counts, 3D array.

n_intndarray

Normalization (monitor/solid-angle) counts, same shape as d_int.

wgtndarray

Per-voxel fit weight, same shape as d_int.

report_fitbool, optional

Print an lmfit fit report after each sweep call. Default False.

bndarray, optional

Measured background counts, same shape as d_int.

mndarray, optional

Measured background monitor counts, same shape as d_int.

Returns:
args_1dlist

[x0, x1, x2, d1d, n1d, w1d, bkg_1d] for the 1D-mode residual/Jacobian calls.

args_2dlist

Analogous list for the 2D-mode calls.

args_3dlist

Analogous list for the 3D-mode calls.

None instead of this tuple if the data shape is too small (any axis < 3) or any mode’s quick_gaussian estimate fails.

estimate_intensity(x0, x1, x2, c, inv_S, y_fit, y, e, mode='3d', bkg_offset=None)[source]#

Background-subtracted box-sum intensity and uncertainty within the 1-sigma peak region (I_1d/I_2d/I_3d).

Sums the background-subtracted normalized intensity over voxels with Mahalanobis distance <= 1, scaled by voxel volume. Distinct from the I_ell and I_prof estimators.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

y_fitndarray

Unused (kept to match the poisson_deviance_fit signature).

yndarray

Observed normalized intensity for this mode.

endarray

Uncertainty on y, same shape as y.

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

bkg_offsetndarray, optional

Additional per-voxel background-rate offset to subtract alongside the fitted B<mode> level.

Returns:
Ifloat

Background-subtracted, volume-scaled summed intensity.

sigfloat

Propagated uncertainty on I.

estimate_peak_strength()[source]#

Fitted 1D-mode amplitude and background values as float arrays.

Returns:
amplitudendarray of shape (3,)

Fitted A1d_0, A1d_1, A1d_2 values.

backgroundndarray of shape (3,)

Fitted B1d_0, B1d_1, B1d_2 values.

extract_amplitude_background()[source]#

Collect the fitted 1D-mode amplitude and background parameters.

Returns:
amplitudesndarray of shape (3,)

A1d_0, A1d_1, A1d_2 lmfit.Parameter objects.

backgroundsndarray of shape (3,)

B1d_0, B1d_1, B1d_2 lmfit.Parameter objects.

extract_intensity(d, n, pk, bkg, vol_frac=1.0, bkg_meas=None)[source]#

Normalization-aware background-subtracted box-sum intensity (the I_ell estimator).

Sums raw counts and normalization over the peak and background-shell regions (positive-normalization voxels only). If bkg_meas is supplied it’s used directly (assumed exact) instead of the shell; otherwise shell counts are scaled by the peak/shell voxel-count ratio, as in extract_raw_intensity but additionally normalized by monitor counts. Final intensity/uncertainty are rescaled by the peak voxel count and by vol_frac to correct for ellipsoid volume clipped by the detector edge.

Parameters:
dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

pkndarray of bool

Peak region mask, same shape as d.

bkgndarray of bool

Background shell region mask, same shape as d.

vol_fracfloat, optional

Fraction (0-1] of the theoretical ellipsoid volume covered by valid voxels; divides intens/sig to correct edge clipping. Default 1.0.

bkg_measndarray, optional

Measured/fitted background counts, same shape as d, used in place of the shell if provided (assumed zero variance).

Returns:
intensfloat

Background-subtracted, volume- and edge-corrected intensity.

sigfloat

Propagated uncertainty; -inf if variance is non-positive.

b, b_errfloat

Estimated background rate (per monitor count) and its uncertainty.

volfloat

Number of voxels in the peak region.

vol_fracfloat

The vol_frac argument, echoed back.

vol_pkint

Same as vol, as an integer.

pk_cnts, pk_normfloat

Summed raw counts and normalization in the peak region (NaN if zero).

bkg_cnts, bkg_normfloat

Summed background counts/normalization used (from bkg_meas or the shell).

ratiofloat

Peak/shell voxel-count ratio used to scale shell counts (1.0 when bkg_meas is used).

extract_raw_intensity(counts, pk, bkg)[source]#

Simple background-subtracted box-sum intensity from raw counts, with Poisson error.

Subtracts the background counts, scaled by the peak/background voxel-count ratio, from the peak-region raw counts.

Parameters:
countsndarray

Raw observed event counts, 3D array. Infinite values treated as NaN.

pkndarray of bool

Peak region mask, same shape as counts.

bkgndarray of bool

Background region mask, same shape as counts.

Returns:
intensfloat

Background-subtracted summed raw intensity.

sigfloat

Propagated Poisson uncertainty; inf if variance is non-positive.

extract_result(args_1d, args_2d, args_3d, xmod)[source]#

Finalize the fit: compute per-mode metrics and populate the post-fit result attributes.

Converts raw counts to normalized intensity/uncertainty for every mode and calls collect_mode_fit_metrics to populate self.reddev, self.intensity, self.sigma, self.fit_metrics, and the missing-support attributes. Also derives the final ellipsoid (center, radii, axes) from S = inv(inv_S), shifts the center back by xmod, and populates self.peak_pos, self.best_fit, self.best_prof, self.best_bkg_prof, and self.best_proj.

Parameters:
args_1dsequence

(x0, x1, x2, d1d, n1d, w1d, bkg_1d, …) where d1d/n1d are 3-tuples of per-axis 1D-projected counts/normalization and bkg_1d is None or a matching 3-tuple of background counts.

args_2dsequence

Analogous to args_1d for the three 2D projections.

args_3dsequence

Analogous to args_1d for the full 3D volume.

xmodsequence of float

(c0mod, c1mod, c2mod) offset subtracted from x0, x1, x2 before fitting; added back to the reported center and coordinates.

Returns:
c0, c1, c2float

Fitted peak center (c0 with xmod added back).

r0, r1, r2float

Principal radii (1-sigma) of the fitted ellipsoid.

v0, v1, v2ndarray of shape (3,)

Principal axis vectors of the fitted ellipsoid.

Returns None instead of this tuple if the optimized inverse covariance is not positive definite.

fit(x0_prof, x1_proj, x2_proj, d, n, dx, xmod, voxel_weights, b=None, m=None)[source]#

Top-level fit entry point: crop data to the valid detector region and fit the ellipsoid.

Converts counts to normalized intensity, rejects too-small or all-invalid input, crops arrays to the bounding box of voxels with positive normalization, then calls estimate_envelope to perform the iterative fit, snapshotting the result via copy_combine on success. Logs elapsed time and outcome for each early-exit path or success.

Parameters:
x0_prof, x1_proj, x2_projndarray

3D meshgrid coordinates along axes 0, 1, 2, original (not peak-centered) frame.

dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

dxobject

Unused (voxel spacing is recomputed via voxels).

xmodsequence of float

(c0mod, c1mod, c2mod) offset subtracted from x0_prof, x1_proj, x2_proj for the fit’s internal coordinates (the ideal Bragg position in the projection frame); stored as self._fit_xmod.

voxel_weightsndarray

Per-voxel fit weight, same shape as d.

bndarray, optional

Measured background counts, same shape as d.

mndarray, optional

Measured background monitor counts, same shape as d.

Returns:
weightstuple or None

The (args_1d, args_2d, args_3d) tuple from estimate_envelope on success; None on any early-exit condition (too small, insufficient valid data, all-invalid crop, an exception in estimate_envelope, or estimate_envelope returning None).

fitted_profile(x0, x1, x2, d, n, c, S, p=0.997, eta=0.5, bkg_meas=None)[source]#

Iteratively matched-filter fit a 1D Gaussian profile along axis 0 (the I_prof estimator).

Alternates matched_filter fits of a Gaussian kernel with re-estimating sigma from the weighted second moment of the net profile, damped by eta and floored at the voxel size. Initial sigma is sqrt(S[0,0]) rescaled from the 99.7%-contour convention using the 1D chi-squared quantile at p.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

csequence of float

Peak center (c0, c1, c2).

Sndarray of shape (3, 3)

Ellipsoid covariance-like matrix (99.7%-contour convention).

pfloat, optional

Confidence-contour probability (0-1) for the initial sigma and the peak/background split. Default 0.997.

etafloat, optional

Damping factor (0-1) for updating sigma between iterations. Default 0.5.

bkg_measndarray, optional

Measured/fitted background counts, 3D array, same shape as d, summed onto the axis-0 profile.

Returns:
I, I_errfloat

Fitted matched-filter intensity and its uncertainty.

bfloat

Fitted background rate.

xndarray

Axis-0 coordinate relative to the center.

y_fitndarray

Fitted profile I*kernel + b (plus background offset if bkg_meas given).

y, endarray

Observed normalized intensity profile and its uncertainty.

gaussian(x0, x1, x2, c, inv_S, mode='3d')[source]#

Unnormalized (peak-value-1) Gaussian profile for a given mode.

inv_S is rescaled from the 99.7%-contour convention to 1-sigma, so the result is 1 at the center and exp(-0.5) at Mahalanobis distance 1.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix (99.7%-contour convention).

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

Returns:
gndarray

Unnormalized Gaussian profile on the mode’s coordinate grid.

gaussian_integral(inv_S, mode='3d')[source]#

Normalization integral sqrt((2*pi)**k * det) of the unit-peak Gaussian for a given mode.

k is the mode’s dimensionality and det the determinant of the 1-sigma covariance from inv_S. Multiplying by amplitude A gives the total fitted intensity for that mode.

Parameters:
inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix (99.7%-contour convention).

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

Returns:
integralfloat

Normalization integral of the unit-peak Gaussian.

gaussian_integral_jac_S(inv_S, d_inv_S, mode='3d')[source]#

Derivative of gaussian_integral with respect to a shape parameter group.

Uses Jacobi’s formula d(det(M))/dx = det(M) * trace(M^-1 dM/dx) on the 1-sigma covariance implied by inv_S.

Parameters:
inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix (99.7%-contour convention).

d_inv_Ssequence of ndarray

Three derivatives of inv_S (e.g. from inv_S_deriv_r, inv_S_deriv_size_shape, or inv_S_deriv_domega), each shape (3, 3).

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

Returns:
dgndarray of shape (3,)

Derivative of gaussian_integral w.r.t. each parameter in d_inv_S.

gaussian_jac_S(x0, x1, x2, c, inv_S, d_inv_S, mode='3d')[source]#

Derivative of the unnormalized gaussian profile with respect to a shape parameter group.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix (99.7%-contour convention).

d_inv_Ssequence of ndarray

Three derivatives of inv_S (e.g. from inv_S_deriv_r, inv_S_deriv_size_shape, or inv_S_deriv_domega), each shape (3, 3).

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

Returns:
dgndarray

Stacked derivatives of the Gaussian profile w.r.t. each parameter in d_inv_S, on the mode’s coordinate grid.

gaussian_jac_c(x0, x1, x2, c, inv_S, mode='3d')[source]#

Derivative of the unnormalized gaussian profile with respect to the center c0, c1, c2.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix (99.7%-contour convention).

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

Returns:
dgndarray

Stacked [dg/dc0, dg/dc1, dg/dc2] on the mode’s coordinate grid; entries for axes unused by mode are zero.

integrate(x0, x1, x2, d, n, c, S, b=None, m=None)[source]#

Final post-fit intensity integration combining box-sum, ellipsoid, and profile estimators.

Builds peak/background masks for both the fitted and estimated-fit ellipsoids (the latter giving an edge-clipping volume-fraction correction, applied to the earlier 3D box-sum intensity), then computes and stores the normalization-aware box-sum estimator (extract_intensity, I_ell), the 1D profile-fit estimator (fitted_profile, I_prof), a properly weighted Poisson profile fit (profile_fit, I_fit) seeded from I_ell, and a raw box-sum estimator for diagnostics. Populates self.intensity, self.sigma, self.weights, self.data_norm_fit, self.peak_background_mask, self.integral, self.profile_iterations, and self.info (all intensity estimators and their uncertainties, plus background/volume diagnostics).

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

csequence of float

Fitted peak center (c0, c1, c2).

Sndarray of shape (3, 3)

Fitted ellipsoid covariance-like matrix.

bndarray, optional

Measured background counts, same shape as d.

mndarray, optional

Measured background monitor counts, same shape as d.

Returns:
intensfloat

Final reported intensity: the weighted Poisson profile fit (profile_fit, I_fit), seeded from the box-sum (I_ell).

sigfloat

Uncertainty on intens; inf if not finite.

inv_S_deriv_domega(r0, r1, r2, domega, delta=1e-06)[source]#

Finite-difference derivatives of inv_S with respect to domega_x, domega_y, domega_z.

Central-differences inv_S = R @ diag(1/r0^2, 1/r1^2, 1/r2^2) @ R.T with R(domega) = self._prior_U @ Exp(domega).

Parameters:
r0, r1, r2float

Principal radii along the rotated axes.

domegasequence of float

Relative-rotation vector (domega_x, domega_y, domega_z).

deltafloat, optional

Finite-difference step size. Default 1e-6.

Returns:
dinv_S_x, dinv_S_y, dinv_S_zndarray of shape (3, 3)

d(inv_S)/d(domega_x), d(domega_y), d(domega_z).

inv_S_deriv_r(r0, r1, r2, u0, u1, u2)[source]#

Analytic derivatives of inv_S with respect to each radius r0, r1, r2.

Parameters:
r0, r1, r2float

Principal radii along the rotated axes.

u0, u1, u2float

Orientation rotation vector components.

Returns:
dinv_S0, dinv_S1, dinv_S2ndarray of shape (3, 3)

d(inv_S)/d(r0), d(r1), d(r2).

inv_S_deriv_size_shape(r0, r1, r2, u0, u1, u2)[source]#

Analytic derivatives of inv_S with respect to log_size, shape_1, shape_2.

Chain-rules inv_S_deriv_r through d(r_i)/d(log_size) = r_i and d(r_i)/d(shape_k) = r_i * SHAPE_BASIS[i, k].

Parameters:
r0, r1, r2float

Principal radii along the rotated axes.

u0, u1, u2float

Orientation rotation vector components.

Returns:
dinv_S_size, dinv_S_shape1, dinv_S_shape2ndarray of shape (3, 3)

d(inv_S)/d(log_size), d(shape_1), d(shape_2).

inv_S_matrix(r0, r1, r2, u0, u1, u2)[source]#

Inverse ellipsoid covariance-like matrix inv_S = U diag(1/r0^2, 1/r1^2, 1/r2^2) U^T.

Parameters:
r0, r1, r2float

Principal radii along the rotated axes.

u0, u1, u2float

Orientation rotation vector components.

Returns:
inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

jacobian(params, args_1d, args_2d, args_3d)[source]#

Full stacked Jacobian matching the residual ordering of residual.

Column-stacks the 1D, 2D, 3D Poisson-deviance Jacobians with the prior Jacobian, sanitizes NaN/inf, and transposes to the (n_residuals, n_vary) convention expected by scipy.optimize.least_squares. Passed to lmfit.Minimizer.minimize in sweep.

Parameters:
paramslmfit.Parameters

Current fit parameters.

args_1dsequence

Positional arguments for jacobian_1d: (x0, x1, x2, ds, ns, ws, bkgs).

args_2dsequence

Positional arguments for jacobian_2d (same layout).

args_3dsequence

Positional arguments for jacobian_3d: (x0, x1, x2, d, n, w, bkg).

Returns:
jacndarray of shape (n_residuals, n_vary)

Full Jacobian matrix, transposed to match residual.

jacobian_1d(params, x0, x1, x2, ds, ns, ws=None, bkgs=None, c=None, inv_S=None, dr=None, du=None)[source]#

Column-stacked Poisson-deviance Jacobians across all three 1D projections.

Parameters:
paramslmfit.Parameters

Current fit parameters.

x0, x1, x2ndarray

3D meshgrid coordinates.

dssequence of ndarray

Raw observed counts for modes (“1d_0”, “1d_1”, “1d_2”).

nssequence of ndarray

Normalization counts, one per element of ds.

wssequence, optional

Per-mode weights; coerced via coerce_weights.

bkgssequence of ndarray, optional

Per-mode additional fixed background counts.

csequence of float, optional

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3), optional

Inverse ellipsoid covariance-like matrix.

drtuple of ndarray, optional

Derivatives of inv_S w.r.t. log_size, shape_1, shape_2.

dutuple of ndarray, optional

Derivatives of inv_S w.r.t. domega_x, domega_y, domega_z.

Returns:
jacndarray

Column-stacked jacobian_mode_poisson outputs for “1d_0”, “1d_1”, “1d_2”.

jacobian_2d(params, x0, x1, x2, ds, ns, ws=None, bkgs=None, c=None, inv_S=None, dr=None, du=None)[source]#

Column-stacked Poisson-deviance Jacobians across all three 2D projections.

Parameters:
paramslmfit.Parameters

Current fit parameters.

x0, x1, x2ndarray

3D meshgrid coordinates.

dssequence of ndarray

Raw observed counts for modes (“2d_0”, “2d_1”, “2d_2”).

nssequence of ndarray

Normalization counts, one per element of ds.

wssequence, optional

Per-mode weights; coerced via coerce_weights.

bkgssequence of ndarray, optional

Per-mode additional fixed background counts.

csequence of float, optional

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3), optional

Inverse ellipsoid covariance-like matrix.

drtuple of ndarray, optional

Derivatives of inv_S w.r.t. log_size, shape_1, shape_2.

dutuple of ndarray, optional

Derivatives of inv_S w.r.t. domega_x, domega_y, domega_z.

Returns:
jacndarray

Column-stacked jacobian_mode_poisson outputs for “2d_0”, “2d_1”, “2d_2”.

jacobian_3d(params, x0, x1, x2, d, n, w=None, bkg=None, c=None, inv_S=None, dr=None, du=None)[source]#

Poisson-deviance Jacobian for the full 3D volume.

Parameters:
paramslmfit.Parameters

Current fit parameters.

x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

wNone, scalar, or ndarray, optional

Per-voxel weight; coerced via coerce_weight.

bkgndarray, optional

Additional fixed background counts added to the model.

csequence of float, optional

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3), optional

Inverse ellipsoid covariance-like matrix.

drtuple of ndarray, optional

Derivatives of inv_S w.r.t. log_size, shape_1, shape_2.

dutuple of ndarray, optional

Derivatives of inv_S w.r.t. domega_x, domega_y, domega_z.

Returns:
jacndarray of shape (n_vary, n_valid)

Jacobian restricted to varying parameters and valid voxels.

jacobian_mode_poisson(params, x0, x1, x2, d, n, w, c, inv_S, dr, du, mode, bkg=None)[source]#

Analytic Jacobian of residual_mode_poisson with respect to the varying fit parameters.

Chain-rules the derivative factor from poisson_deviance_residual_factor through mu = n*(A*gaussian + B) [+ bkg], using the Gaussian derivatives w.r.t. center (gaussian_jac_c) and shape (gaussian_jac_S).

Parameters:
paramslmfit.Parameters

Fit parameters; must contain A<mode>, B<mode>, c0, c1, c2, log_size, shape_1, shape_2, domega_x, domega_y, domega_z.

x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts for this mode.

nndarray

Normalization (monitor/solid-angle) counts for this mode.

wNone, scalar, or ndarray

Per-voxel weight; coerced via coerce_weight.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

drtuple of ndarray

Derivatives of inv_S w.r.t. log_size, shape_1, shape_2 (from inv_S_deriv_size_shape).

dutuple of ndarray

Derivatives of inv_S w.r.t. domega_x, domega_y, domega_z (from inv_S_deriv_domega).

modestr

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”.

bkgndarray, optional

Additional fixed background counts added to the model.

Returns:
jacndarray of shape (n_vary, n_valid)

Jacobian restricted to varying parameters and valid voxels.

matched_filter(d, n, kernel, mask, bkg_meas=None)[source]#

Matched-filter Poisson maximum-likelihood fit of an intensity and background level.

Fits mu = n*(I*kernel + b) + bkg_meas to d by maximizing the Poisson log-likelihood over I and b via L-BFGS-B, seeded from percentile-based estimates. Uncertainty on I comes from the observed Fisher information at the optimum.

Parameters:
dndarray

Raw observed event counts, restricted to mask internally.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

kernelndarray

Matched-filter template (e.g. a normalized Gaussian profile), same shape as d.

maskndarray of bool

Voxels to include in the fit.

bkg_measndarray, optional

Additional fixed measured background counts, same shape as d, added to the model.

Returns:
Ifloat

Fitted matched-filter intensity scale.

I_errfloat

Uncertainty on I; inf if the Fisher information is non-positive.

Afloat

Peak amplitude: I scaled by the kernel’s max within mask.

bfloat

Fitted background rate.

Returns (0.0, inf, 0.0, 0.0) if mask selects no voxels.

missing_support_fraction(x0, x1, x2, c, inv_S, n, mode='3d')[source]#

Fraction of the peak’s total Gaussian mass not captured by valid voxels.

The denominator is the analytic whole-space integral (gaussian_integral), not a box-grid sum, so support extending past the fitting box counts as missing too, same as an internal zero-normalization gap. Diagnostic only: doesn’t affect the fit.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

nndarray

Normalization counts on the same grid; non-finite/non-positive entries mark unmeasured voxels.

modestr, optional

Mode accepted by gaussian. Default “3d”.

Returns:
missing_fractionfloat

missing_support / total_support, or nan if total_support <= 0.

missing_supportfloat

total_support - valid_support, clipped at 0.

total_supportfloat

Analytic whole-space Gaussian integral (box-independent).

mode_model_counts(params, x0, x1, x2, d, n, w, c, inv_S, mode, bkg=None)[source]#

Poisson model counts mu = n*(A*gaussian + B) [+ bkg] for one mode.

Parameters:
paramslmfit.Parameters

Fit parameters; must contain A<mode> and B<mode>.

x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Unused (kept for a uniform call signature).

nndarray

Normalization (monitor/solid-angle) counts for this mode.

wndarray

Unused (kept for a uniform call signature).

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

modestr

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”.

bkgndarray, optional

Additional fixed background counts to add to the model.

Returns:
mundarray

Model (expected) counts.

gndarray

Unnormalized Gaussian profile used in the model.

A, Bfloat

Current A<mode> and B<mode> parameter values.

normalize(x0, x1, x2, d, n, w, mode='3d')[source]#

Project raw data onto a mode and normalize it into a weighted intensity/uncertainty.

Combines profile_project and data_norm, dividing the resulting uncertainty by the projected fit weight.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

wndarray

Per-voxel fit weight, same shape as d.

modestr, optional

Projection mode passed to profile_project. Default “3d”.

Returns:
y_intndarray

Normalized intensity for the given mode.

e_intndarray

Weight-scaled uncertainty for the given mode.

peak_roi(x0, x1, x2, c, S, scale, p=0.997)[source]#

Build peak/background/combined masks and a p-confidence Gaussian kernel around an ellipsoid.

The peak mask is the 1-sigma ellipsoid dilated by one voxel; the combined mask grows by further dilations (up to 3 attempts) until the background shell has at least twice as many voxels as the peak. The kernel y is a Gaussian at the p-confidence inverse covariance, independent of the mask-growing logic.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

csequence of float

Peak center (c0, c1, c2).

Sndarray of shape (3, 3)

Ellipsoid covariance-like matrix.

scalefloat

Unused (recomputed locally from p before use).

pfloat, optional

Confidence-contour probability (0-1) for kernel y. Default 0.997.

Returns:
pkndarray of bool

Peak region mask (1-sigma ellipsoid, dilated by one voxel).

bkgndarray of bool

Background shell mask (mask minus pk).

maskndarray of bool

Combined peak-plus-background mask after dilation growth.

yndarray

Gaussian kernel at the p-confidence-contour inverse covariance, normalized to integrate to 1.

poisson_deviance_fit(x0, x1, x2, c, inv_S, y_fit, y, e, mode='3d')[source]#

Reduced Poisson deviance goodness-of-fit statistic within 1-sigma of the peak center.

Restricts to voxels with Mahalanobis distance d2 <= 2**(2/k), approximates raw counts/normalization from (y, e) via the Poisson error model, and returns the summed deviance 2*(mu - d + d*log(d/mu)) divided by degrees of freedom.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

y_fitndarray

Fitted (model) normalized intensity, same shape as y.

yndarray

Observed normalized intensity for this mode.

endarray

Uncertainty on y, same shape as y.

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

Returns:
reduced_deviancefloat

Summed Poisson deviance over degrees of freedom, or inf if degrees of freedom are not positive.

poisson_deviance_residual_factor(d, mu, w=None, eps=1e-12)[source]#

Signed Poisson deviance residual r = sign(mu-d)*sqrt(2*(mu-d+d*log(d/mu))) and dr/dmu.

Used to chain-rule the Jacobian of a residual-based least-squares fit against a Poisson likelihood (jacobian_mode_poisson). Near mu == d, evaluated via the limiting expression to avoid 0/0. For mu < 0 (unphysical mid-optimization), evaluated at a small positive floor and linearly extrapolated to the true mu, giving a large correctly-signed residual with a bounded Jacobian. Weights apply multiplicatively to both r and dr_dmu.

Parameters:
dndarray

Observed (approximate) raw counts.

mundarray

Model (expected) counts, same shape as d.

wndarray, optional

Per-voxel weight. Default 1 everywhere.

epsfloat, optional

Positive floor for clipping mu. Default 1e-12.

Returns:
rndarray

Signed weighted deviance residual, NaN where invalid.

dr_dmundarray

Weighted derivative of r w.r.t. mu, NaN where invalid.

validndarray of bool

Mask where d, mu, w are finite and d >= 0.

prior_jacobian(params)[source]#

Analytic Jacobian of prior_residual with respect to the varying fit parameters.

Column ordering matches prior_residual (3 center + 1 size + 2 shape + 3 orientation); only the center block mixes across parameters, the rest are diagonal.

Parameters:
paramslmfit.Parameters

Fit parameters; must contain c0, c1, c2, log_size, shape_1, shape_2, domega_x, domega_y, domega_z.

Returns:
jacndarray of shape (n_vary, 9)

Jacobian rows restricted to vary=True parameters.

prior_residual(params)[source]#

Whitened regularization residuals pulling the fit toward the resolution-model prior.

Stacks center (raw offset from the origin), size, shape, and orientation residuals, each scaled by its prior sigma (set by update_estimate).

Parameters:
paramslmfit.Parameters

Fit parameters; must contain c0, c1, c2, log_size, shape_1, shape_2, domega_x, domega_y, domega_z.

Returns:
termsndarray of shape (9,)

Concatenated residuals: 3 center + 1 size + 2 shape + 3 orientation.

profile_fit(d, n, kernel, mask, bkg_meas=None, initial_intensity=None, initial_background=None, max_iterations=100)[source]#

Properly weighted Poisson profile fit of intensity and background from a fixed kernel.

Maximizes the same Poisson log-likelihood as matched_filter (mu = n*(I*kernel + background) + bkg_meas) via L-BFGS-B – convex in (I, background), so unlike the reweighted-least- squares approach this converges to the same optimum regardless of the seed – starting from initial_intensity/ initial_background (the summation-integration estimate) when given, falling back to matched_filter’s percentile-based guess otherwise.

Unlike matched_filter, this tracks the optimizer’s iterate trajectory and, at every iterate, the joint (not marginal) Fisher-information uncertainty on I and background – properly accounting for their correlation, which matters most for weak peaks – so the fit’s convergence can be inspected.

kernel (e.g. the un-normalized Gaussian density from peak_roi) is used exactly as given: its absolute scale fixes the units of I, so it must not be rescaled to sum to one, or the recovered peak intensity would be wrong.

Parameters:
dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

kernelndarray

Fixed peak profile, same shape as d, not renormalized.

maskndarray of bool

Voxels to include in the fit (peak plus background region).

bkg_measndarray, optional

Additional fixed measured background counts, same shape as d, added to the model (as in matched_filter).

initial_intensity, initial_backgroundfloat, optional

Starting guess, e.g. from the summation-integration estimate. Falls back to a percentile-based guess when omitted or non-finite.

max_iterationsint

Maximum number of L-BFGS-B iterations.

Returns:
resulttuple

(intensity, sigma_intensity, background, sigma_background, covariance, correlation, snr, iterations, converged). Sigma values are inf if the information matrix is non-positive at the optimum.

historytuple

(iteration, intensity, sigma_intensity, background, sigma_background) arrays tracing the optimizer’s iterates (including the seed as iteration 0), for diagnostic plotting.

profile_project(x0, x1, x2, d, n, w, mode='3d')[source]#

Project raw counts, normalization, variance proxy, and weight onto a given mode.

“1d_i” sums/averages over the two axes other than i; “2d_i” sums/averages over axis i only; “3d” returns full-array copies. Normalization is rate-averaged (rescaled by the integrated-out voxel size) rather than summed. Invalid input voxels are excluded before projecting; projected voxels with invalid resulting normalization are masked to NaN in all outputs.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

wndarray

Per-voxel fit weight, same shape as d.

modestr, optional

One of “1d_0”, “1d_1”, “1d_2”, “2d_0”, “2d_1”, “2d_2”, “3d”. Default “3d”.

Returns:
d_intndarray

Projected (summed) raw counts.

n_intndarray

Projected (rate-averaged) normalization counts.

v_intndarray

Projected (summed) variance proxy from clip(d, 0, None).

w_intndarray

Projected (averaged) fit weight.

Raises:
ValueError

If mode is not a recognized mode string.

project_background(bkg, mode)[source]#

Sum a 3D background-count array onto a given mode’s projection.

Parameters:
bkgndarray

3D background counts array. Non-finite entries treated as zero.

modestr

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”.

Returns:
bkg_projndarray

bkg summed over axes integrated out by mode (a copy for “3d”).

quick_gaussian(x0, x1, x2, d, n, mode='3d', b=None, m=None)[source]#

Quick amplitude/background estimate for a mode, seeding self.params.

Projects the data onto mode, optionally subtracts a rescaled measured background (b, m), estimates background B from the shell beyond Mahalanobis distance 2 and amplitude A from the 95th percentile within distance 1 above B, then adds non-negative-bounded A<mode>/B<mode> parameters.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

modestr, optional

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”. Default “3d”.

bndarray, optional

Measured background counts, same shape as d.

mndarray, optional

Measured background monitor counts, same shape as d, used with n to rescale b into signal-count units.

Returns:
strongerbool or None

True if estimated amplitude exceeds background. None if fewer than 4 valid voxels (no estimate made).

residual(params, args_1d, args_2d, args_3d)[source]#

Full stacked residual vector combining all modes and the prior.

Concatenates the 1D, 2D, 3D Poisson-deviance residuals with the prior residual, sanitizing NaN/inf so scipy.optimize.least_squares sees a well-behaved cost vector. This is the residual function passed to lmfit.Minimizer in sweep.

Parameters:
paramslmfit.Parameters

Current fit parameters.

args_1dsequence

Positional arguments for residual_1d: (x0, x1, x2, ds, ns, ws, bkgs).

args_2dsequence

Positional arguments for residual_2d (same layout).

args_3dsequence

Positional arguments for residual_3d: (x0, x1, x2, d, n, w, bkg).

Returns:
costndarray

Concatenated residual vector (1D + 2D + 3D + prior terms).

residual_1d(params, x0, x1, x2, ds, ns, ws=None, bkgs=None, c=None, inv_S=None)[source]#

Concatenated Poisson-deviance residuals across all three 1D projections.

Parameters:
paramslmfit.Parameters

Current fit parameters.

x0, x1, x2ndarray

3D meshgrid coordinates.

dssequence of ndarray

Raw observed counts for modes (“1d_0”, “1d_1”, “1d_2”).

nssequence of ndarray

Normalization counts, one per element of ds.

wssequence, optional

Per-mode weights; coerced via coerce_weights.

bkgssequence of ndarray, optional

Per-mode additional fixed background counts.

csequence of float, optional

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3), optional

Inverse ellipsoid covariance-like matrix.

Returns:
resndarray

Concatenated residual_mode_poisson outputs for “1d_0”, “1d_1”, “1d_2”.

residual_2d(params, x0, x1, x2, ds, ns, ws=None, bkgs=None, c=None, inv_S=None)[source]#

Concatenated Poisson-deviance residuals across all three 2D projections.

Parameters:
paramslmfit.Parameters

Current fit parameters.

x0, x1, x2ndarray

3D meshgrid coordinates.

dssequence of ndarray

Raw observed counts for modes (“2d_0”, “2d_1”, “2d_2”).

nssequence of ndarray

Normalization counts, one per element of ds.

wssequence, optional

Per-mode weights; coerced via coerce_weights.

bkgssequence of ndarray, optional

Per-mode additional fixed background counts.

csequence of float, optional

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3), optional

Inverse ellipsoid covariance-like matrix.

Returns:
resndarray

Concatenated residual_mode_poisson outputs for “2d_0”, “2d_1”, “2d_2”.

residual_3d(params, x0, x1, x2, d, n, w=None, bkg=None, c=None, inv_S=None)[source]#

Poisson-deviance residual vector for the full 3D volume.

Parameters:
paramslmfit.Parameters

Current fit parameters.

x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

wNone, scalar, or ndarray, optional

Per-voxel weight; coerced via coerce_weight.

bkgndarray, optional

Additional fixed background counts added to the model.

csequence of float, optional

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3), optional

Inverse ellipsoid covariance-like matrix.

Returns:
resndarray

Flattened signed Poisson deviance residuals for valid voxels.

residual_mode_poisson(params, x0, x1, x2, d, n, w, c, inv_S, mode, bkg=None)[source]#

Flattened Poisson-deviance residual vector for one mode, for use in least-squares fitting.

Parameters:
paramslmfit.Parameters

Fit parameters; must contain A<mode> and B<mode>.

x0, x1, x2ndarray

3D meshgrid coordinates.

dndarray

Raw observed event counts for this mode.

nndarray

Normalization (monitor/solid-angle) counts for this mode.

wNone, scalar, or ndarray

Per-voxel weight; coerced via coerce_weight.

csequence of float

Peak center (c0, c1, c2).

inv_Sndarray of shape (3, 3)

Inverse ellipsoid covariance-like matrix.

modestr

One of “3d”, “2d_0”, “2d_1”, “2d_2”, “1d_0”, “1d_1”, “1d_2”.

bkgndarray, optional

Additional fixed background counts added to the model.

Returns:
resndarray

Flattened signed Poisson deviance residuals for valid voxels.

shape_from_params(params)[source]#

Reconstruct absolute radii and orientation rotation vector from the size/shape/orientation parameters.

Converts log_size, shape_1, shape_2 to radii via SHAPE_BASIS and self._prior_r, and domega_* to an absolute rotation vector via R = self._prior_U @ Exp(domega).

Parameters:
paramslmfit.Parameters

Fit parameters; must contain log_size, shape_1, shape_2, domega_x, domega_y, domega_z.

Returns:
r0, r1, r2float

Ellipsoid principal radii along the rotated axes.

u0, u1, u2float

Absolute orientation rotation vector components.

subtract_profile(d, n, perc=30)[source]#

Subtract a per-axis-0-slice background percentile from the raw counts.

For each index along axis 0, estimates a background level as the perc percentile of the normalized intensity over axes 1 and 2, then subtracts it (rescaled by n) from the raw counts.

Parameters:
dndarray

Raw observed event counts, 3D array.

nndarray

Normalization (monitor/solid-angle) counts, same shape as d.

percfloat, optional

Percentile (0-100) used as the background estimate per axis-0 slice. Default 30.

Returns:
d_subndarray

Background-subtracted raw counts, same shape as d.

nndarray

Unmodified normalization counts (for a convenient d, n = subtract_profile(d, n) pattern).

sweep(args_1d, args_2d, args_3d, protocol=None, n_iter=50, report_fit=False)[source]#

Run one lmfit/scipy.optimize.least_squares minimization pass, accepted only if SNR improves.

Optionally sets which parameters vary via protocol, then minimizes self.residual/self.jacobian with lmfit.Minimizer using the "least_squares" method. The resulting fit is accepted (written back to self.params) only if the post-fit 3D I/sigma (_isig_3d at scales=[1]) is not worse than before the fit; otherwise self.params is left unchanged.

Parameters:
args_1dsequence

Positional arguments for residual_1d/jacobian_1d.

args_2dsequence

Positional arguments for residual_2d/jacobian_2d.

args_3dsequence

Positional arguments for residual_3d/jacobian_3d.

protocolsequence of bool, optional

9-element vary flags for c0, c1, c2, log_size, shape_1, shape_2, domega_x, domega_y, domega_z. If None, current vary settings are left unchanged.

n_iterint, optional

Maximum function evaluations for the solver. Default 50.

report_fitbool, optional

Print an lmfit fit report after minimizing. Default False.

uniform_mode_weights(ys, es, val=1.0)[source]#

Build uniform per-mode weight arrays scaled by the total valid point count.

Each array is filled with sqrt(val) / sqrt(n_valid), where n_valid is the total finite, positive-uncertainty point count across all ys/es pairs, giving each mode a total weight contribution proportional to val regardless of point count.

Parameters:
yssequence of ndarray

Per-mode intensity arrays (used only via es, for counting).

essequence of ndarray

Per-mode uncertainty arrays, one per element of ys.

valfloat, optional

Target weight scale for this group of modes. Default 1.0.

Returns:
weightslist of ndarray

One uniform weight array per element of es.

update_constraints(x0, x1, x2, dx)[source]#

(Re)initialize self.params bounds and starting values from the data extent.

Centers start at zero, bounded to half the coordinate extent per axis; size/shape start at zero; orientation offsets start at zero, bounded to [-pi, pi]. Also resets self.combine_params to None.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

dxfloat

Unused (voxel spacing is recomputed from x0, x1, x2).

update_estimate(shape, sigma_c=None)[source]#

Seed self.params and the resolution-model prior from an initial ellipsoid estimate.

Centers reset to the origin; log_size/shape_*/domega_* reset to zero so the fit starts out reproducing this reference ellipsoid. Also (re)initializes the prior widths and sets self.estimated_fit.

The size/shape/rotation prior widths (prior_cov_sigma, prior_distortion_sigma, prior_rot_sigma) are fixed, fractional tolerances around this reference ellipsoid rather than derived per peak, each set by where its parameter’s linear approximation starts to accumulate more than a given relative error:

  • prior_cov_sigma = 0.45: log_size enters through exp(log_size + SHAPE_BASIS @ shape), so radii are linearized as exp(x) ~ 1 + x (relative error x**2/2); 0.45 is the ~10% error point, since the overall size is the least sensitive of the three to its linearization.

  • prior_distortion_sigma = 0.32: same exp(x) ~ 1 + x bound as size, but for shape_1/shape_2, at the tighter ~5% error point.

  • prior_rot_sigma = 0.24 rad: domega enters through the rotation-vector exponential map, linearized as sin(theta) ~ theta (relative error theta**2/6); 0.24 is the ~1% error point, the tightest of the three since orientation is the most sensitive to its linearization.

Parameters:
shapetuple

(c0, c1, c2, r0, r1, r2, v0, v1, v2): (unused) center, principal radii, and orthonormal principal axis vectors (each shape (3,)).

sigma_cfloat, optional

Center prior width (raw Q units), e.g. from UBModel.get_center_uncertainty. Defaults to 1.0 if not given (e.g. when no UB/lattice information is available).

vech6(M)[source]#

Half-vectorize a symmetric 3x3 matrix into its 6 unique entries.

Parameters:
Mndarray of shape (3, 3)

Symmetric matrix to vectorize.

Returns:
vndarray of shape (6,)

The entries [M00, M11, M22, M12, M02, M01].

voxel_volume(x0, x1, x2)[source]#

Volume of a single voxel, the product of the per-axis voxel spacings.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

Returns:
volumefloat

Product of the three per-axis voxel spacings from voxels.

voxels(x0, x1, x2)[source]#

Voxel spacing along each axis, from adjacent-element differences of the meshgrid.

Parameters:
x0, x1, x2ndarray

3D meshgrid coordinates.

Returns:
dx0, dx1, dx2float

Voxel spacing along each axis.

garnet.reduction.integration module#

class Integration(plan)[source]#

Bases: PeakProjection

Attributes:
timer

Timer for the steps of this reduction.

Methods

append_name(file)

Update filename with identifier name.

cell_centering_name()

Lattice and reflection condition.

extract_peak_counts(peaks_ws, r_cut, n_bins)

Obtain raw event counts in an axis-aligned Q-sample box around each peak, binned straight from the events (without converting to MD), without normalization or ellipsoid fitting.

extract_peak_info(peaks_ws, r_cut[, norm, ...])

Obtain peak information for envelope determination.

get_diagnostic_file(name[, ext])

Diagnostic file.

get_diagnostic_path()

Diagnostic path.

get_file(file[, ws])

Update filename with identifier name and optional workspace name.

get_output_file([ext])

Output file.

get_output_path()

Output path.

get_plot_file(name[, ext])

Plot file.

get_plot_path()

Plot path.

integrate_peaks_normalized(data_ws, ...)

Quick integration of every peak: the events in its resolution ellipsoid less those in the equal-volume shell around it, already normalized by flux and solid angle by IntegratePeaksShapeMD.

model_ellipsoids(peaks_ws, res, r_cut)

Each peak's ellipsoid from a resolution model, as containment radii and Q-lab principal axes (columns), within the region; spheres filling the region if there is no model.

modulation_name()

Modulation vectors.

plot_peak_shape_diagnostics(peaks_ws, res, ...)

Per-peak diagnostic figure: raw counts box (from extract_peak_counts) with the observed (whatever ellipsoid shape is currently stored on the peak, e.g. from PeaksModel.integrate_peaks) and model-predicted (res) ellipse cross-sections overlaid on all three 2D projections.

refine_peak_shapes(peaks, data_ws, peaks_ws, ...)

Observed peak shapes for the resolution model, without binning: starting from each peak's prior ellipsoid, the center and shape of every strong enough peak are refined by a profile fit normalized by flux and solid angle, so they are those of the intensity rather than of the counts.

resolution_name()

Minimum d-spacing and starting radii

save_timing()

Save the timings of this process (in the diagnostics directory, which is cleared at the start of each reduction), for write_timing to combine.

write_timing()

Combine the timings saved by every process into <OutputName>_timing.txt, with each step's total over all runs and its mean per run, and remove them.

add_with_padding

bin_extent

check

cleanup

combine

combine_parallel

create_directories

integrate

integrate_parallel

integrate_peaks

optimize_peaks

optimize_ub

pad_to_shape

predict_add_satellite_peaks

predict_all_peaks

project_ellipsoid_parameters

revert_ellipsoid_parameters

transform_Q

unit_key

update_peak_info

update_peak_offsets

validate_params

voxel_weights

write

add_with_padding(a, b, fill=0)[source]#
append_name(file)[source]#

Update filename with identifier name.

Parameters:
filestr

Original file name.

Returns:
output_filestr

File with updated name for identifier name.

cell_centering_name()[source]#

Lattice and reflection condition.

Returns:
lat_refstr

Underscore separated strings.

combine(files)[source]#
static combine_parallel(plan, files)[source]#
extract_peak_counts(peaks_ws, r_cut, n_bins, peak_indices=None, event_ws='data')[source]#

Obtain raw event counts in an axis-aligned Q-sample box around each peak, binned straight from the events (without converting to MD), without normalization or ellipsoid fitting.

Parameters:
peaks_wsstr

Peaks table.

r_cutfloat or ndarray, shape (n_peaks,)

Box half-width along each Q-sample axis. A scalar applies the same half-width to every peak; a per-peak array (e.g. from ResolutionEllipsoid.predict_roi_radii) sizes each peak’s box individually.

n_binsint

Number of bins along each Q-sample axis.

peak_indicesiterable of int, optional

Only bin these peaks. Defaults to all peaks.

event_wsstr, optional

Name of the events, in momentum. The default is “data”.

Returns:
peak_dictdict

Keyed by peak index, each value a dict with “counts” (raw, weighted event counts), “Q0”, “Q1”, “Q2” (bin centers), “Q”, “R”, “lamda”, “two_theta”, “az_phi”.

extract_peak_info(peaks_ws, r_cut, norm=False, fit=True, bank=None)[source]#

Obtain peak information for envelope determination.

Parameters:
peaks_wsstr

Peaks table.

r_cutlist or float

Cutoff radius parameter(s).

get_file(file, ws='')[source]#

Update filename with identifier name and optional workspace name.

Parameters:
filestr

Original file name.

wsstr, optional

Name of workspace. The default is ‘’.

Returns:
output_filestr

File with updated name for identifier and workspace name.

integrate()[source]#
static integrate_parallel(plan, runs, proc)[source]#
integrate_peaks(data)[source]#
integrate_peaks_normalized(data_ws, peaks_ws, res, r_cut)[source]#

Quick integration of every peak: the events in its resolution ellipsoid less those in the equal-volume shell around it, already normalized by flux and solid angle by IntegratePeaksShapeMD.

Parameters:
data_wsstr

Name of the events, cropped for normalization.

peaks_wsstr

Name of peaks table; its peak shapes are left unchanged.

resResolutionEllipsoid

Resolution model giving each peak’s ellipsoid; if it has no model, spheres filling the radius are used instead.

r_cutfloat

Radius of the region around each peak, which bounds the background shell.

model_ellipsoids(peaks_ws, res, r_cut)[source]#

Each peak’s ellipsoid from a resolution model, as containment radii and Q-lab principal axes (columns), within the region; spheres filling the region if there is no model.

Parameters:
peaks_wsstr

Name of peaks table.

resResolutionEllipsoid or None

Resolution model.

r_cutfloat

Radius of the region around each peak, which bounds the background shell.

Returns:
radiilist

Containment radii of each peak.

directionslist

Principal axes (columns) of each peak in Q-lab.

modulation_name()[source]#

Modulation vectors.

Returns:
modstr

Underscore separated vectors and max order

optimize_peaks(data, md, peaks_ws, centering, cell, run, reindex=False, tol=0.15)[source]#
optimize_ub(data, md, peaks_ws, cell, run)[source]#
pad_to_shape(x, shape, fill=0)[source]#
plot_peak_shape_diagnostics(peaks_ws, res, r_cut, n_bins)[source]#

Per-peak diagnostic figure: raw counts box (from extract_peak_counts) with the observed (whatever ellipsoid shape is currently stored on the peak, e.g. from PeaksModel.integrate_peaks) and model-predicted (res) ellipse cross-sections overlaid on all three 2D projections.

When self.prior_res is set (an instrument prior model was used to size the adaptive per-peak integration radii, see integrate), a second row is prepended comparing that prior model’s prediction against the peak’s stored shape – i.e. two rows per peak: prior vs. observed, then observed vs. the now-updated population model res.

One PNG per peak, saved next to that peak’s other diagnostic files (extract_peak_counts’s “…_counts.nxs”, same peak-named subdirectory). Complements res.plot_diagnostics’s population-level obs-vs-pred scatter with a per-peak, spatial check of the same thing – whether the peak’s stored shape actually matches the counts around it, and whether the box (r_cut) used for the box display looks well-sized.

Parameters:
peaks_wsstr

Peaks table.

resResolutionEllipsoid

Fitted resolution model (fit() already called).

r_cutfloat or ndarray, shape (n_peaks,)

Box half-width(s) for the counts box shown behind the overlay – purely for display, independent of whatever produced the peak’s stored shape.

n_binsint

Number of bins along each Q-sample axis.

predict_add_satellite_peaks(peaks_ws, md_ws, lamda_min, lamda_max)[source]#
predict_all_peaks(data, md, peaks_ws, centering, d_min, lamda_min, lamda_max)[source]#
refine_peak_shapes(peaks, data_ws, peaks_ws, prior_res, r_cut)[source]#

Observed peak shapes for the resolution model, without binning: starting from each peak’s prior ellipsoid, the center and shape of every strong enough peak are refined by a profile fit normalized by flux and solid angle, so they are those of the intensity rather than of the counts. A refined peak is moved to its fitted center and given its fitted ellipsoid (as containment radii, in the sample frame); the others are left without a shape, so the resolution model ignores them.

Parameters:
peaksPeaksModel

Peaks model to integrate with.

data_wsstr

Name of the events, cropped for normalization.

peaks_wsstr

Name of peaks table.

prior_resResolutionEllipsoid or None

Prior resolution model, from the instrument configuration.

r_cutfloat

Radius of the region around each peak.

resolution_name()[source]#

Minimum d-spacing and starting radii

Returns:
res_radstr

Underscore separated strings.

unit_key(v, tol=0.01)[source]#
update_peak_info(peaks_ws, peak_dict)[source]#
update_peak_offsets(peaks_ws, offsets, peak_dict)[source]#
validate_params()[source]#
write(result_file)[source]#

garnet.reduction.normalization module#

class Normalization(plan)[source]#

Bases: SubPlan

Attributes:
timer

Timer for the steps of this reduction.

Methods

append_name(file)

Update filename with identifier name.

binning_name()

Bin size for each dimension.

combine(files)

Merge data and normalization files.

elastic_name()

Elastic channel.

extents_name()

Min/max pairs for each dimensional extents.

get_diagnostic_file(name[, ext])

Diagnostic file.

get_diagnostic_path()

Diagnostic path.

get_file(file[, ws])

Update filename with identifier name and optional workspace name.

get_output_file([ext])

Output file.

get_output_path()

Output path.

get_plot_file(name[, ext])

Plot file.

get_plot_path()

Plot path.

projection_name()

Axes projections.

save_timing()

Save the timings of this process (in the diagnostics directory, which is cleared at the start of each reduction), for write_timing to combine.

symmetry_name()

Laue group name.

write_timing()

Combine the timings saved by every process into <OutputName>_timing.txt, with each step's total over all runs and its mean per run, and remove them.

check

cleanup

combine_parallel

create_directories

normalize

normalize_parallel

validate_params

view

append_name(file)[source]#

Update filename with identifier name.

Parameters:
filestr

Original file name.

Returns:
output_filestr

File with updated name for identifier name.

binning_name()[source]#

Bin size for each dimension.

_N0xN1xN2

Returns:
binsstr

Cross separated integers.

combine(files)[source]#

Merge data and normalization files.

Parameters:
fileslist

Files to be combined.

static combine_parallel(plan, files)[source]#
elastic_name()[source]#

Elastic channel.

Returns:
ccstr

Total or elastic channel.

extents_name()[source]#

Min/max pairs for each dimensional extents.

_[min_0,max_0]_[min_1,max_1]_[min_2,max_2]

Returns:
extentsstr

Underscore separated list.

get_file(file, ws='')[source]#

Update filename with identifier name and optional workspace name.

Parameters:
filestr

Original file name.

wsstr, optional

Name of workspace. The default is ‘’.

Returns:
output_filestr

File with updated name for identifier and workspace name.

normalize()[source]#
static normalize_parallel(plan, runs, proc)[source]#
projection_name()[source]#

Axes projections.

Returns:
projstr

Name of slices.

symmetry_name()[source]#

Laue group name. Spaces are removed and slashes are replaced with underscore.

Returns:
symmetrystr

None or Hermann-Mauguin point group symbol.

validate_params()[source]#
view(result_file)[source]#

garnet.reduction.parallel module#

class ParallelProcessor(n_proc=1)[source]#

Bases: object

Methods

process_dict

process_dict(data, func)[source]#
class ParallelTasks(function, combine=None)[source]#

Bases: object

Methods

run_tasks(plan, n_proc)

Run parallel tasks with processing pool.

safe_function_wrapper

run_tasks(plan, n_proc)[source]#

Run parallel tasks with processing pool.

Parameters:
plandict

Data reduction plan split over each process.

n_procint

Number of processes.

safe_function_wrapper(*args, **kwargs)[source]#
dethread_environment()[source]#
reset_environment()[source]#

garnet.reduction.parametrization module#

class Parametrization(plan)[source]#

Bases: SubPlan

Attributes:
timer

Timer for the steps of this reduction.

Methods

append_name(file)

Update filename with identifier name.

binning_name()

Bin size for each dimension.

combine(files)

Merge data and normalization files.

elastic_name()

Elastic channel.

extents_name()

Min/max pairs for each dimensional extents.

get_diagnostic_file(name[, ext])

Diagnostic file.

get_diagnostic_path()

Diagnostic path.

get_file(file[, ws])

Update filename with identifier name and optional workspace name.

get_output_file([ext])

Output file.

get_output_path()

Output path.

get_plot_file(name[, ext])

Plot file.

get_plot_path()

Plot path.

log_name()

Log name.

projection_name()

Axes projections.

save_timing()

Save the timings of this process (in the diagnostics directory, which is cleared at the start of each reduction), for write_timing to combine.

write_timing()

Combine the timings saved by every process into <OutputName>_timing.txt, with each step's total over all runs and its mean per run, and remove them.

check

cleanup

combine_parallel

create_directories

parametrize

parametrize_parallel

validate_params

view

append_name(file)[source]#

Update filename with identifier name.

Parameters:
filestr

Original file name.

Returns:
output_filestr

File with updated name for identifier name.

binning_name()[source]#

Bin size for each dimension.

_N0xN1xN2

Returns:
binsstr

Cross separated integers.

combine(files)[source]#

Merge data and normalization files.

Parameters:
fileslist

Files to be combined.

static combine_parallel(plan, files)[source]#
elastic_name()[source]#

Elastic channel.

Returns:
ccstr

Total or elastic channel.

extents_name()[source]#

Min/max pairs for each dimensional extents.

_[min_0,max_0]_[min_1,max_1]_[min_2,max_2]

Returns:
extentsstr

Underscore separated list.

get_file(file, ws='')[source]#

Update filename with identifier name and optional workspace name.

Parameters:
filestr

Original file name.

wsstr, optional

Name of workspace. The default is ‘’.

Returns:
output_filestr

File with updated name for identifier and workspace name.

log_name()[source]#

Log name. Spaces are removed and slashes are replaced with underscore.

Returns:
logstr

None or Hermann-Mauguin point group symbol.

parametrize()[source]#
static parametrize_parallel(plan, runs, proc)[source]#
projection_name()[source]#

Axes projections.

Returns:
projstr

Name of slices.

validate_params()[source]#
view(result_file)[source]#

garnet.reduction.peaks module#

class PeakModel(peaks)[source]#

Bases: object

Methods

add_diagnostic_info(no, info)

Log diagnostic info.

clear_peak_shape(no)

Remove the shape of the peak.

get_UB()

UB matrix.

get_angles(no)

Scattering and azimuthal angle of the peak.

get_bank_name(no)

Obtain the bank name.

get_d_from_ub(no)

D-spacing from current oriented lattice for the given peak.

get_d_spacing(no)

Obtain the peak d-spacing.

get_detector_id(no)

Obtain the peak detector id number.

get_goniometer_angles(no)

Goniometer Euler angles of the peak.

get_goniometer_matrix(no)

Goniometer matrix of the peak.

get_hkl(no)

Miller indices.

get_hklmnp(no)

Miller indices.

get_number_peaks()

Total number of peaks in the table.

get_peak_intensity(no)

Obtain the peak intensity.

get_peak_name(no[, merge, d])

Name of peak.

get_peak_shape(no[, r_cut])

Obtain the peak shape parameters.

get_sample_Q(no)

Scattering vector in Q sample coordinates.

get_signal_to_noise(no)

Intensity to uncertainty ratio.

get_wavelength(no)

Wavelength of the peak.

set_hklmnp(no, hklmnp)

Update Miller indices.

set_peak_center(no, c0, c1, c2)

Update the shape of the peak.

set_peak_intensity(no, intens, sig)

Update the peak intensity.

set_peak_shape(no, c0, c1, c2, r0, r1, r2, ...)

Update the shape of the peak.

set_scale_factor(no, scale)

update the peak normalization scale factor.

set_wavelength(no, lamda)

Update the wavelength of the peak.

get_projection_matrix

get_projection_peak_origin

add_diagnostic_info(no, info)[source]#

Log diagnostic info.

Parameters:
noint

Peak index number.

infodict

Mapping of diagnostic key → scalar value. Keys must be a subset of diagnostic_keys (positional run/hkl keys are filled from the peak itself and must not appear in info).

clear_peak_shape(no)[source]#

Remove the shape of the peak.

Parameters:
noint

Peak index number.

get_UB()[source]#

UB matrix.

Returns:
UB2d-array, 3x3

UB-matrix

get_angles(no)[source]#

Scattering and azimuthal angle of the peak.

Parameters:
noint

Peak index number.

Returns:
two_thetafloat

Scattering (polar) angle in degrees.

az_phifloat

Azimuthal angle in degrees.

get_bank_name(no)[source]#

Obtain the bank name.

Parameters:
noint

Peak index number.

Returns:
bankstr

Bank name.

get_d_from_ub(no)[source]#

D-spacing from current oriented lattice for the given peak.

Parameters:
noint

Peak index number.

Returns:
dfloat

D-spacing in Angstroms.

get_d_spacing(no)[source]#

Obtain the peak d-spacing.

Parameters:
noint

Peak index number.

Returns:
dfloat

Interplanar spacing.

get_detector_id(no)[source]#

Obtain the peak detector id number.

Parameters:
noint

Peak index number.

Returns:
det_idint

Detector number.

get_goniometer_angles(no)[source]#

Goniometer Euler angles of the peak.

Parameters:
noint

Peak index number.

Returns:
angleslist

Euler angles (YZY convention) in degrees.

get_goniometer_matrix(no)[source]#

Goniometer matrix of the peak.

Parameters:
noint

Peak index number.

Returns:
R2d-array

Rotation matrix.

get_hkl(no)[source]#

Miller indices.

Parameters:
noint

Peak index number.

Returns:
hkllist

Components of Miller indices.

get_hklmnp(no)[source]#

Miller indices.

Parameters:
noint

Peak index number.

Returns:
hklmnptuple

Components of Miller indices.

get_number_peaks()[source]#

Total number of peaks in the table.

Returns:
nint

Number of peaks

get_peak_intensity(no)[source]#

Obtain the peak intensity.

Parameters:
noint

Peak index number.

Returns:
intensfloat

Intensity.

sigfloat

Uncertainty.

get_peak_name(no, merge=False, d=None)[source]#

Name of peak.

Parameters:
noint

Peak index number.

dfloat, optional

D-spacing to use in the name. When provided the lattice lookup is skipped, so the name is stable even after UB refinement.

Returns:
namestr

Readable name of peak.

get_peak_shape(no, r_cut=inf)[source]#

Obtain the peak shape parameters.

Parameters:
noint

Peak index number.

Returns:
c0, c1, c2float

Peak center.

r0, r1, r2float

Principal radii.

v0, v1, v2list

Principal axis directions.

get_projection_matrix(no)[source]#
get_projection_peak_origin(no)[source]#
get_sample_Q(no)[source]#

Scattering vector in Q sample coordinates.

Parameters:
noint

Peak index number.

Returns:
Qlist

Scattering Q-vector.

get_signal_to_noise(no)[source]#

Intensity to uncertainty ratio.

Parameters:
noint

Peak index number.

Returns:
I_sigfloat

Signal-noise ratio.

get_wavelength(no)[source]#

Wavelength of the peak.

Parameters:
noint

Peak index number.

Returns:
lamdafloat

Wavelength in angstroms.

set_hklmnp(no, hklmnp)[source]#

Update Miller indices.

Parameters:
noint

Peak index number.

hklmnplist

Components of Miller indices.

set_peak_center(no, c0, c1, c2)[source]#

Update the shape of the peak.

Parameters:
noint

Peak index number.

c0, c1, c2float

Peak center.

set_peak_intensity(no, intens, sig)[source]#

Update the peak intensity.

Parameters:
noint

Peak index number.

intensfloat

Intensity.

sigfloat

Uncertainty.

set_peak_shape(no, c0, c1, c2, r0, r1, r2, v0, v1, v2)[source]#

Update the shape of the peak.

Parameters:
noint

Peak index number.

c0, c1, c2float

Peak center.

r0, r1, r2float

Principal radii.

v0, v1, v2list

Principal axis directions.

set_scale_factor(no, scale)[source]#

update the peak normalization scale factor.

Parameters:
noint

Peak index number.

scalefloat

Peak scale factor

set_wavelength(no, lamda)[source]#

Update the wavelength of the peak.

Parameters:
noint

Peak index number.

lamdafloat

Wavelength in angstroms.

class PeaksModel[source]#

Bases: object

Methods

add_peak(peaks, hkl)

Add a peak to an existing table.

centroid_peaks(md, peaks, peak_radius)

Re-center peak locations using centroid within given radius

combine_peaks(peaks, merge)

Merge two peaks workspaces into one.

convert_peaks(peaks)

Remove instrument from peaks.

copy_shaped_peaks(peaks, shaped)

Copy only the peaks that have an ellipsoidal shape.

create_peaks(ws, peaks[, lean])

Create a new peaks table.

delete_peaks(peaks)

Remove peaks.

filter_unique_hkl(peaks, unique_peaks)

Keep one representative per unique (h,k,l,m,n,p) group.

find_peaks(md, peaks, min_Q[, density, ...])

Harvest strong peak locations from Q-sample into a peaks table.

get_UB(ws)

Obtain UB from the oriented lattice.

get_all_goniometer_matrices(ws)

Extract all goniometer matrices.

get_bank_names(peaks)

Obtain the bank names.

get_max_d_spacing(ws)

Obtain the maximum d-spacing from the oriented lattice.

get_peaks_name(peaks)

Name of peaks.

integrate_ellipsoids(data, peaks, peak_radius)

Integrate peaks using ellipsoidal regions.

integrate_peaks(md, peaks, peak_radius[, ...])

Integrate peaks using spherical or ellipsoidal regions.

integrate_peaks_normalized(data, peaks, ...)

Integrate peaks in ellipsoids, normalized by flux and solid angle along each detector's trajectory through them (as MDNorm, but without binning).

integrate_peaks_with_radii(md, peaks, peak_radii)

Integrate peaks using sp ellipsoidal regions.

intensity_profile(md, peaks, peak_radius[, ...])

Integrate peak intensity as profile.

intensity_projection(md, peaks, peak_radius)

Integrate peak intensity as profile.

intensity_vs_radius(md, peaks, peak_radius)

Integrate peak intensity with radius varying from zero to cut off.

load_peaks(filename, peaks)

Load peaks file.

predict_modulated_peaks(peaks, centering, ...)

Predict the modulated peak positions based on main peaks.

predict_peaks(ws, peaks, centering, d_min, ...)

Predict peak Q-sample locations with UB and lattice centering.

predict_satellite_peaks(peaks_ws, data_ws, ...)

Locate satellite peaks from goniometer angles.

remove_duplicate_peaks(peaks)

Omit duplicate peaks from different based on indexing.

remove_peaks_by_d_tolerance(peaks[, tol])

Filter out peaks based on d-spacing tolerance.

remove_unindexed_peaks(peaks)

Filter out unindexes peaks.

remove_weak_peaks(peaks[, sig_noise])

Filter out weak peaks based on signal-to-noise ratio.

renumber_by_size(peaks[, n])

Renumber peak RunNumbers 1..n ordered by observed peak shape size.

renumber_runs_by_index(ws, peaks)

Re-label the runs by index based on goniometer setting.

restore_run_number(peaks)

Restore the run number stashed in bin count.

save_hkl_cw(filename, peaks)

Save peaks in SHELXL CW format sorted by d-spacing.

save_peaks(filename, peaks)

Save peaks file.

scan_threshold(md, peaks, min_Q[, min_found])

Scan peak density threshold from predicted peaks and minimum Q.

set_goniometer(peaks, R)

Update the goniometer on the run.

sort_peaks_by_bank(peaks)

Sort peaks table by ascending bank number

sort_peaks_by_d(peaks)

Sort peaks table by descending d-spacing.

sort_peaks_by_hkl(peaks)

Sort peaks table by descending hkl values.

stash_run_number(peaks)

Park the run number in bin count ahead of a temporary renumbering.

update_scale_factor(peaks, value)

Update counting statistic refrence value for normalization

extract_peaks_roi

get_number_peaks

index_peaks

remove_aluminum_contamination

reset_satellites

add_peak(peaks, hkl)[source]#

Add a peak to an existing table.

Parameters:
peaksstr

Name of peaks table.

hkllist

Miller index.

centroid_peaks(md, peaks, peak_radius)[source]#

Re-center peak locations using centroid within given radius

Parameters:
mdstr

Name of Q-sample.

peaksstr

Name of peaks table.

peak_radiusfloat

Integration region radius.

combine_peaks(peaks, merge)[source]#

Merge two peaks workspaces into one.

Parameters:
peaksstr

Name of peaks table to be added.

mergestr

Name of peaks table to be accumulated.

convert_peaks(peaks)[source]#

Remove instrument from peaks.

Parameters:
peaksstr

Name of peaks table.

copy_shaped_peaks(peaks, shaped)[source]#

Copy only the peaks that have an ellipsoidal shape.

Parameters:
peaksstr

Name of peaks table.

shapedstr

Name of the copy.

Returns:
nint

Number of peaks copied.

create_peaks(ws, peaks, lean=False)[source]#

Create a new peaks table.

wsstr

Name of workspace.

peaksstr

Name of peaks table.

delete_peaks(peaks)[source]#

Remove peaks.

Parameters:
peaksstr

Name of peaks table to be added.

extract_peaks_roi(md, peaks, r_cut, n_bins=21)[source]#
filter_unique_hkl(peaks, unique_peaks)[source]#

Keep one representative per unique (h,k,l,m,n,p) group.

The representative is the peak whose wavelength is closest to the median wavelength of the group (middle-wavelength observation).

Parameters:
peaksstr

Name of input peaks table.

unique_peaksstr

Name of output peaks table with one peak per unique Miller index.

find_peaks(md, peaks, min_Q, density=50, max_peaks=1000)[source]#

Harvest strong peak locations from Q-sample into a peaks table.

Parameters:
mdstr

Name of Q-sample.

peaksstr

Name of peaks table.

min_Qfloat

Minimum Q-spacing enforcing lower limit of peak spacing.

densityint, optional

Threshold density. The default is 1000.

max_peaksint, optional

Maximum number of peaks to find. The default is 50.

get_UB(ws)[source]#

Obtain UB from the oriented lattice.

Parameters:
wsstr

Workspace with UB defined on oriented lattice.

Returns:
UB2d-array

UB matrix.

get_all_goniometer_matrices(ws)[source]#

Extract all goniometer matrices.

Parameters:
wsstr

Name of workspace with goniometer indexing.

Returns:
Rs: list

Goniometer matrices.

get_bank_names(peaks)[source]#

Obtain the bank names.

Parameters:
peaksstr

Name of peaks table.

Returns:
banksstr

Unique bank names.

get_max_d_spacing(ws)[source]#

Obtain the maximum d-spacing from the oriented lattice.

Parameters:
wsstr

Workspace with UB defined on oriented lattice.

Returns:
d_maxfloat

Maximum d-spacing.

get_number_peaks()[source]#
get_peaks_name(peaks)[source]#

Name of peaks.

Returns:
namestr

Readable name of peaks.

index_peaks(peaks, tol=0.15)[source]#
integrate_ellipsoids(data, peaks, peak_radius, cutoff=10)[source]#

Integrate peaks using ellipsoidal regions. Ellipsoid integration adapts itself to the peak distribution.

Parameters:
datastr

Name of event workspace.

peaksstr

Name of peaks table.

peak_radiusfloat

Integration region radius.

cutofffloat

Weak peak threshold. Default is 10.

integrate_peaks(md, peaks, peak_radius, radius_scale=0, background_inner_fact=np.float64(1.2599210498948732), background_outer_fact=np.float64(1.4422495703074083), method='ellipsoid', centroid=True, update=True, integrate_if_on_edge=True, mask_edge_tubes=False)[source]#

Integrate peaks using spherical or ellipsoidal regions. Ellipsoid integration adapts itself to the peak distribution.

Parameters:
mdstr

Name of Q-sample.

peaksstr

Name of peaks table.

peak_radiusfloat

Integration region radius.

radius_scalefloat, optional

Radius scale factor with |Q|. The default is 0.

background_inner_factfloat, optional

Factor of peak radius for background shell. The default is cbrt(2).

background_outer_factfloat, optional

Factor of peak radius for background shell. The default is cbrt(3).

methodstr, optional

Integration method. The default is ‘sphere’.

centroidbool, optional

Shift peak position to centroid. The default is True.

updatebool, optional

Update peak position. The defualt is True.

integrate_if_on_edgebool, optional

Integrate (with a truncated, biased shape) a peak whose region extends past a detector panel/tube boundary instead of excluding it. The default is True.

mask_edge_tubesbool, optional

Exclude peaks landing on a detector’s edge tubes/columns entirely. The default is False.

integrate_peaks_normalized(data, peaks, peak_radii, peak_directions, region_radius, background_inner_fact=np.float64(1.2599210498948732), background_outer_fact=np.float64(1.4422495703074083), solid_angle='sa', flux='flux', profile_fit=False, adjust_center=False, adjust_shape=False, min_signal_to_noise=10.0, max_shape_change=2.0, max_center_shift=3.0)[source]#

Integrate peaks in ellipsoids, normalized by flux and solid angle along each detector’s trajectory through them (as MDNorm, but without binning). The peaks’ own shapes are left unchanged.

By default the events inside each ellipsoid less those in its background shell are counted. With profile_fit, a Gaussian whose standard deviations are the radii is fitted instead, and its center and shape may be refined too.

Parameters:
datastr

Name of events, cropped for normalization.

peaksstr

Name of peaks table.

peak_radiilist

Principal radii of each peak’s ellipsoid.

peak_directionslist

Principal axes (columns) of each peak’s ellipsoid in Q-lab.

region_radiusfloat

Radius of the events considered around each peak.

background_inner_factfloat, optional

Factor of peak radii for background shell. The default is cbrt(2).

background_outer_factfloat, optional

Factor of peak radii for background shell. The default is cbrt(3).

solid_anglestr, optional

Name of solid angle workspace. The default is “sa”.

fluxstr, optional

Name of integrated flux workspace. The default is “flux”.

profile_fitbool, optional

Fit a Gaussian instead of counting. The default is False.

adjust_centerbool, optional

With profile_fit, refine each peak’s center. The default is False.

adjust_shapebool, optional

With profile_fit, refine the shape of each peak with at least min_signal_to_noise. The default is False.

min_signal_to_noisefloat, optional

Smallest intensity over sigma for which a shape is refined. The default is 10.

max_shape_changefloat, optional

Largest factor by which a refined standard deviation may change. The default is 2.

max_center_shiftfloat, optional

With adjust_shape and adjust_center, the largest distance a refined center may move, in standard deviations of the peak’s ellipsoid along each of its axes. The default is 3.

Returns:
diagnosticslist or None

For each peak, counts and normalization in the peak (or fit region) and the shell, whether the peak is only partly covered by the detectors or momentum band, and whether its shape was refined, with the refined standard deviations, Q-lab axes (columns) and center offset, and the time spent integrating it; None if integration failed.

integrate_peaks_with_radii(md, peaks, peak_radii, background_inner_fact=np.float64(1.2599210498948732), background_outer_fact=np.float64(1.4422495703074083), centroid=True, update=True)[source]#

Integrate peaks using sp ellipsoidal regions. Ellipsoid integration adapts itself to the peak distribution.

Parameters:
mdstr

Name of Q-sample.

peaksstr

Name of peaks table.

peak_radiilist

Integration region radii.

background_inner_factfloat, optional

Factor of peak radius for background shell. The default is cbrt(2).

background_outer_factfloat, optional

Factor of peak radius for background shell. The default is cbrt(3).

centroidbool, optional

Shift peak position to centroid. The default is True.

updatebool, optional

Update peak position. The defualt is True.

intensity_profile(md, peaks, peak_radius, background_inner_fact=1, background_outer_fact=1.5)[source]#

Integrate peak intensity as profile.

Parameters:
mdstr

Name of Q-sample.

peaksstr

Name of peaks table.

peak_radiusfloat

Integrat region radius cut off.

background_inner_factfloat, optional

Factor of peak radius for background shell. The default is 1.

background_outer_factfloat, optional

Factor of peak radius for background shell. The default is 1.5.

Returns:
xlist

Peak profile.

ylist

Peak intensity.

lamdalist

Peak wavelength.

intensity_projection(md, peaks, peak_radius, background_inner_fact=1, background_outer_fact=1.5)[source]#

Integrate peak intensity as profile.

Parameters:
mdstr

Name of Q-sample.

peaksstr

Name of peaks table.

peak_radiusfloat

Integrat region radius cut off.

background_inner_factfloat, optional

Factor of peak radius for background shell. The default is 1.

background_outer_factfloat, optional

Factor of peak radius for background shell. The default is 1.5.

Returns:
xlist

Peak projection.

ylist

Peak intensity.

thetalist

Peak angle.

intensity_vs_radius(md, peaks, peak_radius, background_inner_fact=1, background_outer_fact=1.5, steps=51, fix=False)[source]#

Integrate peak intensity with radius varying from zero to cut off.

Parameters:
mdstr

Name of Q-sample.

peaksstr

Name of peaks table.

peak_radiusfloat

Integration region radius cut off.

background_inner_factfloat, optional

Factor of peak radius for background shell. The default is 1

background_outer_factfloat, optional

Factor of peak radius for background shell. The default is 1.5.

stepsint, optional

Number of integration steps. The default is 101.

fixbool, optional

Fix the background shell size

Returns:
radiuslist

Peak radius.

sig_noiselist

Peak signal/noise ratio at lowest threshold.

intenslist

Peak intensity.

load_peaks(filename, peaks)[source]#

Load peaks file.

Parameters:
filenamestr

Name of peaks file with extension .nxs.

peaksstr

Name of peaks table.

predict_modulated_peaks(peaks, centering, d_min, lamda_min, lamda_max, mod_vec_1=[0, 0, 0], mod_vec_2=[0, 0, 0], mod_vec_3=[0, 0, 0], max_order=0, cross_terms=False)[source]#

Predict the modulated peak positions based on main peaks.

Parameters:
wsstr

Name of workspace to predict peaks with UB.

peaksstr

Name of main peaks table.

centeringstr

Lattice centering that provides the reflection condition.

d_minfloat

The lower d-spacing resolution to predict peaks.

lamda_min, lamda_maxfloat

The wavelength band over which to predict peaks.

mod_vec_1, mod_vec_2, mod_vec_3list, optional

Modulation vectors. The default is [0,0,0].

max_orderint, optional

Maximum order greater than zero for satellites. The default is 0.

cross_termsbool, optional

Include modulation cross terms. The default is False.

predict_peaks(ws, peaks, centering, d_min, lamda_min, lamda_max)[source]#

Predict peak Q-sample locations with UB and lattice centering.

Symbol

Reflection condition

P

None

I

\(h+k+l=2n\)

F

\(h,k,l\) unmixed

R

\(-h+k+l=3n\) or \(-h+k+l=3n\)

R(obv)

\(-h+k+l=3n\)

R(rev)

\(h-k+l=3n\)

A

\(k+l=2n\)

B

\(l+h=2n\)

C

\(h+k=2n\)

Parameters:
wsstr

Name of workspace to predict peaks with UB.

peaksstr

Name of peaks table.

centeringstr

Lattice centering that provides the reflection condition.

d_minfloat

The lower d-spacing resolution to predict peaks.

lamda_min, lamda_maxfloat

The wavelength band over which to predict peaks.

predict_satellite_peaks(peaks_ws, data_ws, centering, lamda_min, lamda_max, d_min, mod_vec_1=[0, 0, 0], mod_vec_2=[0, 0, 0], mod_vec_3=[0, 0, 0], max_order=0, cross_terms=False)[source]#

Locate satellite peaks from goniometer angles.

Parameters:
peaks_wsstr

Reference peaks table.

data_wsstr

Q-sample data with goniometer(s).

centeringstr

Lattice centering that provides the reflection condition.

lamda_minfloat

Minimum wavelength.

lamda_maxfloat

Maximum wavelength.

d_minfloat

The lower d-spacing resolution to predict peaks.

mod_vec_1, mod_vec_2, mod_vec_3list, optional

Modulation vectors. The default is [0,0,0].

max_orderint, optional

Maximum order greater than zero for satellites. The default is 0.

cross_termsbool, optional

Include modulation cross terms. The default is False.

remove_aluminum_contamination(peaks, d_min, d_max, delta=0.1)[source]#
remove_duplicate_peaks(peaks)[source]#

Omit duplicate peaks from different based on indexing. Table will be sorted.

Parameters:
peaksstr

Name of peaks table.

remove_peaks_by_d_tolerance(peaks, tol=0.05)[source]#

Filter out peaks based on d-spacing tolerance.

Parameters:
peaksstr

Name of peaks table.

tolfloat, optional

d-spacing tolerance. The default is 0.05.

remove_unindexed_peaks(peaks)[source]#

Filter out unindexes peaks.

Parameters:
peaksstr

Name of peaks table.

remove_weak_peaks(peaks, sig_noise=3)[source]#

Filter out weak peaks based on signal-to-noise ratio.

Parameters:
peaksstr

Name of peaks table.

sig_noisefloat, optional

Minimum signal-to-noise ratio. The default is 3.

renumber_by_size(peaks, n=None)[source]#

Renumber peak RunNumbers 1..n ordered by observed peak shape size.

Bins peaks by their currently-stored peak shape’s largest radius (set by a prior integrate_peaks/integrate_ellipsoids call) into n equal-count bins, smallest (RunNumber=1) to largest (RunNumber=n) – for use with integrate_peaks_with_radii, which selects each radius’s peak group by RunNumber. Call stash_run_number first to preserve the real run numbers, and restore_run_number after re-integrating.

Parameters:
peaksstr

Name of peaks table.

nint or None

Number of size bins. Defaults to one bin per peak.

Returns:
radii_binsndarray, shape (n_bins,)

Mean largest radius of the peaks assigned to each bin, ordered from smallest (bin 1) to largest (bin n).

renumber_runs_by_index(ws, peaks)[source]#

Re-label the runs by index based on goniometer setting.

Parameters:
wsstr

Name of workspace with goniometer indexing.

peaksstr

Name of peaks table.

reset_satellites(peaks)[source]#
restore_run_number(peaks)[source]#

Restore the run number stashed in bin count.

Parameters:
peaksstr

Name of peaks table.

save_hkl_cw(filename, peaks)[source]#

Save peaks in SHELXL CW format sorted by d-spacing.

Parameters:
filenamestr

Name of output file with extension .hkl.

peaksstr

Name of peaks table.

save_peaks(filename, peaks)[source]#

Save peaks file.

Parameters:
filenamestr

Name of peaks file with extension .nxs.

peaksstr

Name of peaks table.

scan_threshold(md, peaks, min_Q, min_found=50)[source]#

Scan peak density threshold from predicted peaks and minimum Q.

Parameters:
mdstr

Name of Q-sample.

peaksstr

Name of predicted peaks table.

min_Qfloat

Minimum Q-spacing enforcing lower limit of peak spacing.

min_foundint, optional

Minimum number of found peaks desired for UB determination. The default is 50.

set_goniometer(peaks, R)[source]#

Update the goniometer on the run.

Parameters:
peaksstr

Name of peaks table.

R2d-array

Goniometer matrix.

sort_peaks_by_bank(peaks)[source]#

Sort peaks table by ascending bank number

Parameters:
peaksstr

Name of peaks table.

sort_peaks_by_d(peaks)[source]#

Sort peaks table by descending d-spacing.

Parameters:
peaksstr

Name of peaks table.

sort_peaks_by_hkl(peaks)[source]#

Sort peaks table by descending hkl values.

Parameters:
peaksstr

Name of peaks table.

stash_run_number(peaks)[source]#

Park the run number in bin count ahead of a temporary renumbering.

Parameters:
peaksstr

Name of peaks table.

update_scale_factor(peaks, value)[source]#

Update counting statistic refrence value for normalization

Parameters:
value: float

Monitor count or proton charge

garnet.reduction.plan module#

class Dumper(stream, default_style=None, default_flow_style=False, canonical=None, indent=None, width=None, allow_unicode=None, line_break=None, encoding=None, explicit_start=None, explicit_end=None, version=None, tags=None, sort_keys=True)[source]#

Bases: SafeDumper

Methods

add_implicit_resolver

add_multi_representer

add_path_resolver

add_representer

analyze_scalar

anchor_node

ascend_resolver

check_empty_document

check_empty_mapping

check_empty_sequence

check_resolver_prefix

check_simple_key

choose_scalar_style

close

descend_resolver

determine_block_hints

dispose

emit

expect_alias

expect_block_mapping

expect_block_mapping_key

expect_block_mapping_simple_value

expect_block_mapping_value

expect_block_sequence

expect_block_sequence_item

expect_document_end

expect_document_root

expect_document_start

expect_first_block_mapping_key

expect_first_block_sequence_item

expect_first_document_start

expect_first_flow_mapping_key

expect_first_flow_sequence_item

expect_flow_mapping

expect_flow_mapping_key

expect_flow_mapping_simple_value

expect_flow_mapping_value

expect_flow_sequence

expect_flow_sequence_item

expect_node

expect_nothing

expect_scalar

expect_stream_start

flush_stream

generate_anchor

ignore_aliases

increase_indent

need_events

need_more_events

open

prepare_anchor

prepare_tag

prepare_tag_handle

prepare_tag_prefix

prepare_version

process_anchor

process_scalar

process_tag

represent

represent_binary

represent_bool

represent_data

represent_date

represent_datetime

represent_dict

represent_float

represent_int

represent_list

represent_mapping

represent_none

represent_scalar

represent_sequence

represent_set

represent_str

represent_undefined

represent_yaml_object

resolve

serialize

serialize_node

write_double_quoted

write_folded

write_indent

write_indicator

write_line_break

write_literal

write_plain

write_single_quoted

write_stream_end

write_stream_start

write_tag_directive

write_version_directive

represent_list(data)[source]#
yaml_representers = {<class 'NoneType'>: <function SafeRepresenter.represent_none>, <class 'bool'>: <function SafeRepresenter.represent_bool>, <class 'bytes'>: <function SafeRepresenter.represent_binary>, <class 'datetime.date'>: <function SafeRepresenter.represent_date>, <class 'datetime.datetime'>: <function SafeRepresenter.represent_datetime>, <class 'dict'>: <function SafeRepresenter.represent_dict>, <class 'float'>: <function SafeRepresenter.represent_float>, <class 'int'>: <function SafeRepresenter.represent_int>, <class 'list'>: <function Dumper.represent_list>, <class 'set'>: <function SafeRepresenter.represent_set>, <class 'str'>: <function SafeRepresenter.represent_str>, <class 'tuple'>: <function SafeRepresenter.represent_list>, None: <function SafeRepresenter.represent_undefined>}[source]#
class ReductionPlan[source]#

Bases: object

Methods

generate_plan(instrument)

Create a template plan.

load_plan(filename)

Load a data reduction plan.

runs_list_to_string(runs)

Convert runs list to string with step notation.

runs_string_to_list(runs_str)

Convert runs string to list.

save_plan(filename[, set_output])

Save a data reduction plan.

set_output(filename)

Change the output directory and name.

template_integration(instrument)

Generate template integration plan.

template_normalization()

Generate template normalization plan.

template_parametrization()

Generate template parametrization plan.

check

validate_file

validate_plan

check(*args, **kwargs)[source]#
generate_plan(instrument)[source]#

Create a template plan.

Parameters:
instrumentstr

Beamline name.

load_plan(filename)[source]#

Load a data reduction plan.

Parameters:
filenamestr

yaml file of reduction plan.

runs_list_to_string(runs)[source]#

Convert runs list to string with step notation.

Parameters:
runslist

Integer run numbers.

Returns:
runs_strstr

Condensed notation for run numbers, including step notation.

runs_string_to_list(runs_str)[source]#

Convert runs string to list.

Parameters:
runs_strstr

Condensed notation for run numbers.

Returns:
runslist

Integer run numbers.

save_plan(filename, set_output=True)[source]#

Save a data reduction plan.

Parameters:
filenamestr

yaml file of reduction plan.

set_output(filename)[source]#

Change the output directory and name.

Parameters:
filenamestr

yaml file of reduction plan.

template_integration(instrument)[source]#

Generate template integration plan.

Parameters:
instrumentstr

Beamline name.

Returns:
paramsdict

Integration plan.

template_normalization()[source]#

Generate template normalization plan.

Parameters:
instrumentstr

Beamline name.

Returns:
paramsdict

Integration plan.

template_parametrization()[source]#

Generate template parametrization plan.

Parameters:
instrumentstr

Beamline name.

Returns:
paramsdict

Integration plan.

validate_file(fname, ext)[source]#
validate_plan()[source]#
class SubPlan(plan)[source]#

Bases: object

Attributes:
timer

Timer for the steps of this reduction.

Methods

append_name(file)

Update filename with identifier name.

get_diagnostic_file(name[, ext])

Diagnostic file.

get_diagnostic_path()

Diagnostic path.

get_output_file([ext])

Output file.

get_output_path()

Output path.

get_plot_file(name[, ext])

Plot file.

get_plot_path()

Plot path.

save_timing()

Save the timings of this process (in the diagnostics directory, which is cleared at the start of each reduction), for write_timing to combine.

write_timing()

Combine the timings saved by every process into <OutputName>_timing.txt, with each step's total over all runs and its mean per run, and remove them.

check

cleanup

create_directories

append_name(file)[source]#

Update filename with identifier name.

Parameters:
filestr

Original file name.

Returns:
output_filestr

File with updated name for identifier name.

check(*args, **kwargs)[source]#
cleanup()[source]#
create_directories()[source]#
get_diagnostic_file(name, ext='.nxs')[source]#

Diagnostic file.

Returns:
namestr

Name to use for file.

diag_filestr

Path name to save diagnostics file.

get_diagnostic_path()[source]#

Diagnostic path.

Returns:
diag_pathstr

Path to save diagnostics file.

get_output_file(ext='.nxs')[source]#

Output file.

Returns:
output_filestr

Output file.

get_output_path()[source]#

Output path.

Returns:
output_pathstr

Output path.

get_plot_file(name, ext='.png')[source]#

Plot file.

Returns:
namestr

Name to use for file.

plot_filestr

Path name to save plots file.

get_plot_path()[source]#

Plot path.

Returns:
plot_pathstr

Path to save plots file.

save_timing()[source]#

Save the timings of this process (in the diagnostics directory, which is cleared at the start of each reduction), for write_timing to combine.

property timer[source]#

Timer for the steps of this reduction.

write_timing()[source]#

Combine the timings saved by every process into <OutputName>_timing.txt, with each step’s total over all runs and its mean per run, and remove them.

check(value, op, other=None, message=None)[source]#

General-purpose assertion with flexible condition types.

Parameters:
valueany

Parameter.

opstr

Conditional operator.

otherTYPE, optional

Reference value. The default is None.

messagestr, optional

Invalid message. The default is None.

Raises:
ValueError

Invalud operator.

TypeError

Invalud conditional.

AssertionError

Invalid value.

Examples

check(3, ‘>’, 0) check(5, ‘in’, [1, 3, 5]) check(None, ‘is’, None) check(3, lambda v: v % 2 == 1)

delete_directory(path)[source]#

Fast removal of nonempty directories.

Parameters:
pathstr

Directory to remove.

load_YAML(filename)[source]#

Load reduction input file.

Parameters:
filenamestr

Output file name.

Returns:
outputdict

Parameters.

save_JSON(output, filename)[source]#

Save reduction output file.

Parameters:
outputstr

Name of file.

filenamestr

Output file name.

save_YAML(output, filename)[source]#

Save reduction output file.

Parameters:
outputdict

Parameters.

filenamestr

Output file name.

garnet.reduction.resolution module#

class PeakShapePlot(axis_labels=None)[source]#

Bases: object

The per-peak box + overlaid ellipse diagnostic figure – a spatial complement to _plot_resolution_diagnostics’s population-level obs-vs-pred scatter, showing the actual counts box each peak’s shape estimate(s) came from.

The figure is built once, for the first peak drawn, and reused for every later peak with the same layout (number of rows and ellipses), updating its images and ellipses rather than building and laying out a new figure each time. Call close when done.

Parameters:
axis_labelslist of (str, str), optional

Per-projection axis labels, one pair per dimension dropped (i.e. projections onto the other two). Defaults to the Q-sample convention [(“Q_1”, “Q_2”), (“Q_0”, “Q_2”), (“Q_0”, “Q_1”)]; pass e.g. local rotated-frame labels when “Q0”/”Q1”/”Q2” aren’t literal Q-sample components.

Methods

draw(samples, filename)

Draw one peak's figure and save it.

close

close()[source]#
draw(samples, filename)[source]#

Draw one peak’s figure and save it.

Parameters:
sampleslist of dict

Per-peak rows with keys “label”, “Q0”/”Q1”/”Q2” (dense grids, ij-indexed), “counts” (matching 3D array), and “ellipses”: a list of dict, each with “center” (3-tuple, same frame as “Q0”/”Q1”/”Q2”), “S” (3x3, containment-scale, same frame), “label”, “color”, and “linestyle” – overlaid on the counts image per sample, in the order given.

filenamestr

Output image path.

class ResolutionEllipsoid(peaks_ws, r_cut=inf, sig_noise_cut=5.0, min_peaks=10, scale_bounds=(0.5, 2.0), mosaic='isotropic', peak_shape_frame='s', lamda_cut_min=None, lamda_cut_max=None)[source]#

Bases: object

Methods

build_design_matrices(peak_indices)

Return per-peak design matrices for the nonlinear optimizer.

estimate_prior_sigmas()

Estimate prior width for S matrix and centroid parameters.

fit([fixed_instrumental, weighting, exclude])

Fit the variance parameters against the workspace's peak shapes.

fit_iterative([n_iter, clip_sigma, ...])

Repeatedly fit(), then hard-exclude peaks whose observed shape disagrees with the just-fitted model by more than clip_sigma robust-scaled residuals, and refit without them.

fit_mosaic_only(model[, weighting, loss, eps])

Fit only the mosaic term(s), holding every instrumental parameter (not just its shape, but its exact calibrated value) fixed at model's values -- unlike fit(fixed_instrumental=...), which still lets one overall scale factor float across all the instrumental columns together.

fit_sn_correction(model[, weighting, loss, ...])

Fit an S/N-dependent correction to an already-calibrated instrument model (e.g. from fit() on a strong, high-S/N peak subset), against self's own (presumably wider, noisier) peak-shape selection -- self.sig_noise_cut set lower than whatever produced model.

fit_trimmed([keep_frac, n_iter, weighting, ...])

Least-Trimmed-Squares-style concentration-step refit: keep only the keep_frac fraction of eligible peaks with the smallest residual under the current fit, hard-refit on exactly that subset, and repeat until the kept set stops changing (the classic Fast-LTS algorithm, Rousseeuw & Van Driessen 1999).

predict_roi_radii(r_cut[, margin, min_frac])

Per-peak box half-width for a second extract_peak_counts pass, sized from this model's predicted ellipsoid instead of the one global r_cut every peak used on the first pass.

renumber_by_size([n])

Renumber peak runNumbers 1..n ordered by model ellipsoid size.

robust_nnls(A, y[, loss, eps])

Bounded least squares (x >= 0) with a single robust (Cauchy-loss) least_squares solve -- no iterative MAD-rescaled reweighting.

robust_nnls_saturation(A, y, lamda[, loss, ...])

Like robust_nnls, but A's column 0 (sigma_gamma_i^2) and column 1 (sigma_nu_i^2) each additionally get scaled by their own x^2 / (1 + x^q)^(2/q) (x = lamda / lambda_0) before every residual evaluation -- lambda_0_gamma_i and lambda_0_nu_i fit independently and nonlinearly (sharing one q, the transition sharpness) alongside the (still nonnegative) linear sigma^2 parameters, in a single robust least_squares solve (no iterative reweighting).

set_variance_parameters(sigmas)

Initialize the model directly from known/prior sigma parameters, bypassing fit() -- e.g. instrument-characteristic divergence parameters (config.instruments.beamlines[...]["DivergenceParams"]) used to seed peak shapes before per-run peaks exist to fit an instrument-specific model against.

set_variance_parameters_deg(sigmas_deg)

Same as set_variance_parameters, but with the angular sigmas given in degrees instead of radians.

write_diagnostics_csv(filename)

Write the per-peak diagnostics table (Q, observed/predicted/error radii along the Stoica-Wilkinson axes, angles, S/N, outlier flag) to a CSV file.

apply

diagnostics

plot_diagnostics

predict_lab_S

predict_sample_S

write_resolution_parameters

apply()[source]#
build_design_matrices(peak_indices)[source]#

Return per-peak design matrices for the nonlinear optimizer.

Parameters:
peak_indiceslist of int
Returns:
dict{peak_index: (6, n_params) ndarray}

Maps each index to the matrix A such that vech6(S_lab) ≈ A @ x.

diagnostics()[source]#
estimate_prior_sigmas()[source]#

Estimate prior width for S matrix and centroid parameters.

Uses the peak shapes stored in the workspace (set by apply()) against the model predictions to compute RMS residuals for both the covariance and centroid.

fit(fixed_instrumental=None, weighting='sn_over_q2', exclude=None)[source]#

Fit the variance parameters against the workspace’s peak shapes.

Parameters:
fixed_instrumentaldict or None

Prior instrumental sigmas (radians, “sigma_dl_mod” dimensionless – same shape as set_variance_parameters expects), e.g. from a beamline’s characterized DivergenceParams. When given, the instrumental columns (gamma_i, nu_i, gamma_f, nu_f, dl_mod, dl_mod_b) are held at this prior shape and only a single overall scale factor plus the mosaic term(s) are fit, rather than letting all instrumental terms float independently – lambda_0 (if present in the prior) is baked into the design matrix directly rather than fit.

weightingstr, optional

How to weight each peak’s row: “sn_over_q2” (default) – signal/noise / Q^2. “sn” – signal/noise alone (no Q^2 down-weighting of the high-Q tail). “sn2” – signal/noise squared – leans further into strong, well-measured peaks than “sn” alone, letting the fit track whatever trend the strongest peaks show while suppressing weak/outlier-prone peaks’ influence, without hard-excluding any of them. “none” – every peak weighted equally. A diagnostic knob to check whether the Q^2 term is masking a real trend at the tails of the Q or wavelength range, where weight is otherwise systematically smaller/larger.

excludeset of int or None, optional

Peak indices (workspace row index, same as diagnostics()’s “i”) to hard-exclude from this fit entirely – e.g. peaks a prior iteration of fit_iterative flagged as sigma-clipped outliers. Unlike the internal Huber down-weighting (which still lets an outlier’s row pull the fit a little), these peaks contribute nothing.

fit_iterative(n_iter=5, clip_sigma=5.0, weighting='sn_over_q2', fixed_instrumental=None)[source]#

Repeatedly fit(), then hard-exclude peaks whose observed shape disagrees with the just-fitted model by more than clip_sigma robust-scaled residuals, and refit without them. Converges when a pass finds no new peaks to exclude (or after n_iter passes).

Unlike the Huber down-weighting already inside fit() (bounded but nonzero influence for every peak, every pass), this drops flagged peaks entirely and for good – useful when a genuine trend in the data is being obscured by a persistent minority of badly-integrated peaks that Huber alone doesn’t fully suppress.

Parameters:
n_iterint, optional

Maximum number of fit-then-clip passes.

clip_sigmafloat, optional

Peaks whose largest-axis relative residual (|obs_r - pred_r| / pred_r, max over the three Stoica-Wilkinson axes) is more than clip_sigma robust (MAD-based) scaled deviations above the median are excluded on the next pass. Large by default (5) – this is a coarse “drop the badly-integrated tail”, not a general-purpose robust estimator (fit() already does that internally).

weighting, fixed_instrumental

Passed through to every fit() call.

Returns:
modeldict

Same shape as fit()’s return, from the final pass.

excludedset of int

Peak indices excluded by the final pass.

fit_mosaic_only(model, weighting='sn', loss='cauchy', eps=1e-12)[source]#

Fit only the mosaic term(s), holding every instrumental parameter (not just its shape, but its exact calibrated value) fixed at model’s values – unlike fit(fixed_instrumental=…), which still lets one overall scale factor float across all the instrumental columns together.

Mosaic’s uncertainty in the ordinary joint fit is contaminated by its correlation with the instrumental terms (sigma_gamma_f in particular, r ~ -0.5 to -0.7 throughout this model’s development) – the joint Jacobian is close to rank-deficient in that direction, so mosaic and sigma_gamma_f can trade off against each other with little cost to the joint residual. Subtracting the instrumental prediction out of the data entirely (rather than fitting a shared scale for it) removes that trade-off by construction: mosaic is the only thing left free, so its resulting uncertainty is a much more trustworthy answer to “is mosaic measurable at all here” – at the cost of not propagating model’s own uncertainty into that number (this is a conditional/profile uncertainty, Var(mosaic | instrumental fixed at its point estimate), not the full marginal one).

Parameters:
modeldict

An already-fit self.model-shaped dict (needs “lambda_0_gamma_i”, “lambda_0_nu_i”, “q_i”, and “variance_parameters”).

weighting, loss, epsas in fit.
Returns:
mosaicdict

Same mosaic-only keys _label_variance_parameters would produce (e.g. “sigma_mosaic” for the isotropic model), plus a “_stderr” entry for each.

residual_normfloat
usedlist of int

Workspace peak indices used, same convention as fit.

fit_sn_correction(model, weighting='sn', loss='cauchy', eps=1e-12, a_bounds=(0.0, 0.999), s0_bounds=None)[source]#

Fit an S/N-dependent correction to an already-calibrated instrument model (e.g. from fit() on a strong, high-S/N peak subset), against self’s own (presumably wider, noisier) peak-shape selection – self.sig_noise_cut set lower than whatever produced model.

sigma_obs = f(sn) * sigma_calc, f(sn) = 1 - a*exp(-sn/s0) – -> 1 (unbiased) at high S/N; bounded below by 1-a as sn -> 0 (a peak-shape-fitting artifact at low counts systematically under- or overestimates the fitted radius, saturating rather than diverging). Equivalently S_obs = f(sn)^2 * S_calc.

model[“variance_parameters”]/[“lambda_0_gamma_i”]/ [“lambda_0_nu_i”] are held fixed rather than refit jointly with (a, s0): only the fixed model side is rescaled by f(sn) here, so there’s no degenerate escape the way there would be if the free linear sigma^2 parameters could also float (see the reverted joint-fit attempt this rescaling scheme replaced) – pushing (a, s0) to an extreme just makes the fixed prediction mismatch the (untouched) observed data worse, not trivially better.

Parameters:
modeldict

An already-fit self.model-shaped dict (needs “lambda_0_gamma_i”, “lambda_0_nu_i”, and “variance_parameters”).

weighting, loss, epsas in fit/robust_nnls_saturation.
a_boundstuple of float, optional

Bounds for a (dimensionless).

s0_boundstuple of float or None, optional

Bounds for s0 (same units as S/N). None (default) computes (1.0, max observed S/N) from this fit’s own retained peaks.

Returns:
afloat
s0float
residual_normfloat
resultOptimizeResult
usedlist of int

Workspace peak indices used, same convention as fit.

fit_trimmed(keep_frac=0.5, n_iter=15, weighting='sn_over_q2', fixed_instrumental=None)[source]#

Least-Trimmed-Squares-style concentration-step refit: keep only the keep_frac fraction of eligible peaks with the smallest residual under the current fit, hard-refit on exactly that subset, and repeat until the kept set stops changing (the classic Fast-LTS algorithm, Rousseeuw & Van Driessen 1999).

Unlike fit_iterative’s open-ended sigma-clipping – which estimates an outlier threshold (MAD-based robust scale) from the data itself, and MAD’s own breakdown point is exactly 50% – this fixes the retained fraction up front, so it tolerates contamination up to 1 - keep_frac by construction rather than by estimating a threshold that itself degrades as contamination approaches 50%.

Parameters:
keep_fracfloat, optional

Fraction of eligible peaks to retain each concentration step (e.g. 0.5 keeps the best half).

n_iterint, optional

Maximum number of concentration steps.

weighting, fixed_instrumental

Passed through to every fit() call.

Returns:
modeldict

Same shape as fit()’s return, from the final step.

keptset of int

Peak indices retained by the final step.

plot_diagnostics(filename)[source]#
predict_lab_S(peak_index)[source]#
predict_roi_radii(r_cut, margin=2.0, min_frac=0.3)[source]#

Per-peak box half-width for a second extract_peak_counts pass, sized from this model’s predicted ellipsoid instead of the one global r_cut every peak used on the first pass.

A single global r_cut is simultaneously too large for compact/ weak peaks (their box is mostly background, biasing the direct centroid/covariance moment low via estimate_peak_shapes’s SNR shrinkage) and potentially too small for genuinely broad ones. Resizing each peak’s box from the model’s own prediction targets both at once. Call after fit().

Parameters:
r_cutfloat

Original (first-pass) box half-width – also the upper bound returned here, since it’s presumably sized to safely avoid absorbing neighboring peaks.

marginfloat

Multiple of the predicted largest containment radius used as the new box half-width, leaving room around the peak for estimate_peak_shapes’s per-box background estimate.

min_fracfloat

Lower bound on the returned radius, as a fraction of r_cut – guards against a pathologically tight box where the model predicts a near-zero radius.

Returns:
radiindarray, shape (n_peaks,)

Per-peak box half-width, in the same order as the peaks workspace.

predict_sample_S(peak_index)[source]#
renumber_by_size(n=None)[source]#

Renumber peak runNumbers 1..n ordered by model ellipsoid size.

Each peak’s predicted largest principal radius is computed from the resolution model. Peaks are sorted from smallest to largest and divided into n equal-count bins; all peaks in the same bin get the same runNumber (1 = smallest, n = largest). When n is None each peak gets its own unique rank (equivalent to n = n_peaks). Call after fit() (and optionally apply()).

Parameters:
nint or None

Number of size bins. Defaults to one bin per peak.

Returns:
radii_binsndarray, shape (n_bins,)

Mean largest radius of the peaks assigned to each bin, ordered from smallest (bin 1) to largest (bin n).

robust_nnls(A, y, loss='cauchy', eps=1e-12)[source]#

Bounded least squares (x >= 0) with a single robust (Cauchy-loss) least_squares solve – no iterative MAD-rescaled reweighting.

Replaces a previous plain-NNLS (Lawson-Hanson) reweighting scheme: NNLS’s active-set solve pins parameters to exactly zero, and which parameters get pinned could flip between iterations as the Huber weights changed, so the fit could oscillate instead of converging. least_squares with method=”trf” enforces x >= 0 continuously instead, and supports a robust loss directly – one pass, no reweighting loop.

robust_nnls_saturation(A, y, lamda, loss='cauchy', eps=1e-12, lambda_0_bounds=(0.01, 100.0), q_bounds=(0.1, 20.0))[source]#

Like robust_nnls, but A’s column 0 (sigma_gamma_i^2) and column 1 (sigma_nu_i^2) each additionally get scaled by their own x^2 / (1 + x^q)^(2/q) (x = lamda / lambda_0) before every residual evaluation – lambda_0_gamma_i and lambda_0_nu_i fit independently and nonlinearly (sharing one q, the transition sharpness) alongside the (still nonnegative) linear sigma^2 parameters, in a single robust least_squares solve (no iterative reweighting).

Parameters:
A, yas in robust_nnls – A’s columns 0/1 must be the

raw (lambda_0=None) sigma_gamma_i/sigma_nu_i columns from _model_design_lab, not yet saturation-scaled.

lamdandarray, shape (A.shape[0],)

Wavelength (Angstrom) for each row of A.

lambda_0_boundstuple of float, optional

Bounds for lambda_0_gamma_i/lambda_0_nu_i (Angstrom). Deliberately wide – a lambda_0 much larger than the data’s wavelength range recovers pure linear growth, and one much smaller recovers a wavelength-independent constant, so the data is free to prefer either limit.

q_boundstuple of float, optional

Bounds for the shared saturation sharpness exponent q.

Returns:
xndarray

Fitted sigma^2 parameters (same meaning/order as robust_nnls, i.e. without lambda_0_gamma_i/ lambda_0_nu_i/q).

lambda_0_gamma_ifloat
lambda_0_nu_ifloat
qfloat
residual_normfloat
resultOptimizeResult

result.x/result.jac have three extra (final) entries/ columns, for lambda_0_gamma_i, lambda_0_nu_i, then q, for _parameter_covariance.

set_variance_parameters(sigmas)[source]#

Initialize the model directly from known/prior sigma parameters, bypassing fit() – e.g. instrument-characteristic divergence parameters (config.instruments.beamlines[…][“DivergenceParams”]) used to seed peak shapes before per-run peaks exist to fit an instrument-specific model against.

Parameters:
sigmasdict

Same keys _label_variance_parameters returns for the current self.mosaic model (e.g. for “isotropic”: sigma_gamma_i, sigma_nu_i, sigma_gamma_f, sigma_nu_f, sigma_dl_mod, sigma_mosaic).

Returns:
modeldict

Same shape fit() sets on self.model, minus the fit-only keys “residual_norm”, “used_peaks”, “robust_weights”.

set_variance_parameters_deg(sigmas_deg)[source]#

Same as set_variance_parameters, but with the angular sigmas given in degrees instead of radians. “sigma_dl_mod”/ “sigma_dl_mod_b” are dimensionless (sigma_lambda/lambda, see _model_design_lab) rather than an angle, so they are passed through unchanged.

Parameters:
sigmas_degdict

Same keys as set_variance_parameters, with every key but “sigma_dl_mod”/”sigma_dl_mod_b”/”lambda_0_gamma_i”/ “lambda_0_nu_i”/”q_i” in degrees. Those five are passed through unchanged – none of them is an angle.

Returns:
modeldict

Same as set_variance_parameters.

write_diagnostics_csv(filename)[source]#

Write the per-peak diagnostics table (Q, observed/predicted/error radii along the Stoica-Wilkinson axes, angles, S/N, outlier flag) to a CSV file. Call after fit().

write_resolution_parameters(filename)[source]#
stoica_wilkinson_transform(two_theta, phi)[source]#

Local instrument-frame resolution axes (n, u, v) in lab (Q-lab) frame.

Parameters:
two_thetafloat

Scattering angle in radians.

phifloat

Azimuthal angle in radians.

Returns:
Tndarray, shape (3, 3)

Rows are the orthonormal n, u, v directions.

garnet.reduction.sample module#

class SampleMaterial(plan, filename)[source]#

Bases: object

Methods

get_material

get_sample_material

get_shape

save_ellipsoid_stl

set_material

set_sample

set_shape

verify_chemical_formula

get_material()[source]#
get_sample_material()[source]#
get_shape()[source]#
save_ellipsoid_stl(filename='/tmp/ellipsoid.stl')[source]#
set_material(material)[source]#
set_sample(sample)[source]#
set_shape(filename)[source]#
verify_chemical_formula(formula)[source]#

garnet.reduction.ub module#

class Optimization(peaks, tol=0.1)[source]#

Bases: object

Methods

get_UB()

Current UB matrux.

get_lattice_parameters()

Current lattice parameters.

get_orientation_angles()

Current orientation angles.

optimize_lattice(cell[, n_cycles, sigma_cut])

Refine the orientation and lattice parameters under constraints.

optimize_lattice_only(cell[, n_cycles, ...])

Refine lattice parameters under cell constraints with U fixed.

residual(x, hkl, Q, fun[, W, weights])

Optimization residual function.

B_matrix

U_matrix

cubic

fixed

hexagonal

monoclinic

orthorhombic

rhombohedral

tetragonal

triclinic

whiten_weight_matrix

B_matrix(a, b, c, alpha, beta, gamma)[source]#
U_matrix(phi, theta, omega)[source]#
cubic(x)[source]#
fixed(x)[source]#
get_UB()[source]#

Current UB matrux.

Returns:
UB2d-array

UB-matrix.

get_lattice_parameters()[source]#

Current lattice parameters.

Returns:
a, b, cfloat

Lattice constants in angstroms.

alpha, beta, gammafloat

Lattice angles in degrees.

get_orientation_angles()[source]#

Current orientation angles.

Returns:
phifloat

Rotation axis azimuthal angle in radians.

thetafloat

Rotation axis polar angle in radians.

omegafloat

Rotation angle in radians.

hexagonal(x)[source]#
monoclinic(x)[source]#
optimize_lattice(cell, n_cycles=5, sigma_cut=3.0)[source]#

Refine the orientation and lattice parameters under constraints.

Iterates least-squares refinement with I/σ weighting and σ-clipping outlier rejection. After each cycle peaks whose unweighted residual norm exceeds sigma_cut times the median are removed, and the next cycle starts from the previous solution.

Parameters:
cellstr

Lattice centering to constrain parameters.

n_cyclesint, optional

Maximum number of sigma-clipping cycles. The default is 5.

sigma_cutfloat, optional

Outlier rejection threshold in units of the median residual norm. The default is 3.0.

optimize_lattice_only(cell, n_cycles=5, sigma_cut=3.0)[source]#

Refine lattice parameters under cell constraints with U fixed.

The crystal orientation (U matrix) is held constant; only the free lattice parameters allowed by cell are optimized.

Parameters:
cellstr

Lattice system to constrain parameters.

n_cyclesint, optional

Maximum number of sigma-clipping cycles. The default is 5.

sigma_cutfloat, optional

Outlier rejection threshold in units of the median residual norm. The default is 3.0.

orthorhombic(x)[source]#
residual(x, hkl, Q, fun, W=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), weights=None)[source]#

Optimization residual function.

Parameters:
xlist

Parameters.

hkllist

Miller indices.

Qlist

Q-sample vectors.

funfunction

Lattice constraint function.

W: 3x3-array

Weight matrix.

weights1d-array, optional

Per-peak scalar weights (e.g. I/σ). The default is None (uniform).

Returns:
residuallist

Least squares residuals.

rhombohedral(x)[source]#
tetragonal(x)[source]#
triclinic(x)[source]#
whiten_weight_matrix(Q)[source]#
class RefineSingleCrystalGoniometer(peaks, tol=0.12, cell='Triclinic', n_iter=1)[source]#

Bases: object

Methods

get_lattice_parameters()

Current lattice parameters.

get_orientation_angles()

Current orientation angles.

optimize_lattice(cell)

Refine the orientation and lattice parameters under constraints.

residual(x, peak_dict, func)

Optimization residual function.

B_matrix

U_matrix

calculate_goniometer

cubic

fixed

hexagonal

monoclinic

orthorhombic

rhombohedral

tetragonal

triclinic

B_matrix(a, b, c, alpha, beta, gamma)[source]#
U_matrix(phi, theta, omega)[source]#
calculate_goniometer(omega, chi, phi)[source]#
cubic(x)[source]#
fixed(x)[source]#
get_lattice_parameters()[source]#

Current lattice parameters.

Returns:
a, b, cfloat

Lattice constants in angstroms.

alpha, beta, gammafloat

Lattice angles in degrees.

get_orientation_angles()[source]#

Current orientation angles.

Returns:
phifloat

Rotation axis azimuthal angle in radians.

thetafloat

Rotation axis polar angle in radians.

omegafloat

Rotation angle in radians.

hexagonal(x)[source]#
monoclinic(x)[source]#
optimize_lattice(cell)[source]#

Refine the orientation and lattice parameters under constraints.

Parameters:
cellstr

Lattice centering to constrain paramters.

orthorhombic(x)[source]#
residual(x, peak_dict, func)[source]#

Optimization residual function.

Parameters:
xlist

Parameters.

peak_dictdictionary

Goniometer angles, Q-lab vectors, Miller indices. .

funcfunction

Lattice constraint function.

Returns:
residuallist

Least squares residuals.

rhombohedral(x)[source]#
tetragonal(x)[source]#
triclinic(x)[source]#
class Reorient(peaks, UB_ref, crystal_system='Triclinic', lattice_system=None)[source]#

Bases: object

Methods

minimize(transforms, P[, tol])

Search the crystal system's full point group (proper and improper) composed with the axis-permutation P, keeping only the combinations that stay a valid, proper (det > 0) transform, for whichever best matches UB_ref.

resolve_axis_ambiguity(crystal_system)

Find the axis permutation (plus, if needed, a compensating single axis flip) of self.UB whose (a, b, c, alpha, beta, gamma) best matches self.UB_ref's, independent of crystal system (this is a general permutation-and-compare search, not specific cell-type logic).

cell_symmetry_matrices

cell_symmetry_matrices(crystal_system, lattice_system)[source]#
minimize(transforms, P, tol=0.12)[source]#

Search the crystal system’s full point group (proper and improper) composed with the axis-permutation P, keeping only the combinations that stay a valid, proper (det > 0) transform, for whichever best matches UB_ref.

P alone may be improper (an odd permutation, det = -1, e.g. a plain swap of two axes) – composing it with the point group’s proper rotations only would then always stay improper, since proper x proper = proper and proper x improper = improper. Point group elements are filtered by the combined transform’s determinant here (not on their own beforehand) so the correct subset – proper or improper – is used in either case.

resolve_axis_ambiguity(crystal_system)[source]#

Find the axis permutation (plus, if needed, a compensating single axis flip) of self.UB whose (a, b, c, alpha, beta, gamma) best matches self.UB_ref’s, independent of crystal system (this is a general permutation-and-compare search, not specific cell-type logic). Only Triclinic/Monoclinic/Orthorhombic actually need it – see __init__.

A plain axis permutation is improper (det < 0) for an odd permutation (e.g. a single a<->b swap) – left uncorrected, that would flip the resulting UB left-handed. Real UB determinations are always right-handed, so an odd permutation is only ever realistic paired with a compensating sign flip on one axis; which axis is tried explicitly (flipping changes 2 of the 3 angles, not just a sign, so the wrong choice would fail the angle match) and only kept if it actually restores det > 0.

Returns:
Pndarray, shape (3, 3)

Proper (det > 0) axis-permutation/sign-flip matrix.

class UBModel(peaks)[source]#

Bases: object

Methods

calculate_hkl()

Calculate hkl values without rounding.

copy_UB(workspace)

Copy UB to another workspace.

determine_UB_from_conventional_cell(a, b, c, ...)

Determine UB with prior known lattice parameters.

determine_UB_with_lattice_parameters(a, b, ...)

Determine UB with prior known lattice parameters.

determine_UB_with_primitive_cell(min_d, max_d)

Determine UB with primitive lattice using min/max lattice constant.

generate_lattice_transforms(cell)

Obtain possible transforms compatabile with a unit cell lattice.

get_UB()

Current UB matrix.

get_center_uncertainty(hkl[, min_frac, max_frac])

Propagate lattice-parameter uncertainties to a peak's Q-space center uncertainty.

get_lattice_parameter_uncertanties()

Get the errors in the lattice constants from the oriented lattice.

get_lattice_parameters()

Current lattice parameters.

get_max_d_spacing()

Obtain the maximum d-spacing from the oriented lattice.

has_UB()

Check if peaks table has a UB determined.

index_peaks([tol, sat_tol, mod_vec_1, ...])

Index the peaks and calculate the lattice parameter uncertainties.

load_UB(filename[, run_number])

Load UB from file.

possible_conventional_cells([max_error, ...])

List possible conventional cells.

refine_UB_with_constraints(cell[, tol])

Refine UB with constraints corresponding to lattice system.

refine_UB_without_constraints([tol, sat_tol])

Refine UB with unconstrained lattice parameters.

refine_U_only(a, b, c, alpha, beta, gamma)

Refine the U orientation only.

save_UB(filename)

Save UB to file.

select_cell(number[, tol])

Transform to conventional cell using form number.

select_type(cell, centering[, tol])

Transform to conventional cell using cell and centering type.

transform_lattice(transform[, tol])

Apply a cell transformation to the lattice.

calculate_transform_extents

centering_matrix

convert_conventional_to_primitive

get_primitive_cell_length_range

shortest_reciprocal_spacing

transform_conventional_to_primitive

transform_primitive_to_conventional

calculate_hkl()[source]#

Calculate hkl values without rounding.

calculate_transform_extents(centering)[source]#
centering_matrix(centering)[source]#
convert_conventional_to_primitive(a, b, c, alpha, beta, gamma, centering)[source]#
copy_UB(workspace)[source]#

Copy UB to another workspace.

Parameters:
workspacefloat

Target workspace to copy the UB to.

determine_UB_from_conventional_cell(a, b, c, alpha, beta, gamma, centering='P', tol=0.2)[source]#

Determine UB with prior known lattice parameters.

Parameters:
a, b, cfloat

Lattice constants in angstroms.

alpha, beta, gammafloat

Lattice angles in degrees.

centeringstr,

Conventional cell reflection condition.

tolfloat, optional

Indexing tolerance. The default is 0.1.

determine_UB_with_lattice_parameters(a, b, c, alpha, beta, gamma, tol=0.2)[source]#

Determine UB with prior known lattice parameters.

Parameters:
a, b, cfloat

Lattice constants in angstroms.

alpha, beta, gammafloat

Lattice angles in degrees.

tolfloat, optional

Indexing tolerance. The default is 0.1.

determine_UB_with_primitive_cell(min_d, max_d, tol=0.15)[source]#

Determine UB with primitive lattice using min/max lattice constant.

Parameters:
min_dfloat

Minimum lattice parameter in ansgroms.

max_dfloat

Maximum lattice parameter in angstroms.

tolfloat, optional

Indexing tolerance. The default is 0.15.

generate_lattice_transforms(cell)[source]#

Obtain possible transforms compatabile with a unit cell lattice.

Parameters:
cellstr

Latttice system.

Returns:
transformsdict

Transform dictionary with symmetry operation as key.

get_UB()[source]#

Current UB matrix.

Returns:
UB2d-array

UB-matrix.

get_center_uncertainty(hkl, min_frac=0.01, max_frac=0.05)[source]#

Propagate lattice-parameter uncertainties to a peak’s Q-space center uncertainty.

Parameters:
hklsequence of float

Miller indices of the peak.

min_fracfloat, optional

Floor on each lattice-parameter uncertainty as a fraction of its value. Default 0.01.

max_fracfloat, optional

Cap on each lattice-parameter uncertainty as a fraction of its value. Default 0.05.

Returns:
sigma_cfloat

Norm of the propagated Q-space center uncertainty (inverse angstroms), or 0.0 if lattice information is unavailable.

get_lattice_parameter_uncertanties()[source]#

Get the errors in the lattice constants from the oriented lattice.

Returns:
errorslist

List of errors in lattice constants.

get_lattice_parameters()[source]#

Current lattice parameters.

Returns:
a, b, cfloat

Lattice constants in angstroms.

alpha, beta, gammafloat

Lattice angles in degrees.

get_max_d_spacing()[source]#

Obtain the maximum d-spacing from the oriented lattice.

Returns:
d_maxfloat

Maximum d-spacing.

get_primitive_cell_length_range(centering)[source]#
has_UB()[source]#

Check if peaks table has a UB determined.

index_peaks(tol=0.15, sat_tol=None, mod_vec_1=[0, 0, 0], mod_vec_2=[0, 0, 0], mod_vec_3=[0, 0, 0], max_order=0, cross_terms=False)[source]#

Index the peaks and calculate the lattice parameter uncertainties.

Parameters:
tolfloat, optional

Indexing tolerance. The default is 0.1.

sat_tolfloat, optional

Satellite indexing tolerance. The default is None.

mod_vec_1, mod_vec_2, mod_vec_3list, optional

Modulation vectors. The default is [0,0,0].

max_orderint, optional

Maximum order greater than zero for satellites. The default is 0.

cross_termsbool, optional

Include modulation cross terms. The default is False.

Returns:
indexinglist

Result of indexing including number indexed and errors.

load_UB(filename, run_number=None)[source]#

Load UB from file.

Parameters:
filenamestr

Name of UB file with extension .mat.

run_numberstr, optional

Run number to replace starred expression.

possible_conventional_cells(max_error=0.2, permutations=True)[source]#

List possible conventional cells.

Parameters:
max_errorfloat, optional

Max scalar error to report form numbers. The default is 0.2.

permutationsbool, optional

Allow permutations of the lattice. The default is True.

Returns:
valslist

List of form results.

refine_UB_with_constraints(cell, tol=0.1)[source]#

Refine UB with constraints corresponding to lattice system.

Lattice system

Lengths

Angles

Cubic

\(a=b=c\)

\(α=β=γ=90\)

Hexagonal

\(a=b\)

\(α=β=90,γ=120\)

Rhombohedral

\(a=b=c\)

\(α=β=γ\)

Tetragonal

\(a=b\)

\(α=β=γ=90\)

Orthorhombic

None

\(α=β=γ=90\)

Monoclinic

None

\(α=γ=90\)

Triclinic

None

None

Parameters:
cellfloat

Lattice system.

tolfloat, optional

Indexing tolerance. The default is 0.1.

refine_UB_without_constraints(tol=0.1, sat_tol=None)[source]#

Refine UB with unconstrained lattice parameters.

Parameters:
tolfloat, optional

Indexing tolerance. The default is 0.1.

sat_tolfloat, optional

Satellite indexing tolerance. The default is None.

refine_U_only(a, b, c, alpha, beta, gamma)[source]#

Refine the U orientation only.

Parameters:
a, b, cfloat

Lattice constants in angstroms.

alpha, beta, gammafloat

Lattice angles in degrees.

save_UB(filename)[source]#

Save UB to file.

Parameters:
filenamestr

Name of UB file with extension .mat.

select_cell(number, tol=0.1)[source]#

Transform to conventional cell using form number.

Parameters:
numberint

Form number.

tolfloat, optional

Indexing tolerance. The default is 0.1.

select_type(cell, centering, tol=0.1)[source]#

Transform to conventional cell using cell and centering type.

Parameters:
cellstr

Cell type.

centeringstr

Centering type.

tolfloat, optional

Indexing tolerance. The default is 0.1.

shortest_reciprocal_spacing(centering, max_index=2)[source]#
transform_conventional_to_primitive(centering)[source]#
transform_lattice(transform, tol=0.1)[source]#

Apply a cell transformation to the lattice.

Parameters:
transform3x3 array-like

Transform to apply to hkl values.

tolfloat, optional

Indexing tolerance. The default is 0.1.

transform_primitive_to_conventional(centering)[source]#
write_ub_info(info_file, run, min_d, opt, ub)[source]#

Write a lattice-refinement diagnostic text file.

Parameters:
info_filestr

Output text file path.

runint

Run number, for the report header.

min_dfloat

Resolution (minimum d-spacing) used for peak prediction.

optOptimization

Lattice optimization (already optimize_lattice()-refined), for the count of peaks used.

ubUBModel

UB model (already refined), for the lattice parameters and their uncertainties.

Module contents#