Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updated docstrings #80

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 14 additions & 0 deletions docs/src/index.md
Original file line number Diff line number Diff line change
@@ -1 +1,15 @@
# LegendSpecFits.jl

This package provides functionality to fit probability distribution functions to histograms.

This package is part of the [LegendJuliaRegistry](https://github.com/legend-exp/LegendJuliaRegistry), meaning that you can download the package as follows:
```julia
import Pkg

# add the LegendJuliaRegistry
Pkg.Registry.add(Pkg.RegistrySpec(url = "https://github.com/legend-exp/LegendJuliaRegistry"))

# add LegendSpecFits
Pkg.add("LegendSpecFits")
```

27 changes: 23 additions & 4 deletions src/aoe_ctc.jl
Original file line number Diff line number Diff line change
Expand Up @@ -5,18 +5,37 @@
"""
## TO DO: Still needs to be updated once finalized

ctc_aoe(aoe_all::Array{T}, ecal_all::Array{T}, qdrift_e_all::Array{T}, compton_bands::Array{T}, peak::T, window::T) where T<:Real
ctc_aoe(aoe_all::Vector{<:Real}, ecal_all::Vector{<:Unitful.RealOrRealQuantity}, qdrift_e_all::Vector{<:Real}, compton_bands::Vector{<:Unitful.RealOrRealQuantity}, 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; aoe_expression::Union{Symbol, String}="aoe", qdrift_expression::Union{Symbol, String} = "qdrift / e", pseudo_prior::NamedTupleDist = NamedTupleDist(empty = true), pseudo_prior_all::NamedTupleDist = NamedTupleDist(empty = true), pol_order::Int=1)

Correct for the drift time dependence of A/E parameter

# Returns (but may change):

# Arguments
* `aoe_all`: All A/E values
* `ecal_all`: All calibrated energies
* `qdrift_e_all`: All charge drift energy
* `compton_bands`: Compton bands
* `peak`: Peak position
* `window`: Window size
* `hist_start`: Histogram start position
* `hist_end`: Histogram end position
* `bin_width`: Bin width

# Keywords
* `aoe_expression`: A/E expression
* `qdrift_expression`: Charge drift expression
* `pseudo_prior`: Pseudo prior
* `peudo_prior_all`: All pseudo priors
* `pol_order`: Polynomial order

# Returns
* `peak`: peak position
* `window`: window size
* `func`: function to correct aoe
* `fct`: correction factor
* `σ_before`: σ before correction
* `σ_after`: σ after correction
* `func`: function to correct aoe
* `func_generic`: generic function to correct aoe
* `converged`: Convergence status of the fit
"""

function ctc_aoe(aoe_all::Vector{<:Real}, ecal_all::Vector{<:Unitful.RealOrRealQuantity}, qdrift_e_all::Vector{<:Real}, compton_bands::Vector{<:Unitful.RealOrRealQuantity}, 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; aoe_expression::Union{Symbol, String}="aoe", qdrift_expression::Union{Symbol, String} = "qdrift / e", pseudo_prior::NamedTupleDist = NamedTupleDist(empty = true), pseudo_prior_all::NamedTupleDist = NamedTupleDist(empty = true), pol_order::Int=1) # deleted m_cal since no calibration
Expand Down
117 changes: 98 additions & 19 deletions src/aoe_cut.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,21 @@

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.

# Arguments
* `aoe_cut`: A/E cut value
* `aoe`: A/E
* `e`: Calibrated energies
* `peak`: Peak position
* `window`: Data window in energy
* `bin_width`: Histogram bin widths
* `result_before`: Survival fraction before A/E cut

# Keywords
* `uncertainty`: Uncertainty of survival fraction
* `fit_func`: Fit function

# Returns
- `sf`: Survival fraction after the cut
* `result.sf`: Survival fraction after the cut
"""
function 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 energy after cut and create histogram
Expand All @@ -25,12 +38,35 @@
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`.

# Arguments
* `aoe`: A/E
* `e`: Calibrated energies

# Keywords
* `dep`: DEP energies
* `window`: Data window in energy
* `dep_sf`: Survival fraction of DEP energies
* `rtol`: Relative tolerance
* `maxiters`: Maximum iterations
* `sigma_high_sided`: Fixed upper limit cut
* `cut_search_interval`: Data interval to search for a cut value
* `bin_width_window`: Specified window around the peak in which the bin algorithm is applied
* `max_e_plot`: Maximum energy plot
* `plot_window`: Plot window
* `fixed_position`: Fixed position of cut
* `fit_func`: Fitted function
* `uncertainty`: Uncertainty

# Returns
- `cut`: AoE cut value
- `n0`: Number of counts before the cut
- `nsf`: Number of counts after the cut
* `lowcut`: low AoE cut value
* `highcut`: high AoE cut value
* `n0`: Number of counts before the cut
* `nsf`: Number of counts after the cut
* `sf`: Survival fraction
"""
function 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,
Expand Down Expand Up @@ -91,9 +127,25 @@
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.

# Arguments
* `aoe`: A/E
* `e`: Calibrated energies
* `peaks`: Data range of several peaks
* `peak_names`: Name of the peak
* `windows`: Energy data window
* `aoe_cut`: A/E cut value

# Keywords
* `uncertainty`: Uncertainty
* `inverted_mode`: If set to false: A/E analysis. If True: LQ analysis
* `bin_width_window`: Specified window around the peak where the binning algorithm is applied to
* `sigma_high_sided`: Fixed upper limit cut
* `fit_funcs`: Fitted functions

# Return
- `result`: Dict of results for each peak
- `report`: Dict of reports for each peak
* `result`: Dict of results for each peak
* `report`: Dict of reports for each peak
"""
function 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}
@assert length(peaks) == length(peak_names) == length(windows) "Length of peaks, peak_names and windows must be equal"
Expand Down Expand Up @@ -125,16 +177,31 @@

"""
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}
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.

# Arguments
* `aoe`: A/E
* `e`: Calibrated energies
* `peak`: Peak position
* `window`: Data window in energy
* `aoe_cut`: A/E cut value

# Keywords
* `Uncertainty`: Uncertainty
* `lq_mode`: Inverts the cut logic
* `low_e_tail`: Low energy tail
* `bin_width_window`: Specified window around the peak to apply the binning algorithm
* `sigma_high_sided`: Fixed upper limit cut

# 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
* `peak`: Peak position
* `fit_func`: fitted function
* `n_before`: Number of counts before the cut
* `n_after`: Number of counts after the cut
* `sf`: Survival fraction
* `gof`: Goodness of fit
"""
function 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}
Expand Down Expand Up @@ -189,17 +256,29 @@


"""
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_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.

# Arguments
* `aoe`: A/E
* `e`: Calibrated energies
* `center`: Center of the fit
* `window`: Data window in energy
* `aoe_cut`: A/E cut value

# Keywords
* `inverted_mode`: Inverted cut logic
* `sigma_high_sided`: Fixed upper limit cut

# 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
* `window`: Window size
* `n_before`: Number of counts before the cut
* `n_after`: Number of counts after the cut
* `sf`: Survival fraction

"""
function 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}
function 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}

Check warning on line 281 in src/aoe_cut.jl

View check run for this annotation

Codecov / codecov/patch

src/aoe_cut.jl#L281

Added line #L281 was not covered by tests
# scale unit
e_unit = u"keV"
# get energy around center
Expand Down
58 changes: 48 additions & 10 deletions src/aoe_filter_optimization.jl
Original file line number Diff line number Diff line change
@@ -1,12 +1,24 @@


"""
prepare_sep_peakhist(e::Array{T}, dep::T,; relative_cut::T=0.5, n_bins_cut::Int=500) where T<:Real
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

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

# Arguments
* `e`: Uncalibrated energies

# Keywords
* `sep`: Single escape peak
* `window`: Energy window
* `relative_cut`: Relative cut
* `n_bins_cut`: Number of cut bins
* `fit_func`: fitted function
* `uncertainty`: Uncertainty
:
# Returns
- `result`: Result of the initial fit
- `report`: Report of the initial fit
* `result`: Result of the initial fit
* `report`: Report of the initial fit
"""
function 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
# get cut window around peak
Expand All @@ -27,18 +39,44 @@ function prepare_sep_peakhist(e::Vector{<:T}; sep::Unitful.Energy{<:Real}=2103.5
end

"""
fit_sf_wl(dep_sep_data, a_grid_wl_sg, optimization_config)
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::T=1592.53u"keV", dep_window::Vector{<:T}=[12.0, 10.0]u"keV",
sep::T=2103.53u"keV", sep_window::Vector{<:T}=[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::Quantity{<:Real}=0.0u"keV^-1", max_aoe_offset::Quantity{<:Real}=0.05u"keV^-1",
dep_cut_search_fit_func::Symbol=:gamma_def, sep_cut_search_fit_func::Symbol=:gamma_def,
uncertainty::Bool = false) where T<:Unitful.Energy{<:Real}

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
* `e_dep`: DEP energy
* `aoe_dep`: Double escape peak A/E
* `e_sep`: SEP energy
* `aoe_sep`: Single escape peak A/E
* `a_grid_wl_sg`: Range of window lengths to sweep through

# Keywords
* `dep`: Double escape peak
* `dep_window`: DEP window
* `sep`: Single escape peak
* `sep_window`: SEP window
* `sep_rel_cut`: SEP relative cut
* `min_aoe_quantile`: Minimum A/E quantile
* `max_aoe_quantile`: Maximum A/E quantile
* `min_aoe_offset`: Minimum A/E offset
* `max_aoe_offset`: Maximum A/E offset
* `dep_cut_search_fit_func`: Symbol describing the function for the DEP cut search
* `sep_cut_search_fit_func`: Symbol describing the function for the SEP cut search
* `uncertainty`: Uncertainty


# Returns
- `result`: optimal window length and corresponding survival fraction
- `report`: report with all window lengths and survival fractions
* `wl`: Window length
* `sf`: Survival fraction
* `n_dep`: Number of DEP energies
* `n_sep`: Number of SEP energies

"""
function 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::T=1592.53u"keV", dep_window::Vector{<:T}=[12.0, 10.0]u"keV",
Expand All @@ -62,7 +100,7 @@ function fit_sf_wl(e_dep::Vector{<:Real}, aoe_dep::ArrayOfSimilarArrays{<:Real},
sep_sfs = Vector{Quantity}(undef, length(a_grid_wl_sg))
wls = Vector{eltype(a_grid_wl_sg)}(undef, length(a_grid_wl_sg))

# for each window lenght, calculate the survival fraction in the SEP
# for each window length, calculate the survival fraction in the SEP
Threads.@threads for i_aoe in eachindex(a_grid_wl_sg)
# get window length
wl = a_grid_wl_sg[i_aoe]
Expand Down
24 changes: 16 additions & 8 deletions src/aoe_fit_calibration.jl
Original file line number Diff line number Diff line change
@@ -1,17 +1,25 @@
@. f_aoe_sigma(x, p) = sqrt(p[1]^2 + p[2]^2/x^2)
f_aoe_mu(x, p) = p[1] .+ p[2].*x
"""
fit_aoe_corrections(e::Array{<:Unitful.Energy{<:Real}}, μ::Array{<:Real}, σ::Array{<:Real})
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")

Fit the corrections for the AoE value of the detector.

# Arguments
* `e`: Calibrated energies
* `μ`: Mean values
* `σ`: Sigma values

# Keywords
* `aoe_expression`: A/E expression
* `e_expression`: Calibrated energy expression

# 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
* `µ_compton`: Mean compton values
* `σ_compton`: Sigma compton values
* `compton_bands`: Compton bands
* `func`: A/E correction function

"""
function 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")
# fit compton band mus with linear function
Expand Down
34 changes: 34 additions & 0 deletions src/aoe_pseudo_prior.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,20 @@
"""
get_aoe_standard_pseudo_prior(h::Histogram, ps::NamedTuple, fit_func::Symbol; fixed_position::Bool=false)

Gets the A/E of a histogram using a standard pseudo prior.

# Arguments
* `h`: Histogram data
* `ps`: Peak statistics
* `fit_func`: Fit function

# Keywords
* `fixed_position`: Determines whether the function uses a fixed position or not

# Returns
* `pprior_base`: Base pseudo prior
"""

function get_aoe_standard_pseudo_prior(h::Histogram, ps::NamedTuple, fit_func::Symbol; fixed_position::Bool=false)
pprior_base = NamedTupleDist(
μ = ifelse(fixed_position, ConstValueDist(ps.peak_pos), Normal(ps.peak_pos, 0.5*ps.peak_sigma)),
Expand Down Expand Up @@ -25,6 +42,23 @@ function get_aoe_standard_pseudo_prior(h::Histogram, ps::NamedTuple, fit_func::S
end
end

"""
get_aoe_pseudo_prior(h::Histogram, ps::NamedTuple, fit_func::Symbol; pseudo_prior::NamedTupleDist=NamedTupleDist(empty = true), kwargs...)

Gets the A/E of a histogram using a pseudo prior.
# Arguments
* `h`: Histogram data
* `ps`: Peak statistics
* `fit_func`: Fit function

# Keywords
* `pseudo_prior`: Initial guess for parameters of histogram

# Returns
* `pseudo_prior`: Initial guess for parameters of histogram

"""

function get_aoe_pseudo_prior(h::Histogram, ps::NamedTuple, fit_func::Symbol; pseudo_prior::NamedTupleDist=NamedTupleDist(empty = true), kwargs...)
standard_pseudo_prior = get_aoe_standard_pseudo_prior(h, ps, fit_func; kwargs...)
# use standard priors in case of no overwrites given
Expand Down
Loading
Loading