-
Notifications
You must be signed in to change notification settings - Fork 122
/
LoadILLPolarizedDiffraction.cpp
480 lines (425 loc) · 18 KB
/
LoadILLPolarizedDiffraction.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
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidDataHandling/LoadILLPolarizedDiffraction.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidGeometry/Instrument/ComponentHelper.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/OptionalBool.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/UnitLabelTypes.h"
#include "MantidKernel/V3D.h"
#include "MantidKernel/VisibleWhenProperty.h"
#include <Poco/Path.h>
namespace Mantid {
namespace DataHandling {
using namespace API;
using namespace Geometry;
using namespace Kernel;
using namespace NeXus;
using Types::Core::DateAndTime;
namespace {
// This defines the number of detector banks in D7
constexpr size_t D7_NUMBER_BANKS = 3;
// This defines the number of physical pixels in D7
constexpr size_t D7_NUMBER_PIXELS = 132;
// This defines the number of pixels per bank in D7
constexpr size_t D7_NUMBER_PIXELS_BANK = 44;
// This defines the number of monitors in the instrument. If there are cases
// where this is no longer one this decleration should be moved.
constexpr size_t NUMBER_MONITORS = 2;
// This defines Time Of Flight measurement mode switch value
constexpr size_t TOF_MODE_ON = 1;
} // namespace
// Register the algorithm into the AlgorithmFactory
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadILLPolarizedDiffraction)
/// Returns confidence. @see IFileLoader::confidence
int LoadILLPolarizedDiffraction::confidence(NexusDescriptor &descriptor) const {
// fields existent only at the ILL Diffraction
if (descriptor.pathExists("/entry0/D7")) {
return 80;
} else {
return 0;
}
}
/// Algorithms name for identification. @see Algorithm::name
const std::string LoadILLPolarizedDiffraction::name() const {
return "LoadILLPolarizedDiffraction";
}
/// Algorithm's version for identification. @see Algorithm::version
int LoadILLPolarizedDiffraction::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string LoadILLPolarizedDiffraction::category() const {
return "DataHandling\\Nexus;ILL\\Diffraction";
}
/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
const std::string LoadILLPolarizedDiffraction::summary() const {
return "Loads ILL D7 instrument polarized diffraction nexus files.";
}
/**
* Constructor
*/
LoadILLPolarizedDiffraction::LoadILLPolarizedDiffraction()
: IFileLoader<NexusDescriptor>() {}
/**
* Initialize the algorithm's properties.
*/
void LoadILLPolarizedDiffraction::init() {
declareProperty(std::make_unique<FileProperty>("Filename", "",
FileProperty::Load, ".nxs"),
"File path of the data file to load");
declareProperty(std::make_unique<WorkspaceProperty<API::WorkspaceGroup>>(
"OutputWorkspace", "", Direction::Output),
"The output workspace.");
const std::vector<std::string> positionCalibrationOptions{"None", "Nexus",
"YIGFile"};
declareProperty(
"PositionCalibration", "None",
std::make_shared<StringListValidator>(positionCalibrationOptions),
"Select the type of pixel position calibration. If None, the pixel "
"positions are read from IDF file. If Nexus, the positions are read from "
"Nexus file. If YIGFile, then the calibration twotheta data is loaded "
"from a user-defined calibration file.");
declareProperty(std::make_unique<FileProperty>(
"YIGFilename", "", FileProperty::OptionalLoad, ".xml"),
"File path of the YIG calibration data file to load");
setPropertySettings("YIGFilename",
std::make_unique<Kernel::EnabledWhenProperty>(
"PositionCalibration", IS_EQUAL_TO, "YIGFile"));
declareProperty("ConvertToScatteringAngle", false,
"Convert the bin edges to scattering angle",
Direction::Input);
declareProperty("TransposeMonochromatic", false,
"Transpose the 2D workspace with monochromatic data",
Direction::Input);
}
std::map<std::string, std::string>
LoadILLPolarizedDiffraction::validateInputs() {
std::map<std::string, std::string> issues;
if (getPropertyValue("PositionCalibration") == "YIGFile" &&
getPropertyValue("YIGFilename") == "") {
issues["PositionCalibration"] =
"YIG-based position calibration of detectors requested but "
"the file was not provided.";
}
return issues;
}
/**
* Executes the algorithm.
*/
void LoadILLPolarizedDiffraction::exec() {
Progress progress(this, 0, 1, 2);
m_fileName = getPropertyValue("Filename");
m_outputWorkspaceGroup = std::make_shared<API::WorkspaceGroup>();
progress.report("Loading the detector polarization analysis data");
loadData();
progress.report("Loading the metadata");
loadMetaData();
setProperty("OutputWorkspace", m_outputWorkspaceGroup);
}
/**
* Loads the polarized detector data, sets up workspaces and labels
* according to the measurement type and data dimensions
*/
void LoadILLPolarizedDiffraction::loadData() {
// open the root entry
NXRoot dataRoot(m_fileName);
// read each entry
for (auto entryNumber = 0;
entryNumber < static_cast<int>((dataRoot.groups().size()));
entryNumber++) {
NXEntry entry = dataRoot.openEntry("entry" + std::to_string(entryNumber));
m_instName = entry.getString("D7/name");
std::string start_time = entry.getString("start_time");
start_time = m_loadHelper.dateTimeInIsoFormat(start_time);
// prepare axes for data
std::vector<double> axis = prepareAxes(entry);
// init the workspace with proper number of histograms and number of
// channels
auto workspace = initStaticWorkspace(entry);
// load data from file
std::string dataName = "data/Detector_data";
NXUInt data = entry.openNXDataSet<unsigned int>(dataName);
data.load();
// Assign detector counts
PARALLEL_FOR_IF(Kernel::threadSafe(*workspace))
for (auto pixel_no = 0; pixel_no < static_cast<int>(D7_NUMBER_PIXELS);
++pixel_no) {
auto &spectrum = workspace->mutableY(pixel_no);
auto &errors = workspace->mutableE(pixel_no);
for (auto channel_no = 0;
channel_no < static_cast<int>(m_numberOfChannels); ++channel_no) {
unsigned int counts = data(pixel_no, 0, channel_no);
spectrum[channel_no] = counts;
errors[channel_no] = std::sqrt(counts);
}
workspace->mutableX(pixel_no) = axis;
}
// load and assign monitor data
for (auto monitor_no = static_cast<int>(D7_NUMBER_PIXELS);
monitor_no < static_cast<int>(D7_NUMBER_PIXELS + NUMBER_MONITORS);
++monitor_no) {
NXUInt monitorData = entry.openNXDataSet<unsigned int>(
"monitor" +
std::to_string(monitor_no + 1 - static_cast<int>(D7_NUMBER_PIXELS)) +
"/data");
monitorData.load();
auto &spectrum = workspace->mutableY(monitor_no);
auto &errors = workspace->mutableE(monitor_no);
for (auto channel_no = 0;
channel_no < static_cast<int>(m_numberOfChannels); channel_no++) {
unsigned int counts = monitorData(0, 0, channel_no);
spectrum[channel_no] = counts;
errors[channel_no] = std::sqrt(counts);
}
workspace->mutableX(monitor_no) = axis;
}
// load the instrument
loadInstrument(workspace);
// rotate detectors to their position during measurement
moveTwoTheta(entry, workspace);
// convert the spectrum axis to scattering angle
if (getProperty("ConvertToScatteringAngle")) {
workspace = convertSpectrumAxis(workspace);
}
// transpose monochromatic data distribution
if (getProperty("TransposeMonochromatic") &&
m_acquisitionMode != TOF_MODE_ON) {
workspace = transposeMonochromatic(workspace);
}
// adds the current entry workspace to the output group
m_outputWorkspaceGroup->addWorkspace(workspace);
entry.close();
}
dataRoot.close();
}
/**
* Dumps the metadata from the file for each entry separately
*/
void LoadILLPolarizedDiffraction::loadMetaData() {
// Open NeXus file
NXhandle nxHandle;
NXstatus nxStat = NXopen(m_fileName.c_str(), NXACC_READ, &nxHandle);
if (nxStat != NX_ERROR) {
for (auto workspaceId = 0;
workspaceId < m_outputWorkspaceGroup->getNumberOfEntries();
++workspaceId) {
MatrixWorkspace_sptr workspace =
std::static_pointer_cast<API::MatrixWorkspace>(
m_outputWorkspaceGroup->getItem(workspaceId));
auto const entryName = std::string("entry" + std::to_string(workspaceId));
m_loadHelper.addNexusFieldsToWsRun(nxHandle, workspace->mutableRun(),
entryName);
}
NXclose(&nxHandle);
}
}
/**
* Initializes the output workspace based on the resolved instrument.
* If there are multiple entries in the file and the current entry
* is not the first one, the returned workspace is a clone
* of the workspace from the first entry
* @param entry : entry linked with the returned workspace
* @return : workspace with the correct data dimensions
*/
API::MatrixWorkspace_sptr
LoadILLPolarizedDiffraction::initStaticWorkspace(const NXEntry &entry) {
const size_t nSpectra = D7_NUMBER_PIXELS + NUMBER_MONITORS;
API::MatrixWorkspace_sptr workspace = WorkspaceFactory::Instance().create(
"Workspace2D", nSpectra, m_numberOfChannels + 1, m_numberOfChannels);
// Set x axis units
NXInt acquisitionMode = entry.openNXInt("acquisition_mode");
acquisitionMode.load();
m_acquisitionMode = acquisitionMode[0];
if (m_acquisitionMode == TOF_MODE_ON) {
auto lblUnit = std::static_pointer_cast<Kernel::Units::Label>(
UnitFactory::Instance().create("Label"));
lblUnit->setLabel("Time", Units::Symbol::Microsecond);
workspace->getAxis(0)->unit() = lblUnit;
} else {
auto lblUnit = std::static_pointer_cast<Kernel::Units::Label>(
UnitFactory::Instance().create("Label"));
lblUnit->setLabel("Wavelength", Units::Symbol::Angstrom);
workspace->getAxis(0)->unit() = lblUnit;
}
// Set y axis unit
workspace->setYUnit("Counts");
// check the polarization direction and set the workspace title
std::string polDirection = entry.getString("D7/POL/actual_state");
std::string flipperState = entry.getString("D7/POL/actual_stateB1B2");
workspace->setTitle(polDirection.substr(0, 1) + "_" + flipperState);
return workspace;
}
/**
* Runs LoadInstrument as child to link the instrument to workspace
* @param workspace : workspace with data from the first entry
*/
void LoadILLPolarizedDiffraction::loadInstrument(
API::MatrixWorkspace_sptr workspace) {
IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");
loadInst->setPropertyValue("Filename", m_instName + "_Definition.xml");
loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", workspace);
loadInst->setProperty("RewriteSpectraMap", OptionalBool(true));
loadInst->execute();
}
/**
* Loads twotheta for each detector pixel from the file
* @param workspace : workspace with loaded instrument
* @param entry : entry from which the pixel 2theta positions will be read
* @param bankId : bank ID for which 2theta positions will be read
* @return : vector of pixel 2theta positions in the chosen bank
*/
std::vector<double> LoadILLPolarizedDiffraction::loadTwoThetaDetectors(
const API::MatrixWorkspace_sptr workspace, const NXEntry &entry,
const int bankId) {
std::vector<double> twoTheta(static_cast<int>(D7_NUMBER_PIXELS_BANK));
if (getPropertyValue("PositionCalibration") == "Nexus") {
NXFloat twoThetaPixels = entry.openNXFloat(
"D7/Detector/bank" + std::to_string(bankId) + "_offset");
twoThetaPixels.load();
float *twoThetaDataStart = twoThetaPixels();
float *twoThetaDataEnd = twoThetaDataStart + D7_NUMBER_PIXELS_BANK;
twoTheta.assign(twoThetaDataStart, twoThetaDataEnd);
} else {
IAlgorithm_sptr loadInst = createChildAlgorithm("LoadParameterFile");
loadInst->setPropertyValue("Filename", getPropertyValue("YIGFilename"));
loadInst->setProperty("Workspace", workspace);
loadInst->execute();
auto instrumentMap = workspace->instrumentParameters();
Instrument_const_sptr instrument = workspace->getInstrument();
IComponent_const_sptr currentBank = instrument->getComponentByName(
std::string("bank" + std::to_string(bankId)));
for (auto pixel_no = 0; pixel_no < static_cast<int>(D7_NUMBER_PIXELS_BANK);
pixel_no++) {
std::string twoThetaRead = currentBank->getParameterAsString(
std::string("twoTheta_pixel_" + std::to_string(pixel_no + 1)));
twoTheta[pixel_no] = std::stod(twoThetaRead);
}
}
return twoTheta;
}
/**
* Rotates each pixel to its corresponding 2theta read from the file
* @param entry : entry from which the 2theta positions will be read
* @param workspace : workspace containing the instrument being moved
*/
void LoadILLPolarizedDiffraction::moveTwoTheta(
const NXEntry &entry, API::MatrixWorkspace_sptr workspace) {
Instrument_const_sptr instrument = workspace->getInstrument();
auto &componentInfo = workspace->mutableComponentInfo();
for (auto bank_no = 0; bank_no < static_cast<int>(D7_NUMBER_BANKS);
++bank_no) {
NXFloat twoThetaBank = entry.openNXFloat(
"D7/2theta/actual_bank" +
std::to_string(bank_no + 2)); // detector bank IDs start at 2
twoThetaBank.load();
if (getPropertyValue("PositionCalibration") == "None") {
Quat rotation(-twoThetaBank[0], V3D(0, 1, 0));
IComponent_const_sptr currentBank = instrument->getComponentByName(
std::string("bank" + std::to_string(bank_no + 2)));
const auto componentIndex =
componentInfo.indexOf(currentBank->getComponentID());
componentInfo.setRotation(componentIndex, rotation);
} else {
std::vector<double> twoThetaPixels =
loadTwoThetaDetectors(workspace, entry, bank_no + 2);
for (auto pixel_no = 0;
pixel_no < static_cast<int>(D7_NUMBER_PIXELS_BANK); ++pixel_no) {
auto const pixelIndex =
bank_no * static_cast<int>(D7_NUMBER_PIXELS_BANK) + pixel_no;
auto const pixel = componentInfo.componentID(pixelIndex);
V3D position = pixel->getPos();
double radius, theta, phi;
position.getSpherical(radius, theta, phi);
position.spherical(radius, twoThetaBank[0] - twoThetaPixels[pixel_no],
phi);
componentInfo.setPosition(pixelIndex, position);
}
}
}
}
/**
* Prepares values for bin edges depending of measurement type
* @param entry : entry from which the number of channels and measurement type
* will be read
* @return : returns vector with bin edges
*/
std::vector<double>
LoadILLPolarizedDiffraction::prepareAxes(const NXEntry &entry) {
// check the mode of measurement and prepare axes for data
std::vector<double> axes;
NXInt acquisitionMode = entry.openNXInt("acquisition_mode");
acquisitionMode.load();
if (acquisitionMode[0] == TOF_MODE_ON) {
NXFloat timeOfFlightInfo = entry.openNXFloat("D7/Detector/time_of_flight");
timeOfFlightInfo.load();
auto channelWidth = static_cast<double>(timeOfFlightInfo[0]);
m_numberOfChannels = size_t(timeOfFlightInfo[1]);
auto tofDelay = timeOfFlightInfo[2];
for (auto channel_no = 0;
channel_no <= static_cast<int>(m_numberOfChannels); channel_no++) {
axes.push_back(static_cast<double>(tofDelay + channel_no * channelWidth));
}
} else {
m_numberOfChannels = 1;
NXFloat wavelength = entry.openNXFloat("D7/monochromator/wavelength");
wavelength.load();
axes.push_back(static_cast<double>(wavelength[0] * 0.99));
axes.push_back(static_cast<double>(wavelength[0] * 1.01));
}
return axes;
}
/**
* Converts the spectrum axis to scattering angle
* @param workspace : workspace to change the
*/
API::MatrixWorkspace_sptr LoadILLPolarizedDiffraction::convertSpectrumAxis(
API::MatrixWorkspace_sptr workspace) {
IAlgorithm_sptr convertSpectrumAxis =
createChildAlgorithm("ConvertSpectrumAxis");
convertSpectrumAxis->initialize();
convertSpectrumAxis->setProperty("InputWorkspace", workspace);
convertSpectrumAxis->setProperty("OutputWorkspace", "__unused_for_child");
convertSpectrumAxis->setProperty("Target", "SignedTheta");
convertSpectrumAxis->setProperty("EMode", "Direct");
convertSpectrumAxis->setProperty("OrderAxis", false);
convertSpectrumAxis->execute();
workspace = convertSpectrumAxis->getProperty("OutputWorkspace");
IAlgorithm_sptr changeSign = createChildAlgorithm("ConvertAxisByFormula");
changeSign->initialize();
changeSign->setProperty("InputWorkspace", workspace);
changeSign->setProperty("OutputWorkspace", "__unused_for_child");
changeSign->setProperty("Axis", "Y");
changeSign->setProperty("Formula", "-y");
changeSign->execute();
return changeSign->getProperty("OutputWorkspace");
}
/**
* Transposes given 2D workspace with monochromatic data
* @param workspace : workspace to be transposed
*/
API::MatrixWorkspace_sptr LoadILLPolarizedDiffraction::transposeMonochromatic(
API::MatrixWorkspace_sptr workspace) {
IAlgorithm_sptr transpose = createChildAlgorithm("Transpose");
transpose->initialize();
transpose->setProperty("InputWorkspace", workspace);
transpose->setProperty("OutputWorkspace", "__unused_for_child");
transpose->execute();
return transpose->getProperty("OutputWorkspace");
}
} // namespace DataHandling
} // namespace Mantid