Mantid
Loading...
Searching...
No Matches
SortHKL.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/Sample.h"
11#include "MantidAPI/TableRow.h"
12#include "MantidAPI/TextAxis.h"
23#include "MantidKernel/Utils.h"
24
25#include <cmath>
26#include <fstream>
27#include <numeric>
28
29using namespace Mantid::Geometry;
30using namespace Mantid::DataObjects;
31using namespace Mantid::Kernel;
32using namespace Mantid::API;
34
35namespace Mantid::Crystal {
36
37// Register the algorithm into the AlgorithmFactory
39
40SortHKL::SortHKL() : m_pointGroups(getAllPointGroups()), m_refConds(getAllReflectionConditions()) {}
41
42SortHKL::~SortHKL() = default;
43
45 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("InputWorkspace", "", Direction::Input),
46 "An input PeaksWorkspace with an instrument.");
47
48 /* TODO: These two properties with string lists keep appearing -
49 * Probably there should be a dedicated Property type or validator. */
50 std::vector<std::string> pgOptions;
51 pgOptions.reserve(2 * m_pointGroups.size() + 5);
52 std::transform(m_pointGroups.cbegin(), m_pointGroups.cend(), std::back_inserter(pgOptions),
53 [](const auto &group) { return group->getSymbol(); });
54 std::transform(m_pointGroups.cbegin(), m_pointGroups.cend(), std::back_inserter(pgOptions),
55 [](const auto &group) { return group->getName(); });
56 // Scripts may have Orthorhombic misspelled from past bug in PointGroupFactory
57 pgOptions.emplace_back("222 (Orthorombic)");
58 pgOptions.emplace_back("mm2 (Orthorombic)");
59 pgOptions.emplace_back("2mm (Orthorombic)");
60 pgOptions.emplace_back("m2m (Orthorombic)");
61 pgOptions.emplace_back("mmm (Orthorombic)");
62 declareProperty("PointGroup", pgOptions[0], std::make_shared<StringListValidator>(pgOptions),
63 "Which point group applies to this crystal?");
64
65 std::vector<std::string> centeringOptions;
66 centeringOptions.reserve(2 * m_refConds.size());
67 std::transform(m_refConds.cbegin(), m_refConds.cend(), std::back_inserter(centeringOptions),
68 [](const auto &condition) { return condition->getSymbol(); });
69 std::transform(m_refConds.cbegin(), m_refConds.cend(), std::back_inserter(centeringOptions),
70 [](const auto &condition) { return condition->getName(); });
71 declareProperty("LatticeCentering", centeringOptions[0], std::make_shared<StringListValidator>(centeringOptions),
72 "Appropriate lattice centering for the peaks.");
73
74 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
75 "Output PeaksWorkspace");
76 declareProperty("OutputChi2", 0.0, "Chi-square is available as output", Direction::Output);
78 std::make_unique<WorkspaceProperty<ITableWorkspace>>("StatisticsTable", "StatisticsTable", Direction::Output),
79 "An output table workspace for the statistics of the peaks.");
80 declareProperty(std::make_unique<PropertyWithValue<std::string>>("RowName", "Overall", Direction::Input),
81 "name of row");
82 declareProperty("Append", false,
83 "Append to output table workspace if true.\n"
84 "If false, new output table workspace (default).");
85 const std::vector<std::string> equivTypes{"Mean", "Median"};
86 declareProperty("EquivalentIntensities", equivTypes.front(), std::make_shared<StringListValidator>(equivTypes),
87 "Replace intensities by mean(default), "
88 "or median.");
89 declareProperty(std::make_unique<PropertyWithValue<double>>("SigmaCritical", 3.0, Direction::Input),
90 "Removes peaks whose intensity deviates more than "
91 "SigmaCritical from the mean (or median).");
92 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("EquivalentsWorkspace", "EquivalentIntensities",
94 "Output Equivalent Intensities");
95 declareProperty("WeightedZScore", false,
96 "Use weighted ZScore if true.\n"
97 "If false, standard ZScore (default).");
98}
99
101 PeaksWorkspace_sptr inputPeaksWorkspace = getProperty("InputWorkspace");
102
103 const std::vector<Peak> &inputPeaks = inputPeaksWorkspace->getPeaks();
104 std::vector<Peak> peaks = getNonZeroPeaks(inputPeaks);
105
106 if (peaks.empty()) {
107 g_log.warning() << "Number of peaks should not be 0 for SortHKL.\n";
108 return;
109 }
110
111 UnitCell cell = inputPeaksWorkspace->sample().getOrientedLattice();
112
113 UniqueReflectionCollection uniqueReflections = getUniqueReflections(peaks, cell);
114 std::string equivalentIntensities = getPropertyValue("EquivalentIntensities");
115 double sigmaCritical = getProperty("SigmaCritical");
116 bool weightedZ = getProperty("WeightedZScore");
117
119 "Workspace2D", uniqueReflections.getReflections().size(), 20, 20);
120 int counter = 0;
121 size_t maxPeaks = 0;
122 auto tAxis = std::make_unique<TextAxis>(uniqueReflections.getReflections().size());
123 UniqWksp->getAxis(0)->unit() = UnitFactory::Instance().create("Wavelength");
124 for (const auto &unique : uniqueReflections.getReflections()) {
125 /* Since all possible unique reflections are explored
126 * there may be 0 observations for some of them.
127 * In that case, nothing can be done.*/
128
129 if (unique.second.count() > 2) {
130 tAxis->setLabel(counter, " " + unique.second.getHKL().toString());
131 auto &UniqX = UniqWksp->mutableX(counter);
132 auto &UniqY = UniqWksp->mutableY(counter);
133 auto &UniqE = UniqWksp->mutableE(counter);
134 counter++;
135 auto wavelengths = unique.second.getWavelengths();
136 auto intensities = unique.second.getIntensities();
137 g_log.debug() << "HKL " << unique.second.getHKL() << "\n";
138 g_log.debug() << "Intensities ";
139 for (const auto &e : intensities)
140 g_log.debug() << e << " ";
141 g_log.debug() << "\n";
142 std::vector<double> zScores;
143 if (!weightedZ) {
144 zScores = Kernel::getZscore(intensities);
145 } else {
146 auto sigmas = unique.second.getSigmas();
147 zScores = Kernel::getWeightedZscore(intensities, sigmas);
148 }
149
150 if (zScores.size() > maxPeaks)
151 maxPeaks = zScores.size();
152 // Possibly remove outliers.
153 auto outliersRemoved = unique.second.removeOutliers(sigmaCritical, weightedZ);
154
155 auto intensityStatistics =
156 Kernel::getStatistics(outliersRemoved.getIntensities(), StatOptions::Mean | StatOptions::Median);
157
158 g_log.debug() << "Mean = " << intensityStatistics.mean << " Median = " << intensityStatistics.median << "\n";
159 // sort wavelengths & intensities
160 for (size_t i = 0; i < wavelengths.size(); i++) {
161 size_t i0 = i;
162 for (size_t j = i + 1; j < wavelengths.size(); j++) {
163 if (wavelengths[j] < wavelengths[i0]) // Change was here!
164 {
165 i0 = j;
166 }
167 }
168 double temp = wavelengths[i0];
169 wavelengths[i0] = wavelengths[i];
170 wavelengths[i] = temp;
171 temp = intensities[i0];
172 intensities[i0] = intensities[i];
173 intensities[i] = temp;
174 }
175 g_log.debug() << "Zscores ";
176 for (size_t i = 0; i < std::min(zScores.size(), static_cast<size_t>(20)); ++i) {
177 UniqX[i] = wavelengths[i];
178 UniqY[i] = intensities[i];
179 if (zScores[i] > sigmaCritical)
180 UniqE[i] = intensities[i];
181 else if (equivalentIntensities == "Mean")
182 UniqE[i] = intensityStatistics.mean - intensities[i];
183 else
184 UniqE[i] = intensityStatistics.median - intensities[i];
185 g_log.debug() << zScores[i] << " ";
186 }
187 g_log.debug() << "\n";
188 }
189 }
190
191 if (counter > 0) {
192 MatrixWorkspace_sptr UniqWksp2 =
193 Mantid::API::WorkspaceFactory::Instance().create("Workspace2D", counter, maxPeaks, maxPeaks);
194 for (int64_t i = 0; i < counter; ++i) {
195 auto &outSpec = UniqWksp2->getSpectrum(i);
196 const auto &inSpec = UniqWksp->getSpectrum(i);
197 outSpec.setHistogram(inSpec.histogram());
198 // Copy the spectrum number/detector IDs
199 outSpec.copyInfoFrom(inSpec);
200 }
201 UniqWksp2->replaceAxis(1, std::move(tAxis));
202 setProperty("EquivalentsWorkspace", UniqWksp2);
203 } else {
204 setProperty("EquivalentsWorkspace", UniqWksp);
205 }
206
207 PeaksStatistics peaksStatistics(uniqueReflections, equivalentIntensities, sigmaCritical, weightedZ);
208
209 // Store the statistics for output.
210 const std::string tableName = getProperty("StatisticsTable");
211 ITableWorkspace_sptr statisticsTable = getStatisticsTable(tableName);
212 insertStatisticsIntoTable(statisticsTable, peaksStatistics);
213
214 // Store all peaks that were used to calculate statistics.
215
216 PeaksWorkspace_sptr outputPeaksWorkspace = getOutputPeaksWorkspace(inputPeaksWorkspace);
217
218 std::vector<Peak> &originalOutputPeaks = outputPeaksWorkspace->getPeaks();
219 originalOutputPeaks.swap(peaksStatistics.m_peaks);
220
221 sortOutputPeaksByHKL(outputPeaksWorkspace);
222
223 setProperty("OutputWorkspace", outputPeaksWorkspace);
224 setProperty("OutputChi2", peaksStatistics.m_chiSquared);
225 setProperty("StatisticsTable", statisticsTable);
226 AnalysisDataService::Instance().addOrReplace(tableName, statisticsTable);
227}
228
231std::vector<Peak> SortHKL::getNonZeroPeaks(const std::vector<Peak> &inputPeaks) const {
232 std::vector<Peak> peaks;
233 peaks.reserve(inputPeaks.size());
234
235 std::remove_copy_if(inputPeaks.begin(), inputPeaks.end(), std::back_inserter(peaks), [](const Peak &peak) {
236 return peak.getIntensity() <= 0.0 || peak.getSigmaIntensity() <= 0.0 || peak.getIntMNP() != V3D(0, 0, 0) ||
237 peak.getHKL() == V3D(0, 0, 0);
238 });
239
240 return peaks;
241}
242
254UniqueReflectionCollection SortHKL::getUniqueReflections(const std::vector<Peak> &peaks, const UnitCell &cell) const {
256 PointGroup_sptr pointGroup = getPointgroup();
257
258 std::pair<double, double> dLimits = getDLimits(peaks, cell);
259
260 UniqueReflectionCollection reflections(cell, dLimits, pointGroup, centering);
261 reflections.addObservations(peaks);
262
263 return reflections;
264}
265
268 ReflectionCondition_sptr centering = std::make_shared<ReflectionConditionPrimitive>();
269
270 const std::string refCondName = getPropertyValue("LatticeCentering");
271 const auto found = std::find_if(m_refConds.crbegin(), m_refConds.crend(),
272 [refCondName](const auto &condition) { return condition->getName() == refCondName; });
273 if (found != m_refConds.crend()) {
274 centering = *found;
275 }
276
277 return centering;
278}
279
283 PointGroup_sptr pointGroup = PointGroupFactory::Instance().createPointGroup("-1");
284
285 std::string pointGroupName = getPropertyValue("PointGroup");
286 size_t pos = pointGroupName.find("Orthorombic");
287 if (pos != std::string::npos) {
288 g_log.warning() << "Orthorhomic is misspelled in your script.\n";
289 pointGroupName.replace(pos, 11, "Orthorhombic");
290 g_log.warning() << "Please correct to " << pointGroupName << ".\n";
291 }
292 const auto found = std::find_if(m_pointGroups.crbegin(), m_pointGroups.crend(),
293 [&pointGroupName](const auto &group) { return group->getName() == pointGroupName; });
294 if (found != m_pointGroups.crend()) {
295 pointGroup = *found;
296 }
297
298 return pointGroup;
299}
300
304std::pair<double, double> SortHKL::getDLimits(const std::vector<Peak> &peaks, const UnitCell &cell) const {
305 auto dLimitIterators = std::minmax_element(peaks.begin(), peaks.end(), [cell](const Peak &lhs, const Peak &rhs) {
306 return cell.d(lhs.getHKL()) < cell.d(rhs.getHKL());
307 });
308
309 return std::make_pair(cell.d((*dLimitIterators.first).getHKL()), cell.d((*dLimitIterators.second).getHKL()));
310}
311
314ITableWorkspace_sptr SortHKL::getStatisticsTable(const std::string &name) const {
315 TableWorkspace_sptr tablews;
316
317 // Init or append to a table workspace
318 bool append = getProperty("Append");
319 if (append && AnalysisDataService::Instance().doesExist(name)) {
320 tablews = AnalysisDataService::Instance().retrieveWS<TableWorkspace>(name);
321 } else {
322 tablews = std::make_shared<TableWorkspace>();
323 tablews->addColumn("str", "Resolution Shell");
324 tablews->addColumn("int", "No. of Unique Reflections");
325 tablews->addColumn("double", "Resolution Min");
326 tablews->addColumn("double", "Resolution Max");
327 tablews->addColumn("double", "Multiplicity");
328 tablews->addColumn("double", "Mean ((I)/sd(I))");
329 tablews->addColumn("double", "Rmerge");
330 tablews->addColumn("double", "Rpim");
331 tablews->addColumn("double", "Data Completeness");
332 }
333
334 return tablews;
335}
336
340 if (!table) {
341 throw std::runtime_error("Can't store statistics into Null-table.");
342 }
343
344 std::string name = getProperty("RowName");
345 double completeness = 0.0;
346 if (name.substr(0, 4) != "bank") {
347 completeness = static_cast<double>(statistics.m_completeness);
348 }
349
350 // append to the table workspace
351 API::TableRow newrow = table->appendRow();
352
353 newrow << name << statistics.m_uniqueReflections << statistics.m_dspacingMin << statistics.m_dspacingMax
354 << statistics.m_redundancy << statistics.m_meanIOverSigma << 100.0 * statistics.m_rMerge
355 << 100.0 * statistics.m_rPim << 100.0 * completeness;
356}
357
360 PeaksWorkspace_sptr outputPeaksWorkspace = getProperty("OutputWorkspace");
361 if (outputPeaksWorkspace != inputPeaksWorkspace) {
362 outputPeaksWorkspace = inputPeaksWorkspace->clone();
363 }
364
365 return outputPeaksWorkspace;
366}
367
369void SortHKL::sortOutputPeaksByHKL(const IPeaksWorkspace_sptr &outputPeaksWorkspace) {
370 // Sort by HKL
371 std::vector<std::pair<std::string, bool>> criteria{{"H", true}, {"K", true}, {"L", true}};
372 outputPeaksWorkspace->sort(criteria);
373}
374
375} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > & rhs
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
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
The PeaksStatistics class is a small helper class that is used in SortHKL.
double m_meanIOverSigma
Average signal to noise ratio in the reflections.
double m_redundancy
Average number of observations for a unique reflection.
double m_dspacingMin
Lower d-spacing limit in the data set, sometimes referred to as upper resolution limit.
double m_rPim
Precision indicating R-factor (R_{p.i.m}).
double m_completeness
Fraction of observed unique reflections in the resolution range defined by d_min and d_max.
double m_dspacingMax
Upper d-spacing limit in the data set.
double m_rMerge
Merging R-factor, R_merge, sometimes also called R_sym.
This class computes all possible unique reflections within the specified d-limits,...
const std::map< Kernel::V3D, UniqueReflection > & getReflections() const
Returns the internally stored reflection map.
void addObservations(const std::vector< DataObjects::Peak > &peaks)
Assigns the supplied peaks to the proper UniqueReflection.
Save a PeaksWorkspace to a Gsas-style ASCII .hkl file.
Definition: SortHKL.h:32
void sortOutputPeaksByHKL(const API::IPeaksWorkspace_sptr &outputPeaksWorkspace)
Sorts the peaks in the workspace by H, K and L.
Definition: SortHKL.cpp:369
std::vector< Mantid::Geometry::ReflectionCondition_sptr > m_refConds
Reflection conditions.
Definition: SortHKL.h:79
API::ITableWorkspace_sptr getStatisticsTable(const std::string &name) const
Create a TableWorkspace for the statistics with appropriate columns or get one from the ADS.
Definition: SortHKL.cpp:314
DataObjects::PeaksWorkspace_sptr getOutputPeaksWorkspace(const DataObjects::PeaksWorkspace_sptr &inputPeaksWorkspace) const
Returns a PeaksWorkspace which is either the input workspace or a clone.
Definition: SortHKL.cpp:359
std::pair< double, double > getDLimits(const std::vector< DataObjects::Peak > &peaks, const Geometry::UnitCell &cell) const
Returns the lowest and highest d-Value in the list.
Definition: SortHKL.cpp:304
void insertStatisticsIntoTable(const API::ITableWorkspace_sptr &table, const PeakStatisticsTools::PeaksStatistics &statistics) const
Inserts statistics the supplied PeaksStatistics-objects into the supplied TableWorkspace.
Definition: SortHKL.cpp:339
Geometry::ReflectionCondition_sptr getCentering() const
Returns the centering extracted from the user-supplied property.
Definition: SortHKL.cpp:267
Geometry::PointGroup_sptr getPointgroup() const
Returns the PointGroup-object constructed from the property supplied by the user.
Definition: SortHKL.cpp:282
std::vector< DataObjects::Peak > getNonZeroPeaks(const std::vector< DataObjects::Peak > &inputPeaks) const
Returns a vector which contains only peaks with I > 0, sigma > 0 and valid HKL.
Definition: SortHKL.cpp:231
const std::string name() const override
Algorithm's name for identification.
Definition: SortHKL.h:38
void init() override
Virtual method - must be overridden by concrete algorithm.
Definition: SortHKL.cpp:44
PeakStatisticsTools::UniqueReflectionCollection getUniqueReflections(const std::vector< DataObjects::Peak > &peaks, const Geometry::UnitCell &cell) const
SortHKL::getUniqueReflections.
Definition: SortHKL.cpp:254
std::vector< Mantid::Geometry::PointGroup_sptr > m_pointGroups
Point Groups possible.
Definition: SortHKL.h:76
void exec() override
Virtual method - must be overridden by concrete algorithm.
Definition: SortHKL.cpp:100
Structure describing a single-crystal peak.
Definition: Peak.h:34
TableWorkspace is an implementation of Workspace in which the data are organised in columns of same s...
Class to implement unit cell of crystals.
Definition: UnitCell.h:44
double d(double h, double k, double l) const
Return d-spacing ( ) for a given h,k,l coordinate.
Definition: UnitCell.cpp:700
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
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.
std::shared_ptr< PointGroup > PointGroup_sptr
Shared pointer to a PointGroup.
Definition: PointGroup.h:67
MANTID_GEOMETRY_DLL const ReflectionConditions & getAllReflectionConditions()
MANTID_GEOMETRY_DLL std::vector< PointGroup_sptr > getAllPointGroups()
Definition: PointGroup.cpp:195
Statistics getStatistics(const std::vector< TYPE > &data, const unsigned int flags=StatOptions::AllStats)
Return a statistics object for the given data set.
Definition: Statistics.cpp:167
std::vector< double > getZscore(const std::vector< TYPE > &data)
Return the Z score values for a dataset.
Definition: Statistics.cpp:81
std::vector< double > getWeightedZscore(const std::vector< TYPE > &data, const std::vector< TYPE > &weights)
There are enough special cases in determining the Z score where it useful to put it in a single funct...
Definition: Statistics.cpp:102
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54