integration: add --reciprocal-profile dial (global reciprocal-space width)

The per-shell profile width is learned in pixels, so it varies ~4x with resolution
(mostly the geometric projection of a near-constant reciprocal-space relrod) and must
be binned per shell -> it starves at high resolution / on sparse data. The new
--reciprocal-profile flag instead learns ONE global width in reciprocal space,
sigma2_q,tan = A + B|q| + C|q|^2: the Jacobian g_tan=cos(2theta) removes the geometric
projection, and C|q|^2 is the crystal mosaicity relrod (variance ~(eta|q|)^2). Applied
per reflection as sigma2_tan,px = (A + B|q| + C|q|^2)/g_tan^2 (B,C clamped >=0;
quadratic->linear->constant fallback).

Off by default. On the sharp HEWL test crystal (mosaicity 0.091deg, so C fits to ~0 and
it reduces to the validated linear form) it is metric-neutral: ISa 16.2->16.3, anomalous
0.92x unchanged, CCref band 90.0->89.9, CC1/2 a touch lower (per-shell isn't starved at
23k spots/shell, and a global fit is less flexible). So: simpler + more transferable at a
small CC1/2 cost, ISa/anomalous held. Its payoff is on MOSAIC crystals (large C|q|^2),
where per-shell starves on the wide weak high-res spots and 6 shells are too coarse; both
lyso test crystals are sharp, so it ships as a dial to try on mosaic data elsewhere.
A separate radial relrod fit was tried and dropped (no gain). See NEXTGEN_INTEGRATOR.md.

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

common/BraggIntegrationSettings.cpp

@@ -77,6 +77,15 @@ 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,6 +20,7 @@ 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);
@@ -28,6 +29,7 @@ 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;
@@ -38,4 +40,5 @@ public:
     [[nodiscard]] float GetDMinLimit_A() const;
 
     [[nodiscard]] float GetMinimumSigmaInRegardsToI() const;
+    [[nodiscard]] bool GetReciprocalProfile() const;
 };

image_analysis/bragg_integration/NEXTGEN_INTEGRATOR.md

@@ -322,6 +322,34 @@ 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)
+
+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|)²`).
+
+`--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.
+
+**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).
+
 ## Lineage: PixelRefine (removed 2026-06-25)
 
 `ProfileIntegrate2D` is the surviving half of an earlier experimental still-integrator, **PixelRefine**,

image_analysis/bragg_integration/ProfileIntegrate2D.cpp

@@ -284,6 +284,64 @@ 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 global profile width (--reciprocal-profile): one global model
+    // sigma2_q,tan = A + B*|q| + C*|q|^2 over all strong spots replaces the per-shell pixel width. The
+    // Jacobian g_tan = cos(2theta) maps the pixel tangential moment into reciprocal space, removing the
+    // geometric projection that makes the pixel width grow ~4x with resolution. The C*|q|^2 term is the
+    // crystal MOSAICITY (relrod variance ~ (eta*|q|)^2): ~0 for a sharp crystal (the fit collapses to
+    // A + B*|q|, tying the per-shell model) but dominant for a mosaic one, where per-shell starves on the
+    // wide, weak high-res spots. Applied per reflection as sigma2_tan,px = (A + B*|q| + C*|q|^2)/g_tan^2.
+    const bool recip_on = settings.GetReciprocalProfile();
+    double recip_A = 0.0, recip_B = 0.0, recip_C = 0.0;
+    bool use_recip = false;
+    if (recip_on && !empirical) {
+        double s1 = 0, sq = 0, sq2 = 0, sq3 = 0, sq4 = 0, sy = 0, sqy = 0, sq2y = 0;
+        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 q = 1.0 / predicted[i].d, yv = cos2t * cos2t * (m2 / m2w);
+            s1 += 1; sq += q; sq2 += q * q; sq3 += q * q * q; sq4 += q * q * q * q;
+            sy += yv; sqy += q * yv; sq2y += q * q * yv;
+        }
+        // Fit quadratic, falling back to linear then constant so B and C stay >= 0 (a relrod cannot shrink
+        // with resolution; a sharp crystal's C comes out as noise/<0 and drops to the linear branch).
+        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);
+        };
+        if (s1 >= 30) {
+            const double D3 = det3(s1, sq, sq2, sq, sq2, sq3, sq2, sq3, sq4), D2 = s1 * sq2 - sq * sq;
+            const double C3 = std::fabs(D3) > 1e-9 ? det3(s1, sq, sy, sq, sq2, sqy, sq2, sq3, sq2y) / D3 : -1.0;
+            const double B3 = std::fabs(D3) > 1e-9 ? det3(s1, sy, sq2, sq, sqy, sq3, sq2, sq2y, sq4) / D3 : 0.0;
+            const double B2 = std::fabs(D2) > 1e-9 ? (s1 * sqy - sq * sy) / D2 : -1.0;
+            if (std::fabs(D3) > 1e-9 && C3 > 0.0 && B3 >= 0.0) {
+                recip_A = det3(sy, sq, sq2, sqy, sq2, sq3, sq2y, sq3, sq4) / D3; recip_B = B3; recip_C = C3;
+            } else if (std::fabs(D2) > 1e-9 && B2 >= 0.0) {
+                recip_A = (sy * sq2 - sqy * sq) / D2; recip_B = B2;
+            } else {
+                recip_A = sy / s1;   // constant width fallback
+            }
+            use_recip = recip_A > 0.0;
+        }
+    }
+
     // --- Pass B: profile-fit each reflection (Kabsch, de-biased variance v = B + I*P; iterate). ---
     std::vector<Reflection> out;
     out.reserve(npredicted);
@@ -296,13 +354,19 @@ std::vector<Reflection> ProfileIntegrateInternal(const DiffractionExperiment &ex
 
         int Rf = R;
         const std::vector<double> *Pvec = &shell_P[sh];   // ProfileEmpirical uses the shared learned grid
-        double s2t = shell_sigma2[sh], s2r = s2t;
+        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) {   // global reciprocal width instead of the per-shell pixel width
+            const double q = 1.0 / std::max(predicted[i].d, 1e-6f), cos2t = 1.0 / std::sqrt(1.0 + tan2t * tan2t);
+            s2t = std::max(0.25, (recip_A + recip_B * q + recip_C * q * q) / (cos2t * cos2t));
+        }
+        double s2r = s2t;
         double ux = 1.0, uy = 0.0;
         bool elong = false;
         if (use_ellipse) {
-            const double rx = predicted[i].predicted_x - beam_x, ry = predicted[i].predicted_y - beam_y;
-            const double Rpx = std::hypot(rx, ry), sbw = bw_sigma * Rpx;
-            const double tan2t = Rpx / F_px;
+            const double sbw = bw_sigma * Rpx;
             const double radial_extra = sbw * sbw + c_radial * tan2t * tan2t;   // bandwidth + parallax + capture
             // Only elongate where the radial streak adds a genuine fraction of a pixel of variance; at
             // low/mid resolution the smear is sub-pixel and elongating just adds noise.

tools/jfjoch_process.cpp

@@ -72,6 +72,7 @@ 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 {
@@ -97,7 +98,8 @@ enum {
     OPT_SCALE_FULLS,
     OPT_CAPTURE_UNCERTAINTY,
     OPT_MOSAICITY,
-    OPT_SMOOTH_G
+    OPT_SMOOTH_G,
+    OPT_RECIPROCAL_PROFILE
 };
 
 static option long_options[] = {
@@ -142,6 +144,7 @@ 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},
     {"reject-outliers", required_argument, nullptr, OPT_REJECT_OUTLIERS},
     {"reject-delta-cchalf", required_argument, nullptr, OPT_REJECT_DELTA_CCHALF},
     {nullptr, 0, nullptr, 0}
@@ -314,6 +317,7 @@ 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
 
@@ -535,6 +539,9 @@ 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;
@@ -773,6 +780,13 @@ 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;