Merge branch 'add-datx-support' into 'main'

Add datx support

See merge request !23

src/hardware.jl

@@ -14,6 +14,14 @@ struct PMT
     t₀::Float64
     status::Union{Int32, Missing}
 end
+function Base.isapprox(lhs::PMT, rhs::PMT; kwargs...)
+    for field in [:id, :status]
+        getfield(lhs, field) == getfield(rhs, field) || return false
+    end
+    for field in [:pos, :dir, :t₀]
+        isapprox(getfield(lhs, field), getfield(rhs, field); kwargs...) || return false
+    end
+end
 
 
 """
@@ -55,6 +63,18 @@ Base.length(d::DetectorModule) = d.n_pmts
 Base.eltype(::Type{DetectorModule}) = PMT
 Base.iterate(d::DetectorModule, state=1) = state > d.n_pmts ? nothing : (d.pmts[state], state+1)
 Base.isless(lhs::DetectorModule, rhs::DetectorModule) = lhs.location < rhs.location
+function Base.isapprox(lhs::DetectorModule, rhs::DetectorModule; kwargs...)
+    for field in [:id, :location, :n_pmts, :status]
+        getfield(lhs, field) == getfield(rhs, field) || return false
+    end
+    for field in [:pos, :q, :t₀]
+        isapprox(getfield(lhs, field), getfield(rhs, field); kwargs...) || return false
+    end
+    for (lhs_pmt, rhs_pmt) in zip(lhs.pmts, rhs.pmts)
+        isapprox(lhs_pmt, rhs_pmt; kwargs...)
+    end
+    true
+end
 """
     Base.getindex(d::DetectorModule, i) = d.pmts[i+1]
 
@@ -320,19 +340,81 @@ end
 
 Create a `Detector` instance from a DETX file.
 """
-function Detector(filename::AbstractString)
-    open(filename, "r") do fobj
-        Detector(fobj)
+function Detector(filename::AbstractString)::Detector
+    _, ext = splitext(filename)
+    if ext == ".detx"
+        return open(filename, "r") do fobj
+            read_detx(fobj)
+        end
+    elseif ext == ".datx"
+        return open(filename, "r") do fobj
+            read_datx(fobj)
+        end
     end
+    error("Unsupported detector file format '$(filename)'")
 end
 
 
 """
-    function Detector(io::IO)
+    function read_detx(io::IO)
+
+Create a `Detector` instance from an ASCII IO stream using the DATX specification.
+"""
+function read_datx(io::IO)
+    comment_marker = [0x23, 0x23, 0x23, 0x23]
+
+    comments = String[]
+    while read(io, 4) == comment_marker
+        push!(comments, _readstring(io))
+    end
+    seek(io, max(0, position(io) - length(comment_marker)))
+    det_id = read(io, Int32)
+    version = parse(Int, _readstring(io)[2:end])
+    validity = DateRange(unix2datetime(read(io, Float64)), unix2datetime(read(io, Float64)))
+    _readstring(io)  # says "UTM", ignoring
+    wgs = _readstring(io)
+    zone = _readstring(io)
+    utm_ref_grid = "$wgs $zone"
+    utm_position = UTMPosition(read(io, Float64), read(io, Float64), read(io, Float64))
+    n_modules = read(io, Int32)
+
+    modules = Dict{Int32, DetectorModule}()
+    locations = Dict{Tuple{Int, Int}, DetectorModule}()
+    strings = Int[]
+    for _ in 1:n_modules
+        module_id = read(io, Int32)
+        location = Location(read(io, Int32), read(io, Int32))
+        if !(location.string in strings)
+            push!(strings, location.string)
+        end
+        module_pos = Position{Float64}(read(io, Float64), read(io, Float64), read(io, Float64))
+        q = Quaternion{Float64}(read(io, Float64), read(io, Float64), read(io, Float64), read(io, Float64))
+        module_t₀ = read(io, Float64)
+        module_status = read(io, Int32)
+        n_pmts = read(io, Int32)
+        pmts = PMT[]
+        for _ in 1:n_pmts
+            pmt_id = read(io, Int32)
+            pmt_pos = Position{Float64}(read(io, Float64), read(io, Float64), read(io, Float64))
+            pmt_dir = Direction{Float64}(read(io, Float64), read(io, Float64), read(io, Float64))
+            pmt_t₀ = read(io, Float64)
+            pmt_status = read(io, Int32)
+            push!(pmts, PMT(pmt_id, pmt_pos, pmt_dir, pmt_t₀, pmt_status))
+        end
+        m = DetectorModule(module_id, module_pos, location, n_pmts, pmts, q, module_status, module_t₀)
+        modules[module_id] = m
+        locations[(location.string, location.floor)] = m
+    end
+    Detector(version, det_id, validity, utm_position, utm_ref_grid, n_modules, modules, locations, strings, comments)
+end
+@inline _readstring(io) = String(read(io, read(io, Int32)))
+
+"""
+    function read_detx(io::IO)
 
-Create a `Detector` instance from an IO stream.
+Create a `Detector` instance from an ASCII IO stream using the DETX specification.
 """
-function Detector(io::IO)
+function read_detx(io::IO)
     lines = readlines(io)
 
     comments = _extract_comments(lines, DETECTOR_COMMENT_PREFIX)
@@ -421,15 +503,7 @@ function Detector(io::IO)
             push!(pmts, PMT(pmt_id, Position(x, y, z), Direction(dx, dy, dz), t0, pmt_status))
         end
 
-        if (ismissing(t₀) || t₀ == 0.0) && floor > 0
-            # t₀ is only available in DETX v4+ and even with supported versions, the value is
-            # sometimes 0 when e.g. the DETX was converted with Jpp from a version which did
-            # not include that informatino (v3 and below). Here, we are using the averaged
-            # PMT t₀s for the module t₀, just like Jpp does nowadays.
-            t₀ = mean([pmt.t₀ for pmt in pmts])
-        end
-
-        # If t₀ is still missing, we default to 0.0
+        # If t₀ is missing, we default to 0.0
         if ismissing(t₀)
             t₀ = 0.0
         end

test/hardware.jl

@@ -81,7 +81,12 @@ const SAMPLES_DIR = joinpath(@__DIR__, "samples")
             @test m.location.floor in 1:4
         end
 
-        @test 478392.31980645156 ≈ d.modules[808992603].t₀
+        # This used to be calculated from the mean PMT positions if the value was missing or 0.0
+        # from KM3io.jl v0.14.9 we will show the value as it is and may implement something like
+        # gett₀(m::DetectorModule) in future to calculate the fallback value
+        #
+        # @test 478392.31980645156 ≈ d.modules[808992603].t₀
+        @test 0.0 ≈ d.modules[808992603].t₀
 
         @test 9 == d.strings[1]
         @test 30 == d.strings[end]
@@ -120,7 +125,7 @@ end
 @testset "DETX writing" begin
     for from_version ∈ 1:5
         for to_version ∈ 1:5
-            out = tempname()
+            out = tempname(;cleanup=false) * ".detx"
             d₀ = Detector(joinpath(SAMPLES_DIR, "v$(from_version).detx"))
             write(out, d₀; version=to_version)
             d = Detector(out)
@@ -133,9 +138,26 @@ end
                     end
                 end
             end
+            rm(out)
         end
     end
 end
+@testset "DATX" begin
+    detx = Detector(datapath("detx", "KM3NeT_00000133_20221025.detx"))
+    datx = Detector(datapath("datx", "KM3NeT_00000133_20221025.datx"))
+    for field in fieldnames(Detector)
+        field == :modules && continue
+        if field == :locations
+            detx_locs = getfield(detx, field)
+            datx_locs = getfield(datx, field)
+            for key in keys(detx_locs)
+                @test isapprox(detx_locs[key], datx_locs[key]; atol=1e-06)
+            end
+            continue
+        end
+        @test getfield(detx, field) == getfield(datx, field)
+    end
+end
 @testset "hydrophones" begin
     hydrophones = read(joinpath(SAMPLES_DIR, "hydrophone.txt"), Hydrophone)
     @test 19 == length(hydrophones)