-
Notifications
You must be signed in to change notification settings - Fork 122
/
EQSANSQ2D.cpp
205 lines (178 loc) · 9.64 KB
/
EQSANSQ2D.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
// 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 "MantidWorkflowAlgorithms/EQSANSQ2D.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidGeometry/Instrument.h"
#include "MantidWorkflowAlgorithms/EQSANSInstrument.h"
#include "Poco/NumberFormatter.h"
namespace Mantid::WorkflowAlgorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(EQSANSQ2D)
using namespace Kernel;
using namespace API;
using namespace Geometry;
void EQSANSQ2D::init() {
auto wsValidator = std::make_shared<WorkspaceUnitValidator>("Wavelength");
declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator),
"Workspace to calculate I(qx,qy) from");
declareProperty("OutputWorkspace", "", Direction::Input);
declareProperty("NumberOfBins", 100, "Number of bins in each dimension of the 2D output", Kernel::Direction::Input);
declareProperty("IQxQyLogBinning", false, "I(qx,qy) log binning when binning is not specified.",
Kernel::Direction::Input);
declareProperty("OutputMessage", "", Direction::Output);
}
/// Returns the value of a run property from a given workspace
/// @param inputWS :: input workspace
/// @param pname :: name of the property to retrieve
double getRunProperty(const MatrixWorkspace_sptr &inputWS, const std::string &pname) {
return inputWS->run().getPropertyValueAsType<double>(pname);
}
/// Execute algorithm
void EQSANSQ2D::exec() {
Progress progress(this, 0.0, 1.0, 3);
progress.report("Setting up I(qx,Qy) calculation");
MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
const int nbins = getProperty("NumberOfBins");
// If the OutputWorkspace property was not given, use the
// name of the input workspace as the base name for the output
std::string outputWSName = getPropertyValue("OutputWorkspace");
if (outputWSName.empty()) {
outputWSName = inputWS->getName();
}
// Determine whether we need frame skipping or not by checking the chopper
// speed
bool frame_skipping = false;
const auto &run = inputWS->run();
if (run.hasProperty("is_frame_skipping")) {
auto prop = run.getProperty("is_frame_skipping");
const auto &typeInfo = *(prop->type_info());
if (typeInfo == typeid(int64_t)) {
frame_skipping = (run.getPropertyValueAsType<int64_t>("is_frame_skipping") == 1);
} else if (typeInfo == typeid(int32_t)) {
frame_skipping = (run.getPropertyValueAsType<int32_t>("is_frame_skipping") == 1);
} else
g_log.warning() << "Unknown property type for is_frame_skipping\n";
}
// Get run properties necessary to calculate the input parameters to Qxy
// Wavelength bandwidths
double wavelength_min = 0.0;
if (inputWS->run().hasProperty("wavelength_min"))
wavelength_min = getRunProperty(inputWS, "wavelength_min");
else if (inputWS->dataX(0).size() > 1)
wavelength_min = (inputWS->dataX(1)[0] + inputWS->dataX(1)[1]) / 2.0;
else if (inputWS->dataX(0).size() == 1)
wavelength_min = inputWS->dataX(1)[0];
else {
g_log.error("Can't determine the minimum wavelength for the input workspace.");
throw std::invalid_argument("Can't determine the minimum wavelength for the input workspace.");
}
double qmax = 0;
if (inputWS->run().hasProperty("qmax")) {
qmax = getRunProperty(inputWS, "qmax");
g_log.debug() << "Using Qmax from run properties = " << qmax << std::endl;
} else {
// This is pointing to the correct parameter in the workspace,
// which has been changed to the right distance in EQSANSLoad.cpp
const double sample_detector_distance = getRunProperty(inputWS, "sample_detector_distance");
const double nx_pixels = inputWS->getInstrument()->getNumberParameter("number-of-x-pixels")[0];
const double ny_pixels = inputWS->getInstrument()->getNumberParameter("number-of-y-pixels")[0];
const double pixel_size_x = inputWS->getInstrument()->getNumberParameter("x-pixel-size")[0];
const double pixel_size_y = inputWS->getInstrument()->getNumberParameter("y-pixel-size")[0];
const double beam_ctr_x = getRunProperty(inputWS, "beam_center_x");
const double beam_ctr_y = getRunProperty(inputWS, "beam_center_y");
double dxmax = pixel_size_x * std::max(beam_ctr_x, nx_pixels - beam_ctr_x);
double dymax = pixel_size_y * std::max(beam_ctr_y, ny_pixels - beam_ctr_y);
double maxdist = std::max(dxmax, dymax);
// This uses the correct parameter in the workspace, which has been
// changed to the right distance in EQSANSLoad.cpp
qmax = 4 * M_PI / wavelength_min * std::sin(0.5 * std::atan(maxdist / sample_detector_distance));
g_log.debug() << "Using calculated Qmax = " << qmax << std::endl;
};
if (frame_skipping) {
// In frame skipping mode, treat each frame separately
const double wavelength_max = getRunProperty(inputWS, "wavelength_max");
const double wavelength_min_f2 = getRunProperty(inputWS, "wavelength_min_frame2");
const double wavelength_max_f2 = getRunProperty(inputWS, "wavelength_max_frame2");
// Frame 1
std::string params =
Poco::NumberFormatter::format(wavelength_min, 2) + ",0.1," + Poco::NumberFormatter::format(wavelength_max, 2);
auto rebinAlg = createChildAlgorithm("Rebin", 0.4, 0.5);
rebinAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWS);
rebinAlg->setPropertyValue("Params", params);
rebinAlg->setProperty("PreserveEvents", false);
rebinAlg->executeAsChildAlg();
const bool log_binning = getProperty("IQxQyLogBinning");
auto qxyAlg = createChildAlgorithm("Qxy", .5, .65);
qxyAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", rebinAlg->getProperty("OutputWorkspace"));
qxyAlg->setProperty<double>("MaxQxy", qmax);
qxyAlg->setProperty<double>("DeltaQ", qmax / nbins);
qxyAlg->setProperty<bool>("SolidAngleWeighting", false);
qxyAlg->setProperty<bool>("IQxQyLogBinning", log_binning);
qxyAlg->executeAsChildAlg();
MatrixWorkspace_sptr qxy_output = qxyAlg->getProperty("OutputWorkspace");
auto replaceAlg = createChildAlgorithm("ReplaceSpecialValues", .65, 0.7);
replaceAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", qxy_output);
replaceAlg->setProperty<double>("NaNValue", 0.0);
replaceAlg->setProperty<double>("NaNError", 0.0);
replaceAlg->executeAsChildAlg();
std::string outputWSName_frame = outputWSName + "_frame1_Iqxy";
declareProperty(
std::make_unique<WorkspaceProperty<>>("OutputWorkspaceFrame1", outputWSName_frame, Direction::Output));
MatrixWorkspace_sptr result = replaceAlg->getProperty("OutputWorkspace");
setProperty("OutputWorkspaceFrame1", result);
// Frame 2
params = Poco::NumberFormatter::format(wavelength_min_f2, 2) + ",0.1," +
Poco::NumberFormatter::format(wavelength_max_f2, 2);
rebinAlg = createChildAlgorithm("Rebin", 0.7, 0.8);
rebinAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWS);
rebinAlg->setPropertyValue("Params", params);
rebinAlg->setProperty("PreserveEvents", false);
rebinAlg->executeAsChildAlg();
qxyAlg = createChildAlgorithm("Qxy", .8, 0.95);
qxyAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", rebinAlg->getProperty("OutputWorkspace"));
qxyAlg->setProperty<double>("MaxQxy", qmax);
qxyAlg->setProperty<double>("DeltaQ", qmax / nbins);
qxyAlg->setProperty<bool>("SolidAngleWeighting", false);
qxyAlg->setProperty<bool>("IQxQyLogBinning", log_binning);
qxyAlg->executeAsChildAlg();
qxy_output = qxyAlg->getProperty("OutputWorkspace");
replaceAlg = createChildAlgorithm("ReplaceSpecialValues", .95, 1.0);
replaceAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", qxy_output);
replaceAlg->setProperty<double>("NaNValue", 0.0);
replaceAlg->setProperty<double>("NaNError", 0.0);
replaceAlg->executeAsChildAlg();
outputWSName_frame = outputWSName + "_frame2_Iqxy";
declareProperty(
std::make_unique<WorkspaceProperty<>>("OutputWorkspaceFrame2", outputWSName_frame, Direction::Output));
result = replaceAlg->getProperty("OutputWorkspace");
setProperty("OutputWorkspaceFrame2", result);
setProperty("OutputMessage", "I(Qx,Qy) computed for each frame");
} else {
// When not in frame skipping mode, simply run Qxy
const bool log_binning = getProperty("IQxQyLogBinning");
auto qxyAlg = createChildAlgorithm("Qxy", .3, 0.9);
qxyAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWS);
qxyAlg->setProperty<double>("MaxQxy", qmax);
qxyAlg->setProperty<double>("DeltaQ", qmax / nbins);
qxyAlg->setProperty<bool>("SolidAngleWeighting", false);
qxyAlg->setProperty<bool>("IQxQyLogBinning", log_binning);
qxyAlg->executeAsChildAlg();
MatrixWorkspace_sptr qxy_output = qxyAlg->getProperty("OutputWorkspace");
auto replaceAlg = createChildAlgorithm("ReplaceSpecialValues", .9, 1.0);
replaceAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", qxy_output);
replaceAlg->setProperty<double>("NaNValue", 0.0);
replaceAlg->setProperty<double>("NaNError", 0.0);
replaceAlg->executeAsChildAlg();
outputWSName += "_Iqxy";
declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspaceFrame1", outputWSName, Direction::Output));
MatrixWorkspace_sptr result = replaceAlg->getProperty("OutputWorkspace");
setProperty("OutputWorkspaceFrame1", result);
setProperty("OutputMessage", "I(Qx,Qy) computed for each frame");
}
}
} // namespace Mantid::WorkflowAlgorithms