24#include <boost/iterator/iterator_facade.hpp>
25#include <boost/math/special_functions/round.hpp>
47const std::string PEAKS{
"Peaks"};
48const std::string HOFFSET{
"Hoffset"};
49const std::string KOFFSET{
"Koffset"};
50const std::string LOFFSET{
"Loffset"};
51const std::string INCLUDEPEAKSINRANGE{
"IncludeAllPeaksInRange"};
52const std::string HMIN{
"Hmin"};
53const std::string HMAX{
"Hmax"};
54const std::string KMIN{
"Kmin"};
55const std::string KMAX{
"Kmax"};
56const std::string LMIN{
"Lmin"};
57const std::string LMAX{
"Lmax"};
58const std::string REFLECTION_COND{
"ReflectionCondition"};
59const std::string ON_DETECTOR{
"RequirePeaksOnDetector"};
60const std::string FRACPEAKS{
"FracPeaks"};
67class PeaksInRangeStrategy {
77 : m_hklGenerator(hklMin, hklMax), m_hklIterator(m_hklGenerator.begin()), m_hklFilter(filter),
78 m_inputPeaks(inputPeaks) {
83 return Progress(alg, 0.0, 1.0, m_hklGenerator.size());
86 void initialHKL(
V3D *hkl, DblMatrix *gonioMatrix,
int *runNumber)
const noexcept {
87 *hkl = *(m_hklGenerator.begin());
88 *gonioMatrix = m_inputPeaks->run().getGoniometer().getR();
89 *runNumber = m_inputPeaks->getPeak(0).getRunNumber();
93 bool nextHKL(
V3D *hkl, DblMatrix *gonioMatrix,
int *runNumber)
noexcept {
95 if (m_hklIterator != m_hklGenerator.end()) {
96 if (!m_hklFilter->isAllowed(*m_hklIterator)) {
97 return nextHKL(hkl, gonioMatrix, runNumber);
99 *hkl = *m_hklIterator;
108 HKLGenerator::const_iterator m_hklIterator;
118class PeaksFromIndexedStrategy {
120 explicit PeaksFromIndexedStrategy(
const PeaksWorkspace *
const inputPeaks)
121 : m_inputPeaks{inputPeaks}, m_currentPeak{0} {}
124 return Progress(alg, 0.0, 1.0, m_inputPeaks->getNumberPeaks());
127 void initialHKL(
V3D *hkl, DblMatrix *gonioMatrix,
int *runNumber)
const noexcept {
128 const auto &initialPeak = m_inputPeaks->getPeak(m_currentPeak);
129 *hkl = initialPeak.getHKL();
130 *gonioMatrix = initialPeak.getGoniometerMatrix();
131 *runNumber = initialPeak.getRunNumber();
135 bool nextHKL(
V3D *hkl, DblMatrix *gonioMatrix,
int *runNumber)
noexcept {
136 bool canContinue{
true};
138 if (m_currentPeak < m_inputPeaks->getNumberPeaks()) {
139 const auto &initialPeak = m_inputPeaks->getPeak(m_currentPeak);
140 *hkl = initialPeak.getHKL();
141 *gonioMatrix = initialPeak.getGoniometerMatrix();
142 *runNumber = initialPeak.getRunNumber();
161std::shared_ptr<Mantid::API::IPeaksWorkspace>
164 auto outPeaks = WorkspaceFactory::Instance().createPeaks();
168 lattice->setMaxOrder(modulationProps.
maxOrder);
169 lattice->setCrossTerm(modulationProps.
crossTerms);
170 const auto &offsets = modulationProps.
offsets;
173 for (
const auto &offset : offsets) {
174 const V3D &modVector{std::get<3>(offset)};
175 const double m{std::get<0>(offset)},
n{std::get<1>(offset)}, p{std::get<2>(offset)};
177 if (m == 1 &&
n == 0 && p == 0)
179 else if (m == 0 &&
n == 1 && p == 0)
181 else if (m == 0 &&
n == 0 && p == 1)
185 lattice->setModVec1(modVector);
188 lattice->setModVec2(modVector);
191 lattice->setModVec3(modVector);
195 outPeaks->mutableSample().setOrientedLattice(std::move(lattice));
213template <
typename SearchStrategy>
215predictFractionalPeaks(
Algorithm *
const alg,
const bool requirePeaksOnDetector,
const PeaksWorkspace &inputPeaks,
219 const InstrumentRayTracer tracer(inputPeaks.
getInstrument());
221 const auto &offsets = modulationProps.
offsets;
222 auto outPeaks = createOutputWorkspace(inputPeaks, modulationProps);
223 using PeakHash = std::array<int, 4>;
224 std::vector<PeakHash> alreadyDonePeaks;
230 searchStrategy.initialHKL(¤tHKL, &gonioMatrix, &runNumber);
231 auto progressReporter = searchStrategy.createProgressReporter(alg);
233 for (
const auto &mnpOffset : offsets) {
234 const V3D currentIntHKL{std::round(currentHKL[0]), std::round(currentHKL[1]), std::round(currentHKL[2])};
235 const V3D candidateHKL{currentIntHKL + std::get<3>(mnpOffset)};
236 const V3D qLab = (gonioMatrix * UB * candidateHKL) * 2 * M_PI * qConvention;
249 peak->setGoniometerMatrix(gonioMatrix);
250 if (requirePeaksOnDetector && peak->getDetectorID() < 0)
253 PeakHash savedPeak{runNumber, boost::math::iround(1000.0 * candidateHKL[0]),
254 boost::math::iround(1000.0 * candidateHKL[1]), boost::math::iround(1000.0 * candidateHKL[2])};
256 auto it = find(alreadyDonePeaks.begin(), alreadyDonePeaks.end(), savedPeak);
257 if (it == alreadyDonePeaks.end())
258 alreadyDonePeaks.emplace_back(savedPeak);
262 peak->setHKL(candidateHKL * qConvention);
263 peak->setIntHKL(currentIntHKL * qConvention);
264 const
double m{std::get<0>(mnpOffset)},
n{std::get<1>(mnpOffset)}, p{std::get<2>(mnpOffset)};
266 peak->setIntMNP(
V3D(m,
n, p));
267 peak->setRunNumber(runNumber);
268 outPeaks->addPeak(*peak);
270 progressReporter.report();
271 if (!searchStrategy.nextHKL(¤tHKL, &gonioMatrix, &runNumber))
293 "Workspace of Peaks with orientation matrix that indexed the peaks and "
294 "instrument loaded");
297 "Offset in the h direction");
299 "Offset in the k direction");
301 "Offset in the h direction");
302 declareProperty(PropertyNames::INCLUDEPEAKSINRANGE,
false,
"If false only offsets from peaks from Peaks are used");
303 declareProperty(PropertyNames::HMIN, -8.0,
"Minimum H value to use during search",
Direction::Input);
304 declareProperty(PropertyNames::HMAX, 8.0,
"Maximum H value to use during search",
Direction::Input);
305 declareProperty(PropertyNames::KMIN, -8.0,
"Minimum K value to use during search",
Direction::Input);
306 declareProperty(PropertyNames::KMAX, 8.0,
"Maximum K value to use during search",
Direction::Input);
307 declareProperty(PropertyNames::LMIN, -8.0,
"Minimum L value to use during search",
Direction::Input);
308 declareProperty(PropertyNames::LMAX, 8.0,
"Maximum L value to use during search",
Direction::Input);
311 std::vector<std::string> propOptions;
312 propOptions.reserve(reflectionConditions.size() + 1);
313 propOptions.emplace_back(
"");
314 std::transform(reflectionConditions.cbegin(), reflectionConditions.cend(), std::back_inserter(propOptions),
315 [](
const auto &condition) { return condition->getName(); });
316 declareProperty(PropertyNames::REFLECTION_COND,
"", std::make_shared<StringListValidator>(propOptions),
317 "If provided, generate a list of possible peaks from this "
318 "reflection condition and use them to predict the fractional "
319 "peaks. This option requires a range of HKL values and "
320 "implies IncludeAllPeaksInRange=true");
322 declareProperty(PropertyNames::ON_DETECTOR,
true,
323 "If true then the predicted peaks are required to hit a "
324 "detector pixel. Default=true",
330 for (
const auto &name : {PropertyNames::HMIN, PropertyNames::HMAX, PropertyNames::KMIN, PropertyNames::KMAX,
331 PropertyNames::LMIN, PropertyNames::LMAX}) {
334 setPropertySettings(name, std::make_unique<Kernel::EnabledWhenProperty>(
335 std::move(includeInRangeEqOne), std::move(reflConditionNotEmpty),
Kernel::OR));
338 for (
const auto &name : {PropertyNames::HOFFSET, PropertyNames::KOFFSET, PropertyNames::LOFFSET}) {
339 setPropertyGroup(name,
"Separate Offsets");
341 for (
const auto &name :
344 setPropertyGroup(name,
"Modulation Vectors");
347 for (
const auto &offsetName : {PropertyNames::HOFFSET, PropertyNames::KOFFSET, PropertyNames::LOFFSET}) {
352 setPropertySettings(offsetName,
353 std::make_unique<Kernel::EnabledWhenProperty>(std::move(modVectorOneAndTwoIsDefault),
359 "Workspace of Peaks with peaks with fractional h,k, and/or l values");
367 std::map<std::string, std::string> helpMessages;
369 if (peaks && peaks->getNumberPeaks() <= 0) {
370 helpMessages[PropertyNames::PEAKS] =
"Input workspace has no peaks.";
373 auto validateRange = [&helpMessages,
this](
const std::string &minName,
const std::string &maxName) {
376 const std::string helpMsg =
"Inconsistent " + minName +
"/" + maxName +
": " + maxName +
" < " + minName;
377 helpMessages[minName] = helpMsg;
378 helpMessages[maxName] = helpMsg;
381 validateRange(PropertyNames::HMIN, PropertyNames::HMAX);
382 validateRange(PropertyNames::KMIN, PropertyNames::KMAX);
383 validateRange(PropertyNames::LMIN, PropertyNames::LMAX);
400 const bool includePeaksInRange =
getProperty(
"IncludeAllPeaksInRange");
405 const std::string reflectionConditionName =
getProperty(PropertyNames::REFLECTION_COND);
406 const bool requirePeakOnDetector =
getProperty(PropertyNames::ON_DETECTOR);
409 if (includePeaksInRange || !reflectionConditionName.empty()) {
412 const auto found = std::find_if(
413 std::cbegin(allConditions), std::cend(allConditions),
414 [&reflectionConditionName](
const auto &condition) {
return condition->getName() == reflectionConditionName; });
417 if (found != std::cend(allConditions)) {
418 filter = std::make_unique<HKLFilterCentering>(*found);
421 filter = std::make_unique<HKLFilterNone>();
423 outPeaks = predictFractionalPeaks(
this, requirePeakOnDetector, *inputPeaks, modulationInfo,
424 PeaksInRangeStrategy(hklMin, hklMax, filter.get(), inputPeaks.get()));
426 outPeaks = predictFractionalPeaks(
this, requirePeakOnDetector, *inputPeaks, modulationInfo,
427 PeaksFromIndexedStrategy(inputPeaks.get()));
443 std::vector<double> hOffsets =
getProperty(PropertyNames::HOFFSET);
444 std::vector<double> kOffsets =
getProperty(PropertyNames::KOFFSET);
445 std::vector<double> lOffsets =
getProperty(PropertyNames::LOFFSET);
446 if (hOffsets.empty())
447 hOffsets.emplace_back(0.0);
448 if (kOffsets.empty())
449 kOffsets.emplace_back(0.0);
450 if (lOffsets.empty())
451 lOffsets.emplace_back(0.0);
452 const bool crossTerms{
false}, saveOnLattice{
false};
#define DECLARE_ALGORITHM(classname)
const std::vector< V3D > modVectors
#define GNU_DIAG_OFF(x)
This is a collection of macros for turning compiler warnings off in a controlled manner.
Base class from which all concrete algorithm classes should be derived.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Geometry::Instrument_const_sptr getInstrument() const
Returns the parameterized instrument.
const Sample & sample() const
Sample accessors.
Helper class for reporting progress from algorithms.
const Geometry::OrientedLattice & getOrientedLattice() const
Get a reference to the sample's OrientedLattice.
A property class for workspaces.
Using a set of offset vectors, either provided as separate lists or as a set of vectors,...
void exec() override
Execute the algorithm.
std::map< std::string, std::string > validateInputs() override
Validate the input once all values are set.
ModulationProperties getModulationInfo()
Structure describing a single-crystal peak.
The class PeaksWorkspace stores information about a set of SCD peaks.
IPeak_uptr createPeak(const Kernel::V3D &QLabFrame, boost::optional< double > detectorDistance=boost::none) const override
Creates an instance of a Peak BUT DOES NOT ADD IT TO THE WORKSPACE.
A class to filter HKLs according to a lattice centering.
This class is responsible for tracking rays and accumulating a list of objects that are intersected a...
Class to implement UB matrix.
const Kernel::DblMatrix & getUB() const
Get the UB matrix.
Support for a property that holds an array of values.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
Manage the lifetime of a class intended to be a singleton.
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
std::vector< MNPOffset > generateOffsetVectors(const std::vector< Kernel::V3D > &modVectors, const int maxOrder, const bool crossTerms)
Calculate a list of HKL offsets from the given modulation vectors.
double qConventionFactor()
convenience overload to pull the convention from the config service
std::vector< Kernel::V3D > validModulationVectors(const std::vector< double > &modVector1, const std::vector< double > &modVector2, const std::vector< double > &modVector3)
Create a list of valid modulation vectors from the input.
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::unique_ptr< Peak > Peak_uptr
std::shared_ptr< ReflectionCondition > ReflectionCondition_sptr
Shared pointer to a ReflectionCondition.
std::unique_ptr< HKLFilter > HKLFilter_uptr
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
MANTID_GEOMETRY_DLL const ReflectionConditions & getAllReflectionConditions()
std::unique_ptr< IPeak > IPeak_uptr
Mantid::Kernel::Matrix< double > DblMatrix
Tie together the names of the properties for the modulation vectors.
static const std::string ModVector3
static const std::string MaxOrder
static const std::string ModVector2
std::vector< MNPOffset > offsets
static ModulationProperties create(const API::IAlgorithm &alg)
Create a ModulationProperties object from an algorithm.
static const std::string CrossTerms
static const std::string ModVector1
static void appendTo(API::IAlgorithm *alg)
Append the common set of properties that relate to modulation vectors to the given algorithm.
Describes the direction (within an algorithm) of a Property.
@ Input
An input workspace.
@ Output
An output workspace.