From 155c53acd8d971a8d017d88437a3a363e25cf06c Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 8 Jun 2026 15:37:27 +0200 Subject: [PATCH] jfjoch_process: First pixel refine integration --- common/DatasetSettings.cpp | 9 ++ common/DatasetSettings.h | 3 + common/DiffractionExperiment.cpp | 9 ++ common/DiffractionExperiment.h | 2 + common/IndexingSettings.h | 2 +- image_analysis/CMakeLists.txt | 2 +- image_analysis/IndexAndRefine.cpp | 115 +++++++++++++++--- image_analysis/IndexAndRefine.h | 26 +++- image_analysis/MXAnalysisWithoutFPGA.cpp | 2 +- .../pixel_refinement/PixelRefine.cpp | 26 ++-- image_analysis/pixel_refinement/PixelRefine.h | 8 +- tools/jfjoch_process.cpp | 34 +++++- 12 files changed, 203 insertions(+), 35 deletions(-) diff --git a/common/DatasetSettings.cpp b/common/DatasetSettings.cpp index 85bb81d5..2bcb144b 100644 --- a/common/DatasetSettings.cpp +++ b/common/DatasetSettings.cpp @@ -389,6 +389,15 @@ std::optional DatasetSettings::GetPolarizationFactor() const { return polarization_factor; } +DatasetSettings &DatasetSettings::BandwidthFWHM(const std::optional &input) { + bandwidth_fwhm = input; + return *this; +} + +std::optional DatasetSettings::GetBandwidthFWHM() const { + return bandwidth_fwhm; +} + float DatasetSettings::GetPoniRot3_rad() const { return poni_rot_3_rad; } diff --git a/common/DatasetSettings.h b/common/DatasetSettings.h index 9ec5672e..3d91c630 100644 --- a/common/DatasetSettings.h +++ b/common/DatasetSettings.h @@ -55,6 +55,7 @@ class DatasetSettings { bool write_nxmx_hdf5_master; std::optional polarization_factor; + std::optional bandwidth_fwhm; // relative X-ray bandwidth, FWHM of dlambda/lambda (e.g. 0.01 for 1%) float poni_rot_1_rad; float poni_rot_2_rad; float poni_rot_3_rad; @@ -99,6 +100,7 @@ public: DatasetSettings& PixelValueLowThreshold(const std::optional &input); DatasetSettings& PixelValueHighThreshold(const std::optional &input); DatasetSettings& PolarizationFactor(const std::optional &input); + DatasetSettings& BandwidthFWHM(const std::optional &input); DatasetSettings& PoniRot1_rad(float input); DatasetSettings& PoniRot2_rad(float input); DatasetSettings& PoniRot3_rad(float input); @@ -150,6 +152,7 @@ public: std::optional IsSaveCalibration() const; std::optional GetPolarizationFactor() const; + std::optional GetBandwidthFWHM() const; float GetPoniRot1_rad() const; float GetPoniRot2_rad() const; float GetPoniRot3_rad() const; diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index ffd5f4ad..6b1a3b66 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -1310,6 +1310,15 @@ std::optional DiffractionExperiment::GetPolarizationFactor() const { return dataset.GetPolarizationFactor(); } +DiffractionExperiment &DiffractionExperiment::BandwidthFWHM(const std::optional &input) { + dataset.BandwidthFWHM(input); + return *this; +} + +std::optional DiffractionExperiment::GetBandwidthFWHM() const { + return dataset.GetBandwidthFWHM(); +} + DiffractionExperiment &DiffractionExperiment::SaveCalibration(const std::optional &input) { dataset.SaveCalibration(input); return *this; diff --git a/common/DiffractionExperiment.h b/common/DiffractionExperiment.h index 131c3bff..a8bbb9a8 100644 --- a/common/DiffractionExperiment.h +++ b/common/DiffractionExperiment.h @@ -281,9 +281,11 @@ public: DiffractionExperiment& ApplySolidAngleCorr(bool input); DiffractionExperiment& PolarizationFactor(const std::optional &input); + DiffractionExperiment& BandwidthFWHM(const std::optional &input); bool GetApplySolidAngleCorr() const; std::optional GetPolarizationFactor() const; + std::optional GetBandwidthFWHM() const; int64_t GetUDPInterfaceCount() const; std::vector GetDetectorModuleConfig(const std::vector& net_config) const; diff --git a/common/IndexingSettings.h b/common/IndexingSettings.h index 6e900440..c428f02e 100644 --- a/common/IndexingSettings.h +++ b/common/IndexingSettings.h @@ -6,7 +6,7 @@ #include enum class IndexingAlgorithmEnum {FFBIDX, FFT, FFTW, Auto, None}; -enum class GeomRefinementAlgorithmEnum {None, OrientationOnly, BeamCenter}; +enum class GeomRefinementAlgorithmEnum {None, OrientationOnly, BeamCenter, PixelRefine}; class IndexingSettings { IndexingAlgorithmEnum algorithm; diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 913870c5..1bf6e448 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -52,4 +52,4 @@ ADD_SUBDIRECTORY(image_preprocessing) ADD_SUBDIRECTORY(azint) ADD_SUBDIRECTORY(pixel_refinement) -TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi) +TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge JFJochPixelRefine gemmi) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index d436f2a7..5dc02d55 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -175,6 +175,10 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In XtalOptimizerRotationOnly(data, msg.spots, 0.05); break; case GeomRefinementAlgorithmEnum::BeamCenter: + case GeomRefinementAlgorithmEnum::PixelRefine: + // PixelRefine still benefits from the classical beam-center + cell + // refinement as a starting point; the pixel-level refinement runs + // later, during integration. if (XtalOptimizer(data, {msg.spots})) { outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl()) .BeamY_pxl(data.geom.GetBeamY_pxl()); @@ -234,7 +238,9 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, const SpotFindingSettings &spot_finding_settings, const CompressedImage &image, BraggPrediction &prediction, - const IndexAndRefine::IndexingOutcome &outcome) { + const IndexAndRefine::IndexingOutcome &outcome, + const AzimuthalIntegrationMapping *mapping, + const AzimuthalIntegrationProfile *profile) { if (!outcome.lattice_candidate) return; @@ -286,14 +292,31 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, .mosaicity_deg = std::fabs(mos_deg) }; - auto pred_start_time = std::chrono::steady_clock::now(); - auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction); - auto pred_end_time = std::chrono::steady_clock::now(); - msg.bragg_prediction_time_s = std::chrono::duration(pred_end_time - pred_start_time).count(); + // Select the integration path: classical 2D integration, or the experimental + // PixelRefine joint geometry/profile/scale refinement (which also produces the + // integrated reflections that flow into the normal save/merge). + const bool use_pixel_refine = + experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() == GeomRefinementAlgorithmEnum::PixelRefine + && !pixel_reference_.empty() && mapping && profile; - auto integration_start_time = std::chrono::steady_clock::now(); - i_outcome.reflections = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), nrefl, msg.number); - msg.integrated_reflections = i_outcome.reflections.size(); + if (use_pixel_refine) { + auto integration_start_time = std::chrono::steady_clock::now(); + PixelRefineIntegrate(msg, image, prediction, outcome, *mapping, *profile, i_outcome); + msg.integrated_reflections = i_outcome.reflections.size(); + auto integration_end_time = std::chrono::steady_clock::now(); + msg.integration_time_s = std::chrono::duration(integration_end_time - integration_start_time).count(); + } else { + auto pred_start_time = std::chrono::steady_clock::now(); + auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction); + auto pred_end_time = std::chrono::steady_clock::now(); + msg.bragg_prediction_time_s = std::chrono::duration(pred_end_time - pred_start_time).count(); + + auto integration_start_time = std::chrono::steady_clock::now(); + i_outcome.reflections = BraggIntegrate2D(outcome.experiment, image, prediction.GetReflections(), nrefl, msg.number); + msg.integrated_reflections = i_outcome.reflections.size(); + auto integration_end_time = std::chrono::steady_clock::now(); + msg.integration_time_s = std::chrono::duration(integration_end_time - integration_start_time).count(); + } constexpr size_t kMaxReflections = 10000; if (i_outcome.reflections.size() > kMaxReflections) { @@ -314,10 +337,10 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, CalcISigma(msg, i_outcome.reflections); CalcWilsonBFactor(msg, i_outcome.reflections); - auto integration_end_time = std::chrono::steady_clock::now(); - msg.integration_time_s = std::chrono::duration(integration_end_time - integration_start_time).count(); - - ScaleImage(msg, i_outcome); + // PixelRefine produces already-scaled reflections; only the classical path + // needs the separate ScaleOnTheFly step. + if (!use_pixel_refine) + ScaleImage(msg, i_outcome); // Copy reflections to outgoing message msg.reflections = i_outcome.reflections; @@ -331,7 +354,9 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, void IndexAndRefine::ProcessImage(DataMessage &msg, const SpotFindingSettings &spot_finding_settings, const CompressedImage &image, - BraggPrediction &prediction) { + BraggPrediction &prediction, + const AzimuthalIntegrationMapping *mapping, + const AzimuthalIntegrationProfile *profile) { if (!indexer_ || !spot_finding_settings.indexing) return; @@ -361,7 +386,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, msg.lattice_type = outcome.symmetry; if (spot_finding_settings.quick_integration) - QuickPredictAndIntegrate(msg, spot_finding_settings, image, prediction, outcome); + QuickPredictAndIntegrate(msg, spot_finding_settings, image, prediction, outcome, mapping, profile); } std::optional IndexAndRefine::FinalizeRotationIndexing() { @@ -377,6 +402,7 @@ std::optional IndexAndRefine::FinalizeRotationIndexing() IndexAndRefine &IndexAndRefine::ReferenceIntensities(std::vector &reference) { scaling_engine = std::make_unique(experiment, reference); + pixel_reference_ = reference; // kept for the experimental PixelRefine path return *this; } @@ -396,6 +422,67 @@ void IndexAndRefine::ScaleImage(DataMessage &msg, IntegrationOutcome& outcome) { msg.image_scale_time_s = std::chrono::duration(scaling_end_time - scaling_start_time).count(); } +bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg, + const CompressedImage &image, + BraggPrediction &prediction, + const IndexAndRefine::IndexingOutcome &outcome, + const AzimuthalIntegrationMapping &mapping, + const AzimuthalIntegrationProfile &profile, + IntegrationOutcome &i_outcome) { + if (!outcome.lattice_candidate) + return false; + + // Build the engine once (lazy: needs the azimuthal mapping, known only here). + std::call_once(pixel_refine_once_, [&] { + pixel_refine_ = std::make_unique(experiment, mapping, pixel_reference_); + }); + if (!pixel_refine_) + return false; + + PixelRefineData prd; + prd.geom = outcome.experiment.GetDiffractionGeometry(); + prd.latt = *outcome.lattice_candidate; + prd.crystal_system = outcome.symmetry.crystal_system; + if (prd.crystal_system == gemmi::CrystalSystem::Trigonal) + prd.crystal_system = gemmi::CrystalSystem::Hexagonal; + prd.centering = outcome.symmetry.centering; + if (const auto bw = experiment.GetBandwidthFWHM()) + prd.bandwidth = bw.value() / 2.3548; // FWHM -> sigma + + std::vector buffer; + const uint8_t *ptr = image.GetUncompressedPtr(buffer); + switch (image.GetMode()) { + case CompressedImageMode::Int8: + pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + case CompressedImageMode::Int16: + pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + case CompressedImageMode::Int32: + pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + case CompressedImageMode::Uint8: + pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + case CompressedImageMode::Uint16: + pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + case CompressedImageMode::Uint32: + pixel_refine_->Run(reinterpret_cast(ptr), profile, prediction, prd); break; + default: + return false; + } + + // PixelRefine output flows into the normal save/merge path: the refined + // geometry/lattice and the already-scaled reflections become the outcome. + i_outcome.reflections = std::move(prd.reflections); + i_outcome.geom = prd.geom; + i_outcome.latt = prd.latt; + i_outcome.image_scale_g = static_cast(prd.scale_factor); + i_outcome.image_scale_b_factor_Ang2 = static_cast(prd.B_factor); + + msg.image_scale_factor = static_cast(prd.scale_factor); + if (prd.B_factor != 0.0) + msg.image_scale_b_factor = static_cast(prd.B_factor); + + return true; +} + ScalingResult IndexAndRefine::ScaleAllImages(const std::vector &reference, size_t nthreads) { ScaleOnTheFly scaling(experiment, reference); scaling.Scale(integration_outcome, nthreads); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 0f03ee97..7674392e 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -8,7 +8,10 @@ #include "../common/DiffractionSpot.h" #include "../common/DiffractionExperiment.h" +#include "../common/AzimuthalIntegrationMapping.h" +#include "../common/AzimuthalIntegrationProfile.h" #include "bragg_prediction/BraggPrediction.h" +#include "pixel_refinement/PixelRefine.h" #include "indexing/IndexerThreadPool.h" #include "lattice_search/LatticeSearch.h" #include "rotation_indexer/RotationIndexer.h" @@ -62,17 +65,36 @@ class IndexAndRefine { const SpotFindingSettings &spot_finding_settings, const CompressedImage &image, BraggPrediction &prediction, - const IndexingOutcome &outcome); + const IndexingOutcome &outcome, + const AzimuthalIntegrationMapping *mapping, + const AzimuthalIntegrationProfile *profile); std::unique_ptr scaling_engine; void ScaleImage(DataMessage &msg, IntegrationOutcome& outcome); + + // Experimental PixelRefine integration path (selected via + // GeomRefinementAlgorithmEnum::PixelRefine). Needs reference intensities; the + // engine is built lazily on first use (when the azimuthal mapping is known) + // and is safe to share across threads (prediction is supplied per call). + std::vector pixel_reference_; + std::unique_ptr pixel_refine_; + std::once_flag pixel_refine_once_; + bool PixelRefineIntegrate(DataMessage &msg, + const CompressedImage &image, + BraggPrediction &prediction, + const IndexingOutcome &outcome, + const AzimuthalIntegrationMapping &mapping, + const AzimuthalIntegrationProfile &profile, + IntegrationOutcome &i_outcome); public: IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer); void AddImageToRotationIndexer(DataMessage &msg); void ForceRotationIndexerLattice(const CrystalLattice& lattice); - void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction); + void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction, + const AzimuthalIntegrationMapping *mapping = nullptr, + const AzimuthalIntegrationProfile *profile = nullptr); IndexAndRefine& ReferenceIntensities(std::vector &reference); ScalingResult ScaleAllImages(const std::vector &reference, size_t nthreads = 0); diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 7bade492..6a638213 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -92,7 +92,7 @@ void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, if (spot_finding_settings.indexing) indexer.ProcessImage(output, spot_finding_settings, CompressedImage(preprocessor_buffer->getBuffer(), experiment.GetXPixelsNum(), experiment.GetYPixelsNum()), - *prediction); + *prediction, &integration, &profile); } output.max_viable_pixel_value = ret.max_value; diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 933edb40..f2e4f93c 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -308,10 +308,8 @@ struct PixelResidual { PixelRefine::PixelRefine(const DiffractionExperiment &experiment, const AzimuthalIntegrationMapping &mapping, - const std::vector &reference, - BraggPrediction &prediction) - : prediction(prediction), - mapping(mapping), + const std::vector &reference) + : mapping(mapping), xpixel(experiment.GetXPixelsNum()), ypixel(experiment.GetYPixelsNum()), experiment(experiment), @@ -324,6 +322,7 @@ PixelRefine::PixelRefine(const DiffractionExperiment &experiment, template void PixelRefine::Run(const T *image, const AzimuthalIntegrationProfile &profile, + BraggPrediction &prediction, PixelRefineData &data) { data.solved = false; data.reflections.clear(); @@ -694,7 +693,12 @@ void PixelRefine::Run(const T *image, r.sigma = static_cast(std::sqrt(1.0 / den)); r.bkg = static_cast(bkg_sum / static_cast(n)); r.observed = true; - r.image_scale_corr = (r.partiality > 0.0f) ? r.rlp / r.partiality : NAN; + // Put I onto a common (merge-ready) scale: I * image_scale_corr is the + // full-spot, scale/DW-corrected intensity. Fold in the refined G, the + // per-reflection B_term and the partiality (rlp = 1 here). + const double B_term = std::exp(-data.B_factor / (4.0 * g.d * g.d)); + const double denom = static_cast(r.partiality) * data.scale_factor * B_term; + r.image_scale_corr = (denom > 0.0) ? static_cast(r.rlp / denom) : NAN; } else { r.I = 0.0f; r.sigma = NAN; @@ -706,9 +710,9 @@ void PixelRefine::Run(const T *image, } // Explicit instantiations for the supported (uncompressed) image pixel types. -template void PixelRefine::Run(const int8_t *, const AzimuthalIntegrationProfile &, PixelRefineData &); -template void PixelRefine::Run(const int16_t *, const AzimuthalIntegrationProfile &, PixelRefineData &); -template void PixelRefine::Run(const int32_t *, const AzimuthalIntegrationProfile &, PixelRefineData &); -template void PixelRefine::Run(const uint8_t *, const AzimuthalIntegrationProfile &, PixelRefineData &); -template void PixelRefine::Run(const uint16_t *, const AzimuthalIntegrationProfile &, PixelRefineData &); -template void PixelRefine::Run(const uint32_t *, const AzimuthalIntegrationProfile &, PixelRefineData &); +template void PixelRefine::Run(const int8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const int16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const int32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const uint8_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const uint16_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); +template void PixelRefine::Run(const uint32_t *, const AzimuthalIntegrationProfile &, BraggPrediction &, PixelRefineData &); diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 07b280cf..fec77c46 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -142,7 +142,6 @@ struct PixelRefineData { }; class PixelRefine { - BraggPrediction &prediction; const AzimuthalIntegrationMapping &mapping; const size_t xpixel, ypixel; const DiffractionExperiment &experiment; @@ -152,11 +151,14 @@ class PixelRefine { public: PixelRefine(const DiffractionExperiment &experiment, const AzimuthalIntegrationMapping &mapping, - const std::vector &reference, - BraggPrediction &prediction); + const std::vector &reference); + // The BraggPrediction is supplied per call (it is mutated): this keeps a + // single PixelRefine instance usable from several threads, each passing its + // own prediction buffer. Only `data` is written; PixelRefine state is const. template void Run(const T *image, const AzimuthalIntegrationProfile &profile, + BraggPrediction &prediction, PixelRefineData &data); }; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index f0f32aab..b01c05b1 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -58,7 +58,7 @@ void print_usage() { std::cout << " -X, --indexing-algorithm Indexing algorithm (FFBIDX|FFT|FFTW|Auto|None)" << std::endl; std::cout << " -S, --space-group Space group number - used for both indexing and scaling" << std::endl; std::cout << " -C, --unit-cell Fix reference unit cell: \"a,b,c,alpha,beta,gamma\"" << std::endl; - std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice)" << std::endl; + std::cout << " -r, --refine Geometry refinement algorithm (none|orientation|beam_and_lattice|pixelrefine)" << std::endl; std::cout << std::endl; std::cout << " Scaling and merging" << std::endl; @@ -73,6 +73,10 @@ void print_usage() { std::cout << " --scaling-iterations Number of scaling iterations with no reference data (default: 3)" << std::endl; std::cout << " --scaling-output Output format for scaling results mtz|cif|txt (default: mtz)" << std::endl; std::cout << " -z, --reference-mtz Reference MTZ file" << std::endl; + std::cout << std::endl; + + std::cout << " Pixel refinement (experimental, select via -r pixelrefine, needs --reference-mtz)" << std::endl; + std::cout << " --bandwidth Relative X-ray bandwidth FWHM (e.g. 0.01 for 1% DMM); default from file or 0" << std::endl; } enum { @@ -87,7 +91,8 @@ enum { OPT_SCALING_OUTPUT, OPT_SINGLE_PASS_ROTATION, OPT_REDO_ROTATION_SPOTS, - OPT_FORCE_ROTATION_LATTICE + OPT_FORCE_ROTATION_LATTICE, + OPT_BANDWIDTH }; static option long_options[] = { @@ -123,6 +128,7 @@ static option long_options[] = { {"scaling-iterations", required_argument, nullptr, OPT_SCALING_ITERATIONS}, {"scaling-high-resolution", required_argument, nullptr, OPT_SCALING_HIGH_RESOLUTION}, {"scaling-output", required_argument, nullptr, OPT_SCALING_OUTPUT}, + {"bandwidth", required_argument, nullptr, OPT_BANDWIDTH}, {nullptr, 0, nullptr, 0} }; @@ -293,6 +299,8 @@ int main(int argc, char **argv) { int64_t scaling_iter = 3; std::optional forced_rotation_lattice; + std::optional bandwidth_fwhm; // relative FWHM of dlambda/lambda + IndexingAlgorithmEnum indexing_algorithm = IndexingAlgorithmEnum::Auto; GeomRefinementAlgorithmEnum refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter; @@ -410,6 +418,8 @@ int main(int argc, char **argv) { refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter; else if (alg == "orientation") refinement_algorithm = GeomRefinementAlgorithmEnum::OrientationOnly; + else if (alg == "pixelrefine") + refinement_algorithm = GeomRefinementAlgorithmEnum::PixelRefine; else { logger.Error("Invalid geom refinement algorithm: {}", alg); print_usage(); @@ -515,6 +525,13 @@ int main(int argc, char **argv) { exit(EXIT_FAILURE); } break; + case OPT_BANDWIDTH: + bandwidth_fwhm = atof(optarg); + if (!(bandwidth_fwhm.value() >= 0.0f)) { + logger.Error("Invalid bandwidth: {}", optarg); + exit(EXIT_FAILURE); + } + break; default: print_usage(); @@ -616,6 +633,19 @@ int main(int argc, char **argv) { logger.Info("Max spot count overridden to {}", max_spot_count_override.value()); } + // X-ray bandwidth: CLI overrides the value carried in the dataset; otherwise + // keep whatever the dataset provided (0 / none -> monochromatic). + if (bandwidth_fwhm) + experiment.BandwidthFWHM(bandwidth_fwhm); + if (experiment.GetBandwidthFWHM()) + logger.Info("X-ray bandwidth FWHM set to {:.4f}", experiment.GetBandwidthFWHM().value()); + + // PixelRefine integration needs reference intensities (the I_true hypothesis). + if (refinement_algorithm == GeomRefinementAlgorithmEnum::PixelRefine && reference_data.empty()) { + logger.Warning("-r pixelrefine needs --reference-mtz; falling back to beam_and_lattice"); + refinement_algorithm = GeomRefinementAlgorithmEnum::BeamCenter; + } + // Configure Indexing IndexingSettings indexing_settings; indexing_settings.Algorithm(indexing_algorithm);