integration: remove --reciprocal-profile (proven per-frame dead-end)

An 11-crystal mosaicity-stratified re-test (/data/rotation_test, off vs on
vs a de-contaminated variant, plus a per-frame dump of the fitted widths)
showed the dial is net-negative and cannot work in the per-frame paradigm:

- The C|q|^2 mosaicity term - the whole point - is unfittable per-frame: the
  fitted curvature a2 comes out ~0 (often negative) on every crystal, with zero
  correlation to the XDS mosaicity (0.09..0.42 deg). Strong spots sit at low q
  where eta^2 q^2 is invisible; the curvature only appears at high q where there
  are ~0 strong spots. The law degenerates to a straight line.
- With a2~0 the high-res width becomes a blind 1/cos^2(2theta) extrapolation,
  2-4x wider than per-shell. The per-shell path's high-res "starvation" (flat
  narrow fallback) is accidentally correct: weak, crowded high-res spots want a
  narrow aperture, not the true wide spot shape.
- The over-wide profile pulls background into weak spots -> R-meas rises, CC1/2
  drops in reliable high-multiplicity shells (pding4_001, pding4_003, MyoB,
  EcwtCQ066). A cap at the widest well-sampled per-shell width recovers the
  regression, confirming over-widening is the harm. No crystal reliably wins;
  the apparent overall-CC gains were all in noise shells (mult 2-3, CC<20%).

Delete the CLI flag, the BraggIntegrationSettings::reciprocal_profile setting,
and the per-frame fit block. Default (per-shell) integration is byte-identical.
NEXTGEN_INTEGRATOR.md records the finding as a dead-end for posterity.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

common/BraggIntegrationSettings.cpp

@@ -77,15 +77,6 @@ IntegratorMode BraggIntegrationSettings::GetIntegrator() const {
     return integrator_mode;
 }
 
-BraggIntegrationSettings &BraggIntegrationSettings::ReciprocalProfile(bool input) {
-    reciprocal_profile = input;
-    return *this;
-}
-
-bool BraggIntegrationSettings::GetReciprocalProfile() const {
-    return reciprocal_profile;
-}
-
 float BraggIntegrationSettings::GetR1() const {
     return r_1;
 }

common/BraggIntegrationSettings.h

@@ -20,7 +20,6 @@ class BraggIntegrationSettings {
     float d_min_limit_A = 1.0;
     std::optional<float> fixed_profile_radius;
     float minimum_sigma_in_regards_to_i = 0.02;
-    bool reciprocal_profile = false;
 
 public:
     BraggIntegrationSettings& R1(float input);
@@ -29,7 +28,6 @@ public:
     BraggIntegrationSettings& DMinLimit_A(float input);
     BraggIntegrationSettings& FixedProfileRadius_recipA(std::optional<float> input);
     BraggIntegrationSettings& Integrator(IntegratorMode input);
-    BraggIntegrationSettings& ReciprocalProfile(bool input);
 
 
     [[nodiscard]] IntegratorMode GetIntegrator() const;
@@ -40,5 +38,4 @@ public:
     [[nodiscard]] float GetDMinLimit_A() const;
 
     [[nodiscard]] float GetMinimumSigmaInRegardsToI() const;
-    [[nodiscard]] bool GetReciprocalProfile() const;
 };

docs/JFJOCH_PROCESS.md

@@ -171,5 +171,4 @@ Integration:
 | --- | --- |
 | `--integrator <txt>` | Spot integrator: `gaussian` (profile-fit, default) \| `empirical` \| `boxsum` (classical fallback) |
 | `--integration-radius <r>` | Signal-box radius `r1`, or `r1,r2,r3` (px). One value ⇒ `r2=r1+2`, `r3=r1+4` |
-| `--reciprocal-profile` | Learn one global reciprocal-space profile width (`A+B·|q|+C·|q|²`) instead of per-shell; helps mosaic/sparse data |
 | `--bandwidth <num>` | Relative X-ray bandwidth FWHM (e.g. `0.01` for a 1% DMM); default from file or 0 (monochromatic) |

image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md

@@ -322,33 +322,41 @@ Open follow-up (deferred, untested): a neighbour mask on the learning accumulati
 alternative to the `r1` window (keeps this spot's genuine [r1,r2] tail; needs a multiplicity mask to
 separate self from neighbour in the overlap annulus).
 
-### Reciprocal-space profile width — `--reciprocal-profile` (experimental dial, default OFF)
+### Reciprocal-space profile width — `--reciprocal-profile` (tried and REMOVED 2026-07-02)
 
-The per-shell width is learned in **pixels**, so it changes ~4× with resolution (mostly the geometric
-projection of a roughly-constant reciprocal-space relrod) and must be binned per shell to have enough
-strong spots — which **starves** at high resolution / on sparse data. Mapping the spot 2nd-moment into
-**reciprocal space** (multiply the pixel tangential moment by the Jacobian `g_tan² = cos²2θ`) removes the
-projection and makes the width far more transferable: across resolution shells the spread collapses from
-CV ≈ 44 % (pixel) to ≈ 11 % radial / 27 % tangential, and the tangential residual is a clean mosaicity
-law `σ²_q,tan = A + B|q| + C|q|²` (`C|q|²` = relrod variance `~(η|q|)²`).
+**The idea.** The per-shell width is learned in **pixels**, so it changes ~4× with resolution (the
+geometric projection of a roughly-constant reciprocal-space relrod) and must be binned per shell, which
+**starves** at high resolution / on sparse data. Mapping the tangential 2nd-moment into **reciprocal
+space** (Jacobian `g_tan² = cos²2θ`) should remove the projection and give a transferable mosaicity law
+`σ²_q,tan = A + B|q| + C|q|²` (`C|q|²` = relrod variance `~(η|q|)²`), fit once per frame — like XDS/DIALS
+integrating in a reciprocal profile coordinate. It was shipped as an off-by-default dial after being
+tuned on the single sharp HEWL crystal, where it was metric-neutral.
 
-`--reciprocal-profile` replaces the per-shell pixel width with **one global** fit of that law (`B`,`C`
-clamped ≥ 0; falls back linear→constant). On the **sharp** HEWL test crystal (η = 0.091°, so the `C|q|²`
-mosaicity term fits to noise/≈0 and it reduces to the validated linear form) it is **metric-neutral**:
-ISa 16.2→16.3, anomalous 0.92×→0.91× (no drop), CCref band 90.0→89.9, CC1/2 a touch lower (1.12 Å
-95→93) — because per-shell isn't starved here (~23 k spots/shell) and a 2–3-parameter global fit is just
-less flexible, while the faithful pure-tangential width is slightly tighter than the (generous) isotropic
-per-shell seed the de-biased fit mildly prefers. So on a sharp crystal it is **simpler + more transferable
-at a small CC1/2 cost, ISa/anomalous unchanged** — kept as an off-by-default dial, not the default.
+**Why it was removed — an 11-crystal mosaicity-stratified re-test (`/data/rotation_test`, off vs on vs a
+de-contaminated variant, plus a per-frame dump of the fitted widths) showed it is net-negative and cannot
+work per-frame.** Three linked, proven facts:
+1. **The `C|q|²` mosaicity term — the whole point — is unfittable per-frame.** The fitted curvature `a2`
+   came out ≈ 0 (often slightly *negative*) on **every** crystal, with **zero correlation** to the XDS
+   mosaicity (0.09→0.42°). Strong spots live at low `q` (q̄ ≈ 0.24–0.47) where `η²q²` is invisible; the
+   curvature only appears at high `q`, where there are ~0 strong spots to constrain it. The law degenerates
+   to a straight line.
+2. With `a2 ≈ 0`, the high-res width becomes a **blind `1/cos²2θ` extrapolation → 2–4× wider than
+   per-shell** (e.g. px²: pding4_001 1.07 flat → 2.0–3.0; cytC 1.10 → 2.7; MyoB 0.89 → 3.4). Meanwhile the
+   per-shell path *starves* at high res (n_strong<30 ⇒ flat global fallback that even **shrinks** vs
+   mid-res) — which turns out to be **accidentally correct**: weak, crowded high-res spots want a *narrow*
+   aperture, not the true (wide) spot shape.
+3. The over-wide profile pulls background into weak spots ⇒ **per-observation scatter (R-meas) rises**,
+   CC1/2 drops in reliable, high-multiplicity shells (pding4_001 mult 25, pding4_003 mult 12, EcwtCQ066,
+   MyoB); the mean intensities survive (CC-overall flat). A confirmatory cap (clamp the reciprocal width at
+   the widest well-sampled per-shell value) **recovers the regression** — direct proof the over-widening is
+   the harm. The only "improvements" seen (EcwtAL500, EP_cs_01-17 overall CC) were entirely in
+   noise-dominated shells (mult 2–3, CC<20 %, R-meas>100 %). **No crystal reliably wins.**
 
-**Where it should pay off, and why it can't be shown here:** a **mosaic** crystal makes the `C|q|²` term
-large and dominant at high resolution, where (a) per-shell starves on the wide, weak high-res spots,
-(b) 6 discrete shells are too coarse for the steep width growth, and (c) the wide spots overflow the
-fixed grid (the ellipse path grows to hold them). This is exactly why DIALS/XDS integrate in a reciprocal
-profile coordinate. Both lyso test crystals are sharp (and the jets are sharp serial stills), so the dial
-can only be **demonstrated on a genuinely mosaic rotation dataset** — the reason it ships as a dial to try
-elsewhere. Tested and dropped along the way: a separate **radial** relrod fit (radial is divergence +
-rotation-smear, not mosaicity → no metric gain on the sharp crystal).
+The per-frame paradigm structurally cannot fit mosaicity, and un-starving the high-res width is net harm,
+so the dial (CLI flag, `BraggIntegrationSettings::reciprocal_profile`, and the per-frame fit block) was
+deleted. A real win would need cross-frame / rolling-window pooling to constrain `a2`, or a two-pass
+mosaicity seed from scaling — both against the per-frame / online + stills philosophy. Also dropped along
+the way: a separate **radial** relrod fit (radial is divergence + rotation-smear, not mosaicity → no gain).
 
 ## Lineage: PixelRefine (removed 2026-06-25)
 

image_analysis/bragg_integration/ProfileIntegrate2D.cpp

@@ -17,10 +17,6 @@ namespace {
 constexpr int N_SHELL = 6;
 constexpr double STRONG_I_OVER_SIGMA = 5.0;
 constexpr int MIN_STRONG_PER_SHELL = 30;
-// --reciprocal-profile width fit: gentle ridge on the curvature coefficient (sharp-crystal prior,
-// so curvature must be earned) and IRLS robust-fit iterations (Huber, to reject outlier spots).
-constexpr double RECIP_RIDGE = 0.05;
-constexpr int RECIP_IRLS_ITERS = 3;
 
 // Radial parallax broadening as the coefficient of tan^2(2theta), i.e. Var(z)/pixel^2 [px^2]. A photon
 // converts at a random depth z in the sensor (exponential with attenuation length L, truncated at the
@@ -288,89 +284,6 @@ std::vector<Reflection> ProfileIntegrateInternal(const DiffractionExperiment &ex
     const double F_px = geom.GetDetectorDistance_mm() / std::max(1e-6f, geom.GetPixelSize_mm());
     const bool use_ellipse = !empirical && (bw_sigma > 0.0 || c_radial > 0.0);
 
-    // Reciprocal-space profile width (--reciprocal-profile): a per-frame model of the tangential
-    // variance in reciprocal space, y(q) = a0 + a1*t + a2*t^2 with t = (q - qbar)/qscale, replacing the
-    // per-shell pixel width. The Jacobian g_tan = cos(2theta) maps the pixel tangential moment into
-    // reciprocal space (removing the ~4x geometric growth with resolution); the t^2 term is the crystal
-    // MOSAICITY (relrod variance ~ (eta*|q|)^2), ~0 for a sharp crystal. Fitting in the CENTERED,
-    // STANDARDIZED variable t rather than raw q keeps the 3x3 normal matrix well-conditioned even when
-    // the strong spots span a narrow q-range (small cell / sparse still) - the raw {1,q,q^2} fit went
-    // near-singular there, letting tiny per-frame jitter swing the curvature into a wild over-wide
-    // profile. The fit is robust (IRLS / Huber: outlier spots can't drag it) with a gentle ridge on the
-    // curvature (sharp crystal = prior). Applied per reflection as sigma2_tan,px = y(q)/g_tan^2, with q
-    // clamped to the fitted strong-spot range (never extrapolated).
-    const bool recip_on = settings.GetReciprocalProfile();
-    double rp_a0 = 0.0, rp_a1 = 0.0, rp_a2 = 0.0, rp_qbar = 0.0, rp_qscale = 1.0, rp_tmin = 0.0, rp_tmax = 0.0;
-    bool use_recip = false;
-    if (recip_on && !empirical) {
-        std::vector<double> qv, yv;
-        qv.reserve(npredicted); yv.reserve(npredicted);
-        for (size_t i = 0; i < npredicted; ++i) {
-            const auto &rh = rough[i];
-            if (!rh.ok || !rh.strong || rh.I <= 0.0 || !(predicted[i].d > 0.0f)) continue;
-            const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
-            const double Rpx = std::hypot(rx, ry);
-            if (Rpx < 1e-6) continue;
-            const double ux = rx / Rpx, uy = ry / Rpx;
-            double m2 = 0.0, m2w = 0.0;
-            for (int dy = -R; dy <= R; ++dy)
-                for (int dx = -R; dx <= R; ++dx) {
-                    if (dx * dx + dy * dy >= r1_sq) continue;
-                    const int64_t x = rh.cx + dx, y = rh.cy + dy;
-                    if (x < 0 || y < 0 || x >= static_cast<int64_t>(xpixel) || y >= static_cast<int64_t>(ypixel)) continue;
-                    const auto px = ptr[y * xpixel + x];
-                    if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue;
-                    const double w = std::max(0.0, (static_cast<double>(px) - rh.bkg) / rh.I);
-                    const double tn = -dx * uy + dy * ux;
-                    m2 += w * tn * tn; m2w += w;
-                }
-            if (m2w <= 0.0) continue;
-            const double tan2t = Rpx / F_px, cos2t = 1.0 / std::sqrt(1.0 + tan2t * tan2t);
-            const double y = cos2t * cos2t * (m2 / m2w);
-            if (!(y > 0.0) || !std::isfinite(y)) continue;
-            qv.push_back(1.0 / predicted[i].d); yv.push_back(y);
-        }
-        const size_t n = qv.size();
-        if (n >= 30) {
-            double qmin = qv[0], qmax = qv[0], qsum = 0.0;
-            for (double q : qv) { qmin = std::min(qmin, q); qmax = std::max(qmax, q); qsum += q; }
-            rp_qbar = qsum / static_cast<double>(n);
-            double var = 0.0; for (double q : qv) var += (q - rp_qbar) * (q - rp_qbar);
-            rp_qscale = std::sqrt(var / static_cast<double>(n));
-            if (!(rp_qscale > 1e-9)) rp_qscale = std::max(1e-6, 0.5 * (qmax - qmin));
-            rp_tmin = (qmin - rp_qbar) / rp_qscale; rp_tmax = (qmax - rp_qbar) / rp_qscale;
-            auto det3 = [](double a, double b, double c, double d, double e, double f, double g, double h, double i) {
-                return a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g);
-            };
-            std::vector<double> wob(n, 1.0);
-            for (int it = 0; it < RECIP_IRLS_ITERS; ++it) {
-                double S0=0,S1=0,S2=0,S3=0,S4=0, T0=0,T1=0,T2=0;
-                for (size_t i = 0; i < n; ++i) {
-                    const double t = (qv[i] - rp_qbar) / rp_qscale, t2 = t*t, w = wob[i];
-                    S0 += w; S1 += w*t; S2 += w*t2; S3 += w*t2*t; S4 += w*t2*t2;
-                    T0 += w*yv[i]; T1 += w*t*yv[i]; T2 += w*t2*yv[i];
-                }
-                const double S4r = S4 + RECIP_RIDGE * S4;   // ridge: shrink curvature toward the sharp prior
-                const double D = det3(S0,S1,S2, S1,S2,S3, S2,S3,S4r);
-                if (!(std::fabs(D) > 1e-12)) { use_recip = false; break; }
-                rp_a0 = det3(T0,S1,S2, T1,S2,S3, T2,S3,S4r) / D;
-                rp_a1 = det3(S0,T0,S2, S1,T1,S3, S2,T2,S4r) / D;
-                rp_a2 = det3(S0,S1,T0, S1,S2,T1, S2,S3,T2) / D;
-                use_recip = true;
-                if (it + 1 < RECIP_IRLS_ITERS) {   // Huber reweight from residuals (scale = 1.4826*MAD)
-                    std::vector<double> res(n);
-                    for (size_t i = 0; i < n; ++i) {
-                        const double t = (qv[i]-rp_qbar)/rp_qscale;
-                        res[i] = std::fabs(yv[i] - (rp_a0 + rp_a1*t + rp_a2*t*t));
-                    }
-                    std::vector<double> tmp(res); std::nth_element(tmp.begin(), tmp.begin()+n/2, tmp.end());
-                    const double mad = std::max(1e-6, 1.4826 * tmp[n/2]);
-                    for (size_t i = 0; i < n; ++i) { const double r = res[i]/(1.5*mad); wob[i] = r <= 1.0 ? 1.0 : 1.0/r; }
-                }
-            }
-        }
-    }
-
     // --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). ---
     std::vector<Reflection> out;
     out.reserve(npredicted);
@@ -386,13 +299,7 @@ std::vector<Reflection> ProfileIntegrateInternal(const DiffractionExperiment &ex
         const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
         const double Rpx = std::hypot(rx, ry);
         const double tan2t = Rpx / F_px;
-        double s2t = shell_sigma2[sh];
-        if (use_recip) {   // per-frame reciprocal width instead of the per-shell pixel width
-            const double q = 1.0 / std::max(predicted[i].d, 1e-6f);
-            const double t = std::clamp((q - rp_qbar) / rp_qscale, rp_tmin, rp_tmax);
-            const double cos2t = 1.0 / std::sqrt(1.0 + tan2t * tan2t);
-            s2t = std::max(0.25, (rp_a0 + rp_a1 * t + rp_a2 * t * t) / (cos2t * cos2t));
-        }
+        const double s2t = shell_sigma2[sh];
         double s2r = s2t;
         double ux = 1.0, uy = 0.0;
         bool elong = false;

tools/jfjoch_process.cpp

@@ -85,7 +85,6 @@ void print_usage() {
     std::cout << "   --bandwidth <num>               Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl;
     std::cout << "   --integration-radius <r>        Signal-box radius r1, or r1,r2,r3 (px). One value => r2=r1+2, r3=r1+4" << std::endl;
     std::cout << "   --integrator <txt>              Spot integrator boxsum|gaussian|empirical (default: gaussian profile-fit; boxsum is the classical fallback)" << std::endl;
-    std::cout << "   --reciprocal-profile            Learn one global reciprocal-space profile width (A+B|q|+C|q|^2) instead of per-shell; helps mosaic/sparse data" << std::endl;
     }
 
 enum {
@@ -112,7 +111,6 @@ enum {
     OPT_CAPTURE_UNCERTAINTY,
     OPT_MOSAICITY,
     OPT_SMOOTH_G,
-    OPT_RECIPROCAL_PROFILE,
     OPT_DETECT_ICE_RINGS,
     OPT_NO_SCALE_FULLS,
     OPT_WRITE_PROCESS_H5,
@@ -166,7 +164,6 @@ static option long_options[] = {
     {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH},
     {"integration-radius", required_argument, nullptr, OPT_INTEGRATION_RADIUS},
     {"integrator", required_argument, nullptr, OPT_INTEGRATOR},
-    {"reciprocal-profile", no_argument, nullptr, OPT_RECIPROCAL_PROFILE},
     {"detect-ice-rings", no_argument, nullptr, OPT_DETECT_ICE_RINGS},
     {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS},
     {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF},
@@ -345,7 +342,6 @@ int main(int argc, char **argv) {
     std::optional<float> d_min_scale_merge;
     std::optional<std::string> integration_radius_arg; // "r1" or "r1,r2,r3"
     std::optional<IntegratorMode> integrator_mode;     // --integrator boxsum|gaussian|empirical
-    bool reciprocal_profile = false;                   // --reciprocal-profile
     std::optional<double> outlier_reject_nsigma;       // merge per-observation outlier rejection
     std::optional<double> delta_cchalf_nsigma;         // per-crystal CC1/2-delta rejection
 
@@ -581,9 +577,6 @@ int main(int argc, char **argv) {
             else if (strcmp(optarg, "empirical") == 0) integrator_mode = IntegratorMode::ProfileEmpirical;
             else { logger.Error("--integrator expects boxsum|gaussian|empirical"); return 1; }
             break;
-        case OPT_RECIPROCAL_PROFILE:
-            reciprocal_profile = true;
-            break;
         case OPT_REJECT_OUTLIERS:
             outlier_reject_nsigma = std::stod(optarg);
             break;
@@ -860,13 +853,6 @@ int main(int argc, char **argv) {
                     : "profile (empirical)");
     }
 
-    if (reciprocal_profile) {
-        BraggIntegrationSettings bis = experiment.GetBraggIntegrationSettings();
-        bis.ReciprocalProfile(true);
-        experiment.ImportBraggIntegrationSettings(bis);
-        logger.Info("Reciprocal-space global profile width enabled (per-shell width replaced)");
-    }
-
     SpotFindingSettings spot_settings;
     spot_settings.enable = true;
     spot_settings.indexing = true;