leonarski_fandClaude Opus 5
f60768d49c
Two changes to the same variance chain; they are in one commit because the second exists to remove an assumption the first was breaking, and separating them leaves a tree that is correct only by luck. The reported sigma was floored at 2% of the intensity, a per-partial I/sigma cap of 50. It applied only to the box-sum seed, never to the profile fit, so the shipped default was unaffected - but the combine back-derives each partial's non-signal variance as sigma^2 - I, and a floored sigma makes that quantity mean nothing. It then read corr^2 * (0.0004 I^2 - I), which is not a background variance. Measured on --integrator boxsum: the reported sigma understated the true scatter by up to 16x at I ~ 21000 counts per partial, and pooled_I amplified a 1 ct/px background drift into an 11.5% intensity error on the strongest reflections. What the floor stood in for - that at high intensity the error is systematic rather than counting - is already carried downstream, twice: the fitted b in v = a*sigma^2 + (b*I)^2, measured from the data rather than assumed, and SigmaWithSystematicFloor on the merged sigma. The floor was that idea applied one level too early with a hardcoded b of 0.02. It arrived without a test or a setter and was unreachable from the CLI, the API and the config. The merge now takes the non-signal variance the integrator actually measured instead of inverting sigma^2 = I + N. That identity is exact for a box sum once the floor is gone and was never exact for a profile fit, whose sigma^2 = 1/den + (wsum/den)^2 * bkg_var is formed against a fitted intensity. The value is carried through BraggFitResult, Reflection and Obs, both engines, both merges, and the process-file round trip; files written before this change are read with the term absent, which is what they had. Battery, 37 crystals, paired: space groups unchanged, reflection sets unchanged, median delta zero on R_meas and CC1/2. --integrator boxsum on the reference crystal goes ISa 8.9 -> 20.2 with a 0.947 -> 1.032. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
88 lines
3.5 KiB
C++
88 lines
3.5 KiB
C++
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "HDF5DataFilePluginReflection.h"
|
|
#include "../include/spdlog/fmt/fmt.h"
|
|
#include "../common/Reflection.h"
|
|
|
|
void HDF5DataFilePluginReflection::OpenFile(HDF5File &data_file, const DataMessage &msg, size_t images_per_file) {
|
|
reflection_group = std::make_unique<HDF5Group>(data_file, "/entry/reflections");
|
|
reflection_group->NXClass("NXcollection");
|
|
}
|
|
|
|
void HDF5DataFilePluginReflection::Write(const DataMessage &msg, uint64_t image_number) {
|
|
if (!reflection_group || msg.reflections.empty())
|
|
return;
|
|
|
|
std::vector<int32_t> h, k, l;
|
|
std::vector<float> I, sigma, d, lp;
|
|
std::vector<float> image, phi, pred_x, pred_y, obs_x, obs_y, bkg, var_bkg, partiality, zeta,
|
|
scale_factor;
|
|
|
|
h.reserve(msg.reflections.size());
|
|
k.reserve(msg.reflections.size());
|
|
l.reserve(msg.reflections.size());
|
|
I.reserve(msg.reflections.size());
|
|
sigma.reserve(msg.reflections.size());
|
|
d.reserve(msg.reflections.size());
|
|
pred_x.reserve(msg.reflections.size());
|
|
pred_y.reserve(msg.reflections.size());
|
|
obs_x.reserve(msg.reflections.size());
|
|
obs_y.reserve(msg.reflections.size());
|
|
bkg.reserve(msg.reflections.size());
|
|
var_bkg.reserve(msg.reflections.size());
|
|
lp.reserve(msg.reflections.size());
|
|
partiality.reserve(msg.reflections.size());
|
|
image.reserve(msg.reflections.size());
|
|
phi.reserve(msg.reflections.size());
|
|
zeta.reserve(msg.reflections.size());
|
|
scale_factor.reserve(msg.reflections.size());
|
|
|
|
for (const auto &refl : msg.reflections) {
|
|
image.emplace_back(refl.image_number);
|
|
h.emplace_back(refl.h);
|
|
k.emplace_back(refl.k);
|
|
l.emplace_back(refl.l);
|
|
I.emplace_back(refl.I);
|
|
sigma.emplace_back(refl.sigma);
|
|
d.emplace_back(refl.d);
|
|
pred_x.emplace_back(refl.predicted_x);
|
|
pred_y.emplace_back(refl.predicted_y);
|
|
obs_x.emplace_back(refl.observed_x);
|
|
obs_y.emplace_back(refl.observed_y);
|
|
bkg.emplace_back(refl.bkg);
|
|
var_bkg.emplace_back(refl.var_bkg);
|
|
lp.emplace_back(1.0/refl.rlp);
|
|
partiality.emplace_back(refl.partiality);
|
|
phi.emplace_back(refl.delta_phi_deg);
|
|
zeta.emplace_back(refl.zeta);
|
|
scale_factor.emplace_back(refl.image_scale_corr);
|
|
}
|
|
|
|
std::string image_group_name = fmt::format("image_{:06d}", image_number);
|
|
HDF5Group image_group(*reflection_group, image_group_name);
|
|
image_group.NXClass("NXreflections");
|
|
image_group.SaveVector("h", h);
|
|
image_group.SaveVector("k", k);
|
|
image_group.SaveVector("l", l);
|
|
image_group.SaveVector("d", d)->Units("Angstrom");
|
|
image_group.SaveVector("delta_phi", phi);
|
|
image_group.SaveVector("predicted_x", pred_x);
|
|
image_group.SaveVector("predicted_y", pred_y);
|
|
image_group.SaveVector("observed_x", obs_x);
|
|
image_group.SaveVector("observed_y", obs_y);
|
|
image_group.SaveVector("int_sum", I);
|
|
image_group.SaveVector("int_err", sigma);
|
|
image_group.SaveVector("background_mean", bkg);
|
|
image_group.SaveVector("background_variance", var_bkg);
|
|
image_group.SaveVector("observed_frame", image);
|
|
image_group.SaveVector("lp", lp);
|
|
image_group.SaveVector("partiality", partiality);
|
|
image_group.SaveVector("zeta", zeta);
|
|
image_group.SaveVector("image_scale_corr", scale_factor);
|
|
}
|
|
|
|
void HDF5DataFilePluginReflection::WriteFinal(HDF5File &data_file) {
|
|
|
|
}
|