One gain mechanism: frames analysed in their own units, the gain applied where a photon count is assumed

The gain of an input in detector units (SMV, marCCD, XTC, the XFEL formats) is measured once, before
anything reads the frames (RugnuxMeasureGain: photon-peak spacing and background variance per unit
mean), or given with --gain, and handed to the integration and the spot finding by one function
(ApplyRugnuxPixelGain). The rotation pre-scan no longer measures it, and the stills reader no longer
converts frames of a detector whose photons are not resolved: SACLA MPCCD, Rayonix and CCD frames are
integrated in ADU at the gain, rotation and stills alike.

Where a photon count was assumed, it is now taken at the gain (bit-identical at gain 1):
- integration: the background clip's one-photon floor, the sigma floor (one photon), the profile fit's
  background floor (0.01 photon), and a still's overload kept as a lower bound, which was left in ADU
  while every other reflection left the integrator in photons - the cause of the earlier loss of the
  engine route on clipped CCD stills (merged low resolution inflated by the gain, the rocking-width
  fit's bandwidth collapsing); CPU and CUDA, Metal through the shared CPU code;
- the adaptive spot threshold (Poisson arm and its one-photon floor), CPU and GPU;
- the pre-scan spot-width measurement (strong-spot SNR, the ring clip) behind the stills bandwidth and
  the rotation signal radius;
- the hot/dead pixel tests, CPU and GPU; the empty-aperture ring clip.
BackgroundPixelGain is gone: one pixel_gain.

A detector whose background resolves single photons (photon peaks: JUNGFRAU in keV, AGIPD) is a photon
counter and its frames are still counted in whole photons as read: analysed in keV at the gain instead,
a photon-counting still set lost 0.008 R-free (fewer spots, a noisier beam-centre trial) - the rounding
is the detector's own discrimination. --gain sets the value, the background still says which kind.

Validation, against rc178: photon-counting and g<=1 rotation p.hkl md5-identical (lyso, thau, cytc,
myob, insu, a small-molecule set, 5reo, 6cee, 6i3j, 2xfw, 6s1u, 8v4j); CCD g>1 at or within noise of rc178 (5uth
R_meas 0.113 -> 0.108). Stills R-free (0-cycle REFMAC): MPCCD sets 0.184 -> 0.185 and 0.199 -> 0.196,
Rayonix identical, photon-resolving sets back on their rc178 numbers.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit da8dbfc8147354b63e375f2d9521dffb80938b6d)
This commit is contained in:
leonarski_f committed 2026-10-23 17:56:51 +02:00
1 parent 15793a6bc2
commit 711f902aa7
30 files changed
+341 -267

No files matched your search

+71
View File
@@ -394,3 +394,74 @@ TEST_CASE("BraggIntegrationEngineCPU_PixelGainGivesPhotons", "[Integration][port
CHECK(out_ph[0].var_bkg == Catch::Approx(out_adu[0].var_bkg / 4.0f));
}
}
// The same frame, once in photons and once in a detector's units at 7 ADU per photon, integrates to the
// same reflections in photons when the engine is given the gain: the background clip, its one-photon
// floor, the variance floors, the profile fit's weights and a still's overload kept as a lower bound
// are all taken at the gain. A sparse background of 0-2 photons per pixel, where the floors bind, weak
// and strong spots, and one spot clipped at the saturation limit.
TEST_CASE("BraggIntegrationEngineCPU_AduIntegratesAsPhotons", "[Integration][portable]") {
constexpr int32_t g = 7, limit = 400;
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();
ImagePreprocessorBuffer photons(npixel);
uint32_t state = 12345;
for (size_t i = 0; i < npixel; ++i) {
state = state * 1664525u + 1013904223u;
const uint32_t u = state >> 24;
photons[i] = u < 200 ? 0 : (u < 245 ? 1 : 2);
}
std::vector<Reflection> predicted;
int n = 0;
for (int cy = 60; cy < 480; cy += 40)
for (int cx = 60; cx < 1000; cx += 40, ++n) {
const double amplitude = n == 7 ? 3000.0 : 2.0 + 15.0 * (n % 23);
for (int dy = -6; dy <= 6; ++dy)
for (int dx = -6; dx <= 6; ++dx) {
int32_t &v = photons[(cy + dy) * width + cx + dx];
v += static_cast<int32_t>(std::lround(amplitude * std::exp(-(dx * dx + dy * dy) / (2 * 1.3 * 1.3))));
if (v >= limit)
v = INT32_MAX;
}
Reflection r{};
r.h = n; r.k = 1; r.l = 1;
r.predicted_x = static_cast<float>(cx);
r.predicted_y = static_cast<float>(cy);
r.d = 2.0f;
r.prescaling_corr = 1.0f;
r.partiality = 1.0f;
predicted.push_back(r);
}
ImagePreprocessorBuffer adu(npixel);
for (size_t i = 0; i < npixel; ++i)
adu[i] = photons[i] == INT32_MAX ? INT32_MAX : g * photons[i];
for (const auto mode : {IntegratorMode::BoxSum, IntegratorMode::ProfileGaussian}) {
BraggIntegrationSettings settings;
settings.Integrator(mode);
experiment.ImportBraggIntegrationSettings(settings);
experiment.Detector().SaturationLimit(limit);
BraggIntegrationEngineCPU in_photons(experiment, PixelMask(experiment));
const auto ref = in_photons.Run(photons, predicted, predicted.size(), 0);
settings.PixelGain(static_cast<float>(g));
experiment.ImportBraggIntegrationSettings(settings);
experiment.Detector().SaturationLimit(g * limit);
BraggIntegrationEngineCPU in_adu(experiment, PixelMask(experiment));
const auto out = in_adu.Run(adu, predicted, predicted.size(), 0);
REQUIRE(out.size() == ref.size());
REQUIRE(ref.size() > 100);
CHECK(in_adu.Counts().overloads_kept == 1);
CHECK(in_photons.Counts().overloads_kept == 1);
for (size_t i = 0; i < ref.size(); ++i) {
CHECK(out[i].h == ref[i].h);
CHECK(out[i].I == Catch::Approx(ref[i].I).epsilon(1e-4).margin(1e-3));
CHECK(out[i].sigma == Catch::Approx(ref[i].sigma).epsilon(1e-4));
CHECK(out[i].bkg == Catch::Approx(ref[i].bkg).epsilon(1e-4).margin(1e-5));
// 1/den - I cancels: float rounding of den shows at 1e-4 of the difference.
CHECK(out[i].var_bkg == Catch::Approx(ref[i].var_bkg).epsilon(1e-3).margin(1e-3));
CHECK(out[i].overloaded == ref[i].overloaded);
}
}
}