diff --git a/Makefile b/Makefile index e92b604e59dd6f0949f81327ece7b726cd3ac980..a0f54d0c9d9b1c76d8fbffe1f4eb62b0df1924c5 100644 --- a/Makefile +++ b/Makefile @@ -10,8 +10,11 @@ test: clean: rm -rf docs/build/ +docwatch: + while true; do fswatch -1 docs src && make doc; done + preview: julia -e 'using LiveServer; serve(dir="docs/build")' -.PHONY: build doc test clean preview +.PHONY: build doc docwatch test clean preview diff --git a/Project.toml b/Project.toml index f9e797d048d158420432af0e7f878540087f3b15..6fdc5ef7c9d4e07396a0fd36c6011ace74cfeaab 100644 --- a/Project.toml +++ b/Project.toml @@ -19,7 +19,7 @@ UnROOT = "3cd96dde-e98d-4713-81e9-a4a1b0235ce9" [compat] DocStringExtensions = "0.8, 0.9" HDF5 = "^0.16.15, ^0.17" -KM3NeTTestData = "^0.4.12" +KM3NeTTestData = "^0.4.14" StaticArrays = "1" UnROOT = "^0.10.26" julia = "1" diff --git a/docs/Project.toml b/docs/Project.toml index 8be78c5c0ebcb4294c6272647962219c57dab232..1841cbae36c44ba32375ba35dcdc01e6bbd2b8ef 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -10,3 +10,4 @@ PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925" [compat] Documenter = "1" FHist = "^0.11" +KM3NeTTestData = "^0.4.14" diff --git a/docs/src/api.md b/docs/src/api.md index 48b315016dcbcb80f4a5b0a8b3f8219df7b9b32c..deb8ec5fc2797b5e47daf90ceeeb8f5b86816c25 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -93,8 +93,11 @@ read(filename::AbstractString, T::Type{AcousticsTriggerParameter}) calibrate calibratetime combine +Orientations +Compass floordist slew +slerp ``` ## Physics diff --git a/docs/src/manual/calibration.md b/docs/src/manual/calibration.md index f0a5deec1050b535a891d3bb405cf42674d53c17..4323e38618540be11b6f5cb431f6d07ed5eadea1 100644 --- a/docs/src/manual/calibration.md +++ b/docs/src/manual/calibration.md @@ -1,5 +1,46 @@ # Calibration -It's implemented but not documented here yet. Check out the docs `calibrate` -function's docstring! +In general, data related to detector calibration are either stored in the +database or in the [Calibration +Archive](https://git.km3net.de/auxiliary_data/calibration). +## Hits + +Hit calibration is implemented but not documented here yet. Check out the docs for [`calibrate()`](@ref). + +## Module Orientations + +Module orientation data is stored in ROOT files which are generated by the Jpp +framework. The corresponding ROOT files, which are the output of the [dynamic +calibration +procedure](https://common.pages.km3net.de/jpp/Position_calibration.PDF), are +stored in the [Calibration +Archive](https://git.km3net.de/auxiliary_data/calibration) under +`orientations/`. + +`KM3io.jl` extends the `Base.read` function with a method which reads the whole +orientations file at once in an object of type [`Orientations`](@ref). This +object can be called to calculate the orientation of a module (as a quaternion) +for a given time, as long as the time is within the time range of the +orientation data. The quaternions from the orientation data are interpolated for +the given time using [`slerp()`](@ref). + +The following example shows how to read the orientation data and obtain the +orientation quaternion for a module at a given time. + +```@example 1 +using KM3io, KM3NeTTestData + +o = read(datapath("calib", "KM3NeT_00000133_D_1.0.0_00017397_00017496_1.orientations.root"), Orientations) + +module_id = 817589211 + +q = o(module_id, 1693408347) +``` + +The quaternions can be converted to [`Compass`](@ref) object which has the +fields `.yaw`, `.pitch` and `.roll`: + +```@example 1 +compass = Compass(q) +``` diff --git a/src/KM3io.jl b/src/KM3io.jl index 4ad54bf04495397f3c59d2c22d3e18e26be3cd99..330cb87be1433cead3e7e6f118d7a0a4b4fae186 100644 --- a/src/KM3io.jl +++ b/src/KM3io.jl @@ -59,6 +59,7 @@ include("hardware.jl") include("root/online.jl") include("root/offline.jl") include("root/root.jl") +include("root/calibration.jl") include("hdf5/hdf5.jl") include("daq.jl") include("acoustics.jl") diff --git a/src/exports.jl b/src/exports.jl index 0446c61cd78f9290feb9754107022a4a332075ca..ac9cab1c17d3d1d7af5b83e4f09f1d29bd49eaf6 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -84,6 +84,8 @@ calibratetime, combine, floordist, slew, +Orientations, +Compass, # Reconstruction RecStageRange, @@ -127,6 +129,7 @@ distance, phi, theta, zenith, +slerp, # Real-time @ip_str, diff --git a/src/root/calibration.jl b/src/root/calibration.jl new file mode 100644 index 0000000000000000000000000000000000000000..e75fd1d960d2570d96d9d0b5168d133ca17b01f1 --- /dev/null +++ b/src/root/calibration.jl @@ -0,0 +1,97 @@ +""" + +A data structure to hold orientations data. This struct should be instantiated by +`Base.read(filename, Orientations)`. + +""" +struct Orientations + module_ids::Set{Int} + times::Dict{Int, Vector{Float64}} + quaternions::Dict{Int, Vector{Quaternion}} +end +function Base.show(io::IO, o::Orientations) + min_t = Inf + max_t = -Inf + for times in values(o.times) + _min_t, _max_t = extrema(times) + if _min_t < min_t + min_t = _min_t + end + if _max_t > max_t + max_t = _max_t + end + end + tâ‚ = unix2datetime(min_t) + tâ‚‚ = unix2datetime(max_t) + print(io, "Orientations of $(length(o.module_ids)) modules ($(tâ‚) to $(tâ‚‚))") +end +function (o::Orientations)(module_id::Integer, time::Real) + times = o.times[module_id] + (time < first(times) || time > last(times)) && error("The requested time is outside of the range of the orientations data.") + + idx2 = searchsortedfirst(times, time) + q2 = o.quaternions[module_id][idx2] + (idx2 >= length(times) || idx2 == 1) && return q2 + idx1 = idx2 - 1 + q1 = o.quaternions[module_id][idx1] + + t1 = times[idx1] + t2 = times[idx2] + Δt = t2 - t1 + + Δt == 0.0 && return q1 + + t = (time - t1) / Δt + + slerp(q1, q2, t) +end +(o::Orientations)(module_id::Integer) = (t=o.times[module_id], q=o.quaternions[module_id]) + + +function Base.read(filename::AbstractString, T::Type{Orientations}) + f = UnROOT.ROOTFile(filename) + module_ids = Set{Int}() + quaternions = Dict{Int, Vector{Quaternion}}() + times = Dict{Int, Vector{Float64}}() + for (module_id, t, a, b, c, d) in zip([UnROOT.LazyBranch(f, "ORIENTATION/ORIENTATION/$(b)") for b in ["id", "t", "JCOMPASS::JQuaternion/a", "JCOMPASS::JQuaternion/b", "JCOMPASS::JQuaternion/c", "JCOMPASS::JQuaternion/d"]]...) + if !(module_id ∈ module_ids) + push!(module_ids, module_id) + quaternions[module_id] = Quaternion[] + times[module_id] = Float64[] + end + push!(quaternions[module_id], Quaternion(a, b, c, d)) + push!(times[module_id], t) + end + T(module_ids, times, quaternions) +end + +""" +A compass with yaw, pitch and roll. +""" +struct Compass + yaw::Float64 + pitch::Float64 + roll::Float64 +end + +""" + +Initialises a [`Compass`](@ref) from a [`Quaternion`](@ref). + +""" +function Compass(q::Quaternion) + yaw = -atan(2.0 * (q.q0 * q.qz + q.qx * q.qy), 1.0 - 2.0 * (q.qy * q.qy + q.qz * q.qz)) + sp = 2.0 * (q.q0 * q.qy - q.qz * q.qx) + + if (sp >= +1.0) + pitch = asin(+1.0) + elseif (sp <= -1.0) + pitch = asin(-1.0) + else + pitch = asin(sp) + end + + roll = -atan(2.0 * (q.q0 * q.qx + q.qy * q.qz), 1.0 - 2.0 * (q.qx * q.qx + q.qy * q.qy)) + + Compass(yaw, pitch, roll) +end diff --git a/src/tools/math.jl b/src/tools/math.jl index 0a41b685165342096f000a7ee5a7d3251ca50b94..a83d243e60476687acb192833d9c69c4790ef50e 100644 --- a/src/tools/math.jl +++ b/src/tools/math.jl @@ -12,3 +12,53 @@ Calculates the disance between two points. """ distance(a::Position, b::Position) = norm(a - b) + + +""" + +Interpolate between two vectors (e.g. quaternions) using the slerp method. `t` +should be between 0 and 1. 0 will produce `qâ‚` and `1` `qâ‚‚`. + +The input vectors `qâ‚` and `qâ‚‚` will be normalised unless `normalized` is +`false`. It is not done by default to shave off a few dozens of nanoseconds. +Make sure to set `normalized=false` if the input vectors are not unit vectors. + +""" +function slerp(qâ‚, qâ‚‚, t::Real; dot_threshold=0.9995, normalized=true) + if !normalized + qâ‚ = normalize(qâ‚) + qâ‚‚ = normalize(qâ‚‚) + end + + dot = qâ‚â‹…qâ‚‚ + + if (dot < 0.0) + qâ‚‚ *= -1 + dot = -dot + end + + sâ‚ = t + sâ‚€ = 1.0 - t + + if dot <= dot_threshold + θ₀ = acos(dot) + θ₠= θ₀ * t + + sâ‚ = sin(θâ‚) / sin(θ₀) + sâ‚€ = cos(θâ‚) - dot * sâ‚ + end + + normalize((sâ‚€ * qâ‚) + (sâ‚ * qâ‚‚)) +end + +# Another implementation which yields slightly different results. +# Further reading: http://number-none.com/product/Understanding%20Slerp,%20Then%20Not%20Using%20It/ +# +# function slerp(qâ‚, qâ‚‚, t::Real; dot_threshold=0.9995) +# dot = acos(qâ‚â‹…qâ‚‚) +# dot > dot_threshold && return normalize(qâ‚ + t*(qâ‚‚ - qâ‚)) +# dot = clamp(dot, -1, 1) +# θ = acos(dot) * t +# q = normalize(qâ‚‚ - qâ‚*dot) +# qâ‚*cos(θ) + q*sin(θ) +# end diff --git a/test/calibration.jl b/test/calibration.jl index f151a58556eb7a768ceaf5ba606255e084c9f335..3296b3b289c74814809fb0dbaddc7db64eaad7c5 100644 --- a/test/calibration.jl +++ b/test/calibration.jl @@ -22,3 +22,31 @@ end det = Detector(datapath("detx", "km3net_offline.detx")) @test 9.61317647058823 ≈ floordist(det) end + + +@testset "orientations" begin + o = read(datapath("calib", "KM3NeT_00000133_D_1.0.0_00017397_00017496_1.orientations.root"), Orientations) + @show o + module_id = 817589211 + min_t, max_t = extrema(o.times[module_id]) + Δt = max_t - min_t + @test [0.8260205110995139, 0.003912907129683348, -0.004395551387888641, -0.5636093359133512] == o(module_id, min_t) + @test [0.8289446524907407, 0.004590185819553083, -0.0007479055911552097, -0.5593113032456739] == o(module_id, max_t) + @test [0.8297266567631056, 0.002991865243189534, -0.004798371006076004, -0.5581412898494195] == o(module_id, min_t + Δt/2) + @test [0.8305219131347711, 0.003947997911424212, -0.0042572917986734805, -0.5569556899628482] == o(module_id, min_t + Δt/3) + + qdata = o(module_id) + @test 5 == length(qdata.t) + @test 5 == length(qdata.q) + @test 1.693407821152e9 == qdata.t[1] + @test [0.8260205110995139, 0.003912907129683348, -0.004395551387888641, -0.5636093359133512] == qdata.q[1] +end + +@testset "Compass" begin + q = Quaternion(0.8260205110995139, 0.003912907129683348, -0.004395551387888641, -0.5636093359133512) + c = Compass(q) + @test c.yaw == 1.1975374212207646 + @test c.pitch == -0.0028509330922497 + @test c.roll == -0.011419325278029469 + +end diff --git a/test/math.jl b/test/math.jl new file mode 100644 index 0000000000000000000000000000000000000000..c4ff0cbd32ad3bbfed24374a85315e532dd60c84 --- /dev/null +++ b/test/math.jl @@ -0,0 +1,29 @@ +import KM3io: slerp +using Test + +@testset "slerp()" begin + q1 = [1.0, 0.0] + q2 = [0.0, 1.0] + @test q1 ≈ slerp(q1, q2, 0) + @test q2 ≈ slerp(q1, q2, 1) + @test [0.9510565162951538, 0.30901699437494745] ≈ slerp(q1, q2, 0.2) + @test [0.70710678, 0.70710678] ≈ slerp(q1, q2, 0.5) + @test [0.45399049973954686, 0.8910065241883678] ≈ slerp(q1, q2, 0.7) + + # should normalise internally + q1 = [0.4, 0.0] + q2 = [0.0, 0.9] + @test [1.0, 0.0] ≈ slerp(q1, q2, 0; normalized=false) + @test [0.0, 1.0] ≈ slerp(q1, q2, 1; normalized=false) + @test [0.9510565162951538, 0.30901699437494745] ≈ slerp(q1, q2, 0.2; normalized=false) + @test [0.70710678, 0.70710678] ≈ slerp(q1, q2, 0.5; normalized=false) + @test [0.45399049973954686, 0.8910065241883678] ≈ slerp(q1, q2, 0.7; normalized=false) + + q1 = [1.0, 0.0, 0.0] + q2 = [0.0, 0.0, 1.0] + @test q1 ≈ slerp(q1, q2, 0) + @test q2 ≈ slerp(q1, q2, 1) + @test [0.9510565162951538, 0.0, 0.30901699437494745] ≈ slerp(q1, q2, 0.2) + @test [0.70710678, 0.0, 0.70710678] ≈ slerp(q1, q2, 0.5) + @test [0.45399049973954686, 0.0, 0.8910065241883678] ≈ slerp(q1, q2, 0.7) +end diff --git a/test/runtests.jl b/test/runtests.jl index 2609cfd9e78f992cb841ad54293f625bda0b53c9..10c19d7fb4100cb7f2de661795435a325ac860fa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,4 +8,5 @@ include("hardware.jl") include("acoustics.jl") include("calibration.jl") include("physics.jl") +include("math.jl") include("controlhost.jl")