Mantid
Loading...
Searching...
No Matches
SaveFullprofResolution.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 +
13
14#include "Poco/File.h"
15#include "boost/algorithm/string.hpp"
16#include "boost/math/special_functions/round.hpp"
17
18#include <fstream>
19#include <iomanip>
20
21using namespace Mantid;
22using namespace Mantid::API;
23using namespace Mantid::Kernel;
24using namespace Mantid::DataObjects;
25using namespace std;
26
27namespace Mantid::DataHandling {
28
30
31//----------------------------------------------------------------------------------------------
35 : API::Algorithm(), m_profileParamMap(), m_profileTableWS(), m_outIrfFilename(), m_bankID(-1),
36 m_fpProfileNumber(-1), m_append(false) {}
37
38//----------------------------------------------------------------------------------------------
42 declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("InputWorkspace", "", Direction::Input),
43 "Input TableWorkspace containing the parameters for .irf file.");
44
45 std::vector<std::string> exts{".irf"};
46 declareProperty(std::make_unique<API::FileProperty>("OutputFilename", "", API::FileProperty::Save, exts),
47 "Name of the output .irf file.");
48
49 std::shared_ptr<BoundedValidator<int>> bankboundval = std::make_shared<BoundedValidator<int>>();
50 bankboundval->setLower(0);
51 declareProperty("Bank", EMPTY_INT(), "Bank number of the parameters belonged to. ");
52
53 vector<string> supportedfunctions;
54 supportedfunctions.emplace_back("Back-to-back exponential convoluted with pseudo-voigt (profile 9)");
55 supportedfunctions.emplace_back("Jason Hodge's function (profile 10)");
56 auto funcvalidator = std::make_shared<StringListValidator>(supportedfunctions);
57 declareProperty("ProfileFunction", "Jason Hodge's function (profile 10)", funcvalidator,
58 "Profile number defined in Fullprof.");
59
60 declareProperty("Append", false,
61 "If true and the output file exists, the "
62 "bank will be appended to the existing "
63 "one.");
64}
65
66//----------------------------------------------------------------------------------------------
70 // Get input parameters
72
73 // Parse the input
75
76 // Generate the string for the file to write
77 std::string filestr;
78 switch (m_fpProfileNumber) {
79 case 9:
80 filestr = toProf9IrfString();
81 break;
82
83 case 10:
84 filestr = toProf10IrfString();
85 break;
86
87 default:
88 throw runtime_error("Profile number is not supported yet.");
89 }
90
91 // Write to file
92 std::ofstream ofile;
93 // Make it work!
94 if (m_append) {
95 ofile.open(m_outIrfFilename.c_str(), std::ofstream::out | std::ofstream::app);
96 g_log.information() << "Opened output file " << m_outIrfFilename << " in append mode. "
97 << "\n";
98 } else {
99 ofile.open(m_outIrfFilename.c_str(), std::ofstream::out | std::ofstream::trunc);
100 g_log.information() << "Opened output file " << m_outIrfFilename << " in new/overwrite mode. "
101 << "\n";
102 }
103 ofile << filestr;
104 ofile.close();
105}
106
107//----------------------------------------------------------------------------------------------
111 // Parameter table
112 m_profileTableWS = getProperty("InputWorkspace");
113
114 // Output file and operation
115 m_outIrfFilename = getPropertyValue("OutputFilename");
116 if (m_outIrfFilename.empty())
117 throw runtime_error("Input file name invalid. ");
118 m_append = getProperty("Append");
119 if (m_append) {
120 // Set append flag to false if file does not exist
121 bool fileexist = Poco::File(m_outIrfFilename).exists();
122 if (!fileexist)
123 m_append = false;
124 }
125
126 // Bank to write out
127 m_bankID = getProperty("Bank");
128
129 // Profile function
130 string proffunction = getProperty("ProfileFunction");
131 if (proffunction == "Back-to-back exponential convoluted with pseudo-voigt (profile 9)")
133 else if (proffunction == "Jason Hodge's function (profile 10)")
135 else {
136 stringstream errmsg;
137 errmsg << "It is impossible to have profile function " << proffunction << " input. ";
138 g_log.error(errmsg.str());
139 throw runtime_error(errmsg.str());
140 }
141}
142
143//----------------------------------------------------------------------------------------------
148 // Check the table workspace
149 std::vector<std::string> colnames = m_profileTableWS->getColumnNames();
150 size_t numcols = colnames.size();
151
152 stringstream dbmsgss("Input table's column names: ");
153 for (auto &colname : colnames) {
154 dbmsgss << setw(20) << colname;
155 }
156 g_log.debug(dbmsgss.str());
157
158 if (colnames[0] != "Name")
159 throw runtime_error("First colunm must be 'Name'");
160
161 // Read out a list of parameter names
162 size_t numpars = m_profileTableWS->rowCount();
163 vector<string> vec_parnames(numpars);
164 int rowbankindex = -1;
165 for (size_t i = 0; i < numpars; ++i) {
166 string parname = m_profileTableWS->cell<string>(i, 0);
167 vec_parnames[i] = parname;
168 if (parname == "BANK")
169 rowbankindex = static_cast<int>(i);
170 }
171
172 // Locate the column number to pass parameters
173 int colindex = -1;
174 if (rowbankindex < 0) {
175 // If there is NO 'BANK', locate first (from left) column starting with
176 // 'Value'
177 for (size_t i = 1; i < numcols; ++i) {
178 if (boost::starts_with(colnames[i], "Value")) {
179 colindex = static_cast<int>(i);
180 break;
181 }
182 }
183 } else {
184 // If there is BANK, Locate first (from left) column starting with 'Value'
185 // and BANK matches
186 for (size_t i = 1; i < numcols; ++i) {
187 if (boost::starts_with(colnames[i], "Value")) {
188 int bankid = boost::math::iround(m_profileTableWS->cell<double>(rowbankindex, i));
189 if (bankid == m_bankID) {
190 colindex = static_cast<int>(i);
191 break;
192 }
193 }
194 }
195 }
196
197 if (colindex < 0) {
198 throw runtime_error("Unable to find column");
199 } else if (colindex >= static_cast<int>(m_profileTableWS->columnCount())) {
200 throw runtime_error("Impossible to have this situation.");
201 }
202
203 // Clear the parameter
204 m_profileParamMap.clear();
205
206 // Parse
207 for (size_t ir = 0; ir < numpars; ++ir) {
208 double parvalue = m_profileTableWS->cell<double>(ir, static_cast<size_t>(colindex));
209 m_profileParamMap.emplace(vec_parnames[ir], parvalue);
210 }
211
212 // Debug output
213 stringstream dbss("Imported Parameter Table: \n");
214 map<string, double>::iterator mit;
215 for (mit = m_profileParamMap.begin(); mit != m_profileParamMap.end(); ++mit)
216 dbss << setw(20) << mit->first << " = " << setprecision(5) << mit->second << '\n';
217 g_log.debug(dbss.str());
218}
219
220//----------------------------------------------------------------------------------------------
224 // Get all parameter values
225 double tofmin = m_profileParamMap["tof-min"];
226 double tofmax = m_profileParamMap["tof-max"];
227 double zero = m_profileParamMap["Zero"];
228 double zerot = m_profileParamMap["Zerot"];
229 double tofstep = m_profileParamMap["step"];
230 double dtt1 = m_profileParamMap["Dtt1"];
231 double dtt1t = m_profileParamMap["Dtt1t"];
232 double dtt2t = m_profileParamMap["Dtt2t"];
233 double xcross = m_profileParamMap["Tcross"];
234 double width = m_profileParamMap["Width"];
235 double sig2 = m_profileParamMap["Sig2"];
236 double sig1 = m_profileParamMap["Sig1"];
237 double sig0 = m_profileParamMap["Sig0"];
238 double gam2 = m_profileParamMap["Gam2"];
239 double gam1 = m_profileParamMap["Gam1"];
240 double gam0 = m_profileParamMap["Gam0"];
241 double alph0 = m_profileParamMap["Alph0"];
242 double alph1 = m_profileParamMap["Alph1"];
243 double alph0t = m_profileParamMap["Alph0t"];
244 double alph1t = m_profileParamMap["Alph1t"];
245 double beta0 = m_profileParamMap["Beta0"];
246 double beta1 = m_profileParamMap["Beta1"];
247 double beta0t = m_profileParamMap["Beta0t"];
248 double beta1t = m_profileParamMap["Beta1t"];
249 int profindex = static_cast<int>(floor(m_profileParamMap["Profile"] + 0.5));
250 double twotheta = m_profileParamMap["twotheta"];
251
252 // Check with profile index
253 if (profindex == 0) {
254 profindex = 10;
255 } else if (profindex != 10) {
256 stringstream errmsg;
257 errmsg << "This column in table has profile number " << profindex << " other than 10.";
258 g_log.error(errmsg.str());
259 throw runtime_error(errmsg.str());
260 }
261
262 stringstream content;
263 content << fixed;
264
265 // Write header
266 if (!m_append) {
267 content << " Instrumental resolution function for POWGEN/SNS ireso: 6"
268 << "\n";
269 content << "! To be used with function NPROF=" << profindex << " in FullProf (Res=6)"
270 << "\n";
271 }
272
273 // Write bank information
274 content << "! ---------------------------------------------- Bank " << m_bankID << " ";
275 if (has_key(m_profileParamMap, "CWL")) {
276 double cwl = m_profileParamMap["CWL"];
277 if (cwl > 0)
278 content << "CWL = " << setprecision(4) << cwl << "A"
279 << "\n";
280 else
281 content << "\n";
282 } else {
283 content << "\n";
284 }
285
286 // Write profile parameter
287 content << "! Type of profile function: back-to-back exponentials * pseudo-Voigt"
288 << "\n";
289 content << "NPROF " << profindex << "\n";
290
291 content << "! Tof-min(us) step Tof-max(us)"
292 << "\n";
293 content << "TOFRG " << setprecision(3) << tofmin << " " << setw(16) << setprecision(5) << tofstep << " " << setw(16)
294 << setprecision(3) << tofmax << "\n";
295
296 content << "! Zero Dtt1"
297 << "\n";
298 content << "ZD2TOF " << setw(16) << setprecision(5) << zero << setw(16) << setprecision(5) << dtt1 << "\n";
299
300 content << "! Zerot Dtt1t Dtt2t x-cross Width"
301 << "\n";
302 content << "ZD2TOT " << setprecision(5) << zerot << setw(16) << setprecision(5) << dtt1t << setw(16)
303 << setprecision(5) << dtt2t << setw(16) << setprecision(10) << xcross << setw(16) << setprecision(5) << width
304 << "\n";
305
306 content << "! TOF-TWOTH of the bank"
307 << "\n";
308 content << "TWOTH " << setprecision(3) << twotheta << "\n";
309
310 // Note that sig0, sig1 and sig2 used in LeBail/Mantid framework are of the
311 // definition in manual.
312 // In .irf file, Sig-0, Sig-1 and Sig-2 are the squared values;
313 content << "! Sig-2 Sig-1 Sig-0"
314 << "\n";
315 content << "SIGMA " << setprecision(6) << sig2 * sig2 << setw(16) << setprecision(6) << sig1 * sig1 << setw(16)
316 << setprecision(6) << sig0 * sig0 << "\n";
317
318 content << "! Gam-2 Gam-1 Gam-0"
319 << "\n";
320 content << "GAMMA " << setw(16) << setprecision(6) << gam2 << " " << setw(16) << setprecision(6) << gam1 << " "
321 << setw(16) << setprecision(6) << gam0 << "\n";
322
323 content << "! alph0 beta0 alph1 beta1"
324 << "\n";
325 content << "ALFBE " << setprecision(6) << alph0 << " " << setw(16) << setprecision(6) << beta0 << " "
326 << setw(16) << setprecision(6) << alph1 << " " << setw(16) << setprecision(6) << beta1 << "\n";
327
328 content << "! alph0t beta0t alph1t beta1t"
329 << "\n";
330 content << "ALFBT " << setprecision(6) << alph0t << " " << setw(16) << setprecision(6) << beta0t << " "
331 << setw(16) << setprecision(6) << alph1t << " " << setw(16) << setprecision(6) << beta1t << "\n";
332 content << "END"
333 << "\n";
334
335 return content.str();
336}
337
338//----------------------------------------------------------------------------------------------
342 // Get all parameter values
343 double tofmin = m_profileParamMap["tof-min"];
344 double tofmax = m_profileParamMap["tof-max"];
345 double zero = m_profileParamMap["Zero"];
346 double tofstep = m_profileParamMap["step"];
347 double dtt1 = m_profileParamMap["Dtt1"];
348 double dtt2 = m_profileParamMap["Dtt2"];
349 double sig2 = m_profileParamMap["Sig2"];
350 double sig1 = m_profileParamMap["Sig1"];
351 double sig0 = m_profileParamMap["Sig0"];
352 double gam2 = m_profileParamMap["Gam2"];
353 double gam1 = m_profileParamMap["Gam1"];
354 double gam0 = m_profileParamMap["Gam0"];
355 double alph0 = m_profileParamMap["Alph0"];
356 double alph1 = m_profileParamMap["Alph1"];
357 double beta0 = m_profileParamMap["Beta0"];
358 double beta1 = m_profileParamMap["Beta1"];
359 int profindex = static_cast<int>(floor(m_profileParamMap["Profile"] + 0.5));
360 double twotheta = m_profileParamMap["twotheta"];
361 if (twotheta < 0)
362 twotheta += 360.;
363
364 // Check with profile index
365 if (profindex == 0) {
366 profindex = 9;
367 } else if (profindex != 9) {
368 stringstream errmsg;
369 errmsg << "This column in table has profile number " << profindex << " other than 9.";
370 g_log.error(errmsg.str());
371 throw runtime_error(errmsg.str());
372 }
373
374 stringstream content;
375 content << fixed;
376
377 // Write header
378 if (!m_append) {
379 content << "Instrumental resolution function for HRPD/ISIS L. Chapon "
380 "12/2003 ireso: 5"
381 << "\n";
382 content << "! To be used with function NPROF=" << profindex << " in FullProf (Res=5)"
383 << "\n";
384 }
385
386 // Write bank information
387 content << "! ---------------------------------------------- Bank " << m_bankID << " ";
388 if (has_key(m_profileParamMap, "CWL")) {
389 double cwl = m_profileParamMap["CWL"];
390 if (cwl > 0.)
391 content << "CWL = " << setprecision(4) << cwl << "A"
392 << "\n";
393 else
394 content << "\n";
395 } else {
396 content << "\n";
397 }
398
399 // Write profile parameters
400 content << "! Type of profile function: back-to-back exponentials * pseudo-Voigt"
401 << "\n";
402 content << "NPROF " << profindex << "\n";
403 content << "! Tof-min(us) step Tof-max(us)"
404 << "\n";
405 content << "TOFRG " << setprecision(3) << tofmin << " " << setw(16) << setprecision(5) << tofstep << " " << setw(16)
406 << setprecision(3) << tofmax << "\n";
407
408 content << "! Dtt1 Dtt2 Zero"
409 << "\n";
410 content << "D2TOF " << setw(16) << setprecision(5) << dtt1 << setw(16) << setprecision(5) << dtt2 << setw(16)
411 << setprecision(5) << zero << "\n";
412
413 content << "! TOF-TWOTH of the bank"
414 << "\n";
415 content << "TWOTH " << setprecision(3) << twotheta << "\n";
416
417 // Note that sig0, sig1 and sig2 used in LeBail/Mantid framework are of the
418 // definition in manual.
419 // In .irf file, Sig-0, Sig-1 and Sig-2 are the squared values;
420 content << "! Sig-2 Sig-1 Sig-0"
421 << "\n";
422 content << "SIGMA " << setprecision(6) << sig2 * sig2 << setw(16) << setprecision(6) << sig1 * sig1 << setw(16)
423 << setprecision(6) << sig0 * sig0 << "\n";
424
425 content << "! Gam-2 Gam-1 Gam-0"
426 << "\n";
427 content << "GAMMA " << setw(16) << setprecision(6) << gam2 << " " << setw(16) << setprecision(6) << gam1 << " "
428 << setw(16) << setprecision(6) << gam0 << "\n";
429
430 content << "! alph0 beta0 alph1 beta1"
431 << "\n";
432 content << "ALFBE " << setprecision(6) << alph0 << " " << setw(16) << setprecision(6) << beta0 << " "
433 << setw(16) << setprecision(6) << alph1 << " " << setw(16) << setprecision(6) << beta1 << "\n";
434
435 content << "END"
436 << "\n";
437
438 return content.str();
439}
440
441//
444bool SaveFullprofResolution::has_key(std::map<std::string, double> profmap, const std::string &key) {
445 map<string, double>::iterator fiter;
446 fiter = profmap.find(key);
447 bool exist = true;
448 if (fiter == profmap.end())
449 exist = false;
450
451 return exist;
452}
453
454} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
@ Save
to specify a file to write to, the file may or may not exist
Definition: FileProperty.h:49
A property class for workspaces.
SaveFullprofResolution : TODO: DESCRIPTION.
void parseTableWorkspace()
Parse input workspace to map of parameters.
std::string m_outIrfFilename
Output Irf file name.
void processProperties()
Write the header information.
std::string toProf10IrfString()
Write out string for profile 10 .irf file.
DataObjects::TableWorkspace_sptr m_profileTableWS
Input table workspace.
std::string toProf9IrfString()
Write out string for profile 9 .irf file.
bool has_key(std::map< std::string, double > profmap, const std::string &key)
Check wether a profile parameter map has the parameter.
std::map< std::string, double > m_profileParamMap
Map containing the name of value of each parameter required by .irf file.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Helper class which provides the Collimation Length for SANS instruments.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53