diff --git a/notebooks/Direct_Light_from_Muon.jl b/notebooks/Direct_Light_from_Muon.jl index 0aa887542f2d86c7c34405b2a786c4c42621acda..36dcb2cdd66da03adab71413c1466eb739b28e7e 100644 --- a/notebooks/Direct_Light_from_Muon.jl +++ b/notebooks/Direct_Light_from_Muon.jl @@ -6,56 +6,71 @@ using InteractiveUtils # ╔═╡ 6f448138-0b23-11ef-1b83-859c53214856 begin - using CairoMakie - using LumenManufaktur + using CairoMakie + using LumenManufaktur end # ╔═╡ 7292e824-7de9-4d1f-ab2d-895f6fadfd88 pkgversion(LumenManufaktur) # ╔═╡ 2a28e349-c9a1-4fb1-a06b-b1b616c20597 -directlightfrommuon(LMParameters(), LumenManufaktur.KM3NeTPMT, 5.23, π, π/2) +directlightfrommuon(LMParameters(), LumenManufaktur.KM3NeTPMT, 5.23, π, π / 2) # ╔═╡ aa28dd35-653e-4433-b0af-4d6050696bc2 -directlightfrommuon(LMParameters(), LumenManufaktur.KM3NeTPMT, 1.15, π/1.5, π/2, 0.01) +directlightfrommuon(LMParameters(), LumenManufaktur.KM3NeTPMT, 1.15, π / 1.5, π / 2, 0.01) # ╔═╡ c89f6755-8dfe-4f8e-a640-d5f52b21c9b6 let - fig = Figure(size = (600, 400)); - ax = Axis( - fig[1, 1], - xlabel = "distance between muon and PMT / m", - ylabel = "npe", - xgridstyle = :dash, - ygridstyle = :dash - ) - Rs = range(1, 10, 1000) - params = LMParameters(dispersion_model=BaileyDispersion(240)) - pmt = LumenManufaktur.KM3NeTPMT - lines!(ax, Rs, [directlightfrommuon(params, pmt, R, π, π/2) for R in Rs], label="θ = π") - lines!(ax, Rs, [directlightfrommuon(params, pmt, R, π/2, π/2) for R in Rs], label="θ = π/2") - axislegend(; position = :rt) - fig + fig = Figure(size = (600, 400)) + ax = Axis( + fig[1, 1], + xlabel = "distance between muon and PMT / m", + ylabel = "npe", + xgridstyle = :dash, + ygridstyle = :dash, + ) + Rs = range(1, 10, 1000) + params = LMParameters(dispersion_model = BaileyDispersion(240)) + pmt = LumenManufaktur.KM3NeTPMT + lines!( + ax, + Rs, + [directlightfrommuon(params, pmt, R, π, π / 2) for R in Rs], + label = "θ = π", + ) + lines!( + ax, + Rs, + [directlightfrommuon(params, pmt, R, π / 2, π / 2) for R in Rs], + label = "θ = π/2", + ) + axislegend(; position = :rt) + fig end # ╔═╡ b6bd1af8-f5af-4536-b903-f894f4bee50f let - fig = Figure(size = (600, 400)); - ax = Axis( - fig[1, 1], - xlabel = "PMT zenith angle [rad]", - ylabel = "npe", - xgridstyle = :dash, - ygridstyle = :dash - ) - θs = range(-π, π, 1000) - params = LMParameters(dispersion_model=BaileyDispersion(240)) - pmt = LumenManufaktur.KM3NeTPMT - for R in 2:10 - lines!(ax, θs, [directlightfrommuon(params, pmt, R, θ, π/2) for θ in θs], label="R = $R m") - end - axislegend(; position = :ct) - fig + fig = Figure(size = (600, 400)) + ax = Axis( + fig[1, 1], + xlabel = "PMT zenith angle [rad]", + ylabel = "npe", + xgridstyle = :dash, + ygridstyle = :dash, + ) + θs = range(-π, π, 1000) + params = LMParameters(dispersion_model = BaileyDispersion(240)) + pmt = LumenManufaktur.KM3NeTPMT + for R = 2:10 + lines!( + ax, + θs, + [directlightfrommuon(params, pmt, R, θ, π / 2) for θ in θs], + label = "R = $R m", + ) + end + axislegend(; position = :ct) + fig end # ╔═╡ a3027c77-7f9f-4309-bca2-7a6288eba83a @@ -68,27 +83,23 @@ LumenManufaktur.Kopelevich() LumenManufaktur.KM3NeTPMT # ╔═╡ a196cb83-6423-4360-b335-0c2fa55183b5 -ANTARESPMT = LumenManufaktur.PMTModel( - 440e-4, - λ -> rand()*20, - λ -> rand(), -) +ANTARESPMT = LumenManufaktur.PMTModel(440e-4, λ -> rand() * 20, λ -> rand()) # ╔═╡ 4d6f0b13-98ea-407a-b23d-23ac47df4a3e let - fig = Figure(size = (600, 400)); - ax = Axis( - fig[1, 1], - xlabel = "distance between muon and PMT / m", - ylabel = "npe", - xgridstyle = :dash, - ygridstyle = :dash - ) - Rs = range(1, 10, 1000) - params = LMParameters(dispersion_model=BaileyDispersion(240)) - pmt = LumenManufaktur.KM3NeTPMT - lines!(ax, Rs, [directlightfrommuon(params, ANTARESPMT, R, π, π/2) for R in Rs]) - fig + fig = Figure(size = (600, 400)) + ax = Axis( + fig[1, 1], + xlabel = "distance between muon and PMT / m", + ylabel = "npe", + xgridstyle = :dash, + ygridstyle = :dash, + ) + Rs = range(1, 10, 1000) + params = LMParameters(dispersion_model = BaileyDispersion(240)) + pmt = LumenManufaktur.KM3NeTPMT + lines!(ax, Rs, [directlightfrommuon(params, ANTARESPMT, R, π, π / 2) for R in Rs]) + fig end # ╔═╡ 00000000-0000-0000-0000-000000000001 diff --git a/src/dispersion.jl b/src/dispersion.jl index df420babc44ac8f43db0eb811dab5463228b98ea..9aed02b498348627b1ccdd66b02f785433f32857 100644 --- a/src/dispersion.jl +++ b/src/dispersion.jl @@ -34,23 +34,23 @@ end @inline function refractionindexgroup(dp::BaileyDispersion, λ) n = refractionindexphase(dp, λ) y = dispersionphase(dp, λ) - n / (1.0 + y*λ/n) + n / (1.0 + y * λ / n) end @inline function dispersionphase(dp::BaileyDispersion, λ) x = 1.0 / λ - -(x^2)*(dp.a2 + x*(2.0*dp.a3 + x*3.0*dp.a4)) + -(x^2) * (dp.a2 + x * (2.0 * dp.a3 + x * 3.0 * dp.a4)) end @inline function dispersiongroup(dp::BaileyDispersion, λ) x = 1.0 / λ - n = refractionindexphase(dp, λ) - np = dispersionphase(dp, λ) - npp = x^3*(2.0*dp.a2 + x*(6.0*dp.a3 + x*12.0*dp.a4)) - ng = n / (1.0 + np*λ/n) + n = refractionindexphase(dp, λ) + np = dispersionphase(dp, λ) + npp = x^3 * (2.0 * dp.a2 + x * (6.0 * dp.a3 + x * 12.0 * dp.a4)) + ng = n / (1.0 + np * λ / n) - ng^2 * (2*np^2 - n*npp) * λ / (n^3); + ng^2 * (2 * np^2 - n * npp) * λ / (n^3) end """ @@ -71,9 +71,9 @@ function wavelength(dp::BaileyDispersion, n, w, eps) v = w while true - y = refractionindexgroup(dp, v) - abs(y - n) < eps && break - v += (n - y) / dispersiongroup(dp, v) + y = refractionindexgroup(dp, v) + abs(y - n) < eps && break + v += (n - y) / dispersiongroup(dp, v) end v diff --git a/src/parameters.jl b/src/parameters.jl index 5ff9073ec5b475c05fb65c958ead1890a54f5061..2c6174916010dd31b9f564995577c7095fb0266d 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -13,7 +13,12 @@ The parameter set for light detection. - `scattering_probability_model`: the scattering probability model (default: p00075) - `absorption_model`: the absorption model """ -Base.@kwdef struct LMParameters{D<:DispersionModel,S<:ScatteringModel,SP<:ScatteringProbabilityModel,A<:AbsorptionModel} +Base.@kwdef struct LMParameters{ + D<:DispersionModel, + S<:ScatteringModel, + SP<:ScatteringProbabilityModel, + A<:AbsorptionModel, +} minimum_distance::Float64 = 1.0e-1 module_radius::Float64 = 0.25 lambda_min::Float64 = 300.0 diff --git a/src/scattering.jl b/src/scattering.jl index 58413f0b798eaabe91304d8fcc26d5e767f0e7cc..7bfe05de17f7a88f79d53902ccdd616a5c00cd80 100644 --- a/src/scattering.jl +++ b/src/scattering.jl @@ -49,7 +49,7 @@ Model specific function to describe light scattering probability in water (p0007 g = 0.924 f = 0.17 - f * rayleigh(x) + (1.0 - f) * henyey_greenstein(g, x) + f * rayleigh(x) + (1.0 - f) * henyey_greenstein(g, x) end """ @@ -61,7 +61,7 @@ Light scattering probability in water (Heneyey-Greenstein). - `x`: cosine scattering angle """ @inline function henyey_greenstein(x) - g = 0.924; + g = 0.924 return henyey_greenstein(g, x) end @@ -76,9 +76,9 @@ Light scattering probability in water (Heneyey-Greenstein). """ @inline function henyey_greenstein(g, x) a0 = (1.0 - g^2) / (4π) - y = 1.0 + g^2 - 2.0*g*x + y = 1.0 + g^2 - 2.0 * g * x - a0 / (y*sqrt(y)) + a0 / (y * sqrt(y)) end """ @@ -101,8 +101,8 @@ Light scattering probability in water (Rayleigh). - `x`: cosine scattering angle """ @inline function rayleigh(g, x) - a0 = 1.0 / (1.0 + g/3.0) / (4π) - a0 * (1.0 + g*x^2) + a0 = 1.0 / (1.0 + g / 3.0) / (4π) + a0 * (1.0 + g * x^2) end @@ -119,21 +119,21 @@ Get the inverse of the attenuation length [m^-1]. """ function inverseattenuationlength(::Scatteringp00075, l_abs, ls, cts) - 1.0 / l_abs + inverseattenuationlengthinterpolator(cts) / ls; + 1.0 / l_abs + inverseattenuationlengthinterpolator(cts) / ls end """ Interpolator for the p00075 model based inverse attenutation calculation. """ const inverseattenuationlengthinterpolator = let - xs = range(-1, 1; length=100000) + xs = range(-1, 1; length = 100000) dx = xs.step.hi xs = collect(xs) ys = Float64[] W = 0.0 for x in xs push!(ys, W) - W += 2π * dx * scatteringprobability(Scatteringp00075(), x+0.5*dx) + W += 2π * dx * scatteringprobability(Scatteringp00075(), x + 0.5 * dx) end # xs[1] = 0.0 # xs[end] = 1.0