diff --git a/Sources/SpectroEDX/Fit/Continuum.swift b/Sources/SpectroEDX/Fit/Continuum.swift index fd1cb8f..3942e90 100644 --- a/Sources/SpectroEDX/Fit/Continuum.swift +++ b/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.. Double { + guard e >= 0.1 else { return 0 } + let table = MassAbsorptionTable.bundled! + let mu = (try? table.mac(element: "Al", energyKeV: e)) ?? 0 + return 9.0 * (beam - e) / e * exp(-pow(0.35 / e, 2.2)) * exp(-massThickness * mu) + } + + /// `incoming` convolved with the detector's Gaussian, sigma from the line-width law at each energy (counts per channel). + private static let continuum: [Double] = { + let raw = energies.map(incoming) + return energies.indices.map { i in + let e = energies[i] + let fwhm = XRayLines.fwhm(resolutionMnKaEV: resolution, atEnergy: max(e, 0.3)) ?? resolution / 1000 + let sigma = fwhm / (8 * log(2.0)).squareRoot() + let reach = Int((5 * sigma / axis.scale).rounded(.up)) + var sum = 0.0, wsum = 0.0 + for j in max(0, i - reach)...min(axis.size - 1, i + reach) { + let w = exp(-0.5 * pow((energies[j] - e) / sigma, 2)) + sum += w * raw[j]; wsum += w + } + return sum / wsum + } + }() + + private static func lineColumn(_ id: String) -> [Double] { + let el = String(id.prefix { $0 != "_" }) + let g = EDSLineModel.build(elements: [el], axis: axis, beamEnergy: beam, resolutionMnKaEV: resolution).groups.first { $0.id == id }! + return EDSLineModel.column(of: g.lines, energies: energies, scale: axis.scale) + } + + private static let alArea = 3e5, siArea = 2e3 + + private static func spectrum(si: Bool) -> [Double] { + let al = lineColumn("Al_Ka"), s = lineColumn("Si_Ka") + return energies.indices.map { continuum[$0] + alArea * al[$0] + (si ? siArea * s[$0] : 0) } + } + + private static func fit(_ counts: [Double], elements: [String], step: ContinuumEdgeStep) throws -> EDSFitResult { + var policy = FitPolicy() + policy.continuumEdgeStep = step + let s = FitSettings.standard(elements: elements, axis: axis, resolutionMnKaEV: resolution, beamEnergy: beam, policy: policy) + return try EDSFit.run(counts: counts, axis: axis, settings: s, computeCovariance: false) + } + + /// Data minus model summed over [lo, hi) keV, and the data there. + private static func band(_ f: EDSFitResult, _ counts: [Double], _ lo: Double, _ hi: Double) -> (residual: Double, data: Double) { + var r = 0.0, d = 0.0 + for (j, c) in f.channels.enumerated() where energies[c] >= lo && energies[c] < hi { r += f.residual[j]; d += counts[c] } + return (r, d) + } + + private static func worstChannel(_ f: EDSFitResult, _ counts: [Double]) -> Double { + var worst = 0.0 + for (j, c) in f.channels.enumerated() where energies[c] >= 1.50 && energies[c] < 1.70 { worst = max(worst, abs(f.residual[j]) / counts[c]) } + return worst + } + + /// Pure Al, no Si. Broadened: the model in 1.65-2.2 keV (under Si K-alpha) within 0.3 % of the data, and no channel in + /// 1.50-1.70 keV off by 4 % or more. The sharp step fails both on the same spectrum (measured: 0.95 % and 5.5 %), which is the + /// defect this form exists for. Measured broadened: 0.06 % and 3.3 % (the rest is the order-5 segment above the edge). + @Test func testTheBroadenedStepFollowsAResolutionBroadenedAlEdgeAndTheSharpDoesNot() throws { + let counts = Self.spectrum(si: false) + let broad = try Self.fit(counts, elements: ["Al"], step: .broadened) + let sharp = try Self.fit(counts, elements: ["Al"], step: .sharp) + let ab = Self.band(broad, counts, 1.65, 2.2), asharp = Self.band(sharp, counts, 1.65, 2.2) + let wb = Self.worstChannel(broad, counts), ws = Self.worstChannel(sharp, counts) + print(String(format: "EDGE pure Al 1.65-2.2 keV: broadened %+.2f %%, sharp %+.2f %%; worst channel 1.50-1.70: broadened %.2f %%, sharp %.2f %%", + ab.residual / ab.data * 100, asharp.residual / asharp.data * 100, wb * 100, ws * 100)) + #expect(abs(ab.residual) / ab.data < 0.003, "broadened 1.65-2.2 keV: \(ab.residual) of \(ab.data)") + #expect(wb < 0.04, "broadened worst channel \(wb)") + #expect(abs(asharp.residual) / asharp.data >= 0.003 || ws >= 0.04, "the sharp step misfits the broadened edge: \(asharp.residual), \(ws)") + } + + /// A weak Si K-alpha (2000 counts against 3e5 Al K-alpha) above the broadened edge: the sharp step reads it 38.7 % low, the + /// broadened one 25.4 % low (measured). Held: the broadened bias is at most 0.75 of the sharp one, and both are pinned within + /// 2 points so a change either way is looked at. The remaining bias is the order-5 segment above the edge (orders (9, 9) read + /// -8 %), stated, not asserted away; results stay unvalidated. + @Test func testAWeakSiLineAboveTheBroadenedEdgeReadsCloserWithTheBroadenedStep() throws { + let counts = Self.spectrum(si: true) + func si(_ step: ContinuumEdgeStep) throws -> Double { + let f = try Self.fit(counts, elements: ["Al", "Si"], step: step) + return f.values[try #require(f.groupIDs.firstIndex(of: "Si_Ka"))] / Self.siArea - 1 + } + let b = try si(.broadened), s = try si(.sharp) + print(String(format: "EDGE Si bias: broadened %+.2f %%, sharp %+.2f %%", b * 100, s * 100)) + #expect(abs(b) <= 0.75 * abs(s), "broadened \(b) vs sharp \(s)") + #expect(abs(b - -0.254) < 0.02, "broadened Si bias \(b)") + #expect(abs(s - -0.387) < 0.02, "sharp Si bias \(s)") + } + + /// `.sharp` is the default and the form without the broadening: the same columns as a call with no resolution. + @Test func testTheSharpStepIsTheDefaultAndTheUnbroadenedForm() { + let e = Array(Self.energies[40..<4000]) + let sharp = ContinuumForm(beamEnergy: Self.beam) + let broad = ContinuumForm(beamEnergy: Self.beam, edgeStep: .broadened) + #expect(sharp.edgeStep == .sharp) + #expect(sharp.columns(energies: e, resolutionMnKaEV: Self.resolution).values == sharp.columns(energies: e).values) + #expect(broad.columns(energies: e).values == sharp.columns(energies: e).values, "no resolution: the sharp form") + #expect(sharp.columns(energies: e, resolutionMnKaEV: Self.resolution).values != broad.columns(energies: e, resolutionMnKaEV: Self.resolution).values) + #expect(sharp != broad) + #expect(!sharp.label.contains("broadened") && broad.label.contains("broadened")) + // Without an edge inside the range the two forms are the same. + let low = Array(Self.energies[40..<300]) // 0.2-1.5 keV + #expect(sharp.columns(energies: low, resolutionMnKaEV: Self.resolution).values == broad.columns(energies: low, resolutionMnKaEV: Self.resolution).values) + } + + /// With every coefficient 1 the Bernstein sums are 1 on both segments, so the broadened continuum equals the sharp one + /// (Kramers / kRef): the weights sum to 1. Every column is non-negative and finite. + @Test func testTheBroadenedWeightsSumToOne() { + let e = Array(Self.energies[40..<4000]) + let broad = ContinuumForm(beamEnergy: Self.beam, edgeStep: .broadened).columns(energies: e, resolutionMnKaEV: Self.resolution).values + let sharp = ContinuumForm(beamEnergy: Self.beam).columns(energies: e).values + for i in e.indices { + let b = broad.reduce(0) { $0 + $1[i] }, s = sharp.reduce(0) { $0 + $1[i] } + #expect(abs(b - s) <= 1e-12 * max(abs(s), 1), "channel \(i): \(b) vs \(s)") + } + #expect(broad.allSatisfy { $0.allSatisfy { $0 >= 0 && $0.isFinite } }) + } + + /// Two close edges (90 eV apart, under three sigma): the weights still sum to 1 (they telescope; a product would over-count between). + @Test func testTheBroadenedWeightsSumToOneWithTwoCloseEdges() { + let e = Array(Self.energies[40..<4000]) + let edges = [1.5596, 1.65] + let broad = ContinuumForm(beamEnergy: Self.beam, edges: edges, orders: [9, 3, 5], edgeStep: .broadened) + .columns(energies: e, resolutionMnKaEV: Self.resolution).values + let sharp = ContinuumForm(beamEnergy: Self.beam, edges: edges, orders: [9, 3, 5]).columns(energies: e).values + for i in e.indices { + let b = broad.reduce(0) { $0 + $1[i] }, s = sharp.reduce(0) { $0 + $1[i] } + #expect(abs(b - s) <= 1e-12 * max(abs(s), 1), "channel \(i): \(b) vs \(s)") + } + #expect(broad.allSatisfy { $0.allSatisfy { $0 >= 0 && $0.isFinite } }) + } + + /// The policy: default `.sharp`, named in the footer when `.broadened`, and a stored policy without the field reads `.sharp`. + @Test func testThePolicyDefaultsToSharpAndOldStoredPoliciesReadIt() throws { + #expect(FitPolicy.default.continuumEdgeStep == .sharp) + // The default policy's canonical JSON is what it was before the field existed: the key is written only when not the default. + #expect(!QuantificationPolicy.default.canonicalJSON.contains("continuumEdgeStep")) + var broad = QuantificationPolicy() + broad.fit.continuumEdgeStep = .broadened + #expect(broad.footerLine?.contains("continuumEdgeStep = broadened") == true, "\(broad.footerLine ?? "nil")") + let round = try JSONDecoder().decode(QuantificationPolicy.self, from: Data(broad.canonicalJSON.utf8)) + #expect(round == broad) + var json = try #require(try JSONSerialization.jsonObject(with: Data(QuantificationPolicy.default.canonicalJSON.utf8)) as? [String: Any]) + var fit = try #require(json["fit"] as? [String: Any]) + fit.removeValue(forKey: "continuumEdgeStep") + json["fit"] = fit + let old = try JSONDecoder().decode(QuantificationPolicy.self, from: try JSONSerialization.data(withJSONObject: json)) + #expect(old == .default) + let settings = FitSettings.standard(elements: ["Al"], axis: Self.axis, resolutionMnKaEV: Self.resolution, beamEnergy: Self.beam, + policy: broad.fit) + guard case .continuum(let form) = settings.background else { Issue.record("no continuum"); return } + #expect(form.edgeStep == .broadened) + } +} diff --git a/docs/parity.md b/docs/parity.md index 314d903..cf5cc4c 100644 --- a/docs/parity.md +++ b/docs/parity.md @@ -275,6 +275,11 @@ the bounds are physics (`familyOvervoltage` 2, `sumPeakFraction` 1, not measured `IdentificationLoopTests` on planted spectra, `RealFileAutoIDTests` (opt-in) on three measured files. The identification fit's Kβ weights are NIST EPQ's (`EDSLineWeights.epqK`), not eXSpy's; the quantification fit keeps eXSpy's, so no parity row moves. +The continuum's step at an edge (`FitPolicy.continuumEdgeStep`): `.sharp` (default) cuts the segments at the edge channel; +`.broadened` (opt-in) weights each segment by the detector's Gaussian step at its edges, sigma from the line-width law at the edge. +No Python counterpart: exSpy's background is a polynomial with no edge. `validation: "none"`; tests `ContinuumEdgeBroadeningTests` +on a generated Al film whose edge is convolved with the detector (no measured spectrum committed). + Detector Si internal fluorescence: eXSpy 0.3.2 has no term for it (`EDSModel.add_family_lines` fits one Gaussian per X-ray line; no dead layer and no detector efficiency either, the k-factor carries the efficiency), so the Si internal-fluorescence counts at Si Kα land in its `Si_Ka` area. SpectroKit's default fit does the same; its `.countsAboveEdge` reference shape, which would remove them, is opt-in. The simulator