Mantid
Loading...
Searching...
No Matches
LoadSpiceXML2DDet.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"
12#include "MantidAPI/Run.h"
18
19#include <boost/algorithm/string.hpp>
20
21#include <Poco/DOM/AutoPtr.h>
22#include <Poco/DOM/DOMParser.h>
23#include <Poco/DOM/Document.h>
24#include <Poco/DOM/NamedNodeMap.h>
25#include <Poco/DOM/Node.h>
26#include <Poco/DOM/NodeFilter.h>
27#include <Poco/DOM/NodeIterator.h>
28#include <Poco/DOM/NodeList.h>
29#include <Poco/SAX/InputSource.h>
30
31#include <algorithm>
32#include <fstream>
33#include <utility>
34
35using namespace std;
36
37namespace Mantid::DataHandling {
38
39using namespace Mantid::API;
40using namespace Mantid::Kernel;
41
43
44const char STRING = 's';
45const char FLOAT32 = 'f';
46const char INT32 = 'i';
47
52SpiceXMLNode::SpiceXMLNode(std::string nodename) : m_name{std::move(nodename)}, m_typechar('s') {}
53
58void SpiceXMLNode::setValue(const std::string &strvalue) { m_value = strvalue; }
59
66void SpiceXMLNode::setParameters(const std::string &nodetype, const std::string &nodeunit,
67 const std::string &nodedescription) {
68 // data type
69 if (nodetype == "FLOAT32") {
70 m_typefullname = nodetype;
72 } else if (nodetype == "INT32") {
73 m_typefullname = nodetype;
75 }
76
77 // unit
78 if (!nodeunit.empty()) {
79 m_unit = nodeunit;
80 }
81
82 // description
83 if (!nodedescription.empty())
84 m_description = nodedescription;
85}
86
89bool SpiceXMLNode::hasUnit() const { return (!m_unit.empty()); }
90
95bool SpiceXMLNode::hasValue() const { return (!m_value.empty()); }
96
101bool SpiceXMLNode::isString() const { return (m_typechar == STRING); }
102
107bool SpiceXMLNode::isInteger() const { return (m_typechar == INT32); }
108
113bool SpiceXMLNode::isDouble() const { return (m_typechar == FLOAT32); }
114
119const std::string SpiceXMLNode::getName() const { return m_name; }
120
125const std::string SpiceXMLNode::getUnit() const { return m_unit; }
126
131const std::string SpiceXMLNode::getDescription() const { return m_description; }
132
137const std::string SpiceXMLNode::getValue() const { return m_value; }
138
142 : m_detXMLFileName(), m_detXMLNodeName(), m_numPixelX(0), m_numPixelY(0), m_loadInstrument(false),
143 m_detSampleDistanceShift(0.0), m_hasScanTable(false), m_ptNumber4Log(0), m_idfFileName() {}
144
148
149const std::string LoadSpiceXML2DDet::name() const { return "LoadSpiceXML2DDet"; }
150
151int LoadSpiceXML2DDet::version() const { return 1; }
152
153const std::string LoadSpiceXML2DDet::category() const { return "DataHandling\\XML"; }
154
155const std::string LoadSpiceXML2DDet::summary() const {
156 return "Load 2-dimensional detector data file in XML format from SPICE. ";
157}
158
163 std::vector<std::string> exts;
164 exts.emplace_back(".xml");
165 exts.emplace_back(".bin");
166 declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::FileAction::Load, exts),
167 "XML file name for one scan including 2D detectors counts from SPICE");
168
169 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
170 "Name of output matrix workspace. "
171 "Output workspace will be an X by Y Workspace2D if instrument "
172 "is not loaded. ");
173
174 declareProperty("DetectorLogName", "Detector",
175 "Log name (i.e., XML node name) for detector counts in XML file."
176 "By default, the name is 'Detector'");
177
178 declareProperty(std::make_unique<ArrayProperty<size_t>>("DetectorGeometry"),
179 "A size-2 unsigned integer array [X, Y] for detector geometry. "
180 "Such that the detector contains X x Y pixels."
181 "If the input data is a binary file, input for DetectorGeometry will be "
182 "overridden "
183 "by detector geometry specified in the binary file");
184
185 declareProperty("LoadInstrument", true,
186 "Flag to load instrument to output workspace. "
187 "HFIR's HB3A will be loaded if InstrumentFileName is not specified.");
188
189 declareProperty(std::make_unique<FileProperty>("InstrumentFilename", "", FileProperty::OptionalLoad, ".xml"),
190 "The filename (including its full or relative path) of an instrument "
191 "definition file. The file extension must either be .xml or .XML when "
192 "specifying an instrument definition file. Note Filename or "
193 "InstrumentName must be specified but not both.");
194
195 declareProperty(std::make_unique<WorkspaceProperty<ITableWorkspace>>("SpiceTableWorkspace", "", Direction::Input,
197 "Name of TableWorkspace loaded from SPICE scan file by LoadSpiceAscii.");
198
199 declareProperty("PtNumber", 0, "Pt. value for the row to get sample log from. ");
200
201 declareProperty("UserSpecifiedWaveLength", EMPTY_DBL(),
202 "User can specify the wave length of the instrument if it is "
203 "drifted from the designed value."
204 "It happens often.");
205
206 declareProperty("ShiftedDetectorDistance", 0.,
207 "Amount of shift of the distance between source and detector centre."
208 "It is used to apply instrument calibration.");
209
210 declareProperty("DetectorCenterXShift", 0.0,
211 "The amount of shift of "
212 "detector center along X "
213 "direction in the unit meter.");
214
215 declareProperty("DetectorCenterYShift", 0.0,
216 "The amount of shift of "
217 "detector center along Y "
218 "direction in the unit meter.");
219}
220
226 m_detXMLNodeName = getPropertyValue("DetectorLogName");
227 std::vector<size_t> vec_pixelgeom = getProperty("DetectorGeometry");
228 if (vec_pixelgeom.size() == 2) {
229 m_numPixelX = vec_pixelgeom[0];
230 m_numPixelY = vec_pixelgeom[1];
231 } else if (vec_pixelgeom.empty()) {
232 m_numPixelX = 0;
233 m_numPixelY = 0;
234 } else {
235 throw std::runtime_error("Input pixels geometry is not correct in format. "
236 "It either has 2 integers or left empty to get "
237 "determined automatically.");
238 }
239 g_log.debug() << "User input poixels numbers: " << m_numPixelX << ", " << m_numPixelY << "\n";
240
241 m_loadInstrument = getProperty("LoadInstrument");
242
243 m_idfFileName = getPropertyValue("InstrumentFilename");
244 m_detSampleDistanceShift = getProperty("ShiftedDetectorDistance");
245
246 // Retreive sample environment data from SPICE scan table workspace
247 std::string spicetablewsname = getPropertyValue("SpiceTableWorkspace");
248 if (!spicetablewsname.empty())
249 m_hasScanTable = true;
250 else
251 m_hasScanTable = false;
252
253 m_ptNumber4Log = getProperty("PtNumber");
254
255 m_userSpecifiedWaveLength = getProperty("UserSpecifiedWaveLength");
256
257 m_detXShift = getProperty("DetectorCenterXShift");
258 m_detYShift = getProperty("DetectorCenterYShift");
259}
260
267 // With given spice scan table, 2-theta is read from there.
268 if (m_hasScanTable) {
269 ITableWorkspace_sptr spicetablews = getProperty("SpiceTableWorkspace");
270 setupSampleLogFromSpiceTable(outws, spicetablews, m_ptNumber4Log);
271 }
272
273 Types::Core::DateAndTime anytime(1000);
274
275 // Process 2theta
276 bool return_true = true;
277 if (!outws->run().hasProperty("2theta") && outws->run().hasProperty("_2theta")) {
278 // Set up 2theta if it is not set up yet
279 double logvalue = std::stod(outws->run().getProperty("_2theta")->value());
280 TimeSeriesProperty<double> *newlogproperty = new TimeSeriesProperty<double>("2theta");
281 newlogproperty->addValue(anytime, logvalue);
282 outws->mutableRun().addProperty(newlogproperty);
283 g_log.information() << "Set 2theta from _2theta (as XML node) with value " << logvalue << "\n";
284 } else if (!outws->run().hasProperty("2theta") && !outws->run().hasProperty("_2theta")) {
285 // Neither 2theta nor _2theta
286 g_log.warning("No 2theta is set up for loading instrument.");
287 return_true = false;
288 }
289
290 // set up the caibrated detector center to beam
292 det_dx->addValue(anytime, m_detXShift);
293 outws->mutableRun().addProperty(det_dx);
294
296 det_dy->addValue(anytime, m_detYShift);
297 outws->mutableRun().addProperty(det_dy);
298
299 // set up Sample-detetor distance calibration
300 double sampledetdistance = m_detSampleDistanceShift;
301 TimeSeriesProperty<double> *distproperty = new TimeSeriesProperty<double>("diffr");
302 distproperty->addValue(anytime, sampledetdistance);
303 outws->mutableRun().addProperty(distproperty);
304
305 return return_true;
306}
307
308//----------------------------------------------------------------------------------------------
313 // Load input
315
316 // check the file end
318 if (m_detXMLFileName.substr(m_detXMLFileName.find_last_of('.') + 1) == "bin") {
319 std::vector<unsigned int> vec_counts = binaryParseIntegers(m_detXMLFileName);
320 outws = createMatrixWorkspace(vec_counts);
321 } else {
322 // Parse detector XML file
323 std::vector<SpiceXMLNode> vec_xmlnode = xmlParseSpice(m_detXMLFileName);
324
325 // Create output workspace
326 if (m_numPixelX * m_numPixelY > 0)
329 else
331 }
332
333 // Set up log for loading instrument
334 bool can_set_instrument = setupSampleLogs(outws);
335
336 if (m_loadInstrument && can_set_instrument) {
338 // set wave length to user specified wave length
339 double wavelength = m_userSpecifiedWaveLength;
340 // if user does not specify wave length then try to get wave length from log
341 // sample _m1 (or m1 as well in future)
342 bool has_wavelength = !(wavelength == EMPTY_DBL());
343 if (!has_wavelength)
344 has_wavelength = getHB3AWavelength(outws, wavelength);
345
346 if (has_wavelength) {
347 setXtoLabQ(outws, wavelength);
348 }
349 }
350
351 setProperty("OutputWorkspace", outws);
352}
353
359std::vector<SpiceXMLNode> LoadSpiceXML2DDet::xmlParseSpice(const std::string &xmlfilename) {
360 // Declare output
361 std::vector<SpiceXMLNode> vecspicenode;
362
363 // Open file
364 std::ifstream ifs;
365 ifs.open(xmlfilename.c_str());
366 if (!ifs.is_open()) {
367 std::stringstream ess;
368 ess << "File " << xmlfilename << " cannot be opened.";
369 throw std::runtime_error(ess.str());
370 }
371
372 // Parse
373 Poco::XML::InputSource src(ifs);
374
375 Poco::XML::DOMParser parser;
376 Poco::AutoPtr<Poco::XML::Document> pDoc = parser.parse(&src);
377
378 // Go though XML
379 Poco::XML::NodeIterator nodeIter(pDoc, Poco::XML::NodeFilter::SHOW_ELEMENT);
380 Poco::XML::Node *pNode = nodeIter.nextNode();
381 while (pNode) {
382 const Poco::XML::XMLString nodename = pNode->nodeName();
383
384 // get number of children
385 Poco::AutoPtr<Poco::XML::NodeList> children = pNode->childNodes();
386 const size_t numchildren = children->length();
387 if (numchildren > 1) {
388 g_log.debug() << "Parent node " << nodename << " has " << numchildren << " children."
389 << "\n";
390 if (nodename == "SPICErack") {
391 // SPICErack is the main parent node. start_time and end_time are there
392 Poco::AutoPtr<Poco::XML::NamedNodeMap> attributes = pNode->attributes();
393 unsigned long numattr = attributes->length();
394 for (unsigned long j = 0; j < numattr; ++j) {
395 std::string attname = attributes->item(j)->nodeName();
396 std::string attvalue = attributes->item(j)->innerText();
397 SpiceXMLNode xmlnode(attname);
398 xmlnode.setValue(attvalue);
399 vecspicenode.emplace_back(xmlnode);
400 g_log.debug() << "SPICErack attribute " << j << " Name = " << attname << ", Value = " << attvalue << "\n";
401 }
402 }
403
404 } else if (numchildren == 1) {
405 std::string innertext = pNode->innerText();
406 Poco::AutoPtr<Poco::XML::NamedNodeMap> attributes = pNode->attributes();
407 unsigned long numattr = attributes->length();
408 g_log.debug() << " Child node " << nodename << "'s attributes: "
409 << "\n";
410
411 SpiceXMLNode xmlnode(nodename);
412 std::string nodetype;
413 std::string nodeunit;
414 std::string nodedescription;
415
416 for (unsigned long j = 0; j < numattr; ++j) {
417 std::string atttext = attributes->item(j)->innerText();
418 std::string attname = attributes->item(j)->nodeName();
419 g_log.debug() << " attribute " << j << " name = " << attname << ", "
420 << "value = " << atttext << "\n";
421 if (attname == "type") {
422 // type
423 nodetype = atttext;
424 } else if (attname == "unit") {
425 // unit
426 nodeunit = atttext;
427 } else if (attname == "description") {
428 // description
429 nodedescription = atttext;
430 }
431 }
432 xmlnode.setParameters(nodetype, nodeunit, nodedescription);
433 xmlnode.setValue(innertext);
434
435 vecspicenode.emplace_back(xmlnode);
436 } else {
437 // An unexpected case but no guarantee for not happening
438 g_log.error("Funny... No child node.");
439 }
440
441 // Move to next node
442 pNode = nodeIter.nextNode();
443 } // ENDWHILE
444
445 // Close file
446 ifs.close();
447
448 return vecspicenode;
449}
450
451//----------------------------------------------------------------------------------------------
453std::vector<unsigned int> LoadSpiceXML2DDet::binaryParseIntegers(std::string &binary_file_name) {
454 // check binary file size
455 ifstream infile(binary_file_name.c_str(), ios::binary);
456 streampos begin, end;
457 begin = infile.tellg();
458 infile.seekg(0, ios::end);
459 end = infile.tellg();
460 g_log.information() << "File size is: " << (end - begin) << " bytes.\n";
461
462 size_t num_unsigned_int = static_cast<size_t>(end - begin) / sizeof(unsigned int);
463 if (num_unsigned_int <= 2)
464 throw std::runtime_error("Input binary file size is too small (<= 2 unsigned int)");
465
466 size_t num_dets = num_unsigned_int - 2;
467 g_log.information() << "File contains " << num_unsigned_int << " unsigned integers and thus " << num_dets
468 << " detectors.\n";
469
470 // define output vector
471 std::vector<unsigned int> vec_counts(num_dets);
472 infile.seekg(0, ios::beg);
473
474 // read each integer... time consuming
475 // int max_count = 0;
476 // char buffer[sizeof(int)];
477 unsigned int buffer;
478 unsigned int total_counts(0);
479
480 // read detector size (row and column)
481 infile.read((char *)&buffer, sizeof(buffer));
482 auto num_rows = static_cast<size_t>(buffer);
483 infile.read((char *)&buffer, sizeof(buffer));
484 auto num_cols = static_cast<size_t>(buffer);
485 if (num_rows * num_cols != num_dets) {
486 g_log.error() << "Input binary file " << binary_file_name << " has inconsistent specification "
487 << "on detector size. "
488 << "First 2 unsigned integers are " << num_rows << ", " << num_cols
489 << ", while the detector number specified in the file is " << num_dets << "\n";
490 throw std::runtime_error("Input binary file has inconsistent specification on detector size.");
491 }
492
493 for (size_t i = 0; i < num_dets; ++i) {
494 // infile.read(buffer, sizeof(int));
495 infile.read((char *)&buffer, sizeof(buffer));
496 vec_counts[i] = buffer;
497 total_counts += buffer;
498 }
499
500 g_log.information() << "For detector " << num_rows << " x " << num_cols << ", total counts = " << total_counts
501 << "\n";
502
503 return vec_counts;
504}
505
506//----------------------------------------------------------------------------------------------
507MatrixWorkspace_sptr LoadSpiceXML2DDet::createMatrixWorkspace(const std::vector<unsigned int> &vec_counts) {
508 // Create matrix workspace
509 size_t numspec = vec_counts.size();
511 std::dynamic_pointer_cast<MatrixWorkspace>(WorkspaceFactory::Instance().create("Workspace2D", numspec, 2, 1));
512
513 g_log.information("Workspace created");
514
515 // set up value
516 for (size_t i = 0; i < numspec; ++i) {
517 outws->mutableX(i)[0] = 0.;
518 outws->mutableX(i)[1] = 1;
519 auto counts = static_cast<double>(vec_counts[i]);
520 outws->mutableY(i)[0] = counts;
521 if (counts > 0.5)
522 outws->mutableE(i)[0] = sqrt(counts);
523 else
524 outws->mutableE(i)[0] = 1.0;
525 }
526
527 return outws;
528}
529
530//-----
541LoadSpiceXML2DDet::xmlCreateMatrixWorkspaceKnownGeometry(const std::vector<SpiceXMLNode> &vecxmlnode,
542 const size_t &numpixelx, const size_t &numpixely,
543 const std::string &detnodename, const bool &loadinstrument) {
544
545 // TODO FIXME - If version 2 works, then this version will be discarded
546
547 // Create matrix workspace
549
550 if (loadinstrument) {
551 size_t numspec = numpixelx * numpixely;
552 outws =
553 std::dynamic_pointer_cast<MatrixWorkspace>(WorkspaceFactory::Instance().create("Workspace2D", numspec, 2, 1));
554 } else {
555 outws = std::dynamic_pointer_cast<MatrixWorkspace>(
556 WorkspaceFactory::Instance().create("Workspace2D", numpixely, numpixelx, numpixelx));
557 }
558
559 // Go through all XML nodes to process
560 size_t numxmlnodes = vecxmlnode.size();
561 bool parsedDet = false;
562 double max_counts = 0.;
563 for (size_t n = 0; n < numxmlnodes; ++n) {
564 // Process node for detector's count
565 const SpiceXMLNode &xmlnode = vecxmlnode[n];
566 if (xmlnode.getName() == detnodename) {
567 // Get node value string (256x256 as a whole)
568 const std::string detvaluestr = xmlnode.getValue();
569
570 // Split
571 std::vector<std::string> vecLines;
572 boost::split(vecLines, detvaluestr, boost::algorithm::is_any_of("\n"));
573 g_log.debug() << "There are " << vecLines.size() << " lines"
574 << "\n";
575
576 // XML file records data in the order of column-major
577 size_t i_col = 0;
578 for (size_t i = 0; i < vecLines.size(); ++i) {
579 std::string &line = vecLines[i];
580
581 // Skip empty line
582 if (line.empty()) {
583 g_log.debug() << "\tFound empty Line at " << i << "\n";
584 continue;
585 }
586
587 // Check whether it exceeds boundary
588 if (i_col == numpixelx) {
589 std::stringstream errss;
590 errss << "Number of non-empty rows (" << i_col + 1 << ") in detector data "
591 << "exceeds user defined geometry size " << numpixelx << ".";
592 throw std::runtime_error(errss.str());
593 }
594
595 // Split line
596 std::vector<std::string> veccounts;
597 boost::split(veccounts, line, boost::algorithm::is_any_of(" \t"));
598
599 // check number of counts per column should not exceeds number of pixels
600 // in Y direction
601 if (veccounts.size() != numpixely) {
602 std::stringstream errss;
603 errss << "[Version 1] Row " << i_col << " contains " << veccounts.size() << " items other than " << numpixely
604 << " counts specified by user.";
605 throw std::runtime_error(errss.str());
606 }
607
608 // scan per column
609 for (size_t j_row = 0; j_row < veccounts.size(); ++j_row) {
610 double counts = std::stod(veccounts[j_row]);
611 size_t rowIndex, columnIndex;
612
613 if (loadinstrument) {
614 // the detector ID and ws index are set up in column-major too!
615 rowIndex = i_col * numpixelx + j_row;
616 columnIndex = 0;
617 } else {
618 rowIndex = j_row;
619 columnIndex = i_col;
620 }
621
622 outws->mutableX(rowIndex)[columnIndex] = static_cast<double>(columnIndex);
623 outws->mutableY(rowIndex)[columnIndex] = counts;
624
625 if (counts > 0)
626 outws->mutableE(rowIndex)[columnIndex] = sqrt(counts);
627 else
628 outws->mutableE(rowIndex)[columnIndex] = 1.0;
629
630 // record max count
631 if (counts > max_counts) {
632 max_counts = counts;
633 }
634 }
635
636 // Update column index (i.e., column number)
637 i_col += 1;
638 } // END-FOR (i-vec line)
639
640 // Set flag
641 parsedDet = true;
642 } else {
643 // Parse to log: because there is no start time. so all logs are single
644 // value type
645 const std::string nodename = xmlnode.getName();
646 const std::string nodevalue = xmlnode.getValue();
647 if (xmlnode.isDouble()) {
648 double dvalue = std::stod(nodevalue);
649 outws->mutableRun().addProperty(new PropertyWithValue<double>(nodename, dvalue));
650 g_log.debug() << "Log name / xml node : " << xmlnode.getName() << " (double) value = " << dvalue << "\n";
651 } else if (xmlnode.isInteger()) {
652 int ivalue = std::stoi(nodevalue);
653 outws->mutableRun().addProperty(new PropertyWithValue<int>(nodename, ivalue));
654 g_log.debug() << "Log name / xml node : " << xmlnode.getName() << " (int) value = " << ivalue << "\n";
655 } else {
656 std::string str_value(nodevalue);
657 if (nodename == "start_time") {
658 // replace 2015-01-17 13:36:45 by 2015-01-17T13:36:45
659 str_value = nodevalue;
660 str_value.replace(10, 1, "T");
661 g_log.debug() << "Replace start_time " << nodevalue << " by Mantid time format " << str_value << "\n";
662 }
663 outws->mutableRun().addProperty(new PropertyWithValue<std::string>(nodename, str_value));
664 }
665 }
666 }
667
668 // Raise exception if no detector node is found
669 if (!parsedDet) {
670 std::stringstream errss;
671 errss << "Unable to find an XML node of name " << detnodename << ". Unable to load 2D detector XML file.";
672 throw std::runtime_error(errss.str());
673 }
674
675 g_log.notice() << "Maximum detector count on it is " << max_counts << "\n";
676
677 return outws;
678}
679
684LoadSpiceXML2DDet::xmlCreateMatrixWorkspaceUnknowGeometry(const std::vector<SpiceXMLNode> &vecxmlnode,
685 const std::string &detnodename, const bool &loadinstrument) {
686
687 // Create matrix workspace
689
690 // Go through all XML nodes to process
691 size_t numxmlnodes = vecxmlnode.size();
692 bool parsedDet = false;
693 double max_counts = 0.;
694
695 // define log value map
696 std::map<std::string, std::string> str_log_map;
697 std::map<std::string, double> dbl_log_map;
698 std::map<std::string, int> int_log_map;
699
700 for (size_t n = 0; n < numxmlnodes; ++n) {
701 // Process node for detector's count
702 const SpiceXMLNode &xmlnode = vecxmlnode[n];
703 if (xmlnode.getName() == detnodename) {
704 // Get node value string (256x256 as a whole)
705 const std::string detvaluestr = xmlnode.getValue();
706
707 outws = this->xmlParseDetectorNode(detvaluestr, loadinstrument, max_counts);
708
709 // Set flag
710 parsedDet = true;
711 } else {
712 // Parse to log: because there is no start time. so all logs are single
713 // value type
714 const std::string nodename = xmlnode.getName();
715 const std::string nodevalue = xmlnode.getValue();
716 if (xmlnode.isDouble()) {
717 double dvalue = std::stod(nodevalue);
718 dbl_log_map.emplace(nodename, dvalue);
719 } else if (xmlnode.isInteger()) {
720 int ivalue = std::stoi(nodevalue);
721 int_log_map.emplace(nodename, ivalue);
722 } else {
723 if (nodename == "start_time") {
724 // replace 2015-01-17 13:36:45 by 2015-01-17T13:36:45
725 std::string str_value(nodevalue);
726 str_value.replace(10, 1, "T");
727 g_log.debug() << "Replace start_time " << nodevalue << " by Mantid time format " << str_value << "\n";
728 str_log_map.emplace(nodename, str_value);
729 } else
730 str_log_map.emplace(nodename, nodevalue);
731 } // END-IF-ELSE (node value type)
732 } // END-IF-ELSE (detector-node or log node)
733 } // END-FOR (xml nodes)
734
735 if (outws) {
736 // Add the property to output workspace
737 for (auto &log_entry : str_log_map) {
738 outws->mutableRun().addProperty(new PropertyWithValue<std::string>(log_entry.first, log_entry.second));
739 }
740 for (auto &log_entry : int_log_map) {
741 outws->mutableRun().addProperty(new PropertyWithValue<int>(log_entry.first, log_entry.second));
742 }
743 for (auto &log_entry : dbl_log_map) {
744 outws->mutableRun().addProperty(new PropertyWithValue<double>(log_entry.first, log_entry.second));
745 }
746 }
747
748 // Raise exception if no detector node is found
749 if (!parsedDet) {
750 std::stringstream errss;
751 errss << "Unable to find an XML node of name " << detnodename << ". Unable to load 2D detector XML file.";
752 throw std::runtime_error(errss.str());
753 }
754
755 g_log.notice() << "Maximum detector count on it is " << max_counts << "\n";
756
757 return outws;
758}
759
760API::MatrixWorkspace_sptr LoadSpiceXML2DDet::xmlParseDetectorNode(const std::string &detvaluestr, bool loadinstrument,
761 double &max_counts) {
762 // Split to lines
763 std::vector<std::string> vecLines;
764 boost::split(vecLines, detvaluestr, boost::algorithm::is_any_of("\n"));
765 g_log.debug() << "There are " << vecLines.size() << " lines"
766 << "\n";
767
768 // determine the number of pixels at X direction (bear in mind that the XML
769 // file records data in column major)
770 size_t num_empty_line = 0;
771 size_t num_weird_line = 0;
772 for (auto &vecLine : vecLines) {
773 if (vecLine.empty())
774 ++num_empty_line;
775 else if (vecLine.size() < 100)
776 ++num_weird_line;
777 }
778 size_t num_pixel_x = vecLines.size() - num_empty_line - num_weird_line;
779 g_log.information() << "There are " << num_empty_line << " lines and " << num_weird_line
780 << " lines are not regular.\n";
781
782 // read the first line to determine the number of pixels at X direction
783 size_t first_regular_line = 0;
784 if (vecLines[first_regular_line].size() < 100)
785 ++first_regular_line;
786 std::vector<std::string> veccounts;
787 boost::split(veccounts, vecLines[first_regular_line], boost::algorithm::is_any_of(" \t"));
788 size_t num_pixel_y = veccounts.size();
789
790 // create output workspace
792
793 if (loadinstrument) {
794 size_t numspec = num_pixel_x * num_pixel_y;
795 outws =
796 std::dynamic_pointer_cast<MatrixWorkspace>(WorkspaceFactory::Instance().create("Workspace2D", numspec, 2, 1));
797 } else {
798 outws = std::dynamic_pointer_cast<MatrixWorkspace>(
799 WorkspaceFactory::Instance().create("Workspace2D", num_pixel_y, num_pixel_x, num_pixel_x));
800 }
801
802 // XML file records data in the order of column-major
803 // FIXME - This may waste the previous result by parsing first line
804 size_t i_col = 0;
805 max_counts = 0;
806 for (size_t i = first_regular_line; i < vecLines.size(); ++i) {
807 std::string &line = vecLines[i];
808
809 // skip empty lines
810 if (line.size() < 100)
811 continue;
812
813 // Skip empty line
814 if (line.empty()) {
815 g_log.debug() << "\tFound empty Line at " << i << "\n";
816 continue;
817 }
818
819 // Check whether it exceeds boundary
820 if (i_col == num_pixel_x) {
821 std::stringstream errss;
822 errss << "Number of non-empty rows (" << i_col + 1 << ") in detector data "
823 << "exceeds user defined geometry size " << num_pixel_x << ".";
824 throw std::runtime_error(errss.str());
825 }
826
827 boost::split(veccounts, line, boost::algorithm::is_any_of(" \t"));
828
829 // check number of counts per column should not exceeds number of pixels
830 // in Y direction
831 if (veccounts.size() != num_pixel_y) {
832 std::stringstream errss;
833 errss << "Row " << i_col << " contains " << veccounts.size() << " items other than " << num_pixel_y
834 << " counts specified by user.";
835 throw std::runtime_error(errss.str());
836 }
837
838 // scan per column
839 for (size_t j_row = 0; j_row < veccounts.size(); ++j_row) {
840 double counts = std::stod(veccounts[j_row]);
841 size_t rowIndex, columnIndex;
842
843 if (loadinstrument) {
844 // the detector ID and ws index are set up in column-major too!
845 rowIndex = i_col * num_pixel_x + j_row;
846 columnIndex = 0;
847 } else {
848 rowIndex = j_row;
849 columnIndex = i_col;
850 }
851
852 outws->mutableX(rowIndex)[columnIndex] = static_cast<double>(columnIndex);
853 outws->mutableY(rowIndex)[columnIndex] = counts;
854
855 if (counts > 0)
856 outws->mutableE(rowIndex)[columnIndex] = sqrt(counts);
857 else
858 outws->mutableE(rowIndex)[columnIndex] = 1.0;
859
860 // record max count
861 if (counts > max_counts) {
862 max_counts = counts;
863 }
864 }
865
866 // Update column index (i.e., column number)
867 i_col += 1;
868 } // END-FOR (i-vec line)
869
870 return outws;
871}
872
881 const ITableWorkspace_sptr &spicetablews, int ptnumber) {
882 size_t numrows = spicetablews->rowCount();
883 std::vector<std::string> colnames = spicetablews->getColumnNames();
884 // FIXME - Shouldn't give a better value?
885 Types::Core::DateAndTime anytime(1000);
886
887 bool foundlog = false;
888 for (size_t ir = 0; ir < numrows; ++ir) {
889 // loop over the table workspace to find the row of the spcified pt number
890 int localpt = spicetablews->cell<int>(ir, 0);
891 if (localpt != ptnumber)
892 continue;
893
894 // set the properties to matrix workspace including all columns
895 for (size_t ic = 1; ic < colnames.size(); ++ic) {
896 double logvalue = spicetablews->cell<double>(ir, ic);
897 std::string &logname = colnames[ic];
898 auto newlogproperty = new TimeSeriesProperty<double>(logname);
899 newlogproperty->addValue(anytime, logvalue);
900 matrixws->mutableRun().addProperty(newlogproperty);
901 }
902
903 // Break as the experiment pointer is found
904 foundlog = true;
905 break;
906 }
907
908 if (!foundlog)
909 g_log.warning() << "Pt. " << ptnumber << " is not found. Log is not loaded to output workspace."
910 << "\n";
911}
912
919bool LoadSpiceXML2DDet::getHB3AWavelength(const MatrixWorkspace_sptr &dataws, double &wavelength) {
920 bool haswavelength(false);
921 wavelength = -1.;
922
923 // FIXME - Now it only search for _m1. In future,
924 // it is better to searc both m1 and _m1
925
926 if (dataws->run().hasProperty("_m1")) {
927 g_log.notice("[DB] Data workspace has property _m1!");
929 dynamic_cast<Kernel::TimeSeriesProperty<double> *>(dataws->run().getProperty("_m1"));
930
931 if (ts && ts->size() > 0) {
932 double m1pos = ts->valuesAsVector()[0];
933 if (fabs(m1pos - (-25.870000)) < 0.2) {
934 wavelength = 1.003;
935 haswavelength = true;
936 } else if (fabs(m1pos - (-39.17)) < 0.2) {
937 wavelength = 1.5424;
938 haswavelength = true;
939 } else {
940 g_log.warning() << "m1 position " << m1pos << " does not have defined mapping to "
941 << "wavelength."
942 << "\n";
943 }
944 } else if (!ts) {
945 g_log.warning("Log _m1 is not TimeSeriesProperty. Treat it as a single "
946 "value property.");
947 double m1pos = std::stod(dataws->run().getProperty("_m1")->value());
948 if (fabs(m1pos - (-25.870000)) < 0.2) {
949 wavelength = 1.003;
950 haswavelength = true;
951 } else if (fabs(m1pos - (-39.17)) < 0.2) {
952 wavelength = 1.5424;
953 haswavelength = true;
954 } else {
955 g_log.warning() << "m1 position " << m1pos << " does not have defined mapping to "
956 << "wavelength."
957 << "\n";
958 }
959 } else {
960 g_log.error("Log _m1 is empty.");
961 }
962 } else {
963 g_log.warning() << "No _m1 log is found."
964 << "\n";
965 }
966
967 if (!haswavelength)
968 g_log.warning("No wavelength is setup!");
969 else
970 g_log.notice() << "[DB] Wavelength = " << wavelength << "\n";
971
972 return haswavelength;
973}
974
980void LoadSpiceXML2DDet::setXtoLabQ(const API::MatrixWorkspace_sptr &dataws, const double &wavelength) {
981
982 size_t numspec = dataws->getNumberHistograms();
983 for (size_t iws = 0; iws < numspec; ++iws) {
984 double ki = 2. * M_PI / wavelength;
985 auto &x = dataws->mutableX(iws);
986 x[0] = ki;
987 x[1] = ki + 0.00001;
988 }
989
990 dataws->getAxis(0)->setUnit("Momentum");
991}
992
998void LoadSpiceXML2DDet::loadInstrument(const API::MatrixWorkspace_sptr &matrixws, const std::string &idffilename) {
999 // load instrument
1000 auto loadinst = createChildAlgorithm("LoadInstrument");
1001 loadinst->initialize();
1002 loadinst->setProperty("Workspace", matrixws);
1003 if (!idffilename.empty()) {
1004 loadinst->setProperty("Filename", idffilename);
1005 } else
1006 loadinst->setProperty("InstrumentName", "HB3A");
1007 loadinst->setProperty("RewriteSpectraMap", Mantid::Kernel::OptionalBool(true));
1008 loadinst->execute();
1009 if (!loadinst->isExecuted())
1010 g_log.error("Unable to load instrument to output workspace");
1011}
1012
1013} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define fabs(x)
Definition: Matrix.cpp:22
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
A property class for workspaces.
LoadSpiceXML2DDet : Load 2D detector data in XML format form SPICE.
std::string m_detXMLNodeName
XML node name in detector counts file.
std::vector< unsigned int > binaryParseIntegers(std::string &binary_file_name)
parse binary integer file
void init() override
Declare properties.
const std::string summary() const override
Summary.
double m_detSampleDistanceShift
shift distance from sample to detector center
API::MatrixWorkspace_sptr xmlCreateMatrixWorkspaceKnownGeometry(const std::vector< SpiceXMLNode > &vecxmlnode, const size_t &numpixelx, const size_t &numpixely, const std::string &detnodename, const bool &loadinstrument)
Create output MatrixWorkspace.
void loadInstrument(const API::MatrixWorkspace_sptr &matrixws, const std::string &idffilename)
Load instrument.
void setupSampleLogFromSpiceTable(const API::MatrixWorkspace_sptr &matrixws, const API::ITableWorkspace_sptr &spicetablews, int ptnumber)
Set up sample logs from table workspace loaded where SPICE data file is loaded.
bool m_hasScanTable
Flag to show whether the SPICE scan table workspace is given.
int m_ptNumber4Log
Pt number for the sample logs to load with presense of Spice scan table workspace.
bool m_loadInstrument
Flag to show whether instrument is required to load.
~LoadSpiceXML2DDet() override
Destructor.
std::string m_detXMLFileName
SPICE detector XML file.
void exec() override
Main execution.
size_t m_numPixelY
Pixel size at Y direction.
std::string m_idfFileName
IDF file name to override Mantid's.
double m_detXShift
shift of detector on X and Y direction
API::MatrixWorkspace_sptr xmlCreateMatrixWorkspaceUnknowGeometry(const std::vector< SpiceXMLNode > &vecxmlnode, const std::string &detnodename, const bool &loadinstrument)
Create output MatrixWorkspace.
void setXtoLabQ(const API::MatrixWorkspace_sptr &dataws, const double &wavelength)
Set output workspace's X-axs as lab-frame Q space.
size_t m_numPixelX
Pixel size at X direction.
API::MatrixWorkspace_sptr createMatrixWorkspace(const std::vector< unsigned int > &vec_counts)
create workspace (good to load instrument) from vector of counts
double m_userSpecifiedWaveLength
User specified wave length.
int version() const override
Algorithm version.
bool getHB3AWavelength(const API::MatrixWorkspace_sptr &dataws, double &wavelength)
Get wavelength from workspace.
const std::string name() const override
Algoriothm name.
std::vector< SpiceXMLNode > xmlParseSpice(const std::string &xmlfilename)
Parse SPICE XML file.
API::MatrixWorkspace_sptr xmlParseDetectorNode(const std::string &detvaluestr, bool loadinstrument, double &max_counts)
const std::string category() const override
Category.
bool setupSampleLogs(const API::MatrixWorkspace_sptr &outws)
Set up sample logs in the output workspace.
bool hasValue() const
Check whether XML node has value set.
bool isDouble() const
Is this node of double type?
bool isInteger() const
Is this node of integer type?
void setValue(const std::string &strvalue)
Set node value in string format.
void setParameters(const std::string &nodetype, const std::string &nodeunit, const std::string &nodedescription)
Set XML node parameters.
const std::string getName() const
Get name of XML node.
bool hasUnit() const
Check whether XML has unit set.
const std::string getValue() const
Get node's value in string.
bool isString() const
Is this node of string type?
const std::string getUnit() const
Get unit of XML node.
const std::string getDescription() const
Get node's description.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
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
OptionalBool : Tri-state bool.
Definition: OptionalBool.h:25
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A specialised Property class for holding a series of time-value pairs.
int size() const override
Returns the number of values at UNIQUE time intervals in the time series.
std::vector< TYPE > valuesAsVector() const
Return the time series's values (unfiltered) as a vector<TYPE>
void addValue(const Types::Core::DateAndTime &time, const TYPE &value)
Add a value to the map using a DateAndTime object.
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
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.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54