Merge pull request #124 from DaGeibl/main

updated LQ sideband selection logic

Author: theHenks

ext/LegendSpecFitsRecipesBaseExt.jl

@@ -1390,9 +1390,9 @@ end
     xlabel := "Drift Time"
     ylabel := "LQ (A.U.)"
     framestyle := :box
-    left_margin := -2Plots.mm
-    bottom_margin := -4Plots.mm
-    top_margin := -3Plots.mm
+    left_margin := (-2, :mm)
+    bottom_margin := (-4, :mm)
+    top_margin := (-3, :mm)
     color := :viridis
     formatter := :plain
     thickness_scaling := 1.6
@@ -1476,15 +1476,15 @@ end
 
 # recipe for the lq_cut report
 
-@recipe function f(report::NamedTuple{(:cut, :fit_result, :temp_hists, :fit_report)}, lq_class::Vector{Float64}, e_cal, plot_type::Symbol)
+@recipe function f(report::NamedTuple{(:cut, :fit_result, :temp_hists, :fit_report, :dep_σ, :edges)}, lq_class::Vector{Float64}, e_cal, plot_type::Symbol)
 
     # Extract cutvalue from the report
     cut_value = Measurements.value.(report.cut)
 
     # Plot configuration for all types
-    left_margin := -2Plots.mm
-    bottom_margin := -4Plots.mm
-    top_margin := -3Plots.mm
+    left_margin := (-2, :mm)
+    bottom_margin := (-4, :mm)
+    top_margin := (-3, :mm)
     thickness_scaling := 1.6
     size := (1200, 900)
     framestyle := :box
@@ -1662,10 +1662,40 @@ end
             report.temp_hists.hist_sb2
         end
 
-
+    elseif plot_type == :energy_spectrum
+        # Plot energy spectrum with DEP and sideband regions
+        left_margin := (-2, :mm)
+        bottom_margin := (-4, :mm)
+        top_margin := (-3, :mm)
+        thickness_scaling := 1.6
+        size := (1200, 900)
+        framestyle := :box
+        formatter := :plain
+        xlabel := "Energy"
+        ylabel := "Counts"
+        
+        @series begin
+            seriestype := :stephist
+            label := "Energy Spectrum (σ: $(round(u"keV", report.dep_σ, sigdigits=3)))"
+            bins := 1000
+            e_cal[1500u"keV" .< e_cal .< 1660u"keV"]
+        end
+    
+        @series begin
+            seriestype := :vline
+            label := "DEP Region"
+            fillcolor := :red
+            [report.edges.DEP_edge_left, report.edges.DEP_edge_right]
+        end
+    
+        @series begin
+            seriestype := :vline
+            label := "Sideband Edges"
+            fillcolor := :green
+            legend := :topleft
+            [report.edges.sb1_edge, report.edges.sb2_edge]
+        end
     end
 end
 
-
-
 end # module LegendSpecFitsRecipesBaseExt

src/lqcut.jl

@@ -159,11 +159,29 @@ function lq_cut(
     dep_µ::Unitful.Energy, dep_σ::Unitful.Energy, e_cal::Vector{<:Unitful.Energy}, lq_classifier::Vector{<:AbstractFloat}; cut_sigma::Float64=3.0, dep_sideband_sigma::Float64=4.5, cut_truncation_sigma::Float64=3.5, uncertainty::Bool=true
     )
 
-    # define sidebands
-    lq_dep = lq_classifier[dep_µ - dep_sideband_sigma * dep_σ .< e_cal .< dep_µ + dep_sideband_sigma * dep_σ]
-    lq_sb1 = lq_classifier[dep_µ -  2 * dep_sideband_sigma * dep_σ .< e_cal .< dep_µ - dep_sideband_sigma * dep_σ]
-    lq_sb2 = lq_classifier[dep_µ + dep_sideband_sigma * dep_σ .< e_cal .< dep_µ + 2 * dep_sideband_sigma * dep_σ]
-    
+    # define sidebands; different for low and high energy resolution detectors to avoid sb reaching into 212-Bi FEP
+    DEP_edge_left  = dep_µ - dep_sideband_sigma * dep_σ
+    DEP_edge_right = dep_µ + dep_sideband_sigma * dep_σ
+
+    if dep_σ < 2.0u"keV"
+        sb1_edge = dep_µ - 2 * dep_sideband_sigma * dep_σ
+        sb2_edge = dep_µ + 2 * dep_sideband_sigma * dep_σ  
+
+        lq_dep = lq_classifier[DEP_edge_left .< e_cal .< DEP_edge_right]
+        lq_sb1 = lq_classifier[sb1_edge .< e_cal .< DEP_edge_left]
+        lq_sb2 = lq_classifier[DEP_edge_right .< e_cal .< sb2_edge]
+    else
+        sb1_edge = dep_µ - 2 * dep_sideband_sigma * dep_σ  
+        sb2_edge = dep_µ - 3 * dep_sideband_sigma * dep_σ
+
+        lq_dep = lq_classifier[DEP_edge_left .< e_cal .< DEP_edge_right]
+        lq_sb1 = lq_classifier[sb1_edge .< e_cal .< DEP_edge_left]
+        lq_sb2 = lq_classifier[sb2_edge .< e_cal .< sb1_edge]
+    end
+
+    # save edges for crosschecks
+    edges_for_crosschecks = (;DEP_edge_left, DEP_edge_right, sb1_edge, sb2_edge)
+
     # generate values for histogram edges
     combined = filter(isfinite,[lq_dep; lq_sb1; lq_sb2])
     ideal_bin_width = get_friedman_diaconis_bin_width(combined)
@@ -206,6 +224,8 @@ function lq_cut(
         fit_result = fit_result,
         temp_hists = temp_hists,
         fit_report = fit_report,
+        dep_σ = dep_σ,
+        edges = edges_for_crosschecks,
     )
 
     return result, report