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 +
14
15#include "boost/math/special_functions/round.hpp"
16
17#include <filesystem>
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 = std::filesystem::exists(m_outIrfFilename);
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 (const 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 const auto it = std::find_if(colnames.begin() + 1, colnames.end(),
178 [](const std::string &colname) { return colname.starts_with("Value"); });
179 if (it != colnames.end()) {
180 colindex = static_cast<int>(std::distance(colnames.begin(), it));
181 }
182 } else {
183 // If there is BANK, Locate first (from left) column starting with 'Value'
184 // and BANK matches
185 for (size_t i = 1; i < numcols; ++i) {
186 if (colnames[i].starts_with("Value")) {
187 int bankid = boost::math::iround(m_profileTableWS->cell<double>(rowbankindex, i));
188 if (bankid == m_bankID) {
189 colindex = static_cast<int>(i);
190 break;
191 }
192 }
193 }
194 }
195
196 if (colindex < 0) {
197 throw runtime_error("Unable to find column");
198 } else if (colindex >= static_cast<int>(m_profileTableWS->columnCount())) {
199 throw runtime_error("Impossible to have this situation.");
200 }
201
202 // Clear the parameter
203 m_profileParamMap.clear();
204
205 // Parse
206 for (size_t ir = 0; ir < numpars; ++ir) {
207 double parvalue = m_profileTableWS->cell<double>(ir, static_cast<size_t>(colindex));
208 m_profileParamMap.emplace(vec_parnames[ir], parvalue);
209 }
210
211 // Debug output
212 stringstream dbss("Imported Parameter Table: \n");
213 map<string, double>::iterator mit;
214 for (mit = m_profileParamMap.begin(); mit != m_profileParamMap.end(); ++mit)
215 dbss << setw(20) << mit->first << " = " << setprecision(5) << mit->second << '\n';
216 g_log.debug(dbss.str());
217}
218
219//----------------------------------------------------------------------------------------------
223 // Get all parameter values
224 double tofmin = m_profileParamMap["tof-min"];
225 double tofmax = m_profileParamMap["tof-max"];
226 double zero = m_profileParamMap["Zero"];
227 double zerot = m_profileParamMap["Zerot"];
228 double tofstep = m_profileParamMap["step"];
229 double dtt1 = m_profileParamMap["Dtt1"];
230 double dtt1t = m_profileParamMap["Dtt1t"];
231 double dtt2t = m_profileParamMap["Dtt2t"];
232 double xcross = m_profileParamMap["Tcross"];
233 double width = m_profileParamMap["Width"];
234 double sig2 = m_profileParamMap["Sig2"];
235 double sig1 = m_profileParamMap["Sig1"];
236 double sig0 = m_profileParamMap["Sig0"];
237 double gam2 = m_profileParamMap["Gam2"];
238 double gam1 = m_profileParamMap["Gam1"];
239 double gam0 = m_profileParamMap["Gam0"];
240 double alph0 = m_profileParamMap["Alph0"];
241 double alph1 = m_profileParamMap["Alph1"];
242 double alph0t = m_profileParamMap["Alph0t"];
243 double alph1t = m_profileParamMap["Alph1t"];
244 double beta0 = m_profileParamMap["Beta0"];
245 double beta1 = m_profileParamMap["Beta1"];
246 double beta0t = m_profileParamMap["Beta0t"];
247 double beta1t = m_profileParamMap["Beta1t"];
248 int profindex = static_cast<int>(floor(m_profileParamMap["Profile"] + 0.5));
249 double twotheta = m_profileParamMap["twotheta"];
250
251 // Check with profile index
252 if (profindex == 0) {
253 profindex = 10;
254 } else if (profindex != 10) {
255 stringstream errmsg;
256 errmsg << "This column in table has profile number " << profindex << " other than 10.";
257 g_log.error(errmsg.str());
258 throw runtime_error(errmsg.str());
259 }
260
261 stringstream content;
262 content << fixed;
263
264 // Write header
265 if (!m_append) {
266 content << " Instrumental resolution function for POWGEN/SNS ireso: 6"
267 << "\n";
268 content << "! To be used with function NPROF=" << profindex << " in FullProf (Res=6)"
269 << "\n";
270 }
271
272 // Write bank information
273 content << "! ---------------------------------------------- Bank " << m_bankID << " ";
274 if (has_key(m_profileParamMap, "CWL")) {
275 double cwl = m_profileParamMap["CWL"];
276 if (cwl > 0)
277 content << "CWL = " << setprecision(4) << cwl << "A"
278 << "\n";
279 else
280 content << "\n";
281 } else {
282 content << "\n";
283 }
284
285 // Write profile parameter
286 content << "! Type of profile function: back-to-back exponentials * pseudo-Voigt"
287 << "\n";
288 content << "NPROF " << profindex << "\n";
289
290 content << "! Tof-min(us) step Tof-max(us)"
291 << "\n";
292 content << "TOFRG " << setprecision(3) << tofmin << " " << setw(16) << setprecision(5) << tofstep << " " << setw(16)
293 << setprecision(3) << tofmax << "\n";
294
295 content << "! Zero Dtt1"
296 << "\n";
297 content << "ZD2TOF " << setw(16) << setprecision(5) << zero << setw(16) << setprecision(5) << dtt1 << "\n";
298
299 content << "! Zerot Dtt1t Dtt2t x-cross Width"
300 << "\n";
301 content << "ZD2TOT " << setprecision(5) << zerot << setw(16) << setprecision(5) << dtt1t << setw(16)
302 << setprecision(5) << dtt2t << setw(16) << setprecision(10) << xcross << setw(16) << setprecision(5) << width
303 << "\n";
304
305 content << "! TOF-TWOTH of the bank"
306 << "\n";
307 content << "TWOTH " << setprecision(3) << twotheta << "\n";
308
309 // Note that sig0, sig1 and sig2 used in LeBail/Mantid framework are of the
310 // definition in manual.
311 // In .irf file, Sig-0, Sig-1 and Sig-2 are the squared values;
312 content << "! Sig-2 Sig-1 Sig-0"
313 << "\n";
314 content << "SIGMA " << setprecision(6) << sig2 * sig2 << setw(16) << setprecision(6) << sig1 * sig1 << setw(16)
315 << setprecision(6) << sig0 * sig0 << "\n";
316
317 content << "! Gam-2 Gam-1 Gam-0"
318 << "\n";
319 content << "GAMMA " << setw(16) << setprecision(6) << gam2 << " " << setw(16) << setprecision(6) << gam1 << " "
320 << setw(16) << setprecision(6) << gam0 << "\n";
321
322 content << "! alph0 beta0 alph1 beta1"
323 << "\n";
324 content << "ALFBE " << setprecision(6) << alph0 << " " << setw(16) << setprecision(6) << beta0 << " "
325 << setw(16) << setprecision(6) << alph1 << " " << setw(16) << setprecision(6) << beta1 << "\n";
326
327 content << "! alph0t beta0t alph1t beta1t"
328 << "\n";
329 content << "ALFBT " << setprecision(6) << alph0t << " " << setw(16) << setprecision(6) << beta0t << " "
330 << setw(16) << setprecision(6) << alph1t << " " << setw(16) << setprecision(6) << beta1t << "\n";
331 content << "END"
332 << "\n";
333
334 return content.str();
335}
336
337//----------------------------------------------------------------------------------------------
341 // Get all parameter values
342 double tofmin = m_profileParamMap["tof-min"];
343 double tofmax = m_profileParamMap["tof-max"];
344 double zero = m_profileParamMap["Zero"];
345 double tofstep = m_profileParamMap["step"];
346 double dtt1 = m_profileParamMap["Dtt1"];
347 double dtt2 = m_profileParamMap["Dtt2"];
348 double sig2 = m_profileParamMap["Sig2"];
349 double sig1 = m_profileParamMap["Sig1"];
350 double sig0 = m_profileParamMap["Sig0"];
351 double gam2 = m_profileParamMap["Gam2"];
352 double gam1 = m_profileParamMap["Gam1"];
353 double gam0 = m_profileParamMap["Gam0"];
354 double alph0 = m_profileParamMap["Alph0"];
355 double alph1 = m_profileParamMap["Alph1"];
356 double beta0 = m_profileParamMap["Beta0"];
357 double beta1 = m_profileParamMap["Beta1"];
358 int profindex = static_cast<int>(floor(m_profileParamMap["Profile"] + 0.5));
359 double twotheta = m_profileParamMap["twotheta"];
360 if (twotheta < 0)
361 twotheta += 360.;
362
363 // Check with profile index
364 if (profindex == 0) {
365 profindex = 9;
366 } else if (profindex != 9) {
367 stringstream errmsg;
368 errmsg << "This column in table has profile number " << profindex << " other than 9.";
369 g_log.error(errmsg.str());
370 throw runtime_error(errmsg.str());
371 }
372
373 stringstream content;
374 content << fixed;
375
376 // Write header
377 if (!m_append) {
378 content << "Instrumental resolution function for HRPD/ISIS L. Chapon "
379 "12/2003 ireso: 5"
380 << "\n";
381 content << "! To be used with function NPROF=" << profindex << " in FullProf (Res=5)"
382 << "\n";
383 }
384
385 // Write bank information
386 content << "! ---------------------------------------------- Bank " << m_bankID << " ";
387 if (has_key(m_profileParamMap, "CWL")) {
388 double cwl = m_profileParamMap["CWL"];
389 if (cwl > 0.)
390 content << "CWL = " << setprecision(4) << cwl << "A"
391 << "\n";
392 else
393 content << "\n";
394 } else {
395 content << "\n";
396 }
397
398 // Write profile parameters
399 content << "! Type of profile function: back-to-back exponentials * pseudo-Voigt"
400 << "\n";
401 content << "NPROF " << profindex << "\n";
402 content << "! Tof-min(us) step Tof-max(us)"
403 << "\n";
404 content << "TOFRG " << setprecision(3) << tofmin << " " << setw(16) << setprecision(5) << tofstep << " " << setw(16)
405 << setprecision(3) << tofmax << "\n";
406
407 content << "! Dtt1 Dtt2 Zero"
408 << "\n";
409 content << "D2TOF " << setw(16) << setprecision(5) << dtt1 << setw(16) << setprecision(5) << dtt2 << setw(16)
410 << setprecision(5) << zero << "\n";
411
412 content << "! TOF-TWOTH of the bank"
413 << "\n";
414 content << "TWOTH " << setprecision(3) << twotheta << "\n";
415
416 // Note that sig0, sig1 and sig2 used in LeBail/Mantid framework are of the
417 // definition in manual.
418 // In .irf file, Sig-0, Sig-1 and Sig-2 are the squared values;
419 content << "! Sig-2 Sig-1 Sig-0"
420 << "\n";
421 content << "SIGMA " << setprecision(6) << sig2 * sig2 << setw(16) << setprecision(6) << sig1 * sig1 << setw(16)
422 << setprecision(6) << sig0 * sig0 << "\n";
423
424 content << "! Gam-2 Gam-1 Gam-0"
425 << "\n";
426 content << "GAMMA " << setw(16) << setprecision(6) << gam2 << " " << setw(16) << setprecision(6) << gam1 << " "
427 << setw(16) << setprecision(6) << gam0 << "\n";
428
429 content << "! alph0 beta0 alph1 beta1"
430 << "\n";
431 content << "ALFBE " << setprecision(6) << alph0 << " " << setw(16) << setprecision(6) << beta0 << " "
432 << setw(16) << setprecision(6) << alph1 << " " << setw(16) << setprecision(6) << beta1 << "\n";
433
434 content << "END"
435 << "\n";
436
437 return content.str();
438}
439
440//
443bool SaveFullprofResolution::has_key(std::map<std::string, double> profmap, const std::string &key) {
444 map<string, double>::iterator fiter;
445 fiter = profmap.find(key);
446 bool exist = true;
447 if (fiter == profmap.end())
448 exist = false;
449
450 return exist;
451}
452
453} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition Algorithm.h:538
Base class from which all concrete algorithm classes should be derived.
Definition Algorithm.h:76
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Kernel::Logger & g_log
Definition Algorithm.h:422
@ Save
to specify a file to write to, the file may or may not exist
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:145
void error(const std::string &msg)
Logs at error level.
Definition Logger.cpp:108
void information(const std::string &msg)
Logs at information level.
Definition Logger.cpp:136
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:24
STL namespace.
@ Input
An input workspace.
Definition Property.h:53