Files
Jungfraujoch/image_analysis/bragg_integration/BraggPrediction.cpp
2025-09-11 14:37:16 +02:00

60 lines
2.3 KiB
C++

// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include "BraggPrediction.h"
std::multimap<float, Reflection> CalcBraggPredictions(const DiffractionExperiment& experiment,
const CrystalLattice& lattice,
float high_res_A,
float ewald_dist_cutoff,
int max_hkl) {
std::multimap<float, Reflection> reflections;
auto geom = experiment.GetDiffractionGeometry();
auto det_width_pxl = static_cast<float>(experiment.GetXPixelsNum());
auto det_height_pxl = static_cast<float>(experiment.GetYPixelsNum());
float one_over_dmax = 1.0f / high_res_A;
float one_over_dmax_sq = one_over_dmax * one_over_dmax;
Coord Astar = lattice.Astar();
Coord Bstar = lattice.Bstar();
Coord Cstar = lattice.Cstar();
for (int h = -max_hkl; h < max_hkl; h++) {
for (int k = -max_hkl; k < max_hkl; k++) {
for (int l = -max_hkl; l < max_hkl; l++) {
// (000) is not too interesting
if (h == 0 && k == 0 && l == 0)
continue;
Coord recip = Astar * static_cast<float>(h) + Bstar * static_cast<float>(k) + Cstar * static_cast<float>(l);
if (recip * recip > one_over_dmax_sq)
continue;
float d = 1.0f / recip.Length();
float dist_ewald_sphere = std::fabs(geom.DistFromEwaldSphere(recip));
if (!std::isfinite(dist_ewald_sphere) || dist_ewald_sphere > ewald_dist_cutoff )
continue;
auto [x,y] = geom.RecipToDector(recip);
if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl))
continue;
reflections.insert(std::make_pair(dist_ewald_sphere, Reflection{
.h = h,
.k = k,
.l = l,
.predicted_x = x,
.predicted_y = y,
.d = d
}));
}
}
}
return reflections;
}