Mantid
Loading...
Searching...
No Matches
PredictFractionalPeaks.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8#include "MantidAPI/Run.h"
9#include "MantidAPI/Sample.h"
23
24#include <boost/iterator/iterator_facade.hpp>
25#include <boost/math/special_functions/round.hpp>
26
43
44namespace {
45
46namespace PropertyNames {
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"};
61} // namespace PropertyNames
62
67class PeaksInRangeStrategy {
68public:
76 PeaksInRangeStrategy(const V3D &hklMin, const V3D &hklMax, HKLFilter *filter, const PeaksWorkspace *const inputPeaks)
77 : m_hklGenerator(hklMin, hklMax), m_hklIterator(m_hklGenerator.begin()), m_hklFilter(filter),
78 m_inputPeaks(inputPeaks) {
79 assert(filter);
80 }
81
82 Progress createProgressReporter(Algorithm *const alg) const noexcept {
83 return Progress(alg, 0.0, 1.0, m_hklGenerator.size());
84 }
85
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();
90 }
91
93 bool nextHKL(V3D *hkl, DblMatrix *gonioMatrix, int *runNumber) noexcept {
94 ++m_hklIterator;
95 if (m_hklIterator != m_hklGenerator.end()) {
96 if (!m_hklFilter->isAllowed(*m_hklIterator)) {
97 return nextHKL(hkl, gonioMatrix, runNumber);
98 } else {
99 *hkl = *m_hklIterator;
100 return true;
101 }
102 } else
103 return false;
104 }
105
106private:
107 HKLGenerator m_hklGenerator;
108 HKLGenerator::const_iterator m_hklIterator;
109 HKLFilter *m_hklFilter;
110 const PeaksWorkspace *const m_inputPeaks;
111}; // namespace
112
118class PeaksFromIndexedStrategy {
119public:
120 explicit PeaksFromIndexedStrategy(const PeaksWorkspace *const inputPeaks)
121 : m_inputPeaks{inputPeaks}, m_currentPeak{0} {}
122
123 Progress createProgressReporter(Algorithm *const alg) const noexcept {
124 return Progress(alg, 0.0, 1.0, m_inputPeaks->getNumberPeaks());
125 }
126
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();
132 }
133
135 bool nextHKL(V3D *hkl, DblMatrix *gonioMatrix, int *runNumber) noexcept {
136 bool canContinue{true};
137 m_currentPeak++;
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();
143 } else {
144 canContinue = false;
145 }
146
147 return canContinue;
148 }
149
150private:
151 const PeaksWorkspace *const m_inputPeaks;
152 int m_currentPeak;
153};
154
161std::shared_ptr<Mantid::API::IPeaksWorkspace>
162createOutputWorkspace(const PeaksWorkspace &inputPeaks, const Mantid::Crystal::ModulationProperties &modulationProps) {
164 auto outPeaks = WorkspaceFactory::Instance().createPeaks();
165 outPeaks->setInstrument(inputPeaks.getInstrument());
166 if (modulationProps.saveOnLattice) {
167 auto lattice = std::make_unique<Mantid::Geometry::OrientedLattice>(inputPeaks.sample().getOrientedLattice());
168 lattice->setMaxOrder(modulationProps.maxOrder);
169 lattice->setCrossTerm(modulationProps.crossTerms);
170 const auto &offsets = modulationProps.offsets;
171 // there should be a maximum of 3 modulation vectors. Store the
172 // order=(1,0,0),(0,1,0), (0,0,1) vectors
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)};
176 int modNum(-1);
177 if (m == 1 && n == 0 && p == 0)
178 modNum = 1;
179 else if (m == 0 && n == 1 && p == 0)
180 modNum = 2;
181 else if (m == 0 && n == 0 && p == 1)
182 modNum = 3;
183 switch (modNum) {
184 case 1:
185 lattice->setModVec1(modVector);
186 break;
187 case 2:
188 lattice->setModVec2(modVector);
189 break;
190 case 3:
191 lattice->setModVec3(modVector);
192 break;
193 }
194 }
195 outPeaks->mutableSample().setOrientedLattice(std::move(lattice));
196 }
197 return outPeaks;
198}
199
213template <typename SearchStrategy>
215predictFractionalPeaks(Algorithm *const alg, const bool requirePeaksOnDetector, const PeaksWorkspace &inputPeaks,
216 const Mantid::Crystal::ModulationProperties &modulationProps, SearchStrategy searchStrategy) {
218
219 const InstrumentRayTracer tracer(inputPeaks.getInstrument());
220 const auto &UB = inputPeaks.sample().getOrientedLattice().getUB();
221 const auto &offsets = modulationProps.offsets;
222 auto outPeaks = createOutputWorkspace(inputPeaks, modulationProps);
223 using PeakHash = std::array<int, 4>;
224 std::vector<PeakHash> alreadyDonePeaks;
225
226 const auto qConvention{Mantid::Crystal::qConventionFactor()};
227 V3D currentHKL;
228 DblMatrix gonioMatrix;
229 int runNumber{0};
230 searchStrategy.initialHKL(&currentHKL, &gonioMatrix, &runNumber);
231 auto progressReporter = searchStrategy.createProgressReporter(alg);
232 while (true) {
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;
237 if (qLab[2] <= 0)
238 continue;
239
240 IPeak_uptr ipeak;
241 try {
242 ipeak = inputPeaks.createPeak(qLab);
243 } catch (...) {
244 // If we can't create a valid peak we have no choice but to skip
245 // it
246 continue;
247 }
248 Peak_uptr peak(static_cast<Peak *>(ipeak.release()));
249 peak->setGoniometerMatrix(gonioMatrix);
250 if (requirePeaksOnDetector && peak->getDetectorID() < 0)
251 continue;
252 GNU_DIAG_OFF("missing-braces")
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])};
255 GNU_DIAG_ON("missing-braces")
256 auto it = find(alreadyDonePeaks.begin(), alreadyDonePeaks.end(), savedPeak);
257 if (it == alreadyDonePeaks.end())
258 alreadyDonePeaks.emplace_back(savedPeak);
259 else
260 continue;
261
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)};
265 if (fabs(m) > 0. || fabs(n) > 0. || fabs(p) > 0.)
266 peak->setIntMNP(V3D(m, n, p));
267 peak->setRunNumber(runNumber);
268 outPeaks->addPeak(*peak);
269 }
270 progressReporter.report();
271 if (!searchStrategy.nextHKL(&currentHKL, &gonioMatrix, &runNumber))
272 break;
273 }
274
275 return outPeaks;
276}
277
278} // namespace
279
280namespace Mantid::Crystal {
281
282DECLARE_ALGORITHM(PredictFractionalPeaks)
283
284
288 using Kernel::Direction;
290
291 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace_sptr::element_type>>(PropertyNames::PEAKS, "",
293 "Workspace of Peaks with orientation matrix that indexed the peaks and "
294 "instrument loaded");
295
296 declareProperty(std::make_unique<Kernel::ArrayProperty<double>>(PropertyNames::HOFFSET, "-0.5,0.0,0.5"),
297 "Offset in the h direction");
298 declareProperty(std::make_unique<Kernel::ArrayProperty<double>>(PropertyNames::KOFFSET, "0"),
299 "Offset in the k direction");
300 declareProperty(std::make_unique<Kernel::ArrayProperty<double>>(PropertyNames::LOFFSET, "-0.5,0.5"),
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);
309
310 const auto &reflectionConditions = getAllReflectionConditions();
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");
321
322 declareProperty(PropertyNames::ON_DETECTOR, true,
323 "If true then the predicted peaks are required to hit a "
324 "detector pixel. Default=true",
327
328 // enable range properties if required
330 for (const auto &name : {PropertyNames::HMIN, PropertyNames::HMAX, PropertyNames::KMIN, PropertyNames::KMAX,
331 PropertyNames::LMIN, PropertyNames::LMAX}) {
332 EnabledWhenProperty includeInRangeEqOne{PropertyNames::INCLUDEPEAKSINRANGE, Kernel::IS_EQUAL_TO, "1"};
333 EnabledWhenProperty reflConditionNotEmpty{PropertyNames::REFLECTION_COND, Kernel::IS_NOT_EQUAL_TO, ""};
334 setPropertySettings(name, std::make_unique<Kernel::EnabledWhenProperty>(
335 std::move(includeInRangeEqOne), std::move(reflConditionNotEmpty), Kernel::OR));
336 }
337 // group offset/modulations options together
338 for (const auto &name : {PropertyNames::HOFFSET, PropertyNames::KOFFSET, PropertyNames::LOFFSET}) {
339 setPropertyGroup(name, "Separate Offsets");
340 }
341 for (const auto &name :
344 setPropertyGroup(name, "Modulation Vectors");
345 }
346 // enable/disable offsets & modulation vectors appropriately
347 for (const auto &offsetName : {PropertyNames::HOFFSET, PropertyNames::KOFFSET, PropertyNames::LOFFSET}) {
351 EnabledWhenProperty modVectorOneAndTwoIsDefault{modVectorOneIsDefault, modVectorTwoIsDefault, Kernel::AND};
352 setPropertySettings(offsetName,
353 std::make_unique<Kernel::EnabledWhenProperty>(std::move(modVectorOneAndTwoIsDefault),
354 std::move(modVectorThreeIsDefault), Kernel::AND));
355 }
356 // Outputs
357 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace_sptr::element_type>>(PropertyNames::FRACPEAKS, "",
359 "Workspace of Peaks with peaks with fractional h,k, and/or l values");
360}
361
366std::map<std::string, std::string> PredictFractionalPeaks::validateInputs() {
367 std::map<std::string, std::string> helpMessages;
368 const PeaksWorkspace_sptr peaks = getProperty(PropertyNames::PEAKS);
369 if (peaks && peaks->getNumberPeaks() <= 0) {
370 helpMessages[PropertyNames::PEAKS] = "Input workspace has no peaks.";
371 }
372
373 auto validateRange = [&helpMessages, this](const std::string &minName, const std::string &maxName) {
374 const double min{getProperty(minName)}, max{getProperty(maxName)};
375 if (max < min) {
376 const std::string helpMsg = "Inconsistent " + minName + "/" + maxName + ": " + maxName + " < " + minName;
377 helpMessages[minName] = helpMsg;
378 helpMessages[maxName] = helpMsg;
379 }
380 };
381 validateRange(PropertyNames::HMIN, PropertyNames::HMAX);
382 validateRange(PropertyNames::KMIN, PropertyNames::KMAX);
383 validateRange(PropertyNames::LMIN, PropertyNames::LMAX);
384
385 // If a modulation vector is provided then maxOrder is needed
390 if (maxOrder == 0 && !modVectors.empty()) {
391 helpMessages[ModulationProperties::MaxOrder] = "Maxorder required when specifying a modulation vector.";
392 }
393 return helpMessages;
394}
395
398 PeaksWorkspace_sptr inputPeaks = getProperty(PropertyNames::PEAKS);
399 auto modulationInfo = getModulationInfo();
400 const bool includePeaksInRange = getProperty("IncludeAllPeaksInRange");
401 const V3D hklMin{getProperty(PropertyNames::HMIN), getProperty(PropertyNames::KMIN),
402 getProperty(PropertyNames::LMIN)};
403 const V3D hklMax{getProperty(PropertyNames::HMAX), getProperty(PropertyNames::KMAX),
404 getProperty(PropertyNames::LMAX)};
405 const std::string reflectionConditionName = getProperty(PropertyNames::REFLECTION_COND);
406 const bool requirePeakOnDetector = getProperty(PropertyNames::ON_DETECTOR);
407
408 IPeaksWorkspace_sptr outPeaks;
409 if (includePeaksInRange || !reflectionConditionName.empty()) {
411 const auto &allConditions = getAllReflectionConditions();
412 const auto found = std::find_if(
413 std::cbegin(allConditions), std::cend(allConditions),
414 [&reflectionConditionName](const auto &condition) { return condition->getName() == reflectionConditionName; });
416 HKLFilter_uptr filter;
417 if (found != std::cend(allConditions)) {
418 filter = std::make_unique<HKLFilterCentering>(*found);
419 } else {
421 filter = std::make_unique<HKLFilterNone>();
422 }
423 outPeaks = predictFractionalPeaks(this, requirePeakOnDetector, *inputPeaks, modulationInfo,
424 PeaksInRangeStrategy(hklMin, hklMax, filter.get(), inputPeaks.get()));
425 } else {
426 outPeaks = predictFractionalPeaks(this, requirePeakOnDetector, *inputPeaks, modulationInfo,
427 PeaksFromIndexedStrategy(inputPeaks.get()));
428 }
429 setProperty(PropertyNames::FRACPEAKS, outPeaks);
430}
431
438 // Input validation ensures that we have either specified offests or
439 // modulation vectors
441
442 if (maxOrder == 0) {
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};
453 return {generateOffsetVectors(hOffsets, kOffsets, lOffsets), maxOrder, crossTerms, saveOnLattice};
454 } else {
455 return ModulationProperties::create(*this);
456 }
457}
458
459} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< V3D > modVectors
Definition: IndexPeaks.cpp:56
const bool crossTerms
Definition: IndexPeaks.cpp:57
const int maxOrder
Definition: IndexPeaks.cpp:55
#define fabs(x)
Definition: Matrix.cpp:22
#define GNU_DIAG_ON(x)
#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.
Definition: Algorithm.h:85
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Geometry::Instrument_const_sptr getInstrument() const
Returns the parameterized instrument.
const Sample & sample() const
Sample accessors.
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
const Geometry::OrientedLattice & getOrientedLattice() const
Get a reference to the sample's OrientedLattice.
Definition: Sample.cpp:154
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.
Structure describing a single-crystal peak.
Definition: Peak.h:34
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.
Definition: ArrayProperty.h:28
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...
Definition: ListValidator.h:29
Manage the lifetime of a class intended to be a singleton.
Class for 3D vectors.
Definition: V3D.h:34
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
Definition: Peak.h:170
std::shared_ptr< ReflectionCondition > ReflectionCondition_sptr
Shared pointer to a ReflectionCondition.
std::unique_ptr< HKLFilter > HKLFilter_uptr
Definition: HKLFilter.h:75
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
Definition: IPeak.h:103
Mantid::Kernel::Matrix< double > DblMatrix
Definition: Matrix.h:206
Tie together the names of the properties for the modulation vectors.
static ModulationProperties create(const API::IAlgorithm &alg)
Create a ModulationProperties object from an algorithm.
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.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54