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

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.