feat: opt-in detector-broadened continuum step at absorption edges

Sources/SpectroEDX/Fit/Continuum.swift

@@ -29,6 +29,22 @@
 //  the model has no nonlinear parameter: n_0 + 1 and n_1 + 1 coefficients, 10 + 6 = 16 by default, all
 //  linear and non-negative.
 //
+//  THE STEP SEEN THROUGH THE DETECTOR (opt-in `ContinuumEdgeStep.broadened`): a detector records the incoming spectrum
+//  convolved with its resolution, so a measured edge is a Gaussian-blurred step (about 75 eV FWHM at 1.56 keV for 128 eV at
+//  Mn K-alpha), never a jump between two channels. The broadened form weights segment s by
+//      w_s(E) = Phi((E - a_s) / sigma(a_s)) - Phi((E - b_s) / sigma(b_s)),
+//  Phi the standard normal CDF, sigma the line-width law's Gaussian sigma at the edge (`XRayLines.fwhm`, the line model's own
+//  resolution), no factor at the fitted range's ends; each segment's polynomial continues past its edge at its end value
+//  (t clamped to [0, 1]). The weights sum to 1 whenever every segment holds channels (a segment with none, two edges in one
+//  channel gap, is skipped and its weight is lost) and are non-negative, so the continuum stays non-negative. This is the step
+//  convolved with the detector Gaussian when the polynomials are flat across the edge's width. An edge outside the fitted range
+//  is not modelled, including the blurred tail of a step just outside it.
+//  Measured: on a pure Al wedge (Titan 200 kV, Super-X; Snelson, Collins, MacArthur, University of Leeds, doi 10.5518/1826,
+//  CC BY 4.0) the sharp step made the model fall 351 counts more than the data between the channels at 1.559 and 1.564 keV,
+//  left the data 1677 counts above the model in 1.53-1.65 keV and the model 39 % above the data in 1.65-2.2 keV; broadened:
+//  42 counts, 907 counts and 19 %. `.sharp` (the default) is the earlier form, value for value: the kit's registered
+//  synthetic generators plant an unbroadened step, so making `.broadened` the default re-measures those results first.
+//
 //  The Al K edge energy 1.5596 keV is the standard tabulated value (X-ray Data Booklet, Table 1-1).
 //  The Si K edge (1.839 keV, a detector dead-layer step) is NOT in the default: add it to `edges` to test it, or model the step
 //  physically with the opt-in dead-layer term (`FitPolicy.continuumDeadLayer`, `DeadLayerFit.swift`), which multiplies `efficiency`.
@@ -60,7 +76,7 @@ public nonisolated struct ContinuumForm: Sendable, Equatable {
 
     public static func == (a: ContinuumForm, b: ContinuumForm) -> Bool {
         a.beamEnergy == b.beamEnergy && a.edges == b.edges && a.orders == b.orders && a.nonNegative == b.nonNegative
-            && sameEfficiency(a.efficiency, b.efficiency)
+            && sameEfficiency(a.efficiency, b.efficiency) && a.edgeStep == b.edgeStep
     }
 
     /// Beam energy E0 in keV (the continuum ends there).
@@ -75,22 +91,27 @@ public nonisolated struct ContinuumForm: Sendable, Equatable {
     /// roll-off is modelled and P_s only carries what epsilon does not. nil: the round-1 form. The model is generic
     /// (`DetectorEfficiency`), not a measured Super-X calibration.
     package var efficiency: (any DetectorEfficiency)?
+    /// How the step at each edge is drawn: sharp at the edge channel (default) or broadened by the detector resolution.
+    package var edgeStep: ContinuumEdgeStep
 
     package static let alKEdge = FitPolicy.alKEdgeKeV
 
     package init(beamEnergy: Double, edges: [Double] = [ContinuumForm.alKEdge], orders: [Int] = [9, 5],
-                 nonNegative: Bool = true, efficiency: (any DetectorEfficiency)? = nil) {
+                 nonNegative: Bool = true, efficiency: (any DetectorEfficiency)? = nil,
+                 edgeStep: ContinuumEdgeStep = .sharp) {
         precondition(orders.count == edges.count + 1, "ContinuumForm needs one order per segment")
         precondition(zip(edges, edges.dropFirst()).allSatisfy { $0 < $1 }, "edges ascend")
         self.beamEnergy = beamEnergy; self.edges = edges; self.orders = orders; self.nonNegative = nonNegative
-        self.efficiency = efficiency
+        self.efficiency = efficiency; self.edgeStep = edgeStep
     }
 
     package var coefficientCount: Int { orders.reduce(0) { $0 + $1 + 1 } }
 
     /// The basis columns over `energies` (ascending, the fitted channels). A coefficient with no
     /// support (an empty segment) is not produced. `names` label them for the result.
-    package func columns(energies: [Double]) -> (values: [[Double]], names: [String]) {
+    /// `resolutionMnKaEV`: the detector resolution that broadens each edge step (the line model's). nil, or `edgeStep == .sharp`,
+    /// or no edge inside the range: the sharp form, value for value the form before broadening existed.
+    package func columns(energies: [Double], resolutionMnKaEV: Double? = nil) -> (values: [[Double]], names: [String]) {
         guard let lo = energies.first, let hi = energies.last, hi > lo else { return ([], []) }
         // Segment boundaries inside the fitted range.
         var bounds = [lo]
@@ -100,6 +121,9 @@ public nonisolated struct ContinuumForm: Sendable, Equatable {
         let eRef = 2.0
         let epsRef = efficiency?.efficiency(energyKeV: eRef) ?? 1
         let kRef = (beamEnergy - eRef) / eRef * (epsRef > 0 ? epsRef : 1)
+        if edgeStep == .broadened, let r = resolutionMnKaEV, bounds.count > 2 {
+            return broadenedColumns(energies: energies, bounds: bounds, segmentIndex: segmentIndex, resolutionMnKaEV: r, kRef: kRef)
+        }
         var values: [[Double]] = [], names: [String] = []
         for s in 0..<(bounds.count - 1) {
             let a = bounds[s], b = bounds[s + 1]
@@ -127,6 +151,47 @@ public nonisolated struct ContinuumForm: Sendable, Equatable {
         return (values, names)
     }
 
+    /// The broadened form (header): B(E) = K(E) eps(E) sum_s w_s(E) P_s(clamp(t_s, 0, 1)) / kRef, w_s = Phi_a - Phi_b.
+    private func broadenedColumns(energies: [Double], bounds: [Double], segmentIndex: [Int], resolutionMnKaEV r: Double,
+                                  kRef: Double) -> (values: [[Double]], names: [String]) {
+        let last = bounds.count - 2
+        func sigma(_ e: Double) -> Double {
+            // A width law with no real width at the edge (a resolution far below the line-width law's floor) falls back to the
+            // Mn K-alpha width, as the proposer does.
+            let fwhm = XRayLines.fwhm(resolutionMnKaEV: r, atEnergy: e) ?? r / 1000
+            return max(fwhm, 1e-9) / (8 * log(2.0)).squareRoot()
+        }
+        func phi(_ z: Double) -> Double { 0.5 * erfc(-z / 2.0.squareRoot()) }
+        var values: [[Double]] = [], names: [String] = []
+        for s in 0...last {
+            let a = bounds[s], b = bounds[s + 1]
+            // The order follows the segment's own channels, as in the sharp form.
+            let own = energies.filter { $0 >= a && (s == last ? $0 <= b : $0 < b) }.count
+            guard own > 0 else { continue }
+            let order = min(orders[segmentIndex[s]], own - 1)
+            let sa = s > 0 ? sigma(a) : 0, sb = s < last ? sigma(b) : 0
+            // Phi(lower edge) - Phi(upper edge): the weights telescope to 1 for any number of edges (a product would over-count
+            // between two close edges); max(0, ...) guards the far tails where the lower edge's sigma is the smaller.
+            let weight = energies.map { e -> Double in
+                max(0, (s > 0 ? phi((e - a) / sa) : 1) - (s < last ? phi((e - b) / sb) : 0))
+            }
+            for k in 0...order {
+                var col = [Double](repeating: 0, count: energies.count)
+                for (i, e) in energies.enumerated() where weight[i] > 0 {
+                    let t = min(max((e - a) / (b - a), 0), 1)
+                    let basis = nonNegative
+                        ? binomial(order, k) * pow(t, Double(k)) * pow(1 - t, Double(order - k))
+                        : pow(t, Double(k))
+                    let eps = efficiency?.efficiency(energyKeV: e) ?? 1
+                    col[i] = ((beamEnergy - e) / e) * eps / kRef * basis * weight[i]
+                }
+                values.append(col)
+                names.append("continuum s\(segmentIndex[s]) t^\(k)")
+            }
+        }
+        return (values, names)
+    }
+
     private func binomial(_ n: Int, _ k: Int) -> Double {
         var r = 1.0
         for i in 0..<k { r = r * Double(n - i) / Double(i + 1) }
@@ -137,7 +202,8 @@ public nonisolated struct ContinuumForm: Sendable, Equatable {
     package var label: String {
         let eps = efficiency.map { "\u{00D7}\u{03B5}(\($0.sourceName))" } ?? ""
         let ord = orders.map(String.init).joined(separator: ",")
-        let split = edges.isEmpty ? "no edge split" : "split at " + edges.map { "\($0) keV" }.joined(separator: ", ")
+        let step = edgeStep == .broadened ? " (step broadened by the detector resolution)" : ""
+        let split = edges.isEmpty ? "no edge split" : "split at " + edges.map { "\($0) keV" }.joined(separator: ", ") + step
         return "continuum: Kramers\(eps)\u{00D7}\(nonNegative ? "Bernstein" : "poly")(\(ord)), \(split); orders chosen on synthetic data"
     }
 
@@ -154,3 +220,11 @@ public nonisolated struct ContinuumForm: Sendable, Equatable {
         + "(generator A, no tail: continuum\u{2013}line projection) and \u{2212}310 \u{00B1} 23 at 882 (generator S); "
         + "\u{2212}402 \u{00B1} 33 at 3000 (A) and \u{2212}310 \u{00B1} 23 at 2899 (S)"
 }
+
+/// How the continuum's step at an absorption edge is drawn (`FitPolicy.continuumEdgeStep`).
+public nonisolated enum ContinuumEdgeStep: String, Sendable, Codable, Equatable {
+    /// Seen through the detector: each segment weighted by the detector's Gaussian step at its edges.
+    case broadened
+    /// A jump between the two channels either side of the edge (the default).
+    case sharp
+}

Sources/SpectroEDX/Fit/EDSFit.swift

@@ -69,7 +69,7 @@ public nonisolated struct FitSettings: Sendable, Equatable {
                     fitTo: defaultFitTo(axis: axis, beamEnergy: beamEnergy, upperLimitKeV: policy.defaultUpperLimitKeV),
                     background: .continuum(ContinuumForm(beamEnergy: beamEnergy, edges: policy.continuumEdges,
                                                          orders: policy.continuumOrders, nonNegative: policy.continuumNonNegative,
-                                                         efficiency: efficiency)),
+                                                         efficiency: efficiency, edgeStep: policy.continuumEdgeStep)),
                     escapePeaks: policy.escapePeaks, lineTail: policy.lineTail)
     }
 }

Sources/SpectroEDX/Fit/LinearDesign.swift

@@ -106,7 +106,7 @@ package nonisolated struct LinearDesign: Sendable {
                 bg.append(Column(kind: .background, name: "polynomial \(k)", values: v, signed: true, supported: true))
             }
         case .continuum(let form):
-            let (cols, names) = form.columns(energies: energies)
+            let (cols, names) = form.columns(energies: energies, resolutionMnKaEV: model.resolutionMnKaEV)
             for (v, nm) in zip(cols, names) {
                 bg.append(Column(kind: .background, name: nm, values: v.map { $0 * axis.scale },
                                  signed: !form.nonNegative, supported: true))

Sources/SpectroEDX/Quant/QuantificationPolicy.swift

@@ -73,6 +73,11 @@ public nonisolated struct FitPolicy: Sendable, Equatable, Codable {
     public var continuumEdges: [Double]
     /// Bernstein basis with non-negative coefficients (the continuum cannot go negative). False: free-sign monomials.
     public var continuumNonNegative: Bool
+    /// How the continuum steps at each edge. Default `.sharp`, the form the registered synthetic results were measured with.
+    /// `.broadened` blurs the step by the detector resolution, as a measured edge is: on a pure Al spectrum (Titan, Super-X) the
+    /// sharp step left the model 39 % above the data in 1.65-2.2 keV, the broadened one 19 %; on a synthetic Al film with a
+    /// broadened edge a weak Si K-alpha read -38.7 % sharp and -25.4 % broadened.
+    public var continuumEdgeStep: ContinuumEdgeStep
     /// A detector efficiency multiplying the whole continuum form (its low-energy roll-off modelled instead of absorbed in the
     /// polynomial). Default nil: the Kramers-only form.
     public var continuumEfficiency: DetectorLayers?
@@ -111,26 +116,51 @@ public nonisolated struct FitPolicy: Sendable, Equatable, Codable {
                 irls: IRLSSettings = IRLSSettings(), escapeRatioWarning: Double = 0.05,
                 axisRefinement: AxisRefinementSettings = AxisRefinementSettings(),
                 lineTail: LineTail? = nil, sumPeaks: SumPeakRule = .whenSignificant, sumPeakSignificance: Double = 3,
-                continuumDeadLayer: ContinuumDeadLayer? = nil) {
+                continuumEdgeStep: ContinuumEdgeStep = .sharp, continuumDeadLayer: ContinuumDeadLayer? = nil) {
         self.continuumDeadLayer = continuumDeadLayer
         self.fitFromKeV = fitFromKeV; self.defaultUpperLimitKeV = defaultUpperLimitKeV
         self.continuumOrders = continuumOrders; self.continuumEdges = continuumEdges
         self.continuumNonNegative = continuumNonNegative; self.continuumEfficiency = continuumEfficiency
         self.escapePeaks = escapePeaks; self.maximumPolynomialOrder = maximumPolynomialOrder
         self.irls = irls; self.escapeRatioWarning = escapeRatioWarning; self.axisRefinement = axisRefinement
         self.lineTail = lineTail; self.sumPeaks = sumPeaks; self.sumPeakSignificance = sumPeakSignificance
+        self.continuumEdgeStep = continuumEdgeStep
     }
 
     /// The fit before sum peaks were the default (one Gaussian per line, no sum peaks): exSpy's line model.
     public static let exspyEquivalent = FitPolicy(sumPeaks: .off)
 
     private enum CodingKeys: String, CodingKey {
         case fitFromKeV, defaultUpperLimitKeV, continuumOrders, continuumEdges, continuumNonNegative, continuumEfficiency, escapePeaks,
-             lineTail, maximumPolynomialOrder, irls, escapeRatioWarning, axisRefinement, sumPeaks, sumPeakSignificance, continuumDeadLayer
+             lineTail, maximumPolynomialOrder, irls, escapeRatioWarning, axisRefinement, sumPeaks, sumPeakSignificance, continuumEdgeStep,
+             continuumDeadLayer
+    }
+
+    /// The synthesized encoding, except that `continuumEdgeStep` is written only when it is not the default, so the default
+    /// policy's canonical JSON stays byte for byte what it was before the field existed.
+    public func encode(to encoder: any Encoder) throws {
+        var c = encoder.container(keyedBy: CodingKeys.self)
+        try c.encode(fitFromKeV, forKey: .fitFromKeV)
+        try c.encode(defaultUpperLimitKeV, forKey: .defaultUpperLimitKeV)
+        try c.encode(continuumOrders, forKey: .continuumOrders)
+        try c.encode(continuumEdges, forKey: .continuumEdges)
+        try c.encode(continuumNonNegative, forKey: .continuumNonNegative)
+        try c.encodeIfPresent(continuumEfficiency, forKey: .continuumEfficiency)
+        try c.encode(escapePeaks, forKey: .escapePeaks)
+        try c.encodeIfPresent(lineTail, forKey: .lineTail)
+        try c.encode(maximumPolynomialOrder, forKey: .maximumPolynomialOrder)
+        try c.encode(irls, forKey: .irls)
+        try c.encode(escapeRatioWarning, forKey: .escapeRatioWarning)
+        try c.encode(axisRefinement, forKey: .axisRefinement)
+        try c.encode(sumPeaks, forKey: .sumPeaks)
+        try c.encode(sumPeakSignificance, forKey: .sumPeakSignificance)
+        if continuumEdgeStep != .sharp { try c.encode(continuumEdgeStep, forKey: .continuumEdgeStep) }
+        try c.encodeIfPresent(continuumDeadLayer, forKey: .continuumDeadLayer)
     }
 
     /// A policy stored before sum peaks were the default has no `sumPeaks`: it reads `.off`, the fit it was stored with (and gets the
-    /// non-default footer line). Every other field is required, as before.
+    /// non-default footer line). One stored before `continuumEdgeStep` existed reads `.sharp`, the form it was run with. Every other
+    /// field is required, as before.
     public init(from decoder: any Decoder) throws {
         let c = try decoder.container(keyedBy: CodingKeys.self)
         self.init(fitFromKeV: try c.decode(Double.self, forKey: .fitFromKeV),
@@ -147,6 +177,7 @@ public nonisolated struct FitPolicy: Sendable, Equatable, Codable {
                   lineTail: try c.decodeIfPresent(LineTail.self, forKey: .lineTail),
                   sumPeaks: try c.decodeIfPresent(SumPeakRule.self, forKey: .sumPeaks) ?? .off,
                   sumPeakSignificance: try c.decodeIfPresent(Double.self, forKey: .sumPeakSignificance) ?? 3,
+                  continuumEdgeStep: try c.decodeIfPresent(ContinuumEdgeStep.self, forKey: .continuumEdgeStep) ?? .sharp,
                   continuumDeadLayer: try c.decodeIfPresent(ContinuumDeadLayer.self, forKey: .continuumDeadLayer))
     }
 

Tests/SpectroCoreTests/Fixtures/spectro-golden/demo-edx-tiny-dm4.txt

@@ -659,6 +659,7 @@ whole.quantification.fitSettings.background.continuum.edges[0] = 3ff8f41f212d773
 whole.quantification.fitSettings.background.continuum.orders: n 2 sum-bits e fnv 54226316c5c52fc9
 whole.quantification.fitSettings.background.continuum.nonNegative = true
 whole.quantification.fitSettings.background.continuum.efficiency = nil
+whole.quantification.fitSettings.background.continuum.edgeStep = .sharp
 whole.quantification.fitSettings.escapePeaks = true
 whole.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325
 whole.quantification.fitSettings.method = .leastSquares
@@ -900,6 +901,7 @@ region-locked-axis.quantification.fitSettings.background.continuum.edges[0] = 3f
 region-locked-axis.quantification.fitSettings.background.continuum.orders: n 2 sum-bits e fnv 54226316c5c52fc9
 region-locked-axis.quantification.fitSettings.background.continuum.nonNegative = true
 region-locked-axis.quantification.fitSettings.background.continuum.efficiency = nil
+region-locked-axis.quantification.fitSettings.background.continuum.edgeStep = .sharp
 region-locked-axis.quantification.fitSettings.escapePeaks = true
 region-locked-axis.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325
 region-locked-axis.quantification.fitSettings.method = .leastSquares

Tests/SpectroCoreTests/Fixtures/spectro-golden/eds-realistic-l-sim.txt

@@ -314,6 +314,7 @@ seed0.seven.quantification.fitSettings.background.continuum.edges[0] = 3ff8f41f2
 seed0.seven.quantification.fitSettings.background.continuum.orders: n 2 sum-bits e fnv 54226316c5c52fc9
 seed0.seven.quantification.fitSettings.background.continuum.nonNegative = true
 seed0.seven.quantification.fitSettings.background.continuum.efficiency = nil
+seed0.seven.quantification.fitSettings.background.continuum.edgeStep = .sharp
 seed0.seven.quantification.fitSettings.escapePeaks = true
 seed0.seven.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325
 seed0.seven.quantification.fitSettings.method = .leastSquares
@@ -685,6 +686,7 @@ seed0.seven-60keV.quantification.fitSettings.background.continuum.edges[0] = 3ff
 seed0.seven-60keV.quantification.fitSettings.background.continuum.orders: n 2 sum-bits e fnv 54226316c5c52fc9
 seed0.seven-60keV.quantification.fitSettings.background.continuum.nonNegative = true
 seed0.seven-60keV.quantification.fitSettings.background.continuum.efficiency = nil
+seed0.seven-60keV.quantification.fitSettings.background.continuum.edgeStep = .sharp
 seed0.seven-60keV.quantification.fitSettings.escapePeaks = true
 seed0.seven-60keV.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325
 seed0.seven-60keV.quantification.fitSettings.method = .leastSquares
@@ -983,6 +985,7 @@ seed0.four.quantification.fitSettings.background.continuum.edges[0] = 3ff8f41f21
 seed0.four.quantification.fitSettings.background.continuum.orders: n 2 sum-bits e fnv 54226316c5c52fc9
 seed0.four.quantification.fitSettings.background.continuum.nonNegative = true
 seed0.four.quantification.fitSettings.background.continuum.efficiency = nil
+seed0.four.quantification.fitSettings.background.continuum.edgeStep = .sharp
 seed0.four.quantification.fitSettings.escapePeaks = true
 seed0.four.quantification.fitSettings.referenceShapes.shapes: n 0 sum 0000000000000000 0 fnv cbf29ce484222325
 seed0.four.quantification.fitSettings.method = .leastSquares