Repository

Example of reading KM3NeT Open Data

This notebook demonstrates read and fill histograms from the KM3NeT/ORCA 433 kt-y oscillations data release using KM3NeT's tools.

Setup and Dependencies

First, we load the necessary packages:

  • KM3io: For reading KM3NeT data formats

  • NuFlux: For neutrino flux calculations

  • KM3OpenOsc: For neutrino oscillation calculations

begin
    using KM3io
    using NuFlux
    using KM3OpenOsc
end

Downloading Open Dataset

ORCA 433 kt-y open dataset can be downloaded from the following dataverse link and later unzipped to have access to the files needed. After unzipping the file into a directory, 5 files will be found including a README among them where the content of each file is carefully explained.

Loading Data Files into the notebook

Here we define paths to two key files:

  • PATH_UNZIPPED_DATASET: Path to unzipped open dataset

  • BINDEF: JSON file containing binning definitions for the analysis

  • OSCFILE: ROOT file containing data from KM3NeT/ORCA

begin
    PATH_UNZIPPED_DATASET = "../../unzipped_ORCA433kt-y_open_dataset/" # to change as a function of the user
    BINDEF = PATH_UNZIPPED_DATASET * "bins_433kt-y_v0.5.json"
    OSCFILE = PATH_UNZIPPED_DATASET * "ORCA6_433kt-y_opendata_v0.5.root"
end
"../../unzipped_ORCA433kt-y_open_dataset/ORCA6_433kt-y_opendata_v0.5.root"

Reading the Oscillation Data File

To read the file we use KM3io.OSCFile which is a function that directly identifies the type of file and what is inside

We load the oscillation data file and extract three data samples:

  • nu: Neutrino Monte Carlo events

  • data: Experimental data events

  • muons: Atmospheric muon background events

f = KM3io.OSCFile(OSCFILE)
OSCFile{OscOpenDataTree of Neutrinos (592099 events), OscOpenDataTree of Data (900 events), OscOpenDataTree of Muons (900 events)}
begin
    nu = f.osc_opendata_nu
    data = f.osc_opendata_data
    muons = f.osc_opendata_muons
end
OscOpenDataTree of Muons (900 events)

Creating Histogram Containers

We initialize empty 2D histogram containers using the binning definition file:

  • hn: For simulated neutrino events

  • hd: For real experimental data

  • hm: For simulated atmospheric muon events

begin
    hn = create_histograms(BINDEF)
    hd = create_histograms(BINDEF)
    hm = create_histograms(BINDEF); nothing
end
fieldnames(typeof(hn))
(:hists_true, :hists_reco)
keys(hn.hists_reco)
KeySet for a Dict{String, Hist2D} with 4 entries. Keys:
  "reco"
  "recoShowers"
  "recoHighPurityTracks"
  "recoLowPurityTracks"
hn.hists_reco["recoShowers"]
2.01000.0-1.01.0
  • edges: ([2.0, 2.95752, 4.37345, 5.3183, 6.46727, 7.86447, 9.56352, 11.6296, 14.1421, 17.1974, 20.9128, 25.4308, 30.9249, 55.6102, 100.0, 1000.0], [-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.3999999999999999, -0.29999999999999993, -0.19999999999999996, -0.09999999999999998 … 0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6000000000000001, 0.7000000000000001, 0.8, 0.9, 1.0])
  • bin counts: [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  • maximum count: 0.0
  • total count: 0.0

Setting Oscillation Parameters

We define the oscillation parameters using the best fit values from the standard oscillations analysis from ORCA 433 kt-y paper JHEP10(2024)206:

  • dm_21, dm_31: Mass-squared differences in eV²

  • theta_12, theta_23, theta_13: Mixing angles

  • dcp: CP-violating phase

These parameters are used to calculate the oscillation matrices (U and H). It uses Neurthino.jl behind the scenes.

begin
    BF = Dict("dm_21" => 7.42e-5,
                       "dm_31" => 2.18e-3,
                       "theta_12" => deg2rad(33.45),
                       "theta_23" => deg2rad(45.57299599919429),
                       "theta_13" => deg2rad(8.62),
                       "dcp" => deg2rad(230))

    U,H = get_oscillation_matrices(BF)
end
(ComplexF64[0.8249422516589652 + 0.0im 0.5449826827686066 + 0.0im -0.09634131252970905 + 0.1148151053225223im; -0.3284406729882627 + 0.06841342277740595im 0.6219810005635141 + 0.04519604930846363im 0.7060759732322109 + 0.0im; 0.449910833254889 + 0.06705856737677285im -0.5586843920183072 + 0.04430098940637267im 0.6920928861894314 + 0.0im], ComplexF64[-0.0007514 + 0.0im, -0.0006772000000000001 + 0.0im, 0.0014286000000000001 + 0.0im])

Loading Neutrino Flux Model

We load the Honda Frejus site flux model, which provides the atmospheric neutrino flux predictions needed for the analysis. Any Honda formatted table can be read using this function. It uses NuFlux.jl behind the scenes.

begin
    NUFLUX_PATH = split(Base.pathof(NuFlux), "src")[1]
    FLUX_DATA_DIR = joinpath(NUFLUX_PATH, "data")
    flux_path = joinpath(FLUX_DATA_DIR, "frj-nu-20-01-000.d") # Get flux of Honda Frejus site from NuFlux
    flux_dict = get_flux_dict(flux_path)
end
Dict{Int32, NuFlux.FluxTable} with 4 entries:
  -12 => FluxTable([-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1  …  0.1,…
  12  => FluxTable([-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1  …  0.1,…
  14  => FluxTable([-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1  …  0.1,…
  -14 => FluxTable([-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1  …  0.1,…

Filling Histograms with Events

We populate our histogram containers with events:

  • Neutrinos with flux weighting and oscillation effects applied

  • Atmospheric muons

  • Experimental data

There is an optional argument which is livetime that can be used to scale up or down montecarlo

begin
    fill_response!(hn, nu, flux_dict, U, H; oscillations=true) # fill neutrinos ,need flux, oscillation parameters
    fill_response!(hm, muons)
    fill_response!(hd, data)
end
hn.hists_true["muon_cc_nu"]
1.010000.0-1.01.0
  • edges: ([1.0, 1.2589254117941673, 1.5848931924611136, 1.9952623149688797, 2.51188643150958, 3.1622776601683795, 3.9810717055349727, 5.011872336272723, 6.309573444801933, 7.943282347242818 … 1258.925411794168, 1584.8931924611143, 1995.2623149688823, 2511.886431509581, 3162.2776601683804, 3981.071705534977, 5011.872336272724, 6309.573444801938, 7943.28234724283, 10000.0], [-1.0, -0.975, -0.95, -0.925, -0.9, -0.875, -0.85, -0.825, -0.8, -0.775 … 0.775, 0.8, 0.8250000000000001, 0.8500000000000001, 0.875, 0.9, 0.925, 0.9500000000000001, 0.9750000000000001, 1.0])
  • bin counts: [0.04085472213947171 0.004977305716097258 … 0.0 0.0; 0.12829694773170977 0.008506359166404558 … 0.0 0.0; … ; 0.05740980482210812 0.057535779183707024 … 0.00011768083104523938 0.0; 0.037587732908178864 0.0380479135922043 … 0.0 0.0]
  • maximum count: 6.423147039212604
  • total count: 2915.5498586995577
hn.hists_reco["recoShowers"]
2.01000.0-1.01.0
  • edges: ([2.0, 2.95752, 4.37345, 5.3183, 6.46727, 7.86447, 9.56352, 11.6296, 14.1421, 17.1974, 20.9128, 25.4308, 30.9249, 55.6102, 100.0, 1000.0], [-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.3999999999999999, -0.29999999999999993, -0.19999999999999996, -0.09999999999999998 … 0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6000000000000001, 0.7000000000000001, 0.8, 0.9, 1.0])
  • bin counts: [3.191829228181372 4.369504546333975 … 0.0 0.0; 11.66977246088045 13.448746078325804 … 0.0 0.0; … ; 19.297699494383636 15.513395665284698 … 0.0 0.0; 12.833286366966924 12.16150406073552 … 0.0 0.0]
  • maximum count: 28.520710756581217
  • total count: 2032.1392944893719

Exporting results to HDF5

You can easily export the filled histograms to HDF5.

Additionally it is possible to export the full response to a single hdf5 file

begin
    export_histograms_hdf5(hn, "neutrino_histograms_from_testdata.h5") # You can easily export the filled histograms to hdf5
    h5f = build_HDF5_file("responses_to_file.h5") # Create h5 file with same structure as responses bins 
    fill_HDF5_file!(h5f, nu, "neutrinos") # Completely export the response as a table in an hdf5 file at a given path 
    fill_HDF5_file!(h5f, muons, "muons")
    fill_HDF5_file!(h5f, data, "data") 
    close(h5f)
end

Visualizing the Reconstruction Channels

First, we load the necessary packages:

  • CairoMakie: For plotting

  • FHist: For histogram operations

This plot shows three different reconstruction channels:

  1. High Purity Tracks: Events classified with track topology with low atmospheric muon contamination

  2. Low Purity Tracks: Events classified with track topology with higher atmospheric muon contamination

  3. Showers: Events classified with a shower topology

For each channel, we display:

  • Stacked histograms of neutrino and muon Monte Carlo events

  • Data points with error bars

The x-axis represents the cosine of the zenith angle (cos θ), where -1 means upgoing events (through the Earth) and 0 means horizontal events.

using CairoMakie, FHist
with_theme(ATLASTHEME) do
    reco_types = ["recoHighPurityTracks", "recoLowPurityTracks", "recoShowers"]
     f = Figure(size = (1200, 400))

    for (i, reco_type) in enumerate(reco_types)
        neutrinos = project(hn.hists_reco[reco_type],:y)
        muons = project(hm.hists_reco[reco_type],:y)
        data = project(hd.hists_reco[reco_type],:y)
        ax = Axis(f[1, i], 
                 xlabel = "Cos theta", 
                 ylabel = i == 1 ? "Event count" : "", 
                 title = reco_type)
        p = stackedhist!(ax, [neutrinos, muons]; error_color=Pattern('/'))
        xlims!(ax, (-1, 0))
        scatter!(ax, bincenters(data), bincounts(data), 
                 color = :black, 
                 markersize = 10)
        if i == 3
            labels = ["Neutrinos", "Muons", "Data"]
            elements = [
                [PolyElement(polycolor = p.attributes.color[][j]) for j in 1:2]..., 
                MarkerElement(color = :black, marker = :circle, markersize = 8)
            ]
        end
    end
    f
end

Examining True Event Categories

This plot shows the distribution of true Monte Carlo event categories:

  1. muon_cc_nu: Charged-current muon neutrino events

  2. muon_cc_nub: Charged-current muon antineutrino events

with_theme(ATLASTHEME) do
    reco_types = ["muon_cc_nu", "muon_cc_nub"]
     f = Figure(size = (1200, 400))
    # Loop through each reconstruction type
    for (i, reco_type) in enumerate(reco_types)
        muon_cc_nu = project(hn.hists_true[reco_type],:y)
        muon_cc_nub = project(hm.hists_true[reco_type],:y)
        ax = Axis(f[1, i], 
                 xlabel = "Cos theta", 
                 ylabel = i == 1 ? "Event count" : "", 
                 title = reco_type)
        p = stackedhist!(ax, [muon_cc_nu, muon_cc_nub]; error_color=Pattern('/'))
        xlims!(ax, (-1, 0))
        if i == 3
            labels = ["muon_cc_nu", "muon_cc_nub"]
            elements = [
                [PolyElement(polycolor = p.attributes.color[][j]) for j in 1:2]..., 
                MarkerElement(color = :black, marker = :circle, markersize = 8)
            ]
        end
    end
    f
end

Built with Julia 1.11.2 and

CairoMakie 0.13.10
FHist 0.11.11
KM3OpenOsc 0.1.9
KM3io 0.18.7
NuFlux 0.1.6