From 87560311bc4c5d8aa5f74040aa71d2d95890bdb7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 22 Oct 2026 20:22:15 +0200 Subject: [PATCH] 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 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- .../bragg_integration/BackgroundGain.cpp | 40 ++++++++++++++++++ .../bragg_integration/BackgroundGain.h | 17 ++++++++ reader/JFJochXFELReader.cpp | 10 +++++ reader/JFJochXFELReader.h | 4 ++ rugnux/rugnux_cli.cpp | 40 +++++++++++++----- tests/BackgroundGainTest.cpp | 42 +++++++++++++++++++ 6 files changed, 142 insertions(+), 11 deletions(-) diff --git a/image_analysis/bragg_integration/BackgroundGain.cpp b/image_analysis/bragg_integration/BackgroundGain.cpp index b6732ec19..cec5b36da 100644 --- a/image_analysis/bragg_integration/BackgroundGain.cpp +++ b/image_analysis/bragg_integration/BackgroundGain.cpp @@ -92,4 +92,44 @@ std::optional Gain(const std::vector &frames) { return static_cast(ratio.back().first); } +void AddToHistogram(const std::vector &image, const std::vector &mask, int64_t saturation, + std::vector &histogram) { + const int64_t end = std::min(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 PhotonPeak(const std::vector &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(h[peak]) - static_cast(h[valley]) + > 5.0 * std::sqrt(static_cast(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(1.0, h[p - 1])); + const double b = std::log(std::max(1.0, h[p])); + const double c = std::log(std::max(1.0, h[p + 1])); + const double curvature = a - 2.0 * b + c; + return static_cast(p + (curvature < 0.0 ? 0.5 * (a - c) / curvature : 0.0)); +} + } // namespace background_gain diff --git a/image_analysis/bragg_integration/BackgroundGain.h b/image_analysis/bragg_integration/BackgroundGain.h index de9c9ae75..bc990ee0a 100644 --- a/image_analysis/bragg_integration/BackgroundGain.h +++ b/image_analysis/bragg_integration/BackgroundGain.h @@ -47,4 +47,21 @@ Levels MeasureFrame(const std::vector &image, size_t width, size_t heig // fewer than MIN_BINS levels holding MIN_PAIRS_PER_BIN pairs. std::optional Gain(const std::vector &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 &image, const std::vector &mask, int64_t saturation, + std::vector &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 PhotonPeak(const std::vector &histogram); + } // namespace background_gain diff --git a/reader/JFJochXFELReader.cpp b/reader/JFJochXFELReader.cpp index f0da13dfa..db84b06ec 100644 --- a/reader/JFJochXFELReader.cpp +++ b/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; } diff --git a/reader/JFJochXFELReader.h b/reader/JFJochXFELReader.h index 868e5210f..063520004 100644 --- a/reader/JFJochXFELReader.h +++ b/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; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e628abd00..e14304e83 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -275,7 +275,7 @@ void print_usage() { std::cout << " --detector-distance 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 Wavelength (A)" << std::endl; std::cout << " --pixel-size Pixel size (um)" << std::endl; - std::cout << " --gain 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 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 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 PONI rotation 1 (rad)" << std::endl; std::cout << " --rot2 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 MeasureGain(JFJochReader &reader) { const auto dataset = reader.GetDataset(); const size_t n = reader.GetNumberOfImages(); const size_t frames = std::min(n, background_gain::GAIN_MAX_FRAMES); std::vector levels; + std::vector histogram(background_gain::PEAK_HISTOGRAM_BINS, 0); JFJochReaderRawImage raw; std::vector buffer; std::vector pixels; @@ -739,7 +741,11 @@ std::optional 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(); diff --git a/tests/BackgroundGainTest.cpp b/tests/BackgroundGainTest.cpp index ab9ac9db9..cf6e28af5 100644 --- a/tests/BackgroundGainTest.cpp +++ b/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 PhotonHistogram(double gain, double sigma, uint32_t seed, double mean = 0.2) { + std::mt19937 rng(seed); + std::poisson_distribution photons(mean); + std::normal_distribution noise(0.0, sigma); + std::vector img(W * H); + for (auto &v: img) + v = static_cast(std::lround(gain * photons(rng) + noise(rng))); + std::vector histogram(background_gain::PEAK_HISTOGRAM_BINS, 0); + background_gain::AddToHistogram(img, std::vector(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 mask; + std::vector 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()); +}