diff --git a/tools/XdsIntegrateParser.cpp b/tools/XdsIntegrateParser.cpp
index 41a908aa..7db549be 100644
--- a/tools/XdsIntegrateParser.cpp
+++ b/tools/XdsIntegrateParser.cpp
@@ -26,24 +26,28 @@ IntegrateMap ParseXdsIntegrateHkl(const std::string& filename) {
         std::istringstream iss(line);
         int32_t h, k, l;
         double I, sigma;
+        double xcal, ycal, zcal, rlp, peak, corr;
+        int32_t maxc;
 
-        if (!(iss >> l >> k >> h >> I >> sigma)) {
+        if (!(iss >> l >> k >> h >> I >> sigma >> xcal >> ycal >> zcal >> rlp >> peak >> corr >> maxc)) {
             continue;
         }
 
         HKLData entry;
-        entry.h = h;
+        entry.h = -h;
         entry.k = k;
         entry.l = l;
         entry.I = I;
         entry.sigma = sigma;
+        entry.rlp = rlp;
+        entry.image_number = zcal;
 
         double v = 0.0;
         while (iss >> v) {
             entry.tail.push_back(v);
         }
 
-        result[hkl_key_16(h, k, l)].push_back(std::move(entry));
+        result[hkl_key_16(-h, k, l)].push_back(std::move(entry));
     }
 
     return result;
@@ -87,30 +91,34 @@ CcHalfByResolutionResult ComputeCcByResolution(
 
         const auto& xds_list = it->second;
 
-        double i_ours = SumIntensity(ours_list);
-        double i_xds = SumIntensity(xds_list);
+        for (const auto &i_x: xds_list) {
+            for (const auto &i_o: ours_list) {
+                if (std::fabs(i_x.image_number - i_o.image_number) > 30.0)
+                    continue;
 
-        int64_t h = ours_list.front().h;
-        int64_t k = ours_list.front().k;
-        int64_t l = ours_list.front().l;
+                int64_t h = i_x.h;
+                int64_t k = i_x.k;
+                int64_t l = i_x.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;
+                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;
+                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;
+                int idx = shell.value();
+                n[idx] += 1;
+                sum_x[idx] += i_o.I;
+                sum_y[idx] += i_x.I;
+                sum_x2[idx] += i_o.I * i_o.I;
+                sum_y2[idx] += i_x.I * i_x.I;
+                sum_xy[idx] += i_o.I * i_x.I;
+            }
+        }
     }
 
     CcHalfByResolutionResult result;
diff --git a/tools/XdsIntegrateParser.h b/tools/XdsIntegrateParser.h
index eab0cef1..776e220b 100644
--- a/tools/XdsIntegrateParser.h
+++ b/tools/XdsIntegrateParser.h
@@ -31,7 +31,11 @@ struct HKLData {
     int64_t h ,k,l;
     double I = 0.0;
     double sigma = 0.0;
-    float last_image = 0;
+    int32_t last_image = 0;
+    int32_t count = 0;
+    float rlp = 0.0;
+    float image_number = -100;
+
     std::vector<double> tail; // any remaining numeric columns
 };
 
diff --git a/tools/jfjoch_extract_hkl.cpp b/tools/jfjoch_extract_hkl.cpp
index 61310e24..da1cad4f 100644
--- a/tools/jfjoch_extract_hkl.cpp
+++ b/tools/jfjoch_extract_hkl.cpp
@@ -86,14 +86,18 @@ int main(int argc, char **argv) {
                         .k = r.k,
                         .l = r.l,
                         .I = r.I,
-                        .last_image = r.image_number
+                        .last_image = i,
+                        .count = 1,
+                        .rlp = r.lp,
+                        .image_number = r.image_number
                     };
                     bool found = false;
                     if (it != agg.end()) {
                         for (auto &val: it->second) {
-                            if (val.last_image == r.image_number - 1) {
+                            if (val.last_image == i - 1) {
                                 val.I += r.I;
-                                val.last_image = r.image_number;
+                                val.last_image = i;
+                                val.count++;
                                 found = true;
                                 break;
                             }
@@ -107,15 +111,25 @@ int main(int argc, char **argv) {
             std::fstream f(output_filename, std::ios::out);
             for (const auto &[key, val]: agg) {
                 if (!val.empty()) {
-                    f << val[0].h << " " << val[0].k << " " << val[0].l << " " << val[0].I;
-                    auto xds_it = xds_result.find(key);
-                    if (xds_it != xds_result.end() && !xds_it->second.empty())
-                        f << " " << xds_it->second[0].I;
-                    f << std::endl;
+
+                    if (xds_result.empty()) {
+                        f << val[0].h << " " << val[0].k << " " << val[0].l << " " << val[0].I;
+                        f << std::endl;
+                    } else {
+                        auto xds_it = xds_result.find(key);
+                        if (xds_it != xds_result.end() && !xds_it->second.empty()) {
+                            f << val[0].h << " " << val[0].k << " " << val[0].l << " " << val[0].I;
+
+                            f << " " << xds_it->second[0].I << " " << 1.0/val[0].rlp << " " << xds_it->second[0].rlp;
+                            f << std::endl;
+                        }
+                    }
+
+
                 }
             }
 
-            auto cc_result = ComputeCcByResolution(latt, agg,xds_result, 1.2, 50.0, 20);
+            auto cc_result = ComputeCcByResolution(latt, agg,xds_result, 1.0, 50.0, 25);
             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;