rugnux: the gain of a stills file that states none is read off its single-photon peak where the background resolves photons

A pedestal-subtracted JUNGFRAU background is mostly empty pixels; the variance per unit mean leaves
out the zeros and read 6.7 keV per photon on a 12 keV run in keV, 2.4 on a 4.6 keV one, so such
files needed --gain. MeasureGain now also histograms the sampled pixel values and takes the gain from
the single-photon peak (background_gain::PhotonPeak): the highest bin past the zero peak's descent,
standing 5 sigma above the valley before it, narrow (below half height a quarter of its value away),
and repeated by a two-photon peak near twice its value. Where there is no such pair the variance per
unit mean decides as before. Float frames are read for the measurement at a sixteenth of their unit,
so a photon of a few keV is resolved rather than rounded to whole keV.

Measured: float32 keV JUNGFRAU files 11.96 (photon energy 12.04 keV) and 4.53 (4.57 keV); the int16
photon copies show no peak past the first bin and keep the 0.55 variance answer, i.e. stay in
photons; the SACLA MPCCD and PAL Rayonix files show humps, not photon peaks, and measure as before.
A gain the file states or --gain gives still wins. Rotation is untouched (the pre-scan's CCD gain
uses only the variance estimator).

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

image_analysis/bragg_integration/BackgroundGain.cpp

@@ -92,4 +92,44 @@ std::optional<float> Gain(const std::vector<Levels> &frames) {
     return static_cast<float>(ratio.back().first);
 }
 
+void AddToHistogram(const std::vector<int32_t> &image, const std::vector<uint32_t> &mask, int64_t saturation,
+                    std::vector<uint64_t> &histogram) {
+    const int64_t end = std::min<int64_t>(saturation, PEAK_HISTOGRAM_BINS);
+    for (size_t i = 0; i < image.size(); i++)
+        if (mask[i] == 0 && image[i] >= 0 && image[i] < end)
+            histogram[image[i]]++;
+}
+
+std::optional<float> PhotonPeak(const std::vector<uint64_t> &h) {
+    // A peak stands out of a valley where it is higher by five times its own counting noise.
+    const auto stands = [&](size_t valley, size_t peak) {
+        return static_cast<double>(h[peak]) - static_cast<double>(h[valley])
+               > 5.0 * std::sqrt(static_cast<double>(h[peak]));
+    };
+    // Down the zero peak to the first valley.
+    size_t v0 = 1;
+    while (v0 + 1 < h.size() && h[v0 + 1] <= h[v0])
+        v0++;
+    if (v0 + 1 >= h.size())
+        return std::nullopt;
+    const size_t p = std::max_element(h.begin() + v0, h.end()) - h.begin();
+    if (!stands(v0, p) || 5 * p / 2 >= h.size())
+        return std::nullopt;
+    // A photon is a narrow peak, its width the noise: below half its height a quarter of its value away.
+    // A CCD's background is a hump as wide as where it lies.
+    if (2 * h[3 * p / 4] >= h[p] || 2 * h[5 * p / 4] >= h[p])
+        return std::nullopt;
+    // The two-photon peak: the histogram rises again, near 2p, past the valley between p and 2p.
+    const size_t v2 = std::min_element(h.begin() + p, h.begin() + 2 * p + 1) - h.begin();
+    const size_t p2 = std::max_element(h.begin() + v2, h.begin() + 5 * p / 2 + 1) - h.begin();
+    if (!stands(v2, p2) || 4 * p2 < 7 * p || 4 * p2 > 9 * p)
+        return std::nullopt;
+    // Within the bin: the vertex of a parabola through the log counts of the peak and its neighbours.
+    const double a = std::log(std::max<double>(1.0, h[p - 1]));
+    const double b = std::log(std::max<double>(1.0, h[p]));
+    const double c = std::log(std::max<double>(1.0, h[p + 1]));
+    const double curvature = a - 2.0 * b + c;
+    return static_cast<float>(p + (curvature < 0.0 ? 0.5 * (a - c) / curvature : 0.0));
+}
+
 }  // namespace background_gain

image_analysis/bragg_integration/BackgroundGain.h

@@ -47,4 +47,21 @@ Levels MeasureFrame(const std::vector<int32_t> &image, size_t width, size_t heig
 // fewer than MIN_BINS levels holding MIN_PAIRS_PER_BIN pairs.
 std::optional<float> Gain(const std::vector<Levels> &frames);
 
+// The gain of a detector whose background resolves single photons (a JUNGFRAU with its pedestal
+// subtracted): the histogram of pixel values has a peak at zero photons, one at a single photon, g,
+// and another at two, 2g. The variance per unit mean above does not see it on a sparse background -
+// the empty pixels it leaves out are most of it - so the quantum is read off the histogram instead.
+constexpr size_t PEAK_HISTOGRAM_BINS = 4096;
+
+// Adds the frame's unmasked pixels with 0 <= value < min(saturation, PEAK_HISTOGRAM_BINS) to histogram
+// (PEAK_HISTOGRAM_BINS long, one bin per value).
+void AddToHistogram(const std::vector<int32_t> &image, const std::vector<uint32_t> &mask, int64_t saturation,
+                    std::vector<uint64_t> &histogram);
+
+// The value of the single-photon peak: the highest bin past the descent of the zero peak, standing out
+// of the valley before it, narrow, and repeated by a two-photon peak near twice its value; nothing where the
+// histogram shows no such pair - integer photon counts, whose single photon is the first bin itself, or
+// a detector whose noise hides the photons.
+std::optional<float> PhotonPeak(const std::vector<uint64_t> &histogram);
+
 }  // namespace background_gain

reader/JFJochXFELReader.cpp

@@ -775,6 +775,16 @@ void JFJochXFELReader::ApplyGain(float adu_per_photon) {
     Place();
 }
 
+void JFJochXFELReader::ClearGain() {
+    for (auto &p: geom_.panels)
+        p.adu_per_photon.reset();
+    Place();
+}
+
+bool JFJochXFELReader::StoredFloat() const {
+    return layout_.stored_float || !layout_.plain_type;
+}
+
 bool JFJochXFELReader::HasPhotonEnergy() const {
     return photon_energy_eV_ > 0;
 }

reader/JFJochXFELReader.h

@@ -159,9 +159,13 @@ public:
 
     // Converts to photons with this gain from here on (detector units per photon, every panel).
     void ApplyGain(float adu_per_photon);
+    // Back to detector units: undoes ApplyGain where nothing read states a gain.
+    void ClearGain();
 
     // False where no gain is known, so the images are in detector units, not photons.
     [[nodiscard]] bool CountsInPhotons() const;
+    // True where the frames are stored as floats, which are rounded to whole detector units as read.
+    [[nodiscard]] bool StoredFloat() const;
     // False where nothing read states the photon energy (a geometry file over frames of no known layout):
     // the caller gives it.
     [[nodiscard]] bool HasPhotonEnergy() const;

rugnux/rugnux_cli.cpp

@@ -275,7 +275,7 @@ void print_usage() {
     std::cout << "   --detector-distance <num>        Detector distance (mm). On stills it is held: neither searched nor refined (the cell absorbs any error). Without it, stills with a reference cell (-C) search the distance and refine it together with the cell, the file's value being one hypothesis" << std::endl;
     std::cout << "   --wavelength <num>               Wavelength (A)" << std::endl;
     std::cout << "   --pixel-size <num>               Pixel size (um)" << std::endl;
-    std::cout << "   --gain <num>                     Detector units per photon, for a file that states none (SACLA, CXI, XTC and --crystfel-geometry input; the frames are converted to photons with it). Default: for SACLA and CXI input, measured on the background of a few frames; XTC frames stay in detector units" << std::endl;
+    std::cout << "   --gain <num>                     Detector units per photon, for a file that states none (SACLA, CXI, XTC and --crystfel-geometry input; the frames are converted to photons with it). Default: for SACLA, CXI and --crystfel-geometry input, measured on a few frames - the single-photon peak of the pixel values where they resolve photons (for a JUNGFRAU file in keV, the photon energy), else the background's variance per unit mean; XTC frames stay in detector units" << std::endl;
     std::cout << "   --saturation <num>               Pixel value from which a pixel is overloaded, in the units the frames are read in (photons where they are converted). Replaces the file's limit and the one stills read off the frames" << std::endl;
     std::cout << "   --rot1 <num>                     PONI rotation 1 (rad)" << std::endl;
     std::cout << "   --rot2 <num>                     PONI rotation 2 (rad)" << std::endl;
@@ -714,14 +714,16 @@ void configure_offline_output(DiffractionExperiment &experiment, const std::stri
 
 namespace {
 
-// The gain of a detector whose file states none - ADU per photon, from the background's variance per
-// unit mean (BackgroundGain.h) on a few frames spread over the run - so a foreign stills file can be
-// converted to photons as it is read.
+// The gain of a detector whose file states none - ADU per photon, on a few frames spread over the run -
+// so a foreign stills file can be converted to photons as it is read: the single-photon peak of the
+// pixel values where the background resolves photons, else the background's variance per unit mean
+// (BackgroundGain.h).
 std::optional<float> MeasureGain(JFJochReader &reader) {
     const auto dataset = reader.GetDataset();
     const size_t n = reader.GetNumberOfImages();
     const size_t frames = std::min<size_t>(n, background_gain::GAIN_MAX_FRAMES);
     std::vector<background_gain::Levels> levels;
+    std::vector<uint64_t> histogram(background_gain::PEAK_HISTOGRAM_BINS, 0);
     JFJochReaderRawImage raw;
     std::vector<uint8_t> buffer;
     std::vector<int32_t> pixels;
@@ -739,7 +741,11 @@ std::optional<float> MeasureGain(JFJochReader &reader) {
             continue;
         levels.push_back(background_gain::MeasureFrame(pixels, w, h, dataset->pixel_mask->GetMask(),
                                                        dataset->experiment.GetSaturationLimit()));
+        background_gain::AddToHistogram(pixels, dataset->pixel_mask->GetMask(),
+                                        dataset->experiment.GetSaturationLimit(), histogram);
     }
+    if (const auto g = background_gain::PhotonPeak(histogram))
+        return g;
     return background_gain::Gain(levels);
 }
 
@@ -2068,17 +2074,29 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     // background clip widened by the gain the pre-scan measures, as for any CCD.
     // A gain at or below one detector unit per photon (a CCD read out with a low gain: one PAL Rayonix
     // set lost most of its indexing when its frames were scaled up by it) is left as it is, as for any CCD.
+    // Float frames are measured at a sixteenth of their unit, as they would otherwise be rounded to whole
+    // units first: a photon of a few keV in a file in keV is a peak only a few units wide.
     if (!gain_per_photon && (input_is_panel_table || input_is_cxi || input_is_crystfel) && !xfel_reader.CountsInPhotons()) {
-        if (const auto g = MeasureGain(reader); g && *g <= 1.0f)
-            logger.Info("Gain: {:.3f} detector units per photon, measured on the background - at or below 1, the "
-                        "frames stay in detector units", *g);
-        else if (g) {
+        const float unit = xfel_reader.StoredFloat() ? 1.0f / 16.0f : 1.0f;
+        if (unit != 1.0f)
+            xfel_reader.ApplyGain(unit);
+        auto g = MeasureGain(reader);
+        if (g)
+            *g *= unit;
+        if (g && *g > 1.0f) {
             logger.Info("Gain: {:.3f} detector units per photon, measured on the background; the frames are "
                         "converted to photons with it", *g);
             xfel_reader.ApplyGain(*g);
-        } else
-            logger.Info("Gain: the file states none and the background does not show it; the frames stay in "
-                        "detector units");
+        } else {
+            if (unit != 1.0f)
+                xfel_reader.ClearGain();
+            if (g)
+                logger.Info("Gain: {:.3f} detector units per photon, measured on the background - at or below 1, "
+                            "the frames stay in detector units", *g);
+            else
+                logger.Info("Gain: the file states none and the background does not show it; the frames stay in "
+                            "detector units");
+        }
     }
 
     const auto dataset = reader.GetDataset();

tests/BackgroundGainTest.cpp

@@ -73,3 +73,45 @@ TEST_CASE("BackgroundGain_settings") {
     CHECK(s.GetBackgroundPixelGain() == 4.5f);
     CHECK_THROWS_AS(s.BackgroundPixelGain(0.5f), JFJochException);
 }
+
+namespace {
+
+// A background of mean photons per pixel, each worth gain units, with Gaussian noise of sigma units, the
+// pedestal subtracted.
+std::vector<uint64_t> PhotonHistogram(double gain, double sigma, uint32_t seed, double mean = 0.2) {
+    std::mt19937 rng(seed);
+    std::poisson_distribution<int> photons(mean);
+    std::normal_distribution<double> noise(0.0, sigma);
+    std::vector<int32_t> img(W * H);
+    for (auto &v: img)
+        v = static_cast<int32_t>(std::lround(gain * photons(rng) + noise(rng)));
+    std::vector<uint64_t> histogram(background_gain::PEAK_HISTOGRAM_BINS, 0);
+    background_gain::AddToHistogram(img, std::vector<uint32_t>(W * H, 0), 1 << 30, histogram);
+    return histogram;
+}
+
+}
+
+TEST_CASE("BackgroundGain_PhotonPeak") {
+    // A photon of 4.57 keV in a file in keV, read at a sixteenth of a keV, 0.3 keV of noise.
+    const auto g = background_gain::PhotonPeak(PhotonHistogram(4.57 * 16, 0.3 * 16, 1));
+    REQUIRE(g.has_value());
+    CHECK_THAT(*g / 16.0, Catch::Matchers::WithinAbs(4.57, 0.05));
+
+    // Whole detector units, 12 per photon.
+    const auto g12 = background_gain::PhotonPeak(PhotonHistogram(12.0, 1.0, 2));
+    REQUIRE(g12.has_value());
+    CHECK_THAT(*g12, Catch::Matchers::WithinAbs(12.0, 0.2));
+
+    // Integer photon counts: the single photon is the first bin, no peak past the zero peak.
+    CHECK_FALSE(background_gain::PhotonPeak(PhotonHistogram(1.0, 0.0, 3)).has_value());
+
+    // 30 photons per pixel at 4 units each: one hump, with no second peak at twice its value.
+    CHECK_FALSE(background_gain::PhotonPeak(PhotonHistogram(4.0, 3.0, 4, 30.0)).has_value());
+
+    // A CCD's background (bias + gain * photons + read noise): a hump, not a pair of photon peaks.
+    std::vector<uint32_t> mask;
+    std::vector<uint64_t> ccd(background_gain::PEAK_HISTOGRAM_BINS, 0);
+    background_gain::AddToHistogram(Frame(4.0, 5, mask), mask, 1 << 30, ccd);
+    CHECK_FALSE(background_gain::PhotonPeak(ccd).has_value());
+}