Initial implementation.

Authors: Lukas Hauertmann <[email protected]> Oliver Schulz <[email protected]>
JuliaPhysics · Dec 1, 2018 · 7554009 · 7554009
1 parent 7a71058
commit 7554009
Show file tree

Hide file tree

Showing 17 changed files with 4,984 additions and 7 deletions.
diff --git a/LICENSE.md b/LICENSE.md
@@ -1,7 +1,9 @@
-The RadiationSpectra.jl package is licensed under the MIT "Expat" License:
+The RadiationSpectra.jl package is licensed under the MIT "Expat" License
+(with the exception of the function `peakfinder`, see below):
 
 > Copyright (c) 2018:
 >
+>    Lukas Hauertmann <[email protected]>
 >    Oliver Schulz <[email protected]>
 > 
 > Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -21,4 +23,24 @@ The RadiationSpectra.jl package is licensed under the MIT "Expat" License:
 > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 > SOFTWARE.
+
+
+The function `peakfinder` is licensed under LGPL 2.1, like the original
+C++ version, part of the CERN ROOT software:
+
+> Original Copyright (c) 1999: Miroslav Morhac
+> Translation to Julia 2018: Lukas Hauertmann <[email protected]>
+> 
+> This library is free software; you can redistribute it and/or
+> modify it under the terms of the GNU Lesser General Public
+> License as published by the Free Software Foundation; either
+> version 2.1 of the License, or (at your option) any later version.
+> 
+> This library is distributed in the hope that it will be useful,
+> but WITHOUT ANY WARRANTY; without even the implied warranty of
+> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+> Lesser General Public License for more details.
 > 
+> You should have received a copy of the GNU Lesser General Public
+> License along with this library; if not, write to the Free Software
+> Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,15 @@
 name = "RadiationSpectra"
-uuid = "52a91a86-ee3f-11e8-1f8c-f56e6b0fe16a"
-authors = ["Oliver Schulz <[email protected]>"]
-version = "0.1.0"
+uuid = "4f207c7e-01da-51d7-a1a0-c8c06dd1d883"
+
+[deps]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

diff --git a/REQUIRE b/REQUIRE
diff --git a/docs/make.jl b/docs/make.jl
@@ -8,11 +8,16 @@ using Documenter
 using RadiationSpectra
 
 makedocs(
-    sitename = "RadiationSpectra",
+    sitename = "RadiationSpectra.jl",
     modules = [RadiationSpectra],
     format = :html,
     pages=[
         "Home" => "index.md",
+        "Manual" => Any[
+            "Peak Finding" => "man/spectrum_deconvolution.md",
+            "Fitting" => "man/fitting.md",
+            "Calibration" => "man/spectrum_calibration.md",
+        ],
         "API" => "api.md",
         "LICENSE" => "LICENSE.md",
     ],

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -1 +1,17 @@
 # RadiationSpectra.jl
+
+RadiationSpectra.jl provides tools to analyse Radiation Spectra.
+
+This includes (for now):
+
+```@contents
+Pages = ["man/spectrum_deconvolution.md"]
+```
+```@contents
+Pages = ["man/fitting.md"]
+```
+```@contents
+Pages = ["man/spectrum_calibration.md"]
+```
+
+
diff --git a/docs/src/man/fitting.md b/docs/src/man/fitting.md
@@ -0,0 +1,128 @@
+# Fitting
+
+```@contents
+Pages = ["fitting.md"]
+```
+
+## `Likelihood (LLH) Fit - 1D-Histogram`
+
+1. Get a spectrum and find a peak to fit
+```@example fitting_hist
+using Plots, RadiationSpectra, StatsBase
+myfont = Plots.font(12) # hide
+pyplot(guidefont=myfont, xtickfont=myfont, ytickfont=myfont, legendfont=myfont, titlefont=myfont) # hide
+
+h_uncal = RadiationSpectra.get_example_spectrum()
+h_decon, peakpos = RadiationSpectra.peakfinder(h_uncal)
+
+strongest_peak_bin_idx = StatsBase.binindex(h_uncal, peakpos[1])
+strongest_peak_bin_width = StatsBase.binvolume(h_uncal, strongest_peak_bin_idx)
+strongest_peak_bin_amplitude = h_uncal.weights[strongest_peak_bin_idx]
+
+plot(h_uncal, st=:step, xlims=[peakpos[1] - strongest_peak_bin_width * 20, peakpos[1] + strongest_peak_bin_width * 20], size=(800,400), ylims=[0, strongest_peak_bin_amplitude * 1.1]) 
+```
+
+2. Write a model function
+```@example fitting_hist
+function model(x::T, par::Vector{T}) where {T}
+    scale::T = par[1]
+    σ::T     = par[2]
+    μ::T     = par[3]
+    cp0::T   = par[4] 
+    cp1::T   = par[5]
+    return scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + cp0 + cp1 * (x - μ)
+end
+function model(x::Vector{T}, par::Vector{T}) where {T} return T[model(v, par) for v in x] end;
+```
+
+3. Set up the fit function [`RadiationSpectra.FitFunction{T}`](@ref)
+```@example fitting_hist
+fitfunc = RadiationSpectra.FitFunction( model ); 
+fitfunc.fitrange = (peakpos[1] - 1000, peakpos[1] + 1000)
+guess_σ = strongest_peak_bin_width * 2
+guess_amplitude = strongest_peak_bin_amplitude * guess_σ
+guess_μ = peakpos[1]
+guess_offset = 0
+guess_bg_slope = 0
+fitfunc.initial_parameters = [ guess_amplitude, guess_σ, guess_μ, guess_offset, guess_bg_slope ]
+fitfunc
+```
+
+4. Performe the fit with the [`RadiationSpectra.llhfit!`](@ref)-function and plot the result
+```@example fitting_hist
+RadiationSpectra.llhfit!(fitfunc, h_uncal)
+
+plot(h_uncal, st=:step, xlims=[peakpos[1] - strongest_peak_bin_width * 20, peakpos[1] + strongest_peak_bin_width * 20], size=(800,400), label="Spectrum", ylims=[0, strongest_peak_bin_amplitude * 1.1])
+plot!(fitfunc, use_initial_parameters=true, lc=:green, label="Guess")
+plot!(fitfunc, lc=:red, label="LLH Fit")
+```
+
+```@example fitting_hist
+fitfunc # hide
+```
+
+## `LSQ Fit - 1D-Histogram`
+
+To perfrom a LSQ Fit on a spectrum repeat the first three steps from the [Likelihood (LLH) Fit - 1D-Histogram](@ref).
+Then, 
+
+4. Performe the fit with the [`RadiationSpectra.lsqfit!`](@ref)-function and plot the result
+```@example fitting_hist
+RadiationSpectra.lsqfit!(fitfunc, h_uncal)
+
+plot(h_uncal, st=:step, xlims=[peakpos[1] - strongest_peak_bin_width * 20, peakpos[1] + strongest_peak_bin_width * 20], size=(800,400), label="Spectrum", ylims=[0, strongest_peak_bin_amplitude * 1.1])
+plot!(fitfunc, use_initial_parameters=true, lc=:green, label="Guess")
+plot!(fitfunc, lc=:red, label="LSQ Fit")
+```
+
+```@example fitting_hist
+fitfunc # hide
+```
+
+## `LSQ Fit - 1D-Data Arrays`
+
+* Write a model function
+```@example fitting_1D_data
+using Plots, RadiationSpectra
+myfont = Plots.font(12) # hide
+pyplot(guidefont=myfont, xtickfont=myfont, ytickfont=myfont, legendfont=myfont, titlefont=myfont) # hide
+function model(x::T, par::Vector{T}) where {T}
+    cp0::T   = par[1] 
+    cp1::T   = par[2]
+    return cp0 + cp1 * x
+end
+function model(x::Vector{T}, par::Vector{T}) where {T} return T[model(v, par) for v in x] end;
+```
+
+* Create some random data 
+```@example fitting_1D_data
+xdata = Float64[1, 2, 3, 6, 8, 12]
+ydata = Float64[model(x, [-0.2, 0.7]) + (rand() - 0.5) for x in xdata]
+xdata_err = Float64[0.5 for i in eachindex(xdata)]
+ydata_err = Float64[1 for i in eachindex(xdata)]
+
+plot(xdata, ydata, xerr=xdata_err, yerr=ydata_err, st=:scatter, size=(800,400), label="Data")
+```
+
+* Set up the fit function [`RadiationSpectra.FitFunction{T}`](@ref)
+```@example fitting_1D_data
+fitfunc = RadiationSpectra.FitFunction( model ); 
+fitfunc.fitrange = (xdata[1], xdata[end])
+guess_offset = 0
+guess_bg_slope = 1
+fitfunc.initial_parameters = Float64[ guess_offset, guess_bg_slope ]
+fitfunc
+```
+
+* Performe the fit and plot the result
+```@example fitting_1D_data
+RadiationSpectra.lsqfit!(fitfunc, xdata, ydata, xdata_err, ydata_err) # xdata_err and ydata_err are optional
+
+plot(xdata, ydata, xerr=xdata_err, yerr=ydata_err, st=:scatter, size=(800,400), label="Data")
+plot!(fitfunc, use_initial_parameters=true, lc=:green, label="Guess")
+plot!(fitfunc, lc=:red, label="LSQ Fit")
+```
+
+```@example fitting_1D_data
+fitfunc # hide
+```
diff --git a/docs/src/man/spectrum_calibration.md b/docs/src/man/spectrum_calibration.md
@@ -0,0 +1,34 @@
+# Spectrum Calibration
+
+The function [`RadiationSpectra.calibrate_spectrum`](@ref) return a calibrated spectrum of an uncalibrated one. 
+As input it needs the uncalibrated spectrum (`StatsBase::Histogram`) and a `Vector` of known peak positions like photon lines.
+
+The calibration is based on the assumption that the calibration function is just ``f(x) = c\cdot x``.
+
+```@example spectrum_calibration
+using Plots, RadiationSpectra 
+myfont = Plots.font(12) # hide
+pyplot(guidefont=myfont, xtickfont=myfont, ytickfont=myfont, legendfont=myfont, titlefont=myfont) # hide
+
+h_uncal = RadiationSpectra.get_example_spectrum()
+photon_lines = [609.312, 911.204, 1120.287, 1460.830, 1764.494] # keV
+h_cal = RadiationSpectra.calibrate_spectrum(h_uncal, photon_lines)
+
+p_uncal = plot(h_uncal, st=:step, label="Uncalibrated spectrum"); 
+p_cal = plot(h_cal, st=:step, label="Calibrated spectrum", xlabel="E / keV"); 
+vline!(photon_lines, lw=0.5, color=:red, label="Photon lines")
+plot(p_uncal, p_cal, size=(800,500), layout=(2, 1)) 
+```
+
+Zoom into one of the peaks:
+```@example spectrum_calibration
+p_cal = plot(h_cal, st=:step, label="Calibrated spectrum", xlabel="E / keV", xlims=[1440, 1480], size=[800, 400]); # hide
+vline!([1460.830], label="Photon line") # hide
+```
+
+## Algorithm 
+
+1. Deconvolution -> Peak finding 
+2. Peak Identification - Which peak corresponds to which photon line? 
+3. Fit identified peaks 
+4. Fit determined peak positions vs 'true' positions (photon lines) to get the calibration constant ``c``.
diff --git a/docs/src/man/spectrum_deconvolution.md b/docs/src/man/spectrum_deconvolution.md
@@ -0,0 +1,17 @@
+# Peak Finding (Spectrum Deconvolution)
+
+See [`RadiationSpectra.peakfinder`](@ref)
+
+```@example
+using RadiationSpectra 
+h_uncal = RadiationSpectra.get_example_spectrum()
+h_decon, peakpos = RadiationSpectra.peakfinder(h_uncal)
+
+using Plots 
+myfont = Plots.font(12) # hide
+pyplot(guidefont=myfont, xtickfont=myfont, ytickfont=myfont, legendfont=myfont, titlefont=myfont) # hide
+p_uncal = plot(h_uncal, st=:step, label="Uncalibrated spectrum", c=1, lw=1.2); 
+p_decon = plot(peakpos, st=:vline, c=:red, label="Peaks", lw=0.3);
+plot!(h_decon, st=:step, label="Deconvoluted spectrum", c=1, lw=1.2); 
+plot(p_uncal, p_decon, size=(800,500), layout=(2, 1)) 
+```