Merge pull request #13 from legend-exp/patch_spms-geds-calibration

Updated `SiPM` calibration and `LAr` cut definition

Author: theHenks

src/calibrate_all.jl

@@ -43,38 +43,37 @@ function calibrate_all(data::LegendData, sel::AnyValiditySelection, datastore::A
 
     trig_e_ch = findall.(ged_events_pre.is_physical_trig)
 
-    trig_e_trap_cal = _fix_vov(getindex.(ged_events_pre.e_trap_cal, trig_e_ch))
-    trig_e_cusp_cal = _fix_vov(getindex.(ged_events_pre.e_cusp_cal, trig_e_ch))
+    trig_e_trap_max_cal = _fix_vov(getindex.(ged_events_pre.e_trap_max_cal, trig_e_ch))
+    trig_e_cusp_max_cal = _fix_vov(getindex.(ged_events_pre.e_cusp_max_cal, trig_e_ch))
     trig_e_trap_ctc_cal = _fix_vov(getindex.(ged_events_pre.e_trap_ctc_cal, trig_e_ch))
     trig_e_cusp_ctc_cal = _fix_vov(getindex.(ged_events_pre.e_cusp_ctc_cal, trig_e_ch))
-    trig_e_short_cal = _fix_vov(getindex.(ged_events_pre.e_313_cal, trig_e_ch))
+    trig_e_535_cal      = _fix_vov(getindex.(ged_events_pre.e_535_cal, trig_e_ch))
     trig_t0 = _fix_vov(getindex.(ged_events_pre.t0, trig_e_ch))
     n_trig = length.(trig_e_ch)
     n_expected_baseline = length.(ged_events_pre.is_baseline) .- length.(trig_e_ch)
 
     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_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_cal = maximum_with_init.(trig_e_trap_max_cal),
+        max_e_cusp_cal = maximum_with_init.(trig_e_cusp_max_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),
+        max_e_short_cal = maximum_with_init.(trig_e_535_cal),
         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_max_cal = trig_e_trap_max_cal,
+        trig_e_cusp_max_cal = trig_e_cusp_max_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,
+        trig_e_535_cal = trig_e_535_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),
@@ -122,7 +121,7 @@ function calibrate_all(data::LegendData, sel::AnyValiditySelection, datastore::A
     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),
+        ged_spm = _build_lar_cut(data, sel, global_events),
     )
 
     result = StructArray(merge(columns(global_events), cross_systems_cols))

src/calibrate_geds.jl

@@ -11,6 +11,7 @@ Apply the calibration specified by `data` and `sel` for the given HPGe
 Also calculates the configured cut/flag values.
 """
 function calibrate_ged_channel_data(data::LegendData, sel::AnyValiditySelection, detector::DetectorIdLike, channel_data::AbstractVector; 
+        e_cal_pars_type::Symbol=:rpars, e_cal_pars_cat::Symbol=:ecal,
         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)
 
@@ -20,7 +21,7 @@ function calibrate_ged_channel_data(data::LegendData, sel::AnyValiditySelection,
     chdata = channel_data[:]
 
     # get energy and psd calibration functions for the detector
-    cal_pf = get_ged_cal_propfunc(data, sel, detector)
+    cal_pf = get_ged_cal_propfunc(data, sel, detector; pars_type=e_cal_pars_type, pars_cat=e_cal_pars_cat)
     psd_pf = get_ged_psd_propfunc(data, sel, detector; pars_type=psd_cal_pars_type, pars_cat=psd_cal_pars_cat)
 
     # get qc labels

src/calibrate_smps.jl

@@ -13,6 +13,8 @@ function calibrate_spm_channel_data(data::LegendData, sel::AnyValiditySelection,
     chdata = channel_data[:]
 
     spmcal_pf = get_spm_cal_propfunc(data, sel, detector)
+    spmdc_sel_pf = get_spm_dc_sel_propfunc(data, sel, detector)
+    spmdc_cal_pf = get_spm_dc_cal_propfunc(data, sel, detector)
 
     # get additional cols to be parsed into the event tier
     chdata_output_pf = if keep_chdata
@@ -24,9 +26,12 @@ function calibrate_spm_channel_data(data::LegendData, sel::AnyValiditySelection,
     cal_output_novv = spmcal_pf.(chdata)
     cal_output = StructArray(map(VectorOfArrays, columns(cal_output_novv)))
 
+    dc_output_novv = NamedTuple{keys(spmdc_sel_pf)}([spmdc_cal_pf[e_type].(spmdc_sel_pf[e_type].(chdata)) for e_type in keys(spmdc_sel_pf)])
+    dc_output = StructArray(map(VectorOfArrays, columns(dc_output_novv)))
+
     chdata_output = chdata_output_pf.(chdata)
 
-    return StructVector(merge(columns(cal_output), columns(chdata_output)))
+    return StructVector(merge(columns(cal_output), columns(dc_output), columns(chdata_output)))
 end
 export calibrate_spm_channel_data
 
@@ -45,34 +50,55 @@ function _single_fiber_esum(
 end
 
 
-# ToDo: Make cut criteria configurable:
 function _lar_cut(
+    colnames::Tuple{Symbol, Symbol, Symbol},
     t_win::AbstractInterval, spm_t::AbstractVector{<:AbstractVector{<:Number}},
-    spmdc::AbstractVector{<:AbstractVector{<:Number}}, spm_pe::AbstractVector{<:AbstractVector{<:Number}}
+    spmdc::AbstractVector{<:AbstractVector{<:Number}}, spm_pe::AbstractVector{<:AbstractVector{<:Number}},
+    pe_ch_threshold::Quantity{<:Real}, pe_sum_threshold::Quantity{<:Real}, multiplicity_threshold::Int
 )
     n_over_thresh::Int = 0
     pe_sum::eltype(eltype(spm_pe)) = zero(eltype(eltype(spm_pe)))
     for i in eachindex(spm_t)
         s_i = _single_fiber_esum(t_win, spm_t[i], spmdc[i], spm_pe[i])
-        if s_i > 0.5
+        if s_i > pe_ch_threshold
             n_over_thresh += 1
         end
         pe_sum += s_i
     end
 
-    lar_cut = n_over_thresh >= 4 || pe_sum >= 4
+    lar_cut = Bool(n_over_thresh >= multiplicity_threshold || pe_sum >= pe_sum_threshold)
 
-    return (lar_cut = lar_cut, spms_win_pe_sum = pe_sum, spms_win_multiplicity = n_over_thresh)
+    return NamedTuple{colnames, Tuple{Bool, eltype(eltype(spm_pe)), Int}}([lar_cut, pe_sum, n_over_thresh])
 end
 
-
-# ToDo: Make time window configurable:
-function _build_lar_cut(global_events::AbstractVector{<:NamedTuple})
+function _build_lar_cut(data::LegendData, sel::AnyValiditySelection, global_events::AbstractVector{<:NamedTuple}, e_filter::Symbol)
     geds_t0 = global_events.geds.t0_start
-    spm_t = global_events.spms.trig_pos
-    spmdc = global_events.spms.trig_is_dc
-    spm_pe = global_events.spms.trig_pe
+
+    dataprod_larcut = get_spms_evt_lar_cut_props(data, sel)
+    dataprod_larcut_filter = dataprod_larcut.energy_types[e_filter]
+
+    spm_t =  getproperty(global_events.spms, Symbol(dataprod_larcut_filter.pos))
+    spmdc =  getproperty(global_events.spms, Symbol(dataprod_larcut_filter.is_dc))
+    spm_pe = getproperty(global_events.spms, e_filter)
+
+    ged_sum_window = dataprod_larcut.ged_sum_window
+
+    t_wins = ClosedInterval.(geds_t0 .+ first(ged_sum_window), geds_t0 .+ last(ged_sum_window))
+
+    pe_ch_threshold = dataprod_larcut.pe_ch_threshold
+    pe_sum_threshold = dataprod_larcut.pe_sum_threshold
+    multiplicity_threshold = dataprod_larcut.multiplicity_threshold
+
+    colnames = Tuple(Symbol.("$(e_filter)_" .* ["lar_cut", "spms_win_pe_sum", "spms_win_multiplicity"]))
+    return StructArray(_lar_cut.(Ref(colnames), t_wins, spm_t, spmdc, spm_pe, Ref(pe_ch_threshold), Ref(pe_sum_threshold), Ref(multiplicity_threshold)))
+end
+
+function _build_lar_cut(data::LegendData, sel::AnyValiditySelection, global_events::AbstractVector{<:NamedTuple})
+    dataprod_larcut = get_spms_evt_lar_cut_props(data, sel)
+    energy_types = keys(dataprod_larcut.energy_types)
+
+    is_valid_lar_propfunc = ljl_propfunc(dataprod_larcut.is_valid_lar)
 
-    t_wins = @. ClosedInterval(geds_t0 - 1u"μs", geds_t0 + 5u"μs")
-    return StructArray(_lar_cut.(t_wins, spm_t, spmdc, spm_pe))
+    lar_cut = StructArray(merge(columns.(_build_lar_cut.(Ref(data), Ref(sel), Ref(global_events), energy_types))...))
+    return StructArray(merge((is_valid_lar = is_valid_lar_propfunc.(lar_cut),), columns(lar_cut)))
 end