LoadGSS.cpp 18.5 KB
Newer Older
1
2
3
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4
5
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
// SPDX - License - Identifier: GPL - 3.0 +
7
8
9
10
//---------------------------------------------------
// Includes
//---------------------------------------------------
#include "MantidDataHandling/LoadGSS.h"
11
#include "MantidAPI/Axis.h"
12
#include "MantidAPI/FileProperty.h"
13
14
#include "MantidAPI/ISpectrum.h"
#include "MantidAPI/MatrixWorkspace.h"
15
#include "MantidAPI/RegisterFileLoader.h"
16
#include "MantidAPI/WorkspaceFactory.h"
17
#include "MantidGeometry/Instrument.h"
18
19
#include "MantidGeometry/Instrument/CompAssembly.h"
#include "MantidGeometry/Instrument/Component.h"
Dimitar Tasev's avatar
Dimitar Tasev committed
20
#include "MantidGeometry/Instrument/Detector.h"
21
#include "MantidKernel/UnitFactory.h"
22

Dimitar Tasev's avatar
Dimitar Tasev committed
23
#include <Poco/File.h>
24
#include <boost/regex.hpp>
25
#include <fstream>
Karl Palmen's avatar
Karl Palmen committed
26
#include <sstream>
Peterson, Peter's avatar
Peterson, Peter committed
27
#include <string>
28
29
30

using namespace Mantid::DataHandling;
using namespace Mantid::API;
31
using namespace Mantid::HistogramData;
32
33
using namespace Mantid::Kernel;

34
35
36
37
38
namespace Mantid {
namespace DataHandling {

DECLARE_FILELOADER_ALGORITHM(LoadGSS)

39
40
41
42
43
44
45
namespace { // anonymous namespace
const boost::regex DET_POS_REG_EXP{"^#.+flight path\\s+([0-9.]+).+"
                                   "tth\\s+([0-9.]+).+"
                                   "DIFC\\s+([0-9.]+)"};
const boost::regex L1_REG_EXP{"^#.+flight path\\s+([0-9.]+)\\s*m"};
} // end of anonymous namespace

46
47
//----------------------------------------------------------------------------------------------
/** Return the confidence with with this algorithm can load the file
LamarMoore's avatar
LamarMoore committed
48
49
 * @param descriptor A descriptor for the file
 * @returns An integer specifying the confidence level. 0 indicates it will not
50
 * be used
LamarMoore's avatar
LamarMoore committed
51
 */
52
int LoadGSS::confidence(Kernel::FileDescriptor &descriptor) const {
Owen Arnold's avatar
Owen Arnold committed
53
54

  if (!descriptor.isAscii() || descriptor.extension() == ".tar")
Zhou, Wenduo's avatar
Zhou, Wenduo committed
55
56
    return 0;

57
58
59
60
61
62
63
64
65
66
67
  std::string str;
  std::istream &file = descriptor.data();
  std::getline(file, str); // workspace title first line
  while (!file.eof()) {
    std::getline(file, str);
    // Skip over empty and comment lines, as well as those coming from files
    // saved with the 'ExtendedHeader' option
    if (str.empty() || str[0] == '#' || str.compare(0, 8, "Monitor:") == 0) {
      continue;
    }
    if (str.compare(0, 4, "BANK") == 0 &&
68
        (str.find("RALF") != std::string::npos || str.find("SLOG") != std::string::npos) &&
69
70
71
72
73
74
75
76
77
        (str.find("FXYE") != std::string::npos)) {
      return 80;
    }
  }
  return 0;
}

//----------------------------------------------------------------------------------------------
/** Initialise the algorithm
LamarMoore's avatar
LamarMoore committed
78
 */
79
void LoadGSS::init() {
80
  const std::vector<std::string> exts{".gsa", ".gss", ".gda", ".txt"};
81
  declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Load, exts),
82
                  "The input filename of the stored data");
83

84
  declareProperty(std::make_unique<API::WorkspaceProperty<>>("OutputWorkspace", "", Kernel::Direction::Output),
85
86
87
88
89
90
91
92
93
                  "Workspace name to load into.");

  declareProperty("UseBankIDasSpectrumNumber", false,
                  "If true, spectrum number corresponding to each bank is to "
                  "be its bank ID. ");
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm
LamarMoore's avatar
LamarMoore committed
94
 */
95
96
97
98
99
100
void LoadGSS::exec() {
  // Process input parameters
  std::string filename = getPropertyValue("Filename");

  bool useBankAsSpectrum = getProperty("UseBankIDasSpectrumNumber");

101
  MatrixWorkspace_sptr outputWorkspace = loadGSASFile(filename, useBankAsSpectrum);
102
103
104
105
106
107

  setProperty("OutputWorkspace", outputWorkspace);
}

//----------------------------------------------------------------------------------------------
/** Main method to load GSAS file
LamarMoore's avatar
LamarMoore committed
108
 */
109
API::MatrixWorkspace_sptr LoadGSS::loadGSASFile(const std::string &filename, bool useBankAsSpectrum) {
110
111
112
113
114
115
116
117
  // Vectors for detector information
  double primaryflightpath = -1;
  std::vector<double> twothetas;
  std::vector<double> difcs;
  std::vector<double> totalflightpaths;
  std::vector<int> detectorIDs;

  // Vectors to store data
118
119
120
121
  std::vector<HistogramData::BinEdges> gsasDataX;
  std::vector<HistogramData::Counts> gsasDataY;
  std::vector<HistogramData::CountStandardDeviations> gsasDataE;

122
123
124
  std::vector<double> vecX, vecY, vecE;

  // progress
Dimitar Tasev's avatar
Dimitar Tasev committed
125
  std::unique_ptr<Progress> prog = nullptr;
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

  // Parameters for reading file
  char currentLine[256];
  std::string wsTitle;
  std::string slogTitle;
  std::string instrumentname = "Generic";
  char filetype = 'x';

  // Gather data
  std::ifstream input(filename.c_str(), std::ios_base::in);
  if (!input.is_open()) {
    // throw exception if file cannot be opened
    std::stringstream errss;
    errss << "Unable to open GSAS file " << filename;
    throw std::runtime_error(errss.str());
Zhou, Wenduo's avatar
Zhou, Wenduo committed
141
142
  }

143
144
145
146
147
148
149
  // First line: Title
  if (!input.eof()) {
    // Get workspace title (should be first line or 2nd line for SLOG)
    input.getline(currentLine, 256);
    wsTitle = currentLine;
  } else {
    throw std::runtime_error("File is empty");
Zhou, Wenduo's avatar
Zhou, Wenduo committed
150
  }
151

152
153
154
155
156
157
158
159
160
161
162
  // Loop all the lines
  bool isOutOfHead = false;
  bool slogtitleset = false;
  bool multiplybybinwidth = false;
  int nSpec = 0;
  bool calslogx0 = true;
  double bc4 = 0;
  double bc3 = 0;

  while (!input.eof() && input.getline(currentLine, 256)) {
    // Initialize progress after NSpec is imported
Dimitar Tasev's avatar
Dimitar Tasev committed
163
    if (nSpec != 0 && !prog) {
164
      prog = std::make_unique<Progress>(this, 0.0, 1.0, nSpec);
165
    }
166

167
168
169
170
    // Set flag to test SLOG
    if (!slogtitleset) {
      slogTitle = currentLine;
      slogtitleset = true;
Zhou, Wenduo's avatar
Zhou, Wenduo committed
171
    }
172

173
174
175
176
177
178
179
180
181
    if (currentLine[0] == '\n' || currentLine[0] == '#') {
      // Comment/information line
      std::string key1, key2;
      std::istringstream inputLine(currentLine, std::ios::in);
      inputLine.ignore(256, ' ');
      inputLine >> key1 >> key2;

      if (key2 == "Histograms") {
        // NSpec (Format: 'nspec HISTOGRAM')
182
        nSpec = std::stoi(key1);
183
        g_log.information() << "Histogram Line:  " << key1 << "  nSpec = " << nSpec << "\n";
184
185
186
187
188
189
190
191
192
193
      } else if (key1 == "Instrument:") {
        // Instrument (Format: 'Instrument XXXX')
        instrumentname = key2;
        g_log.information() << "Instrument    :  " << key2 << "\n";
      } else if (key1 == "with") {
        // Multiply by bin width: (Format: 'with multiplied')
        std::string s1;
        inputLine >> s1;
        if (s1 == "multiplied") {
          multiplybybinwidth = true;
194
          g_log.information() << "Y is multiplied by bin width\n";
195
        } else {
196
          g_log.warning() << "In line '" << currentLine << "', key word " << s1 << " is not allowed!\n";
Zhou, Wenduo's avatar
Zhou, Wenduo committed
197
        }
198
199
      } else if (key1 == "Primary") {
        // Primary flight path ...
200
        boost::smatch result;
201
202
203
        // Have to force a copy of the input or the stack gets corrupted
        // on MSVC when inputLine.str() falls out of scope which then
        // corrupts the value in result
204
        const std::string line = inputLine.str();
205
        if (boost::regex_search(line, result, L1_REG_EXP) && result.size() == 2) {
Peterson, Peter's avatar
Peterson, Peter committed
206
          primaryflightpath = std::stod(std::string(result[1]));
207
208
209

        } else {
          std::stringstream msg;
210
          msg << "Failed to parse primary flight path from line \"" << inputLine.str() << "\"";
211
212
213
214
215
216
          g_log.warning(msg.str());
        }

        std::stringstream msg;
        msg << "L1 = " << primaryflightpath;
        g_log.information(msg.str());
217
218
219
      } else if (key1 == "Total") {
        // Total flight path .... .... including total flying path, difc and
        // 2theta of 1 bank
220
221
222
223
224
225

        double totalpath(0.f);
        double tth(0.f);
        double difc(0.f);

        boost::smatch result;
226
        const std::string line = inputLine.str();
227
        if (boost::regex_search(line, result, DET_POS_REG_EXP) && result.size() == 4) {
Peterson, Peter's avatar
Peterson, Peter committed
228
229
230
          totalpath = std::stod(std::string(result[1]));
          tth = std::stod(std::string(result[2]));
          difc = std::stod(std::string(result[3]));
231
232
        } else {
          std::stringstream msg;
233
          msg << "Failed to parse position from line \"" << inputLine.str() << "\"";
234
235
          g_log.warning(msg.str());
        }
236

237
238
239
        totalflightpaths.emplace_back(totalpath);
        twothetas.emplace_back(tth);
        difcs.emplace_back(difc);
240

241
        std::stringstream msg;
242
        msg << "Bank " << difcs.size() - 1 << ": Total flight path = " << totalpath << "  2Theta = " << tth
243
244
            << "  DIFC = " << difc;
        g_log.information(msg.str());
245
246
247
248
249
250
251
      } // if keys....

    } // ENDIF for Line with #
    else if (currentLine[0] == 'B') {
      // Line start with Bank including file format, X0 information and etc.
      isOutOfHead = true;

252
      // If there is, Save the previous to array and initialize new MantiVec for
253
      // (X, Y, E)
254
      if (!vecX.empty()) {
255
256
257
        gsasDataX.emplace_back(std::move(vecX));
        gsasDataY.emplace_back(std::move(vecY));
        gsasDataE.emplace_back(std::move(vecE));
258
259
260
261
        vecX.clear();
        vecY.clear();
        vecE.clear();

262
        if (prog != nullptr)
263
264
          prog->report();
      }
Zhou, Wenduo's avatar
Zhou, Wenduo committed
265

266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
      // Parse the bank line in format
      // RALF: BANK <SpectraNo> <NBins> <NBins> RALF <BC1> <BC2> <BC1> <BC4>
      // SLOG: BANK <SpectraNo> <NBins> <NBins> SLOG <BC1> <BC2> <BC3> 0>
      // where,
      // BC1 = X[0] * 32
      // BC2 = X[1] * 32 - BC1
      // BC4 = ( X[1] - X[0] ) / X[0]

      int specno, nbin1, nbin2;
      std::istringstream inputLine(currentLine, std::ios::in);

      double bc1 = 0;
      double bc2 = 0;

      inputLine.ignore(256, 'K');
      std::string filetypestring;

      inputLine >> specno >> nbin1 >> nbin2 >> filetypestring;
284
      g_log.debug() << "Bank: " << specno << "  filetypestring = " << filetypestring << '\n';
285

286
      detectorIDs.emplace_back(specno);
287
288
289
290
291
292
293
294
295
296

      if (filetypestring[0] == 'S') {
        // SLOG
        filetype = 's';
        inputLine >> bc1 >> bc2 >> bc3 >> bc4;
      } else if (filetypestring[0] == 'R') {
        // RALF
        filetype = 'r';
        inputLine >> bc1 >> bc2 >> bc1 >> bc4;
      } else {
297
        g_log.error() << "Unsupported GSAS File Type: " << filetypestring << "\n";
298
299
        throw Exception::FileError("Not a GSAS file", filename);
      }
Zhou, Wenduo's avatar
Zhou, Wenduo committed
300

301
302
303
      // Determine x0
      if (filetype == 'r') {
        double x0 = bc1 / 32;
304
        g_log.debug() << "RALF: x0 = " << x0 << "  bc4 = " << bc4 << '\n';
305
        vecX.emplace_back(x0);
306
307
308
309
310
311
312
313
314
315
316
317
318
319
      } else {
        // Cannot calculate x0, turn on the flag
        calslogx0 = true;
      }
    } // Line with B
    else if (isOutOfHead) {
      // Parse data line
      double xValue;
      double yValue;
      double eValue;

      double xPrev;

      // * Get previous X value
320
      if (!vecX.empty()) {
321
322
323
        xPrev = vecX.back();
      } else if (filetype == 'r') {
        // Except if RALF
324
        throw Mantid::Kernel::Exception::NotImplementedError("LoadGSS: File was not in expected format.");
325
326
327
      } else {
        xPrev = -0.0;
      }
Zhou, Wenduo's avatar
Zhou, Wenduo committed
328

329
330
331
332
333
334
335
      // It is different for the definition of X, Y, Z in SLOG and RALF format
      if (filetype == 'r') {
        // RALF
        // LoadGSS produces overlapping columns for some datasets, due to
        // std::setw
        // For this reason we need to read the column values as string and then
        // convert to double
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
        {
          std::string str(currentLine, 15);
          std::istringstream istr(str);
          istr >> xValue;
        }
        {
          std::string str(currentLine + 15, 18);
          std::istringstream istr(str);
          istr >> yValue;
        }
        {
          std::string str(currentLine + 15 + 18, 18);
          std::istringstream istr(str);
          istr >> eValue;
        }

352
353
354
355
        xValue = (2 * xValue) - xPrev;

      } else if (filetype == 's') {
        // SLOG
356
        std::istringstream inputLine(currentLine, std::ios::in);
357
358
359
        inputLine >> xValue >> yValue >> eValue;
        if (calslogx0) {
          // calculation of x0 must use the x'[0]
360
          g_log.debug() << "x'_0 = " << xValue << "  bc3 = " << bc3 << '\n';
361
362

          double x0 = 2 * xValue / (bc3 + 2.0);
363
          vecX.emplace_back(x0);
364
          xPrev = x0;
365
          g_log.debug() << "SLOG: x0 = " << x0 << '\n';
366
          calslogx0 = false;
Zhou, Wenduo's avatar
Zhou, Wenduo committed
367
368
        }

369
370
371
372
        xValue = (2 * xValue) - xPrev;
      } else {
        g_log.error() << "Unsupported GSAS File Type: " << filetype << "\n";
        throw Exception::FileError("Not a GSAS file", filename);
Zhou, Wenduo's avatar
Zhou, Wenduo committed
373
374
      }

375
376
377
      if (multiplybybinwidth) {
        yValue = yValue / (xValue - xPrev);
        eValue = eValue / (xValue - xPrev);
378
      }
Zhou, Wenduo's avatar
Zhou, Wenduo committed
379

380
      // store read in data (x, y, e) to vector
381
382
383
      vecX.emplace_back(std::move(xValue));
      vecY.emplace_back(std::move(yValue));
      vecE.emplace_back(std::move(eValue));
384
385
    } // Date Line
    else {
386
      g_log.warning() << "Line not defined: " << currentLine << '\n';
387
    }
388
389
390
  } // ENDWHILE of reading all lines

  // Get the sizes before using std::move
391
392
393
  auto nHist(static_cast<int>(gsasDataX.size()));
  auto xWidth(static_cast<int>(vecX.size()));
  auto yWidth(static_cast<int>(vecY.size()));
394

395
  // Push the vectors (X, Y, E) of the last bank to gsasData
396
  if (!vecX.empty()) { // Put final spectra into data
397
398
399
400
    gsasDataX.emplace_back(std::move(vecX));
    gsasDataY.emplace_back(std::move(vecY));
    gsasDataE.emplace_back(std::move(vecE));
    ++nHist;
401
402
403
404
405
406
407
408
409
  }
  input.close();

  //********************************************************************************************
  // Construct the workspace for GSS data
  //********************************************************************************************

  // Create workspace & GSS Files data is always in TOF

410
411
  MatrixWorkspace_sptr outputWorkspace = std::dynamic_pointer_cast<MatrixWorkspace>(
      WorkspaceFactory::Instance().create("Workspace2D", nHist, xWidth, yWidth));
412
413
414
415
416
417
418
419
  outputWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");

  // set workspace title
  if (filetype == 'r')
    outputWorkspace->setTitle(wsTitle);
  else
    outputWorkspace->setTitle(slogTitle);

420
  // put data from constructed histograms into outputWorkspace
421
422
423
  if (detectorIDs.size() != static_cast<size_t>(nHist)) {
    // File error is found
    std::ostringstream mess("");
424
    mess << "Number of spectra (" << detectorIDs.size() << ") is not equal to number of histograms (" << nHist << ").";
425
    throw std::runtime_error(mess.str());
Zhou, Wenduo's avatar
Zhou, Wenduo committed
426
427
  }

428
429
  for (int i = 0; i < nHist; ++i) {
    // Move data across
430
431
    outputWorkspace->setHistogram(i, BinEdges(std::move(gsasDataX[i])), Counts(std::move(gsasDataY[i])),
                                  CountStandardDeviations(std::move(gsasDataE[i])));
Zhou, Wenduo's avatar
Zhou, Wenduo committed
432

433
434
    // Reset spectrum number if
    if (useBankAsSpectrum) {
435
      auto specno = static_cast<specnum_t>(detectorIDs[i]);
436
      outputWorkspace->getSpectrum(i).setSpectrumNo(specno);
Zhou, Wenduo's avatar
Zhou, Wenduo committed
437
438
    }
  }
439

440
  // build instrument geometry
Danny Hindson's avatar
Danny Hindson committed
441
442
  createInstrumentGeometry(outputWorkspace, instrumentname, primaryflightpath,
                           detectorIDs, totalflightpaths, twothetas, difcs);
443
444
445
446
447
448

  return outputWorkspace;
}

//----------------------------------------------------------------------------------------------
/** Convert a string containing number and unit to double
LamarMoore's avatar
LamarMoore committed
449
 */
450
double LoadGSS::convertToDouble(std::string inputstring) {
451
  std::string temps;
452
  auto isize = static_cast<int>(inputstring.size());
453
454
  for (int i = 0; i < isize; i++) {
    char thechar = inputstring[i];
455
    if ((thechar <= 'Z' && thechar >= 'A') || (thechar <= 'z' && thechar >= 'a')) {
456
457
458
      break;
    } else {
      temps += thechar;
459
    }
460
  }
461

Peterson, Peter's avatar
Peterson, Peter committed
462
  double rd = std::stod(temps);
463
464
465
466
467
468

  return rd;
}

//----------------------------------------------------------------------------------------------
/** Create the instrument geometry with Instrument
LamarMoore's avatar
LamarMoore committed
469
 */
Danny Hindson's avatar
Danny Hindson committed
470
471
472
473
474
void LoadGSS::createInstrumentGeometry(
    const 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) {
475
  // Check Input
476
  if (detectorids.size() != totalflightpaths.size() || totalflightpaths.size() != twothetas.size()) {
477
478
479
480
481
482
483
484
485
    g_log.warning("Cannot create geometry, because the numbers of L2 and Polar "
                  "are not equal.");
    return;
  }

  // Debug output
  std::stringstream dbss;
  dbss << "L1 = " << primaryflightpath << "\n";
  for (size_t i = 0; i < detectorids.size(); i++) {
486
487
    dbss << "Detector " << detectorids[i] << "  L1+L2 = " << totalflightpaths[i] << "  2Theta = " << twothetas[i]
         << "\n";
488
489
490
491
  }
  g_log.debug(dbss.str());

  // Create a new instrument and set its name
492
  Geometry::Instrument_sptr instrument(new Geometry::Instrument(instrumentname));
493
494

  // Add dummy source and samplepos to instrument
495
  Geometry::Component *samplepos = new Geometry::Component("Sample", instrument.get());
496
497
498
499
  instrument->add(samplepos);
  instrument->markAsSamplePos(samplepos);
  samplepos->setPos(0.0, 0.0, 0.0);

500
  Geometry::ObjComponent *source = new Geometry::ObjComponent("Source", instrument.get());
501
502
503
504
505
506
507
508
509
  instrument->add(source);
  instrument->markAsSource(source);

  double l1 = primaryflightpath;
  source->setPos(0.0, 0.0, -1.0 * l1);

  // Add detectors
  // The L2 and 2-theta values from Raw file assumed to be relative to sample
  // position
510
  const auto numDetector = static_cast<int>(detectorids.size()); // number of detectors
511
512
513
514
515
516
517
518
  // std::vector<int> detID = detectorids;    // detector IDs
  // std::vector<double> angle = twothetas;  // angle between indicent beam and
  // direction from sample to detector (two-theta)

  // Assumption: detector IDs are in the same order of workspace index
  for (int i = 0; i < numDetector; ++i) {
    // a) Create a new detector. Instrument will take ownership of pointer so no
    // need to delete.
519
    Geometry::Detector *detector = new Geometry::Detector("det", detectorids[i], samplepos);
520
521
522
523
524
525
526
527
528
    Kernel::V3D pos;

    // r is L2
    double r = totalflightpaths[i] - l1;
    pos.spherical(r, twothetas[i], 0.0);

    detector->setPos(pos);

    // add copy to instrument, spectrum and mark it
529
530
531
532
533
    auto &spec = workspace->getSpectrum(i);
    spec.clearDetectorIDs();
    spec.addDetectorID(detectorids[i]);
    instrument->add(detector);
    instrument->markAsDetector(detector);
534

535
  } // ENDFOR (i: spectrum)
536
  workspace->setInstrument(instrument);
Zhou, Wenduo's avatar
Zhou, Wenduo committed
537

Danny Hindson's avatar
Danny Hindson committed
538
539
540
541
542
543
  auto &paramMap = workspace->instrumentParameters();
  for (size_t i = 0; i < workspace->getNumberHistograms(); i++) {
    auto detector = workspace->getDetector(i);
    paramMap.addDouble(detector->getComponentID(), "DIFC", difcs[i]);
  }
}
LamarMoore's avatar
LamarMoore committed
544
} // namespace DataHandling
Danny Hindson's avatar
Danny Hindson committed
545

LamarMoore's avatar
LamarMoore committed
546
} // namespace Mantid