diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 80e94cd33..9330bf040 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -78,6 +78,16 @@ high-throughput reduction and analysis of serial femtosecond X-ray diffraction d J. Appl. Cryst. 47, 1118-1131 [doi:10.1107/S1600576714007626](https://doi.org/10.1107/S1600576714007626). +**Per-crystal post-refinement of serial stills** - refining, for each still, the partiality +parameters the spot positions cannot pin (a small reorientation, the crystal's own rocking width and the +radial scale of its reciprocal lattice, which a shot's photon energy or a crystal's cell moves) against +the merged intensities follows the post-refinement of PRIME and the mosaic-block model of Sauter et al.: +M. Uervirojnangkoorn, O. B. Zeldin, A. Y. Lyubimov, J. Hattne, A. S. Brewster, N. K. Sauter, A. T. Brunger +and W. I. Weis, "Enabling X-ray free electron laser crystallography for challenging biological systems +from a limited number of crystals" (2015), eLife 4, e05421; N. K. Sauter, J. Hattne, A. S. Brewster, +N. Echols, P. H. Zwart and P. D. Adams, "Improved crystal orientation and physical properties from +single-shot XFEL stills" (2014), Acta Cryst. D70, 3299-3309 (DOIs not verified here, so not given). + Spot extraction groups strong pixels into spots with the sparse connected-component labelling of the ACTS traccc project: P. Gessinger, H. M. Gray, A. Krasznahorkay, C. Leggett, J. Niermann, A. Salzburger, S. N. Swatman and B. Yeo, "traccc: GPU track reconstruction library for HEP diff --git a/image_analysis/scale_merge/StillsPartialityRefine.cpp b/image_analysis/scale_merge/StillsPartialityRefine.cpp index 6096da7f4..89d5413db 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.cpp +++ b/image_analysis/scale_merge/StillsPartialityRefine.cpp @@ -25,6 +25,10 @@ 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 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; + double SafeInv(double x, double fallback) { if (!std::isfinite(x) || x == 0.0) return fallback; @@ -79,30 +83,6 @@ namespace { } }; - // Analytic partiality for a reflection whose base reciprocal vector is q, tilted by (psi_x, psi_y). - // Mirrors BraggPrediction: dist_ewald = |S| - 1/lambda with S = q_rot + S0, and - // p = exp(-dist_ewald^2 / 2 sigma^2), sigma^2 = gamma0^2 + (gamma_e*d*)^2 + (bw*|q_z|)^2 - three - // independent broadenings added in quadrature: the reciprocal-lattice point's own radius (gamma0, - // ~1/domain size, resolution-INdependent), the mosaic/divergence spread (gamma_e*d*, proportional to - // d*) and the bandwidth smear along the beam. In practice the fit below returns gamma_e ~ 0 and the - // width is essentially gamma0 - see there. - double ComputeP(double qx, double qy, double qz, - double psi_x, double psi_y, - double s0x, double s0y, double s0z, double inv_lambda, - const Width &width) { - double aa[3] = {psi_x, psi_y, 0.0}; - double q[3] = {qx, qy, qz}; - double qr[3]; - ceres::AngleAxisRotatePoint(aa, q, qr); - const double Sx = qr[0] + s0x, Sy = qr[1] + s0y, Sz = qr[2] + s0z; - const double de = std::sqrt(Sx * Sx + Sy * Sy + Sz * Sz) - inv_lambda; - const double dstar = std::sqrt(qr[0] * qr[0] + qr[1] * qr[1] + qr[2] * qr[2]); - const double sig2 = width.Sigma2(dstar, qr[2]); - if (!(sig2 > 0.0)) - return 1.0; - return std::exp(-0.5 * de * de / sig2); - } - // The prediction window cuts the rocking curve. Reflections fill a window of half-width W evenly, so // their plain moment of de is u = W^2 / 3, and the intensity-weighted moment m of a Gaussian of width // s cut at +-W is u * g(a), a = W / s, with g(a) = 3/a^2 (1 - 2 a phi(a) / erf(a / sqrt 2)): 1 for a @@ -170,79 +150,118 @@ namespace { -s * ky, s * kx, c}; } - // ComputeP with the tilt given as its rotation matrix. - struct TiltModel { + // A crystal's model of its partialities: the orientation tilt (psi_x, psi_y) about the two axes + // perpendicular to the beam, a radial scale eta of its reciprocal lattice (q -> (1 + eta) q: a + // wavelength or a cell a little off for this crystal - on an XFEL the photon energy of its shot) and a + // scale exp(ln_width) of the run's rocking width (a crystal's own mosaicity or domain size). + // Per-crystal mosaicity and lattice scale in post-refinement, following Uervirojnangkoorn et al. (2015) + // eLife 4, e05421 and Sauter et al. (2014) Acta Cryst. D70, 3299-3309. + struct CrystalModel { + double psi_x = 0.0, psi_y = 0.0, eta = 0.0, ln_width = 0.0; + }; + + // The partiality of a reflection whose base reciprocal vector is q under a crystal model. Mirrors + // BraggPrediction: dist_ewald = |S| - 1/lambda with S = q' + S0, q' = (1 + eta) R(psi) q, and + // p = exp(-dist_ewald^2 / 2 sigma^2), sigma = exp(ln_width) sigma_run(d*). + + struct PartialityModel { double s0x, s0y, s0z, inv_lambda; Width width; - double P(const std::array &R, const FitObs &o) const { - const double qx = R[0] * o.qx + R[1] * o.qy + R[2] * o.qz; - const double qy = R[3] * o.qx + R[4] * o.qy + R[5] * o.qz; - const double qz = R[6] * o.qx + R[7] * o.qy + R[8] * o.qz; + double At(const std::array &R, const CrystalModel &m, double qx0, double qy0, double qz0) const { + const double f = 1.0 + m.eta; + const double qx = f * (R[0] * qx0 + R[1] * qy0 + R[2] * qz0); + const double qy = f * (R[3] * qx0 + R[4] * qy0 + R[5] * qz0); + const double qz = f * (R[6] * qx0 + R[7] * qy0 + R[8] * qz0); const double de = VecLen(qx + s0x, qy + s0y, qz + s0z) - inv_lambda; - const double sig2 = width.Sigma2(VecLen(qx, qy, qz), qz); + const double sig2 = width.Sigma2(VecLen(qx, qy, qz), qz) * std::exp(2.0 * m.ln_width); if (!(sig2 > 0.0)) return 1.0; return std::exp(-0.5 * de * de / sig2); } + + double P(const std::array &R, const CrystalModel &m, const FitObs &o) const { + return At(R, m, o.qx, o.qy, o.qz); + } }; - // Refine the orientation tilt (psi_x, psi_y) holding the scale G fixed: minimise - // sum_i rho(r_i^2) / 2 + |psi / sigma_prior|^2 / 2, r_i = weight_i (G p_i(psi) Iref_i / corr_i - Iobs_i), - // with the Cauchy rho(s) = k^2 log(1 + s / k^2) - the intensity-space objective of ScaleOnTheFly, and - // the Gaussian prior a proper MAP term since r_i is already in sigma units. Levenberg-Marquardt on the - // normal equations reweighted by rho', the Jacobian by central differences, every trial tilt held inside - // the bound. Returns false when the cost cannot be evaluated, leaving psi as it came in. - bool RefineTilt(const std::vector &obs, const TiltModel &model, double G, double robust_k, - double inv_sigma_prior, double max_tilt, double psi[2]) { + // Refine the crystal model holding the scale G fixed: minimise + // sum_i rho(r_i^2) / 2 + sum_k (x_k / sigma_k)^2 / 2, r_i = weight_i (G p_i(x) Iref_i / corr_i - Iobs_i), + // over x = (psi_x, psi_y, eta, ln_width), with the Cauchy rho(s) = k^2 log(1 + s / k^2) - the intensity- + // space objective of ScaleOnTheFly - and Gaussian priors about the run (x = 0), proper MAP terms since + // r_i is already in sigma units. Levenberg-Marquardt on the normal equations reweighted by rho', the + // Jacobian by central differences, every trial held inside the bounds. Returns false when the cost + // cannot be evaluated, leaving the model as it came in. + bool RefineCrystal(const std::vector &obs, const PartialityModel &model, double G, double robust_k, + const std::array &inv_sigma_prior, const std::array &bound, + CrystalModel &crystal, std::array &variance) { constexpr int MAX_ITERATIONS = 25; - constexpr double STEP = 1e-7; // rad, for the central differences + constexpr std::array STEP = {1e-7, 1e-7, 1e-7, 1e-5}; // for the central differences const double k2 = robust_k * robust_k; - const double prior2 = inv_sigma_prior * inv_sigma_prior; + auto as_model = [](const Eigen::Vector4d &x) { return CrystalModel{x(0), x(1), x(2), x(3)}; }; - auto cost = [&](const double p[2]) { - const auto R = TiltMatrix(p[0], p[1]); - double c = 0.5 * prior2 * (p[0] * p[0] + p[1] * p[1]); + auto cost = [&](const Eigen::Vector4d &x) { + const auto R = TiltMatrix(x(0), x(1)); + const CrystalModel m = as_model(x); + double c = 0.0; + for (int k = 0; k < 4; ++k) + c += 0.5 * inv_sigma_prior[k] * inv_sigma_prior[k] * x(k) * x(k); for (const auto &o: obs) { - const double r = o.weight * (G * model.P(R, o) * o.inv_corr * o.Iref - o.Iobs); + const double r = o.weight * (G * model.P(R, m, o) * o.inv_corr * o.Iref - o.Iobs); c += 0.5 * k2 * std::log1p(r * r / k2); } return c; }; - double current = cost(psi); + Eigen::Vector4d x(crystal.psi_x, crystal.psi_y, crystal.eta, crystal.ln_width); + double current = cost(x); if (!std::isfinite(current)) return false; double lambda = 1e-4; + Eigen::Matrix4d data_A = Eigen::Matrix4d::Identity(); // the normal matrix of the data alone for (int iter = 0; iter < MAX_ITERATIONS; ++iter) { - const auto R = TiltMatrix(psi[0], psi[1]); - const auto Rxp = TiltMatrix(psi[0] + STEP, psi[1]), Rxm = TiltMatrix(psi[0] - STEP, psi[1]); - const auto Ryp = TiltMatrix(psi[0], psi[1] + STEP), Rym = TiltMatrix(psi[0], psi[1] - STEP); - double A00 = prior2, A01 = 0.0, A11 = prior2; - double b0 = prior2 * psi[0], b1 = prior2 * psi[1]; + Eigen::Matrix4d A = Eigen::Matrix4d::Zero(); + Eigen::Vector4d b = Eigen::Vector4d::Zero(); + for (int k = 0; k < 4; ++k) { + A(k, k) = inv_sigma_prior[k] * inv_sigma_prior[k]; + b(k) = A(k, k) * x(k); + } + // The model at x and at x +- one step in each parameter. + std::array, 9> R; + std::array m; + for (int j = 0; j < 9; ++j) { + Eigen::Vector4d xj = x; + if (j > 0) + xj((j - 1) / 2) += (j % 2 ? 1.0 : -1.0) * STEP[(j - 1) / 2]; + R[j] = TiltMatrix(xj(0), xj(1)); + m[j] = as_model(xj); + } for (const auto &o: obs) { const double a = o.weight * G * o.inv_corr * o.Iref; - const double r = a * model.P(R, o) - o.weight * o.Iobs; - const double jx = a * (model.P(Rxp, o) - model.P(Rxm, o)) / (2.0 * STEP); - const double jy = a * (model.P(Ryp, o) - model.P(Rym, o)) / (2.0 * STEP); + const double r = a * model.P(R[0], m[0], o) - o.weight * o.Iobs; + Eigen::Vector4d J; + for (int k = 0; k < 4; ++k) + J(k) = a * (model.P(R[2 * k + 1], m[2 * k + 1], o) - model.P(R[2 * k + 2], m[2 * k + 2], o)) + / (2.0 * STEP[k]); const double w = 1.0 / (1.0 + r * r / k2); // rho' - A00 += w * jx * jx; - A01 += w * jx * jy; - A11 += w * jy * jy; - b0 += w * jx * r; - b1 += w * jy * r; + A += w * J * J.transpose(); + b += w * r * J; } + data_A = A; + for (int k = 0; k < 4; ++k) + data_A(k, k) -= inv_sigma_prior[k] * inv_sigma_prior[k]; bool accepted = false; - double trial[2] = {psi[0], psi[1]}; + Eigen::Vector4d trial = x; double trial_cost = current; while (lambda < 1e10) { - const double M00 = A00 * (1.0 + lambda), M11 = A11 * (1.0 + lambda); - const double det = M00 * M11 - A01 * A01; - if (!(det > 0.0)) + Eigen::Matrix4d M = A; + M.diagonal() *= 1.0 + lambda; + const Eigen::Vector4d dx = M.ldlt().solve(b); + if (!dx.allFinite()) break; - trial[0] = std::clamp(psi[0] - (M11 * b0 - A01 * b1) / det, -max_tilt, max_tilt); - trial[1] = std::clamp(psi[1] - (M00 * b1 - A01 * b0) / det, -max_tilt, max_tilt); + for (int k = 0; k < 4; ++k) + trial(k) = std::clamp(x(k) - dx(k), -bound[k], bound[k]); trial_cost = cost(trial); if (std::isfinite(trial_cost) && trial_cost < current) { accepted = true; @@ -253,14 +272,17 @@ namespace { } if (!accepted) break; - const double step = std::hypot(trial[0] - psi[0], trial[1] - psi[1]); const double decrease = current - trial_cost; - psi[0] = trial[0]; - psi[1] = trial[1]; + const double step = (trial - x).norm(); + x = trial; current = trial_cost; if (decrease <= 1e-6 * current || step <= 1e-10) break; } + crystal = as_model(x); + const Eigen::Matrix4d cov = data_A.inverse(); + for (int k = 0; k < 4; ++k) + variance[k] = cov(k, k); return true; } } @@ -273,12 +295,12 @@ StillsPartialityRefine::StillsPartialityRefine(const DiffractionExperiment &x) min_partiality_(x.GetScalingSettings().GetMinPartiality()), bandwidth_sigma_(x.GetBandwidthFWHM().value_or(0.0f) / 2.3548f) {} -double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, +StillsPartialityRefine::CrystalFit StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, const std::map &reference, const std::optional &run_width, - double radial_scale) const { + double radial_scale, const CrystalPrior &prior) const { if (outcome.reflections.empty()) - return 0.0; + return {}; const Coord Astar = outcome.latt.Astar(); const Coord Bstar = outcome.latt.Bstar(); @@ -337,7 +359,7 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, } if (obs.size() < MIN_FIT_REFLECTIONS || n_de == 0) - return 0.0; + return {}; // Without the run's measured width (too few reflections in the run to measure it), the width comes // from this crystal's own reflections. Solve the 2x2 normal equations for de^2 = A + B d*^2. A degenerate spread in d* (all reflections in @@ -372,44 +394,47 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, .table_dstar = run_width ? &run_width->table_dstar : nullptr, .table_s = run_width ? &run_width->table_s : nullptr}; - double psi[2] = {0.0, 0.0}; + CrystalModel crystal; double G = 1.0; - const TiltModel model{.s0x = S0.x, .s0y = S0.y, .s0z = S0.z, .inv_lambda = inv_lambda, - .width = width}; - const bool refine_tilt = obs.size() >= settings_.min_reflections; - const int inner = refine_tilt ? settings_.inner_iterations : 1; + const PartialityModel model{.s0x = S0.x, .s0y = S0.y, .s0z = S0.z, .inv_lambda = inv_lambda, .width = width}; + const bool refine_crystal = obs.size() >= settings_.min_reflections; + const int inner = refine_crystal ? settings_.inner_iterations : 1; + const std::array inv_sigma_prior = {kRadToDeg / settings_.prior_sigma_deg, + kRadToDeg / settings_.prior_sigma_deg, + 1.0 / prior.sigma_eta, 1.0 / prior.sigma_ln_width}; + std::array variance{NAN, NAN, NAN, NAN}; + const std::array bound = {settings_.max_tilt_rad, settings_.max_tilt_rad, + settings_.max_eta, settings_.max_ln_width}; for (int it = 0; it < inner; ++it) { - // (1) Solve G given the current partialities, and weight the tilt fit by the variances at that model. + // (1) Solve G given the current partialities, and weight the model fit by the variances at that model. + const auto R = TiltMatrix(crystal.psi_x, crystal.psi_y); std::vector coeff(obs.size()); - for (size_t j = 0; j < obs.size(); ++j) { - const double p = ComputeP(obs[j].qx, obs[j].qy, obs[j].qz, psi[0], psi[1], - S0.x, S0.y, S0.z, inv_lambda, width); - coeff[j] = p * obs[j].inv_corr * obs[j].Iref; - } + for (size_t j = 0; j < obs.size(); ++j) + coeff[j] = model.P(R, crystal, obs[j]) * obs[j].inv_corr * obs[j].Iref; G = SolveScaleIRLS(obs, coeff, settings_.robust_k); if (!(G > 0.0) || !std::isfinite(G)) - return 0.0; + return {}; for (size_t j = 0; j < obs.size(); ++j) obs[j].weight = 1.0 / std::sqrt(ModelVariance(obs[j], G * coeff[j])); - if (!refine_tilt) + if (!refine_crystal) break; - // (2) Refine the tilt holding G fixed. A cost that cannot be evaluated keeps the tilt this crystal - // came in with and stops refining it - G alone is still a usable model. - const double inv_sigma_prior = settings_.prior_sigma_deg > 0.0 ? kRadToDeg / settings_.prior_sigma_deg : 0.0; - if (!RefineTilt(obs, model, G, settings_.robust_k, inv_sigma_prior, settings_.max_tilt_rad, psi)) + // (2) Refine the crystal model holding G fixed. A cost that cannot be evaluated keeps the model this + // crystal came in with and stops refining it - G alone is still a usable model. + if (!RefineCrystal(obs, model, G, settings_.robust_k, inv_sigma_prior, bound, crystal, variance)) break; } // 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 * qe_corr * flight_corr / // (partiality * G), the same composition ScaleOnTheFly writes. + const auto R = TiltMatrix(crystal.psi_x, crystal.psi_y); 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, width); + const double p = model.At(R, crystal, q.x, q.y, q.z); r.partiality = static_cast(p); const float corr = r.prescaling_corr * r.qe_corr * r.flight_corr; const double denom = p * G; @@ -428,8 +453,8 @@ double StillsPartialityRefine::RefineOne(IntegrationOutcome &outcome, outcome.image_scale_cc = cc; outcome.image_scale_cc_n = cc_n; - const double tilt_deg = std::sqrt(psi[0] * psi[0] + psi[1] * psi[1]) * kRadToDeg; - return tilt_deg; + return {.refined = refine_crystal, .tilt_deg = std::hypot(crystal.psi_x, crystal.psi_y) * kRadToDeg, + .eta = crystal.eta, .ln_width = crystal.ln_width, .var_eta = variance[2], .var_ln_width = variance[3]}; } double StillsPartialityRefine::MeasureRadialScale(const std::vector &outcomes, @@ -718,6 +743,40 @@ std::optional StillsPartialityRefine::MeasureR return width; } +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). + auto spread = [&](auto value, auto variance, double floor) { + std::vector x, v; + for (const auto &f: fits) + if (f.refined && std::isfinite(variance(f)) && variance(f) > 0.0) { + x.push_back(value(f)); + v.push_back(variance(f)); + } + if (x.size() < 10) + return floor; + double mean = 0.0; + for (double e: x) + mean += e; + mean /= static_cast(x.size()); + double var = 0.0; + for (double e: x) + var += (e - mean) * (e - mean); + var /= static_cast(x.size() - 1); + std::nth_element(v.begin(), v.begin() + v.size() / 2, v.end()); + return std::max(floor, std::sqrt(std::max(0.0, var - v[v.size() / 2]))); + }; + const CrystalPrior prior{ + .sigma_eta = spread([](const CrystalFit &f) { return f.eta; }, + [](const CrystalFit &f) { return f.var_eta; }, MIN_PRIOR_ETA), + .sigma_ln_width = spread([](const CrystalFit &f) { return f.ln_width; }, + [](const CrystalFit &f) { return f.var_ln_width; }, MIN_PRIOR_LN_WIDTH)}; + Logger("StillsPartialityRefine").Info("Stills crystals differ in radial scale by {:.1e} and in rocking width " + "by {:.0f}% (sd over the crystals, net of their own uncertainty)", + prior.sigma_eta, 100.0 * prior.sigma_ln_width); + return prior; +} + double StillsPartialityRefine::Run(std::vector &outcomes, size_t nthreads, const std::optional &reference_cell, const std::optional &given_width, @@ -732,7 +791,6 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si nthreads = AvailableCpus(); nthreads = std::max(1, nthreads); - double last_mean_tilt = 0.0; // Where the frames were predicted through a curve measured before the image pass, the partialities use // that same curve. Measured again here, through windows that follow it, the width reads narrower at low // resolution (its constant term falls to about zero), and the merge it gave moved R-free against the @@ -750,6 +808,13 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si 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 + // uncertainty (method of moments) is how much crystals really differ. On an XFEL each shot has its + // photon energy and the spread is large; on a monochromatic or a stable pink beam it comes out small, + // and the crystals are held to the run's curve. + CrystalPrior prior{.sigma_eta = settings_.max_eta, .sigma_ln_width = settings_.max_ln_width}; + std::vector fits(outcomes.size()); for (int outer = 0; outer < settings_.outer_iterations; ++outer) { // Reference full intensities from the current corrections. const std::vector merged = MergeAll(experiment_, outcomes, nthreads, reference_cell); @@ -757,19 +822,12 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si for (const auto &m: merged) reference[hkl_key_generator_(m)] = m.I; - std::atomic tilt_sum{0.0}; - std::atomic tilt_n{0}; std::atomic next{0}; - auto worker = [&]() { size_t i = next.fetch_add(1); while (i < outcomes.size()) { - const double t = merged_crystal(outcomes[i]) ? RefineOne(outcomes[i], reference, run_width, radial_scale) : 0.0; - if (t > 0.0) { - double prev = tilt_sum.load(); - while (!tilt_sum.compare_exchange_weak(prev, prev + t)) {} - tilt_n.fetch_add(1); - } + fits[i] = merged_crystal(outcomes[i]) ? RefineOne(outcomes[i], reference, run_width, radial_scale, prior) + : CrystalFit{}; i = next.fetch_add(1); } }; @@ -786,8 +844,17 @@ double StillsPartialityRefine::Run(std::vector &outcomes, si f.get(); } - last_mean_tilt = tilt_n > 0 ? tilt_sum.load() / static_cast(tilt_n.load()) : 0.0; + if (outer == 0) + prior = RunPrior(fits); } + double tilt_sum = 0.0; + size_t tilt_n = 0; + for (const auto &f: fits) + if (f.refined) { + tilt_sum += f.tilt_deg; + ++tilt_n; + } + const double last_mean_tilt = tilt_n > 0 ? tilt_sum / static_cast(tilt_n) : 0.0; return last_mean_tilt; } diff --git a/image_analysis/scale_merge/StillsPartialityRefine.h b/image_analysis/scale_merge/StillsPartialityRefine.h index 5b99d6d36..6f644fb8c 100644 --- a/image_analysis/scale_merge/StillsPartialityRefine.h +++ b/image_analysis/scale_merge/StillsPartialityRefine.h @@ -23,10 +23,18 @@ // reshapes the SPATIAL pattern of partialities in a way a single G cannot mimic. That breaks the // degeneracy that killed the scalar-sigma fit. // +// Two more parameters per crystal, each held by a prior about the run: a radial scale eta of its +// reciprocal lattice (q -> (1 + eta) q, which moves the excitation error by about eta lambda d*^2 / 2 - +// on an XFEL the photon energy of the crystal's own shot, which spot positions cannot pin) and a scale of +// the run's rocking width (the crystal's own mosaicity and size). Both change p differently across +// resolution, so neither is degenerate with G. Their priors are not constants: the first cycle refines +// against loose ones, and the run's spread of the estimates, net of their own uncertainty, becomes the +// prior of the next (RunPrior) - large on SASE XFEL data, small on a monochromatic or stable pink beam. +// // Because the integrated intensity I_obs is fixed, refining the tilt only recomputes the partiality // analytically from the stored per-image lattice/geometry (q = A*.h + B*.k + C*.l, then // dist_ewald = |q + S0| - 1/lambda, signed) - NO pixel re-integration. The loop is: merge -> per-crystal -// refine tilt (G profiled out by the same robust IRLS ScaleOnTheFly uses) -> recompute p and +// refine tilt, eta and width (G profiled out by the same robust IRLS ScaleOnTheFly uses) -> recompute p and // image_scale_corr -> re-merge, iterated a few times. Mutates each reflection's `partiality` and // `image_scale_corr` in place; the existing MergeOnTheFly then consumes the improved corrections. The // width of the Gaussian partiality is the run's, measured once from the intensity-weighted excitation @@ -38,8 +46,10 @@ public: int inner_iterations = 3; // (solve G) <-> (refine tilt) alternations per crystal size_t min_reflections = 40; // skip tilt refinement below this (anti-overfit on sparse crystals) double max_tilt_rad = 0.0175; // hard bound on |dpsi| (~1 deg); indexing already refined orientation - double prior_sigma_deg = 0.02; // soft prior pulling dpsi toward 0 (0 = off); tames weak-data overfit, + double prior_sigma_deg = 0.02; // soft prior pulling dpsi toward 0; tames weak-data overfit, // inert on strong data (well-supported tilts overcome it) + double max_eta = 5e-3; // bound on a crystal's radial scale |eta|, and its prior in the first cycle + double max_ln_width = 1.1; // the same for the log of its width scale (a factor of 3) double robust_k = 3.0; // Cauchy loss scale (sigma units) double gamma_e = 0.0; // angular width sigma(d*) = gamma_e*d* (0 = measure from data) }; @@ -84,8 +94,21 @@ private: const double min_partiality_; const float bandwidth_sigma_; - // Refine one crystal against the reference map; returns |dpsi| in degrees (0 if skipped). - double RefineOne(IntegrationOutcome &outcome, - const std::map &reference, - const std::optional &run_width, double radial_scale) const; + // The priors of a crystal's radial scale eta and log width scale (Gaussian, about 0). + struct CrystalPrior { + double sigma_eta, sigma_ln_width; + }; + // One crystal's refined model: |dpsi| in degrees, eta, ln width and the variances of the last two from + // its own data (no prior). refined = false where only G (or nothing) was fitted. + struct CrystalFit { + bool refined = false; + double tilt_deg = 0.0, eta = 0.0, ln_width = 0.0, var_eta = NAN, var_ln_width = NAN; + }; + + // Refine one crystal against the reference map. + CrystalFit RefineOne(IntegrationOutcome &outcome, const std::map &reference, + const std::optional &run_width, double radial_scale, + const CrystalPrior &prior) 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; };