diff --git a/CMakeLists.txt b/CMakeLists.txt index 92e97835..02668e08 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,6 +16,8 @@ SET(JFJOCH_INSTALL_DRIVER_SOURCE OFF CACHE BOOL "Install kernel driver source (i SET(JFJOCH_USE_CUDA ON CACHE BOOL "Compile Jungfraujoch with CUDA") SET(JFJOCH_VIEWER_BUILD OFF CACHE BOOL "Compile Jungfraujoch viewer") +FIND_PACKAGE(ZLIB REQUIRED) + SET (ZLIB_USE_STATIC_LIBS TRUE) OPTION(SLS9 "Build with sls_detector_package v9.2.0" OFF) diff --git a/common/CrystalLattice.h b/common/CrystalLattice.h index c5cce66a..45952cd3 100644 --- a/common/CrystalLattice.h +++ b/common/CrystalLattice.h @@ -4,12 +4,12 @@ #ifndef JUNGFRAUJOCH_CRYSTALLATTICE_H #define JUNGFRAUJOCH_CRYSTALLATTICE_H -#include "../symmetry/gemmi/math.hpp" +#include +#include #include #include #include "Coord.h" #include "UnitCell.h" -#include "../symmetry/gemmi/symmetry.hpp" class CrystalLattice { Coord vec[3]; diff --git a/common/DiffractionExperiment.h b/common/DiffractionExperiment.h index 95cc74c8..11e82816 100644 --- a/common/DiffractionExperiment.h +++ b/common/DiffractionExperiment.h @@ -29,7 +29,7 @@ #include "BraggIntegrationSettings.h" #include "ScalingSettings.h" -#include "../symmetry/gemmi/symmetry.hpp" +#include enum class DetectorMode { Standard, PedestalG0, PedestalG1, PedestalG2, DarkMask diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 024a19bc..2a2477bd 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -22,7 +22,7 @@ #include "CrystalLattice.h" #include "IndexingSettings.h" #include "XrayFluorescenceSpectrum.h" -#include "../symmetry/gemmi/symmetry.hpp" +#include constexpr const uint64_t user_data_release = 6; constexpr const uint64_t user_data_magic_number = 0x52320000UL | user_data_release; diff --git a/gemmi_gph/CMakeLists.txt b/gemmi_gph/CMakeLists.txt index 66d88819..c7353257 100644 --- a/gemmi_gph/CMakeLists.txt +++ b/gemmi_gph/CMakeLists.txt @@ -5,4 +5,4 @@ ADD_LIBRARY(gemmi STATIC symmetry.cpp gz.cpp mtz.cpp sprintf.cpp xds_ascii.cpp gemmi/unitcell.hpp gemmi/math.hpp) TARGET_INCLUDE_DIRECTORIES(gemmi PUBLIC .) -TARGET_LINK_LIBRARIES(gemmi ) \ No newline at end of file +TARGET_LINK_LIBRARIES(gemmi PRIVATE ZLIB::ZLIB) \ No newline at end of file diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 394dbaed..e8bc48a8 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -30,7 +30,9 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC RotationParameters.cpp RotationParameters.h WriteMmcif.cpp - WriteMmcif.h) + WriteMmcif.h + LoadFCalcFromMtz.cpp + LoadFCalcFromMtz.h) FIND_PACKAGE(Eigen3 3.4 REQUIRED NO_MODULE) # provides Eigen3::Eigen diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 9b749e89..6a0e62e1 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -280,15 +280,15 @@ void IndexAndRefine::ScaleImage(size_t n, ScaleOnTheFly &scaling, ScalingResult result.rotation_wedge_deg[n] = res.wedge; } -ScalingResult IndexAndRefine::ScaleAllImages(size_t nthreads) { - auto merge_result = MergeAll(experiment, reflections); - ScaleOnTheFly scaling(merge_result, experiment); +ScalingResult IndexAndRefine::ScaleAllImages(const std::vector &reference, size_t nthreads) { + ScaleOnTheFly scaling(reference, experiment); return scaling.Scale(reflections, mosaicity, nthreads); } -MergeResult IndexAndRefine::Merge() const { +MergeResult IndexAndRefine::Merge(bool calc_statistics) const { MergeResult out; out.merged = MergeAll(experiment, reflections); - out.statistics = MergeStats(experiment, out.merged, reflections); + if (calc_statistics) + out.statistics = MergeStats(experiment, out.merged, reflections); return out; } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index f80335e3..9acca8a1 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -63,9 +63,10 @@ class IndexAndRefine { public: IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer); void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction); + IndexAndRefine& ReferenceIntensities(std::vector &reference); - ScalingResult ScaleAllImages(size_t nthreads = 0); - MergeResult Merge() const; + ScalingResult ScaleAllImages(const std::vector &reference, size_t nthreads = 0); + MergeResult Merge(bool statistics) const; std::optional Finalize(); }; diff --git a/image_analysis/LoadFCalcFromMtz.cpp b/image_analysis/LoadFCalcFromMtz.cpp new file mode 100644 index 00000000..b7153f32 --- /dev/null +++ b/image_analysis/LoadFCalcFromMtz.cpp @@ -0,0 +1,47 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "LoadFCalcFromMtz.h" + +#include +#include +#include +#include + +#include + +#include "../common/Reflection.h" + +std::vector LoadFCalcFromMtz(const std::string& path) { + gemmi::Mtz mtz; + mtz.read_file_gz(path, true); + + const gemmi::Mtz::Column* fc = mtz.column_with_label("F-model", nullptr, 'F'); + if (fc == nullptr) + throw std::runtime_error("MTZ does not contain F-model column"); + + std::vector result; + result.reserve(static_cast(mtz.nreflections)); + + const std::size_t stride = mtz.columns.size(); + + for (int i = 0; i < mtz.nreflections; ++i) { + const std::size_t row = static_cast(i) * stride; + const float f = (*fc)[static_cast(i)]; + + if (std::isnan(f)) + continue; + + MergedReflection r; + r.h = static_cast(mtz.data[row + 0]); + r.k = static_cast(mtz.data[row + 1]); + r.l = static_cast(mtz.data[row + 2]); + r.I = f * f; + r.sigma = NAN; + r.d = 0.0f; + + result.emplace_back(r); + } + + return result; +} diff --git a/image_analysis/LoadFCalcFromMtz.h b/image_analysis/LoadFCalcFromMtz.h new file mode 100644 index 00000000..ad2f29c9 --- /dev/null +++ b/image_analysis/LoadFCalcFromMtz.h @@ -0,0 +1,8 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "../common/Reflection.h" + +std::vector LoadFCalcFromMtz(const std::string& path); diff --git a/image_analysis/WriteMmcif.h b/image_analysis/WriteMmcif.h index 10e05b22..057ee56b 100644 --- a/image_analysis/WriteMmcif.h +++ b/image_analysis/WriteMmcif.h @@ -10,7 +10,7 @@ #include "scale_merge/FrenchWilson.h" #include "../common/UnitCell.h" -#include "../symmetry/gemmi/symmetry.hpp" +#include /// Metadata needed to write a meaningful mmCIF reflection file. struct MmcifMetadata { diff --git a/image_puller/CMakeLists.txt b/image_puller/CMakeLists.txt index 51f875a6..24d01ab1 100644 --- a/image_puller/CMakeLists.txt +++ b/image_puller/CMakeLists.txt @@ -5,4 +5,4 @@ ADD_LIBRARY(JFJochImagePuller ZMQImagePuller.cpp ZMQImagePuller.h TestImagePuller.h TCPImagePuller.cpp TCPImagePuller.h) -TARGET_LINK_LIBRARIES(JFJochImagePuller JFJochZMQ JFJochLogger) +TARGET_LINK_LIBRARIES(JFJochImagePuller JFJochCommon JFJochZMQ JFJochLogger) diff --git a/jungfrau/CMakeLists.txt b/jungfrau/CMakeLists.txt index 072a6900..adc83356 100644 --- a/jungfrau/CMakeLists.txt +++ b/jungfrau/CMakeLists.txt @@ -5,4 +5,6 @@ ADD_LIBRARY(JFCalibration STATIC JFModuleGainCalibration.cpp JFModuleGainCalibration.h JFPedestalCalc.cpp JFPedestalCalc.h) +TARGET_LINK_LIBRARIES(JFCalibration JFJochCommon) + SET_SOURCE_FILES_PROPERTIES(JFPedestalCalc.cpp JFConversionFloatingPoint.cpp PROPERTIES COMPILE_FLAGS -Ofast) diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 62343d27..9c2d5769 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -25,6 +25,7 @@ #include "../image_analysis/IndexAndRefine.h" #include "../receiver/JFJochReceiverPlots.h" #include "../compression/JFJochCompressor.h" +#include "../image_analysis/LoadFCalcFromMtz.h" #include "../image_analysis/scale_merge/FrenchWilson.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" #include "../image_analysis/WriteMmcif.h" @@ -149,6 +150,7 @@ int main(int argc, char **argv) { bool refine_bfactor = false; bool refine_wedge = false; std::optional wedge_for_scaling; + std::string ref_mtz; IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto; @@ -163,11 +165,14 @@ int main(int argc, char **argv) { } int opt; - while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:MP:AD:C:T:t:Bw::")) != -1) { + while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:MP:AD:C:T:t:Bw::z:")) != -1) { switch (opt) { case 'o': output_prefix = optarg; break; + case 'z': + ref_mtz = optarg; + break; case 'N': nthreads = atoi(optarg); break; @@ -319,6 +324,12 @@ int main(int argc, char **argv) { logger.Info("Loaded dataset from {}", input_file); + std::vector reference_data; + if (!ref_mtz.empty()) { + reference_data = LoadFCalcFromMtz(ref_mtz); + logger.Info("Loaded {} reflections from {} MTZ file", reference_data.size(), ref_mtz); + } + // 2. Setup Experiment & Components DiffractionExperiment experiment(dataset->experiment); experiment.BitDepthImage(32).Compression(CompressionAlgorithm::BSHUF_LZ4); @@ -579,10 +590,17 @@ int main(int argc, char **argv) { auto scale_start = std::chrono::steady_clock::now(); for (int i = 0; i < 3; i++) { auto iter_start = std::chrono::steady_clock::now(); - auto scale_result = indexer.ScaleAllImages(); - end_msg.image_scale_factor = scale_result.image_scale_g; - scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); + if (reference_data.empty()) { + auto merge_result = indexer.Merge(false); + auto scale_result = indexer.ScaleAllImages(merge_result.merged); + end_msg.image_scale_factor = scale_result.image_scale_g; + scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); + } else { + auto scale_result = indexer.ScaleAllImages(reference_data); + end_msg.image_scale_factor = scale_result.image_scale_g; + scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); + } auto iter_end = std::chrono::steady_clock::now(); double iter_time = std::chrono::duration(iter_end - iter_start).count(); @@ -592,7 +610,7 @@ int main(int argc, char **argv) { double scale_time = std::chrono::duration(scale_end - scale_start).count(); auto merge_start = std::chrono::steady_clock::now(); - auto merge_result = indexer.Merge(); + auto merge_result = indexer.Merge(true); auto merge_end = std::chrono::steady_clock::now(); double merge_time = std::chrono::duration(merge_end - merge_start).count();