diff --git a/.gitignore b/.gitignore index 99e2ceab4c08aead66c9cd62d6de36d11df14adb..d0f623c2f8cb376944184b7c3e1059c70defe5d1 100644 --- a/.gitignore +++ b/.gitignore @@ -5,4 +5,4 @@ /Manifest.toml /dev/ /docs/build/ -/docs/site/ +/docs/site/ \ No newline at end of file diff --git a/Project.toml b/Project.toml index 3a9a447eb5638cc53416902d0d9e0e181050da0c..6bdb3a8d71ec3fde9c2db24722eedaef1afa21aa 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/api.md b/docs/src/api.md index 540cf7f3491a7ea3154dedc6f27fc172cb592b6c..6331983141aafa297596a060d5dca08f32d5a640 100644 --- a/docs/src/api.md +++ b/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 diff --git a/docs/src/manual/rootfiles.md b/docs/src/manual/rootfiles.md index 90b28157ed37dc9f54667e8953c8503223618333..71ea1302c0d76572f08368eb537ac2fdfe2eb480 100644 --- a/docs/src/manual/rootfiles.md +++ b/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) + +``` diff --git a/src/KM3io.jl b/src/KM3io.jl index 3b6a91b8c7df52fe2cad276b6b92fc334ddbc75c..fd1c11e76af3b22145a90c76dea4032caeecaddb 100644 --- a/src/KM3io.jl +++ b/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") diff --git a/src/exports.jl b/src/exports.jl index 643187f455c8d41b84dc48f01840d424f2cff144..d531fdf6e0152b0c85b66e1bfff6ef14d39fd382 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -80,6 +80,15 @@ Trk, XCalibratedHit, FitInformation, +# Oscillations open data +OscillationsData, +OSCFile, +ResponseMatrixBin, +ResponseMatrixBinNeutrinos, +ResponseMatrixBinMuons, +ResponseMatrixBinData, +OscOpenDataTree, + # Misc I/O tojson, diff --git a/src/root/oscillations.jl b/src/root/oscillations.jl new file mode 100644 index 0000000000000000000000000000000000000000..f98b30113c0f70d5a04c800da64da8dd97610999 --- /dev/null +++ b/src/root/oscillations.jl @@ -0,0 +1,235 @@ +""" +`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 + diff --git a/src/tools/helpers.jl b/src/tools/helpers.jl index be09c1d39b9ccfe93d810ee3540a90df2132dd7e..db0f9dc2e12e77325957df8ed7933924eda8c224 100644 --- a/src/tools/helpers.jl +++ b/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] """ diff --git a/test/oscillations.jl b/test/oscillations.jl new file mode 100644 index 0000000000000000000000000000000000000000..3b2fd3ce4c51ff00741724fa95ef145be96da023 --- /dev/null +++ b/test/oscillations.jl @@ -0,0 +1,37 @@ +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 diff --git a/test/runtests.jl b/test/runtests.jl index 8a277144e786c0f4137035019a4c6900861eddee..4f67d5b79822079ca2194c00a3d9774dfb849b68 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,3 +11,4 @@ include("calibration.jl") include("physics.jl") include("math.jl") include("controlhost.jl") +include("oscillations.jl")