From 0e4e5dc8f112f3b1c763ee7e67f8c2e02895eb13 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 21 Jan 2026 13:47:56 +0100 Subject: [PATCH] jfjoch_extract_hkl: Work in progress -- CC calculation --- tools/XdsIntegrateParser.cpp | 89 +++++++++++++++++++++++++++++++++++- tools/XdsIntegrateParser.h | 19 +++++++- tools/jfjoch_extract_hkl.cpp | 9 +++- 3 files changed, 114 insertions(+), 3 deletions(-) diff --git a/tools/XdsIntegrateParser.cpp b/tools/XdsIntegrateParser.cpp index 8d94dbdb..41a908aa 100644 --- a/tools/XdsIntegrateParser.cpp +++ b/tools/XdsIntegrateParser.cpp @@ -7,6 +7,8 @@ #include #include +#include "../common/ResolutionShells.h" + IntegrateMap ParseXdsIntegrateHkl(const std::string& filename) { std::ifstream in(filename); if (!in.is_open()) { @@ -25,7 +27,7 @@ IntegrateMap ParseXdsIntegrateHkl(const std::string& filename) { int32_t h, k, l; double I, sigma; - if (!(iss >> h >> k >> l >> I >> sigma)) { + if (!(iss >> l >> k >> h >> I >> sigma)) { continue; } @@ -44,5 +46,90 @@ IntegrateMap ParseXdsIntegrateHkl(const std::string& filename) { result[hkl_key_16(h, k, l)].push_back(std::move(entry)); } + return result; +} + + +namespace { +double SumIntensity(const std::vector& v) { + double s = 0.0; + for (const auto& e : v) + s += e.I; + return s; +} +} // namespace + +CcHalfByResolutionResult ComputeCcByResolution( + const CrystalLattice& lattice, + const IntegrateMap& ours, + const IntegrateMap& xds, + float d_min, + float d_max, + int32_t nshells) { + + ResolutionShells shells(d_min, d_max, nshells); + + std::vector sum_x(nshells, 0.0); + std::vector sum_y(nshells, 0.0); + std::vector sum_x2(nshells, 0.0); + std::vector sum_y2(nshells, 0.0); + std::vector sum_xy(nshells, 0.0); + std::vector n(nshells, 0); + + const Coord astar = lattice.Astar(); + const Coord bstar = lattice.Bstar(); + const Coord cstar = lattice.Cstar(); + + for (const auto& [key, ours_list] : ours) { + auto it = xds.find(key); + if (it == xds.end()) + continue; + + const auto& xds_list = it->second; + + double i_ours = SumIntensity(ours_list); + double i_xds = SumIntensity(xds_list); + + int64_t h = ours_list.front().h; + int64_t k = ours_list.front().k; + int64_t l = ours_list.front().l; + + Coord recip = astar * h + bstar * k + cstar * l; + double recip_len = std::sqrt(recip.x * recip.x + recip.y * recip.y + recip.z * recip.z); + if (recip_len <= 0.0) + continue; + + float d = static_cast(1.0 / recip_len); + auto shell = shells.GetShell(d); + if (!shell) + continue; + + int idx = shell.value(); + n[idx] += 1; + sum_x[idx] += i_ours; + sum_y[idx] += i_xds; + sum_x2[idx] += i_ours * i_ours; + sum_y2[idx] += i_xds * i_xds; + sum_xy[idx] += i_ours * i_xds; + } + + CcHalfByResolutionResult result; + result.cc.resize(nshells, std::numeric_limits::quiet_NaN()); + result.pairs = n; + result.shell_mean_one_over_d2 = shells.GetShellMeanOneOverResSq(); + + for (int i = 0; i < nshells; ++i) { + if (n[i] < 2) + continue; + + double denom_x = n[i] * sum_x2[i] - sum_x[i] * sum_x[i]; + double denom_y = n[i] * sum_y2[i] - sum_y[i] * sum_y[i]; + double denom = std::sqrt(denom_x * denom_y); + if (denom <= 0.0) + continue; + + result.cc[i] = (n[i] * sum_xy[i] - sum_x[i] * sum_y[i]) / denom; + } + return result; } \ No newline at end of file diff --git a/tools/XdsIntegrateParser.h b/tools/XdsIntegrateParser.h index 1d5e4c7e..eab0cef1 100644 --- a/tools/XdsIntegrateParser.h +++ b/tools/XdsIntegrateParser.h @@ -10,6 +10,15 @@ #include #include +#include "../common/CrystalLattice.h" + + +struct CcHalfByResolutionResult { + std::vector cc; + std::vector pairs; + std::vector shell_mean_one_over_d2; +}; + inline uint64_t hkl_key_16(int64_t h, int64_t k, int64_t l) { const uint16_t bias = 512; // maps -512..512 -> 0..1024 uint64_t uh = static_cast(std::clamp(h + bias, int64_t(0), int64_t(1024))); @@ -28,4 +37,12 @@ struct HKLData { using IntegrateMap = std::unordered_map>; -IntegrateMap ParseXdsIntegrateHkl(const std::string& filename); \ No newline at end of file +IntegrateMap ParseXdsIntegrateHkl(const std::string& filename); + +CcHalfByResolutionResult ComputeCcByResolution( + const CrystalLattice& lattice, + const IntegrateMap& ours, + const IntegrateMap& xds, + float d_min, + float d_max, + int32_t nshells); diff --git a/tools/jfjoch_extract_hkl.cpp b/tools/jfjoch_extract_hkl.cpp index 126db9fc..61310e24 100644 --- a/tools/jfjoch_extract_hkl.cpp +++ b/tools/jfjoch_extract_hkl.cpp @@ -7,7 +7,6 @@ #include "../reader/JFJochHDF5Reader.h" #include "../common/print_license.h" #include "../common/Logger.h" -#include "XdsIntegrateParser.h" void print_usage(Logger &logger) { logger.Info("Usage ./jfjoch_extract_hkl {} "); @@ -66,6 +65,7 @@ int main(int argc, char **argv) { if (dataset) { if (sum_neighboring) { + CrystalLattice latt; IntegrateMap agg; int64_t total_images = reader.GetNumberOfImages(); @@ -76,6 +76,8 @@ int main(int argc, char **argv) { for (int i = 0; i < total_images; ++i) { auto dataset = reader.LoadImage(i); if (!dataset) continue; + if (dataset->ImageData().indexing_lattice) + latt = dataset->ImageData().indexing_lattice.value(); for (const auto &r: dataset->ImageData().reflections) { int64_t key = hkl_key_16(r.h, r.k, r.l); auto it = agg.find(key); @@ -112,6 +114,11 @@ int main(int argc, char **argv) { f << std::endl; } } + + auto cc_result = ComputeCcByResolution(latt, agg,xds_result, 1.2, 50.0, 20); + for (int i = 0; i < cc_result.cc.size(); ++i) + std::cout << 1/ std::sqrt(cc_result.shell_mean_one_over_d2[i]) << " " << cc_result.cc[i]* 100.0 << " " << cc_result.pairs[i] << std::endl; + } else { auto geom = dataset->Dataset().experiment.GetDiffractionGeometry();