diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 89d5413db..cb69d7497 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -25,6 +26,9 @@ namespace { // Model error of the per-crystal scale fit, as a fraction of the predicted intensity - ScaleOnTheFly's. constexpr double SCALE_MODEL_ERROR = 0.3; + // The width scales BalanceRunWidth searches: from 1/BALANCE_RANGE to BALANCE_RANGE times the table. + constexpr double BALANCE_RANGE = 3.0; + // The smallest priors on a crystal's radial scale and log width scale (RunPrior). constexpr double MIN_PRIOR_ETA = 1e-5; constexpr double MIN_PRIOR_LN_WIDTH = 0.01; @@ -743,6 +747,85 @@ std::optional StillsPartialityRefine::MeasureR return width; } +StillsPartialityRefine::RunWidth StillsPartialityRefine::BalanceRunWidth( + const std::vector &outcomes, const std::map &reference, + const std::optional &reference_cell, const RunWidth &width) const { + // The width table is measured from moments of the excitation error through the prediction window, and + // the tails of the curve are where that is least sure: on rotation sweeps processed as stills, the + // reflections recorded at under a third of their peak came out 30% fainter than the model at low + // resolution and up to 50% brighter at high resolution, against those recorded near their peak. So + // each shell of the table is scaled by the factor s at which its reflections predicted below half + // their peak and those above it put the same scale on the merged intensities: + // sum I / sum p_s I_ref is the same for both, p_s the Gaussian of width s sigma. + if (width.table_dstar.empty()) + return width; + const auto &indexing = experiment_.GetIndexingSettings(); + struct Point { double u, x, ref; }; // eps / sigma, I / G (corrected), the reference intensity + std::vector> shells(width.table_dstar.size()); + const Width w{.gamma0 = width.gamma0, .gamma_e = width.gamma_e, .bw = width.bandwidth_sigma, + .table_dstar = &width.table_dstar, .table_s = &width.table_s}; + for (const auto &o: outcomes) { + if (!o.image_scale_g || !(*o.image_scale_g > 0.0)) + continue; + if (reference_cell && !o.latt.GetUnitCell().is_close(*reference_cell, indexing.GetUnitCellDistTolerance(), + indexing.GetUnitCellAngleTolerance_deg())) + continue; + const Coord S0 = o.geom.GetScatteringVector(); + const double inv_lambda = 1.0 / o.geom.GetWavelength_A(); + for (const auto &r: o.reflections) { + if (r.on_ice_ring || !AcceptReflection(r, d_min_limit_, d_max_limit_) || !std::isfinite(r.I)) + continue; + const auto it = reference.find(hkl_key_generator_(r)); + if (it == reference.end() || !(it->second > 0.0)) + continue; + const Coord q = (o.latt.Astar() * static_cast(r.h) + o.latt.Bstar() * static_cast(r.k) + + o.latt.Cstar() * static_cast(r.l)) * static_cast(width.radial_scale); + const double dstar = VecLen(q.x, q.y, q.z); + const double de = VecLen(q.x + S0.x, q.y + S0.y, q.z + S0.z) - inv_lambda; + const double x = r.I * r.prescaling_corr * r.qe_corr * r.flight_corr / *o.image_scale_g; + if (!std::isfinite(x)) + continue; + // the table shell nearest in d* + size_t k = std::upper_bound(width.table_dstar.begin(), width.table_dstar.end(), dstar) + - width.table_dstar.begin(); + if (k == width.table_dstar.size() || (k > 0 && dstar - width.table_dstar[k - 1] < width.table_dstar[k] - dstar)) + k = k > 0 ? k - 1 : 0; + shells[k].push_back({de / std::sqrt(w.Sigma2(dstar, q.z)), x, it->second}); + } + } + // log(sum I / sum p_s I_ref) over the reflections below half the table's peak, less the same over those + // above it; falls as s grows. + auto gap = [](const std::vector &pts, double s) { + const double half = std::sqrt(2.0 * std::log(2.0)); // |u| at half the peak + double x_lo = 0.0, m_lo = 0.0, x_hi = 0.0, m_hi = 0.0; + for (const auto &p: pts) { + const double ps = std::exp(-0.5 * p.u * p.u / (s * s)); + (std::fabs(p.u) > half ? x_lo : x_hi) += p.x; + (std::fabs(p.u) > half ? m_lo : m_hi) += ps * p.ref; + } + if (!(x_lo > 0.0 && m_lo > 0.0 && x_hi > 0.0 && m_hi > 0.0)) + return std::numeric_limits::quiet_NaN(); + return std::log(x_lo / m_lo) - std::log(x_hi / m_hi); + }; + RunWidth balanced = width; + std::string table; + for (size_t k = 0; k < shells.size(); ++k) { + double lo = 1.0 / BALANCE_RANGE, hi = BALANCE_RANGE; + const double g_lo = gap(shells[k], lo), g_hi = gap(shells[k], hi); + if (shells[k].size() < 400 || !(g_lo > 0.0) || !(g_hi < 0.0)) + continue; + for (int i = 0; i < 40; ++i) { + const double mid = std::sqrt(lo * hi); + (gap(shells[k], mid) > 0.0 ? lo : hi) = mid; + } + balanced.table_s[k] *= std::sqrt(lo * hi); + table += fmt::format(" {:.2f} A x{:.2f}", 1.0 / width.table_dstar[k], std::sqrt(lo * hi)); + } + Logger("StillsPartialityRefine").Info("Stills rocking width balanced on the partials, against resolution:{}", + table.empty() ? " unchanged" : table); + return balanced; +} + StillsPartialityRefine::CrystalPrior StillsPartialityRefine::RunPrior(const std::vector &fits) const { // The spread of a parameter over the crystals: the variance of the estimates less the median variance // they carry from their own data, floored at a small value (a crystal is never held exactly). @@ -805,8 +888,6 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si run_width = MeasureRunWidth(outcomes, reference_cell, radial_scale); } const double radial_scale = run_width ? run_width->radial_scale : 1.0; - if (used_width) - *used_width = run_width; // The priors on a crystal's radial scale and width are the run's own spread of them: the first cycle // refines every crystal against loose priors, and the spread of those estimates less their own @@ -821,6 +902,8 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si std::map reference; for (const auto &m: merged) reference[hkl_key_generator_(m)] = m.I; + if (outer == 0 && run_width) + run_width = BalanceRunWidth(outcomes, reference, reference_cell, *run_width); std::atomic next{0}; auto worker = [&]() { @@ -848,6 +931,9 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si prior = RunPrior(fits); } + if (used_width) + *used_width = run_width; + double tilt_sum = 0.0; size_t tilt_n = 0; for (const auto &f: fits) diff --git a/image_analysis/scale_merge/StillsPartialityRefine.h b/image_analysis/scale_merge/StillsPartialityRefine.h index 6f644fb8c..48c59c76a 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.h +++ b/image_analysis/scale_merge/StillsPartialityRefine.h @@ -109,6 +109,9 @@ private: CrystalFit RefineOne(IntegrationOutcome &outcome, const std::map &reference, const std::optional &run_width, double radial_scale, const CrystalPrior &prior) const; + // The run's width table, each shell rescaled so that its weak and strong partials agree on the scale. + RunWidth BalanceRunWidth(const std::vector &outcomes, const std::map &reference, + const std::optional &reference_cell, const RunWidth &width) const; // The run's spread of eta and ln width over the crystals, net of their own uncertainty. CrystalPrior RunPrior(const std::vector &fits) const; };