From 7c10d62daba739ed6e06507fc795bc07df22b5fd Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 5 Sep 2026 16:28:15 +0200 Subject: [PATCH] integration: the sensor efficiency is carried as its own quantity, not folded into the Lorentz-polarization factor It was multiplied into the per-reflection factor at prediction, so that factor held Lorentz, polarization and efficiency at once and the two spellings that reach a file - the wire key and the reflection dataset - meant something different from what they had meant the day before. The unmerged MTZ had to divide the two apart again at write time to fill its own columns, which is a good sign the wrong thing was being carried. Carry them separately. The prescaling factor is Lorentz and polarization again, what its name and both reference implementations mean by it, and the efficiency is its own field through prediction, integration, serialization and storage. Fifteen sites that want the total now multiply the two - once per reflection, not once per pixel. The efficiency is stored rather than recomputed on read, because the writer has no geometry to recompute it from, and because a file written before the correction existed would have had a radial trend invented for it. Sixty stored files were checked for the one combination that would be ambiguous - the old meaning of the factor beside a stored efficiency - and none carries it. Output does not move. Re-scaling a file written before the efficiency existed is byte-identical, which is a proof rather than a sample, since the stored factor is exactly one there. Where the efficiency is live, one product is reassociated - (L*Q)/P becomes (L/P)*Q - and about a third of the values differ in the last bit or two: every structural column is identical, so no reflection is gained, lost or reindexed, and no intensity in 1.4 million observations moves by as much as 1e-4 of its own sigma. The parity tests now compare the efficiency as well, and their non-vacuity guard watches it rather than the factor it left - which is the same guard that went blind when the efficiency was added to a field it was not watching. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- common/Reflection.h | 27 +++++++------- docs/CBOR.md | 4 +- docs/CHANGELOG.md | 2 +- docs/CPU_DATA_ANALYSIS_INTEGRATION.md | 2 +- docs/HDF5.md | 4 +- image_analysis/UpdateReflectionResolution.cpp | 2 +- image_analysis/WriteReflections.cpp | 37 +++++++++++-------- .../BraggIntegrationEngine.cpp | 2 +- .../bragg_prediction/BraggPrediction.cpp | 2 +- .../bragg_prediction/BraggPredictionGPU.cu | 2 +- .../bragg_prediction/BraggPredictionRot.cpp | 2 +- .../bragg_prediction/BraggPredictionRotGPU.cu | 2 +- .../scale_merge/AnisotropyAnalysis.cpp | 10 +++-- image_analysis/scale_merge/HKLKey.cpp | 6 ++- .../scale_merge/ReindexAmbiguity.cpp | 2 +- .../scale_merge/RotationScaleMerge.cpp | 9 +++-- image_analysis/scale_merge/ScaleOnTheFly.cpp | 7 ++-- .../scale_merge/StillsPartialityRefine.cpp | 11 +++--- tests/CalcBraggPredictionTest.cpp | 20 ++++++---- tests/WriteReflectionsTest.cpp | 12 +++--- 20 files changed, 90 insertions(+), 75 deletions(-) diff --git a/common/Reflection.h b/common/Reflection.h index e1bb0c772..1f78c921f 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -25,24 +25,23 @@ struct Reflection { float var_bkg; // non-signal (background) part of sigma^2, carried to the merge float sigma; float dist_ewald; - // Everything known a priori that multiplies the raw integrated count on its way to a quantity - // proportional to |F|^2: the reciprocal Lorentz factor, the reciprocal polarization factor and - // the sensor's angle-dependent efficiency, whichever of them the path that filled it applies. - // It is deliberately NOT an enumeration of those terms - it is the product of every deterministic - // per-reflection correction, and that set has grown before and will grow again. What it is not is - // a scale: the fitted per-image scale and the partiality stay out of it and are divided in - // separately below. (Named after DIALS's prescaling_correction, which holds the same product.) + // The reciprocal Lorentz factor (rotation only - a still's Lorentz factor is one) times the + // reciprocal polarization factor, and nothing else. This is what LP means everywhere in the + // field, and it is what the CBOR key "rlp" and the HDF5 dataset "lp" store the reciprocal of. + // What it is not is a scale: the fitted per-image scale and the partiality stay out of it and + // are divided in separately below. (Named after DIALS's prescaling_correction.) float prescaling_corr; - // The sensor-efficiency part of prescaling_corr on its own: QE(0)/QE(alpha), always <= 1, and 1 - // where the sensor is opaque or its material and thickness are unknown. The one term singled out - // of the product above, and only so the written files can say which part of it is - // Lorentz-polarization and which is detector response: LP means L/P everywhere else in the - // field, and prescaling_corr / qe_corr is that. Defaulted to 1 so a reflection read from a file - // written before this existed is a no-op rather than a zero. + // The sensor's angle-dependent efficiency, QE(0)/QE(alpha): always <= 1, and exactly 1 where the + // sensor is opaque or its material and thickness are unknown. It is carried BESIDE + // prescaling_corr rather than inside it, because LP and detector response are two different + // things and every file the field reads keeps them apart. The total deterministic correction on + // a reflection is therefore prescaling_corr * qe_corr, and every site that corrects an intensity + // multiplies the two. Defaulted to 1 so a reflection read from a file written before this + // existed is a no-op rather than a zero. float qe_corr = 1.0f; float partiality; // fraction of the reflection recorded in the sampled (rocking) slice float zeta; - float image_scale_corr; // I_true = image_scale_corr * I; = prescaling_corr / (partiality * image_scale) + float image_scale_corr; // I_true = image_scale_corr * I; = prescaling_corr * qe_corr / (partiality * image_scale) bool observed = false; bool on_ice_ring = false; // sits on a hexagonal-ice powder ring: excluded from scaling, kept for merging }; diff --git a/docs/CBOR.md b/docs/CBOR.md index 522f1542d..d827d5aa3 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -167,8 +167,8 @@ See [DECTRIS documentation](https://github.com/dectris/documentation/tree/main/s | - sigma | float | standard deviation, estimated from counting statistics (photons) | | | | - image | float | image number (present for each spot) | | | | - rp | float | Distance to Ewald sphere \[Angstrom^-1\] | | | -| - rlp | float | Prescaling correction: the product of every deterministic per-reflection correction (reciprocal Lorentz, reciprocal polarization, sensor efficiency). The key name is historical and kept for compatibility | | | -| - qe | float | The sensor-efficiency part of `rlp` on its own, QE(0)/QE(alpha); 1 where the sensor is opaque or unknown. Optional | | | +| - rlp | float | Reciprocal Lorentz-polarization factor: the multiplier taking the raw integrated count toward a quantity proportional to \|F\|^2. Lorentz x polarization only - a still has no Lorentz term, so there it is the polarization alone | | | +| - qe | float | Sensor efficiency at the reflection's angle of incidence, QE(0)/QE(alpha); <= 1, and 1 where the sensor is opaque or unknown. Carried beside `rlp`, not inside it: the total correction is `rlp * qe`. Optional | | | | - partiality | float | Partiality of the reflection | | | | - phi | float | phi angle from XDS: difference from middle of current frame, not absolute \[deg\] | | | | - zeta | float | Lorentz zeta factor (reciprocal-space geometry term) | | | diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 80b284329..fe9ba6b5c 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -13,7 +13,7 @@ * `rugnux` corrects each reflection for the sensor's quantum efficiency at the angle the diffracted beam meets the detector. * `rugnux` applies the sensor quantum-efficiency correction to still images processed on a GPU, which were left uncorrected while the same data on a CPU-only machine was corrected. -* The unmerged MTZ keeps the sensor efficiency out of its `LP` column and records it in a new `QE` column, as XDS and DIALS do, so raw counts are `I / LP * QE`; the merged intensities are unchanged, and `_process.h5` gains a matching optional `qe` dataset. +* The unmerged MTZ keeps the sensor efficiency out of its `LP` column and records it in a new `QE` column, as XDS and DIALS do, so raw counts are `I / LP * QE`; `_process.h5` keeps it out of `lp` in the same way, in a new optional `qe` dataset, and the merged intensities are unchanged. * `rugnux` takes the sensor attenuation length from tabulated coefficients rather than a wavelength-cubed approximation, which corrects the parallax term of the spot-width variance on CdTe sensors. * `rugnux --model` treats the model as a hypothesis: it decides the enantiomorph and the indexing only where its R-work beats that of the same model in random orientations, and a model the data reject is still scored, placed and mapped, but leaves the reflection files byte for byte what a run with no model writes. * `rugnux --model` places the model against the data as a rigid body before scoring it, writes sigma_A-weighted 2mFo-DFc and mFo-DFc maps in place of the unweighted 2Fo-Fc and Fo-Fc, and writes the model as it was placed - `_model.cif`, and `_model.pdb` where the PDB format can express the cell - in the cell and space group of the reflection files beside it. diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 15166f456..4820b9220 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -161,7 +161,7 @@ The integrator is selected by `--integrator boxsum|gaussian|empirical` (default ### 9.4 The prescaling correction -Every deterministic per-reflection correction is carried as one multiplicative factor, `prescaling_corr`, which the prediction fills and the integrator adds to. It currently holds the reciprocal Lorentz factor (rotation only), the reciprocal polarization factor from the geometry-based term (§2.2), and the sensor's angle-dependent efficiency (§9.6); it is a product, not a fixed list, and terms have been added to it before. Scaling and merging (§10) use it as the numerator of the reflection's total correction. It is not a scale — the fitted per-image scale and the partiality are separate. The unmerged MTZ splits it back out into an `LP` column (Lorentz x polarization) and a `QE` column (the sensor efficiency, as a divisor normalised to 1 at normal incidence), so that `LP` means what XDS and DIALS mean by it and raw counts are `I / LP * QE`. +The deterministic per-reflection corrections are carried as two multiplicative factors. `prescaling_corr` holds the reciprocal Lorentz factor (rotation only — a still's Lorentz factor is one) times the reciprocal polarization factor from the geometry-based term (§2.2), and nothing else: it is Lorentz x polarization, which is what `LP` means in every format the field reads. `qe_corr` holds the sensor's angle-dependent efficiency (§9.6) beside it, because detector response and beam/crystal geometry are different things. The total deterministic correction on a reflection is the product `prescaling_corr * qe_corr`, and every site that corrects an intensity — the integrator, the scaling fits, the merge ingest, the anisotropy analysis and the unmerged export — multiplies the two. Neither is a scale: the fitted per-image scale and the partiality are separate. The two reach the unmerged MTZ as its `LP` and `QE` columns unchanged (`QE` as a divisor normalised to 1 at normal incidence), so raw counts are `I / LP * QE`. Two terms are deliberately **not** in it. The **per-pixel solid angle**, which the azimuthal profile does divide out (§2.2), is correctly absent here: a Bragg integration sums all the photons in a diff --git a/docs/HDF5.md b/docs/HDF5.md index d3345709f..63694963b 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -435,8 +435,8 @@ mostly onto the standard | `predicted_x`, `predicted_y` | pixel | name standard, units differ | predicted position. NXreflections `predicted_x/_y` are *physical* lengths; the pixel datasets are `predicted_px_x/_y` | | `observed_x`, `observed_y` | pixel | name standard, units differ | observed centroid (pixels; standard pixel form is `observed_px_x/_y`) | | `observed_frame` | | standard | image number of the reflection | -| `lp` | | standard | reciprocal of the prescaling correction (stored as `1/prescaling_corr`). The NXreflections name is `lp`, but the stored product also carries the sensor efficiency term | -| `qe` | | **extension** | the sensor-efficiency part of `lp` on its own, in the same reciprocal convention, so `lp / qe` is the pure Lorentz–polarization reciprocal. Absent in files written before it existed; the reader then takes the whole of `lp` for Lorentz–polarization, which is what those files mean | +| `lp` | | standard | the Lorentz-polarization factor, stored as the reciprocal of the multiplier that is applied. Lorentz x polarization only, which is what the NXreflections name means; the sensor efficiency is `qe`, beside it | +| `qe` | | **extension** | the sensor efficiency, in the same reciprocal convention as `lp`, so the whole correction applied to a raw count is `1/lp * 1/qe`. Absent in files written before it existed; the reader then takes the whole of `lp` for Lorentz–polarization, which is what those files mean | | `partiality` | | standard | recorded fraction of the reflection | | `delta_phi` | deg | **extension** | XDS Δφ: offset from the centre of the current frame | | `zeta` | | **extension** | Lorentz ζ factor (reciprocal-space geometry term) | diff --git a/image_analysis/UpdateReflectionResolution.cpp b/image_analysis/UpdateReflectionResolution.cpp index 3f3bbddbb..9a2e4af84 100644 --- a/image_analysis/UpdateReflectionResolution.cpp +++ b/image_analysis/UpdateReflectionResolution.cpp @@ -30,7 +30,7 @@ ResolutionStats UpdateReflectionResolution(const UnitCell &cell, std::vector SumRockingEvents(const std::vector & for (size_t m = i; m < j; ++m) { const Reflection &r = *parts[m].r; const double p = r.partiality; + const float pc = r.prescaling_corr * r.qe_corr; // the whole correction, LP x QE sum_p += p; - sum_I += static_cast(r.I) * r.prescaling_corr; - sum_var += static_cast(r.sigma) * r.sigma * r.prescaling_corr * r.prescaling_corr; - sum_var_bkg += static_cast(r.var_bkg) * r.prescaling_corr * r.prescaling_corr; - p_corr += p * r.prescaling_corr; + sum_I += static_cast(r.I) * pc; + sum_var += static_cast(r.sigma) * r.sigma * pc * pc; + sum_var_bkg += static_cast(r.var_bkg) * pc * pc; + p_corr += p * pc; p_qe += p * r.qe_corr; p_frame += p * r.image_number; p_x += p * DetectorX(r); @@ -689,14 +690,17 @@ std::vector SumRockingEvents(const std::vector & if (sum_p < min_partiality || sum_p < min_captured_fraction) continue; - // The prescaling factor is applied by the writer, which multiplies I by prescaling_corr, so divide the - // event's own factor back out of the sums here. It is the same geometry for every part of - // one event to a median 2e-4, so the file's I / LP * QE is still the raw count sum. - full.prescaling_corr = static_cast(p_corr / sum_p); + // The prescaling factor is applied by the writer, which multiplies I by the whole correction, so + // divide the event's own factor back out of the sums here. It is the same geometry for every part + // of one event to a median 2e-4, so the file's I / LP * QE is still the raw count sum. The event's + // mean is taken on the whole correction and on its QE part; the Lorentz-polarization half is what + // is left when the second is divided out of the first. + const float mean_corr = static_cast(p_corr / sum_p); full.qe_corr = static_cast(p_qe / sum_p); - full.I = static_cast(sum_I / full.prescaling_corr); - full.sigma = static_cast(std::sqrt(sum_var) / full.prescaling_corr); - full.var_bkg = static_cast(sum_var_bkg / (static_cast(full.prescaling_corr) * full.prescaling_corr)); + full.prescaling_corr = mean_corr / full.qe_corr; + full.I = static_cast(sum_I / mean_corr); + full.sigma = static_cast(std::sqrt(sum_var) / mean_corr); + full.var_bkg = static_cast(sum_var_bkg / (static_cast(mean_corr) * mean_corr)); full.partiality = static_cast(sum_p); full.image_number = static_cast(p_frame / sum_p); full.observed_x = static_cast(p_x / sum_p); @@ -819,15 +823,16 @@ void WriteUnmergedMtzReflections(const std::vector &outcomes mtz.data.push_back(static_cast(hkl[2])); mtz.data.push_back(static_cast((partials ? 256 : 0) + isym)); mtz.data.push_back(static_cast(batch)); - mtz.data.push_back(r.I * r.prescaling_corr); - mtz.data.push_back(r.sigma * r.prescaling_corr); + const float corr = r.prescaling_corr * r.qe_corr; + mtz.data.push_back(r.I * corr); + mtz.data.push_back(r.sigma * corr); mtz.data.push_back(r.partiality); mtz.data.push_back(DetectorX(r)); mtz.data.push_back(DetectorY(r)); mtz.data.push_back(phi_start_deg(r.image_number) + wedge_deg / 2.0f); - // prescaling_corr is the whole prescaling product; qe_corr is its sensor-efficiency part. LP is what - // is left when that part is divided out, and QE is that part in the divide-by direction. - mtz.data.push_back(r.prescaling_corr / r.qe_corr); + // prescaling_corr is Lorentz x polarization, which is what LP means; qe_corr is the sensor + // efficiency beside it, written in the divide-by direction. + mtz.data.push_back(r.prescaling_corr); mtz.data.push_back(1.0f / r.qe_corr); mtz.data.push_back(0.0f); // FLAG: nothing here is a rejected observation mtz.data.push_back(r.delta_phi_deg); diff --git a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp index b1dad1171..f24b6e5b5 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngine.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngine.cpp @@ -233,7 +233,7 @@ std::vector BraggIntegrationEngine::Finalize(const std::vector(image_number); out.push_back(refl); } diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index 2d29b573a..be62c953c 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -210,7 +210,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal .observed_y = NAN, .d = d, .dist_ewald = dist_ewald_sphere, - .prescaling_corr = qe_corr, + .prescaling_corr = 1.0f, .qe_corr = qe_corr, .partiality = 1.0f, .zeta = 1.0, diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 28eb59ed7..9a87749d9 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -142,7 +142,7 @@ namespace { out.observed_y = NAN; out.d = 1.0f / sqrtf(recip_sq); out.dist_ewald = dist_ewald; - out.prescaling_corr = qe_corr; + out.prescaling_corr = 1.0f; out.qe_corr = qe_corr; out.partiality = 1.0f; out.zeta = 1.0f; diff --git a/image_analysis/bragg_prediction/BraggPredictionRot.cpp b/image_analysis/bragg_prediction/BraggPredictionRot.cpp index 9b9643da3..525784466 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRot.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRot.cpp @@ -195,7 +195,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys .observed_y = NAN, .d = d, .dist_ewald = dist_ewald_sphere, - .prescaling_corr = lorentz_reciprocal * qe_corr, + .prescaling_corr = lorentz_reciprocal, .qe_corr = qe_corr, .partiality = partiality, .zeta = zeta_abs, diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index 513059635..2994a38eb 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -190,7 +190,7 @@ namespace { out[count].observed_y = NAN; out[count].d = 1.0f / sqrtf(p0_sq); out[count].dist_ewald = dist_ewald; - out[count].prescaling_corr = lorentz * qe_corr; + out[count].prescaling_corr = lorentz; out[count].qe_corr = qe_corr; out[count].partiality = partiality; out[count].zeta = zeta_abs; diff --git a/image_analysis/scale_merge/AnisotropyAnalysis.cpp b/image_analysis/scale_merge/AnisotropyAnalysis.cpp index 8b3f437dc..20e9e7a6e 100644 --- a/image_analysis/scale_merge/AnisotropyAnalysis.cpp +++ b/image_analysis/scale_merge/AnisotropyAnalysis.cpp @@ -1042,8 +1042,9 @@ std::vector ScaledObservations(const std::vector 0.0f && r.d > 0.0f - && std::isfinite(r.prescaling_corr) && r.prescaling_corr > 0.0f && std::isfinite(r.partiality); + && std::isfinite(pc) && pc > 0.0f && std::isfinite(r.partiality); }; if (!rotation) { @@ -1055,7 +1056,7 @@ std::vector ScaledObservations(const std::vector= min_partiality)) continue; - const double corr = r.prescaling_corr / (r.partiality * g[i]); + const double corr = (r.prescaling_corr * r.qe_corr) / (r.partiality * g[i]); out.push_back({r.h, r.k, r.l, static_cast(r.I * corr), static_cast(r.sigma * corr), r.d}); } @@ -1132,8 +1133,9 @@ std::vector ScaledObservations(const std::vector(r.I) * r.prescaling_corr; - sum_var += static_cast(r.sigma) * r.sigma * r.prescaling_corr * r.prescaling_corr; + const float pc = r.prescaling_corr * r.qe_corr; + sum_I += static_cast(r.I) * pc; + sum_var += static_cast(r.sigma) * r.sigma * pc * pc; p_g += r.partiality * g[parts[m].image]; } const Reflection &first = *parts[i].r; diff --git a/image_analysis/scale_merge/HKLKey.cpp b/image_analysis/scale_merge/HKLKey.cpp index 3dc4b4d2d..0ca60d009 100644 --- a/image_analysis/scale_merge/HKLKey.cpp +++ b/image_analysis/scale_merge/HKLKey.cpp @@ -81,7 +81,8 @@ bool AcceptReflection(const Reflection &r, std::optional d_min_limit, st return false; if (d_max_limit && r.d > d_max_limit) return false; - if (!std::isfinite(r.prescaling_corr) || r.prescaling_corr == 0.0f) + const float corr = r.prescaling_corr * r.qe_corr; + if (!std::isfinite(corr) || corr == 0.0f) return false; if (!std::isfinite(r.sigma) || r.sigma <= 0.0) return false; @@ -97,7 +98,8 @@ bool AcceptReflection(const Reflection &r, double d_min_limit, double d_max_limi return false; if (d_max_limit > 0.0 && r.d > d_max_limit) return false; - if (!std::isfinite(r.prescaling_corr) || r.prescaling_corr == 0.0f) + const float corr = r.prescaling_corr * r.qe_corr; + if (!std::isfinite(corr) || corr == 0.0f) return false; if (!std::isfinite(r.sigma) || r.sigma <= 0.0) return false; diff --git a/image_analysis/scale_merge/ReindexAmbiguity.cpp b/image_analysis/scale_merge/ReindexAmbiguity.cpp index a58b6962f..f28b60dfa 100644 --- a/image_analysis/scale_merge/ReindexAmbiguity.cpp +++ b/image_analysis/scale_merge/ReindexAmbiguity.cpp @@ -173,7 +173,7 @@ double ReindexAmbiguityResolver::ReferenceCC(const std::vector &refl const auto it = reference_data.find(hkl_key_generator(h[0], h[1], h[2])); if (it == reference_data.end()) continue; - const double x = static_cast(r.I) * r.prescaling_corr / r.partiality; + const double x = static_cast(r.I) * (r.prescaling_corr * r.qe_corr) / r.partiality; const double y = it->second; if (!std::isfinite(x) || !std::isfinite(y)) continue; diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 09cefa06c..39610d3b1 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -608,10 +608,11 @@ void RotationScaleMerge::BuildInRangeObservations(const std::vector &ke Obs &obs = partials[at]; obs.h = r.h; obs.k = r.k; obs.l = r.l; obs.I = r.I; obs.sigma = r.sigma; obs.d = r.d; - // Divide the incident flux out of prescaling_corr here: it is the numerator of corr, so - // the fits below see only what is left. The flux is the one term in the product that is - // measured from the data rather than computed from the geometry. - obs.prescaling_corr = static_cast(r.prescaling_corr / frame_flux[o]); + // The whole deterministic correction - Lorentz x polarization x sensor efficiency - is the + // numerator of corr, so the two halves are multiplied back together here. The incident flux is + // divided out at the same time so the fits below see only what is left; it is the one term in + // the product that is measured from the data rather than computed from the geometry. + obs.prescaling_corr = static_cast(r.prescaling_corr * r.qe_corr / frame_flux[o]); obs.partiality = r.partiality; obs.zeta = r.zeta; obs.delta_phi = r.delta_phi_deg; obs.bkg = r.bkg; obs.var_bkg = r.var_bkg; obs.px = r.predicted_x; obs.py = r.predicted_y; obs.image_number = r.image_number; diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index 8da6ca1bf..d6a789616 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -114,7 +114,7 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { const auto it = reference_data.find(hkl_key_generator(r)); if (it == reference_data.end()) continue; - const double coeff = r.partiality * SafeInv(r.prescaling_corr, 1.0) * it->second; + const double coeff = r.partiality * SafeInv(r.prescaling_corr * r.qe_corr, 1.0) * it->second; obs.push_back({coeff, static_cast(r.I), SafeInv(r.sigma, 1.0)}); } @@ -127,9 +127,10 @@ void ScaleOnTheFly::Scale(IntegrationOutcome &integration_outcome) const { } for (auto &r: integration_outcome.reflections) { + const float corr = r.prescaling_corr * r.qe_corr; const double denom = r.partiality * result.G; - r.image_scale_corr = (std::isfinite(r.prescaling_corr) && std::isfinite(denom) && denom > 0.0) - ? static_cast(r.prescaling_corr / denom) + r.image_scale_corr = (std::isfinite(corr) && std::isfinite(denom) && denom > 0.0) + ? static_cast(corr / denom) : NAN; } diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index bef2884fd..2ae8ee39d 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -32,7 +32,7 @@ namespace { double qx, qy, qz; double Iref; double Iobs; - double inv_corr; // 1 / prescaling_corr + double inv_corr; // 1 / (prescaling_corr * qe_corr) double weight; // 1 / sigma }; @@ -186,7 +186,7 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, .qx = q.x, .qy = q.y, .qz = q.z, .Iref = it->second, .Iobs = static_cast(r.I), - .inv_corr = SafeInv(r.prescaling_corr, 1.0), + .inv_corr = SafeInv(r.prescaling_corr * r.qe_corr, 1.0), .weight = SafeInv(r.sigma, 1.0), }); @@ -306,16 +306,17 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, before.emplace_back(r.partiality, r.image_scale_corr); // Write the refined partiality + scale correction onto every reflection of the crystal (not only the - // fit subset), so the merge sees a consistent model. image_scale_corr = prescaling_corr / + // fit subset), so the merge sees a consistent model. image_scale_corr = prescaling_corr * qe_corr / // (partiality * G), the same composition ScaleOnTheFly writes. for (auto &r: outcome.reflections) { const Coord q = base_q(r); const double p = ComputeP(q.x, q.y, q.z, psi[0], psi[1], S0.x, S0.y, S0.z, inv_lambda, gamma0, gamma_e, bw); r.partiality = static_cast(p); + const float corr = r.prescaling_corr * r.qe_corr; const double denom = p * G; - r.image_scale_corr = (std::isfinite(r.prescaling_corr) && std::isfinite(denom) && denom > 0.0) - ? static_cast(r.prescaling_corr / denom) + r.image_scale_corr = (std::isfinite(corr) && std::isfinite(denom) && denom > 0.0) + ? static_cast(corr / denom) : NAN; } outcome.image_scale_g = static_cast(G); diff --git a/tests/CalcBraggPredictionTest.cpp b/tests/CalcBraggPredictionTest.cpp index a1f2df524..8bc490ba5 100644 --- a/tests/CalcBraggPredictionTest.cpp +++ b/tests/CalcBraggPredictionTest.cpp @@ -412,12 +412,13 @@ TEST_CASE("BraggPrediction_CPU_GPU_consistency_tilted") { CHECK(cpu_r.predicted_x == Catch::Approx(gpu_r.predicted_x).margin(0.1)); CHECK(cpu_r.predicted_y == Catch::Approx(gpu_r.predicted_y).margin(0.1)); CHECK(cpu_r.d == Catch::Approx(gpu_r.d).margin(0.01)); - // The prescaling factor is part of the prediction, not a downstream product: the - // sensor-efficiency term lived on the CPU path alone for a while because nothing here - // compared it. + // Both halves of the correction are part of the prediction, not a downstream product: + // the sensor-efficiency term lived on the CPU path alone for a while because nothing here + // compared it, and it is now a field of its own - so it is compared as a field of its own. CHECK(cpu_r.prescaling_corr == Catch::Approx(gpu_r.prescaling_corr).epsilon(1e-4)); + CHECK(cpu_r.qe_corr == Catch::Approx(gpu_r.qe_corr).epsilon(1e-4)); CHECK(cpu_r.image_scale_corr == Catch::Approx(gpu_r.image_scale_corr).epsilon(1e-4)); - min_corr = std::min(min_corr, cpu_r.prescaling_corr); + min_corr = std::min(min_corr, cpu_r.qe_corr); matched++; } } @@ -425,14 +426,16 @@ TEST_CASE("BraggPrediction_CPU_GPU_consistency_tilted") { // Most reflections should match (allow for some numerical differences at boundaries) CHECK(matched > cpu_count * 0.95); // ... and the comparison above must not be vacuous: on this geometry (320 um Si at 12 keV, - // reflections out to 2 A) the sensor correction is several per cent, so a factor that stayed - // at 1 on both sides would mean the correction had been dropped from BOTH paths. + // reflections out to 2 A) the sensor correction is several per cent, so a qe_corr that stayed + // at 1 on both sides would mean the correction had been dropped from BOTH paths. It is checked + // on qe_corr and not on prescaling_corr, which no longer carries it. CHECK(min_corr < 0.99f); } TEST_CASE("BraggPredictionRot_CPU_GPU_consistency_tilted") { // The rotation counterpart of the test above. Same purpose: every per-reflection quantity the - // two implementations both produce has to agree, the prescaling factor included. + // two implementations both produce has to agree - the Lorentz-polarization factor and the sensor + // efficiency each in their own field, so neither can hide behind the other in a product. DiffractionExperiment experiment(DetJF4M()); experiment.DetectorDistance_mm(100.0).BeamX_pxl(1500.0).BeamY_pxl(1000.0) .PoniRot1_rad(0.04).PoniRot2_rad(-0.025) @@ -473,7 +476,8 @@ TEST_CASE("BraggPredictionRot_CPU_GPU_consistency_tilted") { CHECK(cpu_r.predicted_y == Catch::Approx(gpu_r.predicted_y).margin(0.1)); CHECK(cpu_r.d == Catch::Approx(gpu_r.d).margin(0.01)); CHECK(cpu_r.prescaling_corr == Catch::Approx(gpu_r.prescaling_corr).epsilon(1e-3)); - min_corr = std::min(min_corr, cpu_r.prescaling_corr); + CHECK(cpu_r.qe_corr == Catch::Approx(gpu_r.qe_corr).epsilon(1e-3)); + min_corr = std::min(min_corr, cpu_r.qe_corr); matched++; } } diff --git a/tests/WriteReflectionsTest.cpp b/tests/WriteReflectionsTest.cpp index e458f836b..845354287 100644 --- a/tests/WriteReflectionsTest.cpp +++ b/tests/WriteReflectionsTest.cpp @@ -82,7 +82,7 @@ TEST_CASE("Unmerged MTZ: LP is Lorentz-polarization and QE carries the sensor ef IntegrationOutcome outcome; const float raw[3] = {1000.0f, 250.0f, 40.0f}; - const float corr[3] = {1.75f, 2.50f, 0.90f}; // the whole prescaling product + const float lp[3] = {1.75f, 2.50f, 0.90f}; // Lorentz x polarization, and nothing else const float qe[3] = {0.9375f, 0.8125f, 1.0f}; // 1.0 = the sensor said nothing to correct for (int i = 0; i < 3; ++i) { Reflection r{}; @@ -91,7 +91,7 @@ TEST_CASE("Unmerged MTZ: LP is Lorentz-polarization and QE carries the sensor ef r.d = 5.0f + i; r.I = raw[i]; // the writer is what applies the factor r.sigma = std::sqrt(raw[i]); - r.prescaling_corr = corr[i]; + r.prescaling_corr = lp[i]; r.qe_corr = qe[i]; r.partiality = 1.0f; r.predicted_x = 100.0f + i; r.predicted_y = 200.0f + i; @@ -117,14 +117,14 @@ TEST_CASE("Unmerged MTZ: LP is Lorentz-polarization and QE carries the sensor ef const float LP = mtz.data[i * mtz.columns.size() + c_lp->idx]; const float QE = mtz.data[i * mtz.columns.size() + c_qe->idx]; INFO("row " << i); - // LP holds Lorentz x polarization alone: the sensor term has been divided out of it. - CHECK(LP == Catch::Approx(corr[i] / qe[i]).epsilon(1e-5)); + // LP holds Lorentz x polarization alone: the sensor term was never inside it. + CHECK(LP == Catch::Approx(lp[i]).epsilon(1e-5)); // QE is a divisor normalised to 1 at normal incidence, so it never drops below 1. CHECK(QE == Catch::Approx(1.0f / qe[i]).epsilon(1e-5)); CHECK(QE >= 1.0f); // ... and the two together put the raw counts back. CHECK(I / LP * QE == Catch::Approx(raw[i]).epsilon(1e-4)); - // The intensity itself did not move: it is still the fully corrected value. - CHECK(I == Catch::Approx(raw[i] * corr[i]).epsilon(1e-5)); + // The intensity itself is the fully corrected value - both halves applied. + CHECK(I == Catch::Approx(raw[i] * lp[i] * qe[i]).epsilon(1e-5)); } }