Mantid
Loading...
Searching...
No Matches
LoadGSS.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 +
7//---------------------------------------------------
8// Includes
9//---------------------------------------------------
11#include "MantidAPI/Axis.h"
13#include "MantidAPI/ISpectrum.h"
22
23#include <Poco/File.h>
24#include <boost/regex.hpp>
25#include <fstream>
26#include <sstream>
27#include <string>
28
29using namespace Mantid::DataHandling;
30using namespace Mantid::API;
31using namespace Mantid::HistogramData;
32using namespace Mantid::Kernel;
33
34namespace Mantid::DataHandling {
35
37
38namespace { // anonymous namespace
39const boost::regex DET_POS_REG_EXP{"^#.+flight path\\s+([0-9.]+).+"
40 "tth\\s+([0-9.]+).+"
41 "DIFC\\s+([0-9.]+)"};
42const boost::regex L1_REG_EXP{"^#.+flight path\\s+([0-9.]+)\\s*m"};
43} // end of anonymous namespace
44
45//----------------------------------------------------------------------------------------------
52
53 if (!descriptor.isAscii() || descriptor.extension() == ".tar")
54 return 0;
55
56 std::string str;
57 std::istream &file = descriptor.data();
58 std::getline(file, str); // workspace title first line
59 while (!file.eof()) {
60 std::getline(file, str);
61 // Skip over empty and comment lines, as well as those coming from files
62 // saved with the 'ExtendedHeader' option
63 if (str.empty() || str[0] == '#' || str.compare(0, 8, "Monitor:") == 0) {
64 continue;
65 }
66 if (str.compare(0, 4, "BANK") == 0 &&
67 (str.find("RALF") != std::string::npos || str.find("SLOG") != std::string::npos) &&
68 (str.find("FXYE") != std::string::npos)) {
69 return 80;
70 }
71 }
72 return 0;
73}
74
75//----------------------------------------------------------------------------------------------
79 const std::vector<std::string> exts{".gsa", ".gss", ".gda", ".txt"};
80 declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Load, exts),
81 "The input filename of the stored data");
82
83 declareProperty(std::make_unique<API::WorkspaceProperty<>>("OutputWorkspace", "", Kernel::Direction::Output),
84 "Workspace name to load into.");
85
86 declareProperty("UseBankIDasSpectrumNumber", false,
87 "If true, spectrum number corresponding to each bank is to "
88 "be its bank ID. ");
89}
90
91//----------------------------------------------------------------------------------------------
95 // Process input parameters
96 std::string filename = getPropertyValue("Filename");
97
98 bool useBankAsSpectrum = getProperty("UseBankIDasSpectrumNumber");
99
100 MatrixWorkspace_sptr outputWorkspace = loadGSASFile(filename, useBankAsSpectrum);
101
102 setProperty("OutputWorkspace", outputWorkspace);
103}
104
105//----------------------------------------------------------------------------------------------
108API::MatrixWorkspace_sptr LoadGSS::loadGSASFile(const std::string &filename, bool useBankAsSpectrum) {
109 // Vectors for detector information
110 double primaryflightpath = -1;
111 std::vector<double> twothetas;
112 std::vector<double> difcs;
113 std::vector<double> totalflightpaths;
114 std::vector<int> detectorIDs;
115
116 // Vectors to store data
117 std::vector<HistogramData::BinEdges> gsasDataX;
118 std::vector<HistogramData::Counts> gsasDataY;
119 std::vector<HistogramData::CountStandardDeviations> gsasDataE;
120
121 std::vector<double> vecX, vecY, vecE;
122
123 // progress
124 std::unique_ptr<Progress> prog = nullptr;
125
126 // Parameters for reading file
127 char currentLine[256];
128 std::string wsTitle;
129 std::string slogTitle;
130 std::string instrumentname = "Generic";
131 char filetype = 'x';
132
133 // Gather data
134 std::ifstream input(filename.c_str(), std::ios_base::in);
135 if (!input.is_open()) {
136 // throw exception if file cannot be opened
137 std::stringstream errss;
138 errss << "Unable to open GSAS file " << filename;
139 throw std::runtime_error(errss.str());
140 }
141
142 // First line: Title
143 if (!input.eof()) {
144 // Get workspace title (should be first line or 2nd line for SLOG)
145 input.getline(currentLine, 256);
146 wsTitle = currentLine;
147 } else {
148 throw std::runtime_error("File is empty");
149 }
150
151 // Loop all the lines
152 bool isOutOfHead = false;
153 bool slogtitleset = false;
154 bool multiplybybinwidth = false;
155 int nSpec = 0;
156 bool calslogx0 = true;
157 double bc4 = 0;
158 double bc3 = 0;
159
160 while (!input.eof() && input.getline(currentLine, 256)) {
161 // Initialize progress after NSpec is imported
162 if (nSpec != 0 && !prog) {
163 prog = std::make_unique<Progress>(this, 0.0, 1.0, nSpec);
164 }
165
166 // Set flag to test SLOG
167 if (!slogtitleset) {
168 slogTitle = currentLine;
169 slogtitleset = true;
170 }
171
172 if (currentLine[0] == '\n' || currentLine[0] == '#') {
173 // Comment/information line
174 std::string key1, key2;
175 std::istringstream inputLine(currentLine, std::ios::in);
176 inputLine.ignore(256, ' ');
177 inputLine >> key1 >> key2;
178
179 if (key2 == "Histograms") {
180 // NSpec (Format: 'nspec HISTOGRAM')
181 nSpec = std::stoi(key1);
182 g_log.information() << "Histogram Line: " << key1 << " nSpec = " << nSpec << "\n";
183 } else if (key1 == "Instrument:") {
184 // Instrument (Format: 'Instrument XXXX')
185 instrumentname = key2;
186 g_log.information() << "Instrument : " << key2 << "\n";
187 } else if (key1 == "with") {
188 // Multiply by bin width: (Format: 'with multiplied')
189 std::string s1;
190 inputLine >> s1;
191 if (s1 == "multiplied") {
192 multiplybybinwidth = true;
193 g_log.information() << "Y is multiplied by bin width\n";
194 } else {
195 g_log.warning() << "In line '" << currentLine << "', key word " << s1 << " is not allowed!\n";
196 }
197 } else if (key1 == "Primary") {
198 // Primary flight path ...
199 boost::smatch result;
200 // Have to force a copy of the input or the stack gets corrupted
201 // on MSVC when inputLine.str() falls out of scope which then
202 // corrupts the value in result
203 const std::string line = inputLine.str();
204 if (boost::regex_search(line, result, L1_REG_EXP) && result.size() == 2) {
205 primaryflightpath = std::stod(std::string(result[1]));
206
207 } else {
208 std::stringstream msg;
209 msg << "Failed to parse primary flight path from line \"" << inputLine.str() << "\"";
210 g_log.warning(msg.str());
211 }
212
213 std::stringstream msg;
214 msg << "L1 = " << primaryflightpath;
215 g_log.information(msg.str());
216 } else if (key1 == "Total") {
217 // Total flight path .... .... including total flying path, difc and
218 // 2theta of 1 bank
219
220 double totalpath(0.f);
221 double tth(0.f);
222 double difc(0.f);
223
224 boost::smatch result;
225 const std::string line = inputLine.str();
226 if (boost::regex_search(line, result, DET_POS_REG_EXP) && result.size() == 4) {
227 totalpath = std::stod(std::string(result[1]));
228 tth = std::stod(std::string(result[2]));
229 difc = std::stod(std::string(result[3]));
230 } else {
231 std::stringstream msg;
232 msg << "Failed to parse position from line \"" << inputLine.str() << "\"";
233 g_log.warning(msg.str());
234 }
235
236 totalflightpaths.emplace_back(totalpath);
237 twothetas.emplace_back(tth);
238 difcs.emplace_back(difc);
239
240 std::stringstream msg;
241 msg << "Bank " << difcs.size() - 1 << ": Total flight path = " << totalpath << " 2Theta = " << tth
242 << " DIFC = " << difc;
243 g_log.information(msg.str());
244 } // if keys....
245
246 } // ENDIF for Line with #
247 else if (currentLine[0] == 'B') {
248 // Line start with Bank including file format, X0 information and etc.
249 isOutOfHead = true;
250
251 // If there is, Save the previous to array and initialize new MantiVec for
252 // (X, Y, E)
253 if (!vecX.empty()) {
254 gsasDataX.emplace_back(std::move(vecX));
255 gsasDataY.emplace_back(std::move(vecY));
256 gsasDataE.emplace_back(std::move(vecE));
257 vecX.clear();
258 vecY.clear();
259 vecE.clear();
260
261 if (prog != nullptr)
262 prog->report();
263 }
264
265 // Parse the bank line in format
266 // RALF: BANK <SpectraNo> <NBins> <NBins> RALF <BC1> <BC2> <BC1> <BC4>
267 // SLOG: BANK <SpectraNo> <NBins> <NBins> SLOG <BC1> <BC2> <BC3> 0>
268 // where,
269 // BC1 = X[0] * 32
270 // BC2 = X[1] * 32 - BC1
271 // BC4 = ( X[1] - X[0] ) / X[0]
272
273 int specno, nbin1, nbin2;
274 std::istringstream inputLine(currentLine, std::ios::in);
275
276 double bc1 = 0;
277 double bc2 = 0;
278
279 inputLine.ignore(256, 'K');
280 std::string filetypestring;
281
282 inputLine >> specno >> nbin1 >> nbin2 >> filetypestring;
283 g_log.debug() << "Bank: " << specno << " filetypestring = " << filetypestring << '\n';
284
285 detectorIDs.emplace_back(specno);
286
287 if (filetypestring[0] == 'S') {
288 // SLOG
289 filetype = 's';
290 inputLine >> bc1 >> bc2 >> bc3 >> bc4;
291 } else if (filetypestring[0] == 'R') {
292 // RALF
293 filetype = 'r';
294 inputLine >> bc1 >> bc2 >> bc1 >> bc4;
295 } else {
296 g_log.error() << "Unsupported GSAS File Type: " << filetypestring << "\n";
297 throw Exception::FileError("Not a GSAS file", filename);
298 }
299
300 // Determine x0
301 if (filetype == 'r') {
302 double x0 = bc1 / 32;
303 g_log.debug() << "RALF: x0 = " << x0 << " bc4 = " << bc4 << '\n';
304 vecX.emplace_back(x0);
305 } else {
306 // Cannot calculate x0, turn on the flag
307 calslogx0 = true;
308 }
309 } // Line with B
310 else if (isOutOfHead) {
311 // Parse data line
312 double xValue;
313 double yValue;
314 double eValue;
315
316 double xPrev;
317
318 // * Get previous X value
319 if (!vecX.empty()) {
320 xPrev = vecX.back();
321 } else if (filetype == 'r') {
322 // Except if RALF
323 throw Mantid::Kernel::Exception::NotImplementedError("LoadGSS: File was not in expected format.");
324 } else {
325 xPrev = -0.0;
326 }
327
328 // It is different for the definition of X, Y, Z in SLOG and RALF format
329 if (filetype == 'r') {
330 // RALF
331 // LoadGSS produces overlapping columns for some datasets, due to
332 // std::setw
333 // For this reason we need to read the column values as string and then
334 // convert to double
335 {
336 std::string str(currentLine, 15);
337 std::istringstream istr(str);
338 istr >> xValue;
339 }
340 {
341 std::string str(currentLine + 15, 18);
342 std::istringstream istr(str);
343 istr >> yValue;
344 }
345 {
346 std::string str(currentLine + 15 + 18, 18);
347 std::istringstream istr(str);
348 istr >> eValue;
349 }
350
351 xValue = (2 * xValue) - xPrev;
352
353 } else if (filetype == 's') {
354 // SLOG
355 std::istringstream inputLine(currentLine, std::ios::in);
356 inputLine >> xValue >> yValue >> eValue;
357 if (calslogx0) {
358 // calculation of x0 must use the x'[0]
359 g_log.debug() << "x'_0 = " << xValue << " bc3 = " << bc3 << '\n';
360
361 double x0 = 2 * xValue / (bc3 + 2.0);
362 vecX.emplace_back(x0);
363 xPrev = x0;
364 g_log.debug() << "SLOG: x0 = " << x0 << '\n';
365 calslogx0 = false;
366 }
367
368 xValue = (2 * xValue) - xPrev;
369 } else {
370 g_log.error() << "Unsupported GSAS File Type: " << filetype << "\n";
371 throw Exception::FileError("Not a GSAS file", filename);
372 }
373
374 if (multiplybybinwidth) {
375 yValue = yValue / (xValue - xPrev);
376 eValue = eValue / (xValue - xPrev);
377 }
378
379 // store read in data (x, y, e) to vector
380 vecX.emplace_back(xValue);
381 vecY.emplace_back(yValue);
382 vecE.emplace_back(eValue);
383 } // Date Line
384 else {
385 g_log.warning() << "Line not defined: " << currentLine << '\n';
386 }
387 } // ENDWHILE of reading all lines
388
389 // Get the sizes before using std::move
390 auto nHist(static_cast<int>(gsasDataX.size()));
391 auto xWidth(static_cast<int>(vecX.size()));
392 auto yWidth(static_cast<int>(vecY.size()));
393
394 // Push the vectors (X, Y, E) of the last bank to gsasData
395 if (!vecX.empty()) { // Put final spectra into data
396 gsasDataX.emplace_back(std::move(vecX));
397 gsasDataY.emplace_back(std::move(vecY));
398 gsasDataE.emplace_back(std::move(vecE));
399 ++nHist;
400 }
401 input.close();
402
403 //********************************************************************************************
404 // Construct the workspace for GSS data
405 //********************************************************************************************
406
407 // Create workspace & GSS Files data is always in TOF
408
409 MatrixWorkspace_sptr outputWorkspace = std::dynamic_pointer_cast<MatrixWorkspace>(
410 WorkspaceFactory::Instance().create("Workspace2D", nHist, xWidth, yWidth));
411 outputWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
412
413 // set workspace title
414 if (filetype == 'r')
415 outputWorkspace->setTitle(wsTitle);
416 else
417 outputWorkspace->setTitle(slogTitle);
418
419 // put data from constructed histograms into outputWorkspace
420 if (detectorIDs.size() != static_cast<size_t>(nHist)) {
421 // File error is found
422 std::ostringstream mess("");
423 mess << "Number of spectra (" << detectorIDs.size() << ") is not equal to number of histograms (" << nHist << ").";
424 throw std::runtime_error(mess.str());
425 }
426
427 for (int i = 0; i < nHist; ++i) {
428 // Move data across
429 outputWorkspace->setHistogram(i, BinEdges(std::move(gsasDataX[i])), Counts(std::move(gsasDataY[i])),
430 CountStandardDeviations(std::move(gsasDataE[i])));
431
432 // Reset spectrum number if
433 if (useBankAsSpectrum) {
434 auto specno = static_cast<specnum_t>(detectorIDs[i]);
435 outputWorkspace->getSpectrum(i).setSpectrumNo(specno);
436 }
437 }
438
439 // build instrument geometry
440 createInstrumentGeometry(outputWorkspace, instrumentname, primaryflightpath, detectorIDs, totalflightpaths, twothetas,
441 difcs);
442
443 return outputWorkspace;
444}
445
446//----------------------------------------------------------------------------------------------
449double LoadGSS::convertToDouble(std::string inputstring) {
450 std::string temps;
451 auto isize = static_cast<int>(inputstring.size());
452 for (int i = 0; i < isize; i++) {
453 char thechar = inputstring[i];
454 if ((thechar <= 'Z' && thechar >= 'A') || (thechar <= 'z' && thechar >= 'a')) {
455 break;
456 } else {
457 temps += thechar;
458 }
459 }
460
461 double rd = std::stod(temps);
462
463 return rd;
464}
465
466//----------------------------------------------------------------------------------------------
469void LoadGSS::createInstrumentGeometry(const MatrixWorkspace_sptr &workspace, const std::string &instrumentname,
470 const double &primaryflightpath, const std::vector<int> &detectorids,
471 const std::vector<double> &totalflightpaths,
472 const std::vector<double> &twothetas, const std::vector<double> &difcs) {
473 // Check Input
474 if (detectorids.size() != totalflightpaths.size() || totalflightpaths.size() != twothetas.size()) {
475 g_log.warning("Cannot create geometry, because the numbers of L2 and Polar "
476 "are not equal.");
477 return;
478 }
479
480 // Debug output
481 std::stringstream dbss;
482 dbss << "L1 = " << primaryflightpath << "\n";
483 for (size_t i = 0; i < detectorids.size(); i++) {
484 dbss << "Detector " << detectorids[i] << " L1+L2 = " << totalflightpaths[i] << " 2Theta = " << twothetas[i]
485 << "\n";
486 }
487 g_log.debug(dbss.str());
488
489 // Create a new instrument and set its name
490 Geometry::Instrument_sptr instrument(new Geometry::Instrument(instrumentname));
491
492 // Add dummy source and samplepos to instrument
493 Geometry::Component *samplepos = new Geometry::Component("Sample", instrument.get());
494 instrument->add(samplepos);
495 instrument->markAsSamplePos(samplepos);
496 samplepos->setPos(0.0, 0.0, 0.0);
497
498 Geometry::ObjComponent *source = new Geometry::ObjComponent("Source", instrument.get());
499 instrument->add(source);
500 instrument->markAsSource(source);
501
502 double l1 = primaryflightpath;
503 source->setPos(0.0, 0.0, -1.0 * l1);
504
505 // Add detectors
506 // The L2 and 2-theta values from Raw file assumed to be relative to sample
507 // position
508 const auto numDetector = static_cast<int>(detectorids.size()); // number of detectors
509 // std::vector<int> detID = detectorids; // detector IDs
510 // std::vector<double> angle = twothetas; // angle between indicent beam and
511 // direction from sample to detector (two-theta)
512
513 // Assumption: detector IDs are in the same order of workspace index
514 for (int i = 0; i < numDetector; ++i) {
515 // a) Create a new detector. Instrument will take ownership of pointer so no
516 // need to delete.
517 Geometry::Detector *detector = new Geometry::Detector("det", detectorids[i], samplepos);
518 Kernel::V3D pos;
519
520 // r is L2
521 double r = totalflightpaths[i] - l1;
522 pos.spherical(r, twothetas[i], 0.0);
523
524 detector->setPos(pos);
525
526 // add copy to instrument, spectrum and mark it
527 auto &spec = workspace->getSpectrum(i);
528 spec.clearDetectorIDs();
529 spec.addDetectorID(detectorids[i]);
530 instrument->add(detector);
531 instrument->markAsDetector(detector);
532
533 } // ENDFOR (i: spectrum)
534 workspace->setInstrument(instrument);
535
536 auto &paramMap = workspace->instrumentParameters();
537 for (size_t i = 0; i < workspace->getNumberHistograms(); i++) {
538 auto detector = workspace->getDetector(i);
539 paramMap.addDouble(detector->getComponentID(), "DIFC", difcs[i]);
540 }
541}
542} // namespace Mantid::DataHandling
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
#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
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
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
A property class for workspaces.
Loads a file as saved by SaveGSS.
Definition: LoadGSS.h:24
void createInstrumentGeometry(const API::MatrixWorkspace_sptr &workspace, const std::string &instrumentname, const double &primaryflightpath, const std::vector< int > &detectorids, const std::vector< double > &totalflightpaths, const std::vector< double > &twothetas, const std::vector< double > &difcs)
Create an instrument geometry.
Definition: LoadGSS.cpp:469
double convertToDouble(std::string inputstring)
Convert a string (value+unit) to double (value)
Definition: LoadGSS.cpp:449
void exec() override
Execution code.
Definition: LoadGSS.cpp:94
int confidence(Kernel::FileDescriptor &descriptor) const override
Returns a confidence value that this algorithm can load a file.
Definition: LoadGSS.cpp:51
void init() override
Initialisation code.
Definition: LoadGSS.cpp:78
API::MatrixWorkspace_sptr loadGSASFile(const std::string &filename, bool useBankAsSpectrum)
Main method to load GSAS.
Definition: LoadGSS.cpp:108
Component is a wrapper for a Component which can modify some of its parameters, e....
Definition: Component.h:41
void setPos(double, double, double) override
Set the IComponent position, x, y, z respective to parent (if present)
Definition: Component.cpp:204
This class represents a detector - i.e.
Definition: Detector.h:30
Base Instrument Class.
Definition: Instrument.h:47
Object Component class, this class brings together the physical attributes of the component to the po...
Definition: ObjComponent.h:33
Records the filename and the description of failure.
Definition: Exception.h:98
Marks code as not implemented yet.
Definition: Exception.h:138
Defines a wrapper around an open file.
static bool isAscii(const std::string &filename, const size_t nbytes=256)
Returns true if the file is considered ascii.
std::istream & data()
Access the open file stream.
const std::string & extension() const
Access the file extension.
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 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...
Class for 3D vectors.
Definition: V3D.h:34
void spherical(const double R, const double theta, const double phi) noexcept
Sets the vector position based on spherical coordinates.
Definition: V3D.cpp:57
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
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.
std::shared_ptr< Instrument > Instrument_sptr
Shared pointer to an instrument object.
int32_t specnum_t
Typedef for a spectrum Number.
Definition: IDTypes.h:16
@ Output
An output workspace.
Definition: Property.h:54