64 lines
2.0 KiB
C++
64 lines
2.0 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "ScaleOnTheFly.h"
|
|
|
|
#include <ceres/ceres.h>
|
|
|
|
double SafeInv(double x, double fallback) {
|
|
if (!std::isfinite(x) || x == 0.0)
|
|
return fallback;
|
|
return 1.0 / x;
|
|
}
|
|
|
|
struct IntensityRotationResidual {
|
|
IntensityRotationResidual(const Reflection &r, double Itrue, double sigma)
|
|
: Iobs(r.I),
|
|
Itrue(Itrue),
|
|
weight(SafeInv(sigma, 1.0)),
|
|
delta_phi_deg(r.delta_phi_deg),
|
|
lp(SafeInv(r.rlp, 1.0)),
|
|
c1(r.zeta / std::sqrt(2.0)),
|
|
b_resolution_coeff(SafeInv(-r.d * r.d / 4.0, 0.0)) {
|
|
}
|
|
|
|
template<typename T>
|
|
bool operator()(const T *const G,
|
|
const T *const B,
|
|
const T *const mosaicity,
|
|
const T *const wedge,
|
|
T *residual) const {
|
|
const T partiality = T(1.0);
|
|
if (mosaicity > 0) {
|
|
const T half_wedge = wedge[0] / T(2.0);
|
|
const T arg_plus = T(delta_phi_deg + half_wedge) * T(c1) / mosaicity[0];
|
|
const T arg_minus = T(delta_phi_deg - half_wedge) * T(c1) / mosaicity[0];
|
|
partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0);
|
|
}
|
|
const T B_term = ceres::exp(T(B[0]) * b_resolution_coeff);
|
|
residual[0] = (G[0] * partiality * T(lp) * Itrue - T(Iobs)) * T(weight);
|
|
return true;
|
|
}
|
|
|
|
double Iobs;
|
|
double Itrue;
|
|
double weight;
|
|
double delta_phi_deg;
|
|
double lp;
|
|
double c1;
|
|
double b_resolution_coeff;
|
|
};
|
|
|
|
ScaleOnTheFly::ScaleOnTheFly(const std::vector<MergedReflection> &ref,
|
|
std::optional<gemmi::SpaceGroup> sg,
|
|
bool merge_friedel) : sg(sg), merge_friedel(merge_friedel) {
|
|
for (const auto &r : ref) {
|
|
const auto key = CanonicalHKL(r, merge_friedel, sg);
|
|
reference_data[key] = r.I;
|
|
}
|
|
}
|
|
|
|
std::optional<ScalingResult> ScaleOnTheFly::Scale(std::vector<Reflection> &reflections) {
|
|
|
|
}
|