Mantid
Loading...
Searching...
No Matches
LoadPDFgetNFile.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/Axis.h"
13#include "MantidKernel/Unit.h"
15
16#include <boost/algorithm/string.hpp>
17#include <boost/algorithm/string/split.hpp>
18#include <boost/algorithm/string/trim.hpp>
19#include <fstream>
20#include <iomanip>
21
22using namespace Mantid;
23using namespace Mantid::Kernel;
24using namespace Mantid::API;
25using namespace Mantid::DataObjects;
26using namespace std;
27using namespace boost;
28
29// FIXME Add label and unit to output workspace
30// FIXME Consider to output multiple workspaces if there are multiple column
31// data (X, Y1, E1, Y2, E2)
32
33namespace Mantid::DataHandling {
34
36
37//----------------------------------------------------------------------------------------------
44int LoadPDFgetNFile::confidence(Kernel::FileDescriptor &descriptor) const {
45 // check the file extension
46 const std::string &extn = descriptor.extension();
47 // Only allow known file extensions
48 if (extn != "sq" && extn != "sqa" && extn != "sqb" && extn != "gr" && extn != "ain" && extn != "braw" &&
49 extn != "bsmo") {
50 return 0;
51 }
52
53 if (!descriptor.isAscii())
54 return 0;
55
56 auto &file = descriptor.data();
57 std::string str;
58 std::getline(file, str); // workspace title first line
59 while (!file.eof()) {
60 std::getline(file, str);
61 if (startsWith(str, "#L")) {
62 return 80;
63 }
64 }
65
66 return 0;
67}
68
69//----------------------------------------------------------------------------------------------
73 const std::vector<std::string> exts{".sq", ".sqa", ".sqb", ".gr", ".ain", ".braw", ".bsmo"};
74 auto fileproperty =
75 std::make_unique<FileProperty>("Filename", "", FileProperty::Load, exts, Kernel::Direction::Input);
76 this->declareProperty(std::move(fileproperty), "The input filename of the stored data");
77
78 // auto wsproperty = new WorkspaceProperty<Workspace2D>("OutputWorkspace",
79 // "Anonymous", Kernel::Direction::Output);
80 // this->declareProperty(wsproperty, "Name of output workspace. ");
81 declareProperty(std::make_unique<API::WorkspaceProperty<>>("OutputWorkspace", "", Kernel::Direction::Output),
82 "Workspace name to load into.");
83}
84
85//----------------------------------------------------------------------------------------------
89 // 1. Parse input file
90 std::string inpfilename = getProperty("Filename");
91
92 parseDataFile(inpfilename);
93
94 // 2. Generate output workspace
96
97 setProperty("OutputWorkspace", outWS);
98}
99
100//----------------------------------------------------------------------------------------------
105void LoadPDFgetNFile::parseDataFile(const std::string &filename) {
106 // 1. Open file
107 std::ifstream ifile;
108 ifile.open((filename.c_str()));
109
110 if (!ifile.is_open()) {
111 stringstream errmsg;
112 errmsg << "Unable to open file " << filename << ". Quit!";
113 g_log.error() << errmsg.str() << '\n';
114 throw std::runtime_error(errmsg.str());
115 } else {
116 g_log.notice() << "Open PDFgetN File " << filename << '\n';
117 }
118
119 // 2. Parse
120 bool readdata = false;
121
122 char line[256];
123 while (ifile.getline(line, 256)) {
124 string sline(line);
125
126 if (!readdata && startsWith(sline, "#L")) {
127 // a) Find header line for the data segment in the file
128 parseColumnNameLine(sline);
129 readdata = true;
130
131 // Set up the data structure
132 size_t numcols = mColumnNames.size();
133 for (size_t i = 0; i < numcols; ++i) {
134 std::vector<double> tempvec;
135 mData.emplace_back(tempvec);
136 }
137
138 } else if (readdata) {
139 // b) Parse data
140 parseDataLine(sline);
141 } else {
142 // c) Do nothing
143 ;
144 }
145 } // ENDWHILE
146
147 if (!readdata) {
148 stringstream errmsg;
149 errmsg << "Unable to find a line staring with #L as the indicator of data "
150 "segment. ";
151 g_log.error() << errmsg.str() << '\n';
152 throw std::runtime_error(errmsg.str());
153 }
154}
155
156//----------------------------------------------------------------------------------------------
159bool LoadPDFgetNFile::startsWith(const std::string &s, const std::string &header) const {
160 bool answer = true;
161
162 if (s.size() < header.size()) {
163 answer = false;
164 } else {
165 size_t numchars = header.size();
166 for (size_t i = 0; i < numchars; ++i) {
167 char c0 = s[i];
168 char c1 = header[i];
169 if (c0 != c1) {
170 answer = false;
171 break;
172 }
173 }
174 }
175
176 return answer;
177}
178
182 // 1. Split
183 vector<string> terms;
184 boost::split(terms, line, is_any_of(" \t\n"), token_compress_on);
185
186 // 2. Validate
187 if (terms.empty()) {
188 throw std::runtime_error("There is nothing in the input line!");
189 }
190
191 string header = terms[0];
192 if (header != "#L") {
193 stringstream errmsg;
194 errmsg << "Expecting header as #L. Input line has header as " << header << ". Unable to proceed. ";
195 g_log.error() << errmsg.str() << '\n';
196 throw std::runtime_error(errmsg.str());
197 }
198
199 // 3. Parse
200 size_t numcols = terms.size() - 1;
201 stringstream msgss;
202 msgss << "Column Names: ";
203 for (size_t i = 0; i < numcols; ++i) {
204 this->mColumnNames.emplace_back(terms[i + 1]);
205 msgss << setw(-3) << i << ": " << setw(-10) << mColumnNames[i];
206 }
207 g_log.information() << msgss.str() << '\n';
208}
209
213 // 1. Trim (stripg) and Split
214 boost::trim(line);
215 vector<string> terms;
216 boost::split(terms, line, is_any_of(" \t\n"), token_compress_on);
217
218 // 2. Validate
219 size_t numcols = mData.size();
220 if (line[0] == '#') {
221 // Comment/information line to indicate the start of another section of data
222 return;
223 } else if (terms.size() != numcols) {
224 // Data line with incorrect number of columns
225 stringstream warnss;
226 warnss << "Line (" << line << ") has incorrect number of columns other than " << numcols << " as expected. ";
227 g_log.warning(warnss.str());
228 return;
229 }
230
231 // 3. Parse
232 for (size_t i = 0; i < numcols; ++i) {
233 string temps = terms[i];
234 double tempvalue;
235 if (temps == "NaN") {
236 // FIXME: Need to discuss with Peter about how to treat NaN value
237 // tempvalue = DBL_MAX-1.0;
238 tempvalue = 0.0;
239 } else if (temps == "-NaN") {
240 // tempvalue = -DBL_MAX+1.0;
241 // FIXME: Need to discuss with Peter about how to treat NaN value
242 tempvalue = 0.0;
243 } else {
244 tempvalue = std::stod(temps);
245 }
246
247 mData[i].emplace_back(tempvalue);
248 }
249}
250
251//----------------------------------------------------------------------------------------------
253 // 1. Set X
254 string xcolname = mColumnNames[0];
255
256 if (xcolname == "Q") {
257 string unit = "MomentumTransfer";
258 ws->getAxis(0)->setUnit(unit);
259 } else if (xcolname == "r") {
260 ws->getAxis(0)->setUnit("AtomicDistance");
261 } else {
262 stringstream errss;
263 errss << "X axis " << xcolname << " is not supported for unit. \n";
264 g_log.warning() << errss.str() << '\n';
265 }
266
267 // 2. Set Y
268 string ycolname = mColumnNames[1];
269 string ylabel;
270 if (ycolname == "G(r)") {
271 ylabel = "PDF";
272 } else if (ycolname == "S") {
273 ylabel = "S";
274 } else {
275 ylabel = "Intensity";
276 }
277 ws->setYUnitLabel(ylabel);
278}
279
280size_t calcVecSize(const std::vector<double> &data0, std::vector<size_t> &numptsvec, size_t &numsets, bool xascend) {
281 size_t vecsize = 1;
282 auto prex = data0[0];
283 for (size_t i = 1; i < data0.size(); ++i) {
284 double curx = data0[i];
285 if (((xascend) && (curx < prex)) || ((!xascend) && (curx > prex))) {
286 // X in ascending order and hit the end of one set of data
287 // X in descending order and hit the end of one set of data
288 // Record the current data set information and start the next data set
289 numsets += 1;
290 numptsvec.emplace_back(vecsize);
291 vecsize = 1;
292 } else {
293 // In the middle of a set of data
294 ++vecsize;
295 }
296 // Loop variable udpate
297 prex = curx;
298 } // ENDFOR
299
300 return vecsize;
301}
302
303void LoadPDFgetNFile::checkSameSize(const std::vector<size_t> &numptsvec, size_t numsets) {
304 bool samesize = true;
305 for (size_t i = 0; i < numsets; ++i) {
306 if (i > 0) {
307 if (numptsvec[i] != numptsvec[i - 1]) {
308 samesize = false;
309 }
310 }
311 g_log.information() << "Set " << i << ": Number of Points = " << numptsvec[i] << '\n';
312 }
313 if (!samesize) {
314 stringstream errmsg;
315 errmsg << "Multiple bank (number of banks = " << numsets
316 << ") have different size of data array. Unable to handle this "
317 "situation.";
318 g_log.error() << errmsg.str() << '\n';
319 throw std::runtime_error(errmsg.str());
320 }
321}
322
327 // 0. Check
328 if (mData.empty()) {
329 throw runtime_error("Data set has not been initialized. Quit!");
330 }
331
332 // 1. Figure out direction of X and number of data set
333 bool xascend = true;
334 if (mData.size() >= 2 && mData[0][1] < mData[0][0]) {
335 xascend = false;
336 }
337
338 size_t numsets = 0;
339 vector<size_t> numptsvec;
340 size_t arraysize = mData[0].size();
341 if (arraysize <= 1) {
342 throw runtime_error("Number of columns in data is less and equal to 1. It "
343 "is unphysically too small.");
344 }
345
346 // Record the last data set information
347 ++numsets;
348 numptsvec.emplace_back(calcVecSize(mData[0], numptsvec, numsets, xascend));
349
350 checkSameSize(numptsvec, numsets);
351
352 size_t size = numptsvec[0];
353
354 // 2. Generate workspace2D object and set the unit
355 outWS = std::dynamic_pointer_cast<Workspace2D>(
356 API::WorkspaceFactory::Instance().create("Workspace2D", numsets, size, size));
357
358 setUnit(outWS);
359
360 // 3. Set number
361 size_t numspec = outWS->getNumberHistograms();
362 for (size_t i = 0; i < numspec; ++i) {
363 auto &X = outWS->mutableX(i);
364 auto &Y = outWS->mutableY(i);
365 auto &E = outWS->mutableE(i);
366
367 size_t baseindex = i * size;
368 for (size_t j = 0; j < size; ++j) {
369 size_t index;
370 if (xascend)
371 index = j;
372 else
373 index = (size - 1) - j;
374
375 X[index] = mData[0][baseindex + j];
376 Y[index] = mData[1][baseindex + j];
377 E[index] = mData[2][baseindex + j];
378 }
379 }
380}
381
382} // namespace Mantid::DataHandling
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define DECLARE_FILELOADER_ALGORITHM(classname)
DECLARE_FILELOADER_ALGORITHM should be used in place of the standard DECLARE_ALGORITHM macro when wri...
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
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
A property class for workspaces.
LoadPDFgetNFile : TODO: DESCRIPTION.
void parseDataFile(const std::string &filename)
Parse PDFgetN data file.
void init() override
Implement abstract Algorithm methods.
void setUnit(const DataObjects::Workspace2D_sptr &ws)
Set X and Y axis unit and lebel.
void generateDataWorkspace()
Generate output workspace.
std::vector< std::string > mColumnNames
Names of the columns of the data.
DataObjects::Workspace2D_sptr outWS
Output data workspace.
void checkSameSize(const std::vector< size_t > &numptsvec, size_t numsets)
void parseDataLine(std::string line)
Parse data line.
bool startsWith(const std::string &s, const std::string &header) const
Check whether a string starts from a specified sub-string.
void parseColumnNameLine(std::string line)
Parse column name line staring with #L.
std::vector< std::vector< double > > mData
Data structure to hold input: Size = Number of columns in input file.
void exec() override
Implement abstract Algorithm methods.
Defines a wrapper around an open file.
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 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
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Kernel::Logger g_log("ExperimentInfo")
static logger object
bool startsWith(const string &str, const string &prefix)
Returns true if str starts with prefix.
size_t calcVecSize(const std::vector< double > &data0, std::vector< size_t > &numptsvec, size_t &numsets, bool xascend)
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
Helper class which provides the Collimation Length for SANS instruments.
Definition: NDArray.h:49
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54