Merge pull request #8 from legend-exp/dev

Bug Fixes

Author: oschulz

.github/workflows/ci.yml

@@ -25,7 +25,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9'
+          - '1.10'
           - '1'
           - 'nightly'
         os:

Project.toml

@@ -23,4 +23,4 @@ PropertyFunctions = "0.2.1"
 StructArrays = "0.4, 0.5, 0.6"
 Tables = "1"
 Unitful = "1.6"
-julia = "1.9"
+julia = "1.10"

src/LegendEventAnalysis.jl

@@ -12,7 +12,7 @@ using LegendDataManagement
 
 using IntervalSets: AbstractInterval, ClosedInterval, leftendpoint, rightendpoint
 using ProgressMeter: @showprogress
-using PropertyFunctions: PropertyFunction, @pf, filterby
+using PropertyFunctions: PropertyFunction, @pf, filterby, PropSelFunction
 
 using Tables: columns
 using Unitful: RealOrRealQuantity as RealQuantity
@@ -23,7 +23,7 @@ include("flatten_over_channels.jl")
 include("build_global_events.jl")
 include("calibrate_geds.jl")
 include("calibrate_smps.jl")
-include("calibrate_puls.jl")
+include("calibrate_aux.jl")
 include("calibrate_all.jl")
 
 end # module

src/calibrate_all.jl

@@ -7,25 +7,32 @@
 Calibrate all channels in the given datastore, using the metadata
 processing configuration for `data` and `sel`.
 """
-function calibrate_all(data::LegendData, sel::AnyValiditySelection, datastore::AbstractDict; tier::DataTierLike=:jldsp)
+function calibrate_all(data::LegendData, sel::AnyValiditySelection, datastore::AbstractDict, tier::DataTierLike=:jldsp)
     ds = datastore
 
+    @debug "Calibrating all channels in for `ValiditySelection` $(sel) in `DataTier` $(tier)"
     chinfo = channelinfo(data, sel)
-    geds_channels::Vector{ChannelId} = filterby(@pf $system == :geds && $processable && $usability != :off && $is_blinded)(chinfo).channel
-    spms_channels::Vector{ChannelId} = filterby(@pf $system == :spms && $processable)(chinfo).channel
-    puls_channels::Vector{ChannelId} = filterby(@pf $system in [:puls, :bsln])(chinfo).channel
-
+    geds_channels::Vector{ChannelId} = filterby(get_ged_evt_chsel_propfunc(data, sel))(chinfo).channel
+    @debug "Loaded $(length(geds_channels)) HPGe channels"
+    hitgeds_channels::Vector{ChannelId} = filterby(get_ged_evt_hitchsel_propfunc(data, sel))(chinfo).channel
+    @debug "Loaded $(length(hitgeds_channels)) HPGe hit channels"
+    spms_channels::Vector{ChannelId} = filterby(get_spms_evt_chsel_propfunc(data, sel))(chinfo).channel
+    @debug "Loaded $(length(spms_channels)) SiPM channels"
+    aux_channels::Vector{ChannelId} = filterby(get_aux_evt_chsel_propfunc(data, sel))(chinfo).channel
+    @debug "Loaded auxiliary channels: $(join(string.(filterby(get_aux_evt_chsel_propfunc(data, sel))(chinfo).detector), ", "))"
 
     # HPGe:
-
-    ged_caldata = Dict([
-        let detector = channelinfo(data, sel, channel).detector,
-            chdata = ds[channel, tier][:]
-            channel => calibrate_ged_channel_data(data, sel, detector, chdata)
+    @debug "Calibrating HPGe channels"
+    ged_kwargs = get_ged_evt_kwargs(data, sel)
+    ged_caldata_v = Vector{StructVector}(undef, length(geds_channels))
+    Threads.@threads for i in eachindex(geds_channels)
+        let detector = channelinfo(data, sel, geds_channels[i]).detector, chdata = ds[geds_channels[i], tier][:]
+            ged_caldata_v[i] = calibrate_ged_channel_data(data, sel, detector, chdata; ged_kwargs...)
         end
-        for channel in geds_channels
-    ])
+    end
+    ged_caldata = Dict(geds_channels .=> ged_caldata_v)
 
+    @debug "Building global events for HPGe channels"
     ged_events_pre = build_global_events(ged_caldata, geds_channels)
 
     min_t0(t0::AbstractVector{<:Number}) = isempty(t0) ? eltype(t0)(NaN) : minimum(t0)
@@ -47,23 +54,29 @@ function calibrate_all(data::LegendData, sel::AnyValiditySelection, datastore::A
 
     maximum_with_init(A) = maximum(A, init=zero(eltype((A))))
 
+    is_valid_hit(trig_chs::AbstractVector{<:Int}, hit_channels::AbstractVector{<:Int}) = all(x -> x in hit_channels, trig_chs)
+    # Main.@infiltrate
+
     ged_additional_cols = (
         t0_start = min_t0.(trig_t0),
         trig_t0 = trig_t0,
         multiplicity = n_trig,
-        max_e_ch = max_e_ch,
+        max_e_ch_idxs = max_e_ch,
+        max_e_ch = only.(getindex.(ged_events_pre.channel, max_e_ch)),
         max_e_trap_cal = maximum_with_init.(trig_e_trap_cal),
         max_e_cusp_cal = maximum_with_init.(trig_e_cusp_cal),
         max_e_trap_ctc_cal = maximum_with_init.(trig_e_trap_ctc_cal),
         max_e_cusp_ctc_cal = maximum_with_init.(trig_e_cusp_ctc_cal),
         max_e_short_cal = maximum_with_init.(trig_e_short_cal),
-        trig_e_ch = trig_e_ch,
+        trig_e_ch_idxs = trig_e_ch,
+        trig_e_ch = getindex.(ged_events_pre.channel, trig_e_ch),
         trig_e_trap_cal = trig_e_trap_cal,
         trig_e_cusp_cal = trig_e_cusp_cal,
         trig_e_trap_ctc_cal = trig_e_trap_ctc_cal,
         trig_e_cusp_ctc_cal = trig_e_cusp_ctc_cal,
         trig_e_short_cal = trig_e_short_cal,
         is_valid_qc = count.(ged_events_pre.is_baseline) .== n_expected_baseline,
+        is_valid_hit = is_valid_hit.(getindex.(ged_events_pre.channel, trig_e_ch), Ref(Int.(hitgeds_channels))),
         is_discharge_recovery = any.(ged_events_pre.is_discharge_recovery_ml),
         is_saturated = any.(ged_events_pre.is_saturated),
         is_discharge = any.(ged_events_pre.is_discharge),
@@ -72,41 +85,41 @@ function calibrate_all(data::LegendData, sel::AnyValiditySelection, datastore::A
 
 
     # SiPM:
-
-    spm_caldata = Dict([
-        let detector = channelinfo(data, sel, channel).detector,
-            chdata = ds[channel, tier][:]
-            channel => calibrate_spm_channel_data(data, sel, detector, chdata)
+    @debug "Calibrating SiPM channels"
+    spm_caldata_v = Vector{StructVector}(undef, length(spms_channels))
+    Threads.@threads for i in eachindex(spms_channels)
+        let detector = channelinfo(data, sel, spms_channels[i]).detector, chdata = ds[spms_channels[i], tier][:]
+            spm_caldata_v[i] = calibrate_spm_channel_data(data, sel, detector, chdata)
         end
-        for channel in spms_channels
-    ])
+    end
+    spm_caldata = Dict(spms_channels .=> spm_caldata_v)
+    @debug "Building global events for SiPM channels"
     spm_events_novov = build_global_events(spm_caldata, spms_channels)
     spm_events = StructArray(map(_fix_vov, columns(spm_events_novov)))
 
-
-    # Pulser:
-
-    pls_caldata = Dict([
-        let detector = channelinfo(data, sel, channel).detector,
-            chdata = ds[channel, tier][:]
-            channel => calibrate_pls_channel_data(data, sel, detector, chdata)
-        end
-        for channel in puls_channels if haskey(ds, string(channel))
-    ])
-    pls_events = build_global_events(pls_caldata, puls_channels)
 
+    @debug "Calibrating auxiliary channels"
+    # aux & Forced Trigger
+    aux_caldata = 
+        [Dict(
+            let detector = channelinfo(data, sel, channel).detector,
+                chdata = ds[channel, tier][:]
+                channel => calibrate_aux_channel_data(data, sel, detector, chdata)
+            end
+            ) for channel in aux_channels]
+    @debug "Building global events for auxiliary channels"
+    aux_events = NamedTuple{Tuple(get_aux_evt_levelname_propfunc.(Ref(data), Ref(sel), reduce(vcat, collect.(keys.(aux_caldata)))))}(build_global_events.(aux_caldata))
 
     # Cross-system:
-
-    system_events = (
+    @debug "Building cross-system events"
+    system_events = merge((
         geds = ged_events,
         spms = spm_events,
-        puls = pls_events,
-    )
+    ), aux_events)
 
     global_events_pre = build_cross_system_events(system_events)
-    single_pls_col = StructVector(map(Broadcast.BroadcastFunction(only), columns(global_events_pre.puls)))
-    global_events = StructVector(merge(columns(global_events_pre), (puls = single_pls_col,)))
+    aux_cols = NamedTuple{keys(aux_events)}([StructVector(map(Broadcast.BroadcastFunction(only), columns(getproperty(global_events_pre, k)))) for k in keys(aux_events)])
+    global_events = StructVector(merge(Base.structdiff(columns(global_events_pre), NamedTuple{keys(aux_events)}), (aux = StructArray(aux_cols),)))
 
     cross_systems_cols = (
         ged_spm = _build_lar_cut(global_events),

src/calibrate_puls.jl src/calibrate_aux.jlsrc/calibrate_puls.jl renamed to src/calibrate_aux.jl

@@ -1,28 +1,28 @@
 # This file is a part of LegendEventAnalysis.jl, licensed under the MIT License (MIT).
 
 """
-    calibrate_pls_channel_data(data::LegendData, sel::ValiditySelection, detector::DetectorId, channel_data::AbstractVector)
+    calibrate_aux_channel_data(data::LegendData, sel::ValiditySelection, detector::DetectorId, channel_data::AbstractVector)
 
-Apply the calibration specified by `data` and `sel` for pulser referred to
+Apply the calibration specified by `data` and `sel` for aux channel referred to
 by the `detector` ID to the single-channel `channel_data` for that detector.
 
 Also calculates the configured cut/flag values.
 """
-function calibrate_pls_channel_data(data::LegendData, sel::AnyValiditySelection, detector::DetectorId, channel_data::AbstractVector)
+function calibrate_aux_channel_data(data::LegendData, sel::AnyValiditySelection, detector::DetectorId, channel_data::AbstractVector)
     chdata = channel_data[:]
 
     # get calibration function for the detector
-    pulsercal_pf = get_pulser_cal_propfunc(data, sel, detector)
+    auxcal_pf = get_aux_cal_propfunc(data, sel, detector)
 
     # get additional cols to be parsed into the event tier
-    chdata_output_pf = get_ged_evt_chdata_propfunc(data, sel, detector)
+    chdata_output_pf = get_aux_evt_chdata_propfunc(data, sel, detector)
 
     # get additional columns
     chdata_output = chdata_output_pf.(chdata)
 
     # apply calibrations
-    cal_output = pulsercal_pf.(chdata)
+    cal_output = auxcal_pf.(chdata)
 
     return StructVector(merge(columns(chdata_output), columns(cal_output)))
 end
-export calibrate_pls_channel_data
+export calibrate_aux_channel_data

src/calibrate_geds.jl

@@ -1,18 +1,25 @@
 # This file is a part of LegendEventAnalysis.jl, licensed under the MIT License (MIT).
 
 """
-    calibrate_ged_channel_data(data::LegendData, sel::ValiditySelection, detector::DetectorId, channel_data::AbstractVector)
+    calibrate_ged_channel_data(data::LegendData, sel::AnyValiditySelection, detector::DetectorId, channel_data::AbstractVector; 
+        psd_cal_pars_type::Symbol=:ppars, psd_cal_pars_cat::Symbol=:aoe, psd_cut_pars_type::Symbol=:ppars, psd_cut_pars_cat::Symbol=:aoe,
+        keep_chdata::Bool=false)
 
 Apply the calibration specified by `data` and `sel` for the given HPGe
 `detector` to the single-channel `channel_data` for that detector.
 
 Also calculates the configured cut/flag values.
 """
-function calibrate_ged_channel_data(data::LegendData, sel::AnyValiditySelection, detector::DetectorId, channel_data::AbstractVector)
+function calibrate_ged_channel_data(data::LegendData, sel::AnyValiditySelection, detector::DetectorId, channel_data::AbstractVector; 
+        psd_cal_pars_type::Symbol=:ppars, psd_cal_pars_cat::Symbol=:aoe, psd_cut_pars_type::Symbol=:ppars, psd_cut_pars_cat::Symbol=:aoe,
+        keep_chdata::Bool=false)
+    
+    # get all chdata
     chdata = channel_data[:]
 
     # get energy and psd calibration functions for the detector
     cal_pf = get_ged_cal_propfunc(data, sel, detector)
+    psd_pf = get_ged_psd_propfunc(data, sel, detector; pars_type=psd_cal_pars_type, pars_cat=psd_cal_pars_cat)
 
     # get qc labels
     cut_pf = get_ged_qc_cuts_propfunc(data, sel, detector)
@@ -23,29 +30,29 @@ function calibrate_ged_channel_data(data::LegendData, sel::AnyValiditySelection,
     cut_is_trig_pf = get_ged_qc_is_trig_propfunc(data, sel, detector)
 
     # get additional cols to be parsed into the event tier
-    chdata_output_pf = get_ged_evt_chdata_propfunc(data, sel, detector)
+    chdata_output_pf = if keep_chdata
+        PropSelFunction{propertynames(chdata)}()
+    else
+        get_ged_evt_chdata_propfunc(data, sel)
+    end
 
     # get postcal psd flags
-    postcal_pf = let aoe_window = dataprod_pars_aoe_window(data, sel, detector)
-        @pf (
-            aeo_low_cut = $aoe_classifier > leftendpoint(aoe_window),
-            aoe_ds_cut = $aoe_classifier in aoe_window,
-        )
-    end
+    postcal_pf = get_ged_aoe_cut_propfunc(data, sel, detector; pars_type=psd_cut_pars_type, pars_cat=psd_cut_pars_cat)
 
     # apply calibrations 
     cal_output = cal_pf.(chdata)
-    cal_chdata = StructArray(merge(columns(cal_output), columns(chdata)))
+    psd_output = psd_pf.(StructArray(merge(columns(cal_output), columns(chdata))))
+    cal_chdata = StructArray(merge(columns(cal_output), columns(psd_output), columns(chdata)))
 
     # get cut labels
     cut_output = cut_pf.(cal_chdata)
 
-    # get additional columns
-    chdata_output = chdata_output_pf.(chdata)
-
     # get postcal data
     postcal_data = postcal_pf.(cal_chdata)
 
+    # get additional columns 
+    chdata_output = chdata_output_pf.(chdata)
+
     # get cut flags
     is_physical = cut_is_physical_pf.(cut_output)
     is_baseline = cut_is_baseline_pf.(cut_output)
@@ -57,6 +64,6 @@ function calibrate_ged_channel_data(data::LegendData, sel::AnyValiditySelection,
         is_physical_trig = is_physical_trig,
     )
 
-    return StructVector(merge(columns(chdata_output), columns(cal_output), columns(postcal_data), columns(cut_output), additional_cols))
+    return StructVector(merge(columns(chdata_output), columns(cal_output), columns(psd_output), columns(postcal_data), columns(cut_output), additional_cols))
 end
 export calibrate_ged_channel_data

src/calibrate_smps.jl

@@ -8,21 +8,25 @@ Apply the calibration specified by `data` and `sel` for the given SiPM
 
 Also calculates the configured cut/flag values.
 """
-function calibrate_spm_channel_data(data::LegendData, sel::AnyValiditySelection, detector::DetectorId, channel_data::AbstractVector)
+function calibrate_spm_channel_data(data::LegendData, sel::AnyValiditySelection, detector::DetectorId, channel_data::AbstractVector; 
+    keep_chdata::Bool=false)
     chdata = channel_data[:]
 
     spmcal_pf = get_spm_cal_propfunc(data, sel, detector)
 
-    # ToDo: Make channel output configurable:
-    chdata_output = (
-        timestamp = chdata.timestamp,
-        trig_pos = VectorOfArrays(chdata.trig_pos),
-    )
+    # get additional cols to be parsed into the event tier
+    chdata_output_pf = if keep_chdata
+        PropSelFunction{propertynames(chdata)}()
+    else
+        get_spms_evt_chdata_propfunc(data, sel)
+    end
 
     cal_output_novv = spmcal_pf.(chdata)
     cal_output = StructArray(map(VectorOfArrays, columns(cal_output_novv)))
 
-    return StructVector(merge(chdata_output, columns(cal_output)))
+    chdata_output = chdata_output_pf.(chdata)
+
+    return StructVector(merge(columns(cal_output), columns(chdata_output)))
 end
 export calibrate_spm_channel_data
 
@@ -58,7 +62,7 @@ function _lar_cut(
 
     lar_cut = n_over_thresh >= 4 || pe_sum >= 4
 
-    return (lar_cut = lar_cut, smps_win_pe_sum = pe_sum)
+    return (lar_cut = lar_cut, spms_win_pe_sum = pe_sum, spms_win_multiplicity = n_over_thresh)
 end