LegendSpecFits

Fitting routines for the analysis of germanium detector spectra in the LEGEND experiment.

The package provides gamma peak shapes and binned Poisson likelihood fits for energy calibration and resolution curves, charge-trapping correction, A/E and LQ pulse shape discrimination, SiPM calibration, filter parameter optimization, and goodness-of-fit tests. Most fit functions return a tuple (result, report) of the fit results and the data needed to plot them with LegendMakie.lplot.

BinnedPoissonProcess(intensity::AbstractMeasure, edges::AbstractVector{<:Real}) <: MeasureBase.AbstractMeasure

Binned view of a Poisson point process with the given intensity measure: its variates are vectors of non-negative integer bin counts for the given bin edges (logdensityof assumes valid counts, binned_poisson_likelihood validates them). The count in each bin is Poisson-distributed, with the intensity mass in the bin as expectation, computed via analytic MeasureBase.smf differences where available (all spectral shapes, distributed over weighted measures and superpositions) and via Simpson's rule otherwise. Distribution intensities are converted to measures automatically.

CalFuncDict

Alias for Dict{Int,PolCalFunc{2,Float64}}: linear calibration functions (PolCalFunc with two coefficients) keyed by an integer identifier.

ExpDecay(μ, τ) <: SpectralShape

Exponential shape with density exp(-(x - μ) / τ), unit height at μ and decaying towards larger values. Infinite total mass over ℝ.

GaussPeak(μ, σ) <: SpectralShape

Gaussian peak shape with mean μ, standard deviation σ > 0 and unit total mass.

GaussianStep(μ, σ) <: SpectralShape

Gaussian-smoothed step shape with unit plateau height below μ, erfc((x - μ) / (√2 σ)) / 2. Infinite total mass over ℝ.

HighETailPeak(μ, σ, θ) <: SpectralShape

Exponentially modified Gaussian peak shape with unit total mass, Gaussian mean μ and standard deviation σ > 0, and an exponential tail of scale θ >= 0 toward higher energies, the mirror image of LowETailPeak. The shape underlying highEtail_peakshape.

HistLogLike(h::Histogram{<:Real,1})
(ll::HistLogLike)(f_fit)

Calculates the binned Poisson log-likelihood of the histogram h, with all terms that do not depend on the model computed once.

The constructor stores the bin centers, bin widths, counts and the terms loggamma(k + 1) of the counts k. Calling ll(f_fit) returns the same value as hist_loglike(f_fit, h). A HistLogLike is constructed once and reused when the likelihood is evaluated repeatedly, for example inside an optimizer.

The expected count in each bin is bin_width * f_fit(bin_center). f_fit must accept all bin centers of h.

Arguments

• h::Histogram{<:Real,1}: Histogram whose counts are compared with the model.
• f_fit: Callable that returns the density at a position x.

Returns

The log-likelihood as a number. It is -Inf if a bin count is negative or not an integer. A negative expected count throws a DomainError, as Poisson(λ) does.

LowETailPeak(μ, σ, θ) <: SpectralShape

Exponentially modified Gaussian peak shape with unit total mass, Gaussian mean μ and standard deviation σ > 0, and an exponential tail of scale θ >= 0 toward lower energies. The shape underlying lowEtail_peakshape.

PolCalFunc(p0, p1, ...)

Callable polynomial x -> p0 + p1 * x + p2 * x^2 + ... with the given coefficients, ordered from the constant term upward. The coefficients are stored in params and can be of any Number type, including Measurement and Unitful.Quantity.

Examples

f = PolCalFunc(1.0, 2.0)  # 1 + 2x
f(3.0)                    # => 7.0

PolynomialShape(p0, p1, ...; x0 = 0) <: SpectralShape

Polynomial shape with density p0 + p1 (x - x0) + p2 (x - x0)^2 + .... The coefficients must keep the density non-negative on the fit region. An expansion point x0 near the region of interest improves the numerical conditioning of the coefficients.

abstract type SpectralShape <: BAT.BATMeasure

Supertype of the spectral shape measures.

_fit_fwhm_ft(e_grid, e_grid_ft, rt, min_e, max_e, rel_cut_fit; n_bins, peak, window, ft_fwhm_tol)

Implements the method of fit_fwhm_ft without qdrift: the energies are not corrected.

_get_model_counts(f_fit::Base.Callable, v_ml, bin_centers, bin_widths)

Returns the expected counts in each bin: the bin width times the model f_fit(x, v_ml) evaluated at the bin center.

_prepare_data(h::Histogram{<:Real,1})

Returns the tuple (counts, bin_widths, bin_centers) of the one-dimensional histogram h.

aoe_compton_background_peakshape(
    x::Real, μ::Real, σ::Real,
    background::Real, δ::Real
)

Evaluates the background part of aoe_compton_peakshape at x: an exponentially modified Gaussian (mirrored, tail toward lower A/E) with amplitude background, mean μ, standard deviation σ and scale δ.

aoe_compton_peakshape(
    x::Real, μ::Real, σ::Real, n::Real,
    background::Real, δ::Real
)

Evaluates the A/E distribution of events in a Compton region at x: a Gaussian single-site-event (SSE) peak with mean μ, standard deviation σ and amplitude n, plus a multi-site-event (MSE) background background modeled as an exponentially modified Gaussian with scale δ that extends toward lower A/E.

aoe_compton_peakshape_components(fit_func::Symbol; background_center = nothing)
aoe_compton_peakshape_components(fit_func::Symbol, v::NamedTuple; background_center = v.μ)

Returns the component functions of the A/E peak shape fit_func (:aoe_one_bck or :aoe_two_bck) as a NamedTuple with the fields

• funcs: component functions f(x, v): f_sig and f_bck for :aoe_one_bck; f_sig, f_bck_one and f_bck_two for :aoe_two_bck. In the method with v, these are functions of x only, with the parameters v applied.
• labels, colors, linestyles: plot label, color and line style of each component.

The components are used in fit reports and plot recipes.

aoe_compton_signal_peakshape(
    x::Real, μ::Real, σ::Real, n::Real
)

Evaluates the signal part of aoe_compton_peakshape at x: a Gaussian single-site-event peak with mean μ, standard deviation σ and amplitude n.

array_to_tuple(a::AbstractArray, as_nt::NamedTuple)

Returns a NamedTuple with the keys of as_nt and the values of a, in order.

autocal_energy(
    e::AbstractArray{<:Real}, photon_lines::Vector{<:Unitful.Energy{<:Real}};
    mode::Symbol = :fit, min_e::Real = 100, max_e::Real = maximum(e),
    max_e_binning_quantile::Real = 0.5, σ::Real = 2.0, threshold::Real = 50.0,
    min_n_peaks::Int = length(photon_lines), max_n_peaks::Int = 4 * length(photon_lines),
    α::Real = 0.01, rtol::Real = 5e-3
)

Determines a linear energy calibration constant from uncalibrated energies e by matching peaks in the spectrum to the known photon energies photon_lines.

The uncalibrated spectrum is histogrammed with the Freedman-Diaconis bin width of the values between min_e and the max_e_binning_quantile quantile of e. The peak matching is done by RadiationSpectra.calibrate_spectrum (mode = :fit) or RadiationSpectra.determine_calibration_constant_through_peak_ratios (mode = :ratio).

Keywords

• mode::Symbol = :fit: Matching method, :fit or :ratio. Other values throw an ArgumentError.
• min_e::Real = 100: Lower end of the histogrammed range (uncalibrated units).
• max_e::Real = maximum(e): Upper end of the histogrammed range.
• max_e_binning_quantile::Real = 0.5: Quantile of e that limits the data used to choose the bin width.
• σ::Real = 2.0, threshold::Real = 50.0: Expected peak width (in bins) and detection threshold of the peak finder.
• min_n_peaks::Int, max_n_peaks::Int: Minimum and maximum number of peaks to search for.
• α::Real = 0.01, rtol::Real = 5e-3: Significance level and relative tolerance of the peak matching.

Returns

A tuple (result, report) in which both are the same NamedTuple with the fields

• f_calib: function that multiplies uncalibrated values by the calibration constant, giving values in keV.
• h_cal: histogram of the calibrated energies from 0 to 3000 keV in 0.5 keV bins.
• h_uncal: histogram of the uncalibrated energies.
• c: calibration constant in keV per uncalibrated unit.
• peak_positions: positions of the identified peaks in uncalibrated units.
• threshold: peak finder threshold that was used.

background_peakshape(
    x::Real, μ::Real, σ::Real,
    step_amplitude::Real, background::Real;
    background_slope::Real = 0.0, background_exp = 0.0, background_center::Real = μ
)

Evaluates the background part of a gamma peak at x:

\[A\,\mathrm{S}(x; \mu, \sigma) + b\,e^{-d\,(x - x_c)} + s\,(x - x_c)\,,\]

where $A$ is step_amplitude, $b$ is background, $d$ is background_exp, $s$ is background_slope and $x_c$ is background_center. The step $\mathrm{S}$ is a Gaussian-smoothed step with width σ that falls from 1 to 0 across μ with increasing x (more background below the peak). With the default keywords the exponential term reduces to the constant background.

binned_poisson_likelihood(model, h::Histogram{<:Real,1})

Constructs the binned Poisson likelihood of the observed histogram h, i.e. the likelihood of BinnedPoissonProcess(model(p), h.edges) given the bin counts h.weights. model is a Markov kernel that maps parameters to an intensity measure (e.g. a @pf property function combining shape measures, see GaussPeak).

A fixed intensity measure instead of model yields a parameter-independent likelihood.

The bin counts of h must be non-negative integers: a weighted histogram is not Poisson-distributed data and requires a different statistical treatment.

chi2fit(f_fit::Function, x, y; pull_t, v_init, lower_bound, upper_bound, pseudo_prior, uncertainty, correlated)
chi2fit(f_fit, x, y, yerr; kwargs...)
chi2fit(f_fit, x, y, yerr, xerr; kwargs...)
chi2fit(n_poly::Int, x, y; kwargs...)
chi2fit(f_outer::Function, n_poly::Int, x, y; kwargs...)

Fits the model f_fit to the data (x, y) by minimizing a chi-squared function with the Minuit minimizer (migrad).

x and y are vectors of numbers or Measurements. The uncertainty of each x is propagated through the model with forward-mode automatic differentiation and added in quadrature to the uncertainty of y. If all uncertainties are zero, all points have equal weight. Separate yerr (and xerr) vectors are converted to Measurements.

Arguments

• f_fit: Model function f(x, p1, p2, ...). The number of fit parameters is first(methods(f_fit)).nargs - 2, so f_fit needs a single method with explicit parameter arguments.
• n_poly::Int: Fits a polynomial of degree n_poly instead, with parameters ordered from the constant term upward (see PolCalFunc). With f_outer, the model is f_outer(poly(x)).

Keywords

• pull_t::Vector{<:NamedTuple}: One entry per parameter. A NamedTuple (mean, std) adds a Gaussian pull term ((p - mean) / std)^2 to the chi-squared function; an empty NamedTuple() adds none. The default adds no pull terms.
• v_init::Vector{<:Real}: Initial parameter values. Default is one for every parameter.
• lower_bound::Vector{<:Real}, upper_bound::Vector{<:Real}: Parameter limits. Default is no limit.
• pseudo_prior::Union{ContinuousMultivariateDistribution, Nothing} = nothing: If given, the minimization runs in a transformed space where one unit is one prior standard deviation, which keeps parameters in the support of the prior.
• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit from the Hessian of the chi-squared function.
• correlated::Bool = true: If true, the returned parameters carry their full covariance, so that uncertainties of derived quantities, including report.f_fit, account for the parameter correlations. If false, the parameters are independent Measurements with the same marginal uncertainties, which overestimates the uncertainty of derived quantities when the parameters are correlated.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• par: best-fit parameters as Measurements. The uncertainty is NaN if uncertainty is false or the fit did not converge.
• converged: whether the minimizer converged. A warning is issued if not.
• gof: goodness of fit, only present if the uncertainties were calculated, as a NamedTuple with converged, pvalue, chi2min, dof, covmat and residuals_norm. The number of degrees of freedom counts each pull term as one data point; pvalue is NaN if it is zero.

report contains par, f_fit (the fitted model as a function of x), x, y and gof, for plotting with LegendMakie.lplot.

ctc_aoe(
    aoe_all::Vector{<:Real}, ecal_all::Vector{<:Unitful.Energy{<:Real}},
    qdrift_e_all::Vector{<:Real}, compton_bands::Vector{<:Unitful.Energy{<:Real}},
    peak::Real = 0.0, window::Tuple{<:Real, <:Real} = (50.0, 8.0),
    hist_start::Real = -20.0, hist_end::Real = 5.0, bin_width::Real = 0.05;
    compton_window::Unitful.Energy{<:Real} = 20.0u"keV",
    aoe_expression = "aoe", qdrift_expression = "qdrift / e",
    pseudo_prior::NamedTupleDist = NamedTupleDist(empty = true),
    pseudo_prior_all::NamedTupleDist = NamedTupleDist(empty = true),
    fit_func::Symbol = :aoe_two_bck, pol_order::Int = 1
)

Corrects the drift time dependence of the A/E parameter.

Events with a calibrated energy ecal_all in one of the Compton bands [band, band + compton_window] for band in compton_bands are selected, and those with A/E within window around peak are used. The corrected value is aoe + Σ_i fct[i] * qdrift_e^i for i = 1:pol_order. The coefficients fct minimize the Gaussian width of the A/E peak, found with a bounded derivative-free optimization (NLopt LN_BOBYQA) limited by the time limit of set_timelimit. The corrected A/E is finally shifted and scaled to mean 0 and standard deviation 1.

Arguments

• aoe_all: A/E values of all events.
• ecal_all: Calibrated energies of all events, with an energy unit.
• qdrift_e_all: Drift-time parameter divided by the energy of all events.
• compton_bands: Lower edges of the energy bands that are used, with an energy unit.
• peak: Center of the A/E peak.
• window: Widths (lower, upper) of the A/E range around peak that is used.
• hist_start, hist_end, bin_width: Range and bin width of the A/E histograms.

Keywords

• compton_window::Unitful.Energy = 20.0u"keV": Width of each Compton band.
• aoe_expression, qdrift_expression: Names of the A/E and drift-time variables in the returned func.
• pseudo_prior: Replaces fields of the default prior of the A/E peak fits.
• pseudo_prior_all: Accepted but not used.
• fit_func::Symbol = :aoe_two_bck: Shape of the A/E peak fits (see fit_single_aoe_compton).
• pol_order::Int = 1: Order of the polynomial in qdrift_e.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• peak, window: the input values.
• func: the correction and normalization as a Julia expression string in aoe_expression and qdrift_expression.
• fct: the coefficients of the correction.
• σ_start, σ_optimal: the optimized Gaussian width of the A/E peak without and with the correction.
• σ_before, σ_after, σ_after_norm: the widths of the A/E peak before the correction, after it, and after the final normalization.
• before, after, after_norm: the corresponding fit results.
• converged: whether the optimizer converged.

report contains peak, window, fct, bin_width, bin_width_qdrift, the selected events (aoe_peak, aoe_ctc, aoe_ctc_norm, qdrift_peak), the histograms h_before, h_after, h_after_norm, the widths and the reports of the three fits, for plotting with LegendMakie.lplot.

ctc_energy(
    e::Vector{<:Unitful.Energy{<:Real}}, qdrift::Vector{<:Real}, peak::Unitful.Energy{<:Real},
    window::Tuple{<:Unitful.Energy{<:Real}, <:Unitful.Energy{<:Real}},
    m_cal_simple::Unitful.Energy{<:Real} = 1.0u"keV";
    e_expression::Union{Symbol, String} = "e", pol_order::Int = 1
)

Corrects the dependence of the energy on the drift time (charge-trapping correction) in the peak at peak.

The events within window = (lower, upper) around peak are used. The corrected energy is e + Σ_i fct[i] * qdrift^i for i = 1:pol_order. The coefficients fct minimize the logarithm of the ratio of the FWHM to the height of the fitted peak, found with a bounded derivative-free optimization (NLopt LN_BOBYQA) limited by the time limit of set_timelimit. The peak is fitted with fit_single_peak_th228.

Arguments

• e: Calibrated energies.
• qdrift: Drift-time parameter of each event.
• peak: Peak energy.
• window: Widths (lower, upper) of the energy window around peak.
• m_cal_simple: Simple calibration constant (energy per uncalibrated unit). fct in the result is divided by it, so that the correction func applies to the uncalibrated energy.

Keywords

• e_expression: Name of the energy variable in the returned func.
• pol_order::Int = 1: Order of the polynomial in qdrift.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• peak, window, m_cal_simple: the input values.
• func: the correction as a Julia expression string, e_expression + fct[1] * qdrift^1 + ....
• fct: the coefficients of the correction (divided by m_cal_simple).
• fwhm_before, fwhm_after: FWHM of the peak before and after the correction, in keV.
• converged: whether the optimization converged.

report contains peak, window, fct, bin_width, bin_width_qdrift, the selected events (e_peak, e_ctc, qdrift_peak), the histograms h_before and h_after, the FWHM values and the reports of both fits (report_before, report_after), for plotting with LegendMakie.lplot.

cut_single_peak(
    x::AbstractVector{<:Unitful.RealOrRealQuantity}, min_x::T, max_x::T;
    n_bins::Int = -1, relative_cut::Float64 = 0.5, n_tries::Int = 5
) where T<:Unitful.RealOrRealQuantity

Determines a window around the single peak in the values x that lie strictly between min_x and max_x.

The values are histogrammed and the window is the range around the highest bin where the counts exceed relative_cut times the maximum count. If the lower edge of the window does not lie below the peak maximum, the histogram is rebuilt with i-times finer bins for i = 2, ..., n_tries, stopping at the first valid window. If none is found, a warning is logged and the window of the last histogram is returned. The output has the unit of min_x.

Keywords

• n_bins::Int = -1: Number of histogram bins. A negative value selects the Freedman-Diaconis bin width.
• relative_cut::Float64 = 0.5: Fraction of the maximum bin content that defines the window edges.
• n_tries::Int = 5: Maximum number of histograms.

Returns

A NamedTuple (low, high, max) with the left bin edge where the counts first reach the cut level on the rising side, the left bin edge where they fall to the cut level on the falling side, and the left edge of the maximum bin.

double_gaussian(
    x::Real, μ1::Real, σ1::Real, n1::Real, 
    μ2::Real, σ2::Real, n2::Real
)

Evaluates the sum of two Gaussian densities at x, with means μ1, μ2, standard deviations σ1, σ2 and amplitudes n1, n2.

estimate_combined_peak_stats(peakstats::StructArray; calib_type::Symbol = :th228)

Estimates the starting values for a combined fit of several peaks from the single-peak statistics peakstats. Only calib_type = :th228 is accepted; other values throw an error.

estimate_fwhm(v::NamedTuple)

Calculates the FWHM of a peak from the fit parameters v.

If v has the field skew_fraction, the FWHM is determined numerically from the Gaussian signal plus low-energy tail. Otherwise it is 2√(2 ln 2) σ. Returns NaN if the half-maximum positions cannot be found or v has no field σ.

estimate_fwhm_aoe_compton(v::NamedTuple, f_fit::Function)

Calculates the FWHM of the peak described by the fit function f_fit(x) with the parameters v. The maximum is searched on a 0.001 grid between v.μ - v.σ and v.μ + v.σ. Returns NaN if the half-maximum positions cannot be found.

estimate_single_peak_stats(h::Histogram; calib_type::Symbol = :th228)
estimate_single_peak_stats(e::Vector{<:Real}, bin_width::Real; calib_type::Symbol = :th228)

Estimates the position, width, size and background level of a single peak in the histogram h, which is used as a starting point for peak fits.

h must contain a single peak with a Gaussian-like shape. The second method, available for calib_type = :th228 only, takes the unbinned values e, builds a histogram with the given bin_width and halves the bin width until the estimated FWHM is non-zero.

Keywords

• calib_type::Symbol = :th228: Estimation variant, one of :th228 (peak on a background with a step, as in the 228Th calibration spectrum), :psd (peak in a pulse shape discrimination parameter, e.g. A/E) or :simple (peak position, width and counts only). Any other value throws an error.

Returns

A NamedTuple with the fields

• peak_pos: estimated position of the peak, from the maximum and the midpoint of the half-maximum range.
• peak_fwhm: full width at half maximum (FWHM) of the peak.
• peak_sigma: estimated standard deviation of the peak.
• peak_counts: estimated number of counts in the peak.
• mean_background: estimated mean background per unit of the horizontal axis (:th228 and :psd only; for :psd it is the number of counts below three standard deviations of the peak).
• bin_width, mean_background_step, mean_background_std: bin width of h, estimated height of the background step across the peak and standard deviation of the background per unit of the horizontal axis (:th228 only).

If the half-maximum range is empty, the width is estimated from the quarter-maximum range, and a fallback of peak_sigma = 1 is used if that also vanishes (:th228 and :simple).

evaluate_energy_point(
    e::AbstractVector{<:Real}, qdrift::AbstractVector{<:Real},
    peak::Unitful.Energy, window::Tuple, min_e::Real, max_e::Real, rel_cut_fit::Real;
    n_bins::Int = -1, apply_ctc::Bool = true, containment::Real = 0.9
)

Evaluates the shape of a gamma peak for one filter setting.

Events with finite, positive e and qdrift are kept. The peak mode is found with cut_single_peak on (min_e, max_e), and the energies are calibrated to keV with m_cal_simple = peak / mode. Events within 15 half-widths of the mode are kept for the fit. With apply_ctc = true, ctc_energy fits the linear charge-trapping correction, and the metrics are computed from the peak fit after the correction. Otherwise the uncorrected peak is fitted with fit_single_peak_th228. The metrics come from peak_width_metrics, which includes the low-energy tail of the peak.

Arguments

• e: Uncalibrated energies (ADC) of the events.
• qdrift: Drift-time parameter of every event, aligned with e.
• peak::Unitful.Energy: Energy of the gamma line.
• window::Tuple: Energy range (left, right) around peak for the correction and its fit.
• min_e, max_e: Range of e (ADC) that brackets the peak.
• rel_cut_fit: relative_cut of cut_single_peak.

Keywords

• n_bins::Int = -1: Number of histogram bins for cut_single_peak; a negative value selects it automatically.
• apply_ctc::Bool = true: Apply the charge-trapping correction before the peak fit.
• containment::Real = 0.9: Fraction contained by w90, see peak_width_metrics.

Returns

A NamedTuple with the fields

• fwhm_before: FWHM in keV of the uncorrected peak fit (equal to fwhm for apply_ctc = false).
• fwhm, w90, tail_scale, centroid: metrics of peak_width_metrics in keV, taken from the fit after the correction and scaled by peak / μ with the fitted peak position μ.
• asym, tail_fraction: unitless metrics of peak_width_metrics.
• fct: first-order coefficient of the correction, in ADC per qdrift unit (0 for apply_ctc = false).
• ct_spread: FWHM-equivalent spread in keV of the drift-time term that the correction removes, 2 * sqrt(2 * log(2)) * |fct| * m_cal_simple * std(qdrift) scaled by peak / μ. std(qdrift) is taken over the events inside the correction window (0 for apply_ctc = false).
• m_cal_simple: calibration factor in keV per ADC.
• ctc_converged: convergence flag of the correction (true for apply_ctc = false).
• n_events: number of events within 15 half-widths of the mode.

Throws an ErrorException if fewer than 100 events lie in (min_e, max_e), if the peak-statistics guess has a non-positive entry, or if the peak fit returns a non-finite or non-positive σ.

ex_gauss_pdf(x::Real, μ::Real, σ::Real, θ::Real)

Evaluates the probability density function of an exponentially modified Gaussian distribution with Gaussian mean μ and standard deviation σ and exponential scale θ at x.

It is the density of the random variable rand(Normal(μ, σ)) + rand(Exponential(θ)). The function throws an AssertionError if the result is not finite.

ex_step_gauss(x::Real, l::Real, k::Real, t::Real, d::Real)

Evaluates (exp(k * (x - l)) + d) / (exp((x - l) / t) + l) at x.

This is the energy dependence of the A/E resolution model of GERDA (Eq. (1) in EPJC 73 (2013) 2583).

expand_vars(v::NamedTuple)::StructArray

Broadcasts all fields of v (scalars or arrays) to a common array size and returns them as a StructArray.

Scalar fields are repeated lazily (without copying) along the common size.

exponential_decay(x::Real, amplitude::Real, decay::Real, offset::Real)

Evaluates amplitude * exp(-decay * x) + offset at x.

fit_aoe_compton(
    peakhists::Vector{<:Histogram}, peakstats::StructArray, compton_bands::Array{T};
    uncertainty::Bool = false, fit_func::Symbol = :aoe_one_bck
) where T<:Unitful.Energy{<:Real}

Fits the A/E distribution in each Compton band with fit_single_aoe_compton. The fits run in parallel on the available threads.

Arguments

• peakhists: A/E histogram of each band, e.g. from generate_aoe_compton_bands.
• peakstats: Initial peak statistics of each band.
• compton_bands: Energies of the bands; they are used as keys of the result.

Keywords

• uncertainty::Bool = false: Calculate parameter uncertainties and goodness of fit.
• fit_func::Symbol = :aoe_one_bck: :aoe_one_bck or :aoe_two_bck (see get_aoe_fit_functions).

Returns

A tuple (result, report) of two OrderedDicts keyed by compton_bands, with the result and report NamedTuples of fit_single_aoe_compton as values.

fit_aoe_compton_combined(
    peakhists::Vector{<:Histogram}, peakstats::StructArray, compton_bands::Array{T},
    result_corrections::NamedTuple;
    fit_func::Symbol = :aoe_one_bck, aoe_expression = "a / e", e_expression = "e",
    uncertainty::Bool = false
) where T<:Unitful.Energy{<:Real}

Fits the A/E distributions of all Compton bands simultaneously, with the peak position and width given by smooth functions of the band energy E: μ(E) = μA + μB E and σ(E) = √(σA² + σB² / E²). The four parameters μA, μB, σA, σB are fitted by minimizing the summed negative log-likelihood of all bands (NLopt LN_BOBYQA), where the other parameters of each band are fitted for every trial.

Arguments

• peakhists: A/E histogram of each band, e.g. from generate_aoe_compton_bands.
• peakstats: Initial peak statistics of each band.
• compton_bands: Energies of the bands.
• result_corrections: Result of fit_aoe_corrections, which gives the start values and prior widths of the four parameters.

Keywords

• fit_func::Symbol = :aoe_one_bck: :aoe_one_bck or :aoe_two_bck (see get_aoe_fit_functions).
• aoe_expression, e_expression: Names of the A/E and energy variables in the returned func.
• uncertainty::Bool = false: Calculate uncertainties from the Hessian of the summed likelihood, and the goodness of fit.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• μA, μB, σA, σB: the fitted parameters (Measurements if uncertainty is true).
• gof: converged and, if the uncertainties were calculated, pvalue, chi2, dof (summed over the converged band fits), covmat and the mean, median and standard deviation of the normalized residuals.
• μ_compton, σ_compton: NamedTuples with par and func in the format of the result of fit_aoe_corrections.
• compton_bands: NamedTuple with the unitless band energies e.
• func: the normalized A/E as a Julia expression string.

report contains the best-fit values v, report_μ and report_σ (values of μ(E) and σ(E) at the band energies with plot labels) and report_bands, an OrderedDict of the reports of the individual band fits keyed by band energy, for plotting with LegendMakie.lplot.

fit_aoe_corrections(
    e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Real}, σ::Array{<:Real};
    aoe_expression::Union{String, Symbol} = "a / e", e_expression::Union{String, Symbol} = "e"
)

Fits the energy dependence of the mean μ and standard deviation σ of the A/E distribution in energy bands (Compton bands) centered at e, to obtain a correction that normalizes A/E to mean 0 and standard deviation 1.

μ and σ are vectors of Measurements. Points whose uncertainty is zero are dropped. The mean is fitted with a linear function μ(E) = p1 + p2 E and the standard deviation with σ(E) = √(p1² + p2² / E²), both with chi2fit.

Keywords

• aoe_expression: Name of the A/E variable in the returned func.
• e_expression: Name of the energy variable in the returned expressions.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• μ_compton: chi2fit result of the mean fit, extended by par (with units), func (expression string) and μ.
• σ_compton: the same for the standard deviation fit, with σ instead of μ.
• compton_bands: NamedTuple with the energy values e (unitless, in the unit of the input).
• func: the normalized A/E, (a - μ(E)) / σ(E), as a Julia expression string.

report contains report_μ and report_σ, for plotting with LegendMakie.lplot.

fit_binned_double_gauss(h::Histogram, ps::NamedTuple; uncertainty::Bool = true)

Fits the sum of two Gaussians (double_gaussian) to the histogram h by maximizing the binned Poisson likelihood.

The priors of the first Gaussian are centered on the peak statistics ps (peak_pos, peak_sigma, peak_counts, e.g. from estimate_single_peak_stats); the second Gaussian has a wide prior with a mean between 0 and 1200.

Keywords

• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit.

Returns

A tuple (result, report). result is a NamedTuple with the parameters μ1, σ1, n1, μ2, σ2, n2 (Measurements if uncertainty is true and the fit converged) and, in that case, gof. report contains the best-fit values v, the histogram h, the total fit function f_fit and the two components f_gauss_1 and f_gauss_2.

fit_binned_half_trunc_gauss(
    h_nocut::Histogram, cuts::NamedTuple{(:low, :high, :max)} = (low = NaN, high = NaN, max = NaN);
    left::Bool = false, uncertainty::Bool = true
)

Fits one half of a Gaussian with amplitude n to the histogram h_nocut by maximizing the binned Poisson likelihood. Mean and standard deviation are fitted.

The right half (between cuts.max and cuts.high) is used by default; with left = true, the left half (between cuts.low and cuts.max). The units of the cuts must agree.

Keywords

• left::Bool = false: Fit the left half instead of the right half.
• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit.

Returns

A NamedTuple (result, report) with the same structure as fit_binned_trunc_gauss.

fit_binned_trunc_gauss(
    h_nocut::Histogram, cuts::NamedTuple{(:low, :high, :max)} = (low = NaN, high = NaN, max = NaN);
    uncertainty::Bool = true
)

Fits a Gaussian with amplitude n, truncated to [cuts.low, cuts.high], to the histogram h_nocut by maximizing the binned Poisson likelihood. A NaN bound is replaced by the edge of the histogram. The units of the cuts must agree.

Keywords

• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit.

Returns

A NamedTuple (result, report). result has the fields μ, σ and n (number of counts; Measurements if uncertainty is true and the fit converged) and, in that case, gof with pvalue, chi2, dof, covmat, residuals, residuals_norm and bin_centers. report contains the fit function f_fit (counts per bin), the histogram h, μ, σ and gof, for plotting with LegendMakie.lplot.

fit_calibration(
    pol_order::Int, µ::AbstractVector{<:Union{Unitful.RealOrRealQuantity, Measurement{<:Unitful.RealOrRealQuantity}}},
    peaks::AbstractVector{<:Quantity};
    e_expression::Union{Symbol, String} = "e", m_cal_simple = nothing, uncertainty::Bool = true
)

Fits a polynomial calibration curve that maps the uncalibrated peak positions µ to the known peak energies peaks, using chi2fit.

pol_order = 1 fits a linear function and pol_order = 2 a quadratic function. The peak energies peaks define the energy unit of the result. Uncertainties of µ enter the fit through error propagation. If µ carries a unit, it is converted to that energy unit and stripped, with a warning.

Keywords

• e_expression::Union{Symbol, String} = "e": Name of the uncalibrated energy variable in func and func_err.
• m_cal_simple = nothing: Slope of a simple calibration used as the starting point and prior center of the fit. If nothing, the ratio of the highest peak energy to its position µ is used.
• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• par: best-fit coefficients with energy units; entry i multiplies e^(i-1).
• converged, gof: convergence flag and goodness of fit as returned by chi2fit.
• func: the calibration curve as a Julia expression string in e_expression.
• func_err: the same expression with the uncertainty of the curve written out as value ± error, including the correlations between the fit parameters.
• µ, peaks: the input peak positions and energies.

report is the report of chi2fit extended by e_unit, par and type = :cal, for plotting with LegendMakie.lplot.

fit_enc_sigmas(
    enc_grid::Matrix{T}, enc_grid_rt::StepRangeLen,
    min_enc::T, max_enc::T, nbins::Int, rel_cut_fit::T
) where T<:Real

Finds the rise time (RT) of a trapezoidal filter that minimizes the equivalent noise charge (ENC).

For each RT, the ENC values of that row of enc_grid are fitted with a truncated Gaussian (see fit_single_trunc_gauss) after restricting them to a window found with cut_single_peak. The Gaussian standard deviation is the ENC for that RT. Rows that are all zero or all non-finite are skipped. An error is thrown if no RT yields a valid fit.

Arguments

• enc_grid: ENC values; row i holds the values for enc_grid_rt[i].
• enc_grid_rt: Range of RT values (with time units) of the rows of enc_grid.
• min_enc, max_enc: Range of ENC values considered in the fit.
• nbins: Number of histogram bins for each fit.
• rel_cut_fit: Relative cut level for cut_single_peak.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• rt: RT with the smallest ENC, with the step of enc_grid_rt as uncertainty.
• min_enc: corresponding ENC (a Measurement).

report contains rt, min_enc, the RT values enc_grid_rt and the ENC values enc of all successful fits, for plotting with LegendMakie.lplot.

fit_fwhm(
    pol_order::Int, peaks::Vector{<:Unitful.Energy{<:Real}}, fwhm::Vector{<:Unitful.Energy{<:Real}};
    e_type_cal::Symbol = :e_cal, e_expression::Union{Symbol, String} = "e",
    uncertainty::Bool = true, correlated::Bool = true, pull_fano::Bool = false
)
fit_fwhm(peaks, fwhm; kwargs...)

Fits the energy resolution curve fwhm(E) = √(enc + fano·E [+ ct·E²]) to the full widths at half maximum fwhm of the peaks at energies peaks, using chi2fit. The polynomial under the square root has degree pol_order (1 or 2). The method without pol_order uses pol_order = 1.

The energies are fitted in MeV and the squared widths in keV². The Fano factor is derived from the linear coefficient with a pair creation energy of 2.96 eV for germanium. The method evaluates the FWHM at Qββ = 2039.061 keV (the qbb field).

Keywords

• e_type_cal::Symbol = :e_cal: Name of the calibrated energy variable in func_cal and func_cal_err.
• e_expression::Union{Symbol, String} = "e": Name of the energy variable in func and func_err.
• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit.
• correlated::Bool = true: Return parameters that carry their full covariance (see chi2fit).
• pull_fano::Bool = false: For pol_order = 2, add a Gaussian pull term that ties the linear coefficient to the Fano statistics of germanium (Fano factor 0.112 ± 0.015), so that the E² term absorbs the other contributions. Without it, and for pol_order = 1, the linear term is free.

Returns

A tuple (result, report). result is a NamedTuple with the fields of the chi2fit result (converged, gof) and

• par: enc, the FWHM √enc in keV that the curve approaches at zero energy; fano, the dimensionless Fano factor; and, for pol_order = 2, the dimensionless coefficient ct of E².
• qbb: FWHM at 2039.061 keV, in keV.
• func, func_cal: the resolution curve as Julia expression strings in e_expression and e_type_cal.
• func_err, func_cal_err: the same expressions with the uncertainty of the curve written out as value ± error, including the correlations between the fit parameters. Evaluating them at an energy that is itself a Measurement nests the two uncertainties.
• peaks, fwhm: the input values.

report is the report of chi2fit with f_fit replaced by the fitted FWHM curve in keV, extended by e_unit, par, qbb and type = :fwhm, for plotting with LegendMakie.lplot.

fit_fwhm_ft(
    e_grid::Matrix, e_grid_ft::StepRangeLen, qdrift::Vector{<:Real}, rt,
    min_e::Real, max_e::Real, rel_cut_fit::Real, apply_ctc::Bool = true;
    n_bins::Int = -1, peak = 2614.5u"keV", window = (35.0u"keV", 25.0u"keV"),
    objective::Symbol = :fwhm, ctc_rel_uncertainty = nothing, containment::Real = 0.9, rel_tol::Real = 0.005
)
fit_fwhm_ft(
    e_grid::Matrix, e_grid_ft::StepRangeLen, rt, min_e, max_e, rel_cut_fit;
    n_bins::Int = -1, peak = 2614.5u"keV", window = (35.0u"keV", 25.0u"keV"), ft_fwhm_tol = 0.1u"keV"
)

Determines the flat-top time (FT) of a trapezoidal filter that gives the best peak shape of a gamma line.

Row i of e_grid holds the energies for the flat-top time e_grid_ft[i].

With qdrift, each row is evaluated in parallel with evaluate_energy_point. The flat-top time with the smallest objective (see objective_value) is selected; flat-top times within rel_tol of the smallest value tie, and the longest one wins (select_optimum). A row whose evaluation throws is kept in the report with its error text in status and NaN metrics, and a warning is logged.

Without qdrift, each row is restricted to a window around the peak and fitted without correction with fit_single_peak_th228. The FWHM values are converted to keV with one factor that places the mean peak position of all rows at peak. Among the flat-top times with an FWHM within ft_fwhm_tol of the smallest one, the shortest wins. Rows with fewer than 100 energies in (min_e, max_e) or an invalid fit start are skipped.

Arguments

• e_grid: Energies (ADC), one row per flat-top time.
• e_grid_ft::StepRangeLen: Flat-top times of the rows, with time units.
• qdrift::Vector{<:Real}: Drift-time parameter of every event, for the charge-trapping correction.
• rt: Rise time of the scan. It is not used in the computation.
• min_e, max_e: Range of e_grid (ADC) that brackets the peak.
• rel_cut_fit: relative_cut of cut_single_peak.
• apply_ctc::Bool = true: Apply the charge-trapping correction before the peak fit (with qdrift only).

Keywords

• n_bins::Int = -1: Number of histogram bins for cut_single_peak; a negative value selects it automatically.
• peak = 2614.5u"keV": Energy of the gamma line.
• window = (35.0u"keV", 25.0u"keV"): Energy range (left, right) around peak for the correction and its fit (with qdrift only).
• objective::Symbol = :fwhm: Objective of objective_value (with qdrift only).
• ctc_rel_uncertainty = nothing: Required for the objective :fwhm_ct, which also needs apply_ctc = true.
• containment::Real = 0.9: Fraction contained by w90, see peak_width_metrics.
• rel_tol::Real = 0.005: Relative tolerance of the tie.
• ft_fwhm_tol = 0.1u"keV": Absolute FWHM tolerance of the tie (without qdrift only).

Returns

A tuple (result, report). With qdrift, result is a NamedTuple with the fields

• ft: selected flat-top time, with step(e_grid_ft) as uncertainty.
• min_fwhm: FWHM in keV of the selected flat-top time.
• min_objective: objective value of the selected flat-top time.
• objective: the objective used.

report, for plotting with LegendMakie.lplot, contains the fields of result and

• table: TypedTables.Table with one row per flat-top time: the columns ft and status ("ok" or the error text), followed by the fields of evaluate_energy_point.
• objective_values: objective for every row of table.

Without qdrift, result has the fields ft and min_fwhm, and report contains ft, min_fwhm, the flat-top times e_grid_ft and the FWHM values fwhm of all successful fits.

Throws an ErrorException if no flat-top time is evaluated successfully. If the rows of e_grid do not match e_grid_ft, the method with qdrift throws a DimensionMismatch and the method without qdrift an AssertionError. The method with qdrift throws an ArgumentError for an invalid objective or ctc_rel_uncertainty, or for :fwhm_ct with apply_ctc = false.

fit_half_centered_trunc_gauss(
    x::Vector{<:Unitful.RealOrRealQuantity}, μ::T,
    cuts::NamedTuple{(:low, :high, :max), Tuple{T, T, T}};
    left::Bool = false, uncertainty::Bool = true
) where T<:Unitful.RealOrRealQuantity

Fits one half of a Gaussian with fixed mean μ to the values x by unbinned maximum likelihood. Only the standard deviation is fitted.

The right half (x >= μ, up to cuts.high) is used by default; with left = true, the left half (x <= μ, down to cuts.low). cuts can be obtained from cut_single_peak.

Keywords

• left::Bool = false: Fit the left half instead of the right half.
• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit.

Returns

A NamedTuple (result, report) with the same structure as fit_single_trunc_gauss.

fit_half_trunc_gauss(
    x::Vector{<:Unitful.RealOrRealQuantity}, cuts::NamedTuple{(:low, :high, :max), Tuple{T, T, T}};
    left::Bool = false, uncertainty::Bool = true
) where T<:Unitful.RealOrRealQuantity

Fits one half of a Gaussian to the values x by unbinned maximum likelihood. Mean and standard deviation are fitted.

The right half is used by default, i.e. the values between cuts.max (the peak maximum) and cuts.high; with left = true, the left half between cuts.low and cuts.max. cuts can be obtained from cut_single_peak.

Keywords

• left::Bool = false: Fit the left half instead of the right half.
• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit.

Returns

A NamedTuple (result, report) with the same structure as fit_single_trunc_gauss.

fit_lq_compton(
    peakhists::Vector{<:Histogram}, peakstats::StructArray, compton_bands::Array{T};
    pars_lq::NamedTuple{(:μ, :μ_err, :σ, :σ_err)} = (μ = nothing, μ_err = nothing, σ = nothing, σ_err = nothing),
    uncertainty::Bool = false
) where T<:Unitful.Energy{<:Real}

Fits the LQ distribution in each Compton band with fit_single_lq_compton.

The model is a Gaussian peak plus an exponentially modified Gaussian background (aoe_compton_peakshape). If pars_lq.μ is given, μ and σ in each band's peak statistics are replaced by the energy-dependent functions of the parameters pars_lq.μ (linear) and pars_lq.σ (exponential decay). Bands whose fit throws an error are skipped with a warning.

Arguments

• peakhists: LQ histogram of each band.
• peakstats: Initial peak statistics of each band.
• compton_bands: Energies of the bands; they are used as keys of the result.

Returns

A tuple (result, report) of two OrderedDicts keyed by compton_bands, with the result and report NamedTuples of fit_single_lq_compton as values.

fit_peaks(
    peakhists::Array, peakstats::StructArray, th228_lines::Vector;
    calib_type::Symbol, kwargs...
)

Fits a gamma peak shape to each histogram in peakhists, one per calibration line in th228_lines, using fit_single_peak_th228 with the initial values in peakstats. The fits run in parallel on the available threads.

Arguments

• peakhists: One histogram per line.
• peakstats: StructArray of initial peak statistics, one entry per line, as returned by estimate_single_peak_stats (:th228 variant).
• th228_lines: Calibration line energies; they are used as keys of the result.

Keywords

• calib_type::Symbol: Required. Only :th228 is supported; other values throw an error.
• e_unit::Union{Nothing, Unitful.EnergyUnits} = nothing: Energy unit attached to the fitted μ, σ, fwhm and centroid. Without it, they are unitless.
• uncertainty::Bool = true: Calculate parameter uncertainties and goodness of fit.
• low_e_tail::Bool = true: Fit the low-energy tail.
• iterative_fit::Bool = false: If true and the covariance matrix of a fit is not positive definite, repeat that fit without low-energy tail.
• fit_func::Vector{Symbol}: Fit function of each line (see fit_single_peak_th228). Default is :gamma_def for all lines.
• pseudo_prior::NamedTupleDist = NamedTupleDist(empty = true): Replaces fields of the default prior of every fit (see fit_single_peak_th228).
• m_cal_simple = 1.0: Simple calibration constant; μ, σ, fwhm and centroid are divided by it, which converts them back to the scale before the simple calibration.

Returns

A tuple (result, report) of two OrderedDicts keyed by the line energies in th228_lines. The values of result are the result NamedTuples of fit_single_peak_th228 and the values of report are its report NamedTuples, for plotting with LegendMakie.lplot.

fit_peaks_combined(
    peakhists::Array, peakstats::StructArray, th228_lines::Array;
    calib_type::Symbol = :th228, uncertainty::Bool = true, fixed_position::Bool = false
)

Fits several peaks simultaneously by maximum likelihood with parameters shared between the peaks.

The peak position is a linear function of the line energy e, μ = μ_slope * e + μ_intercept, and the low-energy tail fraction is a linear function of e, skew_fraction = skew_fraction_slope * e + skew_fraction_intercept. The width σ, amplitude n, step_amplitude, skew_width and background are fitted separately for each peak. Priors are built from peakstats.

Arguments

• peakhists: One histogram per line.
• peakstats: StructArray of initial peak statistics, one entry per line.
• th228_lines: Line energies e (as real numbers).

Keywords

• calib_type::Symbol = :th228: Only :th228 is supported; other values throw an error.
• uncertainty::Bool = true: Accepted for interface compatibility; no uncertainties are calculated.
• fixed_position::Bool = false: Fix the peak positions to the positions in peakstats.

Returns

A tuple (v_ml, opt_r) with the best-fit parameters as a NamedTuple and the optimizer result.

fit_pmt_spe(x::AbstractVector{<:Real}; kwargs...)

Fits the single-photoelectron (1 p.e.) pulse-height spectrum of a photomultiplier tube (PMT) with an exponential and the 1 p.e. and 2 p.e. peaks.

The model, in counts per unit of x, is

f(x) = n_exp / τ * exp(-(x - x_min) / τ) + n_1 * N(x; μ1, σ1) + r2 * n_1 * N(x; 2 μ1, √2 σ1)

where N is a normalized Gaussian. n_exp and n_1 are the numbers of entries of the exponential (above x_min) and of the 1 p.e. peak, and r2 is the ratio of 2 p.e. to 1 p.e. entries, bounded by 1 so that the fit cannot place the 1 p.e. peak at half its position.

The pedestal is a sharp, non-Gaussian peak at the low end and is not fitted. The fit starts at the valley between the pedestal and the 1 p.e. peak. The histogram of x in hist_range with bins of bin_width is smoothed with a 3-bin moving average. Starting at first(peak_search_range), the valley is the first bin after which the smoothed spectrum rises for two bins. The start value of μ1 is the maximum of the smoothed spectrum between the valley and last(peak_search_range). The fit range starts at x_min = valley - x_min_offset (at least first(hist_range)), and the binned Poisson likelihood of the histogram above x_min is maximized, with L-BFGS-B and, if that does not converge, with BOBYQA. The bounds of the parameters are μ1 in peak_search_range, σ1 between bin_width / 2 and the width of peak_search_range, τ between bin_width / 2 and the fit range, r2 in [0, 1], and positive n_exp and n_1.

Arguments

• x::AbstractVector{<:Real}: Pulse heights. Non-finite values are ignored.

Keywords

• hist_range = (0.0, 100.0): Range of the histogram, two numbers (a tuple or a vector) in the unit of x. It must be a whole number of bins.
• bin_width::Real = 1.0: Bin width.
• peak_search_range = (6.0, 60.0): Range in which the valley and the 1 p.e. peak are searched, two numbers within hist_range.
• x_min_offset::Real = 2.0: Distance of the fit start below the valley.
• σ1_rel_init::Real = 0.33: Start value of σ1 relative to the start value of μ1.
• τ_init::Real = 5.0: Start value of τ.
• r2_init::Real = 0.1: Start value of r2.
• min_counts::Real = 50: Minimum number of entries above x_min.
• uncertainty::Bool = true: Compute the parameter uncertainties and the goodness of fit.

The defaults are for pulse heights in ADC units with the 1 p.e. peak at about 10 to 20 units.

An ArgumentError is thrown if a range does not have two elements or peak_search_range does not lie within hist_range, if hist_range is not a whole number of bins, if peak_search_range covers fewer than three bins, if the smoothed spectrum does not rise again in peak_search_range (no 1 p.e. peak above the pedestal), or if there are fewer than min_counts entries above x_min. A warning is emitted if the fit does not converge.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• μ1, σ1: position and width of the 1 p.e. peak.
• τ: decay constant of the exponential.
• n_exp, n_1: numbers of entries of the exponential and of the 1 p.e. peak.
• r2: ratio of 2 p.e. to 1 p.e. entries.
• x_min: start of the fit range.
• gof: goodness of fit. With uncertainty, it has converged, pvalue, chi2, dof, covmat, residuals, residuals_norm and bin_centers; without, only converged.

With uncertainty, the parameters are Measurements. Their uncertainties are NaN if the fit did not converge or the Hessian matrix is not finite.

report contains the histogram h of all of hist_range, the model f_fit (NaN below x_min), f_components (a NamedTuple with the functions exponential, pe1 and pe2), x_min, μ1 and gof, for plotting with LegendMakie.lplot.

fit_sf_wl(
    e_dep::Vector{<:Real}, aoe_dep::ArrayOfSimilarArrays{<:Real},
    e_sep::Vector{<:Real}, aoe_sep::ArrayOfSimilarArrays{<:Real}, a_grid_wl_sg::StepRangeLen;
    dep = 1592.53u"keV", dep_window = [12.0, 10.0]u"keV",
    sep = 2103.53u"keV", sep_window = [25.0, 25.0]u"keV", sep_rel_cut::Real = 0.5,
    min_aoe_quantile::Real = 0.1, max_aoe_quantile::Real = 0.99,
    min_aoe_offset = 0.0u"keV^-1", max_aoe_offset = 0.05u"keV^-1",
    dep_cut_search_fit_func::Symbol = :gamma_def, sep_cut_search_fit_func::Symbol = :gamma_def,
    uncertainty::Bool = false
)

Finds the window length of the A/E filter that minimizes the survival fraction of the single escape peak (SEP) when the A/E cut keeps a fixed fraction of the double escape peak (DEP).

For each window length in a_grid_wl_sg, the A/E cut is determined on the DEP with get_low_aoe_cut and applied to the SEP, whose survival fraction is determined with get_peak_survival_fraction. Only survival fractions between 1 % and 100 % are considered. The energies are calibrated with a simple calibration constant from the SEP fit (fit_single_peak_th228). An error is thrown if no window length gives a valid survival fraction.

Arguments

• e_dep, e_sep: Uncalibrated energies of the DEP and SEP events.
• aoe_dep, aoe_sep: A/E values of those events for each window length; element i of the flattened array corresponds to a_grid_wl_sg[i].
• a_grid_wl_sg: Range of window lengths.

Keywords

• dep, dep_window, sep, sep_window: Energies of the peaks and lower and upper widths of the fit windows.
• sep_rel_cut: Relative cut level for locating the SEP.
• min_aoe_quantile, max_aoe_quantile: Quantiles of the DEP A/E values that bound the cut search.
• min_aoe_offset, max_aoe_offset: Offsets added to those bounds.
• dep_cut_search_fit_func, sep_cut_search_fit_func: Peak shapes for the DEP and SEP fits (see fit_single_peak_th228).
• uncertainty::Bool = false: Calculate fit uncertainties.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• wl: optimal window length, with the step of a_grid_wl_sg as uncertainty.
• sf: SEP survival fraction at that window length, in percent.
• n_dep, n_sep: number of DEP and SEP events.

report contains wl, min_sf, the window lengths a_grid_wl_sg and the survival fractions sfs of all successful window lengths, for plotting with LegendMakie.lplot.

fit_single_aoe_compton(
    h::Histogram, ps::NamedTuple;
    uncertainty::Bool = true, pseudo_prior::NamedTupleDist = NamedTupleDist(empty = true),
    fit_func::Symbol = :aoe_one_bck, background_center = ps.peak_pos, fixed_position::Bool = false
)

Fits the A/E distribution of the histogram h by maximizing the binned Poisson likelihood, starting from the peak statistics ps.

The model is a Gaussian single-site-event peak plus a background of multi-site events modeled as one (:aoe_one_bck) or two (:aoe_two_bck) exponentially modified Gaussians (see get_aoe_fit_functions). ps is a result of estimate_single_peak_stats (:psd variant).

Keywords

• uncertainty::Bool = true: Calculate parameter uncertainties (from the Hessian of the likelihood) and goodness of fit.
• pseudo_prior::NamedTupleDist: Replaces fields of the default prior.
• fit_func::Symbol = :aoe_one_bck: :aoe_one_bck or :aoe_two_bck.
• background_center: Passed to the peak components in the report.
• fixed_position::Bool = false: Fix the peak position μ to ps.peak_pos.

Returns

A tuple (result, report). result is a NamedTuple with the fit parameters (μ, σ, n, B, δ, and μ2, σ2, B2f, δ2 for :aoe_two_bck) as Measurements, the fwhm of the fitted distribution, fit_func and gof. gof holds converged and, if the uncertainties were calculated, pvalue, chi2, dof (from p_value), covmat and the mean, median and standard deviation of the normalized residuals. report contains the best-fit values v, the histogram h, f_fit, f_components and gof, for plotting with LegendMakie.lplot.

fit_single_aoe_compton_with_fixed_μ_and_σ(
    h::Histogram, μ::Number, σ::Number, ps::NamedTuple;
    fit_func::Symbol = :aoe_one_bck, background_center = μ, uncertainty::Bool = false
)

Fits the A/E distribution of one Compton band like fit_single_aoe_compton, with the peak position μ and width σ fixed to the given values. Only the remaining parameters are fitted and receive uncertainties.

Returns

A tuple (result, report) with the same structure as fit_single_aoe_compton (without fwhm and fit_func).

fit_single_lq_compton(h::Histogram, ps::NamedTuple; uncertainty::Bool = true)

Fits the LQ distribution of the histogram h with aoe_compton_peakshape by maximizing the binned Poisson likelihood. The uniform priors of the parameters are centered on the peak statistics ps (peak_pos, peak_sigma, peak_counts, mean_background).

Returns

A tuple (result, report). result is a NamedTuple with the parameters μ, σ, n, B, δ (Measurements if uncertainty is true and the fit converged, plain numbers otherwise) and, in the former case, gof with pvalue, chi2, dof (from p_value), covmat, residuals, residuals_norm and bin_centers. report contains the best-fit values v, the histogram h and the functions f_fit, f_sig and f_bck of the total, signal and background shapes.

fit_single_peak_th228(
    h::Histogram, ps::NamedTuple;
    uncertainty::Bool = true, low_e_tail::Bool = true, fixed_position::Bool = false,
    pseudo_prior::NamedTupleDist = NamedTupleDist(empty = true),
    fit_func::Symbol = :gamma_def, background_center = ps.peak_pos, m_cal_simple::Real = 1.0
)

Fits a gamma peak shape to the histogram h by maximizing the binned Poisson likelihood, starting from the peak statistics ps, and calculates the FWHM of the fitted shape.

ps is the NamedTuple returned by estimate_single_peak_stats for calib_type = :th228, with the fields peak_pos, peak_fwhm, peak_sigma, peak_counts, bin_width, mean_background, mean_background_step and mean_background_std. It sets the location and width of the default prior of each parameter. The FWHM and its uncertainty are determined by Monte Carlo error propagation with get_peak_fwhm_th228.

Keywords

• uncertainty::Bool = true: Calculate parameter uncertainties (from the Hessian of the likelihood) and goodness of fit. The uncertainties are NaN if false or if the fit did not converge.
• low_e_tail::Bool = true: Fit the low- and high-energy tails. If false, the tail fractions are fixed to zero.
• fixed_position::Bool = false: Fix the peak position μ to ps.peak_pos.
• pseudo_prior::NamedTupleDist = NamedTupleDist(empty = true): Replaces fields of the default prior. Its keys must be parameters of fit_func.
• fit_func::Symbol = :gamma_def: Peak shape, built from gamma_peakshape:
  • :gamma_def: Gaussian signal, low-energy tail, flat background and step background.
  • :gamma_tails: :gamma_def plus a high-energy tail (skew_fraction_highE, skew_width_highE).
  • :gamma_bckSlope: :gamma_def plus a linear background slope (background_slope).
  • :gamma_bckExp: :gamma_def plus an exponential background (background_exp).
  • :gamma_bckFlat: :gamma_def without the step background (no step_amplitude).
  • :gamma_tails_bckFlat: :gamma_tails without the step background.
  • :gamma_minimal: Gaussian signal and flat background only (μ, σ, n, background).
• background_center::Union{Real, Nothing} = ps.peak_pos: Reference position of the slope and exponential background terms.
• m_cal_simple::Real = 1.0: Simple calibration constant; μ, σ, fwhm and centroid are divided by it after the fit.

Returns

A tuple (result, report). result is a NamedTuple with

• the fit parameters of fit_func (μ, σ, n, skew_fraction, skew_width, background, step_amplitude, ...) as Measurements;
• fwhm: FWHM of the fitted signal including the low-energy tail;
• centroid: mean position of the peak including the tails (see peak_centroid);
• fit_func: the peak shape that was used;
• gof: NamedTuple with converged and, if the uncertainties were calculated, pvalue, chi2, dof (Poisson likelihood-ratio based, see p_value_poissonll), covmat, and mean_residuals, median_residuals, std_residuals of the normalized residuals.

report contains the best-fit values v, the histogram h, the fit function f_fit, the peak components f_components and gof (including residuals and residuals_norm), for plotting with LegendMakie.lplot.

fit_single_trunc_gauss(
    x::Vector{<:Unitful.RealOrRealQuantity},
    cuts::NamedTuple{(:low, :high, :max)} = (low = NaN, high = NaN, max = NaN);
    uncertainty::Bool = true
)

Fits a Gaussian, truncated to [cuts.low, cuts.high], to the values x by unbinned maximum likelihood. A NaN bound is replaced by the minimum or maximum of x.

cuts can be obtained from cut_single_peak. Only cuts.low and cuts.high are used. The goodness of fit is determined from a histogram with the Freedman-Diaconis bin width.

Keywords

• uncertainty::Bool = true: Calculate parameter uncertainties (from the Hessian of the likelihood) and goodness of fit.

Returns

A NamedTuple (result, report). result has the fields μ and σ (mean and standard deviation of the Gaussian, in the unit of x; Measurements if uncertainty is true and the fit converged) and, in that case, gof with pvalue, chi2, dof, covmat, residuals, residuals_norm and bin_centers. report contains the fitted density f_fit, a histogram h normalized to a probability density, μ, σ and gof, for plotting with LegendMakie.lplot.

fit_sipm_spectrum(
    pe_cal::Vector{<:Real}, min_pe::Real = 0.5, max_pe::Real = 3.5;
    n_mixtures::Int = ceil(Int, (max_pe - min_pe) * 4), nIter::Int = 25, nInit::Int = 50,
    method::Symbol = :kmeans, kind = :diag, Δpe_peak_assignment::Real = 0.3,
    f_uncal::Function = identity, uncertainty::Bool = true
)

Fits a Gaussian mixture model to the calibrated SiPM amplitudes pe_cal (in photoelectron units) between min_pe and max_pe, and derives the position and resolution of each photoelectron peak.

The mixture is fitted with the expectation-maximization algorithm of GaussianMixtures.jl. The position of the peak at integer PE value pe is the weighted mean of the Gaussian centers within pe ± Δpe_peak_assignment; its resolution is the FWHM of the sum of these Gaussians. Uncertainties are determined from the Hessian of the binned likelihood.

Arguments

• pe_cal: Calibrated amplitudes in PE units.
• min_pe, max_pe: Fit range.

Keywords

• n_mixtures::Int: Number of Gaussians.
• nIter::Int = 25, nInit::Int = 50: Number of EM iterations and of initializations.
• method::Symbol = :kmeans: Initialization method of the mixture.
• kind = :diag: Covariance type of the mixture.
• Δpe_peak_assignment::Real = 0.3: Distance in PE within which a Gaussian is assigned to a peak.
• f_uncal::Function = identity: Function that converts calibrated values to uncalibrated ones, to express positions and resolutions in uncalibrated units.
• uncertainty::Bool = true: Calculate uncertainties and goodness of fit.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• μ, σ, w: means, standard deviations and weights of the Gaussians (Measurements).
• n_mixtures, n_pos_mixtures: number of Gaussians in total and assigned to each peak.
• peaks: integer PE values of the peaks.
• positions_cal, resolutions_cal: peak positions and FWHM in PE units.
• positions, resolutions: the same in uncalibrated units (via f_uncal).
• gof: pvalue, chi2, dof (from p_value_poissonll), covmat and mean, median and standard deviation of the normalized residuals (empty if uncertainty is false).

report contains the histogram h_cal, the fit function f_fit, the components f_fit_components(x, i), min_pe, max_pe, bin_width, n_mixtures, n_pos_mixtures, peaks, positions, μ and gof, for plotting with LegendMakie.lplot.

fit_sipm_threshold(
    thresholds::Vector{<:Real}, min_cut::Real = minimum(thresholds), max_cut::Real = maximum(thresholds);
    n_bins::Int = -1, relative_cut::Real = 0.2, fit_thresholds::Bool = true, uncertainty::Bool = true
)

Determines the typical trigger threshold from the vector thresholds between min_cut and max_cut. The vector thresholds is modified: values outside the range are removed.

If fit_thresholds is true, a truncated Gaussian is fitted to the histogram of the thresholds with fit_binned_trunc_gauss, using a window from cut_single_peak.

Keywords

• n_bins::Int = -1: Number of histogram bins. A value below 1 selects the Freedman-Diaconis bin width.
• relative_cut::Real = 0.2: Relative cut level of cut_single_peak.
• fit_thresholds::Bool = true: Fit the histogram. Otherwise μ and σ are the mean and standard deviation.
• uncertainty::Bool = true: Calculate fit uncertainties.

Returns

A tuple (result, report). result is the result of fit_binned_trunc_gauss (or μ and σ only without the fit) extended by μ_simple and σ_simple, the mean and standard deviation of the thresholds. report is the report of the fit, or the histogram without the fit.

fit_sipm_wl(
    trig_max_grid::VectorOfVectors{<:Real}, e_grid_wl::StepRangeLen,
    thresholds::Vector{<:Real} = zeros(length(e_grid_wl));
    min_pe_peak::Int = 1, max_pe_peak::Int = 5, n_fwhm_noise_cut::Real = 2.0,
    peakfinder_threshold::Real = 5.0, initial_max_amp::Real = 50.0,
    initial_max_bin_width_quantile::Real = 0.9999, peakfinder_rtol::Real = 0.1,
    peakfinder_α::Real = 0.1, peakfinder_σ::Real = -1.0,
    min_pe_fit::Real = 0.6, max_pe_fit::Real = 3.5, Δpe_peak_assignment::Real = 0.3
)

Finds the window length of the SiPM trigger filter with the best 1 PE separation.

For each window length in e_grid_wl, the trigger amplitudes are calibrated with sipm_simple_calibration and fitted with fit_sipm_spectrum. The gain is the distance between the first PE peak and the noise peak in uncalibrated units. The window length that minimizes √(resolution of the 1 PE peak) * √threshold / gain is selected (the threshold factor is 1 if all thresholds are zero). Window lengths that fail are skipped with a warning; an error is thrown if none succeeds.

Arguments

• trig_max_grid: Trigger amplitudes; element i belongs to e_grid_wl[i].
• e_grid_wl: Range of window lengths.
• thresholds: Trigger threshold for each window length, used as initial_min_amp of the calibration.

Keywords

The keywords min_pe_peak, max_pe_peak, n_fwhm_noise_cut, peakfinder_threshold, initial_max_amp, initial_max_bin_width_quantile, peakfinder_rtol, peakfinder_α and peakfinder_σ are passed to sipm_simple_calibration. min_pe_fit, max_pe_fit and Δpe_peak_assignment are passed to fit_sipm_spectrum as min_pe, max_pe and Δpe_peak_assignment.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• wl: optimal window length, with the step of e_grid_wl as uncertainty.
• obj: value of the objective at the optimum.
• res_1pe, pos_1pe: FWHM and position of the 1 PE peak, in uncalibrated units.
• gain: distance between the 1 PE peak and the noise peak.
• threshold: threshold at the optimum.

report contains wl, min_obj, the values gain, res_1pe and pos_1pe for all window lengths (undefined for failed ones), threshold, a_grid_wl_sg and obj for the successful window lengths, and report_simple and report_fit of the optimal window length, for plotting with LegendMakie.lplot.

fit_subpeaks_th228(
    h_survived::Histogram, h_cut::Histogram, h_result::NamedTuple;
    uncertainty::Bool = false, low_e_tail::Bool = true, fix_σ::Bool = true,
    fix_skew_fraction::Bool = true, fix_skew_width::Bool = true,
    pseudo_prior::NamedTupleDist = NamedTupleDist(empty = true),
    fit_func::Symbol = :gamma_def, background_center::Real = h_result.μ
)

Fits the two histograms h_survived and h_cut, which together form a histogram that was already fitted with fit_single_peak_th228 (result h_result), simultaneously. The fit determines the signal survival fraction sf and the background survival fraction bsf (the fraction of the peak and background of h_result that remains in h_survived).

The position, amplitude and background level are fixed to the values in h_result. fix_σ, fix_skew_fraction and fix_skew_width fix the peak width and the tail parameters of both histograms to the values in h_result; if false, they are fitted separately for each histogram.

Returns

A tuple (result, report). result holds the fit parameters as Measurements, in particular the signal survival fraction sf and the background survival fraction bsf (and sasf, the survival fraction of the step amplitude, if h_result has a step_amplitude), the fwhm, the fit_func and gof with converged and the entries survived and cut (pvalue, chi2, dof, covmat). report has the entries survived and cut (each with v, h, f_fit, f_components, gof) and sf and bsf.

gamma_peakshape(
    x::Real, μ::Real, σ::Real, n::Real,
    step_amplitude::Real, skew_fraction::Real, skew_width::Real,
    background::Real;
    skew_fraction_highE::Real = 0.0, skew_width_highE::Real = 0.0,
    background_kwargs...
)

Evaluates the shape of a gamma peak in an energy spectrum at x. The shape is the sum of signal_peakshape, lowEtail_peakshape, highEtail_peakshape and background_peakshape.

Arguments

• x: Energy at which the shape is evaluated.
• μ, σ: Mean and standard deviation of the Gaussian signal peak.
• n: Total number of counts in the peak (Gaussian part plus tails).
• step_amplitude: Amplitude of the step-like background from Compton-scattered gammas.
• skew_fraction, skew_width: Fraction of n in the low-energy tail and its exponential scale relative to μ (the scale is skew_width * μ).
• background: Energy-independent background level.

Keywords

• skew_fraction_highE = 0.0, skew_width_highE = 0.0: Fraction and relative exponential scale of the high-energy tail. The default disables the tail.
• background_kwargs...: Passed to background_peakshape (background_slope, background_exp, background_center).

LegendSpecFits.gauss_pdf(x::Real, μ::Real, σ::Real)

Evaluates the probability density function of a normal distribution with mean μ and standard deviation σ at x. Equivalent to pdf(Normal(μ, σ), x).

generate_aoe_compton_bands(
    aoe::Vector{<:Real}, e::Vector{<:T}, compton_bands::Vector{<:T}, compton_window::T
) where T<:Unitful.Energy{<:Real}

Generates A/E histograms for energy bands in the Compton region and estimates the peak parameters in each.

For each lower band edge c in compton_bands, the events with positive A/E and c < e < c + compton_window are selected and histogrammed around the A/E peak with Freedman-Diaconis bin widths. The peak statistics are determined with estimate_single_peak_stats (:psd variant). The histograms are then rebuilt in the range from 15 standard deviations below to 3 standard deviations above the peak position. Simple least-squares fits give the energy dependence of the peak position (linear) and width (√(p1 + p2 / E²)).

Returns

A NamedTuple with the fields

• peakhists: histogram of each band.
• peakstats: StructArray of peak statistics of each band.
• min_aoe, max_aoe: A/E range of the histogram of each band.
• e_unit: energy unit of the fits.
• mean_peak_pos, std_peak_pos: mean and standard deviation of the peak positions of all bands.
• mean_peak_sigma, std_peak_sigma: mean and standard deviation of the peak widths of bands 20 and above.
• simple_pars_aoe_μ, simple_pars_error_aoe_μ: parameters and standard errors of the linear fit of the peak position over energy.
• simple_pars_aoe_σ, simple_pars_error_aoe_σ: parameters and standard errors of the fit of the peak width over energy.

get_aoe_fit_functions(; background_center = nothing)

Returns a NamedTuple of the A/E Compton-region fit functions f(x, v), where v is a NamedTuple of the fit parameters:

• aoe_one_bck: aoe_compton_peakshape with parameters μ, σ, n, B (background amplitude), δ.
• aoe_two_bck: two_emg_aoe_compton_peakshape with the additional parameters μ2, σ2, B2f (fraction of the background in the second component), δ2.

background_center is not used.

get_centered_gaussian_window_cut(
    x::Vector{<:Unitful.RealOrRealQuantity}, min_x::T, max_x::T, n_σ::Real;
    center::T = zero(min_x), n_bins::Int = 500, relative_cut::Float64 = 0.2,
    left::Bool = false, fixed_center::Bool = true
) where T<:Unitful.RealOrRealQuantity

Determines a symmetric cut window of ±n_σ standard deviations around a peak in the values x between min_x and max_x.

The peak region is found with cut_single_peak. A one-sided truncated Gaussian is fitted to that region, with fit_half_centered_trunc_gauss (mean fixed to center) if fixed_center is true and with fit_half_trunc_gauss (free mean) otherwise. The fitted mean and standard deviation define the window.

Keywords

• center::T = zero(min_x): Mean of the Gaussian if fixed_center is true.
• n_bins::Int = 500: Number of histogram bins for cut_single_peak.
• relative_cut::Float64 = 0.2: Relative cut level for cut_single_peak.
• left::Bool = false: Fit the left side of the peak instead of the right side.
• fixed_center::Bool = true: Fix the mean of the Gaussian to center.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• low_cut, high_cut: lower and upper edge of the window, μ ∓ n_σ * σ.
• μ, σ: fitted mean and standard deviation.
• gof: goodness-of-fit result of the Gaussian fit.
• low_cut_fit, high_cut_fit: edges of the region used for the fit.
• max_cut_fit: position of the maximum bin.

All quantities carry the unit of min_x. report holds the fit function f_fit, the histogram h and the values μ, σ, gof, low_cut and high_cut, for plotting with LegendMakie.lplot.

Deprecated alias of get_continuum_survival_fraction.

get_continuum_survival_fraction(e, center, window, survival_flag; low_flag = survival_flag)
get_continuum_survival_fraction(cut_parameter, e, center, window; low_cut, high_cut, selection)

Counts the events in a continuum window before and after an event selection and calculates the survival fraction.

The window is center ± window. Unlike for peaks, no fit is made.

Arguments

• e: Energies of all events, as Unitful.Energy values.
• center, window: Center and half-width of the energy window.
• survival_flag::AbstractVector{Bool}: true for each event that survives.
• cut_parameter: Instead of survival_flag, a parameter (such as A/E) from which the selection is built: an event survives if low_cut < cut_parameter < high_cut and selection is true.

Keywords

• low_flag::AbstractVector{Bool} = survival_flag: Events that pass the lower cut only; used for the report histogram h_after_low.
• low_cut, high_cut: Exclusive bounds on cut_parameter. Default is no bound. h_after_low applies low_cut only, h_after_ds applies the final selection including high_cut and selection.
• selection::AbstractVector{Bool}: Additional selection. Default selects all events.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• window: center ± window.
• n_before, n_after: event counts with Poisson uncertainty √n.
• sf: survival fraction in percent, with binomial uncertainty.

report contains the histograms h_before, h_after_low and h_after_ds (in keV) and window, n_before, n_after, sf, e_unit and bin_width, for plotting with LegendMakie.lplot.

get_distribution_transform(d::Distribution, pprior::Prior)

Returns a DistributionTransform for the given Distribution and Prior.

This package defines no methods for this function.

get_friedman_diaconis_bin_width(x::AbstractVector{<:Real})
get_friedman_diaconis_bin_width(x::AbstractVector{<:Quantity{<:Real}})

Returns the bin width $2\,\mathrm{IQR}(x)/\sqrt[3]{n}$ of the Freedman-Diaconis rule, where IQR is the interquartile range of x and n is the number of values. For values with units, the result carries the same unit.

get_low_aoe_cut(
    aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T};
    dep::T = 1592.53u"keV", window::Vector{<:T} = [12.0, 10.0]u"keV", dep_sf::Float64 = 0.9,
    rtol::Float64 = 0.001, maxiters::Int = 300, sigma_high_sided::Float64 = Inf,
    cut_search_interval::Tuple = (-25.0 * unit(first(aoe)), 1.0 * unit(first(aoe))),
    bin_width_window::T = 3.0u"keV", max_e_plot::T = 3000.0u"keV",
    plot_window::Vector{<:T} = [12.0, 50.0]u"keV",
    fixed_position::Bool = true, fit_func::Symbol = :gamma_def, uncertainty::Bool = true
) where T<:Unitful.Energy{<:Real}

Determines the lower A/E cut value at which a fraction dep_sf of the events in the double escape peak (DEP) survives.

The DEP at energy dep is fitted with fit_single_peak_th228 in the window [dep - window[1], dep + window[2]]. A bisection root search in cut_search_interval with relative tolerance rtol then finds the cut value for which the survival fraction from get_sf_after_aoe_cut equals dep_sf.

Keywords

• dep::T = 1592.53u"keV": Energy of the double escape peak.
• window::Vector{<:T} = [12.0, 10.0]u"keV": Lower and upper width of the fit window around dep.
• dep_sf::Float64 = 0.9: Target survival fraction of the peak.
• rtol::Float64 = 0.001, maxiters::Int = 300: Relative tolerance and maximum number of iterations of the root search.
• sigma_high_sided::Float64 = Inf: Upper A/E cut value (in the unit of aoe); Inf applies no upper cut.
• cut_search_interval: Interval of A/E values searched for the lower cut.
• bin_width_window::T = 3.0u"keV": Half-width of the energy window around dep that is used to determine the Freedman-Diaconis bin width.
• max_e_plot::T = 3000.0u"keV", plot_window::Vector{<:T} = [12.0, 50.0]u"keV": Energy range of the full-spectrum histograms and lower and upper width of the DEP histograms in report.
• fixed_position::Bool = true: Fix the fitted peak position to dep.
• fit_func::Symbol = :gamma_def: Peak shape (see fit_single_peak_th228).
• uncertainty::Bool = true: Calculate fit uncertainties.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• lowcut: lower A/E cut value, with rtol * lowcut as uncertainty.
• highcut: upper A/E cut value (sigma_high_sided, with the unit of aoe).
• n0: number of counts in the peak before the cut (from the fit).
• nsf: number of counts in the peak after the cut.
• sf: survival fraction after the cut, in percent.

report contains the histograms h_before, h_after_low (lower cut only) and h_after_ds (lower and upper cut) of the full spectrum and dep_h_before, dep_h_after_low, dep_h_after_ds around the DEP, together with sf, n0, lowcut, highcut, e_unit and bin_width, for plotting with LegendMakie.lplot.

get_mc_value_shapes(v::NamedTuple, v_err::NamedTuple, n::Integer)

Draws n samples of the parameters in v, each parameter from an independent Normal distribution with mean v[k] and standard deviation v_err[k]. Returns a vector of n NamedTuples with the keys of v.

get_mc_value_shapes(v::NamedTuple, v_err::Matrix, n::Integer)

Draws n samples of the fit parameters from a multivariate normal distribution with mean v and covariance matrix v_err.

If v_err is not positive definite, the nearest positive definite matrix (see nearestSPD) is used. Samples are discarded unless the third parameter (amplitude) is positive, the fifth parameter (skew fraction) lies in (0, 0.25) and the sixth parameter (skew width) is positive, so the returned vector can be shorter than n. These checks assume the parameter order of the single-peak fit models.

get_number_of_bins(x::AbstractArray; method::Symbol = :sqrt)

Returns the number of histogram bins for the data x as an Int.

Keywords

• method::Symbol = :sqrt: Rule for the bin count. One of :sqrt, :sturges, :rice, :scott, :doane or :fd (Freedman-Diaconis). Any other value throws an AssertionError.

get_peak_fwhm_aoe_compton(v_ml::NamedTuple, v_ml_err::Union{Matrix, NamedTuple}, f_fit::Function, uncertainty::Bool = true)

Calculates the FWHM of an A/E peak from the best-fit parameters v_ml, with the fit function f_fit(x), and its uncertainty by Monte Carlo error propagation.

v_ml_err is either the covariance matrix of the parameters (10 000 correlated samples) or a NamedTuple of standard deviations (1000 independent samples). With uncertainty = false, no sampling is done and the uncertainty is NaN.

Returns

A tuple (fwhm, fwhm_err).

get_peak_fwhm_th228(v_ml::NamedTuple, v_ml_err::Union{Matrix, NamedTuple}, uncertainty::Bool = true)

Calculates the FWHM of a gamma peak from the best-fit parameters v_ml and its uncertainty by Monte Carlo error propagation.

v_ml_err is either the covariance matrix of the parameters (10 000 correlated samples) or a NamedTuple of standard deviations (1000 independent samples). With uncertainty = false, no sampling is done and the uncertainty is NaN.

Returns

A tuple (fwhm, fwhm_err) of the FWHM and its standard deviation over the finite samples.

Deprecated alias of get_peak_survival_fraction.

get_peak_survival_fraction(h_before, h_survived, h_cut; uncertainty, fit_func)
get_peak_survival_fraction(e, survival_flag; uncertainty, fit_func)
get_peak_survival_fraction(e, peak, window, survival_flag; uncertainty, bin_width_window, fit_func)
get_peak_survival_fraction(cut_parameter, e, peak, window; low_cut, high_cut, selection, kwargs...)

Fits a peak before and after an event selection and calculates the fraction of the peak counts that survives.

The peak is first fitted with fit_single_peak_th228 before the selection and then the surviving and rejected events are fitted simultaneously with fit_subpeaks_th228.

Arguments

• h_before, h_survived, h_cut: One-dimensional histograms (unitless energy) of all, surviving and rejected events with the same binning. The latter two must add up to the first.
• e: Energies of the events in the peak region, as Unitful.Energy values. Without peak and window, the whole range of e is histogrammed and fitted.
• survival_flag::AbstractVector{Bool}: true for each event that survives.
• peak, window: Peak energy and a tuple (lower, upper) of window widths below and above it.
• cut_parameter: Instead of survival_flag, a parameter (such as A/E) from which the selection is built: an event survives if low_cut < cut_parameter < high_cut and selection is true.

Keywords

• uncertainty::Bool = true: Calculate fit uncertainties.
• fit_func::Symbol = :gamma_def: Peak shape (see fit_single_peak_th228).
• bin_width_window::T = 2.0u"keV": Half-width of the energy window around peak that is used to determine the Freedman-Diaconis bin width.
• low_cut, high_cut: Exclusive bounds on cut_parameter. Default is no bound.
• selection::AbstractVector{Bool}: Additional selection. Default selects all events.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• peak: the peak energy (methods with peak only).
• fit_func: the peak shape that was used.
• n_before, n_after: number of counts in the peak before and after the selection.
• sf: survival fraction in percent.
• gof: goodness-of-fit results before and after.

report contains peak (methods with peak only), n_before, n_after, sf and the reports before and after of the fits, for plotting with LegendMakie.lplot.

Deprecated alias of get_peaks_survival_fractions.

get_peaks_survival_fractions(e, peaks, peak_names, windows, survival_flag; uncertainty, bin_width_window, fit_funcs)
get_peaks_survival_fractions(e, peaks, peak_names, left_window_sizes, right_window_sizes, survival_flag; kwargs...)
get_peaks_survival_fractions(cut_parameter, e, peaks, peak_names, windows; low_cut, high_cut, selection, kwargs...)
get_peaks_survival_fractions(cut_parameter, e, peaks, peak_names, left_window_sizes, right_window_sizes; kwargs...)

Fits the survival fraction of several peaks after an event selection, with get_peak_survival_fraction for each peak. The peaks are processed in parallel on the available threads.

Arguments

• e: Energies of all events, as Unitful.Energy values.
• peaks: Peak energies.
• peak_names::AbstractVector{Symbol}: Names used as keys of the result.
• windows: One tuple (lower, upper) of window widths below and above each peak. Alternatively, the vectors left_window_sizes and right_window_sizes hold the widths.
• survival_flag::AbstractVector{Bool}: true for each event that survives the selection.
• cut_parameter: Instead of survival_flag, a parameter (such as A/E) from which the selection is built: an event survives if low_cut < cut_parameter < high_cut and selection is true.

Keywords

• uncertainty::Bool = true: Calculate fit uncertainties.
• bin_width_window::T = 2.0u"keV": Half-width of the energy window around each peak that is used to determine the Freedman-Diaconis bin width.
• fit_funcs::AbstractVector{Symbol}: Peak shape for each peak (see fit_single_peak_th228). Default is :gamma_def for all peaks.
• low_cut, high_cut: Exclusive bounds on cut_parameter (cut-parameter methods). Default is no bound.
• selection::AbstractVector{Bool}: Additional selection (cut-parameter methods). Default selects all events.

Returns

A tuple (result, report) of two OrderedDicts keyed by peak_names. The values are the result and report NamedTuples of get_peak_survival_fraction.

get_residuals(f_fit::Base.Callable, h::Histogram{<:Real,1}, v_ml::Union{NamedTuple, AbstractVector})

Calculates bin-wise residuals and normalized residuals of the histogram h with respect to the model f_fit(x, v_ml), and a bin-wise p-value based on the Poisson distribution of each bin. Only bins with positive expected counts are included.

Returns

A tuple (residuals, residuals_norm, p_value_binwise, bin_centers) with

• residuals: model counts minus data counts in each bin.
• residuals_norm: residuals divided by the square root of the model counts.
• p_value_binwise: probability, for a Poisson distribution with the model counts as mean, of a deviation from the mean at least as large as observed.
• bin_centers: centers of the included bins.

get_sf_after_aoe_cut(
    aoe_cut::Unitful.RealOrRealQuantity, aoe::Vector{<:Unitful.RealOrRealQuantity},
    e::Vector{<:T}, peak::T, window::Vector{T}, bin_width::T, result_before::NamedTuple;
    uncertainty::Bool = true, fit_func::Symbol = :gamma_def
) where T<:Unitful.Energy{<:Real}

Calculates the signal survival fraction of a peak after the A/E cut value aoe_cut (events with aoe >= aoe_cut survive).

The energies e in the window [peak - window[1], peak + window[2]] are histogrammed with the bin width bin_width separately for the surviving and the cut events. Both histograms are fitted simultaneously with fit_subpeaks_th228, using the fit result result_before of the peak before the cut as fixed parameters.

Returns

The survival fraction sf as a fraction between 0 and 1 (a Measurement if uncertainty is true).

get_th228_fit_functions(; background_center::Union{Real, Nothing} = nothing)

Returns a NamedTuple of the gamma peak shape functions f(x, v) used in the calibration fits, where v is a NamedTuple of the fit parameters. Each is built from gamma_peakshape:

• gamma_def: Gaussian signal, low-energy tail, flat background and step background.
• gamma_tails: gamma_def plus a high-energy tail.
• gamma_sigWithTail: signal and low-energy tail only, without background.
• gamma_bckSlope: gamma_def plus a linear background slope.
• gamma_bckExp: gamma_def plus an exponential background.
• gamma_bckFlat: gamma_def without the step background.
• gamma_tails_bckFlat: gamma_tails without the step background.
• gamma_minimal: Gaussian signal and flat background only.

background_center is the reference position of the slope and exponential background terms. If nothing, the peak position v.μ is used.

heaviside(x)

Returns zero(x) for x < 0 and one(x) otherwise (including x == 0).

highEtail_peakshape(
    x::Real, μ::Real, σ::Real, n::Real,
    skew_fraction_h::Real, skew_width_h::Real
)

Evaluates the high-energy tail of a gamma peak at x: n * skew_fraction_h times an exponentially modified Gaussian that extends toward higher energies, with exponential scale skew_width_h * μ. Returns zero if σ is zero.

hist_loglike(f_fit::Base.Callable, h::Histogram{<:Real,1})

Calculates the Poisson log-likelihood of the bin counts of the histogram h given the density f_fit(x). f_fit must accept all values x on the horizontal axis of h.

The expected count in each bin is f_fit at the bin center times the bin width (midpoint rule).

When the same histogram is evaluated many times, construct a HistLogLike once instead.

likelihood_ratio(f_fit::Base.Callable, h::Histogram{<:Real,1})

Calculates the likelihood-ratio statistic -2 (log L_model - log L_saturated) of the histogram h, where L_model is the Poisson likelihood for the expected counts bin_width * f_fit(bin_center) and L_saturated is the likelihood with the expected counts set equal to the observed counts.

linear_function(x::Real, slope::Real, intercept::Real)

Evaluates the linear function slope * x + intercept at x.

lowEtail_peakshape(
    x::Real, μ::Real, σ::Real, n::Real,
    skew_fraction::Real, skew_width::Real
)

Evaluates the low-energy tail of a gamma peak at x: n * skew_fraction times an exponentially modified Gaussian that extends toward lower energies, with exponential scale skew_width * μ. Returns zero if σ is zero.

lq_ctc_correction(
    lq::Vector{<:AbstractFloat}, dt_eff::Vector{<:Unitful.RealOrRealQuantity},
    e_cal::Vector{<:Unitful.Energy{<:Real}}, dep_μ::Unitful.AbstractQuantity, dep_σ::Unitful.AbstractQuantity;
    ctc_dep_edgesigma::Float64 = 3.0, ctc_lq_precut_relative_cut::Float64 = 0.25,
    lq_outlier_sigma::Float64 = 2.0, ctc_driftime_cutoff_method::Symbol = :percentile,
    dt_eff_outlier_sigma::Float64 = 2.0, lq_e_corr_expression = "lq / e",
    dt_eff_expression = "qdrift / e", ctc_dt_eff_low_quantile::Float64 = 0.15,
    ctc_dt_eff_high_quantile::Float64 = 0.95, pol_fit_order::Int = 1, uncertainty::Bool = false
)

Corrects the drift time dependence of the LQ parameter with the double escape peak (DEP).

The events within dep_μ ± ctc_dep_edgesigma * dep_σ are selected. Outliers in LQ and in the effective drift time are removed, and a polynomial is fitted to LQ as a function of drift time with chi2fit. The correction subtracts this polynomial from lq.

Arguments

• lq: Energy-corrected LQ parameter.
• dt_eff: Effective drift time.
• e_cal: Calibrated energy.
• dep_μ, dep_σ: Mean and standard deviation of the DEP.

Keywords

• ctc_dep_edgesigma: Half-width of the DEP selection in units of dep_σ.
• ctc_lq_precut_relative_cut: Relative cut level of cut_single_peak for the LQ precut.
• lq_outlier_sigma: LQ values farther than this many standard deviations from the fitted LQ mean are removed.
• ctc_driftime_cutoff_method: Drift time range selection. :percentile uses the quantiles ctc_dt_eff_low_quantile and ctc_dt_eff_high_quantile; :gaussian uses a truncated Gaussian fit with dt_eff_outlier_sigma standard deviations (not for detectors with a double-peaked drift time distribution); :double_gaussian uses the positions where a fitted double Gaussian falls to 10 % of its maximum.
• dt_eff_outlier_sigma: Used by :gaussian.
• lq_e_corr_expression, dt_eff_expression: Names of the LQ and drift-time variables in the returned func.
• ctc_dt_eff_low_quantile, ctc_dt_eff_high_quantile: Used by :percentile.
• pol_fit_order: Order of the polynomial.
• uncertainty: Calculate fit uncertainties.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• func: the corrected LQ as a Julia expression string.
• func_generic: the same with generic parameter names.
• fit_result: result of the polynomial fit (chi2fit).
• box_constraints: NamedTuple (lq_lower, lq_upper, t_lower, t_upper) of the LQ and drift-time ranges used.

report contains lq_report and drift_report (reports of the fits of LQ and drift time; drift_report is nothing for :percentile), lq_box, drift_time_func (the fitted polynomial), dep_left and dep_right.

lq_norm(
    dep_μ::Unitful.Energy, dep_σ::Unitful.Energy, e_cal::Vector{<:Unitful.Energy},
    lq_classifier::Vector{<:AbstractFloat};
    dep_sideband_sigma::Float64 = 4.5, cut_truncation_sigma::Float64 = 3.5,
    uncertainty::Bool = true, lq_class_expression = "lq / e  - (slope * qdrift / e + y_inter)"
)

Normalizes the charge-trapping-corrected LQ classifier to mean 0 and standard deviation 1 using the double escape peak (DEP).

The LQ histograms of two sidebands next to the DEP (bounds at multiples of dep_sideband_sigma * dep_σ from dep_μ; both sides if dep_σ < 2 keV, both on the low-energy side otherwise) are subtracted from the DEP histogram, negative bins are set to zero, and a truncated Gaussian is fitted to the result with fit_binned_trunc_gauss.

Arguments

• dep_μ, dep_σ: Mean and standard deviation of the DEP.
• e_cal: Calibrated energies.
• lq_classifier: LQ classifier values, typically after lq_ctc_correction.

Keywords

• dep_sideband_sigma: DEP region and sideband width in units of dep_σ.
• cut_truncation_sigma: Fit range around the peak of the subtracted histogram, in standard deviations.
• uncertainty: Calculate fit uncertainties.
• lq_class_expression: Name of the LQ classifier in the returned func.

Returns

A tuple (result, report). result is a NamedTuple with fit_result (result of the Gaussian fit) and func, the normalized LQ as a Julia expression string. report contains fit_result, fit_report, temp_hists (histograms hist_dep, hist_sb1, hist_sb2, hist_subtracted, hist_corrected), dep_σ and the energy edges of the regions.

muncert(x)

Returns the uncertainty of a Measurements.Measurement. Applied element-wise to arrays and field-wise to NamedTuples; strings and symbols are returned unchanged.

mvalue(x)

Returns the central value of a Measurements.Measurement. Applied element-wise to arrays and field-wise to NamedTuples; strings and symbols are returned unchanged.

nearestSPD(A::Matrix{<:Real}; n_iter::Signed = 1)

Returns the symmetric positive definite matrix nearest to A.

The calculation follows the singular value decomposition method of Higham (1988). The step is repeated recursively, with n_iter counting the repetitions, until the result is positive definite or n_iter exceeds 5.

neg_log_likelihood_single_aoe_compton_with_fixed_μ_and_σ(
    h::Histogram, μ::Real, σ::Real, ps::NamedTuple;
    fit_func::Symbol = :aoe_one_bck, background_center = μ, optimize::Bool = true
)

Returns the negative log-likelihood of the A/E histogram h after fitting all parameters except the fixed μ and σ, as in fit_single_aoe_compton_with_fixed_μ_and_σ. With optimize = false, the likelihood is evaluated at the parameters (μ, σ) merged with ps instead.

objective_value(metrics::NamedTuple, objective::Symbol; ctc_rel_uncertainty = nothing)

Computes the value of an optimization objective from the metrics of one filter setting.

Arguments

• metrics::NamedTuple: Row with the fields of peak_width_metrics, as returned by it or by evaluate_energy_point. :fwhm_ct also needs the field ct_spread.
• objective::Symbol: One of
  • :w90: w90.
  • :fwhm: fwhm.
  • :asym: asym.
  • :fwhm_asym: fwhm * asym.
  • :fwhm_ct: sqrt(fwhm^2 + (k * ct_spread)^2) with k = ctc_rel_uncertainty.

Keywords

• ctc_rel_uncertainty = nothing: Relative uncertainty k of the charge-trapping correction, a finite number >= 0. Required for :fwhm_ct, ignored by the other objectives.

Returns

The value of the objective, which is minimized by the optimization. Throws an ArgumentError for an unknown objective and for :fwhm_ct without a valid ctc_rel_uncertainty.

optimize_rtft(
    f_energy, qdrift::AbstractVector{<:Real}, rt0::Unitful.Time, ft0::Unitful.Time;
    rt_bounds::Tuple, ft_bounds::Tuple, feasible = (rt, ft) -> true, step0::Tuple, lattice::Unitful.Time,
    objective::Symbol = :fwhm, ctc_rel_uncertainty = nothing, rel_tol::Real = 0.005, max_evals::Int = 60,
    peak, window, min_e, max_e, rel_cut_fit, n_bins = -1, apply_ctc = true, containment = 0.9
)

Searches the rise time and flat-top time of a trapezoidal filter jointly for the minimum of an objective.

Every candidate (rt, ft) lies on the lattice of multiples of lattice and is evaluated at most once with evaluate_energy_point. The search is a compass pattern search with diagonal moves, starting at (rt0, ft0) rounded to the lattice. The admissible neighbors at the current steps are evaluated in parallel with Threads.@threads. A neighbor is admissible if it lies inside the bounds, satisfies feasible, and has not been evaluated. The search moves to the best neighbor if its objective is below (1 - rel_tol) times the current one. Otherwise it halves both steps, down to one lattice step. It stops when both steps equal the lattice and no neighbor improves, or after max_evals evaluations.

The optimum is selected with select_optimum among the final point and its evaluated neighbors: points within rel_tol of the best objective tie, and the tie goes to the shortest rt, then the longest ft. An evaluation that throws is recorded with its error text in status and NaN metrics, and a warning is logged.

Arguments

• f_energy: Function f_energy(rt, ft) that returns the energies (ADC) of the waveform set filtered with this setting, aligned with qdrift. It is called from several threads concurrently.
• qdrift: Drift-time parameter of every event.
• rt0, ft0: Seed rise time and flat-top time.

Keywords

• rt_bounds::Tuple, ft_bounds::Tuple: Allowed range (min, max) of the rise time and the flat-top time. Required.
• feasible = (rt, ft) -> true: Function that returns false for settings that must not be evaluated.
• step0::Tuple: Initial steps (rt step, ft step), rounded to multiples of lattice (at least one). Required.
• lattice::Unitful.Time: Grid spacing of all candidates. Required.
• objective::Symbol = :fwhm: Objective of objective_value.
• ctc_rel_uncertainty = nothing: Required for :fwhm_ct (which also needs apply_ctc = true), see objective_value.
• rel_tol::Real = 0.005: Relative improvement needed to move, and relative tolerance of the tie.
• max_evals::Int = 60: Maximum number of evaluations.
• peak, window, min_e, max_e, rel_cut_fit: Required, see evaluate_energy_point.
• n_bins = -1, apply_ctc = true, containment = 0.9: See evaluate_energy_point.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• rt, ft: selected rise time and flat-top time.
• fwhm, w90, asym, fct, ct_spread: fields of evaluate_energy_point at the optimum.
• min_objective: objective value at the optimum.
• objective: the objective used.
• n_evals: number of evaluations, including failed ones.
• converged: true if both steps reached the lattice and no evaluated neighbor improved; false if max_evals was reached or all admissible neighbors failed.

report, for plotting with LegendMakie.lplot, contains

• rt, ft, objective: as in result.
• table: TypedTables.Table with one row per evaluation: the columns rt, ft, stage (:seed or :neighbor), step_rt, step_ft and status ("ok" or the error text), followed by the fields of evaluate_energy_point.
• path: accepted (rt, ft) positions, starting with the seed.
• objective_values: objective for every row of table.

Throws an ArgumentError for invalid arguments or a seed that is infeasible or out of bounds after rounding, and an ErrorException if the seed evaluation fails.

p_value(f_fit::Base.Callable, h::Histogram{<:Real,1}, v_ml::Union{NamedTuple, AbstractVector})

Calculates the p-value of a fit with a chi-squared test that assumes Gaussian uncertainties, the model counts being the variance of each bin.

Only bins with at least one count enter the test. f_fit(x, v_ml) is the fit function and v_ml the best-fit parameters. For bins with few counts the Gaussian assumption is poor; see p_value_poissonll for a test that stays valid there.

Returns

A tuple (pval, chi2, dof) with the p-value, the chi-squared value and the degrees of freedom (number of non-empty bins minus the number of parameters). pval is NaN if dof is negative.

p_value_MC(
    f_fit::Base.Callable, h::Histogram{<:Real,1}, ps::NamedTuple, v_ml::NamedTuple;
    n_samples::Int64 = 1000
)

Calculates a p-value with Monte Carlo sampling. This is computationally more expensive than p_value and p_value_poissonll.

n_samples histograms are drawn with Poisson-distributed bin counts around the best-fit model f_fit(x, v_ml). Each sample is fitted with fit_single_peak_th228, starting from the peak statistics ps, and the negative log-likelihood of the sample fit is calculated. The p-value is the fraction of samples whose negative log-likelihood is at most that of h.

p_value_poissonll(f_fit::Base.Callable, h::Histogram{<:Real,1}, v_ml::Union{NamedTuple, AbstractVector})

Calculates the p-value of a maximum-likelihood fit with Poisson statistics with the moment-corrected deviance test. The Poisson deviance D of poisson_deviance is compared with its mean and variance under the fitted model, the mean reduced by the number of fitted parameters, and the p-value is the upper tail of the normal distribution at z = (D − (E[D] − n_par)) / sqrt(Var[D]). Unlike the χ² approximation of the deviance this stays valid for bins with few or no counts (Kaastra 2017, A&A 605, A51, for the equivalent C statistic).

Returns

A tuple (pval, chi2, dof) with chi2 = D and dof the number of bins minus the number of parameters. pval is 0 if D is not finite.

peak_centroid(v::NamedTuple)

Calculates the centroid (mean position) of a gamma peak from the fit parameters v: the Gaussian position μ shifted by the low-energy tail, -skew_fraction * μ * skew_width, and by the high-energy tail, +skew_fraction_highE * μ * skew_width_highE. A tail is only included if its parameters are fields of v.

peak_search_gamma(
    e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}};
    peakfinder_σ::Real = 2.0, peakfinder_threshold::Real = 10.0,
    peak_quantile::ClosedInterval{<:Real} = 0.0..1.0,
    bin_quantile::ClosedInterval{<:Real} = peak_quantile, quantile_perc::Float64 = NaN
)

Searches for the most prominent peak in the uncalibrated spectrum e_uncal and returns a simple calibration constant that assigns it to the line in gamma_lines.

Keywords

• peakfinder_σ::Real = 2.0: Expected standard deviation of a peak, in bins (see RadiationSpectra.peakfinder).
• peakfinder_threshold::Real = 10.0: Threshold in the deconvoluted spectrum above which a bin is identified as a peak if the previous bin was not (see RadiationSpectra.peakfinder).
• peak_quantile = 0.0..1.0: Quantile range of e_uncal that is searched for peaks. All gamma_lines must lie within this window.
• bin_quantile = peak_quantile: Quantile range of e_uncal used to choose the bin width.
• quantile_perc::Float64 = NaN: If not NaN, the peak is set to this quantile of e_uncal and no peak search is done.

If more than one line is given, the number of peaks found must equal the number of lines, otherwise an error is thrown.

Returns

A NamedTuple with the fields

• h_uncal: histogram of e_uncal with the bin width used in the peak search.
• c: calibration constant (energy of the line divided by peak_guess).
• peak_guess: position of the identified peak in units of e_uncal.

peak_width_metrics(v::NamedTuple; containment::Real = 0.9)

Computes width and shape metrics of the signal part of a gamma peak.

The signal part is the Gaussian plus low-energy tail of the peak shape, without step and background. The metrics are evaluated on the normalized trapezoid-rule cumulative of this shape on a uniform grid with a step of at most σ / 50.

Arguments

• v::NamedTuple: Peak parameters μ, σ, n, skew_fraction and skew_width, for example the fit result of fit_single_peak_th228. Plain numbers and Measurements are accepted; uncertainties are ignored.

Keywords

• containment::Real = 0.9: Fraction of the signal shape that the interval w90 contains, in (0, 1).

Returns

A NamedTuple with the fields

• fwhm: full width at half maximum of the signal shape, as LegendSpecFits.estimate_fwhm.
• w90: length of the shortest interval that contains the fraction containment of the shape.
• asym: integral of the shape below μ divided by the integral above μ.
• tail_fraction: skew_fraction.
• tail_scale: decay length of the tail, skew_width * μ.
• centroid: mean of the shape, as peak_centroid.

All widths and positions have the units of v. Throws an ArgumentError for non-finite or out-of-range parameters or containment, and an ErrorException if the FWHM cannot be determined.

peakhists_gamma(
    e::Vector{<:Unitful.Energy{<:Real}}, gamma_lines::Vector{<:Unitful.Energy{<:Real}},
    window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}};
    binning_peak_window::Unitful.Energy{<:Real} = 10.0u"keV"
)
peakhists_gamma(e, gamma_lines, left_window_sizes, right_window_sizes; binning_peak_window)

Creates histograms of the energies e around each of the gamma_lines and estimates the peak statistics in each window.

Arguments

• e: Calibrated energies.
• gamma_lines: Gamma line energies.
• window_sizes: One tuple (left, right) per line with the window sizes below and above the line. Alternatively, vectors of left and right window sizes can be passed separately.

Keywords

• binning_peak_window::Unitful.Energy = 10.0u"keV": Energy window around each line used to choose the Freedman-Diaconis bin width.

Returns

A tuple (peakhists, peakstats, h, bin_widths) with

• peakhists: vector of histograms around the lines (values in e_unit, keV).
• peakstats: StructArray of the results of estimate_single_peak_stats for each window.
• h: histogram of all of e from 0 to 3000 keV with the median bin width.
• bin_widths: bin width of each peak histogram with energy unit.

peakshape_components(fit_func::Symbol; background_center::Union{Real, Nothing} = nothing)
peakshape_components(fit_func::Symbol, v::NamedTuple; background_center = v.μ)

Returns the component functions of the peak shape fit_func (one of the names in get_th228_fit_functions except gamma_sigWithTail) as a NamedTuple with the fields

• funcs: component functions f(x, v) (f_sig, f_lowEtail, f_bck and, for the shapes with a high-energy tail, f_highEtail). In the method with v, these are functions of x only, with the parameters v applied.
• labels, colors, linestyles: plot label, color and line style of each component.

The components are used in fit reports and plot recipes.

poisson_deviance(f_fit::Base.Callable, h::Histogram{<:Real,1}, v_ml::Union{NamedTuple, AbstractVector})

Calculates the Poisson deviance of the histogram h under the model f_fit(x, v_ml), D = 2 Σ_i [λ_i − k_i + k_i log(k_i / λ_i)], with k_i the counts and λ_i the bin width times the model at the bin center (the likelihood-ratio statistic of Baker and Cousins). Also calculates its mean E[D] and variance Var[D] under the model, summed from the exact Poisson moments of every bin's term. A bin whose expectation vanishes contributes nothing when it is empty and Inf otherwise.

Returns

A tuple (D, E, V) of the deviance, its mean and its variance.

prepare_sep_peakhist(
    e::Vector{<:T}; sep::Unitful.Energy{<:Real} = 2103.53u"keV",
    window::Vector{<:Unitful.Energy{<:Real}} = [25.0, 25.0]u"keV", relative_cut::T = 0.5,
    n_bins_cut::Int = -1, fit_func::Symbol = :gamma_def, uncertainty::Bool = true
) where T<:Real

Fits the single escape peak (SEP) at energy sep in the uncalibrated energies e and derives a calibration constant from its position.

A simple calibration from cut_single_peak selects the histogram range sep ± window, and the peak is fitted with fit_single_peak_th228.

Returns

A tuple (result, report) of the fit result and report, where result is extended by m_calib, the calibration constant sep / μ.

qc_window_cut(tab::Table, config::AbstractDict)
qc_window_cut(x::Vector{<:Unitful.RealOrRealQuantity}, min::T, max::T, sigma::Real; col_expression = "col", kwargs...)

Determines quality-control window cuts on one or more columns, each from a centered Gaussian fit with get_centered_gaussian_window_cut.

For a table, config maps each column name (a Symbol) to a dictionary with the entries :min, :max (range of the data used), :sigma (half-width of the window in standard deviations) and :kwargs (keyword arguments of get_centered_gaussian_window_cut as a dictionary or NamedTuple). The columns are processed in parallel on the available threads. For a vector x, these values are passed directly and kwargs are forwarded.

Keywords

• col_expression: Name of the variable in the returned cut string (vector method).

Returns

A tuple (result, report).

• Table method: result is a NamedTuple with the field qc, a string of the conditions of all columns joined by &&, and one entry per column holding its result; report is an OrderedDict keyed by column name.
• Vector method: result is the result of get_centered_gaussian_window_cut extended by qc, the cut condition string <low_cut> < col_expression && col_expression < <high_cut> (with units if present); report is its report (f_fit, h, μ, σ, gof, low_cut, high_cut).

select_optimum(objective_values::AbstractVector{<:Real}, keys::AbstractVector{<:Tuple}; rel_tol::Real)

Selects the optimum among candidates with tied objective values.

All candidates whose objective is at most minimum(objective_values) * (1 + rel_tol) tie. Among the tied candidates, the one with the lexicographically smallest tuple in keys is selected. A negated entry in the tuple prefers the larger value of that quantity. NaN objectives are never selected.

Arguments

• objective_values: Objective of every candidate.
• keys: Tie-breaking tuple of every candidate, with the same indices as objective_values.

Keywords

• rel_tol::Real: Relative tolerance of the tie, in [0, Inf). Required.

Returns

The index of the selected candidate in eachindex(objective_values, keys). Throws an ArgumentError if all objectives are NaN, if the smallest objective is not positive and finite, or if rel_tol is outside [0, Inf).

set_memlimit(gig::Float64)

Sets the module-level memory limit, in GB, for the whole current Julia process, regardless of how many fits run in parallel. The default is 10 GB. The value is stored in LegendSpecFits.optim_max_mem.

set_timelimit(sec::Float64)

Sets the time limit, in seconds, for a single optimization step that uses LegendSpecFits.optim_time_limit as maxtime (the NLopt-based fits in ctc_energy and ctc_aoe). The default is 20 s.

signal_peakshape(
    x::Real, μ::Real, σ::Real, n::Real,
    skew_fraction::Real; skew_fraction_highE::Real = 0.0
)

Evaluates the Gaussian signal part of a gamma peak at x: n * (1 - skew_fraction - skew_fraction_highE) times the density of Normal(μ, σ). Returns zero if σ is zero.

simple_calibration(
    e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}},
    window_sizes::Vector{<:Tuple{Unitful.Energy{<:Real}, Unitful.Energy{<:Real}}};
    calib_type::Symbol, kwargs...
)
simple_calibration(
    e_uncal::Vector{<:Real}, gamma_lines::Vector{<:Unitful.Energy{<:Real}},
    left_window_sizes::Vector{<:Unitful.Energy{<:Real}},
    right_window_sizes::Vector{<:Unitful.Energy{<:Real}}; calib_type::Symbol, kwargs...
)

Determines a first, linear energy calibration of the uncalibrated energies e_uncal from the position of one prominent peak, and builds histograms and peak statistics around all gamma_lines.

The calibration constant is the ratio of the known energy of the selected line to the peak position found in the uncalibrated spectrum. The peak search histogram uses the Freedman-Diaconis bin width.

Arguments

• e_uncal: Uncalibrated energies.
• gamma_lines: Known gamma line energies, e.g. from 228Th.
• window_sizes: One tuple (left, right) per line giving the window below and above the line that is histogrammed. Alternatively, vectors of left and right window sizes can be passed separately.

Keywords

• calib_type::Symbol: Required. :th228 searches for the 2614.5 keV line and sets the peak-finder parameters peakfinder_σ = 5.0, peakfinder_threshold = 10.0, peak_quantile = 0.9..1.0 and bin_quantile = 0.05..0.5, overriding any user values. :gamma uses the keywords below for generic gamma lines. Other values throw an error.
• peaksearch_type::Symbol = :most_prominent: (:gamma only) :most_prominent uses the highest peak found, which requires the peak finder to find all gamma_lines; :single uses only gamma_lines[peakfinder_gamma].
• peakfinder_gamma::Int = 1: Index of the line used for peaksearch_type = :single.
• binning_peak_window::Unitful.Energy = 10.0u"keV": Energy window around each line used to choose the bin width of the peak histograms.
• peakfinder_σ::Real = 2.0, peakfinder_threshold::Real = 10.0: Expected peak width (in bins) and detection threshold of RadiationSpectra.peakfinder.
• peak_quantile, bin_quantile: Quantile ranges (ClosedInterval, default 0.0..1.0) of e_uncal for the peak search and for choosing the bin width.
• quantile_perc::Float64 = NaN: If not NaN, the peak is set to this quantile of e_uncal instead of being searched for.

Returns

A tuple (result, report). result is a NamedTuple with the fields

• h_calsimple: histogram of the roughly calibrated energies (0 to 3000 keV).
• h_uncal: histogram of the uncalibrated energies used in the peak search.
• c: calibration constant with energy unit (energy per uncalibrated unit).
• unit: energy unit of the calibrated values.
• bin_width, peakbinwidths: median and per-line bin width of the peak histograms.
• peak_guess: position of the identified peak in uncalibrated units.
• peakhists: vector of histograms around the lines (in unit).
• peakstats: StructArray of peak statistics from estimate_single_peak_stats.

report contains h_calsimple, h_uncal, c, peak_guess, peakhists and peakstats, for plotting with LegendMakie.lplot.

sipm_simple_calibration(pe_uncal::Vector{<:Real}; kwargs...)
sipm_simple_calibration(pe_uncal_vov::VectorOfVectors{<:Real}; kwargs...)

Determines a linear calibration of SiPM pulse amplitudes to photoelectron (PE) units from the positions of the 1 PE and 2 PE peaks found by a peak finder.

The calibration is pe = c * amplitude + offset, with c the inverse of the distance between the first two peaks and offset chosen such that the first peak lies at min_pe_peak. The noise peak is located with cut_single_peak. The threshold between noise and the 1 PE peak is placed at the center of the first valley after the noise peak; if the valley is too short, at the noise peak maximum plus 3 times its half-width, and if none is found, at the end of the search range. The trigger threshold is the lowest minimum of the Savitzky-Golay smoothed calibrated spectrum between the noise threshold and the 1 PE peak. Peaks are searched on a Savitzky-Golay smoothed histogram with RadiationSpectra.determine_calibration_constant_through_peak_ratios; if that fails, local maxima are used. Neighboring peaks closer than 60 % of the largest peak distance are merged.

The method for VectorOfVectors takes the trigger amplitudes of each waveform. It determines the threshold from the finite triggers (of waveforms with at most cut_pool_max_mult triggers if single_trigger_only is true), keeps the triggers above it (only waveforms with exactly one such trigger if single_trigger_only is true), and calibrates these with the vector method.

Keywords

• initial_min_amp::Real = 0.0, initial_max_amp::Real = 50.0: Amplitude range for the noise peak and peak search.
• relative_cut_noise_cut::Real = 0.5: Relative cut level of cut_single_peak for the noise peak.
• n_fwhm_noise_cut::Real = 5.0: Upper limit of the valley search in half-widths of the noise peak above its center. 0.0 skips the valley search and uses initial_min_amp as threshold. A negative value places the threshold at the noise peak maximum plus this factor times the noise peak half-width.
• min_pe_peak::Int = 1, max_pe_peak::Int = 5: First and last PE peak considered (vector method).
• initial_max_bin_width_quantile::Real = 0.9: Quantile of the amplitudes up to which the bin width is determined (vector method).
• peakfinder_σ::Real = -1.0: Expected peak width in bins; a non-positive value estimates it from the noise peak (must be an integer otherwise).
• peakfinder_threshold::Real = 10.0, peakfinder_rtol::Real = 0.1, peakfinder_α::Real = 0.05: Threshold, relative tolerance and significance level of the peak finder.
• single_trigger_only::Bool = true, cut_pool_max_mult::Int = 3: Waveform selection (VectorOfVectors method).

Returns

A tuple (result, report). result is a NamedTuple with the fields

• peakpos: peak positions in uncalibrated units.
• c, offset: slope and offset of the calibration, x_cal = x * c + offset.
• noisepeakpos, noisepeakwidth: position and width of the noise peak.
• noise_threshold, noise_threshold_cal: threshold between noise and the 1 PE peak, uncalibrated and calibrated.
• trigger_threshold, trigger_threshold_cal: the lowest minimum of the smoothed spectrum between the noise threshold and the 1 PE peak, uncalibrated and calibrated; equal to the noise threshold if there is no minimum below the spectrum at the noise threshold.

report contains peakpos, peakpos_cal, the histograms h_uncal and h_calsimple (uncalibrated and calibrated), h_uncal_full and h_calsimple_full (spectra of all triggers before the waveform selection; equal to h_uncal and h_calsimple for the vector method), the noise and trigger thresholds, valley_found, the calibrated amplitudes pe_simple_cal and the calibration function f_simple_calib with its inverse f_simple_uncal, for plotting with LegendMakie.lplot.

LegendSpecFits.stabilize_dist(d::Distribution)

Returns a numerically stable version of the distribution d.

For a Weibull distribution with infinite log-density at zero, the result is truncated to start at a small positive value. Other distributions are returned unchanged.

step_gauss(x::Real, μ::Real, σ::Real)

Evaluates the convolution of a Heaviside step function with the PDF of Normal(μ, σ) at x, i.e. the cumulative distribution function of Normal(μ, σ). It rises from 0 to 1 around μ.

The result is not a probability density because it is not normalizable.

subhist(h::Histogram{<:Any, 1}, r::Tuple{<:Real,<:Real})
subhist(h, i::Interval)
subhist(x::Array{<:Real})
subhist(x::Array, r::Tuple{<:Real,<:Real})

Returns a histogram restricted to a range.

• For a one-dimensional Histogram h, the result contains all bins from the bin that contains r[1] to the bin that contains r[2], clamped to the edges of h. An IntervalSets.Interval i is converted to the tuple (i.left, i.right).
• For a data vector x, the histogram is built from the values of x with the Freedman-Diaconis bin width, spanning minimum(x) to maximum(x). If r is given, only values strictly between r[1] and r[2] are used.

tuple_to_array(nt::NamedTuple)

Returns a Vector with the values of the fields of nt in field order.

two_emg_aoe_compton_background_peakshape(
    x::Real,
    μ::Real, σ::Real, background::Real, δ::Real,
    μ2::Real, σ2::Real, background2_fraction::Real, δ2::Real
)

Evaluates the background part of two_emg_aoe_compton_peakshape at x.

two_emg_aoe_compton_peakshape(
    x::Real,
    μ::Real, σ::Real, n::Real, background::Real, δ::Real,
    μ2::Real, σ2::Real, background2_fraction::Real, δ2::Real
)

Evaluates the A/E distribution of a Compton region with two background components at x: a Gaussian signal peak (μ, σ, n) plus background times a mixture of two mirrored exponentially modified Gaussians. The first has parameters μ, σ, δ and weight 1 - background2_fraction, the second has μ2, σ2, δ2 and weight background2_fraction. Returns zero if σ or σ2 is zero.

weibull_from_mx(m::Real, x::Real, p_x::Real = 0.6827)::Weibull

Constructs a Weibull distribution with median m whose p_x-quantile is x. The default p_x = 0.6827 corresponds to a one-sigma quantile.

The result is suited as a prior for positive quantities.