Mantid
Loading...
Searching...
No Matches
PredictPeaks.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 +
10#include "MantidAPI/Run.h"
11#include "MantidAPI/Sample.h"
29
30#include <fstream>
32
33namespace Mantid::Crystal {
34
35// Register the algorithm into the AlgorithmFactory
36DECLARE_ALGORITHM(PredictPeaks)
37
38using namespace Mantid::API;
39using namespace Mantid::DataObjects;
40using namespace Mantid::Geometry;
41using namespace Mantid::Kernel;
42
46 : m_refConds(getAllReflectionConditions()), m_runNumber(-1), m_inst(), m_pw(), m_sfCalculator(),
47 m_qConventionFactor(qConventionFactor()) {}
48
52 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("InputWorkspace", "", Direction::Input),
53 "An input workspace (MatrixWorkspace, MDEventWorkspace, or "
54 "PeaksWorkspace) containing:\n"
55 " - The relevant Instrument (calibrated as needed).\n"
56 " - A sample with a UB matrix.\n"
57 " - The goniometer rotation matrix.");
58
59 declareProperty(std::make_unique<PropertyWithValue<double>>("WavelengthMin", 0.1, Direction::Input),
60 "Minimum wavelength limit at which to start looking for single-crystal "
61 "peaks.");
62 declareProperty(std::make_unique<PropertyWithValue<double>>("WavelengthMax", 100.0, Direction::Input),
63 "Maximum wavelength limit at which to stop looking for single-crystal "
64 "peaks.");
65
66 declareProperty(std::make_unique<PropertyWithValue<double>>("MinDSpacing", 1.0, Direction::Input),
67 "Minimum d-spacing of peaks to consider. Default = 1.0");
68 declareProperty(std::make_unique<PropertyWithValue<double>>("MaxDSpacing", 100.0, Direction::Input),
69 "Maximum d-spacing of peaks to consider.");
70
71 declareProperty("CalculateGoniometerForCW", false,
72 "This will calculate the goniometer rotation (around y-axis "
73 "only) for a constant wavelength.");
74
75 auto nonNegativeDbl = std::make_shared<BoundedValidator<double>>();
76 nonNegativeDbl->setLower(0);
77 declareProperty("Wavelength", DBL_MAX, nonNegativeDbl, "Wavelength to use when calculating goniometer angle");
78 setPropertySettings("Wavelength", std::make_unique<EnabledWhenProperty>("CalculateGoniometerForCW", IS_NOT_DEFAULT));
79
80 declareProperty("InnerGoniometer", false,
81 "Whether the goniometer to be calculated is the most inner "
82 "(phi) or most outer (omega)");
83 setPropertySettings("InnerGoniometer",
84 std::make_unique<EnabledWhenProperty>("CalculateGoniometerForCW",
86
87 declareProperty("FlipX", false,
88 "Used when calculating goniometer angle if the q_lab x value "
89 "should be negative, hence the detector of the other side "
90 "(right) of the beam");
91 setPropertySettings("FlipX", std::make_unique<EnabledWhenProperty>(
92 "CalculateGoniometerForCW", Mantid::Kernel::ePropertyCriterion::IS_NOT_DEFAULT));
93
94 declareProperty(std::make_unique<PropertyWithValue<double>>("MinAngle", -180, Direction::Input),
95 "Minimum goniometer rotation angle");
96 setPropertySettings("MinAngle", std::make_unique<EnabledWhenProperty>("CalculateGoniometerForCW", IS_NOT_DEFAULT));
97
98 declareProperty(std::make_unique<PropertyWithValue<double>>("MaxAngle", 180, Direction::Input),
99 "Maximum goniometer rotation angle");
100 setPropertySettings("MaxAngle", std::make_unique<EnabledWhenProperty>("CalculateGoniometerForCW", IS_NOT_DEFAULT));
101
102 // Build up a list of reflection conditions to use
103 std::vector<std::string> propOptions;
104 propOptions.reserve(m_refConds.size());
105 std::transform(m_refConds.cbegin(), m_refConds.cend(), std::back_inserter(propOptions),
106 [](const auto &condition) { return condition->getName(); });
107 declareProperty("ReflectionCondition", "Primitive", std::make_shared<StringListValidator>(propOptions),
108 "Which reflection condition applies to this crystal, "
109 "reducing the number of expected HKL peaks?");
110
111 declareProperty("CalculateStructureFactors", false,
112 "Calculate structure factors for the predicted peaks. This "
113 "option only works if the sample of the input workspace has "
114 "a crystal structure assigned.");
115
116 declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("HKLPeaksWorkspace", "", Direction::Input,
118 "Optional: An input PeaksWorkspace with the HKL of the peaks "
119 "that we should predict. \n"
120 "The WavelengthMin/Max and Min/MaxDSpacing parameters are "
121 "unused if this is specified.");
122
123 declareProperty("RoundHKL", true,
124 "When using HKLPeaksWorkspace, this will round the HKL "
125 "values in the HKLPeaksWorkspace to the nearest integers if "
126 "checked.\n"
127 "Keep unchecked to use the original values");
128
129 setPropertySettings("RoundHKL", std::make_unique<EnabledWhenProperty>("HKLPeaksWorkspace", IS_NOT_DEFAULT));
130
131 // Disable some props when using HKLPeaksWorkspace
132 auto makeSet = [] {
133 std::unique_ptr<IPropertySettings> set = std::make_unique<EnabledWhenProperty>("HKLPeaksWorkspace", IS_DEFAULT);
134 return set;
135 };
136 setPropertySettings("WavelengthMin", makeSet());
137 setPropertySettings("WavelengthMax", makeSet());
138 setPropertySettings("MinDSpacing", makeSet());
139 setPropertySettings("MaxDSpacing", makeSet());
140 setPropertySettings("ReflectionCondition", makeSet());
141
142 std::vector<std::string> peakTypes = {"Peak", "LeanElasticPeak"};
143 declareProperty("OutputType", "Peak", std::make_shared<StringListValidator>(peakTypes),
144 "Type of Peak in OutputWorkspace");
146 "CalculateWavelength", true,
147 "When OutputType is LeanElasticPeak you can choose to calculate the wavelength of the peak using the instrument "
148 "and check it is in the valid range or alternatively just accept every peak while not setting the goniometer "
149 "(Q-lab will be incorrect).");
150 setPropertySettings("CalculateWavelength", std::make_unique<EnabledWhenProperty>("OutputType", IS_NOT_DEFAULT));
151
152 declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
153 "An output PeaksWorkspace.");
154
155 declareProperty("PredictPeaksOutsideDetectors", false,
156 "Use an extended detector space (if defined for the"
157 " instrument) to predict peaks which do not fall onto any"
158 "detector. This may produce a very high number of results.");
159
160 auto nonNegativeInt = std::make_shared<BoundedValidator<int>>();
161 nonNegativeInt->setLower(0);
162 declareProperty("EdgePixels", 0, nonNegativeInt, "Remove peaks that are at pixels this close to edge. ");
163}
164
168 // Get the input properties
169 Workspace_sptr rawInputWorkspace = getProperty("InputWorkspace");
170 m_edge = this->getProperty("EdgePixels");
171 m_leanElasticPeak = (getPropertyValue("OutputType") == "LeanElasticPeak");
172 bool usingInstrument = !(m_leanElasticPeak && !getProperty("CalculateWavelength"));
173
174 ExperimentInfo_sptr inputExperimentInfo = std::dynamic_pointer_cast<ExperimentInfo>(rawInputWorkspace);
175
176 MatrixWorkspace_sptr matrixWS = std::dynamic_pointer_cast<MatrixWorkspace>(rawInputWorkspace);
177 PeaksWorkspace_sptr peaksWS = std::dynamic_pointer_cast<PeaksWorkspace>(rawInputWorkspace);
178 MultipleExperimentInfos_sptr mdWS = std::dynamic_pointer_cast<MultipleExperimentInfos>(rawInputWorkspace);
179
180 std::vector<DblMatrix> gonioVec;
181 if (matrixWS) {
182 // Retrieve the goniometer rotation matrix
183 try {
184 auto goniometerMatrix = matrixWS->run().getGoniometerMatrices();
185 std::move(goniometerMatrix.begin(), goniometerMatrix.end(), back_inserter(gonioVec));
186 } catch (std::runtime_error &e) {
187 // If there is no goniometer matrix, use identity matrix instead.
188 g_log.error() << "Error getting the goniometer rotation matrix from the "
189 "InputWorkspace.\n"
190 << e.what() << '\n';
191 g_log.warning() << "Using identity goniometer rotation matrix instead.\n";
192 }
193 } else if (peaksWS) {
194 // Sort peaks by run number so that peaks with equal goniometer matrices are
195 // adjacent
196 std::vector<std::pair<std::string, bool>> criteria;
197 criteria.emplace_back("RunNumber", true);
198
199 peaksWS->sort(criteria);
200
201 // Get all goniometer matrices
202 DblMatrix lastGoniometerMatrix = Matrix<double>(3, 3, false);
203 for (int i = 0; i < static_cast<int>(peaksWS->getNumberPeaks()); ++i) {
204 IPeak &p = peaksWS->getPeak(i);
205 DblMatrix currentGoniometerMatrix = p.getGoniometerMatrix();
206 if (!(currentGoniometerMatrix == lastGoniometerMatrix)) {
207 gonioVec.emplace_back(currentGoniometerMatrix);
208 lastGoniometerMatrix = currentGoniometerMatrix;
209 }
210 }
211
212 } else if (mdWS) {
213 if (mdWS->getNumExperimentInfo() <= 0)
214 throw std::invalid_argument("Specified a MDEventWorkspace as InputWorkspace but it does not have "
215 "any ExperimentInfo associated. Please choose a workspace with a "
216 "full instrument and sample.");
217
218 inputExperimentInfo = mdWS->getExperimentInfo(0);
219
220 // Retrieve the goniometer rotation matrices for each experiment info
221 for (uint16_t i = 0; i < mdWS->getNumExperimentInfo(); ++i) {
222 try {
223 auto goniometerMatrix = mdWS->getExperimentInfo(i)->mutableRun().getGoniometerMatrices();
224 std::move(goniometerMatrix.begin(), goniometerMatrix.end(), back_inserter(gonioVec));
225 } catch (std::runtime_error &e) {
226 // If there is no goniometer matrix, use identity matrix instead.
227 gonioVec.emplace_back(DblMatrix(3, 3, true));
228
229 g_log.error() << "Error getting the goniometer rotation matrix from the "
230 "InputWorkspace.\n"
231 << e.what() << '\n';
232 g_log.warning() << "Using identity goniometer rotation matrix instead.\n";
233 }
234 }
235 }
236
237 // If there's no goniometer matrix at this point, push back an identity
238 // matrix.
239 if (gonioVec.empty()) {
240 gonioVec.emplace_back(DblMatrix(3, 3, true));
241 }
242
243 if (usingInstrument) {
244 setInstrumentFromInputWorkspace(inputExperimentInfo);
245 setRunNumberFromInputWorkspace(inputExperimentInfo);
248 }
249
250 // Create the output
251 if (m_leanElasticPeak) {
252 m_pw = std::make_shared<LeanElasticPeaksWorkspace>();
253 } else {
254 m_pw = std::make_shared<PeaksWorkspace>();
255 }
256 // Copy instrument, sample, etc.
257 m_pw->copyExperimentInfoFrom(inputExperimentInfo.get());
258
259 const Sample &sample = inputExperimentInfo->sample();
260
261 // Retrieve the OrientedLattice (UnitCell) from the workspace
262 const OrientedLattice &orientedLattice = sample.getOrientedLattice();
263
264 // Get the UB matrix from it
265 Matrix<double> ub(3, 3, true);
266 ub = orientedLattice.getUB();
267
268 std::vector<V3D> possibleHKLs;
269 IPeaksWorkspace_sptr possibleHKLWorkspace = getProperty("HKLPeaksWorkspace");
270
271 if (!possibleHKLWorkspace) {
272 fillPossibleHKLsUsingGenerator(orientedLattice, possibleHKLs);
273 } else {
274 fillPossibleHKLsUsingPeaksWorkspace(possibleHKLWorkspace, possibleHKLs);
275 }
276
278
279 /* The wavelength filtering can not be done before because it depends
280 * on q being correctly oriented, so an additional filtering step is required.
281 */
282 double lambdaMin = getProperty("WavelengthMin");
283 double lambdaMax = getProperty("WavelengthMax");
284
285 Progress prog(this, 0.0, 1.0, possibleHKLs.size() * gonioVec.size());
286 prog.setNotifyStep(0.01);
287
288 if (usingInstrument)
289 m_detectorCacheSearch = std::make_unique<DetectorSearcher>(m_inst, m_pw->detectorInfo());
290
291 if (!usingInstrument) {
292 for (auto &possibleHKL : possibleHKLs) {
294 }
295 } else if (getProperty("CalculateGoniometerForCW")) {
296 size_t allowedPeakCount = 0;
297
298 double wavelength = getProperty("Wavelength");
299 if (wavelength == DBL_MAX) {
300 if (m_inst->hasParameter("wavelength")) {
301 wavelength = m_inst->getNumberParameter("wavelength").at(0);
302 } else {
303 throw std::runtime_error("Could not get wavelength, neither Wavelength algorithm property "
304 "set nor instrument wavelength parameter");
305 }
306 }
307 double angleMin = getProperty("MinAngle");
308 double angleMax = getProperty("MaxAngle");
309 bool innerGoniometer = getProperty("InnerGoniometer");
310 bool flipX = getProperty("FlipX");
311 for (auto &possibleHKL : possibleHKLs) {
312 Geometry::Goniometer goniometer(gonioVec.front());
313 V3D q_sample = ub * possibleHKL * (2.0 * M_PI * m_qConventionFactor);
314 goniometer.calcFromQSampleAndWavelength(q_sample, wavelength, flipX, innerGoniometer);
315 std::vector<double> angles = goniometer.getEulerAngles("YZY");
316 double angle = innerGoniometer ? angles[2] : angles[0];
317 DblMatrix orientedUB = goniometer.getR() * ub;
318
319 // NOTE: use standard formula
320 // qLab = goniometerMatirx * qSample
321 // = goniometerMatirx * (2pi * UB * hkl * signConvention)
322 V3D q_lab = goniometer.getR() * q_sample;
323 // NOTE: use standard formula
324 // 4pi 4pi |Q^lab_z|
325 // lambda = --------- sin(theta) = ----------------
326 // Q^lab |Q^lab|^2
327 double lambda = (4.0 * M_PI * std::abs(q_lab.Z())) / q_lab.norm2();
328
329 if (!std::isfinite(angle) || angle < angleMin || angle > angleMax)
330 continue;
331
332 if (std::abs(wavelength - lambda) < 0.01) {
333 g_log.information() << "Found goniometer rotation to be in YZY convention [" << angles[0] << ", " << angles[1]
334 << ". " << angles[2] << "] degrees for Q sample = " << q_sample << "\n";
335 calculateQAndAddToOutput(possibleHKL, orientedUB, goniometer.getR());
336 ++allowedPeakCount;
337 }
338 prog.report();
339 }
340
341 logNumberOfPeaksFound(allowedPeakCount);
342
343 } else {
344 for (auto &goniometerMatrix : gonioVec) {
345 // Final transformation matrix (HKL to Q in lab frame)
346 DblMatrix orientedUB = goniometerMatrix * ub;
347
348 /* Because of the additional filtering step it's better to keep track of
349 * the allowed peaks with a counter. */
350 HKLFilterWavelength lambdaFilter(orientedUB, lambdaMin, lambdaMax);
351
352 size_t allowedPeakCount = 0;
353
354 bool useExtendedDetectorSpace = getProperty("PredictPeaksOutsideDetectors");
355 if (useExtendedDetectorSpace && !m_inst->getComponentByName("extended-detector-space")) {
356 g_log.warning() << "Attempting to find peaks outside of detectors but "
357 "no extended detector space has been defined\n";
358 }
359
360 for (auto &possibleHKL : possibleHKLs) {
361 if (lambdaFilter.isAllowed(possibleHKL)) {
362 calculateQAndAddToOutput(possibleHKL, orientedUB, goniometerMatrix);
363 ++allowedPeakCount;
364 }
365 prog.report();
366 }
367
368 logNumberOfPeaksFound(allowedPeakCount);
369 }
370 }
371
372 // Sort peaks by run number so that peaks with equal goniometer matrices are
373 // adjacent
374 std::vector<std::pair<std::string, bool>> criteria;
375 criteria.emplace_back("RunNumber", true);
377 criteria.emplace_back("BankName", true);
378 m_pw->sort(criteria);
379
380 for (int i = 0; i < static_cast<int>(m_pw->getNumberPeaks()); ++i) {
381 m_pw->getPeak(i).setPeakNumber(i);
382 }
383 setProperty<IPeaksWorkspace_sptr>("OutputWorkspace", m_pw);
384}
385
391void PredictPeaks::logNumberOfPeaksFound(size_t allowedPeakCount) const {
392 if (auto pw = std::dynamic_pointer_cast<PeaksWorkspace>(m_pw)) {
393 const bool usingExtendedDetectorSpace = getProperty("PredictPeaksOutsideDetectors");
394 const auto &peaks = pw->getPeaks();
395 size_t offDetectorPeakCount = 0;
396 size_t onDetectorPeakCount = 0;
397
398 for (const auto &peak : peaks) {
399 if (peak.getDetectorID() == -1) {
400 offDetectorPeakCount++;
401 } else {
402 onDetectorPeakCount++;
403 }
404 }
405
406 g_log.notice() << "Out of " << allowedPeakCount << " allowed peaks within parameters, " << onDetectorPeakCount
407 << " were found to hit a detector";
408 if (usingExtendedDetectorSpace) {
409 g_log.notice() << " and " << offDetectorPeakCount << " were found in "
410 << "extended detector space.";
411 }
412
413 g_log.notice() << "\n";
414 }
415}
416
419 // Check that there is an input workspace that has a sample.
420 if (!inWS || !inWS->getInstrument())
421 throw std::invalid_argument("Did not specify a valid InputWorkspace with a "
422 "full instrument.");
423
424 m_inst = inWS->getInstrument();
425}
426
429 if (!inWS) {
430 throw std::runtime_error("Failed to get run number");
431 }
432
433 m_runNumber = inWS->getRunNumber();
434}
435
438 const auto sample = m_inst->getSample();
439 if (!sample)
440 throw std::runtime_error("Instrument sample position has not been set");
441 const V3D samplePos = sample->getPos();
442
443 // L1 path and direction
444 V3D beamDir = m_inst->getSource()->getPos() - samplePos;
445
446 if ((fabs(beamDir.X()) > 1e-2) || (fabs(beamDir.Y()) > 1e-2)) // || (beamDir.Z() < 0))
447 throw std::invalid_argument("Instrument must have a beam direction that "
448 "is only in the +Z direction for this "
449 "algorithm to be valid..");
450}
451
455 std::vector<V3D> &possibleHKLs) const {
456 const double dMin = getProperty("MinDSpacing");
457 const double dMax = getProperty("MaxDSpacing");
458
459 // --- Reflection condition ----
460 // Use the primitive by default
461 ReflectionCondition_sptr refCond = std::make_shared<ReflectionConditionPrimitive>();
462 // Get it from the property
463 const std::string refCondName = getPropertyValue("ReflectionCondition");
464 const auto found = std::find_if(m_refConds.crbegin(), m_refConds.crend(), [&refCondName](const auto &condition) {
465 return condition->getName() == refCondName;
466 });
467 if (found != m_refConds.crend()) {
468 refCond = *found;
469 }
470
471 HKLGenerator gen(orientedLattice, dMin);
472 auto filter =
473 std::make_shared<HKLFilterCentering>(refCond) & std::make_shared<HKLFilterDRange>(orientedLattice, dMin, dMax);
474
475 V3D hklMin = *(gen.begin());
476
477 g_log.information() << "HKL range for d_min of " << dMin << " to d_max of " << dMax << " is from " << hklMin << " to "
478 << hklMin * -1.0 << ", a total of " << gen.size() << " possible HKL's\n";
479
480 if (gen.size() > 10000000000)
481 throw std::invalid_argument("More than 10 billion HKLs to search. Is "
482 "your d_min value too small?");
483
484 possibleHKLs.clear();
485 possibleHKLs.reserve(gen.size());
486 std::remove_copy_if(gen.begin(), gen.end(), std::back_inserter(possibleHKLs), (~filter)->fn());
487}
488
491 std::vector<V3D> &possibleHKLs) const {
492 possibleHKLs.clear();
493 possibleHKLs.reserve(peaksWorkspace->getNumberPeaks());
494
495 bool roundHKL = getProperty("RoundHKL");
496
497 /* Q is at the end multiplied with the factor determined in the
498 * constructor (-1 for crystallography, 1 otherwise). So to avoid
499 * "flippling HKLs" when it's not required, the HKLs of the input
500 * workspace are also multiplied by the factor that is appropriate
501 * for the convention stored in the workspace.
502 */
503 double peaks_q_convention_factor = qConventionFactor(peaksWorkspace->getConvention());
504
505 for (int i = 0; i < static_cast<int>(peaksWorkspace->getNumberPeaks()); ++i) {
506 IPeak &p = peaksWorkspace->getPeak(i);
507 // Get HKL from that peak
508 V3D hkl = p.getHKL() * peaks_q_convention_factor;
509
510 if (roundHKL)
511 hkl.round();
512
513 possibleHKLs.emplace_back(hkl);
514 } // for each hkl in the workspace
515}
516
531 bool calculateStructureFactors = getProperty("CalculateStructureFactors");
532
533 if (calculateStructureFactors && sample.hasCrystalStructure()) {
534 CrystalStructure crystalStructure = sample.getCrystalStructure();
535 crystalStructure.setCell(sample.getOrientedLattice());
536
537 m_sfCalculator = StructureFactorCalculatorFactory::create<StructureFactorCalculatorSummation>(crystalStructure);
538 }
539}
540
554void PredictPeaks::calculateQAndAddToOutput(const V3D &hkl, const DblMatrix &orientedUB,
555 const DblMatrix &goniometerMatrix) {
556 // The q-vector direction of the peak is = goniometer * ub * hkl_vector
557 // This is in inelastic convention: momentum transfer of the LATTICE!
558 // Also, q does have a 2pi factor = it is equal to 2pi/wavelength.
559 const auto q = orientedUB * hkl * (2.0 * M_PI * m_qConventionFactor);
560 const auto params = getPeakParametersFromQ(q);
561 const auto detectorDir = std::get<0>(params);
562 const auto wl = std::get<1>(params);
563
564 const bool useExtendedDetectorSpace = getProperty("PredictPeaksOutsideDetectors");
565 const auto result = m_detectorCacheSearch->findDetectorIndex(q);
566 const auto hitDetector = std::get<0>(result);
567 const auto index = std::get<1>(result);
568
569 if (!hitDetector && !useExtendedDetectorSpace) {
570 return;
571 }
572
573 const auto &detInfo = m_pw->detectorInfo();
574 const auto &det = detInfo.detector(index);
575 std::unique_ptr<Peak> peak;
576
577 if (hitDetector) {
578 // peak hit a detector to add it to the list
579 peak = std::make_unique<Peak>(m_inst, det.getID(), wl);
580 if (!peak->getDetector()) {
581 return;
582 }
583 } else if (useExtendedDetectorSpace) {
584 // use extended detector space to try and guess peak position
585 const auto returnedComponent = m_inst->getComponentByName("extended-detector-space");
586 // Check that the component is valid
587 const auto component = std::dynamic_pointer_cast<const ObjComponent>(returnedComponent);
588 if (!component)
589 throw std::runtime_error("PredictPeaks: user requested use of a extended "
590 "detector space to predict peaks but there is no"
591 "definition in the IDF");
592
593 // find where this Q vector should intersect with "extended" space
594 Geometry::Track track(detInfo.samplePosition(), detectorDir);
595 if (!component->interceptSurface(track))
596 return;
597
598 // The exit point is the vector to the place that we hit a detector
599 const auto magnitude = track.back().exitPoint.norm();
600 peak = std::make_unique<Peak>(m_inst, q, boost::optional<double>(magnitude));
601 }
602
603 if (m_edge > 0 && edgePixel(m_inst, peak->getBankName(), peak->getCol(), peak->getRow(), m_edge))
604 return;
605
606 // Only add peaks that hit the detector
607 peak->setGoniometerMatrix(goniometerMatrix);
608 // Save the run number found before.
609 peak->setRunNumber(m_runNumber);
610 peak->setHKL(hkl * m_qConventionFactor);
611 peak->setIntHKL(hkl * m_qConventionFactor);
612
613 if (m_sfCalculator) {
614 peak->setIntensity(m_sfCalculator->getFSquared(hkl));
615 }
616
617 // Add it to the workspace
618 m_pw->addPeak(*peak);
619}
620
632 // The q-vector direction of the peak is = goniometer * ub * hkl_vector
633 // This is in inelastic convention: momentum transfer of the LATTICE!
634 // Also, q does have a 2pi factor = it is equal to 2pi/wavelength.
635 const auto q = UB * hkl * (2.0 * M_PI * m_qConventionFactor);
636 auto peak = std::make_unique<LeanElasticPeak>(q);
637
638 // Save the run number found before.
639 peak->setRunNumber(m_runNumber);
640 peak->setHKL(hkl * m_qConventionFactor);
641 peak->setIntHKL(hkl * m_qConventionFactor);
642
643 if (m_sfCalculator) {
644 peak->setIntensity(m_sfCalculator->getFSquared(hkl));
645 }
646
647 // Add it to the workspace
648 m_pw->addPeak(*peak);
649}
650
656std::tuple<V3D, double> PredictPeaks::getPeakParametersFromQ(const V3D &q) const {
657 double norm_q = q.norm();
658 // Default for ki-kf has -q
659 const double qBeam = q.scalar_prod(m_refBeamDir) * m_qConventionFactor;
660 double one_over_wl = (norm_q * norm_q) / (2.0 * qBeam);
661 double wl = (2.0 * M_PI) / one_over_wl;
662 // Default for ki-kf has -q
663 V3D detectorDir = q * -m_qConventionFactor;
664 detectorDir[m_refFrame->pointingAlongBeam()] = one_over_wl - qBeam;
665 detectorDir.normalize();
666 return std::make_tuple(detectorDir, wl);
667}
668
672 m_refFrame = m_inst->getReferenceFrame();
673 m_refBeamDir = m_refFrame->vecPointingAlongBeam();
674}
675
676} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > * lambda
#define fabs(x)
Definition: Matrix.cpp:22
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
This class stores information about the sample used in particular run.
Definition: Sample.h:33
const Geometry::OrientedLattice & getOrientedLattice() const
Get a reference to the sample's OrientedLattice.
Definition: Sample.cpp:154
bool hasCrystalStructure() const
Returns true if the sample actually holds a CrystalStructure.
Definition: Sample.cpp:195
const Geometry::CrystalStructure & getCrystalStructure() const
Definition: Sample.cpp:181
A property class for workspaces.
void calculateQAndAddToOutputLeanElastic(const Kernel::V3D &hkl, const Kernel::DblMatrix &UB)
Calculates Q from HKL and adds a peak to the output workspace.
void setInstrumentFromInputWorkspace(const API::ExperimentInfo_sptr &inWS)
Tries to set the internally stored instrument from an ExperimentInfo-object.
int m_runNumber
Run number of input workspace.
Definition: PredictPeaks.h:90
void checkBeamDirection() const
Checks that the beam direction is +Z, throws exception otherwise.
int m_edge
Number of edge pixels with no peaks.
Definition: PredictPeaks.h:83
void setRunNumberFromInputWorkspace(const API::ExperimentInfo_sptr &inWS)
Sets the run number from the supplied ExperimentInfo or throws an exception.
void exec() override
Run the algorithm.
std::tuple< Kernel::V3D, double > getPeakParametersFromQ(const Kernel::V3D &q) const
Get the predicted detector direction from Q.
void logNumberOfPeaksFound(size_t allowedPeakCount) const
Log the number of peaks found to fall on and off detectors.
std::shared_ptr< const Geometry::ReferenceFrame > m_refFrame
Reference frame for the instrument.
Definition: PredictPeaks.h:94
void fillPossibleHKLsUsingPeaksWorkspace(const API::IPeaksWorkspace_sptr &peaksWorkspace, std::vector< Kernel::V3D > &possibleHKLs) const
Fills possibleHKLs with all HKLs from the supplied PeaksWorkspace.
Geometry::Instrument_const_sptr m_inst
Instrument reference.
Definition: PredictPeaks.h:92
void calculateQAndAddToOutput(const Kernel::V3D &hkl, const Kernel::DblMatrix &orientedUB, const Kernel::DblMatrix &goniometerMatrix)
Calculates Q from HKL and adds a peak to the output workspace.
std::unique_ptr< API::DetectorSearcher > m_detectorCacheSearch
Detector search cache for fast look-up of detectors.
Definition: PredictPeaks.h:88
std::vector< Mantid::Geometry::ReflectionCondition_sptr > m_refConds
Reflection conditions possible.
Definition: PredictPeaks.h:86
void setReferenceFrameAndBeamDirection()
Cache the reference frame and beam direction from the instrument.
Mantid::API::IPeaksWorkspace_sptr m_pw
Output peaks workspace.
Definition: PredictPeaks.h:98
Geometry::StructureFactorCalculator_sptr m_sfCalculator
Definition: PredictPeaks.h:99
void fillPossibleHKLsUsingGenerator(const Geometry::OrientedLattice &orientedLattice, std::vector< Kernel::V3D > &possibleHKLs) const
Fills possibleHKLs with all HKLs that are allowed within d- and lambda-limits.
void init() override
Initialise the properties.
void setStructureFactorCalculatorFromSample(const API::Sample &sample)
Assigns a StructureFactorCalculator if available in sample.
Kernel::V3D m_refBeamDir
Direction of the beam for this instrument.
Definition: PredictPeaks.h:96
Three components are required to describe a crystal structure:
void setCell(const UnitCell &cell)
Assigns a new unit cell.
Class to represent a particular goniometer setting, which is described by the rotation matrix.
Definition: Goniometer.h:55
std::vector< double > getEulerAngles(const std::string &convention="YZX")
Return Euler angles acording to a convention.
Definition: Goniometer.cpp:282
const Kernel::DblMatrix & getR() const
Return global rotation matrix.
Definition: Goniometer.cpp:80
void calcFromQSampleAndWavelength(const Mantid::Kernel::V3D &Q, double wavelength, bool flip_x=false, bool inner=false)
Calculate goniometer for rotation around y-asix for constant wavelength from Q Sample.
Definition: Goniometer.cpp:197
bool isAllowed(const Kernel::V3D &hkl) const noexcept override
Returns true if lambda of the reflection is within the limits.
const const_iterator & end() const
Returns an iterator which "points at" one element past the end.
Definition: HKLGenerator.h:146
const const_iterator & begin() const
Returns an iterator to the beginning of the sequence.
Definition: HKLGenerator.h:143
size_t size() const
Returns the number of HKLs to be generated.
Definition: HKLGenerator.h:140
Structure describing a single-crystal peak.
Definition: IPeak.h:26
virtual Mantid::Kernel::Matrix< double > getGoniometerMatrix() const =0
virtual Mantid::Kernel::V3D getHKL() const =0
Class to implement UB matrix.
const Kernel::DblMatrix & getUB() const
Get the UB matrix.
Defines a track as a start point and a direction.
Definition: Track.h:165
LType::reference back()
Returns a reference to the last link.
Definition: Track.h:213
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
void setNotifyStep(double notifyStepPct)
Override the frequency at which notifications are sent out.
The concrete, templated class for properties.
Class for 3D vectors.
Definition: V3D.h:34
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
Definition: V3D.h:274
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
double norm() const noexcept
Definition: V3D.h:263
constexpr double norm2() const noexcept
Vector length squared.
Definition: V3D.h:265
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
void round() noexcept
Round each component to the nearest integer.
Definition: V3D.cpp:140
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< ExperimentInfo > ExperimentInfo_sptr
Shared pointer to ExperimentInfo.
std::shared_ptr< MultipleExperimentInfos > MultipleExperimentInfos_sptr
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
double qConventionFactor()
convenience overload to pull the convention from the config service
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< ReflectionCondition > ReflectionCondition_sptr
Shared pointer to a ReflectionCondition.
MANTID_GEOMETRY_DLL bool edgePixel(const Geometry::Instrument_const_sptr &inst, const std::string &bankName, int col, int row, int Edge)
Function to find peaks near detector edge.
Definition: EdgePixel.cpp:22
MANTID_GEOMETRY_DLL const ReflectionConditions & getAllReflectionConditions()
Mantid::Kernel::Matrix< double > DblMatrix
Definition: Matrix.h:206
Peak indexing algorithm, which works by assigning multiple possible HKL values to each peak and then ...
Definition: IndexSXPeaks.h:29
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54