Mantid
Loading...
Searching...
No Matches
SaveGDA.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
9#include "MantidAPI/Axis.h"
15#include "MantidKernel/Unit.h"
16
17#include <boost/optional.hpp>
18
19#include <cmath>
20#include <fstream>
21#include <iomanip>
22#include <sstream>
23
24namespace Mantid::DataHandling {
25
26using namespace API;
27
28namespace { // helper functions
29
30const int POINTS_PER_LINE = 4;
31
32double mean(const std::vector<double> &values) {
33 return std::accumulate(values.begin(), values.end(), 0.0) / static_cast<double>(values.size());
34}
35
36// Compute the mean resolution of the x axis of the input workspace
37// Resolution is calculated as the difference between adjacent pairs of values,
38// normalised by the second of the two
39double computeAverageDeltaTByT(const HistogramData::HistogramX &tValues) {
40 std::vector<double> deltaTByT;
41 deltaTByT.reserve(tValues.size() - 1);
42
43 std::adjacent_difference(tValues.begin(), tValues.end(), std::back_inserter(deltaTByT),
44 [](const double previous, const double current) { return (previous - current) / current; });
45 // First element is just first element of tValues, so remove it
46 deltaTByT.erase(deltaTByT.begin());
47 return mean(deltaTByT);
48}
49
50std::string generateBankHeader(int bank, int minT, size_t numberBins, double deltaTByT) {
51 std::stringstream stream;
52 const auto numberLines = static_cast<size_t>(std::ceil(static_cast<double>(numberBins) / POINTS_PER_LINE));
53
54 stream << std::setprecision(2) << "BANK " << bank << " " << numberBins << " " << numberLines << " RALF " << minT
55 << " 96 " << minT << " " << deltaTByT << " ALT";
56 return stream.str();
57}
58
59boost::optional<std::vector<std::string>> getParamLinesFromGSASFile(const std::string &paramsFilename) {
60 // ICONS signifies that a line contains TOF to D conversion factors
61 const static std::string paramLineDelimiter = "ICONS";
62 std::ifstream paramsFile;
63 paramsFile.open(paramsFilename);
64
65 if (paramsFile.is_open()) {
66 std::vector<std::string> paramLines;
67 std::string line;
68 while (std::getline(paramsFile, line)) {
69 if (line.find(paramLineDelimiter) != std::string::npos) {
70 paramLines.emplace_back(line);
71 }
72 }
73 return paramLines;
74 } else {
75 return boost::none;
76 }
77}
78
79} // anonymous namespace
80
81DECLARE_ALGORITHM(SaveGDA)
82
83SaveGDA::CalibrationParams::CalibrationParams(const double _difc, const double _difa, const double _tzero)
84 : difa(_difa), difc(_difc), tzero(_tzero) {}
85
86const std::string SaveGDA::name() const { return "SaveGDA"; }
87
88const std::string SaveGDA::summary() const {
89 return "Save a group of focused banks to the MAUD three-column GDA format";
90}
91
92int SaveGDA::version() const { return 1; }
93
94const std::vector<std::string> SaveGDA::seeAlso() const { return {"SaveBankScatteringAngles", "AlignDetectors"}; }
95
96const std::string SaveGDA::category() const { return "DataHandling\\Text;Diffraction\\DataHandling"; }
97
98const std::string SaveGDA::PROP_OUTPUT_FILENAME = "OutputFilename";
99
100const std::string SaveGDA::PROP_INPUT_WS = "InputWorkspace";
101
102const std::string SaveGDA::PROP_PARAMS_FILENAME = "GSASParamFile";
103
104const std::string SaveGDA::PROP_GROUPING_SCHEME = "GroupingScheme";
105
108 "A GroupWorkspace where every sub-workspace is a "
109 "single-spectra focused run corresponding to a particular "
110 "bank");
111
112 const static std::vector<std::string> outExts{".gda"};
113 declareProperty(std::make_unique<FileProperty>(PROP_OUTPUT_FILENAME, "", FileProperty::Save, outExts),
114 "The name of the file to save to");
115
116 const static std::vector<std::string> paramsExts{".ipf", ".prm", ".parm", ".iprm"};
117 declareProperty(std::make_unique<FileProperty>(PROP_PARAMS_FILENAME, "", FileProperty::Load, paramsExts),
118 "GSAS calibration file containing conversion factors from D to TOF");
119
121 "An array of bank IDs, where the value at element i is the "
122 "ID of the bank in " +
123 PROP_PARAMS_FILENAME + " to associate spectrum i with");
124}
125
127 const std::string filename = getProperty(PROP_OUTPUT_FILENAME);
128 std::ofstream outFile(filename.c_str());
129
130 if (!outFile) {
131 throw Kernel::Exception::FileError("Unable to create file: ", filename);
132 }
133
134 outFile << std::fixed << std::setprecision(0) << std::setfill(' ');
135
137 const auto calibParams = parseParamsFile();
138 const std::vector<int> groupingScheme = getProperty(PROP_GROUPING_SCHEME);
139
140 for (int i = 0; i < inputWS->getNumberOfEntries(); ++i) {
141 const auto ws = inputWS->getItem(i);
142 const auto matrixWS = std::dynamic_pointer_cast<MatrixWorkspace>(ws);
143
144 auto x = matrixWS->dataX(0);
145 const size_t bankIndex(groupingScheme[i] - 1);
146 if (bankIndex >= calibParams.size()) {
147 throw Kernel::Exception::IndexError(bankIndex, calibParams.size(), "Bank number out of range");
148 }
149 const auto &bankCalibParams = calibParams[bankIndex];
150
151 // For historic reasons, TOF is scaled by 32 in MAUD
152 const static double tofScale = 32;
153 std::vector<double> tofScaled;
154 tofScaled.reserve(x.size());
155 Kernel::Units::dSpacing dSpacingUnit;
156 std::vector<double> yunused;
157 dSpacingUnit.toTOF(x, yunused, 0., Kernel::DeltaEMode::Elastic,
158 {{Kernel::UnitParams::difa, bankCalibParams.difa},
159 {Kernel::UnitParams::difc, bankCalibParams.difc},
160 {Kernel::UnitParams::tzero, bankCalibParams.tzero}});
161 std::transform(x.begin(), x.end(), std::back_inserter(tofScaled),
162 [](const double tofVal) { return tofVal * tofScale; });
163 const auto averageDeltaTByT = computeAverageDeltaTByT(tofScaled);
164
165 const auto &intensity = matrixWS->y(0);
166 const auto &error = matrixWS->e(0);
167 const auto numPoints = std::min({tofScaled.size(), intensity.size(), error.size()});
168
169 const auto header =
170 generateBankHeader(i + 1, static_cast<int>(std::round(tofScaled[0])), numPoints, averageDeltaTByT);
171
172 outFile << std::left << std::setw(80) << header << '\n' << std::right;
173
174 for (size_t j = 0; j < numPoints; ++j) {
175 outFile << std::setw(8) << tofScaled[j] << std::setw(7) << intensity[j] * 1000 << std::setw(5) << error[j] * 1000;
176
177 if (j % POINTS_PER_LINE == POINTS_PER_LINE - 1) {
178 // new line every 4 points
179 outFile << '\n';
180 } else if (j == numPoints - 1) {
181 // make sure line is 80 characters long
182 outFile << std::string(80 - (i % POINTS_PER_LINE + 1) * 20, ' ') << '\n';
183 }
184 }
185 }
186}
187
188std::map<std::string, std::string> SaveGDA::validateInputs() {
189 std::map<std::string, std::string> issues;
190 boost::optional<std::string> inputWSIssue;
191
193 for (const auto &ws : *inputWS) {
194 const auto matrixWS = std::dynamic_pointer_cast<MatrixWorkspace>(ws);
195 if (matrixWS) {
196 if (matrixWS->getNumberHistograms() != 1) {
197 inputWSIssue = "The workspace " + matrixWS->getName() +
198 " has the wrong number of histograms. It "
199 "should contain data for a single focused "
200 "spectra";
201 } else if (matrixWS->getAxis(0)->unit()->unitID() != "dSpacing") {
202 inputWSIssue = "The workspace " + matrixWS->getName() +
203 " has incorrect units. SaveGDA "
204 "expects input workspaces with "
205 "units of D-spacing";
206 }
207 } else { // not matrixWS
208 inputWSIssue = "The workspace " + ws->getName() + " is of the wrong type. It should be a MatrixWorkspace";
209 }
210 }
211 if (inputWSIssue) {
212 issues[PROP_INPUT_WS] = *inputWSIssue;
213 }
214
215 const std::vector<int> groupingScheme = getProperty(PROP_GROUPING_SCHEME);
216 const auto numSpectraInGroupingScheme = groupingScheme.size();
217 const auto numSpectraInWS = static_cast<size_t>(inputWS->getNumberOfEntries());
218 if (numSpectraInGroupingScheme != numSpectraInWS) {
219 issues[PROP_GROUPING_SCHEME] = "The grouping scheme must contain one entry for every focused spectrum "
220 "in the input workspace. " +
221 PROP_GROUPING_SCHEME + " has " + std::to_string(numSpectraInGroupingScheme) +
222 " entries whereas " + PROP_INPUT_WS + " has " + std::to_string(numSpectraInWS);
223 }
224
225 return issues;
226}
227
228std::vector<SaveGDA::CalibrationParams> SaveGDA::parseParamsFile() const {
229 const std::string paramsFilename = getProperty(PROP_PARAMS_FILENAME);
230 const auto paramLines = getParamLinesFromGSASFile(paramsFilename);
231 if (!paramLines) {
232 g_log.error(strerror(errno));
233 throw Kernel::Exception::FileError("Could not read GSAS parameter file", paramsFilename);
234 }
235 std::vector<CalibrationParams> calibParams;
236 for (const auto &paramLine : *paramLines) {
237 std::vector<std::string> lineItems;
238 boost::algorithm::split(lineItems, paramLine, boost::is_any_of("\t "), boost::token_compress_on);
239 calibParams.emplace_back(std::stod(lineItems[3]), std::stod(lineItems[4]), std::stod(lineItems[5]));
240 }
241 return calibParams;
242}
243
244} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double error
Definition: IndexPeaks.cpp:133
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
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
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
A property class for workspaces.
const std::string name() const override
function to return a name of the algorithm, must be overridden in all algorithms
Definition: SaveGDA.cpp:86
int version() const override
function to return a version of the algorithm, must be overridden in all algorithms
Definition: SaveGDA.cpp:92
static const std::string PROP_GROUPING_SCHEME
Definition: SaveGDA.h:43
static const std::string PROP_PARAMS_FILENAME
Definition: SaveGDA.h:42
const std::string category() const override
function to return a category of the algorithm.
Definition: SaveGDA.cpp:96
std::map< std::string, std::string > validateInputs() override
Perform validation of ALL the input properties of the algorithm.
Definition: SaveGDA.cpp:188
std::vector< CalibrationParams > parseParamsFile() const
Definition: SaveGDA.cpp:228
const std::string summary() const override
function returns a summary message that will be displayed in the default GUI, and in the help.
Definition: SaveGDA.cpp:88
static const std::string PROP_INPUT_WS
Definition: SaveGDA.h:41
void exec() override
Virtual method - must be overridden by concrete algorithm.
Definition: SaveGDA.cpp:126
const std::vector< std::string > seeAlso() const override
Function to return all of the seeAlso (these are not validated) algorithms related to this algorithm....
Definition: SaveGDA.cpp:94
static const std::string PROP_OUTPUT_FILENAME
Definition: SaveGDA.h:40
void init() override
Virtual method - must be overridden by concrete algorithm.
Definition: SaveGDA.cpp:106
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Records the filename and the description of failure.
Definition: Exception.h:98
Exception for index errors.
Definition: Exception.h:284
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void toTOF(std::vector< double > &xdata, std::vector< double > &ydata, const double &_l1, const int &_emode, std::initializer_list< std::pair< const UnitParams, double > > params)
Convert from the concrete unit to time-of-flight.
Definition: Unit.cpp:147
d-Spacing in Angstrom
Definition: Unit.h:383
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53