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)); } }