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.");
60 "Minimum wavelength limit at which to start looking for single-crystal "
63 "Maximum wavelength limit at which to stop looking for single-crystal "
67 "Minimum d-spacing of peaks to consider. Default = 1.0");
69 "Maximum d-spacing of peaks to consider.");
72 "This will calculate the goniometer rotation (around y-axis "
73 "only) for a constant wavelength.");
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");
81 "Whether the goniometer to be calculated is the most inner "
82 "(phi) or most outer (omega)");
84 std::make_unique<EnabledWhenProperty>(
"CalculateGoniometerForCW",
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");
95 "Minimum goniometer rotation angle");
99 "Maximum goniometer rotation angle");
103 std::vector<std::string> 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?");
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.");
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.");
124 "When using HKLPeaksWorkspace, this will round the HKL "
125 "values in the HKLPeaksWorkspace to the nearest integers if "
127 "Keep unchecked to use the original values");
133 std::unique_ptr<IPropertySettings> set = std::make_unique<EnabledWhenProperty>(
"HKLPeaksWorkspace",
IS_DEFAULT);
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).");
153 "An output PeaksWorkspace.");
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.");
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. ");
174 ExperimentInfo_sptr inputExperimentInfo = std::dynamic_pointer_cast<ExperimentInfo>(rawInputWorkspace);
177 PeaksWorkspace_sptr peaksWS = std::dynamic_pointer_cast<PeaksWorkspace>(rawInputWorkspace);
180 std::vector<DblMatrix> gonioVec;
184 auto goniometerMatrix = matrixWS->run().getGoniometerMatrices();
185 std::move(goniometerMatrix.begin(), goniometerMatrix.end(), back_inserter(gonioVec));
186 }
catch (std::runtime_error &e) {
188 g_log.
error() <<
"Error getting the goniometer rotation matrix from the "
191 g_log.
warning() <<
"Using identity goniometer rotation matrix instead.\n";
193 }
else if (peaksWS) {
196 std::vector<std::pair<std::string, bool>> criteria;
197 criteria.emplace_back(
"RunNumber",
true);
199 peaksWS->sort(criteria);
203 for (
int i = 0; i < static_cast<int>(peaksWS->getNumberPeaks()); ++i) {
204 IPeak &p = peaksWS->getPeak(i);
206 if (!(currentGoniometerMatrix == lastGoniometerMatrix)) {
207 gonioVec.emplace_back(currentGoniometerMatrix);
208 lastGoniometerMatrix = currentGoniometerMatrix;
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.");
218 inputExperimentInfo = mdWS->getExperimentInfo(0);
221 for (uint16_t i = 0; i < mdWS->getNumExperimentInfo(); ++i) {
223 auto goniometerMatrix = mdWS->getExperimentInfo(i)->mutableRun().getGoniometerMatrices();
224 std::move(goniometerMatrix.begin(), goniometerMatrix.end(), back_inserter(gonioVec));
225 }
catch (std::runtime_error &e) {
227 gonioVec.emplace_back(
DblMatrix(3, 3,
true));
229 g_log.
error() <<
"Error getting the goniometer rotation matrix from the "
232 g_log.
warning() <<
"Using identity goniometer rotation matrix instead.\n";
239 if (gonioVec.empty()) {
240 gonioVec.emplace_back(
DblMatrix(3, 3,
true));
243 if (usingInstrument) {
252 m_pw = std::make_shared<LeanElasticPeaksWorkspace>();
254 m_pw = std::make_shared<PeaksWorkspace>();
257 m_pw->copyExperimentInfoFrom(inputExperimentInfo.get());
259 const Sample &sample = inputExperimentInfo->sample();
266 ub = orientedLattice.
getUB();
268 std::vector<V3D> possibleHKLs;
271 if (!possibleHKLWorkspace) {
285 Progress prog(
this, 0.0, 1.0, possibleHKLs.size() * gonioVec.size());
291 if (!usingInstrument) {
292 for (
auto &possibleHKL : possibleHKLs) {
295 }
else if (
getProperty(
"CalculateGoniometerForCW")) {
296 size_t allowedPeakCount = 0;
299 if (wavelength == DBL_MAX) {
300 if (
m_inst->hasParameter(
"wavelength")) {
301 wavelength =
m_inst->getNumberParameter(
"wavelength").at(0);
303 throw std::runtime_error(
"Could not get wavelength, neither Wavelength algorithm property "
304 "set nor instrument wavelength parameter");
309 bool innerGoniometer =
getProperty(
"InnerGoniometer");
311 for (
auto &possibleHKL : possibleHKLs) {
316 double angle = innerGoniometer ? angles[2] : angles[0];
322 V3D q_lab = goniometer.
getR() * q_sample;
327 double lambda = (4.0 * M_PI * std::abs(q_lab.
Z())) / q_lab.
norm2();
329 if (!std::isfinite(angle) || angle < angleMin || angle > angleMax)
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";
344 for (
auto &goniometerMatrix : gonioVec) {
346 DblMatrix orientedUB = goniometerMatrix * ub;
352 size_t allowedPeakCount = 0;
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";
360 for (
auto &possibleHKL : possibleHKLs) {
361 if (lambdaFilter.
isAllowed(possibleHKL)) {
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);
380 for (
int i = 0; i < static_cast<int>(
m_pw->getNumberPeaks()); ++i) {
381 m_pw->getPeak(i).setPeakNumber(i);
383 setProperty<IPeaksWorkspace_sptr>(
"OutputWorkspace",
m_pw);
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;
398 for (
const auto &peak : peaks) {
399 if (peak.getDetectorID() == -1) {
400 offDetectorPeakCount++;
402 onDetectorPeakCount++;
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.";
420 if (!inWS || !inWS->getInstrument())
421 throw std::invalid_argument(
"Did not specify a valid InputWorkspace with a "
424 m_inst = inWS->getInstrument();
430 throw std::runtime_error(
"Failed to get run number");
438 const auto sample =
m_inst->getSample();
440 throw std::runtime_error(
"Instrument sample position has not been set");
441 const V3D samplePos = sample->getPos();
444 V3D beamDir =
m_inst->getSource()->getPos() - samplePos;
446 if ((
fabs(beamDir.
X()) > 1e-2) || (
fabs(beamDir.
Y()) > 1e-2))
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..");
455 std::vector<V3D> &possibleHKLs)
const {
464 const auto found = std::find_if(
m_refConds.crbegin(),
m_refConds.crend(), [&refCondName](
const auto &condition) {
465 return condition->getName() == refCondName;
473 std::make_shared<HKLFilterCentering>(refCond) & std::make_shared<HKLFilterDRange>(orientedLattice, dMin, dMax);
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";
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?");
484 possibleHKLs.clear();
485 possibleHKLs.reserve(gen.
size());
486 std::remove_copy_if(gen.
begin(), gen.
end(), std::back_inserter(possibleHKLs), (~filter)->fn());
491 std::vector<V3D> &possibleHKLs)
const {
492 possibleHKLs.clear();
493 possibleHKLs.reserve(peaksWorkspace->getNumberPeaks());
503 double peaks_q_convention_factor =
qConventionFactor(peaksWorkspace->getConvention());
505 for (
int i = 0; i < static_cast<int>(peaksWorkspace->getNumberPeaks()); ++i) {
506 IPeak &p = peaksWorkspace->getPeak(i);
508 V3D hkl = p.
getHKL() * peaks_q_convention_factor;
513 possibleHKLs.emplace_back(hkl);
531 bool calculateStructureFactors =
getProperty(
"CalculateStructureFactors");
537 m_sfCalculator = StructureFactorCalculatorFactory::create<StructureFactorCalculatorSummation>(crystalStructure);
561 const auto detectorDir = std::get<0>(params);
562 const auto wl = std::get<1>(params);
564 const bool useExtendedDetectorSpace =
getProperty(
"PredictPeaksOutsideDetectors");
566 const auto hitDetector = std::get<0>(result);
567 const auto index = std::get<1>(result);
569 if (!hitDetector && !useExtendedDetectorSpace) {
573 const auto &detInfo =
m_pw->detectorInfo();
574 const auto &det = detInfo.detector(
index);
575 std::unique_ptr<Peak> peak;
579 peak = std::make_unique<Peak>(
m_inst, det.getID(), wl);
580 if (!peak->getDetector()) {
583 }
else if (useExtendedDetectorSpace) {
585 const auto returnedComponent =
m_inst->getComponentByName(
"extended-detector-space");
587 const auto component = std::dynamic_pointer_cast<const ObjComponent>(returnedComponent);
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");
595 if (!component->interceptSurface(track))
599 const auto magnitude = track.
back().exitPoint.norm();
600 peak = std::make_unique<Peak>(
m_inst, q, boost::optional<double>(magnitude));
607 peak->setGoniometerMatrix(goniometerMatrix);
618 m_pw->addPeak(*peak);
636 auto peak = std::make_unique<LeanElasticPeak>(q);
648 m_pw->addPeak(*peak);
657 double norm_q = q.
norm();
660 double one_over_wl = (norm_q * norm_q) / (2.0 * qBeam);
661 double wl = (2.0 * M_PI) / one_over_wl;
664 detectorDir[
m_refFrame->pointingAlongBeam()] = one_over_wl - qBeam;
666 return std::make_tuple(detectorDir, wl);
#define DECLARE_ALGORITHM(classname)
const std::vector< double > * lambda
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Helper class for reporting progress from algorithms.
This class stores information about the sample used in particular run.
const Geometry::OrientedLattice & getOrientedLattice() const
Get a reference to the sample's OrientedLattice.
bool hasCrystalStructure() const
Returns true if the sample actually holds a CrystalStructure.
const Geometry::CrystalStructure & getCrystalStructure() const
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.
double m_qConventionFactor
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.
void checkBeamDirection() const
Checks that the beam direction is +Z, throws exception otherwise.
int m_edge
Number of edge pixels with no peaks.
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.
PredictPeaks()
Constructor.
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.
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.
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.
std::vector< Mantid::Geometry::ReflectionCondition_sptr > m_refConds
Reflection conditions possible.
void setReferenceFrameAndBeamDirection()
Cache the reference frame and beam direction from the instrument.
Mantid::API::IPeaksWorkspace_sptr m_pw
Output peaks workspace.
Geometry::StructureFactorCalculator_sptr m_sfCalculator
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.
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.
std::vector< double > getEulerAngles(const std::string &convention="YZX")
Return Euler angles acording to a convention.
const Kernel::DblMatrix & getR() const
Return global rotation matrix.
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.
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.
const const_iterator & begin() const
Returns an iterator to the beginning of the sequence.
size_t size() const
Returns the number of HKLs to be generated.
Structure describing a single-crystal peak.
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.
LType::reference back()
Returns a reference to the last link.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void notice(const std::string &msg)
Logs at notice level.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
void setNotifyStep(double notifyStepPct)
Override the frequency at which notifications are sent out.
The concrete, templated class for properties.
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
constexpr double X() const noexcept
Get x.
double normalize()
Make a normalized vector (return norm value)
constexpr double Y() const noexcept
Get y.
double norm() const noexcept
constexpr double norm2() const noexcept
Vector length squared.
constexpr double Z() const noexcept
Get z.
void round() noexcept
Round each component to the nearest integer.
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
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.
MANTID_GEOMETRY_DLL const ReflectionConditions & getAllReflectionConditions()
Mantid::Kernel::Matrix< double > DblMatrix
Peak indexing algorithm, which works by assigning multiple possible HKL values to each peak and then ...
@ Input
An input workspace.
@ Output
An output workspace.