jfjoch_extract_hkl: Work in progress -- CC calculation

tools/XdsIntegrateParser.cpp

@@ -7,6 +7,8 @@
 #include <sstream>
 #include <stdexcept>
 
+#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<HKLData>& 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<double> sum_x(nshells, 0.0);
+    std::vector<double> sum_y(nshells, 0.0);
+    std::vector<double> sum_x2(nshells, 0.0);
+    std::vector<double> sum_y2(nshells, 0.0);
+    std::vector<double> sum_xy(nshells, 0.0);
+    std::vector<int32_t> 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<float>(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<double>::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;
 }

tools/XdsIntegrateParser.h

@@ -10,6 +10,15 @@
 #include <unordered_map>
 #include <cstdint>
 
+#include "../common/CrystalLattice.h"
+
+
+struct CcHalfByResolutionResult {
+    std::vector<double> cc;
+    std::vector<int32_t> pairs;
+    std::vector<float> 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<uint16_t>(std::clamp(h + bias, int64_t(0), int64_t(1024)));
@@ -28,4 +37,12 @@ struct HKLData {
 
 using IntegrateMap = std::unordered_map<uint64_t, std::vector<HKLData>>;
 
-IntegrateMap ParseXdsIntegrateHkl(const std::string& filename);
+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);

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 {<options>} <path to master file>");
@@ -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();