SaveGDA.cpp 9.76 KB
Newer Older
1
2
3
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4
5
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
// SPDX - License - Identifier: GPL - 3.0 +
Joseph Ramsay's avatar
Joseph Ramsay committed
7
8
9
10
11
12
13
#include "MantidDataHandling/SaveGDA.h"

#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidAPI/WorkspaceProperty.h"
14
#include "MantidKernel/ArrayProperty.h"
Joseph Ramsay's avatar
Joseph Ramsay committed
15
#include "MantidKernel/Unit.h"
16

Joseph Ramsay's avatar
Joseph Ramsay committed
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <boost/optional.hpp>

#include <cmath>
#include <fstream>
#include <iomanip>
#include <sstream>

namespace Mantid {
namespace DataHandling {

using namespace API;

namespace { // helper functions

31
32
const int POINTS_PER_LINE = 4;

Joseph Ramsay's avatar
Joseph Ramsay committed
33
double mean(const std::vector<double> &values) {
34
  return std::accumulate(values.begin(), values.end(), 0.0) / static_cast<double>(values.size());
Joseph Ramsay's avatar
Joseph Ramsay committed
35
36
}

37
38
39
// Compute the mean resolution of the x axis of the input workspace
// Resolution is calculated as the difference between adjacent pairs of values,
// normalised by the second of the two
Joseph Ramsay's avatar
Joseph Ramsay committed
40
41
42
43
double computeAverageDeltaTByT(const HistogramData::HistogramX &tValues) {
  std::vector<double> deltaTByT;
  deltaTByT.reserve(tValues.size() - 1);

44
45
  std::adjacent_difference(tValues.begin(), tValues.end(), std::back_inserter(deltaTByT),
                           [](const double previous, const double current) { return (previous - current) / current; });
Joseph Ramsay's avatar
Joseph Ramsay committed
46
47
48
49
50
  // First element is just first element of tValues, so remove it
  deltaTByT.erase(deltaTByT.begin());
  return mean(deltaTByT);
}

51
std::string generateBankHeader(int bank, int minT, size_t numberBins, double deltaTByT) {
Joseph Ramsay's avatar
Joseph Ramsay committed
52
  std::stringstream stream;
53
  const auto numberLines = static_cast<size_t>(std::ceil(static_cast<double>(numberBins) / POINTS_PER_LINE));
54

55
56
  stream << std::setprecision(2) << "BANK " << bank << " " << numberBins << "  " << numberLines << " RALF  " << minT
         << "  96  " << minT << " " << deltaTByT << " ALT";
Joseph Ramsay's avatar
Joseph Ramsay committed
57
58
59
  return stream.str();
}

60
boost::optional<std::vector<std::string>> getParamLinesFromGSASFile(const std::string &paramsFilename) {
61
62
  // ICONS signifies that a line contains TOF to D conversion factors
  const static std::string paramLineDelimiter = "ICONS";
63
64
65
66
  std::ifstream paramsFile;
  paramsFile.open(paramsFilename);

  if (paramsFile.is_open()) {
67
    std::vector<std::string> paramLines;
68
69
    std::string line;
    while (std::getline(paramsFile, line)) {
70
      if (line.find(paramLineDelimiter) != std::string::npos) {
71
72
73
        paramLines.emplace_back(line);
      }
    }
74
    return paramLines;
75
76
77
78
79
  } else {
    return boost::none;
  }
}

Joseph Ramsay's avatar
Joseph Ramsay committed
80
81
82
83
} // anonymous namespace

DECLARE_ALGORITHM(SaveGDA)

84
SaveGDA::CalibrationParams::CalibrationParams(const double _difa, const double _difc, const double _tzero)
85
86
    : difa(_difa), difc(_difc), tzero(_tzero) {}

Joseph Ramsay's avatar
Joseph Ramsay committed
87
88
89
90
91
92
93
94
const std::string SaveGDA::name() const { return "SaveGDA"; }

const std::string SaveGDA::summary() const {
  return "Save a group of focused banks to the MAUD three-column GDA format";
}

int SaveGDA::version() const { return 1; }

95
const std::vector<std::string> SaveGDA::seeAlso() const { return {"SaveBankScatteringAngles", "AlignDetectors"}; }
Joseph Ramsay's avatar
Joseph Ramsay committed
96

97
const std::string SaveGDA::category() const { return "DataHandling\\Text;Diffraction\\DataHandling"; }
Joseph Ramsay's avatar
Joseph Ramsay committed
98

99
const std::string SaveGDA::PROP_OUTPUT_FILENAME = "OutputFilename";
Joseph Ramsay's avatar
Joseph Ramsay committed
100
101
102

const std::string SaveGDA::PROP_INPUT_WS = "InputWorkspace";

103
104
105
106
const std::string SaveGDA::PROP_PARAMS_FILENAME = "GSASParamFile";

const std::string SaveGDA::PROP_GROUPING_SCHEME = "GroupingScheme";

Joseph Ramsay's avatar
Joseph Ramsay committed
107
void SaveGDA::init() {
108
  declareProperty(std::make_unique<WorkspaceProperty<WorkspaceGroup>>(PROP_INPUT_WS, "", Kernel::Direction::Input),
Joseph Ramsay's avatar
Joseph Ramsay committed
109
110
111
112
                  "A GroupWorkspace where every sub-workspace is a "
                  "single-spectra focused run corresponding to a particular "
                  "bank");

113
  const static std::vector<std::string> outExts{".gda"};
114
  declareProperty(std::make_unique<FileProperty>(PROP_OUTPUT_FILENAME, "", FileProperty::Save, outExts),
Joseph Ramsay's avatar
Joseph Ramsay committed
115
                  "The name of the file to save to");
116

117
118
119
120
121
122
123
124
  const static std::vector<std::string> paramsExts{".ipf", ".prm", ".parm", ".iprm"};
  declareProperty(std::make_unique<FileProperty>(PROP_PARAMS_FILENAME, "", FileProperty::Load, paramsExts),
                  "GSAS calibration file containing conversion factors from D to TOF");

  declareProperty(std::make_unique<Kernel::ArrayProperty<int>>(PROP_GROUPING_SCHEME),
                  "An array of bank IDs, where the value at element i is the "
                  "ID of the bank in " +
                      PROP_PARAMS_FILENAME + " to associate spectrum i with");
Joseph Ramsay's avatar
Joseph Ramsay committed
125
126
127
}

void SaveGDA::exec() {
128
  const std::string filename = getProperty(PROP_OUTPUT_FILENAME);
Joseph Ramsay's avatar
Joseph Ramsay committed
129
130
131
132
133
134
135
136
137
  std::ofstream outFile(filename.c_str());

  if (!outFile) {
    throw Kernel::Exception::FileError("Unable to create file: ", filename);
  }

  outFile << std::fixed << std::setprecision(0) << std::setfill(' ');

  const API::WorkspaceGroup_sptr inputWS = getProperty(PROP_INPUT_WS);
138
139
140
  const auto calibParams = parseParamsFile();
  const std::vector<int> groupingScheme = getProperty(PROP_GROUPING_SCHEME);

141
  for (int i = 0; i < inputWS->getNumberOfEntries(); ++i) {
Joseph Ramsay's avatar
Joseph Ramsay committed
142
    const auto ws = inputWS->getItem(i);
143
    const auto matrixWS = std::dynamic_pointer_cast<MatrixWorkspace>(ws);
Joseph Ramsay's avatar
Joseph Ramsay committed
144

145
    auto &x = matrixWS->dataX(0);
146
147
    const size_t bankIndex(groupingScheme[i] - 1);
    if (bankIndex >= calibParams.size()) {
148
      throw Kernel::Exception::IndexError(bankIndex, calibParams.size(), "Bank number out of range");
149
150
    }
    const auto &bankCalibParams = calibParams[bankIndex];
151

152
153
    // For historic reasons, TOF is scaled by 32 in MAUD
    const static double tofScale = 32;
154
    std::vector<double> tofScaled;
155
156
    tofScaled.reserve(x.size());
    Kernel::Units::dSpacing dSpacingUnit;
157
158
    std::vector<double> yunused;
    dSpacingUnit.toTOF(x, yunused, 0., Kernel::DeltaEMode::Elastic,
159
160
161
162
163
164
                       {{Kernel::UnitParams::difa, bankCalibParams.difa},
                        {Kernel::UnitParams::difa, bankCalibParams.difc},
                        {Kernel::UnitParams::tzero, bankCalibParams.tzero}});
    std::transform(x.begin(), x.end(), std::back_inserter(tofScaled),
                   [](const double tofVal) { return tofVal * tofScale; });

165
    const auto averageDeltaTByT = computeAverageDeltaTByT(tofScaled);
Joseph Ramsay's avatar
Joseph Ramsay committed
166
167
168

    const auto &intensity = matrixWS->y(0);
    const auto &error = matrixWS->e(0);
169
    const auto numPoints = std::min({tofScaled.size(), intensity.size(), error.size()});
Joseph Ramsay's avatar
Joseph Ramsay committed
170

171
    const auto header =
172
        generateBankHeader(i + 1, static_cast<int>(std::round(tofScaled[0])), numPoints, averageDeltaTByT);
173
174
175

    outFile << std::left << std::setw(80) << header << '\n' << std::right;

Joseph Ramsay's avatar
Joseph Ramsay committed
176
    for (size_t j = 0; j < numPoints; ++j) {
177
      outFile << std::setw(8) << tofScaled[j] << std::setw(7) << intensity[j] * 1000 << std::setw(5) << error[j] * 1000;
Joseph Ramsay's avatar
Joseph Ramsay committed
178

179
      if (j % POINTS_PER_LINE == POINTS_PER_LINE - 1) {
Joseph Ramsay's avatar
Joseph Ramsay committed
180
181
        // new line every 4 points
        outFile << '\n';
182
      } else if (j == numPoints - 1) {
Joseph Ramsay's avatar
Joseph Ramsay committed
183
        // make sure line is 80 characters long
184
        outFile << std::string(80 - (i % POINTS_PER_LINE + 1) * 20, ' ') << '\n';
Joseph Ramsay's avatar
Joseph Ramsay committed
185
186
187
188
189
190
191
192
193
194
195
      }
    }
  }
}

std::map<std::string, std::string> SaveGDA::validateInputs() {
  std::map<std::string, std::string> issues;
  boost::optional<std::string> inputWSIssue;

  const API::WorkspaceGroup_sptr inputWS = getProperty(PROP_INPUT_WS);
  for (const auto &ws : *inputWS) {
196
    const auto matrixWS = std::dynamic_pointer_cast<MatrixWorkspace>(ws);
Joseph Ramsay's avatar
Joseph Ramsay committed
197
198
199
200
201
202
    if (matrixWS) {
      if (matrixWS->getNumberHistograms() != 1) {
        inputWSIssue = "The workspace " + matrixWS->getName() +
                       " has the wrong number of histograms. It "
                       "should contain data for a single focused "
                       "spectra";
203
      } else if (matrixWS->getAxis(0)->unit()->unitID() != "dSpacing") {
Joseph Ramsay's avatar
Joseph Ramsay committed
204
205
206
        inputWSIssue = "The workspace " + matrixWS->getName() +
                       " has incorrect units. SaveGDA "
                       "expects input workspaces with "
207
                       "units of D-spacing";
Joseph Ramsay's avatar
Joseph Ramsay committed
208
209
      }
    } else { // not matrixWS
210
      inputWSIssue = "The workspace " + ws->getName() + " is of the wrong type. It should be a MatrixWorkspace";
Joseph Ramsay's avatar
Joseph Ramsay committed
211
212
213
214
215
    }
  }
  if (inputWSIssue) {
    issues[PROP_INPUT_WS] = *inputWSIssue;
  }
216
217
218

  const std::vector<int> groupingScheme = getProperty(PROP_GROUPING_SCHEME);
  const auto numSpectraInGroupingScheme = groupingScheme.size();
219
  const auto numSpectraInWS = static_cast<size_t>(inputWS->getNumberOfEntries());
220
  if (numSpectraInGroupingScheme != numSpectraInWS) {
221
222
223
224
    issues[PROP_GROUPING_SCHEME] = "The grouping scheme must contain one entry for every focused spectrum "
                                   "in the input workspace. " +
                                   PROP_GROUPING_SCHEME + " has " + std::to_string(numSpectraInGroupingScheme) +
                                   " entries whereas " + PROP_INPUT_WS + " has " + std::to_string(numSpectraInWS);
225
226
  }

Joseph Ramsay's avatar
Joseph Ramsay committed
227
228
229
  return issues;
}

230
231
232
233
234
std::vector<SaveGDA::CalibrationParams> SaveGDA::parseParamsFile() const {
  const std::string paramsFilename = getProperty(PROP_PARAMS_FILENAME);
  const auto paramLines = getParamLinesFromGSASFile(paramsFilename);
  if (!paramLines) {
    g_log.error(strerror(errno));
235
    throw Kernel::Exception::FileError("Could not read GSAS parameter file", paramsFilename);
236
237
238
239
  }
  std::vector<CalibrationParams> calibParams;
  for (const auto &paramLine : *paramLines) {
    std::vector<std::string> lineItems;
240
241
    boost::algorithm::split(lineItems, paramLine, boost::is_any_of("\t "), boost::token_compress_on);
    calibParams.emplace_back(std::stod(lineItems[3]), std::stod(lineItems[4]), std::stod(lineItems[5]));
242
243
244
245
  }
  return calibParams;
}

LamarMoore's avatar
LamarMoore committed
246
247
} // namespace DataHandling
} // namespace Mantid