/
LoadSPE.cpp
246 lines (215 loc) · 7.33 KB
/
LoadSPE.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
//---------------------------------------------------
// Includes
//---------------------------------------------------
#include "MantidDataHandling/LoadSPE.h"
#include "MantidDataHandling/SaveSPE.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/BinEdgeAxis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/Histogram1D.h"
#include "MantidKernel/UnitFactory.h"
#include <cstdio>
#include <limits>
#include <fstream>
namespace Mantid {
namespace DataHandling {
using namespace Kernel;
using namespace API;
DECLARE_FILELOADER_ALGORITHM(LoadSPE)
/**
* Return the confidence with with this algorithm can load the file
* @param descriptor A descriptor for the file
* @returns An integer specifying the confidence level. 0 indicates it will not
* be used
*/
int LoadSPE::confidence(Kernel::FileDescriptor &descriptor) const {
if (!descriptor.isAscii())
return 0;
auto &file = descriptor.data();
std::string fileline;
// First line - expected to be a line 2 columns which is histogram & bin
// numbers
std::getline(file, fileline);
std::istringstream is(fileline);
unsigned int dummy(0);
is >> dummy >> dummy;
if (is.fail()) {
return 0; // Couldn't read 2 numbers so fail
}
// Trying to read another should produce eof
is >> dummy;
if (!is.eof())
return 0;
// Next line should be comment line: "### Phi Grid" or "### Q Grid"
std::getline(file, fileline);
if (fileline.find("Phi Grid") != std::string::npos ||
fileline.find("Q Grid") != std::string::npos) {
return 80;
} else
return 0;
}
//---------------------------------------------------
// Private member functions
//---------------------------------------------------
/**
* Initialise the algorithm
*/
void LoadSPE::init() {
declareProperty(
make_unique<FileProperty>("Filename", "", FileProperty::Load, ".spe"),
"The name of the SPE file to load.");
declareProperty(make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
Direction::Output),
"The name to use for the output workspace");
}
/**
* Execute the algorithm
*/
void LoadSPE::exec() {
// Retrieve filename and try to open the file
m_filename = getPropertyValue("Filename");
FILE *speFile;
speFile = fopen(m_filename.c_str(), "r");
if (!speFile) {
g_log.error("Failed to open file: " + m_filename);
throw Exception::FileError("Failed to open file: ", m_filename);
}
// The first two numbers are the number of histograms and the number of bins
size_t nhist = 0, nbins = 0;
unsigned int nhistTemp = 0, nbinsTemp = 0;
int retval = fscanf(speFile, "%8u%8u\n", &nhistTemp, &nbinsTemp);
if (retval != 2)
reportFormatError("Header line");
// Cast from temp values to size_t values
nhist = static_cast<size_t>(nhistTemp);
nbins = static_cast<size_t>(nbinsTemp);
// Next line should be comment line: "### Phi Grid" or "### Q Grid"
char comment[100];
fgets(comment, 100, speFile);
if (comment[0] != '#')
reportFormatError(std::string(comment));
// Create the axis that will hold the phi values
auto *phiAxis = new BinEdgeAxis(nhist + 1);
// Look at previously read comment field to see what unit vertical axis should
// have
if (comment[4] == 'Q' || comment[4] == 'q') {
phiAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
} else {
phiAxis->unit() = boost::make_shared<Units::Phi>();
}
// Read in phi grid
for (size_t i = 0; i <= nhist; ++i) {
double phi;
retval = fscanf(speFile, "%10le", &phi);
if (retval != 1) {
std::stringstream ss;
ss << "Reading phi value" << i;
reportFormatError(ss.str());
}
phiAxis->setValue(i, phi);
}
// Read to EOL
fgets(comment, 100, speFile);
// Next line should be comment line: "### Energy Grid"
fgets(comment, 100, speFile);
if (comment[0] != '#')
reportFormatError(std::string(comment));
// Now the X bin boundaries
MantidVecPtr XValues;
MantidVec &X = XValues.access();
X.resize(nbins + 1);
for (size_t i = 0; i <= nbins; ++i) {
retval = fscanf(speFile, "%10le", &X[i]);
if (retval != 1) {
std::stringstream ss;
ss << "Reading energy value" << i;
reportFormatError(ss.str());
}
}
// Read to EOL
fgets(comment, 100, speFile);
// Now create the output workspace
MatrixWorkspace_sptr workspace = WorkspaceFactory::Instance().create(
"Workspace2D", nhist, nbins + 1, nbins);
workspace->getAxis(0)->unit() = UnitFactory::Instance().create("DeltaE");
workspace->setDistribution(true); // It should be a distribution
workspace->setYUnitLabel("S(Phi,Energy)");
// Replace the default spectrum axis with the phi values one
workspace->replaceAxis(1, phiAxis);
// Now read in the data spectrum-by-spectrum
Progress progress(this, 0, 1, nhist);
for (size_t j = 0; j < nhist; ++j) {
// Set the common X vector
workspace->setX(j, XValues);
// Read in the Y & E data
readHistogram(speFile, workspace, j);
progress.report();
}
// Close the file
fclose(speFile);
// Set the output workspace property
setProperty("OutputWorkspace", workspace);
}
/** Reads in the data corresponding to a single spectrum
* @param speFile :: The file handle
* @param workspace :: The output workspace
* @param index :: The index of the current spectrum
*/
void LoadSPE::readHistogram(FILE *speFile, API::MatrixWorkspace_sptr workspace,
size_t index) {
// First, there should be a comment line
char comment[100];
fgets(comment, 100, speFile);
if (comment[0] != '#')
reportFormatError(std::string(comment));
// Then it's the Y values
MantidVec &Y = workspace->dataY(index);
const size_t nbins = workspace->blocksize();
int retval;
for (size_t i = 0; i < nbins; ++i) {
retval = fscanf(speFile, "%10le", &Y[i]);
// g_log.error() << Y[i] << std::endl;
if (retval != 1) {
std::stringstream ss;
ss << "Reading data value" << i << " of histogram " << index;
reportFormatError(ss.str());
}
// -10^30 is the flag for not a number used in SPE files (from
// www.mantidproject.org/images/3/3d/Spe_file_format.pdf)
if (Y[i] == SaveSPE::MASK_FLAG) {
Y[i] = std::numeric_limits<double>::quiet_NaN();
}
}
// Read to EOL
fgets(comment, 100, speFile);
// Another comment line
fgets(comment, 100, speFile);
if (comment[0] != '#')
reportFormatError(std::string(comment));
// And then the error values
MantidVec &E = workspace->dataE(index);
for (size_t i = 0; i < nbins; ++i) {
retval = fscanf(speFile, "%10le", &E[i]);
if (retval != 1) {
std::stringstream ss;
ss << "Reading error value" << i << " of histogram " << index;
reportFormatError(ss.str());
}
}
// Read to EOL
fgets(comment, 100, speFile);
return;
}
/** Called if the file is not formatted as expected
* @param what :: A string describing where the problem occurred
* @throw Mantid::Kernel::Exception::FileError terminating the algorithm
*/
void LoadSPE::reportFormatError(const std::string &what) {
g_log.error("Unexpected formatting in file " + m_filename + " : " + what);
throw Exception::FileError("Unexpected formatting in file: ", m_filename);
}
} // namespace DataHandling
} // namespace Mantid