LoadSpiceXML2DDet.cpp 32.6 KB
Newer Older
1
#include "MantidDataHandling/LoadSpiceXML2DDet.h"
2
#include "MantidAPI/Axis.h"
3
#include "MantidAPI/FileProperty.h"
4
#include "MantidAPI/ITableWorkspace.h"
Hahn, Steven's avatar
Hahn, Steven committed
5
#include "MantidAPI/MatrixWorkspace.h"
6
#include "MantidAPI/Run.h"
7
#include "MantidAPI/WorkspaceFactory.h"
Hahn, Steven's avatar
Hahn, Steven committed
8
#include "MantidAPI/WorkspaceProperty.h"
9
#include "MantidKernel/ArrayProperty.h"
Hahn, Steven's avatar
Hahn, Steven committed
10
#include "MantidKernel/OptionalBool.h"
11
#include "MantidKernel/TimeSeriesProperty.h"
12

13
#include <boost/algorithm/string.hpp>
14

15
16
17
18
19
20
21
22
23
24
#include <Poco/DOM/AutoPtr.h>
#include <Poco/DOM/Document.h>
#include <Poco/DOM/DOMParser.h>
#include <Poco/DOM/NamedNodeMap.h>
#include <Poco/DOM/Node.h>
#include <Poco/DOM/NodeFilter.h>
#include <Poco/DOM/NodeIterator.h>
#include <Poco/DOM/NodeList.h>
#include <Poco/SAX/InputSource.h>

25
26
27
#include <algorithm>
#include <fstream>

28
29
30
31
32
33
namespace Mantid {
namespace DataHandling {

using namespace Mantid::API;
using namespace Mantid::Kernel;

34
35
DECLARE_ALGORITHM(LoadSpiceXML2DDet)

36
37
38
39
const char STRING = 's';
const char FLOAT32 = 'f';
const char INT32 = 'i';

40
41
42
43
/** Constructor for SpiceXMLNode
 * @brief SpiceXMLNode::SpiceXMLNode
 * @param nodename
 */
44
SpiceXMLNode::SpiceXMLNode(const std::string &nodename)
45
    : m_name{nodename}, m_typechar('s') {}
46

47
48
49
50
/** Set node value in string format
 * @brief SpiceXMLNode::setValue
 * @param strvalue
 */
51
52
void SpiceXMLNode::setValue(const std::string &strvalue) { m_value = strvalue; }

53
54
55
56
57
58
59
60
61
/** Set XML node parameters
 * @brief SpiceXMLNode::setValues
 * @param nodetype
 * @param nodeunit
 * @param nodedescription
 */
void SpiceXMLNode::setParameters(const std::string &nodetype,
                                 const std::string &nodeunit,
                                 const std::string &nodedescription) {
62
  // data type
63
  if (nodetype == "FLOAT32") {
64
65
    m_typefullname = nodetype;
    m_typechar = FLOAT32;
66
  } else if (nodetype == "INT32") {
67
68
69
70
71
    m_typefullname = nodetype;
    m_typechar = INT32;
  }

  // unit
72
  if (!nodeunit.empty()) {
73
74
75
76
    m_unit = nodeunit;
  }

  // description
77
  if (!nodedescription.empty())
78
79
80
    m_description = nodedescription;
}

81
82
/** Check whether XML has unit set
 */
83
bool SpiceXMLNode::hasUnit() const { return (!m_unit.empty()); }
84

85
86
87
88
/** Check whether XML node has value set
 * @brief SpiceXMLNode::hasValue
 * @return
 */
89
bool SpiceXMLNode::hasValue() const { return (!m_value.empty()); }
90

91
92
93
94
/** Is this node of string type?
 * @brief SpiceXMLNode::isString
 * @return
 */
95
bool SpiceXMLNode::isString() const { return (m_typechar == STRING); }
96

97
98
99
100
/** Is this node of integer type?
 * @brief SpiceXMLNode::isInteger
 * @return
 */
101
bool SpiceXMLNode::isInteger() const { return (m_typechar == INT32); }
102

103
104
105
106
/** Is this node of double type?
 * @brief SpiceXMLNode::isDouble
 * @return
 */
107
bool SpiceXMLNode::isDouble() const { return (m_typechar == FLOAT32); }
108

109
110
111
112
/** Get name of XML node
 * @brief SpiceXMLNode::getName
 * @return
 */
113
const std::string SpiceXMLNode::getName() const { return m_name; }
114
115
116
117
118

/** Get unit of XML node
 * @brief SpiceXMLNode::getUnit
 * @return
 */
119
const std::string SpiceXMLNode::getUnit() const { return m_unit; }
120
121
122
123
124

/** Get node's description
 * @brief SpiceXMLNode::getDescription
 * @return
 */
125
const std::string SpiceXMLNode::getDescription() const { return m_description; }
126
127
128
129
130

/** Get node's value in string
 * @brief SpiceXMLNode::getValue
 * @return
 */
131
132
const std::string SpiceXMLNode::getValue() const { return m_value; }

133
134
/** Constructor
 */
135
136
137
138
LoadSpiceXML2DDet::LoadSpiceXML2DDet()
    : m_detXMLFileName(), m_detXMLNodeName(), m_numPixelX(0), m_numPixelY(0),
      m_loadInstrument(false), m_detSampleDistanceShift(0.0),
      m_hasScanTable(false), m_ptNumber4Log(0), m_idfFileName() {}
139
140
141

/** Destructor
 */
142
LoadSpiceXML2DDet::~LoadSpiceXML2DDet() = default;
143

144
145
146
147
148
149
const std::string LoadSpiceXML2DDet::name() const {
  return "LoadSpiceXML2DDet";
}

int LoadSpiceXML2DDet::version() const { return 1; }

Nick Draper's avatar
Nick Draper committed
150
151
152
const std::string LoadSpiceXML2DDet::category() const {
  return "DataHandling\\XML";
}
153
154
155
156
157
158

const std::string LoadSpiceXML2DDet::summary() const {
  return "Load 2-dimensional detector data file in XML format from SPICE. ";
}

/** Declare properties
159
160
161
162
 * @brief LoadSpiceXML2DDet::init
 */
void LoadSpiceXML2DDet::init() {
  declareProperty(
163
164
      make_unique<FileProperty>("Filename", "", FileProperty::FileAction::Load,
                                ".xml"),
165
166
      "XML file name for one scan including 2D detectors counts from SPICE");

167
  declareProperty(
168
169
      make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "",
                                                      Direction::Output),
170
171
172
      "Name of output matrix workspace. "
      "Output workspace will be an X by Y Workspace2D if instrument "
      "is not loaded. ");
173

Zhou, Wenduo's avatar
Zhou, Wenduo committed
174
175
176
177
  declareProperty(
      "DetectorLogName", "Detector",
      "Log name (i.e., XML node name) for detector counts in XML file."
      "By default, the name is 'Detector'");
178
179

  declareProperty(
180
      make_unique<ArrayProperty<size_t>>("DetectorGeometry"),
181
182
      "A size-2 unsigned integer array [X, Y] for detector geometry. "
      "Such that the detector contains X x Y pixels.");
183

184
185
186
187
  declareProperty(
      "LoadInstrument", true,
      "Flag to load instrument to output workspace. "
      "HFIR's HB3A will be loaded if InstrumentFileName is not specified.");
188
189

  declareProperty(
190
191
      make_unique<FileProperty>("InstrumentFilename", "",
                                FileProperty::OptionalLoad, ".xml"),
192
193
194
195
196
      "The filename (including its full or relative path) of an instrument "
      "definition file. The file extension must either be .xml or .XML when "
      "specifying an instrument definition file. Note Filename or "
      "InstrumentName must be specified but not both.");

197
  declareProperty(
198
      make_unique<WorkspaceProperty<ITableWorkspace>>(
199
          "SpiceTableWorkspace", "", Direction::Input, PropertyMode::Optional),
200
      "Name of TableWorkspace loaded from SPICE scan file by LoadSpiceAscii.");
201

202
203
  declareProperty("PtNumber", 0,
                  "Pt. value for the row to get sample log from. ");
204

205
  declareProperty("UserSpecifiedWaveLength", EMPTY_DBL(),
206
207
                  "User can specify the wave length of the instrument if it is "
                  "drifted from the designed value."
208
                  "It happens often.");
209

Zhou, Wenduo's avatar
Zhou, Wenduo committed
210
211
212
213
  declareProperty(
      "ShiftedDetectorDistance", 0.,
      "Amount of shift of the distance between source and detector centre."
      "It is used to apply instrument calibration.");
214

215
216
217
  declareProperty("DetectorCenterXShift", 0.0, "The amount of shift of "
                                               "detector center along X "
                                               "direction in the unit meter.");
218

219
220
221
  declareProperty("DetectorCenterYShift", 0.0, "The amount of shift of "
                                               "detector center along Y "
                                               "direction in the unit meter.");
222
223
}

224
225
/** Process inputs arguments
 * @brief processInputs
226
 */
Zhou, Wenduo's avatar
Zhou, Wenduo committed
227
void LoadSpiceXML2DDet::processInputs() {
228
229
  m_detXMLFileName = getPropertyValue("Filename");
  m_detXMLNodeName = getPropertyValue("DetectorLogName");
230
  std::vector<size_t> vec_pixelgeom = getProperty("DetectorGeometry");
231
  if (vec_pixelgeom.size() == 2) {
232
233
    m_numPixelX = vec_pixelgeom[0];
    m_numPixelY = vec_pixelgeom[1];
234
  } else if (vec_pixelgeom.empty()) {
235
236
    m_numPixelX = 0;
    m_numPixelY = 0;
237
238
239
240
  } else {
    throw std::runtime_error("Input pixels geometry is not correct in format. "
                             "It either has 2 integers or left empty to get "
                             "determined automatically.");
241
  }
242
243
  g_log.debug() << "User input poixels numbers: " << m_numPixelX << ", "
                << m_numPixelY << "\n";
244

245
  m_loadInstrument = getProperty("LoadInstrument");
246

247
248
  m_idfFileName = getPropertyValue("InstrumentFilename");
  m_detSampleDistanceShift = getProperty("ShiftedDetectorDistance");
249

250
  // Retreive sample environment data from SPICE scan table workspace
251
  std::string spicetablewsname = getPropertyValue("SpiceTableWorkspace");
252
  if (!spicetablewsname.empty())
253
254
255
    m_hasScanTable = true;
  else
    m_hasScanTable = false;
256

257
  m_ptNumber4Log = getProperty("PtNumber");
258
259

  m_userSpecifiedWaveLength = getProperty("UserSpecifiedWaveLength");
260
261
262

  m_detXShift = getProperty("DetectorCenterXShift");
  m_detYShift = getProperty("DetectorCenterYShift");
263
264
265
266
}

/** Set up sample logs especially 2theta and diffr for loading instrument
 * @brief LoadSpiceXML2DDet::setupSampleLogs
267
 * @param outws :: workspace to have sample logs to set up
268
269
 * @return
 */
Zhou, Wenduo's avatar
Zhou, Wenduo committed
270
bool LoadSpiceXML2DDet::setupSampleLogs(API::MatrixWorkspace_sptr outws) {
271
  // With given spice scan table, 2-theta is read from there.
Zhou, Wenduo's avatar
Zhou, Wenduo committed
272
  if (m_hasScanTable) {
273
274
    ITableWorkspace_sptr spicetablews = getProperty("SpiceTableWorkspace");
    setupSampleLogFromSpiceTable(outws, spicetablews, m_ptNumber4Log);
275
276
  }

277
278
  Kernel::DateAndTime anytime(1000);

279
280
  // Process 2theta
  bool return_true = true;
Zhou, Wenduo's avatar
Zhou, Wenduo committed
281
282
  if (!outws->run().hasProperty("2theta") &&
      outws->run().hasProperty("_2theta")) {
283
    // Set up 2theta if it is not set up yet
Peterson, Peter's avatar
Peterson, Peter committed
284
    double logvalue = std::stod(outws->run().getProperty("_2theta")->value());
285
286
287
288
    TimeSeriesProperty<double> *newlogproperty =
        new TimeSeriesProperty<double>("2theta");
    newlogproperty->addValue(anytime, logvalue);
    outws->mutableRun().addProperty(newlogproperty);
Zhou, Wenduo's avatar
Zhou, Wenduo committed
289
290
291
292
    g_log.information() << "Set 2theta from _2theta (as XML node) with value "
                        << logvalue << "\n";
  } else if (!outws->run().hasProperty("2theta") &&
             !outws->run().hasProperty("_2theta")) {
293
294
295
    // Neither 2theta nor _2theta
    g_log.warning("No 2theta is set up for loading instrument.");
    return_true = false;
296
  }
297

298
  // set up the caibrated detector center to beam
299
300
301
  TimeSeriesProperty<double> *det_dx = new TimeSeriesProperty<double>("deltax");
  det_dx->addValue(anytime, m_detXShift);
  outws->mutableRun().addProperty(det_dx);
302

303
304
305
  TimeSeriesProperty<double> *det_dy = new TimeSeriesProperty<double>("deltay");
  det_dy->addValue(anytime, m_detYShift);
  outws->mutableRun().addProperty(det_dy);
306

307
  // set up Sample-detetor distance calibration
308
  double sampledetdistance = m_detSampleDistanceShift;
309
310
311
312
313
  TimeSeriesProperty<double> *distproperty =
      new TimeSeriesProperty<double>("diffr");
  distproperty->addValue(anytime, sampledetdistance);
  outws->mutableRun().addProperty(distproperty);

314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
  return return_true;
}

/** Main execution
 * @brief LoadSpiceXML2DDet::exec
 */
void LoadSpiceXML2DDet::exec() {
  // Load input
  processInputs();

  // Parse detector XML file
  std::vector<SpiceXMLNode> vec_xmlnode = parseSpiceXML(m_detXMLFileName);

  // Create output workspace
  MatrixWorkspace_sptr outws;
329
330
331
  if (m_numPixelX * m_numPixelY > 0)
    outws = createMatrixWorkspace(vec_xmlnode, m_numPixelX, m_numPixelY,
                                  m_detXMLNodeName, m_loadInstrument);
332
  else
333
334
    outws = createMatrixWorkspaceVersion2(vec_xmlnode, m_detXMLNodeName,
                                          m_loadInstrument);
335
336
337
338
339
340

  // Set up log for loading instrument
  bool can_set_instrument = setupSampleLogs(outws);

  if (m_loadInstrument && can_set_instrument) {
    loadInstrument(outws, m_idfFileName);
341
342
    // set wave length to user specified wave length
    double wavelength = m_userSpecifiedWaveLength;
343
344
    // if user does not specify wave length then try to get wave length from log
    // sample _m1 (or m1 as well in future)
345
    bool has_wavelength = !(wavelength == EMPTY_DBL());
346
347
348
    if (!has_wavelength)
      has_wavelength = getHB3AWavelength(outws, wavelength);

349
350
    if (has_wavelength) {
      setXtoLabQ(outws, wavelength);
351
352
    }
  }
353
354
355
356

  setProperty("OutputWorkspace", outws);
}

357
358
/** Parse SPICE XML file for one Pt./measurement
 * @brief LoadSpiceXML2DDet::parseSpiceXML
359
 * @param xmlfilename :: name of the XML file to parse
360
 * @return vector of SpiceXMLNode containing information in XML file
361
 */
Zhou, Wenduo's avatar
Zhou, Wenduo committed
362
363
std::vector<SpiceXMLNode>
LoadSpiceXML2DDet::parseSpiceXML(const std::string &xmlfilename) {
364
365
366
  // Declare output
  std::vector<SpiceXMLNode> vecspicenode;

367
  // Open file
368
369
  std::ifstream ifs;
  ifs.open(xmlfilename.c_str());
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
  if (!ifs.is_open()) {
    std::stringstream ess;
    ess << "File " << xmlfilename << " cannot be opened.";
    throw std::runtime_error(ess.str());
  }

  // Parse
  Poco::XML::InputSource src(ifs);

  Poco::XML::DOMParser parser;
  Poco::AutoPtr<Poco::XML::Document> pDoc = parser.parse(&src);

  // Go though XML
  Poco::XML::NodeIterator nodeIter(pDoc, Poco::XML::NodeFilter::SHOW_ELEMENT);
  Poco::XML::Node *pNode = nodeIter.nextNode();
  while (pNode) {
    const Poco::XML::XMLString nodename = pNode->nodeName();
387
388
389
390
391
    // const Poco::XML::XMLString nodevalue = pNode->nodeValue();

    // get number of children
    size_t numchildren = pNode->childNodes()->length();
    if (numchildren > 1) {
392
393
394
      g_log.debug() << "Parent node " << nodename << " has " << numchildren
                    << " children."
                    << "\n";
395
      if (nodename == "SPICErack") {
396
        // SPICErack is the main parent node.  start_time and end_time are there
Zhou, Wenduo's avatar
Zhou, Wenduo committed
397
398
        unsigned long numattr = pNode->attributes()->length();
        for (unsigned long j = 0; j < numattr; ++j) {
399
400
401
402
403
404
405
406
407
408
          std::string attname = pNode->attributes()->item(j)->nodeName();
          std::string attvalue = pNode->attributes()->item(j)->innerText();
          SpiceXMLNode xmlnode(attname);
          xmlnode.setValue(attvalue);
          vecspicenode.push_back(xmlnode);
          g_log.debug() << "SPICErack attribute " << j << " Name = " << attname
                        << ", Value = " << attvalue << "\n";
        }
      }

409
    } else if (numchildren == 1) {
410
      std::string innertext = pNode->innerText();
411
412
413
      unsigned long numattr = pNode->attributes()->length();
      g_log.debug() << "  Child node " << nodename << "'s attributes: "
                    << "\n";
414
415

      SpiceXMLNode xmlnode(nodename);
416
417
418
      std::string nodetype;
      std::string nodeunit;
      std::string nodedescription;
419

420
      for (unsigned long j = 0; j < numattr; ++j) {
421
422
        std::string atttext = pNode->attributes()->item(j)->innerText();
        std::string attname = pNode->attributes()->item(j)->nodeName();
423
424
        g_log.debug() << "     attribute " << j << " name = " << attname << ", "
                      << "value = " << atttext << "\n";
425
        if (attname == "type") {
426
427
          // type
          nodetype = atttext;
428
        } else if (attname == "unit") {
429
430
          // unit
          nodeunit = atttext;
431
        } else if (attname == "description") {
432
433
434
          // description
          nodedescription = atttext;
        }
435
      }
436
      xmlnode.setParameters(nodetype, nodeunit, nodedescription);
437
438
      xmlnode.setValue(innertext);

439
      vecspicenode.push_back(xmlnode);
440
    } else {
441
      // An unexpected case but no guarantee for not happening
442
443
444
      g_log.error("Funny... No child node.");
    }

445
    // Move to next node
446
    pNode = nodeIter.nextNode();
447
  } // ENDWHILE
448
449
450
451

  // Close file
  ifs.close();

452
  return vecspicenode;
453
454
}

455
456
457
458
459
460
/** Create MatrixWorkspace from Spice XML file
 * @brief LoadSpiceXML2DDet::createMatrixWorkspace
 * @param vecxmlnode :: vector of SpiceXMLNode obtained from XML file
 * @param numpixelx :: number of pixel in x-direction
 * @param numpixely :: number of pixel in y-direction
 * @param detnodename :: the XML node's name for detector counts.
461
 * @param loadinstrument :: flag to load instrument to output workspace or not.
462
463
 * @return
 */
464
MatrixWorkspace_sptr LoadSpiceXML2DDet::createMatrixWorkspace(
465
    const std::vector<SpiceXMLNode> &vecxmlnode, const size_t &numpixelx,
466
467
    const size_t &numpixely, const std::string &detnodename,
    const bool &loadinstrument) {
468

469
470
  // TODO FIXME - If version 2 works, then this version will be discarded

471
  // Create matrix workspace
472
473
474
475
476
  MatrixWorkspace_sptr outws;

  if (loadinstrument) {
    size_t numspec = numpixelx * numpixely;
    outws = boost::dynamic_pointer_cast<MatrixWorkspace>(
477
        WorkspaceFactory::Instance().create("Workspace2D", numspec, 2, 1));
478
479
480
481
482
  } else {
    outws = boost::dynamic_pointer_cast<MatrixWorkspace>(
        WorkspaceFactory::Instance().create("Workspace2D", numpixely, numpixelx,
                                            numpixelx));
  }
483

484
485
486
  // Go through all XML nodes to process
  size_t numxmlnodes = vecxmlnode.size();
  bool parsedDet = false;
487
  double max_counts = 0.;
488
489
490
  for (size_t n = 0; n < numxmlnodes; ++n) {
    // Process node for detector's count
    const SpiceXMLNode &xmlnode = vecxmlnode[n];
491
    if (xmlnode.getName() == detnodename) {
492
493
      // Get node value string (256x256 as a whole)
      const std::string detvaluestr = xmlnode.getValue();
494

495
496
497
      // Split
      std::vector<std::string> vecLines;
      boost::split(vecLines, detvaluestr, boost::algorithm::is_any_of("\n"));
498
499
      g_log.debug() << "There are " << vecLines.size() << " lines"
                    << "\n";
500

501
      // XML file records data in the order of column-major
502
      size_t i_col = 0;
503
504
      for (size_t i = 0; i < vecLines.size(); ++i) {
        std::string &line = vecLines[i];
505

506
        // Skip empty line
507
        if (line.empty()) {
508
          g_log.debug() << "\tFound empty Line at " << i << "\n";
509
510
511
512
          continue;
        }

        // Check whether it exceeds boundary
513
        if (i_col == numpixelx) {
514
          std::stringstream errss;
515
          errss << "Number of non-empty rows (" << i_col + 1
516
                << ") in detector data "
517
                << "exceeds user defined geometry size " << numpixelx << ".";
518
          throw std::runtime_error(errss.str());
519
        }
520

521
522
523
524
        // Split line
        std::vector<std::string> veccounts;
        boost::split(veccounts, line, boost::algorithm::is_any_of(" \t"));

Zhou, Wenduo's avatar
Zhou, Wenduo committed
525
526
        // check number of counts per column should not exceeds number of pixels
        // in Y direction
527
        if (veccounts.size() != numpixely) {
528
          std::stringstream errss;
529
530
          errss << "[Version 1] Row " << i_col << " contains "
                << veccounts.size() << " items other than " << numpixely
531
532
533
534
                << " counts specified by user.";
          throw std::runtime_error(errss.str());
        }

535
        // scan per column
536
        for (size_t j_row = 0; j_row < veccounts.size(); ++j_row) {
Peterson, Peter's avatar
Peterson, Peter committed
537
          double counts = std::stod(veccounts[j_row]);
538
          size_t rowIndex, columnIndex;
539
540

          if (loadinstrument) {
541
            // the detector ID and ws index are set up in column-major too!
542
            rowIndex = i_col * numpixelx + j_row;
543
            columnIndex = 0;
544
          } else {
545
546
            rowIndex = j_row;
            columnIndex = i_col;
547
          }
548

549
550
          outws->mutableX(rowIndex)[columnIndex] =
              static_cast<double>(columnIndex);
551
552
553
          outws->mutableY(rowIndex)[columnIndex] = counts;

          if (counts > 0)
554
            outws->mutableE(rowIndex)[columnIndex] = sqrt(counts);
555
          else
556
            outws->mutableE(rowIndex)[columnIndex] = 1.0;
557

558
559
560
561
562
          // record max count
          if (counts > max_counts) {
            max_counts = counts;
          }
        }
563

564
        // Update column index (i.e., column number)
565
        i_col += 1;
566
      } // END-FOR (i-vec line)
567

568
569
570
571
572
573
574
575
      // Set flag
      parsedDet = true;
    } else {
      // Parse to log: because there is no start time.  so all logs are single
      // value type
      const std::string nodename = xmlnode.getName();
      const std::string nodevalue = xmlnode.getValue();
      if (xmlnode.isDouble()) {
Peterson, Peter's avatar
Peterson, Peter committed
576
        double dvalue = std::stod(nodevalue);
577
578
        outws->mutableRun().addProperty(
            new PropertyWithValue<double>(nodename, dvalue));
579
580
        g_log.debug() << "Log name / xml node : " << xmlnode.getName()
                      << " (double) value = " << dvalue << "\n";
581
      } else if (xmlnode.isInteger()) {
582
        int ivalue = std::stoi(nodevalue);
583
584
        outws->mutableRun().addProperty(
            new PropertyWithValue<int>(nodename, ivalue));
585
586
        g_log.debug() << "Log name / xml node : " << xmlnode.getName()
                      << " (int) value = " << ivalue << "\n";
587
      } else {
Zhou, Wenduo's avatar
Zhou, Wenduo committed
588
        std::string str_value(nodevalue);
589
        if (nodename == "start_time") {
590
          // replace 2015-01-17 13:36:45 by  2015-01-17T13:36:45
Zhou, Wenduo's avatar
Zhou, Wenduo committed
591
592
          str_value = nodevalue;
          str_value.replace(10, 1, "T");
593
594
          g_log.debug() << "Replace start_time " << nodevalue
                        << " by Mantid time format " << str_value << "\n";
Zhou, Wenduo's avatar
Zhou, Wenduo committed
595
        }
596
        outws->mutableRun().addProperty(
597
            new PropertyWithValue<std::string>(nodename, str_value));
598
      }
599
    }
600
601
602
603
  }

  // Raise exception if no detector node is found
  if (!parsedDet) {
604
605
606
607
608
    std::stringstream errss;
    errss << "Unable to find an XML node of name " << detnodename
          << ". Unable to load 2D detector XML file.";
    throw std::runtime_error(errss.str());
  }
609

610
611
  g_log.notice() << "Maximum detector count on it is " << max_counts << "\n";

612
613
614
  return outws;
}

615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
/** create the output matrix workspace without knowledge of detector geometry
 *
 */
MatrixWorkspace_sptr LoadSpiceXML2DDet::createMatrixWorkspaceVersion2(
    const std::vector<SpiceXMLNode> &vecxmlnode, const std::string &detnodename,
    const bool &loadinstrument) {

  // Create matrix workspace
  MatrixWorkspace_sptr outws;

  // Go through all XML nodes to process
  size_t numxmlnodes = vecxmlnode.size();
  bool parsedDet = false;
  double max_counts = 0.;

  // define log value map
  std::map<std::string, std::string> str_log_map;
  std::map<std::string, double> dbl_log_map;
  std::map<std::string, int> int_log_map;

  for (size_t n = 0; n < numxmlnodes; ++n) {
    // Process node for detector's count
    const SpiceXMLNode &xmlnode = vecxmlnode[n];
638
    if (xmlnode.getName() == detnodename) {
639
640
641
642
643
644
645
646
647
648
649
650
651
      // Get node value string (256x256 as a whole)
      const std::string detvaluestr = xmlnode.getValue();

      outws = this->parseDetectorNode(detvaluestr, loadinstrument, max_counts);

      // Set flag
      parsedDet = true;
    } else {
      // Parse to log: because there is no start time.  so all logs are single
      // value type
      const std::string nodename = xmlnode.getName();
      const std::string nodevalue = xmlnode.getValue();
      if (xmlnode.isDouble()) {
Peterson, Peter's avatar
Peterson, Peter committed
652
        double dvalue = std::stod(nodevalue);
653
        dbl_log_map.emplace(nodename, dvalue);
654
      } else if (xmlnode.isInteger()) {
655
        int ivalue = std::stoi(nodevalue);
656
        int_log_map.emplace(nodename, ivalue);
657
      } else {
658
        if (nodename == "start_time") {
659
660
661
          // replace 2015-01-17 13:36:45 by  2015-01-17T13:36:45
          std::string str_value(nodevalue);
          str_value.replace(10, 1, "T");
662
663
          g_log.debug() << "Replace start_time " << nodevalue
                        << " by Mantid time format " << str_value << "\n";
664
          str_log_map.emplace(nodename, str_value);
665
        } else
666
667
          str_log_map.emplace(nodename, nodevalue);
      } // END-IF-ELSE (node value type)
668
669
    }   // END-IF-ELSE (detector-node or log node)
  }     // END-FOR (xml nodes)
670

671
  // Add the property to output workspace
672
673
674
675
  for (std::map<std::string, std::string>::iterator miter = str_log_map.begin();
       miter != str_log_map.end(); ++miter) {
    outws->mutableRun().addProperty(
        new PropertyWithValue<std::string>(miter->first, miter->second));
676
  }
677
678
679
680
  for (std::map<std::string, int>::iterator miter = int_log_map.begin();
       miter != int_log_map.end(); ++miter) {
    outws->mutableRun().addProperty(
        new PropertyWithValue<int>(miter->first, miter->second));
681
  }
682
683
684
685
  for (std::map<std::string, double>::iterator miter = dbl_log_map.begin();
       miter != dbl_log_map.end(); ++miter) {
    outws->mutableRun().addProperty(
        new PropertyWithValue<double>(miter->first, miter->second));
686
687
  }

688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
  // Raise exception if no detector node is found
  if (!parsedDet) {
    std::stringstream errss;
    errss << "Unable to find an XML node of name " << detnodename
          << ". Unable to load 2D detector XML file.";
    throw std::runtime_error(errss.str());
  }

  g_log.notice() << "Maximum detector count on it is " << max_counts << "\n";

  return outws;
}

/**
 */
703
704
705
706
707
708
709
710
711
712
713
API::MatrixWorkspace_sptr
LoadSpiceXML2DDet::parseDetectorNode(const std::string &detvaluestr,
                                     bool loadinstrument, double &max_counts) {
  // Split to lines
  std::vector<std::string> vecLines;
  boost::split(vecLines, detvaluestr, boost::algorithm::is_any_of("\n"));
  g_log.debug() << "There are " << vecLines.size() << " lines"
                << "\n";

  // determine the number of pixels at X direction (bear in mind that the XML
  // file records data in column major)
714
715
  size_t num_empty_line = 0;
  size_t num_weird_line = 0;
716
  for (size_t iline = 0; iline < vecLines.size(); ++iline) {
717
    if (vecLines[iline].empty())
718
      ++num_empty_line;
719
    else if (vecLines[iline].size() < 100)
720
      ++num_weird_line;
721
722
723
724
  }
  size_t num_pixel_x = vecLines.size() - num_empty_line - num_weird_line;
  g_log.information() << "There are " << num_empty_line << " lines and "
                      << num_weird_line << " lines are not regular.\n";
725
726

  // read the first line to determine the number of pixels at X direction
727
728
  size_t first_regular_line = 0;
  if (vecLines[first_regular_line].size() < 100)
729
    ++first_regular_line;
730
  std::vector<std::string> veccounts;
731
732
  boost::split(veccounts, vecLines[first_regular_line],
               boost::algorithm::is_any_of(" \t"));
733
734
735
736
  size_t num_pixel_y = veccounts.size();

  // create output workspace
  MatrixWorkspace_sptr outws;
737

738
739
740
741
742
743
744
745
746
  if (loadinstrument) {
    size_t numspec = num_pixel_x * num_pixel_y;
    outws = boost::dynamic_pointer_cast<MatrixWorkspace>(
        WorkspaceFactory::Instance().create("Workspace2D", numspec, 2, 1));
  } else {
    outws = boost::dynamic_pointer_cast<MatrixWorkspace>(
        WorkspaceFactory::Instance().create("Workspace2D", num_pixel_y,
                                            num_pixel_x, num_pixel_x));
  }
747

748
749
750
751
  // XML file records data in the order of column-major
  // FIXME - This may waste the previous result by parsing first line
  size_t i_col = 0;
  max_counts = 0;
752
  for (size_t i = first_regular_line; i < vecLines.size(); ++i) {
753
    std::string &line = vecLines[i];
754

755
756
757
758
    // skip empty lines
    if (line.size() < 100)
      continue;

759
760
761
762
763
    // Skip empty line
    if (line.empty()) {
      g_log.debug() << "\tFound empty Line at " << i << "\n";
      continue;
    }
764

765
766
767
768
769
770
771
772
    // Check whether it exceeds boundary
    if (i_col == num_pixel_x) {
      std::stringstream errss;
      errss << "Number of non-empty rows (" << i_col + 1
            << ") in detector data "
            << "exceeds user defined geometry size " << num_pixel_x << ".";
      throw std::runtime_error(errss.str());
    }
773

774
775
776
777
778
779
780
781
782
783
784
785
786
    // Split line
    std::vector<std::string> veccounts;
    boost::split(veccounts, line, boost::algorithm::is_any_of(" \t"));

    // check number of counts per column should not exceeds number of pixels
    // in Y direction
    if (veccounts.size() != num_pixel_y) {
      std::stringstream errss;
      errss << "Row " << i_col << " contains " << veccounts.size()
            << " items other than " << num_pixel_y
            << " counts specified by user.";
      throw std::runtime_error(errss.str());
    }
787

788
789
    // scan per column
    for (size_t j_row = 0; j_row < veccounts.size(); ++j_row) {
Peterson, Peter's avatar
Peterson, Peter committed
790
      double counts = std::stod(veccounts[j_row]);
791
      size_t rowIndex, columnIndex;
792

793
794
795
796
797
798
799
800
      if (loadinstrument) {
        // the detector ID and ws index are set up in column-major too!
        rowIndex = i_col * num_pixel_x + j_row;
        columnIndex = 0;
      } else {
        rowIndex = j_row;
        columnIndex = i_col;
      }
801

802
803
      outws->mutableX(rowIndex)[columnIndex] = static_cast<double>(columnIndex);
      outws->mutableY(rowIndex)[columnIndex] = counts;
804

805
806
807
808
      if (counts > 0)
        outws->mutableE(rowIndex)[columnIndex] = sqrt(counts);
      else
        outws->mutableE(rowIndex)[columnIndex] = 1.0;
809

810
811
812
813
814
      // record max count
      if (counts > max_counts) {
        max_counts = counts;
      }
    }
815

816
817
818
    // Update column index (i.e., column number)
    i_col += 1;
  } // END-FOR (i-vec line)
819

820
  return outws;
821
822
}

823
824
825
826
827
828
829
830
831
832
/** Set up sample logs from table workspace loaded where SPICE data file is
 * loaded
 * @brief LoadSpiceXML2DDet::setupSampleLogFromSpiceTable
 * @param matrixws
 * @param spicetablews
 * @param ptnumber
 */
void LoadSpiceXML2DDet::setupSampleLogFromSpiceTable(
    MatrixWorkspace_sptr matrixws, ITableWorkspace_sptr spicetablews,
    int ptnumber) {
833
834
835
836
837
  size_t numrows = spicetablews->rowCount();
  std::vector<std::string> colnames = spicetablews->getColumnNames();
  // FIXME - Shouldn't give a better value?
  Kernel::DateAndTime anytime(1000);

838
  bool foundlog = false;
839
  for (size_t ir = 0; ir < numrows; ++ir) {
840
    // loop over the table workspace to find the row of the spcified pt number
841
842
843
844
    int localpt = spicetablews->cell<int>(ir, 0);
    if (localpt != ptnumber)
      continue;

845
    // set the properties to matrix workspace including all columns
846
    for (size_t ic = 1; ic < colnames.size(); ++ic) {
847
848
      double logvalue = spicetablews->cell<double>(ir, ic);
      std::string &logname = colnames[ic];
849
      auto newlogproperty = new TimeSeriesProperty<double>(logname);
850
851
852
853
      newlogproperty->addValue(anytime, logvalue);
      matrixws->mutableRun().addProperty(newlogproperty);
    }

854
855
    // Break as the experiment pointer is found
    foundlog = true;
856
857
858
    break;
  }

859
860
861
862
  if (!foundlog)
    g_log.warning() << "Pt. " << ptnumber
                    << " is not found.  Log is not loaded to output workspace."
                    << "\n";
863
864
}

865
866
867
868
869
870
/** Get wavelength if the instrument is HB3A
 * @brief LoadSpiceXML2DDet::getHB3AWavelength
 * @param dataws
 * @param wavelength
 * @return
 */
871
872
873
bool LoadSpiceXML2DDet::getHB3AWavelength(MatrixWorkspace_sptr dataws,
                                          double &wavelength) {
  bool haswavelength(false);
874
  wavelength = -1.;
875

876
877
878
  // FIXME - Now it only search for _m1.  In future,
  //         it is better to searc both m1 and _m1

879
  if (dataws->run().hasProperty("_m1")) {
880
    g_log.notice("[DB] Data workspace has property _m1!");
881
882
883
    Kernel::TimeSeriesProperty<double> *ts =
        dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
            dataws->run().getProperty("_m1"));
884

885
886
887
888
889
    if (ts && ts->size() > 0) {
      double m1pos = ts->valuesAsVector()[0];
      if (fabs(m1pos - (-25.870000)) < 0.2) {
        wavelength = 1.003;
        haswavelength = true;
890
891
892
      } else if (fabs(m1pos - (-39.17)) < 0.2) {
        wavelength = 1.5424;
        haswavelength = true;
893
894
895
896
897
898
899
900
901
      } else {
        g_log.warning() << "m1 position " << m1pos
                        << " does not have defined mapping to "
                        << "wavelength."
                        << "\n";
      }
    } else if (!ts) {
      g_log.warning("Log _m1 is not TimeSeriesProperty.  Treat it as a single "
                    "value property.");
Peterson, Peter's avatar
Peterson, Peter committed
902
      double m1pos = std::stod(dataws->run().getProperty("_m1")->value());
903
904
905
      if (fabs(m1pos - (-25.870000)) < 0.2) {
        wavelength = 1.003;
        haswavelength = true;
906
907
908
      } else if (fabs(m1pos - (-39.17)) < 0.2) {
        wavelength = 1.5424;
        haswavelength = true;
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
      } else {
        g_log.warning() << "m1 position " << m1pos
                        << " does not have defined mapping to "
                        << "wavelength."
                        << "\n";
      }
    } else {
      g_log.error("Log _m1 is empty.");
    }
  } else {
    g_log.warning() << "No _m1 log is found."
                    << "\n";
  }

  if (!haswavelength)
    g_log.warning("No wavelength is setup!");
  else
    g_log.notice() << "[DB] Wavelength = " << wavelength << "\n";

  return haswavelength;
}

931
932
933
934
935
/** Set x axis to momentum (lab frame Q)
 * @brief LoadSpiceXML2DDet::setXtoLabQ
 * @param dataws
 * @param wavelength
 */
936
937
938
939
940
941
void LoadSpiceXML2DDet::setXtoLabQ(API::MatrixWorkspace_sptr dataws,
                                   const double &wavelength) {

  size_t numspec = dataws->getNumberHistograms();
  for (size_t iws = 0; iws < numspec; ++iws) {
    double ki = 2. * M_PI / wavelength;
942
943
944
    auto &x = dataws->mutableX(iws);
    x[0] = ki;
    x[1] = ki + 0.00001;
945
946
947
948
949
  }

  dataws->getAxis(0)->setUnit("Momentum");
}

950
951
952
953
954
955
956
/** Load instrument
 * @brief LoadSpiceXML2DDet::loadInstrument
 * @param matrixws
 * @param idffilename
 */
void LoadSpiceXML2DDet::loadInstrument(API::MatrixWorkspace_sptr matrixws,
                                       const std::string &idffilename) {
957
958
959
960
  // load instrument
  API::IAlgorithm_sptr loadinst = createChildAlgorithm("LoadInstrument");
  loadinst->initialize();
  loadinst->setProperty("Workspace", matrixws);
961
  if (!idffilename.empty()) {
962
    loadinst->setProperty("Filename", idffilename);
963
  } else
964
    loadinst->setProperty("InstrumentName", "HB3A");
965
  loadinst->setProperty("RewriteSpectraMap",
966
                        Mantid::Kernel::OptionalBool(true));
967
968
969
970
971
972
973
  loadinst->execute();
  if (loadinst->isExecuted())
    matrixws = loadinst->getProperty("Workspace");
  else
    g_log.error("Unable to load instrument to output workspace");
}

974
975
} // namespace DataHandling
} // namespace Mantid