LegendSpecFits

Template for Julia packages.

DEP_µ::Unitful.Energy, DEP_σ::Unitful.Energy, e_cal::Vector{<:Unitful.Energy}, lq_classifier::Vector{Float64}; cut_sigma::Float64=3.0, truncation_sigma::Float64=2.0)

Evaluates the cutoff value for the LQ cut. The function performs a binned gaussian fit on the sidebandsubtracted LQ histogram and evaluates the cutoff value difined at 3σ of the fit.

Returns

* `result`: NamedTuple of the cutoff value
* `report`: NamedTuple of the fit result, fit report and temporary histograms

_fit_fwhm_ft(e_grid::Matrix, e_grid_ft::StepRangeLen{Quantity{<:T}, Base.TwicePrecision{Quantity{<:T}}, Base.TwicePrecision{Quantity{<:T}}, Int64}, rt::Unitful.RealOrRealQuantity, min_e::T, max_e::T, nbins::Int64, rel_cut_fit::T; default_ft::Quantity{T}=3.0u"µs", ft_fwhm_tol::Unitful.Energy{<:Real} = 0.1u"keV") where {T <:Real}

Fit the FWHM values in e_grid for each FT in e_grid_ft with a Gamma Peakshape and return the optimal FT and the corresponding FWHM value. The cut values cut for each flat-top time a window for better histogramming.

Arguments

• e_grid: 2D array of energy values for each FT in e_grid_ft
• e_grid_ft: 1D array of FT values for which the FWHM values in e_grid are calculated
• rt: RT value for which the FWHM values in e_grid are calculated
• min_e: minimum energy value to consider for the fit
• max_e: maximum energy value to consider for the fit
• nbins: number of bins to use for the histogram of energy values
• rel_cut_fit: relative cut value to use for the fit
• ft_fwhm_tol: search for lowest "optimal" ft within minimum(fwhm) + ft_fwhm_tol to avoid artificially large ft

Returns

• ft: optimal FT value
• min_fwhm: corresponding FWHM value

_fit_fwhm_ft_ctc(e_grid::Matrix, e_grid_ft::StepRangeLen, qdrift::Vector{<:Real}, rt::Unitful.RealOrRealQuantity, min_e::T, max_e::T, nbins::Int64, rel_cut_fit::T; default_ft::Quantity{T}=3.0u"µs", peak::Unitful.Energy{<:Real}=2614.5u"keV", window::Tuple{<:Unitful.Energy{<:Real}, <:Unitful.Energy{<:Real}}=(35.0u"keV", 25.0u"keV"), ft_fwhm_tol::Unitful.Energy{<:Real} = 0.1u"keV") where {T <:Real}

Fit the FWHM values in e_grid for each FT in e_grid_ft with a Gamma Peakshape and return the optimal FT and the corresponding FWHM value. The cut values cut for each flat-top time a window for better histogramming.

Arguments

• e_grid: 2D array of energy values for each FT in e_grid_ft
• e_grid_ft: 1D array of FT values for which the FWHM values in e_grid are calculated
• qdrift: drift time values for each energy value in e_grid
• rt: RT value for which the FWHM values in e_grid are calculated
• min_e: minimum energy value to consider for the fit
• max_e: maximum energy value to consider for the fit
• nbins: number of bins to use for the histogram of energy values
• rel_cut_fit: relative cut value to use for the fit
• ft_fwhm_tol: search for lowest "optimal" ft within minimum(fwhm) + ft_fwhm_tol to avoid artificially large ft

Returns

• ft: optimal FT value
• min_fwhm: corresponding FWHM value

_get_model_counts(f_fit::Base.Callable,v_ml::Union{NamedTuple, AbstractVector},bin_centers::StepRangeLen,bin_widths::StepRangeLen)

aux. function to get modelled peakshape based on histogram binning and best-fit parameter

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

aux. function to convert histogram data into bin edges, bin width and bin counts

advanced_time_and_memory_control(x::Optim.OptimizationState, start_time::Float64, time_to_setup::Float64; time_limit::Float64=60.0, mem_limit::Float64=30.0)

Control function to stop optimization based on time and memory usage.

Arguments

• x::Optim.OptimizationState: optimization state
• start_time::Float64: start time
• time_to_setup::Float64: time to setup
• time_limit::Float64: time limit
• mem_limit::Float64: memory limit

Return

• Bool: true if optimization should stop

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

Describes the background shape of a typical A/E Compton peak in a detector as a step like background for MSE events.

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

Describes the shape of a typical A/E Compton peak in a detector as a gaussian SSE peak and a step like background for MSE events.

aoe_compton_peakshape_components(fit_func::Symbol; background_center::Real)

This function defines the components (signal, low/high-energy tail, backgrounds) of the fit function used in gamma specfits. These component functions are used in the fit-report and in plot receipes

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

Describes the signal shape of a typical A/E Compton peak in a detector as a gaussian SSE peak.

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

Return a NamedTuple with the values of a in the order given by fieldnames(as_nt).

autocal_energy(e::AbstractArray{<:Real}, photon_lines::Vector{<:Unitful.RealOrRealQuantity}; min_e::Real=100, 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)

Compute an energy calibration from raw reconstructed energy deposition values based on a given number of known photon lines which are contained in the spectrum

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

)

Describes the background part of the shape of a typical gamma peak in a detector: components:

• step-function scaled with step_amplitude`
• energy-independent background: background`
• linear slope: background_slope (optional)
• exponential decay: background_exp (optional)

baseline_qc(data, qc_config)

Perform simple Gaussian fits on the baseline disitrbutions for a given data set.

Returns

• result: Namedtuple containing the cut and fit results for each baseline variable
• report: Namedtuple containing plotable objects.

fit_chisq(x::AbstractVector{<:Real},y::AbstractVector{<:Real},yerr::AbstractVector{<:Real}, f_fit::Function;pull_t::Vector{<:NamedTuple} = fill(NamedTuple(), first(methods(f_fit)).nargs - 2), v_init::Vector = [])

Least square fit with chi2 minimization

Input:

• x : x-values
• y : y-values
• yerr : 1 sigma uncertainty on y
• ffit : fit/model function. e.g. for a linear function: flin(x,p1,p2) = p1 .* x .+ p2

The numer of fit parameter is determined with first(methods(f_fit)).nargs - 2. That's why it's important that ffit has the synthax f(x,arg1,arg2,arg3,...) pullt : pull term, a vector of NamedTuple with fields mean and std. A Gaussian pull term is added to the chi2 function to account for systematic uncertainties. If left blank, no pull term is used. v_init : initial value for fit parameter optimization. If left blank, the initial value is set to 1 or guessed roughly for all fit parameters

Return:

• result : NamedTuple with the optimized fit parameter and the fit function
• report:

ctc_energy(e::Array{T}, qdrift::Array{T}, peak::T, window::T) where T<:Real

Correct for the drift time dependence of the energy by minimizing the ratio of the FWHM and the peak height of the peak around peak in e with a cut window of window. The drift time dependence is given by qdrift.

Returns

* `peak`: peak position
* `window`: window size
* `fct`: correction factor
* `fwhm_before`: FWHM before correction
* `fwhm_after`: FWHM after correction
* `func`: function to correct energy
* `func_generic`: generic function to correct energy

cut_single_peak(x::Array, min_x::Float64, max_x::Float64,; n_bins::Int=-1, relative_cut::Float64=0.5)

Cut out a single peak from the array x between min_x and max_x. The number of bins is the number of bins to use for the histogram. The relative cut is the fraction of the maximum counts to use for the cut.

Returns

* `max`: maximum position of the peak
* `low`: lower edge of the cut peak
* `high`: upper edge of the cut peak

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

Evaluates the sum of two gaussians at x with parameters μ1, σ1, n1, μ2, σ2, n2.

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

Estimate the peak position, FWHM, sigma, counts and background of a peak from a histogram.

estimate_fwhm(v::NamedTuple, v_err::NamedTuple)

Get the FWHM of a peak from the fit parameters.

Returns

* `fwhm`: the FWHM of the peak

estimate_single_peak_stats(h::Histogram, calib_type::Symbol=:th228)

Estimate statistics/parameters for a single peak in the given histogram h.

h must only contain a single peak. The peak should have a Gaussian-like shape. calib_type specifies the calibration type. Currently :th228, :psd and :simple is implemented..

Returns

NamedTuple with the fields * peak_pos: estimated position of the peak (in the middle of the peak) * 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 value

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

The PDF of an Exponentially modified Gaussian distribution with Gaussian parameters μ, σ and exponential scale θ at x.

It is the PDF of the distribution that descibes the random process rand(Normal(μ, σ)) + rand(Exponential(θ)).

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

Evaluates an extended step gauss model at x with parameters l, k, t and d.

expand_vars(v::NamedTuple)::StructArray

Expand all fields in v (scalars or arrays) to same array size and return a StructArray.

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

Evaluates an exponential decay function at x with parameters amplitude, decay and offset.

fit_aoe_compton(peakhists::Array, peakstats::StructArray, compton_bands::Array{T}) where T<:Real

Fit the A/E Compton bands using the f_aoe_compton function consisting of a gaussian SSE peak and a step like background for MSE events.

Returns

* `result`: Dict of NamedTuples of the fit results containing values and errors for each compton band
* `report`: Dict of NamedTuples of the fit report which can be plotted for each compton band

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

Performed a combined fit over all A/E Compton band using the f_aoe_compton function consisting of a gaussian SSE peak and a step like background for MSE events, assuming f_aoe_mu for μ and f_aoe_sigma for σ.

Returns

* `v_ml`: The fit result from the maximum-likelihood fit.
* `report_μ`: Report to plot the combined fit result for the enery-dependence of `μ`.
* `report_σ`: Report to plot the combined fit result for the enery-dependence of `σ`.
* `report_bands`: Dict of NamedTuples of the fit report which can be plotted for each compton band

fitaoecorrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Real}, σ::Array{<:Real})

Fit the corrections for the AoE value of the detector.

Returns

• e: Energy values
• μ: Mean values
• σ: Sigma values
• μ_scs: Fit result for the mean values
• f_μ_scs: Fit function for the mean values
• σ_scs: Fit result for the sigma values
• f_σ_scs: Fit function for the sigma values

fit_binned_double_gauss(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background, :μ, :σ), NTuple{7, T}}; uncertainty::Bool=true) where T<:Real

Perform a binned fit of the peakshape to the data in h using the initial values in ps while using the f_double_gauss function consisting of a double gaussian peak. The priors for the first gaussian peak are given by the ps tuple. For the priors of the second gaussian peak a wide window around the first peak is used.

Returns

* `result`: NamedTuple of the fit results containing values and errors
* `report`: NamedTuple of the fit report which can be plotted

fit_binned_trunc_gauss(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background, :μ, :σ), NTuple{7, T}}; uncertainty::Bool=true) where T<:Real

Perform a binned fit of the peakshape to the data in h using the initial values in ps while using the f_gauss function consisting of a gaussian peak multiplied with an amplitude n.

Returns

* `result`: NamedTuple of the fit results containing values and errors
* `report`: NamedTuple of the fit report which can be plotted

fit_calibration(pol_order::Int, µ::AbstractVector{<:Union{Real,Measurement{<:Real}}}, peaks::AbstractVector{<:Union{Real,Measurement{<:Real}}}; pull_t::Vector{<:NamedTuple}=fill(NamedTuple(), pol_order+1), v_init::Vector = [], uncertainty::Bool=true )

Fit the calibration lines with polynomial function of polorder order polorder == 1 -> linear function pol_order == 2 -> quadratic function

Returns

* `result`: NamedTuple with the following fields
    * `par`: best-fit parameters
    * `gof`: godness of fit
* `report`:

fit_enc_sigmas(enc_grid::Matrix{T}, enc_grid_rt::StepRangeLen{Quantity{<:T}, Base.TwicePrecision{Quantity{<:T}}, Base.TwicePrecision{Quantity{<:T}}, Int64}, min_enc::T, max_enc::T, nbins::Int64, rel_cut_fit::T) where T<:Real

Fit the ENC values in enc_grid for each RT in enc_grid_rt with a Gaussian and return the optimal RT and the corresponding ENC value.

Arguments

• enc_grid: 2D array of ENC values for each RT in enc_grid_rt
• enc_grid_rt: 1D array of RT values for which the ENC values in enc_grid are calculated
• min_enc: minimum ENC value to consider for the fit
• max_enc: maximum ENC value to consider for the fit
• nbins: number of bins to use for the histogram of ENC values
• rel_cut_fit: relative cut value to use for the fit

Returns

• rt: optimal RT value
• min_enc: corresponding ENC value

fitFWHM(fit_fwhm(peaks::Vector{T}, fwhm::Vector{T}) where T<:Real

Fit the FWHM of the peaks to a quadratic function.

Returns

* `qbb`: the FWHM at 2039 keV
* `err`: the uncertainties of the fit parameters
* `v`: the fit result parameters
* `f_fit`: the fitted function

fit_fwhm_ft(e_grid::Matrix, e_grid_ft::StepRangeLen, rt::Unitful.RealOrRealQuantity, min_e::T, max_e::T, rel_cut_fit::T, apply_ctc::Bool=true; kwargs...)
fit_fwhm_ft(e_grid::Matrix, e_grid_ft::StepRangeLen, rt::Unitful.RealOrRealQuantity, min_e, max_e, rel_cut_fit; kwargs...)
fit_fwhm_ft(e_grid::Matrix, e_grid_ft::StepRangeLen, qdrift::Vector{<:Real}, rt::Unitful.RealOrRealQuantity, min_e, max_e, rel_cut_fit; kwargs...)

Fit the FWHM values in e_grid for each FT in e_grid_ft with a Gamma Peakshape and return the optimal FT and the corresponding FWHM value. The cut values cut for each flat-top time a window for better histogramming. If the apply_ctc flag is set to true, the CTC correction is applied to the energy values. Othwise, if a qdrift vector is provided, the CTC correction is applied to the energy values.

Arguments

• e_grid: 2D array of energy values for each FT in e_grid_ft
• e_grid_ft: 1D array of FT values for which the FWHM values in e_grid are calculated
• rt: RT value for which the FWHM values in e_grid are calculated
• min_e: minimum energy value to consider for the fit
• max_e: maximum energy value to consider for the fit
• rel_cut_fit: relative cut value to use for the fit

Returns

• ft: optimal FT value
• min_fwhm: corresponding FWHM value

fit_half_centered_trunc_gauss(x::Array, cuts::NamedTuple{(:low, :high, :max), Tuple{Float64, Float64, Float64}})

Fit a single truncated Gaussian to the data x between cut.low and cut.high. The peak center is fixed at μ and the peak is cut in half either in the left or right half.

Returns report and result` with:

* `f_fit`: fitted function
* `μ`: mean of the Gaussian
* `σ`: standard deviation of the Gaussian

fit_half_trunc_gauss(x::Array, cuts::NamedTuple{(:low, :high, :max), Tuple{Float64, Float64, Float64}})

Fit a single truncated Gaussian to the data x between cut.low and cut.high. The peak center is fixed at μ and the peak is cut in half either in the left or right half.

Returns report and result with:

* `f_fit`: fitted function
* `μ`: mean of the Gaussian
* `σ`: standard deviation of the Gaussian

fit_lq_compton(peakhists::Array, peakstats::StructArray, compton_bands::Array{T}) where T<:Real

Fit the A/E Compton bands using the f_lq_compton function consisting of a gaussian SSE peak and a step like background for MSE events.

Returns

* `result`: Dict of NamedTuples of the fit results containing values and errors for each compton band
* `report`: Dict of NamedTuples of the fit report which can be plotted for each compton band

fit_peaks(peakhists::Array, peakstats::StructArray, th228_lines::Array,; calib_type::Symbol=:th228, uncertainty::Bool=true, low_e_tail::Bool=true)

Perform a fit of the peakshape to the data in peakhists using the initial values in peakstats to the calibration lines in th228_lines.

Returns

* `peak_fit_plots`: array of plots of the peak fits
* `return_vals`: dictionary of the fit results

fit_peaks_combined(peakhists::Array, peakstats::StructArray, th228_lines::Array{T},; calib_type::Symbol=:th228, uncertainty::Bool=true, fixed_position::Bool=false) where T<:Real

Fit the peaks of a histogram to a combined peakshape function while sharring parameters between peaks.

fit_sf_wl(dep_sep_data, a_grid_wl_sg, optimization_config)

Fit a A/E filter window length for the SEP data and return the optimal window length and the corresponding survival fraction.

Arguments

• dep_sep_data: NamedTuple with the DEP and SEP data
• a_grid_wl_sg: range of window lengths to sweep through
• optimization_config: configuration dictionary

Returns

• result: optimal window length and corresponding survival fraction
• report: report with all window lengths and survival fractions

fit_single_aoe_compton(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background, :μ, :σ), NTuple{7, T}}; uncertainty::Bool=true) where T<:Real

Perform a fit of the peakshape to the data in h using the initial values in ps while using the f_aoe_compton function consisting of a gaussian SSE peak and a step like background for MSE events.

Returns

* `result`: NamedTuple of the fit results containing values and errors
* `report`: NamedTuple of the fit report which can be plotted

fit_single_aoe_compton_with_fixed_μ_and_σ(h::Histogram, μ::Number, σ::Number, ps::NamedTuple; uncertainty::Bool=true)

Fit a single A/E Compton band using the f_aoe_compton function consisting of a gaussian SSE peak and a step like background for MSE events using fixed values for μ and σ.

Returns

* `neg_log_likelihood`: The negative log-likelihood of the likelihood fit
* `report`: Dict of NamedTuples of the fit report which can be plotted for each compton band

fit_single_lq_compton(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background, :μ, :σ), NTuple{7, T}}; uncertainty::Bool=true) where T<:Real

Perform a fit of the peakshape to the data in h using the initial values in ps while using the f_lq_compton function consisting of a gaussian SSE peak and a step like background for MSE events.

Returns

* `result`: NamedTuple of the fit results containing values and errors
* `report`: NamedTuple of the fit report which can be plotted

fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background), NTuple{5, T}};, uncertainty::Bool=true, fixed_position::Bool=false, low_e_tail::Bool=true) where T<:Real

Perform a fit of the peakshape to the data in h using the initial values in ps while using the gamma_peakshape with low-E tail. Also, FWHM is calculated from the fitted peakshape with MC error propagation. The peak position can be fixed to the value in ps by setting fixed_position=true. If the low-E tail should not be fitted, it can be disabled by setting low_e_tail=false.

Returns

* `result`: NamedTuple of the fit results containing values and errors
* `report`: NamedTuple of the fit report which can be plotted

fit_single_trunc_gauss(x::Array, cuts::NamedTuple{(:low, :high, :max), Tuple{Float64, Float64, Float64}})

Fit a single truncated Gaussian to the data x between min_x and max_x. Returns report and resultwith: *ffit: fitted function *μ: mean of the Gaussian *μerr: error of the mean *σ: standard deviation of the Gaussian *σ_err: error of the standard deviation *n`: number of counts in the peak

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=50, nInit::Int=50, 
    method::Symbol=:kmeans, kind=:diag, Δpe_peak_assignment::Real=0.3, f_uncal::Function=identity, uncertainty::Bool=true)

Fit a Gaussian Mixture Model to the given pe calibration data and return the fit parameters.

Arguments

• pe_cal::Vector{<:Real}: the pe calibration data
• min_pe::Real=0.5: the minimum pe to consider
• max_pe::Real=3.5: the maximum pe to consider
• n_mixtures::Int=ceil(Int, (max_pe - min_pe) * 4): the number of mixtures to fit
• nIter::Int=50: the number of iterations for the EM algorithm
• nInit::Int=50: the number of initializations for the EM algorithm
• method::Symbol=:kmeans: the method to use for initialization
• kind::Symbol=:diag: the kind of covariance matrix to use
• Δpe_peak_assignment::Real=0.3: the range to consider for peak assignment
• f_uncal::Function=identity: the function to use for uncalibration
• uncertainty::Bool=true: whether to calculate the uncertainty

Returns

• result: a tuple with the fit parameters
• report: a tuple with the fit report which can be plotted via a recipe

fit_single_peak_th228(h::Histogram, ps::NamedTuple{(:peak_pos, :peak_fwhm, :peak_sigma, :peak_counts, :mean_background), NTuple{5, T}};, uncertainty::Bool=true, fixed_position::Bool=false, low_e_tail::Bool=true) where T<:Real

Perform a simultaneous fit of two peaks (h_survived and h_cut) that together would form a histogram h, from which the result h_result was already determined using fit_single_peak_th228. Also, FWHM is calculated from the fitted peakshape with MC error propagation. The peak position can be fixed to the value in ps by setting fixed_position=true. If the low-E tail should not be fitted, it can be disabled by setting low_e_tail=false.

Returns

* `result`: NamedTuple of the fit results containing values and errors, in particular the signal survival fraction `sf` and the background survival frachtion `bsf`.
* `report`: NamedTuple of the fit report which can be plotted

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... 
)

Standard gamma peakshape: Describes the shape of a typical gamma peak in a detector. Components:

• Gaussian signal peak with μ, σ, n - skew_fraction - skew_fraction_highE
• low-energy tail: skew_fraction, skew_width
• high-energy tail: skew_fraction_highE, skew_width_highE (optional, default off)
• background:
  • energy-independent background
  • step-function scaled with step_amplitude from Compton scattered gammas
  • linear slope: background_slope (optional, default off)
  • exponential decay: background_exp (optional, default off)

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

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}

Generate histograms for the A/E Compton bands and estimate peak parameters. The compton bands are cutted out of the A/E spectrum and then binned using the Freedman-Diaconis Rule. For better performance the binning is only done in the area around the peak. The peak parameters are estimated using the estimate_single_peak_stats_psd function.

Returns

* `peakhists`: Array of histograms for each compton band
* `peakstats`: StructArray of peak parameters for each compton band
* `min_aoe`: Array of minimum A/E values for each compton band
* `max_aoe`: Array of maximum A/E values for each compton band
* `mean_peak_pos`: Mean peak position of all compton bands
* `std_peak_pos`: Standard deviation of the peak position of all compton bands
* `simple_pars_aoe_μ`: Simple curve fit parameters for the peak position energy depencence
* `simple_pars_error_aoe_μ`: Simple curve fit parameter errors for the peak position energy depencence
* `simple_pars_aoe_σ`: Simple curve fit parameters for the peak sigma energy depencence
* `simple_pars_error_aoe_σ`: Simple curve fit parameter errors for the peak sigma energy depencence

get_centered_gaussian_window_cut(x::Array, min_x::Float64, max_x::Float64, n_σ::Real, center::Float64=0.0, n_bins_cut::Int=500, relative_cut::Float64=0.2, left::Bool=false)

Cut out a single peak from the array x between min_x and max_x by fitting a truncated one-sided Gaussian and extrapolating a window cut with n_σ standard deviations. The center and side of the fit can be specified with left and center variable.

Returns

* `low_cut`: lower edge of the cut peak
* `high_cut`: upper edge of the cut peak
* `center`: center of the peak
* `σ`: standard deviation of the Gaussian
* `low_cut_fit`: lower edge of the cut peak from the fit
* `high_cut_fit`: upper edge of the cut peak from the fit
* `err`: error of the fit parameters

get_continuum_survival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, center::T, window::T, aoe_cut::Unitful.RealOrRealQuantity,; inverted_mode::Bool=false, sigma_high_sided::Unitful.RealOrRealQuantity=Inf*unit(first(aoe))) where T<:Unitful.Energy{<:Real}

Get the survival fraction of a continuum after a AoE cut value aoe_cut for a given center and window size.

Returns

• center: Center of the continuum
• window: Window size
• n_before: Number of counts before the cut
• n_after: Number of counts after the cut
• sf: Survival fraction

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

Return a DistributionTransform for the given Distribution and Prior.

get_friedman_diaconis_bin_width(x::AbstractArray)

Return the bin width for the given data x using the Friedman-Diaconis rule.

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{<:Unitful.RealOrRealQuantity, <:Unitful.RealOrRealQuantity}=(-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}

Get the AoE cut value for a given dep and window size while performing a peak fit with fixed position. The AoE cut value is determined by finding the cut value for which the number of counts after the cut is equal to dep_sf times the number of counts before the cut. The algorhithm utilizes a root search algorithm to find the cut value with a relative tolerance of rtol.

Returns

• cut: AoE cut value
• n0: Number of counts before the cut
• nsf: Number of counts after the cut

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

Return a NamedTuple with the same fields as v and v_err but with Normal distributions for each field.

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

Generate n random samples of fit parameters using their respective best-fit values v and covariance matrix v_err

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

Return the number of bins for the given data x using the given method.

get_peak_fwhm_th228(v_ml::NamedTuple, v_ml_err::NamedTuple)

Get the FWHM of a peak from the fit parameters while performing a MC error propagation.

Returns

* `fwhm`: the FWHM of the peak
* `fwhm_err`: the uncertainty of the FWHM of the peak

get_peak_survival_fraction(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peak::T, window::Vector{T}, aoe_cut::Unitful.RealOrRealQuantity,; 
uncertainty::Bool=true, inverted_mode::Bool=false, bin_width_window::T=2.0u"keV", sigma_high_sided::Unitful.RealOrRealQuantity=Inf*unit(first(aoe)), fit_func::Symbol=:gamma_def) where T<:Unitful.Energy{<:Real}

Get the survival fraction of a peak after a AoE cut value aoe_cut for a given peak and window size while performing a peak fit with fixed position.

Returns

• peak: Peak position
• n_before: Number of counts before the cut
• n_after: Number of counts after the cut
• sf: Survival fraction
• err: Uncertainties

get_peakhists_th228(e::Array, th228_lines::Array, window_sizes::Array, e_unit::String="keV", proxy_binning_peak::Float64=2103.5, proxy_binning_peak_window::Float64=10.0)

Create histograms around the calibration lines and return the histograms and the peak statistics.

Returns

* `peakhists`: array of histograms around the calibration lines
* `peakstats`: array of statistics for the calibration line fits

get_peaks_survival_fractions(aoe::Vector{<:Unitful.RealOrRealQuantity}, e::Vector{<:T}, peaks::Vector{<:T}, peak_names::Vector{Symbol}, windows::Vector{<:Tuple{T, T}}, aoe_cut::Unitful.RealOrRealQuantity,; uncertainty::Bool=true, inverted_mode::Bool=false, bin_width_window::T=2.0u"keV", sigma_high_sided::Unitful.RealOrRealQuantity=Inf*unit(first(aoe)), fit_funcs::Vector{Symbol}=fill(:gamma_def, length(peaks))) where T<:Unitful.Energy{<:Real}

Get the survival fraction of a peak after a AoE cut value aoe_cut for a given peak and window size while performing a peak fit with fixed position.

Return

• result: Dict of results for each peak
• report: Dict of reports for each peak

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

Calculate bin-wise residuals and normalized residuals. Calcualte bin-wise p-value based on poisson distribution for each bin.

Input:

• f_fitfunction handle of fit function (peakshape)
• h histogram of data
• v_ml best-fit parameters

Returns:

• residuals difference: model - data (histogram bin count)
• residuals_norm normalized residuals: model - data / sqrt(model)
• p_value_binwise p-value for each bin based on poisson distribution

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}

Get the survival fraction after a AoE cut value aoe_cut for a given peak and window size from a combined fit to the survived and cut histograms.

Returns

• sf: Survival fraction after the cut

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

This function defines the gamma peakshape fit functions used in the calibration specfits.

• gamma_def: "default" gamma peakshape with gaussian signal, low-energy tail, and background (flat + step)
• gamma_tails: default gamma peakshape + high-energy tail
• gamma_bckSlope: default gamma peakshape + linear background slope
• gamma_bckExp: default gamma peakshape + exponential background
• gamma_bckFlat: default gamma peakshape - step background (only flat component!)
• gammatailsbckFlat: default gamma peakshape + high-energy tail - step background (only flat component!)

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

Describes the high-E signal tail part of the shape

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

Calculate the Poisson log-likelihood of a fit function f_fit(x) and a histogram h. f_fit must accept all values x on the horizontal axis of the histogram.

Currently uses a simple midpoint-rule integration of f_fit over the bins of h.

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

Evaluates a linear function at x with parameters slope and intercept.

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

Describes the low-E signal tail part of the shape of a typical gamma peak in a detector.

lq_norm::Vector{Float64}, dt_eff::Vector{<:Unitful.RealOrRealQuantity}, e_cal::Vector{<:Unitful.Energy{<:Real}}, DEP_µ::Unitful.AbstractQuantity, DEP_σ::Unitful.AbstractQuantity; 
DEP_edgesigma::Float64=3.0 , mode::Symbol=:percentile, drift_cutoff_sigma::Float64 = 2.0, prehist_sigma::Float64=2.5, e_expression::Union{String,Symbol}="e", dt_eff_low_quantile::Float64=0.15, dt_eff_high_quantile::Float64=0.95)

Perform the drift time correction on the LQ data using the DEP peak. The function cuts outliers in lq and drift time, then performs a linear fit on the remaining data. The data is Corrected by subtracting the linear fit from the lq data.

Returns

* `result`: NamedTuple of the function used for lq classifier construction
* `report`: NamedTuple of the histograms used for the fit, the cutoff values and the DEP edges

nearestSPD(A::Matrix{<:Real})

Returns the nearest positive definite matrix to A Calculation is based on matrix factorization techniques described in https://www.sciencedirect.com/science/article/pii/0024379588902236

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

calculate p-value based on least-squares, assuming gaussian uncertainty baseline method to get goodness-of-fit (gof)

input:

• f_fitfunction handle of fit function (peakshape)
• h histogram of data
• v_ml best-fit parameters

returns:

• pval p-value of chi2 test
• chi2 chi2 value
• dof degrees of freedom

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

alternative p-value calculation via Monte Carlo sampling. Warning: computational more expensive than pvaule() and pvalue_LogLikeRatio()

Input:

• f_fitfunction handle of fit function (peakshape)
• h histogram of data
• ps best-fit parameters
• v_ml best-fit parameters
• n_samples number of samples

Performed Steps:

• Create n_samples randomized histograms. For each bin, samples are drawn from a Poisson distribution with λ = model peak shape (best-fit parameter)
• Each sample histogram is fit using the model function f_fit
• For each sample fit, the max. loglikelihood fit is calculated

Returns

• % p value –> comparison of sample max. loglikelihood and max. loglikelihood of best-fit

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

p-value via poisson likelihood ratio: baseline for ML fits using Poisson statistics and bins with low number of counts

peak_centroid(v::NamedTuple)

calculate centroid of gamma peak from fit parameters

peakshape_components(fit_func::Symbol; background_center::Real)

This function defines the components (signal, low/high-energy tail, backgrounds) of the fit function used in gamma specfits. These component functions are used in the fit-report and in plot receipes

prepare_sep_peakhist(e::Array{T}, dep::T,; relative_cut::T=0.5, n_bins_cut::Int=500) where T<:Real

Prepare an array of uncalibrated SEP energies for parameter extraction and calibration.

Returns

• result: Result of the initial fit
• report: Report of the initial fit

pulser_cal_qc(data, pulser_config; n_pulser_identified=100)

Perform simple QC cuts on the data and return the data for energy calibration.

Returns

- pulser_idx: indices of the pulser events

set_memlimit(gig::Float64)
Set memory limit in GB for the whole current process ignoring if things could run in parallel

set_timelimit(sec::Float64)
Set time limit in seconds a single Optim.jl optimization step

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

Describes the signal part of the shape of a typical gamma peak in a detector.

simple_calibration(e_uncal::Array, th228_lines::Array, window_size::Float64=25.0, n_bins::Int=15000, calib_type::String="th228")

Perform a simple calibration for the uncalibrated energy array e_uncal using the calibration type calib_type and the calibration lines th228_lines. The window size is the size of the window around the calibration line to use for the calibration. The number of bins is the number of bins to use for the histogram.

Returns * h_calsimple: histogram of the calibrated energy array * h_uncal: histogram of the uncalibrated energy array * c: calibration factor * fep_guess: estimated full energy peak (FEP) * peakhists: array of histograms around the calibration lines * peakstats: array of statistics for the calibration line fits

sipm_simple_calibration(pe_uncal::Array)

Perform a simple calibration for the uncalibrated p.e. spectrum array pe_uncal using just the 1 p.e. and 2 p.e. peak positions estimated by a peakfinder.

Inputs: * pe_uncal: array of uncalibrated peak amplitudes kwargs: * initial_min_amp: uncalibrated amplitude value as a left boundary to build the uncalibrated histogram where the peak search is performed on. For the peak search with noise peak, this value is consecutively increased i.o.t exclude the noise peak from the histogram. * initial_max_quantile: quantile of the uncalibrated amplitude array to used as right boundary to build the uncalibrated histogram * peakfinder_σ: sigma value in number of bins for peakfinder * peakfinder_threshold: threshold value for peakfinder

Returns * pe_simple_cal: array of the calibrated pe array with the simple calibration * func: function to use for the calibration (pe_simple_cal = pe_uncal .* c .+ offset) * c: calibration factor * offset: calibration offset * peakpos: 1 p.e. and 2 p.e. peak positions in uncalibrated amplitude * peakpos_cal: 1 p.e. and 2 p.e. peak positions in calibrated amplitude * h_uncal: histogram of the uncalibrated pe array * h_calsimple: histogram of the calibrated pe array with the simple calibration

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

Evaluates the convulution of a Heaviside step function and the PDF of Normal(μ, σ) at x.

The result does not correspond to a PDF as it is not normalizable.

subhist(h::Histogram, r::Tuple{<:Real,<:Real})

Return a new Histogram with the bins in the range r.

tuple_to_array(nt::NamedTuple, fields::Vector{Symbol})

Return an array with the values of the fields in nt in the order given by fields.

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

Construct a Weibull distribution with a given median m and a given p_x-quantile x.

Useful to construct priors for positive quantities.