Mantid
Loading...
Searching...
No Matches
ConvertCWSDMDtoHKL.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
12#include "MantidAPI/Run.h"
13#include "MantidAPI/Sample.h"
15
22
30
31#include <fstream>
32
33namespace Mantid::MDAlgorithms {
34
35using namespace Mantid::API;
36using namespace Mantid::Kernel;
37using namespace Mantid::MDAlgorithms;
38using namespace Mantid::DataObjects;
39using namespace Mantid::Geometry;
40
42
43
45void ConvertCWSDMDtoHKL::init() {
46 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
47 "Name of the input MDEventWorkspace that stores detectors "
48 "counts from a constant-wave powder diffraction experiment.");
49
50 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace", "", Direction::Input,
52 "Input Peaks Workspace");
53
54 declareProperty(std::make_unique<ArrayProperty<double>>("UBMatrix"),
55 "A comma seperated list of doubles for UB matrix from (0,0), (0,1)"
56 "... (2,1),(2,2)");
57
58 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("OutputWorkspace", "", Direction::Output),
59 "Name of the output MDEventWorkspace in HKL-space.");
60
61 declareProperty(std::make_unique<FileProperty>("QSampleFileName", "", API::FileProperty::OptionalSave),
62 "Name of file for sample sample.");
63
64 declareProperty(std::make_unique<FileProperty>("HKLFileName", "", API::FileProperty::OptionalSave),
65 "Name of file for HKL.");
66}
67
71 // Get inputs
72 IMDEventWorkspace_sptr inputWS = getProperty("InputWorkspace");
73 if (inputWS->getSpecialCoordinateSystem() != Mantid::Kernel::QSample) {
74 std::stringstream errmsg;
75 errmsg << "Input MDEventWorkspace's coordinate system is not QSample but " << inputWS->getSpecialCoordinateSystem()
76 << ".";
77 throw std::invalid_argument(errmsg.str());
78 }
79
81
82 // Test indexing. Will be delete soon
83 if (false) {
84 Kernel::V3D qsample; // [1.36639,-2.52888,-4.77349]
85 qsample.setX(1.36639);
86 qsample.setY(-2.52888);
87 qsample.setZ(-4.77349);
88
89 std::vector<Kernel::V3D> vec_qsample;
90 vec_qsample.emplace_back(qsample);
91 std::vector<Kernel::V3D> vec_mindex(0);
92
93 convertFromQSampleToHKL(vec_qsample, vec_mindex);
94 g_log.notice() << "[DB Number of returned Miller Index = " << vec_mindex.size() << "\n";
95 g_log.notice() << "[DB] Output HKL = " << vec_mindex[0].toString() << "\n";
96 }
97
98 // Get events information for future processing
99 std::vector<Kernel::V3D> vec_event_qsample;
100 std::vector<signal_t> vec_event_signal;
101 std::vector<detid_t> vec_event_det;
102 exportEvents(inputWS, vec_event_qsample, vec_event_signal, vec_event_det);
103
104 std::string qsamplefilename = getPropertyValue("QSampleFileName");
105 // Abort with an empty string
106 if (!qsamplefilename.empty())
107 saveEventsToFile(qsamplefilename, vec_event_qsample, vec_event_signal, vec_event_det);
108
109 // Convert to HKL
110 std::vector<Kernel::V3D> vec_event_hkl;
111 convertFromQSampleToHKL(vec_event_qsample, vec_event_hkl);
112
113 // Get file name
114 std::string hklfilename = getPropertyValue("HKLFileName");
115 // Abort mission if no file name is given
116 if (!hklfilename.empty())
117 saveEventsToFile(hklfilename, vec_event_hkl, vec_event_signal, vec_event_det);
118
119 // Create output workspace
120 m_outputWS = createHKLMDWorkspace(vec_event_hkl, vec_event_signal, vec_event_det);
121 // Experiment info
122 ExperimentInfo_sptr expinfo = std::make_shared<ExperimentInfo>();
123 expinfo->setInstrument(inputWS->getExperimentInfo(0)->getInstrument());
124 expinfo->mutableRun().setGoniometer(inputWS->getExperimentInfo(0)->run().getGoniometer(), false);
125 expinfo->mutableRun().addProperty("run_number", 1);
126 m_outputWS->addExperimentInfo(expinfo);
127
128 setProperty("OutputWorkspace", m_outputWS);
129}
130
132 std::string peakwsname = getPropertyValue("PeaksWorkspace");
133 if (!peakwsname.empty() && AnalysisDataService::Instance().doesExist(peakwsname)) {
134 // Get from peak works
135 DataObjects::PeaksWorkspace_sptr peakws = getProperty("PeaksWorkspace");
136 m_UB = Kernel::Matrix<double>(peakws->sample().getOrientedLattice().getUB());
137 } else {
138 // Get from ...
139 std::vector<double> ub_array = getProperty("UBMatrix");
140 if (ub_array.size() != 9)
141 throw std::invalid_argument("Input UB matrix must have 9 elements");
142
144 for (size_t i = 0; i < 3; ++i) {
145 for (size_t j = 0; j < 3; ++j) {
146 m_UB[i][j] = ub_array[i * 3 + j];
147 }
148 }
149 }
150}
151
152//--------------------------------------------------------------------------
157void ConvertCWSDMDtoHKL::exportEvents(const IMDEventWorkspace_sptr &mdws, std::vector<Kernel::V3D> &vec_event_qsample,
158 std::vector<signal_t> &vec_event_signal, std::vector<detid_t> &vec_event_det) {
159 // Set the size of the output vectors
160 size_t numevents = mdws->getNEvents();
161 g_log.information() << "Number of events = " << numevents << "\n";
162
163 vec_event_qsample.resize(numevents);
164 vec_event_signal.resize(numevents);
165 vec_event_det.resize(numevents);
166
167 // Go through to get value
168 auto mditer = mdws->createIterator();
169 size_t nextindex = 1;
170 bool scancell = true;
171 size_t currindex = 0;
172 while (scancell) {
173 size_t numevent_cell = mditer->getNumEvents();
174 for (size_t iev = 0; iev < numevent_cell; ++iev) {
175 // Check
176 if (currindex >= vec_event_qsample.size())
177 throw std::runtime_error("Logic error in event size!");
178
179 float tempx = mditer->getInnerPosition(iev, 0);
180 float tempy = mditer->getInnerPosition(iev, 1);
181 float tempz = mditer->getInnerPosition(iev, 2);
182 signal_t signal = mditer->getInnerSignal(iev);
183 detid_t detid = mditer->getInnerDetectorID(iev);
184 Kernel::V3D qsample(tempx, tempy, tempz);
185 vec_event_qsample[currindex] = qsample;
186 vec_event_signal[currindex] = signal;
187 vec_event_det[currindex] = detid;
188
189 ++currindex;
190 }
191
192 // Advance to next cell
193 if (mditer->next()) {
194 // advance to next cell
195 mditer->jumpTo(nextindex);
196 ++nextindex;
197 } else {
198 // break the loop
199 scancell = false;
200 }
201 }
202}
203
204//--------------------------------------------------------------------------
207void ConvertCWSDMDtoHKL::saveMDToFile(const std::vector<std::vector<coord_t>> &vecEventQsample,
208 const std::vector<float> &vecEventSignal) {
209 // Get file name
210 std::string filename = getPropertyValue("QSampleFileName");
211
212 // Abort with an empty string
213 if (filename.empty())
214 return;
215 if (vecEventQsample.size() != vecEventSignal.size())
216 throw std::runtime_error("Input vectors of Q-sample and signal have different sizes.");
217
218 // Write to file
219 std::ofstream ofile;
220 ofile.open(filename.c_str());
221
222 size_t numevents = vecEventQsample.size();
223 for (size_t i = 0; i < numevents; ++i) {
224 ofile << vecEventQsample[i][0] << ", " << vecEventQsample[i][1] << ", " << vecEventQsample[i][2] << ", "
225 << vecEventSignal[i] << "\n";
226 }
227 ofile.close();
228}
229
230//--------------------------------------------------------------------------
233void ConvertCWSDMDtoHKL::saveEventsToFile(const std::string &filename, std::vector<Kernel::V3D> &vecEventPos,
234 const std::vector<signal_t> &vecEventSignal,
235 const std::vector<detid_t> &vecEventDetid) {
236 // Check
237 if (vecEventDetid.size() != vecEventPos.size() || vecEventPos.size() != vecEventSignal.size())
238 throw std::invalid_argument("Input vectors for HKL, signal and detector ID have different size.");
239
240 std::ofstream ofile;
241 ofile.open(filename.c_str());
242
243 for (size_t i = 0; i < vecEventPos.size(); ++i) {
244 ofile << vecEventPos[i].X() << ", " << vecEventPos[i].Y() << ", " << vecEventPos[i].Z() << ", " << vecEventSignal[i]
245 << ", " << vecEventDetid[i] << "\n";
246 }
247 ofile.close();
248}
249
250//--------------------------------------------------------------------------
253void ConvertCWSDMDtoHKL::convertFromQSampleToHKL(const std::vector<V3D> &q_vectors, std::vector<V3D> &miller_indices) {
254
255 Matrix<double> tempUB(m_UB);
256
257 int original_indexed = 0;
258 double original_error = 0;
259 double tolerance = 0.55; // to make sure no output is invalid
260 original_indexed =
261 IndexingUtils::CalculateMillerIndices(tempUB, q_vectors, tolerance, miller_indices, original_error);
262
263 g_log.notice() << "[DB] " << original_indexed << " peaks are indexed."
264 << "\n";
265}
266
267//----------------------------------------------------------------------------------------------
273 const std::vector<signal_t> &vec_signal,
274 const std::vector<detid_t> &vec_detid) {
275 // Check
276 if (vec_hkl.size() != vec_signal.size() || vec_signal.size() != vec_detid.size())
277 throw std::invalid_argument("Input vectors for HKL, signal and detector "
278 "IDs are of different size!");
279
280 // Create workspace in Q_sample with dimenion as 3
281 size_t nDimension = 3;
282 IMDEventWorkspace_sptr mdws = MDEventFactory::CreateMDWorkspace(nDimension, "MDEvent");
283
284 // Extract Dimensions and add to the output workspace.
285 std::vector<std::string> vec_ID(3);
286 vec_ID[0] = "H";
287 vec_ID[1] = "K";
288 vec_ID[2] = "L";
289
290 std::vector<std::string> dimensionNames(3);
291 dimensionNames[0] = "H";
292 dimensionNames[1] = "K";
293 dimensionNames[2] = "L";
294
296
297 // Add dimensions
298 std::vector<double> m_extentMins(3);
299 std::vector<double> m_extentMaxs(3);
300 std::vector<size_t> m_numBins(3, 100);
301 getRange(vec_hkl, m_extentMins, m_extentMaxs);
302
303 // Get MDFrame of HKL type with RLU
304 auto unitFactory = makeMDUnitFactoryChain();
305 auto unit = unitFactory->create(Units::Symbol::RLU.ascii());
306 Mantid::Geometry::HKL frame(unit);
307
308 for (size_t i = 0; i < nDimension; ++i) {
309 std::string id = vec_ID[i];
310 std::string name = dimensionNames[i];
311 // std::string units = "A^-1";
313 id, name, frame, static_cast<coord_t>(m_extentMins[i]), static_cast<coord_t>(m_extentMaxs[i]), m_numBins[i])));
314 }
315
316 // Set coordinate system
317 mdws->setCoordinateSystem(coordinateSystem);
318
319 // Creates a new instance of the MDEventInserter to output workspace
320 MDEventWorkspace<MDEvent<3>, 3>::sptr mdws_mdevt_3 = std::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3>>(mdws);
321 MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> inserter(mdws_mdevt_3);
322
323 // Go though each spectrum to conver to MDEvent
324 for (size_t iq = 0; iq < vec_hkl.size(); ++iq) {
325 Kernel::V3D hkl = vec_hkl[iq];
326 std::vector<Mantid::coord_t> millerindex(3);
327 millerindex[0] = static_cast<float>(hkl.X());
328 millerindex[1] = static_cast<float>(hkl.Y());
329 millerindex[2] = static_cast<float>(hkl.Z());
330
331 signal_t signal = vec_signal[iq];
332 signal_t error = std::sqrt(signal);
333 uint16_t runnumber = 1;
334 detid_t detid = vec_detid[iq];
335
336 // Insert
337 inserter.insertMDEvent(static_cast<float>(signal), static_cast<float>(error * error),
338 static_cast<uint16_t>(runnumber), 0, detid, millerindex.data());
339 }
340
341 return mdws;
342}
343
344void ConvertCWSDMDtoHKL::getRange(const std::vector<Kernel::V3D> &vec_hkl, std::vector<double> &extentMins,
345 std::vector<double> &extentMaxs) {
346 assert(extentMins.size() == 3);
347 assert(extentMaxs.size() == 3);
348
349 for (size_t i = 0; i < 3; ++i) {
350 double minvalue = vec_hkl[0][i];
351 double maxvalue = vec_hkl[0][i];
352 for (size_t j = 1; j < vec_hkl.size(); ++j) {
353 double thisvalue = vec_hkl[j][i];
354 if (thisvalue < minvalue)
355 minvalue = thisvalue;
356 else if (thisvalue > maxvalue)
357 maxvalue = thisvalue;
358 }
359 extentMins[i] = minvalue;
360 extentMaxs[i] = maxvalue;
361 }
362}
363
364} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double error
Definition: IndexPeaks.cpp:133
double tolerance
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
@ OptionalSave
to specify a file to write to but an empty string is
Definition: FileProperty.h:50
A property class for workspaces.
static API::IMDEventWorkspace_sptr CreateMDWorkspace(size_t nd, const std::string &eventType="MDLeanEvent", const Mantid::API::MDNormalization &preferredNormalization=Mantid::API::MDNormalization::VolumeNormalization, const Mantid::API::MDNormalization &preferredNormalizationHisto=Mantid::API::MDNormalization::VolumeNormalization)
Create a MDEventWorkspace of the given type.
MDEventInserter : Helper class that provides a generic interface for adding events to an MDWorkspace ...
Templated class for the multi-dimensional event workspace.
HKL : HKL MDFrame.
Definition: HKL.h:21
static int CalculateMillerIndices(const Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &q_vectors, double tolerance, std::vector< Kernel::V3D > &miller_indices, double &ave_error)
Given a UB, get list of Miller indices for specifed Qs and tolerance.
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.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Numerical Matrix class.
Definition: Matrix.h:42
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
static const UnitLabel RLU
Reciprocal lattice units.
Class for 3D vectors.
Definition: V3D.h:34
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
void setZ(const double zz) noexcept
Set is z position.
Definition: V3D.h:230
void setX(const double xx) noexcept
Set is x position.
Definition: V3D.h:218
void setY(const double yy) noexcept
Set is y position.
Definition: V3D.h:224
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
ConvertCWSDMDtoHKL : TODO: DESCRIPTION.
void saveEventsToFile(const std::string &filename, std::vector< Kernel::V3D > &vecEventPos, const std::vector< signal_t > &vecEventSignal, const std::vector< detid_t > &vecEventDetid)
Save HKL to file for 3D visualization.
const std::string name() const override
Algorithm's name.
void getRange(const std::vector< Kernel::V3D > &vec_hkl, std::vector< double > &extentMins, std::vector< double > &extentMaxs)
API::IMDEventWorkspace_sptr createHKLMDWorkspace(const std::vector< Kernel::V3D > &vec_hkl, const std::vector< signal_t > &vec_signal, const std::vector< detid_t > &vec_detid)
Create output workspace.
void saveMDToFile(const std::vector< std::vector< coord_t > > &vecEventQsample, const std::vector< float > &vecEventSignal)
Save Q-sample to file.
API::IMDEventWorkspace_sptr m_outputWS
void convertFromQSampleToHKL(const std::vector< Kernel::V3D > &q_vectors, std::vector< Kernel::V3D > &miller_indices)
Convert from Q-sample to HKL.
void exportEvents(const API::IMDEventWorkspace_sptr &mdws, std::vector< Kernel::V3D > &vec_event_qsample, std::vector< signal_t > &vec_event_signal, std::vector< detid_t > &vec_event_det)
Export events from an MDEventWorkspace for future processing It is a convenient algorithm if number o...
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< ExperimentInfo > ExperimentInfo_sptr
Shared pointer to ExperimentInfo.
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< MDHistoDimension > MDHistoDimension_sptr
Shared pointer to a MDHistoDimension.
SpecialCoordinateSystem
Special coordinate systems for Q3D.
MDUnitFactory_uptr MANTID_KERNEL_DLL makeMDUnitFactoryChain()
Convience method. Pre-constructed builder chain.
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54