diff --git a/src/LumenManufaktur.jl b/src/LumenManufaktur.jl
index 0b5cb89c6ad843e7aeccb5cd6bab51d1b99df64f..c01106e10538337a69756dc74f0f0f1d172188f3 100644
--- a/src/LumenManufaktur.jl
+++ b/src/LumenManufaktur.jl
@@ -12,7 +12,9 @@ include("scattering.jl")
include("absorption.jl")
include("pmt.jl")
include("parameters.jl")
+include("utils.jl")
include("muons.jl")
+include("deltarays.jl")
include("exports.jl")
end
diff --git a/src/deltarays.jl b/src/deltarays.jl
new file mode 100644
index 0000000000000000000000000000000000000000..312a7663495ed56e9809a546995395368adb58c4
--- /dev/null
+++ b/src/deltarays.jl
@@ -0,0 +1,268 @@
+"""
+ deltayprobability(x)
+
+Emission profile of photons from delta-rays.
+
+Profile is taken from reference ANTARES-SOFT-2002-015, J. Brunner (fig. 3).
+
+Returns the probability.
+
+# Arguments
+- `x`: cosine emission angle
+"""
+function deltarayprobability(x)
+ 0.188 * exp(-1.25 * abs(x - 0.90)^1.30)
+end
+
+
+"""
+ directlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
+
+Probability density function for direct light from delta-rays with a
+closest distance of `R` [m] to the PMT, the angles `θ` \\[rad\\] (zenith) and `ϕ`
+\\[rad\\] (azimuth) with respect to the PMT axis and a time difference `Δt` [ns]
+relative to direct Cherenkov light. Returns d^2P/dt/dE [npe/ns ⋅ m/GeV].
+
+# Arguments
+
+- `params`: parameters of the setup
+- `pmt`: the PMT model
+- `R`: (closest) distance [m] between muon and PMT
+- `ϕ`: zenith angle \\[rad\\] which is 0 when the PMT points away from the muon
+ track (in x-direction) and rotates counter clockwise to the y-axis when
+ viewed from above, while the z-axis points upwards.
+- `θ`: azimuth angle \\[rad\\] which is 0 when the PMT points upwards (along the
+ z-axis) and π/2 when pointing downwards. Between 0 and π/2, the PMT points to
+ the z-axis
+- `Δt`: time difference [ns] relative to the Cherenkov light
+"""
+function directlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
+ value = 0.0
+
+ R = max(R, params.minimum_distance)
+ A = pmt.photocathode_area
+ t = R * tanthetac(params.n) / C + Δt # time [ns]
+ D = 2.0 * sqrt(A / π)
+
+ px = sin(θ) * cos(ϕ)
+ pz = cos(θ)
+
+ n0 = refractionindexgroup(params.dispersion_model, params.lambda_max)
+ n1 = refractionindexgroup(params.dispersion_model, params.lambda_min)
+ ni = sqrt(R * R + C * t * C * t) / R # maximal index of refraction
+
+ n0 >= ni && return value
+
+ nj = min(n1, ni)
+
+ w = params.lambda_max
+
+ for (m_x, m_y) in zip(params.legendre_coefficients...)
+
+ ng = 0.5 * (nj + n0) + m_x * 0.5 * (nj - n0)
+ dn = m_y * 0.5 * (nj - n0)
+
+ w = wavelength(params.dispersion_model, ng, w, 1.0e-5)
+
+ dw = dn / abs(dispersiongroup(params.dispersion_model, w))
+
+ n = refractionindexphase(params.dispersion_model, w)
+
+ l_abs = absorptionlength(params.absorption_model, w)
+ ls = scatteringlength(params.scattering_model, w)
+
+ npe = cherenkov(w, n) * dw * pmt.quantum_efficiency(w)
+
+ root = Root(R, ng, t) # square root
+
+ !root.isvalid && continue
+
+ for z in (root.most_upstream, root.most_downstream)
+
+ C * t <= z && continue
+
+ d = sqrt(z * z + R * R) # distance traveled by photon
+
+ ct0 = -z / d
+ st0 = R / d
+
+ dz = D / st0 # average track length
+ ct = st0 * px + ct0 * pz # cosine angle of incidence on PMT
+
+ U = pmt.angular_acceptance(ct) # PMT angular acceptance
+ V = exp(-d / l_abs - d / ls) # absorption & scattering
+ W = A / (d * d) # solid angle
+
+ Ja = deltarayprobability(ct0) # d^2N/dcos/dϕ
+ Jd = (1.0 - ng * ct0) / C # d t/ dz
+ Je = ng * st0 * st0 * st0 / (R * C) # d^2t/(dz)^2
+
+ value += npe * geanc() * U * V * W * Ja / (abs(Jd) + 0.5 * Je * dz)
+ end
+
+ end
+ return value
+end
+
+"""
+ scatteredlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
+
+Probability density function for scattered light from delta-rays with a
+closest distance of `R` [m] to the PMT, the angles `θ` \\[rad\\] (zenith) and `ϕ`
+\\[rad\\] (azimuth) with respect to the PMT axis and a time difference `Δt` [ns]
+relative to direct Cherenkov light. Returns d^2P/dt/dE [npe/ns ⋅ m/GeV].
+
+# Arguments
+
+- `params`: parameters of the setup
+- `pmt`: the PMT model
+- `R`: (closest) distance [m] between muon and PMT
+- `ϕ`: zenith angle \\[rad\\] which is 0 when the PMT points away from the muon
+ track (in x-direction) and rotates counter clockwise to the y-axis when
+ viewed from above, while the z-axis points upwards.
+- `θ`: azimuth angle \\[rad\\] which is 0 when the PMT points upwards (along the
+ z-axis) and π/2 when pointing downwards. Between 0 and π/2, the PMT points to
+ the z-axis
+- `Δt`: time difference [ns] relative to the Cherenkov light
+"""
+function scatteredlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
+ value = 0.0
+
+ R = max(R, params.minimum_distance)
+ A = pmt.photocathode_area
+ t = R * tanthetac(params.n) / C + Δt # time [ns]
+
+ px = sin(θ) * cos(ϕ)
+ py = sin(θ) * sin(ϕ)
+ pz = cos(θ)
+
+ n0 = refractionindexgroup(params.dispersion_model, params.lambda_max)
+ n1 = refractionindexgroup(params.dispersion_model, params.lambda_min)
+ ni = sqrt(R^2 + C^2 * t^2) / R # maximal index of refraction
+
+ n0 >= ni && return value
+
+ nj = min(ni, n1)
+
+ w = params.lambda_max
+
+
+ for (m_x, m_y) in zip(params.legendre_coefficients...)
+
+ ng = 0.5 * (nj + n0) + m_x * 0.5 * (nj - n0)
+ dn = m_y * 0.5 * (nj - n0)
+
+ w = wavelength(params.dispersion_model, ng, w, 1.0e-5)
+
+ dw = dn / abs(dispersiongroup(params.dispersion_model, w))
+
+ n = refractionindexphase(params.dispersion_model, w)
+
+ l_abs = absorptionlength(params.absorption_model, w)
+ ls = scatteringlength(params.scattering_model, w)
+
+
+ npe = cherenkov(w, n) * dw * pmt.quantum_efficiency(w)
+
+ npe <= 0 && continue
+
+ Jc = 1.0 / ls # dN/dx
+
+ root = Root(R, ng, t) # square root
+
+ !root.isvalid && continue
+
+ zmin = root.most_upstream
+ zmax = root.most_downstream
+
+ zap = 1.0 / l_abs
+
+ xmin = exp(zap * zmax)
+ xmax = exp(zap * zmin)
+
+ n_coefficients = length(params.legendre_coefficients[1])
+
+ for (p_x, p_y) in zip(params.legendre_coefficients...)
+
+ x = 0.5 * (xmax + xmin) + p_x * 0.5 * (xmax - xmin)
+ dx = p_y * 0.5 * (xmax - xmin)
+
+ z = log(x) / zap
+ dz = -dx / (zap * x)
+
+ D = sqrt(z^2 + R^2)
+ cd = -z / D
+ sd = R / D
+
+ qx = cd * px + 0 - sd * pz
+ qy = 1 * py
+ qz = sd * px + 0 + cd * pz
+
+ d = (C * t - z) / ng # photon path
+
+ ds = 2.0 / (n_coefficients + 1)
+
+ sb = 0.5ds
+ while sb < 1.0 - 0.25ds
+
+ for k in (-1, 1)
+
+ cb = k * sqrt((1.0 + sb) * (1.0 - sb))
+ dcb = k * ds * sb / cb
+
+ v = 0.5 * (d + D) * (d - D) / (d - D * cb)
+ u = d - v
+
+ u <= 0 && continue
+ v <= 0 && continue
+
+ cts = (D * cb - v) / u # cosine scattering angle
+
+ V = exp(
+ -d * inverseattenuationlength(
+ params.scattering_probability_model,
+ l_abs,
+ ls,
+ cts,
+ ),
+ )
+
+ cts < 0.0 &&
+ v * sqrt((1.0 + cts) * (1.0 - cts)) < params.module_radius &&
+ continue
+
+ W = min(A / (v * v), 2.0 * π) # solid angle
+ Ja = scatteringprobability(params.scattering_probability_model, cts) # d^2P/dcos/dϕ
+ Jd = ng * (1.0 - cts) / C # dt/du
+
+ ca = (D - v * cb) / u
+ sa = v * sb / u
+
+ dp = π / length(params.integration_points)
+ dom = dcb * dp * v^2 / (u^2)
+
+ for (cp, sp) in params.integration_points.xy
+ ct0 = cd * ca - sd * sa * cp
+
+ vx = sb * cp * qx
+ vy = sb * sp * qy
+ vz = cb * qz
+
+ U =
+ pmt.angular_acceptance(vx + vy + vz) +
+ pmt.angular_acceptance(vx - vy + vz) # PMT angular acceptance
+
+ Jb = deltarayprobability(ct0) # d^2N/dcos/dϕ
+
+ value +=
+ dom * npe * geanc() * dz * U * V * W * Ja * Jb * Jc / abs(Jd)
+ end
+ end
+ sb += ds
+
+ end
+ end
+ end
+
+ value
+end
diff --git a/src/muons.jl b/src/muons.jl
index 88a5a4f47ab9036f9900e2f62c0ede5ac844eb23..b45bc4c241cf3ecefb8bcfef31610e9f5cee838c 100644
--- a/src/muons.jl
+++ b/src/muons.jl
@@ -15,24 +15,6 @@ function cherenkov(λ, n)
end
-"""
- deltayprobability(x)
-
-Emission profile of photons from delta-rays.
-
-Profile is taken from reference ANTARES-SOFT-2002-015, J. Brunner (fig. 3).
-
-Returns the probability.
-
-# Arguments
-- `x`: cosine emission angle
-"""
-function deltarayprobability(x)
- 0.188 * exp(-1.25 * abs(x - 0.90)^1.30)
-end
-
-
-
"""
geanc()
@@ -315,307 +297,3 @@ function scatteredlightfrommuon(params::LMParameters, pmt::PMTModel, D, cd, θ,
return value
end
-
-
-
-"""
- directlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
-
-Probability density function for direct light from delta-rays with a
-closest distance of `R` [m] to the PMT, the angles `θ` \\[rad\\] (zenith) and `ϕ`
-\\[rad\\] (azimuth) with respect to the PMT axis and a time difference `Δt` [ns]
-relative to direct Cherenkov light. Returns d^2P/dt/dE [npe/ns ⋅ m/GeV].
-
-# Arguments
-
-- `params`: parameters of the setup
-- `pmt`: the PMT model
-- `R`: (closest) distance [m] between muon and PMT
-- `ϕ`: zenith angle \\[rad\\] which is 0 when the PMT points away from the muon
- track (in x-direction) and rotates counter clockwise to the y-axis when
- viewed from above, while the z-axis points upwards.
-- `θ`: azimuth angle \\[rad\\] which is 0 when the PMT points upwards (along the
- z-axis) and π/2 when pointing downwards. Between 0 and π/2, the PMT points to
- the z-axis
-- `Δt`: time difference [ns] relative to the Cherenkov light
-"""
-function directlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
- value = 0.0
-
- R = max(R, params.minimum_distance)
- A = pmt.photocathode_area
- t = R * tanthetac(params.n) / C + Δt # time [ns]
- D = 2.0 * sqrt(A / π)
-
- px = sin(θ) * cos(ϕ)
- pz = cos(θ)
-
- n0 = refractionindexgroup(params.dispersion_model, params.lambda_max)
- n1 = refractionindexgroup(params.dispersion_model, params.lambda_min)
- ni = sqrt(R * R + C * t * C * t) / R # maximal index of refraction
-
- n0 >= ni && return value
-
- nj = min(n1, ni)
-
- w = params.lambda_max
-
- for (m_x, m_y) in zip(params.legendre_coefficients...)
-
- ng = 0.5 * (nj + n0) + m_x * 0.5 * (nj - n0)
- dn = m_y * 0.5 * (nj - n0)
-
- w = wavelength(params.dispersion_model, ng, w, 1.0e-5)
-
- dw = dn / abs(dispersiongroup(params.dispersion_model, w))
-
- n = refractionindexphase(params.dispersion_model, w)
-
- l_abs = absorptionlength(params.absorption_model, w)
- ls = scatteringlength(params.scattering_model, w)
-
- npe = cherenkov(w, n) * dw * pmt.quantum_efficiency(w)
-
- root = Root(R, ng, t) # square root
-
- !root.isvalid && continue
-
- for z in (root.most_upstream, root.most_downstream)
-
- C * t <= z && continue
-
- d = sqrt(z * z + R * R) # distance traveled by photon
-
- ct0 = -z / d
- st0 = R / d
-
- dz = D / st0 # average track length
- ct = st0 * px + ct0 * pz # cosine angle of incidence on PMT
-
- U = pmt.angular_acceptance(ct) # PMT angular acceptance
- V = exp(-d / l_abs - d / ls) # absorption & scattering
- W = A / (d * d) # solid angle
-
- Ja = deltarayprobability(ct0) # d^2N/dcos/dϕ
- Jd = (1.0 - ng * ct0) / C # d t/ dz
- Je = ng * st0 * st0 * st0 / (R * C) # d^2t/(dz)^2
-
- value += npe * geanc() * U * V * W * Ja / (abs(Jd) + 0.5 * Je * dz)
- end
-
- end
- return value
-end
-
-"""
- scatteredlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
-
-Probability density function for scattered light from delta-rays with a
-closest distance of `R` [m] to the PMT, the angles `θ` \\[rad\\] (zenith) and `ϕ`
-\\[rad\\] (azimuth) with respect to the PMT axis and a time difference `Δt` [ns]
-relative to direct Cherenkov light. Returns d^2P/dt/dE [npe/ns ⋅ m/GeV].
-
-# Arguments
-
-- `params`: parameters of the setup
-- `pmt`: the PMT model
-- `R`: (closest) distance [m] between muon and PMT
-- `ϕ`: zenith angle \\[rad\\] which is 0 when the PMT points away from the muon
- track (in x-direction) and rotates counter clockwise to the y-axis when
- viewed from above, while the z-axis points upwards.
-- `θ`: azimuth angle \\[rad\\] which is 0 when the PMT points upwards (along the
- z-axis) and π/2 when pointing downwards. Between 0 and π/2, the PMT points to
- the z-axis
-- `Δt`: time difference [ns] relative to the Cherenkov light
-"""
-function scatteredlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
- value = 0.0
-
- R = max(R, params.minimum_distance)
- A = pmt.photocathode_area
- t = R * tanthetac(params.n) / C + Δt # time [ns]
-
- px = sin(θ) * cos(ϕ)
- py = sin(θ) * sin(ϕ)
- pz = cos(θ)
-
- n0 = refractionindexgroup(params.dispersion_model, params.lambda_max)
- n1 = refractionindexgroup(params.dispersion_model, params.lambda_min)
- ni = sqrt(R^2 + C^2 * t^2) / R # maximal index of refraction
-
- n0 >= ni && return value
-
- nj = min(ni, n1)
-
- w = params.lambda_max
-
-
- for (m_x, m_y) in zip(params.legendre_coefficients...)
-
- ng = 0.5 * (nj + n0) + m_x * 0.5 * (nj - n0)
- dn = m_y * 0.5 * (nj - n0)
-
- w = wavelength(params.dispersion_model, ng, w, 1.0e-5)
-
- dw = dn / abs(dispersiongroup(params.dispersion_model, w))
-
- n = refractionindexphase(params.dispersion_model, w)
-
- l_abs = absorptionlength(params.absorption_model, w)
- ls = scatteringlength(params.scattering_model, w)
-
-
- npe = cherenkov(w, n) * dw * pmt.quantum_efficiency(w)
-
- npe <= 0 && continue
-
- Jc = 1.0 / ls # dN/dx
-
- root = Root(R, ng, t) # square root
-
- !root.isvalid && continue
-
- zmin = root.most_upstream
- zmax = root.most_downstream
-
- zap = 1.0 / l_abs
-
- xmin = exp(zap * zmax)
- xmax = exp(zap * zmin)
-
- n_coefficients = length(params.legendre_coefficients[1])
-
- for (p_x, p_y) in zip(params.legendre_coefficients...)
-
- x = 0.5 * (xmax + xmin) + p_x * 0.5 * (xmax - xmin)
- dx = p_y * 0.5 * (xmax - xmin)
-
- z = log(x) / zap
- dz = -dx / (zap * x)
-
- D = sqrt(z^2 + R^2)
- cd = -z / D
- sd = R / D
-
- qx = cd * px + 0 - sd * pz
- qy = 1 * py
- qz = sd * px + 0 + cd * pz
-
- d = (C * t - z) / ng # photon path
-
- ds = 2.0 / (n_coefficients + 1)
-
- sb = 0.5ds
- while sb < 1.0 - 0.25ds
-
- for k in (-1, 1)
-
- cb = k * sqrt((1.0 + sb) * (1.0 - sb))
- dcb = k * ds * sb / cb
-
- v = 0.5 * (d + D) * (d - D) / (d - D * cb)
- u = d - v
-
- u <= 0 && continue
- v <= 0 && continue
-
- cts = (D * cb - v) / u # cosine scattering angle
-
- V = exp(
- -d * inverseattenuationlength(
- params.scattering_probability_model,
- l_abs,
- ls,
- cts,
- ),
- )
-
- cts < 0.0 &&
- v * sqrt((1.0 + cts) * (1.0 - cts)) < params.module_radius &&
- continue
-
- W = min(A / (v * v), 2.0 * π) # solid angle
- Ja = scatteringprobability(params.scattering_probability_model, cts) # d^2P/dcos/dϕ
- Jd = ng * (1.0 - cts) / C # dt/du
-
- ca = (D - v * cb) / u
- sa = v * sb / u
-
- dp = π / length(params.integration_points)
- dom = dcb * dp * v^2 / (u^2)
-
- for (cp, sp) in params.integration_points.xy
- ct0 = cd * ca - sd * sa * cp
-
- vx = sb * cp * qx
- vy = sb * sp * qy
- vz = cb * qz
-
- U =
- pmt.angular_acceptance(vx + vy + vz) +
- pmt.angular_acceptance(vx - vy + vz) # PMT angular acceptance
-
- Jb = deltarayprobability(ct0) # d^2N/dcos/dϕ
-
- value +=
- dom * npe * geanc() * dz * U * V * W * Ja * Jb * Jc / abs(Jd)
- end
- end
- sb += ds
-
- end
- end
- end
-
- value
-end
-
-"""
- Root(R, n, t)
-
-Helper struct to find solution(s) to ``z`` of the square root expression:
-
-```math
-\begin{aligned}
-ct(z=0) & = & z + n \\sqrt{z^2 + R^2}
-\end{aligned}
-```
-
-where ``n = 1/\\cos(\\theta_{c})`` is the index of refraction.
-
-# Arguments
-
-- `R`: minimal distance of approach [m]
-- `n`: index of refraction
-- `t`: time at z = 0 [ns]
-
-# Fields
-
-- `most_upstream`: most upstream solution
-- `most_downstream`: most downstream solution
-- `isvalid`: validity of the solution
-```
-"""
-struct Root
- most_upstream::Float64
- most_downstream::Float64
- isvalid::Bool
-
- function Root(R, n, t)
- a = n^2 - 1.0
- b = 2 * C * t
- c = R^2 * n^2 - C^2 * t^2
-
- q = b * b - 4 * a * c
-
- isvalid = false
- if q >= 0.0
-
- first = (-b - sqrt(q)) / (2a)
- second = (-b + sqrt(q)) / (2a)
-
- isvalid = C * t > second
- end
- new(first, second, isvalid)
- end
-end
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 0000000000000000000000000000000000000000..67b4b243028a5578895a6133ca761987ca9acaa2
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,49 @@
+"""
+ Root(R, n, t)
+
+Helper struct to find solution(s) to ``z`` of the square root expression:
+
+```math
+\begin{aligned}
+ct(z=0) & = & z + n \\sqrt{z^2 + R^2}
+\end{aligned}
+```
+
+where ``n = 1/\\cos(\\theta_{c})`` is the index of refraction.
+
+# Arguments
+
+- `R`: minimal distance of approach [m]
+- `n`: index of refraction
+- `t`: time at z = 0 [ns]
+
+# Fields
+
+- `most_upstream`: most upstream solution
+- `most_downstream`: most downstream solution
+- `isvalid`: validity of the solution
+```
+"""
+struct Root
+ most_upstream::Float64
+ most_downstream::Float64
+ isvalid::Bool
+
+ function Root(R, n, t)
+ a = n^2 - 1.0
+ b = 2 * C * t
+ c = R^2 * n^2 - C^2 * t^2
+
+ q = b * b - 4 * a * c
+
+ isvalid = false
+ if q >= 0.0
+
+ first = (-b - sqrt(q)) / (2a)
+ second = (-b + sqrt(q)) / (2a)
+
+ isvalid = C * t > second
+ end
+ new(first, second, isvalid)
+ end
+end