Mantid
Loading...
Searching...
No Matches
LoadGSASInstrumentFile.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 +
11#include "MantidAPI/TableRow.h"
19
20#include <boost/algorithm/string.hpp>
21#include <boost/algorithm/string/finder.hpp>
22#include <boost/algorithm/string/iter_find.hpp>
23#include <boost/algorithm/string/predicate.hpp>
24#include <boost/algorithm/string/trim.hpp>
25
26#include <Poco/DOM/AutoPtr.h>
27#include <Poco/DOM/DOMWriter.h>
28#include <Poco/DOM/Element.h>
29
30#include <algorithm>
31#include <fstream>
32
33using namespace Mantid;
34using namespace Mantid::API;
35using namespace Mantid::DataObjects;
36using namespace Mantid::Kernel;
37using namespace std;
38using namespace Poco::XML;
39
42
43namespace Mantid::DataHandling {
44
46
47//----------------------------------------------------------------------------------------------
51 // Input file name
52 declareProperty(std::make_unique<FileProperty>("Filename", "", FileProperty::Load, ".prm"),
53 "Path to an GSAS file to load.");
54
55 // Output table workspace
56 auto wsprop = std::make_unique<WorkspaceProperty<API::ITableWorkspace>>("OutputTableWorkspace", "", Direction::Output,
58 declareProperty(std::move(wsprop), "Name of the output TableWorkspace containing "
59 "instrument parameter information read from file. ");
60
61 // Use bank numbers as given in file
62 declareProperty(std::make_unique<PropertyWithValue<bool>>("UseBankIDsInFile", true, Direction::Input),
63 "Use bank IDs as given in file rather than ordinal number of bank. "
64 "If the bank IDs in the file are not unique, it is advised to set this "
65 "to false.");
66
67 // Bank to import
68 declareProperty(std::make_unique<ArrayProperty<int>>("Banks"), "ID(s) of specified bank(s) to load, "
69 "The IDs are as specified by UseBankIDsInFile. "
70 "Default is all banks contained in input .prm file.");
71
72 // Workspace to put parameters into. It must be a workspace group with one
73 // worskpace per bank from the prm file
74 declareProperty(
76 "A workspace group with the instrument to which we add the "
77 "parameters from the GSAS instrument file, with one "
78 "workspace for each bank of the .prm file");
79
80 // Workspaces for each bank
81 declareProperty(std::make_unique<ArrayProperty<int>>("WorkspacesForBanks"),
82 "For each bank,"
83 " the ID of the corresponding workspace in same order as the "
84 "banks are specified. "
85 "ID=1 refers to the first workspace in the workspace group, "
86 "ID=2 refers to the second workspace and so on. "
87 "Default is all workspaces in numerical order. "
88 "If default banks are specified, they too are taken to be in "
89 "numerical order");
90}
91
92//----------------------------------------------------------------------------------------------
96 // Get input
97 string datafile = getProperty("Filename");
98
99 // Import data
100 vector<string> lines;
101 loadFile(datafile, lines);
102
103 // Check Histogram type - only PNTR is currently supported
104 std::string histType = getHistogramType(lines);
105 if (histType != "PNTR") {
106 throw std::runtime_error("Error on checking histogram type: " + histType + "\n");
107 }
108
109 size_t numBanks = getNumberOfBanks(lines);
110 g_log.debug() << numBanks << "banks in file \n";
111
112 // Examine bank information
113 vector<size_t> bankStartIndex;
114 scanBanks(lines, bankStartIndex);
115
116 if (bankStartIndex.empty()) {
117 throw std::runtime_error("No nanks found in file. \n");
118 }
119
120 if (numBanks != bankStartIndex.size()) { // The stated number of banks does
121 // not equal the number of banks
122 // found
123 g_log.warning() << "The number of banks found" << bankStartIndex.size()
124 << "is not equal to the number of banks stated" << numBanks << ".\n";
125 g_log.warning() << "Number of banks found is used.";
126 numBanks = bankStartIndex.size();
127 }
128
129 // Parse banks and export profile parameters
130 map<size_t, map<string, double>> bankparammap;
131 for (size_t i = 0; i < numBanks; ++i) {
132 size_t bankid = i + 1;
133 g_log.debug() << "Parse bank " << bankid << " of total " << numBanks << ".\n";
134 map<string, double> parammap;
135 parseBank(parammap, lines, bankid, bankStartIndex[bankid - 1]);
136 bankparammap.emplace(bankid, parammap);
137 g_log.debug() << "Bank starts at line" << bankStartIndex[i] + 1 << "\n";
138 }
139
140 // Get Workspace property
141 WorkspaceGroup_sptr wsg = getProperty("Workspace");
142 // Generate output table workspace
143 API::ITableWorkspace_sptr outTabWs = genTableWorkspace(bankparammap);
144 if (!getPropertyValue("OutputTableWorkspace").empty()) {
145 // Output the output table workspace
146 setProperty("OutputTableWorkspace", outTabWs);
147 }
148 if (wsg) {
149 vector<int> bankIds = getProperty("Banks");
150 vector<int> workspaceIds = getProperty("WorkspacesForBanks");
151 map<int, size_t> workspaceOfBank;
152
153 // Deal with bankIds
154 if (!bankIds.empty()) {
155 // If user provided a list of banks, check that they exist in the .prm
156 // file
157 for (auto bankId : bankIds) {
158 if (!bankparammap.count(bankId)) {
159 std::stringstream errorString;
160 errorString << "Bank " << bankId << " not found in .prm file";
161 throw runtime_error(errorString.str());
162 }
163 }
164 } else {
165 // Else, use all available banks
166 bankIds.reserve(bankparammap.size());
167 std::transform(bankparammap.cbegin(), bankparammap.cend(), std::back_inserter(bankIds),
168 [](const auto &bank) { return static_cast<int>(bank.first); });
169 }
170
171 // Generate workspaceOfBank
172 LoadFullprofResolution::createBankToWorkspaceMap(bankIds, workspaceIds, workspaceOfBank);
173 // Put parameters into workspace group
175 for (size_t i = 0; i < bankIds.size(); ++i) {
176 int bankId = bankIds[i];
177 size_t wsId = workspaceOfBank[bankId];
178 Workspace_sptr wsi = wsg->getItem(wsId - 1);
179 auto workspace = std::dynamic_pointer_cast<MatrixWorkspace>(wsi);
180 // Get column from table workspace
181 API::Column_const_sptr OutTabColumn = outTabWs->getColumn(i + 1);
182 std::string parameterXMLString;
184 OutTabColumn, workspace, static_cast<int>(bankparammap[i]["NPROF"]), parameterXMLString);
185 // Load the string into the workspace
186 Algorithm_sptr loadParamAlg = createChildAlgorithm("LoadParameterFile");
187 loadParamAlg->setProperty("ParameterXML", parameterXMLString);
188 loadParamAlg->setProperty("Workspace", workspace);
189 loadParamAlg->execute();
190 }
191 }
192}
193
194//----------------------------------------------------------------------------------------------
199void LoadGSASInstrumentFile::loadFile(const string &filename, vector<string> &lines) {
200 string line;
201
202 // the variable of type ifstream:
203 ifstream myfile(filename.c_str());
204
205 // check to see if the file is opened:
206 if (myfile.is_open()) {
207 // while there are still lines in the
208 // file, keep reading:
209 while (!myfile.eof()) {
210 // place the line from myfile into the
211 // line variable:
212 getline(myfile, line);
213
214 // display the line we gathered:
215 boost::algorithm::trim(line);
216 if (!line.empty())
217 lines.emplace_back(line);
218 }
219
220 // close the stream:
221 myfile.close();
222 } else {
223 stringstream errmsg;
224 errmsg << "Input .prm file " << filename << " cannot be open. ";
225 g_log.error(errmsg.str());
226 throw runtime_error(errmsg.str());
227 }
228}
229
230/* Get the histogram type
231 * @param lines :: vector of strings for each non-empty line in .irf file
232 * @return Histogram type code
233 */
234std::string LoadGSASInstrumentFile::getHistogramType(const vector<string> &lines) {
235 // We assume there is just one HTYPE line, look for it from beginning and
236 // return its value.
237 std::string lookFor = "INS HTYPE";
238 for (size_t i = 0; i < lines.size(); ++i) {
239 if (lines[i].substr(0, lookFor.size()) == lookFor) {
240 if (lines[i].size() < lookFor.size() + 7) {
241 // line too short
242 return "HTYPE line too short";
243 }
244 return lines[i].substr(lookFor.size() + 3, 4); // Found
245 }
246 }
247 return "HTYPE line not found";
248}
249
250/* Get number of banks
251 * @param lines :: vector of strings for each non-empty line in .irf file
252 * @return number of banks
253 */
254size_t LoadGSASInstrumentFile::getNumberOfBanks(const vector<string> &lines) {
255 // We assume there is just one BANK line, look for it from beginning and
256 // return its value.
257 std::string lookFor = "INS BANK";
258 for (size_t i = 0; i < lines.size(); ++i) {
259 if (lines[i].substr(0, lookFor.size()) == lookFor) {
260 if (lines[i].size() < lookFor.size() + 3) {
261 // line too short
262 return 0;
263 }
264 return boost::lexical_cast<size_t>(lines[i].substr(lookFor.size() + 2, 1)); // Found
265 }
266 }
267 return 0;
268}
269
275void LoadGSASInstrumentFile::scanBanks(const std::vector<std::string> &lines, std::vector<size_t> &bankStartIndex) {
276 // We look for each line that contains 'BNKPAR' and take it to be the first
277 // line of a bank.
278 // We currently ignore the bank number and assume they are numbered according
279 // to their order in the file.
280 for (size_t i = 0; i < lines.size(); ++i) {
281 string line = lines[i];
282 if (line.substr(0, 3) == "INS") { // Ignore all lines that don't begin with INS
283 if (line.find("BNKPAR") != string::npos) { // We've found start of a new bank
284 bankStartIndex.emplace_back(i);
285 }
286 } // INS
287 } // for(i)
288}
289
290//----------------------------------------------------------------------------------------------
298void LoadGSASInstrumentFile::parseBank(std::map<std::string, double> &parammap, const std::vector<std::string> &lines,
299 size_t bankid, size_t startlineindex) {
300 double param1, param2, param3, param4;
301
302 // We ignore the first lines of the bank
303 // The first useful line starts with "INS nPRCF", where n is the bank number
304 // From this line we get the profile function and number of parameters
305 size_t currentLineIndex = findINSPRCFLine(lines, startlineindex, param1, param2, param3, param4);
306
307 parammap["NPROF"] = param2;
308
309 // Ikeda-Carpenter PV
310 // Then read 15 parameters from the next four INS lines.
311 currentLineIndex = findINSPRCFLine(lines, currentLineIndex + 1, param1, param2, param3, param4);
312 parammap["Alph0"] = param1;
313 parammap["Alph1"] = param2;
314 parammap["Beta0"] = param3;
315 parammap["Beta1"] = param4; // Kappa
316
317 currentLineIndex = findINSPRCFLine(lines, currentLineIndex + 1, param1, param2, param3, param4);
318 parammap["Sig0"] = param1;
319 parammap["Sig1"] = param2;
320 parammap["Sig2"] = param3;
321 parammap["Gam0"] = param4;
322
323 findINSPRCFLine(lines, currentLineIndex + 1, param1, param2, param3, param4);
324 parammap["Gam1"] = param1;
325 parammap["Gam2"] = param2;
326 if (param3 != 0.0) {
327 g_log.warning() << "Bank" << bankid << "stec not 0, but " << param3;
328 }
329 if (param4 != 0.0) {
330 g_log.warning() << "Bank" << bankid << "ptec not 0, but " << param4;
331 }
332}
333
334//----------------------------------------------------------------------------------------------
345size_t LoadGSASInstrumentFile::findINSPRCFLine(const std::vector<std::string> &lines, size_t lineIndex, double &param1,
346 double &param2, double &param3, double &param4) {
347 for (size_t i = lineIndex; i < lines.size(); ++i) {
348 string line = lines[i];
349 if ((line.substr(0, 3) == "INS") && (line.substr(6, 4) == "PRCF")) {
350 std::istringstream paramLine;
351 paramLine.str(lines[i].substr(15));
352 paramLine >> param1 >> param2 >> param3 >> param4;
353 return i;
354 }
355 }
356 throw std::runtime_error("Unexpected end of file reached while searching for INS line. \n");
357}
358
359//----------------------------------------------------------------------------------------------
365TableWorkspace_sptr LoadGSASInstrumentFile::genTableWorkspace(map<size_t, map<string, double>> bankparammap) {
366 g_log.notice() << "Start to generate table workspace ...."
367 << ".\n";
368
369 // Retrieve some information
370 size_t numbanks = bankparammap.size();
371 if (numbanks == 0)
372 throw runtime_error("Unable to generate a table from an empty map!");
373
374 auto bankmapiter = bankparammap.begin();
375 size_t numparams = bankmapiter->second.size();
376
377 // vector of all parameter name
378 vector<string> vec_parname;
379 vector<size_t> vec_bankids;
380
381 map<string, double>::iterator parmapiter;
382 for (parmapiter = bankmapiter->second.begin(); parmapiter != bankmapiter->second.end(); ++parmapiter) {
383 string parname = parmapiter->first;
384 vec_parname.emplace_back(parname);
385 }
386
387 for (bankmapiter = bankparammap.begin(); bankmapiter != bankparammap.end(); ++bankmapiter) {
388 size_t bankid = bankmapiter->first;
389 vec_bankids.emplace_back(bankid);
390 }
391
392 g_log.debug() << "[DBx240] Number of imported parameters is " << numparams
393 << ", Number of banks = " << vec_bankids.size() << "."
394 << "\n";
395
396 // Create TableWorkspace
397 auto tablews = std::make_shared<TableWorkspace>();
398
399 // set columns :
400 // Any 2 columns cannot have the same name.
401 tablews->addColumn("str", "Name");
402 for (size_t i = 0; i < numbanks; ++i) {
403 stringstream colnamess;
404 size_t bankid = vec_bankids[i];
405 colnamess << "Value_" << bankid;
406 tablews->addColumn("double", colnamess.str());
407 }
408
409 g_log.debug() << "Number of column = " << tablews->columnCount() << ".\n";
410
411 // add BANK ID row
412 TableRow newrow = tablews->appendRow();
413 newrow << "BANK";
414 for (size_t i = 0; i < numbanks; ++i)
415 newrow << static_cast<double>(vec_bankids[i]);
416
417 g_log.debug() << "Number of row now = " << tablews->rowCount() << ".\n";
418
419 // add profile parameter rows
420 for (size_t i = 0; i < numparams; ++i) {
421 newrow = tablews->appendRow();
422
423 string parname = vec_parname[i];
424 newrow << parname;
425
426 for (size_t j = 0; j < numbanks; ++j) {
427 size_t bankid = vec_bankids[j];
428
429 // Locate map of bank 'bankid'
430 map<size_t, map<string, double>>::iterator bpmapiter;
431 bpmapiter = bankparammap.find(bankid);
432 if (bpmapiter == bankparammap.end()) {
433 throw runtime_error("Bank cannot be found in map.");
434 }
435
436 // Locate parameter
437 parmapiter = bpmapiter->second.find(parname);
438 if (parmapiter == bpmapiter->second.end()) {
439 throw runtime_error("Parameter cannot be found in a bank's map.");
440 } else {
441 double pvalue = parmapiter->second;
442 newrow << pvalue;
443 }
444
445 } // END(j)
446 } // END(i)
447
448 return tablews;
449}
450
451} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
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
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
static void getTableRowNumbers(const API::ITableWorkspace_sptr &tablews, std::map< std::string, size_t > &parammap)
Get row numbers of the parameters in the table workspace.
static void putParametersIntoWorkspace(const API::Column_const_sptr &, const API::MatrixWorkspace_sptr &ws, int nProf, std::string &parameterXMLString)
Put parameters into a matrix workspace.
static void createBankToWorkspaceMap(const std::vector< int > &banks, const std::vector< int > &workspaces, std::map< int, size_t > &workspaceOfBank)
Create Bank to Workspace Correspondence.
static std::map< std::string, size_t > m_rowNumbers
Place to store the row numbers.
LoadGSASInstrumentFile : Load GSAS instrument file to TableWorkspace(s)
std::string getHistogramType(const std::vector< std::string > &lines)
Get Histogram type.
size_t findINSPRCFLine(const std::vector< std::string > &lines, size_t lineIndex, double &param1, double &param2, double &param3, double &param4)
Find first INS line at or after lineIndex.
size_t getNumberOfBanks(const std::vector< std::string > &lines)
Get Number of banks.
DataObjects::TableWorkspace_sptr genTableWorkspace(std::map< size_t, std::map< std::string, double > > bankparammap)
Generate output workspace.
void scanBanks(const std::vector< std::string > &lines, std::vector< size_t > &bankStartIndex)
Scan imported file for bank information.
void exec() override
Implement abstract Algorithm methods.
void loadFile(const std::string &filename, std::vector< std::string > &lines)
Load file to a vector of strings.
void parseBank(std::map< std::string, double > &parammap, const std::vector< std::string > &lines, size_t bankid, size_t startlineindex)
Parse bank in file to a map.
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
The concrete, templated class for properties.
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< const Column > Column_const_sptr
Definition: Column.h:229
std::shared_ptr< Algorithm > Algorithm_sptr
Typedef for a shared pointer to an Algorithm.
Definition: Algorithm.h:61
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::shared_ptr< Instrument > Instrument_sptr
Shared pointer to an instrument object.
Helper class which provides the Collimation Length for SANS instruments.
STL namespace.
@ InOut
Both an input & output workspace.
Definition: Property.h:55
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54