Bragg integration: a CCD rotation sweep's sigmas at the measured gain, intensities in photons

The pre-scan measures the gain of an SMV/marCCD detector (background variance per unit mean) but the
integrator used it only for the background clip: sigma, var_bkg, the strong-spot flag (I/sigma >= 5)
and the profile-vs-box-sum guard (10 box sigmas) all treated one ADU as one photon, so on a Saturn
(g 4.7-5.8) the guard sat at ~4.5 true sigma. The new BraggIntegrationSettings::PixelGain (set by
rugnux from the same measurement, rotation only, at max(1, g)) takes every variance the engine forms
at the gain, and Finalize converts I, sigma, bkg and var_bkg to photons, which is what the merge's
variance model (var_bkg + I) assumes. CPU and CUDA engines; Metal runs the CPU code for this. Gain 1
(every photon-counting input, ADSC/marCCD whose measured ratio is below 1, and stills) is bit-identical.

Battery --only, 12 open-arm CCD/background-dominated sets vs the base build (987e812e2), GPU:
the six Saturn sets (g 4.7-5.8) all improve - guard fallbacks 4.1-17.3% -> 1.2-6.1%,
anomalous peak height 3mc4 x1.10, 3p85 x1.02, 3r6o x1.06, 5cc8 x1.03, 5uth x1.05, 5vml x1.04;
R_meas 0.880-0.993x on all six; R_model shell-scaled better on 5/6 (5cc8 +0.0005); ISa 5uth +1.5%,
5vml +15%, 3p85 +3%, the twinned 3mc4/3r6o -5%/-1%. 6cee, 6i3j, 6v2r, 7kcn, 8v4j, 9jzo unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu

common/BraggIntegrationSettings.cpp

@@ -177,6 +177,17 @@ float BraggIntegrationSettings::GetBackgroundPixelGain() const {
     return bkg_pixel_gain;
 }
 
+BraggIntegrationSettings &BraggIntegrationSettings::PixelGain(float input) {
+    check_finite("Pixel gain", input);
+    check_min("Pixel gain", input, 1.0);
+    pixel_gain = input;
+    return *this;
+}
+
+float BraggIntegrationSettings::GetPixelGain() const {
+    return pixel_gain;
+}
+
 BraggIntegrationSettings &BraggIntegrationSettings::BackgroundRadialCorrection(std::optional<bool> input) {
     bkg_radial_correction = input;
     return *this;

common/BraggIntegrationSettings.h

@@ -103,6 +103,13 @@ class BraggIntegrationSettings {
     // there the photon clip cuts the background's upper tail at about 2 sigma and adds a positive offset
     // to every reflection. Set by rugnux from the pre-scan for SMV and marCCD input (BackgroundGain.h).
     float bkg_pixel_gain = 1.0f;
+    // Detector units per photon of the integrated intensities. Above 1 the variances the integrator
+    // forms - the box sum's, the profile fit's and the background's - are taken at this many ADU per
+    // photon, so the strong-spot selection and the profile/summation guard compare in photons, and
+    // each reflection leaves the integrator converted to photons (I, sigma, bkg, var_bkg), which is
+    // what the merge's own variance model (var_bkg + I) assumes. Set by rugnux for rotation sweeps of a
+    // CCD from the same measurement as bkg_pixel_gain; 1 for a detector that counts photons.
+    float pixel_gain = 1.0f;
     // Radial background curvature correction. The signal disk and the background annulus are
     // concentric, so for ANY background linear in position their means are equal - a plane fit buys
     // nothing and the leading error is the CURVATURE of the radial background, which the flat annulus
@@ -181,6 +188,7 @@ public:
     BraggIntegrationSettings& BackgroundTrimFraction(float input);
     BraggIntegrationSettings& BackgroundClipNSigma(float input);
     BraggIntegrationSettings& BackgroundPixelGain(float input);
+    BraggIntegrationSettings& PixelGain(float input);
     BraggIntegrationSettings& BackgroundRadialCorrection(std::optional<bool> input);
     BraggIntegrationSettings& MaxHKL(std::optional<int> input);
     BraggIntegrationSettings& Overlap(OverlapMode input);
@@ -202,6 +210,7 @@ public:
     [[nodiscard]] float GetBackgroundTrimFraction() const;
     [[nodiscard]] float GetBackgroundClipNSigma() const;
     [[nodiscard]] float GetBackgroundPixelGain() const;
+    [[nodiscard]] float GetPixelGain() const;
     // Unset = auto (gate per image on the smooth-ice score); see bkg_radial_correction.
     [[nodiscard]] std::optional<bool> GetBackgroundRadialCorrection() const;
     [[nodiscard]] std::optional<int> GetMaxHKL() const;

image_analysis/bragg_integration/BraggIntegrationEngine.cpp

@@ -118,6 +118,7 @@ BraggIntegrationEngine::BraggIntegrationEngine(const DiffractionExperiment &expe
     // longer overrides it for broadband data.
     bkg_clip_nsigma = settings.GetBackgroundClipNSigma();
     bkg_pixel_gain = settings.GetBackgroundPixelGain();
+    pixel_gain = settings.GetPixelGain();
     bkg_trim = bkg_clip_nsigma > 0.0f ? 0.0f : settings.GetBackgroundTrimFraction();
 
     // Overlap treatment. Ownership is decided out to the fit grid's half size, which is where the
@@ -265,10 +266,11 @@ std::vector<Reflection> BraggIntegrationEngine::Finalize(const std::vector<Refle
         if (!fr.ok)
             continue;
         Reflection refl = predicted[i];
-        refl.I = fr.I;
-        refl.sigma = fr.sigma;
-        refl.bkg = fr.bkg;
-        refl.var_bkg = fr.var_bkg;
+        // In photons: the engine measured in ADU, with variances already taken at the gain.
+        refl.I = fr.I / pixel_gain;
+        refl.sigma = fr.sigma / pixel_gain;
+        refl.bkg = fr.bkg / pixel_gain;
+        refl.var_bkg = fr.var_bkg / (pixel_gain * pixel_gain);
         if (fr.has_observed) {
             refl.observed_x = fr.observed_x;
             refl.observed_y = fr.observed_y;

image_analysis/bragg_integration/BraggIntegrationEngine.h

@@ -181,6 +181,9 @@ protected:
     double bw_sigma;             // bandwidth sigma [dimensionless, * Rpx -> px]
     float  bkg_clip_nsigma;      // high-outlier background sigma-clip multiplier (0 = no clip)
     float  bkg_pixel_gain = 1.0f;   // background pixel variance per unit mean the clip scales by
+    // ADU per photon (BraggIntegrationSettings::pixel_gain): every variance the engine forms is taken
+    // at it, and Finalize converts the reflections to photons with it. 1 changes nothing, bit for bit.
+    float  pixel_gain = 1.0f;
     bool   use_ellipse;          // radially elongate the per-reflection Gaussian
 
     double c_radial;             // radial variance coefficient of tan^2(2theta): parallax + capture

image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp

@@ -428,9 +428,13 @@ void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection
         // error enters n_inner times over: var(I) = I_sum + n_inner^2 * bkg/n_bkg_used. Leaving
         // the second term out understates sigma by sqrt(1 + n_inner/n_bkg) - 1.109x at the
         // default r1=4/r2=6/r3=10 stencil, on every reflection of every dataset.
+        // bkg_var is in counted units, as the profile fit reads it; every variance formed from it is
+        // taken at the gain, so a CCD's sigma - and the strong flag and the guard read from it - is in
+        // photons rather than in ADU counted as photons.
+        const double g = pixel_gain;
         out.bkg_var = out.bkg / n_bkg_used;
-        out.var_bkg = static_cast<double>(s.n_inner_valid) * out.bkg
-                      + static_cast<double>(s.n_inner_valid) * s.n_inner_valid * out.bkg_var;
+        out.var_bkg = g * (static_cast<double>(s.n_inner_valid) * out.bkg
+                           + static_cast<double>(s.n_inner_valid) * s.n_inner_valid * out.bkg_var);
         out.I_sum = s.I_sum;
         out.n_inner = static_cast<int>(s.n_inner_valid);
         out.n_disk = s.n_disk;
@@ -438,7 +442,7 @@ void BraggIntegrationEngineCPU::FinishPassA(PassA &a, size_t i, const Reflection
         const double var_bkg_term = static_cast<double>(s.n_inner_valid) * s.n_inner_valid * out.bkg_var;
         out.sigma = 1.0;
         if (s.I_sum > 0) {
-            out.sigma = std::max(out.sigma, std::sqrt(static_cast<double>(s.I_sum) + var_bkg_term));
+            out.sigma = std::max(out.sigma, std::sqrt(g * (static_cast<double>(s.I_sum) + var_bkg_term)));
             // The centroid of the SIGNAL, not of the disk. Weighting by the raw counts weights by
             // signal plus background, and the background is flat over a disk centred on the
             // PREDICTION - so its own centroid is the prediction exactly, and it pulls the answer
@@ -850,15 +854,17 @@ std::vector<Reflection> BraggIntegrationEngineCPU::ProfilePasses(const Sampler &
         // 1/den is the profile-fit variance with the background taken as exact. The fit is
         // I = sum(P*(px-bkg)/v) / sum(P^2/v), so dI/dbkg = -wsum/den and the background estimate's
         // own error adds (wsum/den)^2 * var(bkg) - the same term the box sum was missing.
-        double sigma = std::sqrt(1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var);
+        // Both variances at the gain (FinishPassA), so the guard below compares like with like.
+        const double g = pixel_gain;
+        double sigma = std::sqrt(g * (1.0 / den + (wsum / den) * (wsum / den) * rh.bkg_var));
         // var_bkg is the NON-SIGNAL part of that variance, and 1/den is the fit's variance taken at the
         // fitted intensity, so the signal part to remove is I itself - not max(0, I). Clamping it leaves
         // a down-fluctuated reflection's own (deflated) variance standing as its background variance,
         // which is 2-3x too small; the merge then weights exactly the down-fluctuated observations up.
         // The whole point of a separate var_bkg is that it does not move with the observation's own
         // fluctuation, and 1/den - I is what does not (1/den grows with I one for one).
-        double var_bkg = std::max(0.0, 1.0 / den - I
-                                       + (wsum / den) * (wsum / den) * rh.bkg_var);
+        double var_bkg = g * std::max(0.0, 1.0 / den - I
+                                           + (wsum / den) * (wsum / den) * rh.bkg_var);
         // The seed is a sum over the disk the box sum actually read, so when pixels are missing from
         // both - excluded to a neighbour, or unreadable - the fit's full-profile intensity has to be
         // scaled down to that same disk before the two are comparable. Nothing dropped gives 1.

image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu

@@ -20,6 +20,7 @@ struct BraggGpuParams {
     int   R, G, GG;
     float bkg_clip_nsigma;   // high-side background sigma-clip multiplier (0 = no clip)
     float bkg_pixel_gain;    // background pixel variance per unit mean (1 = counted photons)
+    float pixel_gain;        // ADU per photon every variance is taken at (see the CPU engine)
     int   empirical;      // ProfileEmpirical vs ProfileGaussian
     int   use_ellipse;
     float bw_sigma;
@@ -421,7 +422,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
     double sigma = 1.0;
     uint8_t hasobs = 0; double ox = 0.0, oy = 0.0;
     if (Isum > 0) {
-        sigma = fmax(sigma, sqrt((double) Isum + var_bkg_term));
+        sigma = fmax(sigma, sqrt((double) p.pixel_gain * ((double) Isum + var_bkg_term)));
         // The centroid of the SIGNAL. Weighting by raw counts weights by signal plus background,
         // and the background is flat over a disk centred on the PREDICTION, so its centroid is the
         // prediction and it pulls the answer there - the displacement comes out shrunk by
@@ -441,7 +442,7 @@ __global__ void boxsum(const float *px_x, const float *px_y, const BraggPanel *p
     cy_o[i] = (int) lroundf(cy);
     I_o[i] = (float) I; sigma_o[i] = (float) sigma; bkg_o[i] = (float) bkg;
     bkgvar_o[i] = (float) bkg_var;
-    varbkg_o[i] = (float) ((double) s_ninner_valid * bkg + var_bkg_term);
+    varbkg_o[i] = (float) ((double) p.pixel_gain * ((double) s_ninner_valid * bkg + var_bkg_term));
     isum_o[i] = (float) Isum;
     ninner_o[i] = s_ninner_valid;
     rbin_o[i] = min(max((int) lroundf(s_r0), 0), n_rad > 0 ? n_rad - 1 : 0);
@@ -828,9 +829,10 @@ __global__ void fit(const int32_t *img, const uint32_t *owner, const float *px_x
             // 1/s_den takes the background as exact; dI/dbkg = -s_wsum/s_den adds the background
             // estimate's own error (see the CPU engine).
             const float wr = s_wsum / s_den;
-            float I = s_I, sigma = sqrtf(1.0f / s_den + wr * wr * bkgvar_a[i]);
+            // Both at the gain, as the seed's (see the CPU engine).
+            float I = s_I, sigma = sqrtf(p.pixel_gain * (1.0f / s_den + wr * wr * bkgvar_a[i]));
             // The signal part to remove is I itself, not max(0, I) - see the CPU engine.
-            float var_bkg = fmaxf(0.0f, 1.0f / s_den - I + wr * wr * bkgvar_a[i]);
+            float var_bkg = p.pixel_gain * fmaxf(0.0f, 1.0f / s_den - I + wr * wr * bkgvar_a[i]);
             // Guard against profile-fit runaways (see the CPU engine): fall back to the summation seed
             // when the profile result diverges from it. Pixels missing from both - excluded to a
             // neighbour, or unreadable - scale the full-profile intensity back to the disk the seed
@@ -1012,6 +1014,7 @@ std::vector<Reflection> BraggIntegrationEngineGPU::Run(const ImagePreprocessorBu
         .R = R, .G = G, .GG = GG,
         .bkg_clip_nsigma = mode != IntegratorMode::BoxSum ? bkg_clip_nsigma : 0.0f,
         .bkg_pixel_gain = bkg_pixel_gain,
+        .pixel_gain = pixel_gain,
         .empirical = empirical ? 1 : 0,
         .use_ellipse = use_ellipse ? 1 : 0,
         .bw_sigma = static_cast<float>(bw_sigma), .c_radial = static_cast<float>(c_radial),

rugnux/RugnuxPreScan.cpp

@@ -606,6 +606,10 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         if (g) {
             BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings();
             bis.BackgroundPixelGain(std::max(1.0f, *g));
+            // On a rotation sweep the same gain sets the integrated intensities' variances (and their
+            // unit, photons). Stills keep their sigmas in ADU for now.
+            if (experiment_.IsRotationIndexing())
+                bis.PixelGain(std::max(1.0f, *g));
             experiment_.ImportBraggIntegrationSettings(bis);
         }
         if (!g)
@@ -614,8 +618,10 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru
         else
             logger.Info("Background gain: a background pixel scatters {:.2f}x its mean in ADU ({} frames){}",
                         *g, (ordinals.size() + gain_frame_stride - 1) / gain_frame_stride,
-                        *g > 1.0f ? fmt::format(" => background clip at mean + {:.1f}*sqrt({:.2f}*mean)",
-                                                experiment_.GetBraggIntegrationSettings().GetBackgroundClipNSigma(), *g)
+                        *g > 1.0f ? fmt::format(" => background clip at mean + {:.1f}*sqrt({:.2f}*mean){}",
+                                                experiment_.GetBraggIntegrationSettings().GetBackgroundClipNSigma(), *g,
+                                                experiment_.IsRotationIndexing()
+                                                        ? "; intensities and sigmas in photons at this gain" : "")
                                   : std::string(" - at or below 1, the photon clip is kept"));
     }
 

tests/BraggIntegrationEngineCPUTest.cpp

@@ -351,3 +351,46 @@ TEST_CASE("BraggIntegrationEngineCPU_StillKeepsOnlySaturatedOverloads", "[Integr
     CHECK(out[1].overloaded);
     CHECK(still.Counts().overloads_kept == 1);
 }
+
+// A detector that reports gain ADU per photon: the same spot read at gain 4 leaves the integrator in
+// photons - a quarter of the ADU intensity, at the Poisson sigma of that many photons, i.e. half the
+// sigma the ADU counted as photons gave.
+TEST_CASE("BraggIntegrationEngineCPU_PixelGainGivesPhotons", "[Integration][portable]") {
+    DiffractionExperiment experiment(DetJF(2));
+    experiment.DetectorDistance_mm(100.0f).IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(400.0f).BeamY_pxl(400.0f);
+    const size_t width = experiment.GetXPixelsNum(), npixel = experiment.GetPixelsNum();
+
+    const int cx = 600, cy = 300;
+    ImagePreprocessorBuffer image(npixel);
+    for (size_t i = 0; i < npixel; ++i)
+        image[i] = 40;
+    for (int dy = -6; dy <= 6; ++dy)
+        for (int dx = -6; dx <= 6; ++dx)
+            image[(cy + dy) * width + cx + dx] += static_cast<int32_t>(std::lround(800.0 * std::exp(-(dx * dx + dy * dy) / (2 * 1.3 * 1.3))));
+
+    Reflection r{};
+    r.h = 1; r.k = 2; r.l = 3;
+    r.predicted_x = cx; r.predicted_y = cy;
+    r.d = 2.0f;
+    r.prescaling_corr = 1.0f;
+    r.partiality = 1.0f;
+    const std::vector<Reflection> predicted{r};
+
+    for (const auto mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) {
+        BraggIntegrationSettings settings;
+        settings.Integrator(mode);
+        experiment.ImportBraggIntegrationSettings(settings);
+        BraggIntegrationEngineCPU adu(experiment, PixelMask(experiment));
+        const auto out_adu = adu.Run(image, predicted, 1, 0);
+        settings.PixelGain(4.0f);
+        experiment.ImportBraggIntegrationSettings(settings);
+        BraggIntegrationEngineCPU photons(experiment, PixelMask(experiment));
+        const auto out_ph = photons.Run(image, predicted, 1, 0);
+        REQUIRE(out_adu.size() == 1);
+        REQUIRE(out_ph.size() == 1);
+        CHECK(out_ph[0].I == Catch::Approx(out_adu[0].I / 4.0f));
+        CHECK(out_ph[0].sigma == Catch::Approx(out_adu[0].sigma / 2.0f));
+        CHECK(out_ph[0].bkg == Catch::Approx(out_adu[0].bkg / 4.0f));
+        CHECK(out_ph[0].var_bkg == Catch::Approx(out_adu[0].var_bkg / 4.0f));
+    }
+}

tests/BraggIntegrationEngineGPUTest.cpp

@@ -150,13 +150,15 @@ double CompareCpuVsGpu(IntegratorMode mode, std::optional<float> bandwidth_fwhm,
                      float r1 = 0.0f, float r2 = 0.0f, float r3 = 0.0f,
                      OverlapMode overlap = OverlapMode::Off, float companion_dx = 0.0f,
                      bool clip_spots = false, float odd_partiality = 1.0f,
-                     const SpotFootprint *footprint = nullptr) {
+                     const SpotFootprint *footprint = nullptr, float pixel_gain = 1.0f) {
     DiffractionExperiment experiment =
             MakeExperiment(mode, bandwidth_fwhm, clip_nsigma, radial, DetJF(2), stencil_k, r1, r2, r3,
                            overlap);
-    if (footprint) {
+    {
         BraggIntegrationSettings settings = experiment.GetBraggIntegrationSettings();
-        settings.Footprint(*footprint);
+        if (footprint)
+            settings.Footprint(*footprint);
+        settings.PixelGain(pixel_gain);
         experiment.ImportBraggIntegrationSettings(settings);
     }
     const size_t width = experiment.GetXPixelsNum();
@@ -324,6 +326,16 @@ TEST_CASE("BraggIntegrationEngineGPU_MatchesCPU") {
         CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 8, 0.0f,
                         0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 0.01f);
     }
+    // A CCD's gain: every variance at it, so the strong flag and the guard decide in photons, and the
+    // reflections leave in photons - computed independently in the two engines.
+    SECTION("ProfileGaussian pixel gain") {
+        CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 4.0f, false, 60, 0.0f,
+                        0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 1.0f, nullptr, 4.0f);
+    }
+    SECTION("BoxSum pixel gain") {
+        CompareCpuVsGpu(IntegratorMode::BoxSum, std::nullopt, 4.0f, false, 60, 0.0f,
+                        0.0f, 0.0f, 0.0f, OverlapMode::Off, 0.0f, false, 1.0f, nullptr, 4.0f);
+    }
     SECTION("ProfileGaussian mono trim") { CompareCpuVsGpu(IntegratorMode::ProfileGaussian, std::nullopt, 0.0f); }
     // Unreadable pixels inside the signal disks themselves: the MINPK rescue keeps the reflection and
     // fits it over what is left, and the peak-loss rule throws back the ones that lost the profile's