Merge branch '31-add-support-to-oscillations-open-data-release' into 'main'

Resolve "Add support to oscillations open data release"

Closes #31

See merge request !47

.gitignore

@@ -5,4 +5,4 @@
 /Manifest.toml
 /dev/
 /docs/build/
-/docs/site/
+/docs/site/
\ No newline at end of file

Project.toml

@@ -4,6 +4,7 @@ authors = ["Tamas Gal", "Johannes Schumann"]
 version = "0.17.10"
 
 [deps]
+Corpuscles = "e78a0372-c628-4773-9c8d-eb17d149bf93"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
@@ -25,10 +26,11 @@ KM3DB = "a9013879-bb44-4449-9e5b-40f9ac008ab0"
 KM3ioKM3DBExt = "KM3DB"
 
 [compat]
+Corpuscles = "^2.1.2"
 DocStringExtensions = "0.8, 0.9"
 HDF5 = "^0.16.15, ^0.17"
 KM3DB = "0.2.3"
-KM3NeTTestData = "^0.4.17"
+KM3NeTTestData = "^0.4.18"
 StaticArrays = "1"
 UnROOT = "^0.10.26"
 julia = "1"

docs/src/api.md

@@ -38,6 +38,17 @@ SummarysliceHeader
 SummaryFrame
 ```
 
+## Oscillations Open Data
+```@docs
+OscillationsData
+OSCFile
+ResponseMatrixBin
+ResponseMatrixBinNeutrinos
+ResponseMatrixBinMuons
+ResponseMatrixBinData
+OscOpenDataTree
+```
+
 ## HDF5
 ```@docs
 H5File
@@ -172,4 +183,4 @@ hashistory
 ```@docs
 angle
 distance
-```
+```
\ No newline at end of file

docs/src/manual/rootfiles.md

@@ -301,3 +301,44 @@ ROOTFile{OnlineTree (136335 events, 107632 summaryslices)}
 
 Now you can use it as if it was on your local filesystem. `UnROOT.jl` will take
 care of loading only the needed data from the server.
+
+## [Oscillations Open Dataformat](@id oscillations dataformat)
+
+The [oscillations
+dataformat](https://git.km3net.de/vcarretero/prepare-opendata-orca6-433kty/-/tree/main?ref_type=heads)
+is used to store the responses from a particular oscillations analysis data release. The
+`OSCFile` type represents an actual ROOT file and it is essentially a
+vector of Response like entries (`Vector{ResponseMatrixBin}`) . Depending on what is stored in the initial ROOT file, neutrinos, data and muons response trees are accessible via the `.osc_opendata_nu`, `.osc_opendata_data` and `.osc_opendata_muons`  fields of the `ROOTFile` type respectively.
+
+### ResponseMatrixBin
+
+The `ResponseMatrixBin` stores individual directions of a bin in order to fill a histogram.
+
+``` julia-repl
+julia> using KM3io, KM3NeTTestData
+
+julia> f = OSCFile(datapath("oscillations", "ORCA6_433kt-y_opendata_v0.4_testdata.root"))
+OSCFile{OscOpenDataTree of Neutrinos (59301 events), OscOpenDataTree of Data (106 events), OscOpenDataTree of Muons (99 events)}
+
+julia> f.osc_opendata_nu
+OscOpenDataTree (59301 events)
+
+julia> f.osc_opendata_nu[1]
+KM3io.ResponseMatrixBinNeutrinos(10, 1, 30, 18, -12, 1, 1, 52.25311519561337, 2730.388047646041)
+
+julia> dump(f.osc_opendata_nu[1])
+KM3io.ResponseMatrixBinNeutrinos
+  E_reco_bin: Int64 10
+  Ct_reco_bin: Int64 1
+  E_true_bin: Int64 30
+  Ct_true_bin: Int64 18
+  Flav: Int16 -12
+  IsCC: Int16 1
+  AnaClass: Int16 1
+  W: Float64 52.25311519561337
+  Werr: Float64 2730.388047646041
+
+julia> f.osc_opendata_data[1]
+KM3io.ResponseMatrixBinData(2, 6, 1, 2.0)
+
+```

src/KM3io.jl

@@ -27,6 +27,7 @@ using StaticArrays: FieldVector, @SArray, SVector
 import UnROOT
 
 using HDF5
+using Corpuscles
 
 
 include("exports.jl")
@@ -63,6 +64,7 @@ include("types.jl")
 include("hardware.jl")
 include("root/online.jl")
 include("root/offline.jl")
+include("root/oscillations.jl")
 include("root/root.jl")
 include("root/calibration.jl")
 include("hdf5/hdf5.jl")

src/exports.jl

@@ -80,6 +80,15 @@ Trk,
 XCalibratedHit,
 FitInformation,
 
+# Oscillations open data
+OscillationsData,
+OSCFile,
+ResponseMatrixBin,
+ResponseMatrixBinNeutrinos,
+ResponseMatrixBinMuons,
+ResponseMatrixBinData,
+OscOpenDataTree,
+
 # Misc I/O
 tojson,
 

src/root/oscillations.jl

+"""
+`OscillationsData` is an abstract type representing the data in an oscillation open data file.
+
+"""
+abstract type OscillationsData end
+
+"""
+
+`OSCFile` is a structure representing an oscillation open data file. Depending on the trees inside the root file it will have different fields (neutrino, muons, data).
+
+"""
+struct OSCFile
+    _fobj::Union{UnROOT.ROOTFile,Dict}
+    rawroot::Union{UnROOT.ROOTFile,Nothing}
+    osc_opendata_nu::Union{OscillationsData,Nothing}
+    osc_opendata_data::Union{OscillationsData,Nothing}
+    osc_opendata_muons::Union{OscillationsData,Nothing}
+
+    function OSCFile(filename::AbstractString)
+        if endswith(filename, ".root")
+            fobj = UnROOT.ROOTFile(filename)
+            osc_opendata_nu = ROOT.TTREE_OSC_OPENDATA_NU ∈ keys(fobj) ? OscOpenDataTree(fobj, ROOT.TTREE_OSC_OPENDATA_NU) : nothing
+            osc_opendata_data = ROOT.TTREE_OSC_OPENDATA_DATA ∈ keys(fobj) ? OscOpenDataTree(fobj, ROOT.TTREE_OSC_OPENDATA_DATA) : nothing
+            osc_opendata_muons = ROOT.TTREE_OSC_OPENDATA_MUONS ∈ keys(fobj) ? OscOpenDataTree(fobj, ROOT.TTREE_OSC_OPENDATA_MUONS) : nothing
+            return new(fobj, fobj, osc_opendata_nu, osc_opendata_data, osc_opendata_muons)
+        end
+    end
+end
+Base.close(f::OSCFile) = close(f._fobj)
+function Base.show(io::IO, f::OSCFile)
+    if isa(f._fobj, UnROOT.ROOTFile)
+        s = String[]
+        !isnothing(f.osc_opendata_nu) && push!(s, "$(f.osc_opendata_nu)")
+        !isnothing(f.osc_opendata_data) && push!(s, "$(f.osc_opendata_data)")
+        !isnothing(f.osc_opendata_muons) && push!(s, "$(f.osc_opendata_muons)")
+        info = join(s, ", ")
+        print(io, "OSCFile{$info}")
+    else
+        print(io, "Empty OSCFile")
+    end
+
+end
+
+const NUE_PDGID = Particle("nu(e)0").pdgid.value
+const ANUE_PDGID = Particle("~nu(e)0").pdgid.value
+const NUMU_PDGID = Particle("nu(mu)0").pdgid.value
+const ANUMU_PDGID = Particle("~nu(mu)0").pdgid.value
+const NUTAU_PDGID = Particle("nu(tau)0").pdgid.value
+const ANUTAU_PDGID = Particle("~nu(tau)0").pdgid.value
+
+"""
+
+A `ResponseMatrixBin` is an abstract type representing a bin in a response matrix.
+
+"""
+abstract type ResponseMatrixBin end
+
+"""
+
+A concrete type representing a response matrix bin for neutrino events.
+
+"""
+struct ResponseMatrixBinNeutrinos <: ResponseMatrixBin
+    E_reco_bin::Int64
+    Ct_reco_bin::Int64
+    E_true_bin::Int64
+    Ct_true_bin::Int64
+    Flav::Int16
+    IsCC::Int16
+    AnaClass::Int16
+    W::Float64
+    Werr::Float64
+end
+
+"""
+
+A concrete type representing a response matrix bin for muon events. There is no true quantities for muon events.
+
+"""
+struct ResponseMatrixBinMuons <: ResponseMatrixBin
+    E_reco_bin::Int64
+    Ct_reco_bin::Int64
+    AnaClass::Int16
+    W::Float64
+    Werr::Float64
+end
+
+"""
+
+A concrete type representing a response matrix bin for data events. There is no true quantities for data events.
+
+"""
+struct ResponseMatrixBinData <: ResponseMatrixBin
+    E_reco_bin::Int64
+    Ct_reco_bin::Int64
+    AnaClass::Int16
+    W::Float64
+end
+
+
+"""
+
+Function to get the PDG ID for a given neutrino flavor and neutrino/antineutrino flag.
+
+"""
+function _getpdgnumber(flav::Integer, isNB::Integer)
+    flav == 0 && isNB == 0 && return NUE_PDGID # nu(e)0
+    flav == 0 && isNB == 1 && return ANUE_PDGID # ~nu(e)0
+    flav == 1 && isNB == 0 && return NUMU_PDGID # nu(mu)0
+    flav == 1 && isNB == 1 && return ANUMU_PDGID # ~nu(mu)0
+    flav == 2 && isNB == 0 && return NUTAU_PDGID # nu(tau)0
+    flav == 2 && isNB == 1 && return ANUTAU_PDGID # ~nu(tau)0
+
+    error("Invalid flavor: $flav($isNB)")
+end
+
+"""
+
+Function to get the name of the analysis class based on its identifier.
+
+"""
+function _getanaclassname(fClass::Integer)
+    fClass == 1 && return "HighPurityTracks"
+    fClass == 2 && return "Showers"
+    fClass == 3 && return "LowPurityTracks"
+    error("Invalid class: $fClass)")
+end
+
+"""
+
+`OscOpenDataTree` is a structure representing an oscillation open data tree, it will be represented as response functions.
+
+"""
+struct OscOpenDataTree{T} <: OscillationsData
+    _fobj::UnROOT.ROOTFile
+    #header::Union{MCHeader, Missing} # no header for now, subject to change
+    _bin_lookup_map::Dict{Tuple{Int,Int,Int},Int} # Not implemented for now
+    _t::T  # carry the type to ensure type-safety
+    tpath::String
+
+    function OscOpenDataTree(fobj::UnROOT.ROOTFile, tpath::String)
+        if tpath == ROOT.TTREE_OSC_OPENDATA_NU
+            branch_paths = ["E_reco_bin",
+                "Ct_reco_bin",
+                "Flav",
+                "IsCC",
+                "IsNB",
+                "E_true_bin",
+                "Ct_true_bin",
+                "W",
+                "WE",
+                "Class",
+            ]
+
+        elseif tpath == ROOT.TTREE_OSC_OPENDATA_DATA
+            branch_paths = ["E_reco_bin",
+            "Ct_reco_bin",
+            "W",
+            "Class",
+            ]
+        elseif tpath == ROOT.TTREE_OSC_OPENDATA_MUONS 
+            branch_paths = ["E_reco_bin",
+            "Ct_reco_bin",
+            "W",
+            "WE",
+            "Class",
+            ]
+        end
+
+
+        t = UnROOT.LazyTree(fobj, tpath, branch_paths)
+
+        new{typeof(t)}(fobj, Dict{Tuple{Int,Int,Int},Int}(), t, tpath)
+    end
+end
+
+"""
+
+Construct an `OscOpenDataTree` from a ROOT file and a tree path.
+
+"""
+OscOpenDataTree(filename::AbstractString, tpath::String) = OscOpenDataTree(UnROOT.ROOTFile(filename), tpath)
+
+Base.close(f::OscOpenDataTree) = close(f._fobj)
+Base.length(f::OscOpenDataTree) = length(f._t)
+Base.firstindex(f::OscOpenDataTree) = 1
+Base.lastindex(f::OscOpenDataTree) = length(f)
+Base.eltype(::OscOpenDataTree) = ResponseMatrixBin
+function Base.iterate(f::OscOpenDataTree, state=1)
+    state > length(f) ? nothing : (f[state], state + 1)
+end
+function Base.show(io::IO, f::OscOpenDataTree)
+    data_name = f.tpath == ROOT.TTREE_OSC_OPENDATA_NU ? "Neutrinos" :
+                f.tpath == ROOT.TTREE_OSC_OPENDATA_DATA ? "Data" :
+                f.tpath == ROOT.TTREE_OSC_OPENDATA_MUONS ? "Muons" : "Unknown"
+    print(io, "OscOpenDataTree of $(data_name) ($(length(f)) events)")
+end
+
+Base.getindex(f::OscOpenDataTree, r::UnitRange) = [f[idx] for idx ∈ r]
+Base.getindex(f::OscOpenDataTree, mask::BitArray) = [f[idx] for (idx, selected) ∈ enumerate(mask) if selected]
+function Base.getindex(f::OscOpenDataTree, idx::Integer)
+    if idx > length(f)
+        throw(BoundsError(f, idx))
+    end
+    idx > length(f) && throw(BoundsError(f, idx))
+    e = f._t[idx]  # the event as NamedTuple: struct of arrays
+
+    if f.tpath == ROOT.TTREE_OSC_OPENDATA_NU
+        ResponseMatrixBinNeutrinos(
+            e.E_reco_bin,
+            e.Ct_reco_bin,
+            e.E_true_bin,
+            e.Ct_true_bin,
+            _getpdgnumber(e.Flav, e.IsNB),
+            e.IsCC,
+            e.Class,
+            e.W,
+            e.WE)
+    elseif f.tpath == ROOT.TTREE_OSC_OPENDATA_MUONS
+        ResponseMatrixBinMuons(
+            e.E_reco_bin,
+            e.Ct_reco_bin,
+            e.Class,
+            e.W,
+            e.WE)
+    elseif f.tpath == ROOT.TTREE_OSC_OPENDATA_DATA
+        ResponseMatrixBinData(
+            e.E_reco_bin,
+            e.Ct_reco_bin,
+            e.Class,
+            e.W)
+    end
+
+end
+

src/tools/helpers.jl

@@ -22,6 +22,7 @@ end
 
 countevents(tree::OfflineTree) = length(tree)
 countevents(tree::OnlineTree) = length(tree.events)
+countevents(tree::OscillationsData) = length(tree)
 triggercounterof(e::Evt) = e.trigger_counter
 frameindexof(e::Evt) = e.frame_index
 triggercounterof(e::DAQEvent) = e.header.trigger_counter
@@ -58,6 +59,7 @@ function getevent(tree::T, frame_index, trigger_counter) where T<:Union{OnlineTr
 end
 getevent(tree::OfflineTree, idx) = tree[idx]
 getevent(tree::OnlineTree, idx) = tree.events[idx]
+getevent(tree::OscillationsData, idx) = tree[idx]
 
 
 """

test/oscillations.jl

+using KM3io
+import UnROOT
+using KM3NeTTestData
+using Test
+
+const OSCFILE = datapath("oscillations", "ORCA6_433kt-y_opendata_v0.4_testdata.root")
+
+@testset "Oscillations open data files" begin
+    f = OSCFile(OSCFILE)
+    nu = f.osc_opendata_nu
+    data = f.osc_opendata_data
+    muons = f.osc_opendata_muons
+    @test 59301 == length(nu)
+    @test 1 == nu[1].AnaClass
+    @test 1 == nu[1].Ct_reco_bin
+    @test 18 == nu[1].Ct_true_bin
+    @test 10 == nu[1].E_reco_bin
+    @test 30 == nu[1].E_true_bin
+    @test -12 == nu[1].Flav
+    @test 1 == nu[1].IsCC
+    @test isapprox(nu[1].W, 52.25311519561337)
+    @test isapprox(nu[1].Werr, 2730.388047646041)
+
+    @test 106 == length(data)
+    @test 1 == data[1].AnaClass
+    @test 6 == data[1].Ct_reco_bin
+    @test 2 == data[1].E_reco_bin
+    @test isapprox(data[1].W, 2.0)
+
+    @test 99 == length(muons)
+    @test 1 == muons[1].AnaClass
+    @test 4 == muons[1].Ct_reco_bin
+    @test 1 == muons[1].E_reco_bin
+    @test isapprox(data[1].W, 2.0)
+
+    close(f)
+end
\ No newline at end of file

test/runtests.jl

@@ -11,3 +11,4 @@ include("calibration.jl")
 include("physics.jl")
 include("math.jl")
 include("controlhost.jl")
+include("oscillations.jl")