Skip to content
Snippets Groups Projects
Verified Commit ab75cde8 authored by Tamas Gal's avatar Tamas Gal :speech_balloon:
Browse files

Implement scattered light from deltarays

parent 02232575
No related branches found
No related tags found
No related merge requests found
Pipeline #51519 passed
Stage: test
Stage: docs
Hide whitespace changes
Inline Side-by-side
src/exports.jl
+ 1
− 0
View file @ ab75cde8
......@@ -4,6 +4,7 @@ export
directlightfromdeltarays,
scatteredlightfromdeltarays,
LMParameters,
PDFIntegrationPoints,
PMTModel,
absorptionlength,
scatteringlength,
......
src/muons.jl
+ 120
− 114
View file @ ab75cde8
......@@ -386,8 +386,8 @@ function directlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ
d = sqrt(z * z + R * R) # distance traveled by photon
ct0 = -z / d;
st0 = R / d;
ct0 = -z / d
st0 = R / d
dz = D / st0 # average track length
ct = st0 * px + ct0 * pz # cosine angle of incidence on PMT
......@@ -400,171 +400,177 @@ function directlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ
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);
value += npe * geanc() * U * V * W * Ja / (abs(Jd) + 0.5 * Je * dz)
end
end
return value
end
# /**
# * Probability density function for scattered light from delta-rays.
# *
# * \param R_m distance between muon and PMT [m]
# * \param θ zenith angle orientation PMT [rad]
# * \param ϕ azimuth angle orientation PMT [rad]
# * \param Δt time difference relative to direct Cherenkov light [ns]
# * \return d^2P/dt/dx [npe/ns x m/GeV]
# */
# double getScatteredLightFromDeltaRays(const double R_m, const double θ,
# const double ϕ,
# const double Δt) const {
# double value = 0;
"""
directlightfromdeltarays(params::LMParameters, pmt::PMTModel, R, θ, ϕ, Δt)
# const double R = std::max(R_m, getRmin());
# const double t = R * getTanΘC() / C + Δt; // time [ns]
# const double A = pmt.photocathode_area;
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].
# const double px = sin(θ) * cos(ϕ);
# const double py = sin(θ) * sin(ϕ);
# const double pz = cos(θ);
# Arguments
# const double n0 = refractionindexgroup(params.dispersion_model, params.lambda_max);
# const double n1 = refractionindexgroup(params.dispersion_model, params.lambda_min);
# const double ni =
# sqrt(R * R + C * t * C * t) / R; // maximal index of refraction
- `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
# if (n0 >= ni) {
# return value;
# }
R = max(R, params.minimum_distance)
A = pmt.photocathode_area
t = R * tanthetac(params.n) / C + Δt # time [ns]
# const double nj = std::min(ni, n1);
px = sin(θ) * cos(ϕ)
py = sin(θ) * sin(ϕ)
pz = cos(θ)
# double w = params.lambda_max;
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
# for (const_iterator i = begin(); i != end(); ++i) {
n0 >= ni && return value
# const double ng = 0.5 * (nj + n0) + i->getX() * 0.5 * (nj - n0);
# const double dn = i->getY() * 0.5 * (nj - n0);
nj = min(ni, n1)
# w = getWavelength(ng, w, 1.0e-5);
w = params.lambda_max
# const double dw = dn / fabs(getDispersionGroup(w));
# const double n = getIndexOfRefractionPhase(w);
for (m_x, m_y) in zip(params.legendre_coefficients...)
# const double l_abs = getAbsorptionLength(w);
# const double ls = getScatteringLength(w);
ng = 0.5 * (nj + n0) + m_x * 0.5 * (nj - n0)
dn = m_y * 0.5 * (nj - n0)
# const double npe = cherenkov(w, n) * dw * getQE(w);
w = wavelength(params.dispersion_model, ng, w, 1.0e-5)
# if (npe <= 0) {
# continue;
# }
dw = dn / abs(dispersiongroup(params.dispersion_model, w))
# const double Jc = 1.0 / ls; // dN/dx
n = refractionindexphase(params.dispersion_model, w)
# JRoot rz(R, ng, t); // square root
l_abs = absorptionlength(params.absorption_model, w)
ls = scatteringlength(params.scattering_model, w)
# if (!rz.is_valid) {
# continue;
# }
# const double zmin = rz.first;
# const double zmax = rz.second;
npe = cherenkov(w, n) * dw * pmt.quantum_efficiency(w)
# const double zap = 1.0 / l_abs;
npe <= 0 && continue
# const double xmin = exp(zap * zmax);
# const double xmax = exp(zap * zmin);
Jc = 1.0 / ls # dN/dx
root = Root(R, ng, t) # square root
!root.isvalid && continue
zmin = root.most_upstream
zmax = root.most_downstream
# for (const_iterator j = begin(); j != end(); ++j) {
zap = 1.0 / l_abs
# const double x = 0.5 * (xmax + xmin) + j->getX() * 0.5 * (xmax - xmin);
# const double dx = j->getY() * 0.5 * (xmax - xmin);
xmin = exp(zap * zmax)
xmax = exp(zap * zmin)
# const double z = log(x) / zap;
# const double dz = -dx / (zap * x);
n_coefficients = length(params.legendre_coefficients[1])
# const double D = sqrt(z * z + R * R);
# const double cd = -z / D;
# const double sd = R / D;
for (p_x, p_y) in zip(params.legendre_coefficients...)
# const double qx = cd * px + 0 - sd * pz;
# const double qy = 1 * py;
# const double qz = sd * px + 0 + cd * pz;
x = 0.5 * (xmax + xmin) + p_x * 0.5 * (xmax - xmin)
dx = p_y * 0.5 * (xmax - xmin)
# const double d = (C * t - z) / ng; // photon path
z = log(x) / zap
dz = -dx / (zap * x)
# // const double V = exp(-d/l_abs); //
# // absorption
D = sqrt(z^2 + R^2)
cd = -z / D
sd = R / D
# const double ds = 2.0 / (size() + 1);
qx = cd * px + 0 - sd * pz
qy = 1 * py
qz = sd * px + 0 + cd * pz
# for (double sb = 0.5 * ds; sb < 1.0 - 0.25 * ds; sb += ds) {
d = (C * t - z) / ng # photon path
# for (int buffer[] = {-1, +1, 0}, *k = buffer; *k != 0; ++k) {
ds = 2.0 / (n_coefficients + 1)
# const double cb = (*k) * sqrt((1.0 + sb) * (1.0 - sb));
# const double dcb = (*k) * ds * sb / cb;
# const double v = 0.5 * (d + D) * (d - D) / (d - D * cb);
# const double u = d - v;
sb = 0.5ds
# for (double sb = 0.5 * ds; sb < 1.0 - 0.25 * ds; sb += ds) {
while sb < 1.0 - 0.25ds
# if (u <= 0) {
# continue;
# }
# if (v <= 0) {
# continue;
# }
for k in (-1, 1)
# const double cts = (D * cb - v) / u; // cosine scattering angle
cb = k * sqrt((1.0 + sb) * (1.0 - sb))
dcb = k * ds * sb / cb
# const double V =
# exp(-d * getInverseAttenuationLength(l_abs, ls, cts));
v = 0.5 * (d + D) * (d - D) / (d - D * cb)
u = d - v
# if (cts < 0.0 &&
# v * sqrt((1.0 + cts) * (1.0 - cts)) < MODULE_RADIUS_M) {
# continue;
# }
u <= 0 && continue
v <= 0 && continue
# const double W = std::min(A / (v * v), 2.0 * PI); // solid angle
# const double Ja = getScatteringProbability(cts); // d^2P/dcos/dϕ
# const double Jd = ng * (1.0 - cts) / C; // dt/du
cts = (D * cb - v) / u # cosine scattering angle
# const double ca = (D - v * cb) / u;
# const double sa = v * sb / u;
V = exp(
-d * inverseattenuationlength(
params.scattering_probability_model,
l_abs,
ls,
cts,
),
)
# const double dp = PI / phd.size();
# const double dom = dcb * dp * v * v / (u * u);
cts < 0.0 &&
v * sqrt((1.0 + cts) * (1.0 - cts)) < params.module_radius &&
continue
# for (const_iterator l = phd.begin(); l != phd.end(); ++l) {
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
# const double cp = l->getX();
# const double sp = l->getY();
ca = (D - v * cb) / u
sa = v * sb / u
# const double ct0 = cd * ca - sd * sa * cp;
dp = π / length(params.integration_points)
dom = dcb * dp * v^2 / (u^2)
# const double vx = sb * cp * qx;
# const double vy = sb * sp * qy;
# const double vz = cb * qz;
for (cp, sp) in params.integration_points.xy
ct0 = cd * ca - sd * sa * cp
# const double U =
# pmt.angular_acceptance(vx + vy + vz) +
# pmt.angular_acceptance(vx - vy + vz); // PMT angular acceptance
vx = sb * cp * qx
vy = sb * sp * qy
vz = cb * qz
# const double Jb = deltarayprobability(ct0); // d^2N/dcos/dϕ
U =
pmt.angular_acceptance(vx + vy + vz) +
pmt.angular_acceptance(vx - vy + vz) # PMT angular acceptance
# value += dom * npe * geanc() * dz * U * V * W * Ja * Jb * Jc /
# fabs(Jd);
# }
# }
# }
# }
# }
Jb = deltarayprobability(ct0) # d^2N/dcos/dϕ
# return value;
# }
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)
......
src/parameters.jl
+ 10
− 0
View file @ ab75cde8
struct PDFIntegrationPoints
xy::Vector{Tuple{Float64, Float64}}
function PDFIntegrationPoints(n::Int)
new([(cos(x), sin(x)) for x = (0.5π/n):(π/n):π])
end
end
Base.length(p::PDFIntegrationPoints) = length(p.xy)
"""
The parameter set for light detection.
......@@ -23,6 +31,7 @@ Base.@kwdef struct LMParameters{
module_radius::Float64 = 0.25
lambda_min::Float64 = 300.0
lambda_max::Float64 = 700.0
integration_points::PDFIntegrationPoints = PDFIntegrationPoints(20)
n::Float64 = 1.3800851282
legendre_coefficients::Tuple{Vector{Float64},Vector{Float64}} = gausslegendre(5)
dispersion_model::D = BaileyDispersion()
......@@ -35,6 +44,7 @@ function Base.show(io::IO, p::LMParameters)
println(io, " minimum distance = $(p.minimum_distance) m")
println(io, " module_radius = $(p.module_radius) m")
println(io, " lambda min / max = $(p.lambda_min) nm / $(p.lambda_max) nm")
println(io, " number of integrations points = $(length(p.integration_points))")
println(
io,
" degree of Legendre polynomials = $(length(first(p.legendre_coefficients)))",
......
test/runtests.jl
+ 15
− 5
View file @ ab75cde8
......@@ -2,7 +2,8 @@ using LumenManufaktur
using Test
@testset "direct light from muon" begin
params = LMParameters(dispersion_model = DispersionORCA)
params =
LMParameters(dispersion_model = DispersionORCA, integration_points = PDFIntegrationPoints(5))
pmt = LumenManufaktur.KM3NeTPMT
# The relative tolerances are based on comparisons between
......@@ -35,11 +36,20 @@ using Test
scatteredlightfrommuon(params, pmt, 2.0, 0.5, π / 2, π, 1.23),
rtol = 0.000001,
)
#
@test isapprox(
0.00432646,
directlightfromdeltarays(params, pmt, 1.5, π/2, π/3, 1.23),
rtol = 0.000001
directlightfromdeltarays(params, pmt, 1.5, π / 2, π / 3, 1.23),
rtol = 0.000001,
)
@test isapprox(
0.00283604,
scatteredlightfromdeltarays(params, pmt, 1.5, π / 2, π / 3, 1.23),
rtol = 0.00001,
)
@test isapprox(
0.000202688,
scatteredlightfromdeltarays(params, pmt, 2.5, π / 3, π / 4, 2.46),
rtol = 0.000001,
)
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment