Mantid
Loading...
Searching...
No Matches
SumOverlappingTubes.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
20#include "MantidHistogramData/LinearGenerator.h"
28#include "MantidKernel/Unit.h"
31
32#include <boost/math/special_functions/round.hpp>
33
34namespace Mantid::Algorithms {
35
36DECLARE_ALGORITHM(SumOverlappingTubes)
37
38using namespace API;
39using namespace Geometry;
40using namespace HistogramData;
41using namespace DataObjects;
42using namespace Kernel;
43
45 declareProperty(std::make_unique<ArrayProperty<std::string>>("InputWorkspaces", std::make_shared<ADSValidator>()),
46 "The names of the input workspaces as a list. You may also "
47 "group workspaces using the GUI or [[GroupWorkspaces]], and "
48 "specify the name of the group instead.");
49 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
50 "Name of the output workspace.");
51 std::vector<std::string> outputTypes{"2DTubes", "2D", "1D"};
52 declareProperty("OutputType", "2D", std::make_shared<StringListValidator>(outputTypes),
53 "Whether to have the output in raw 2D, with no "
54 "Debye-Scherrer cone correction, 2D or 1D.");
55 declareProperty(std::make_unique<ArrayProperty<double>>("ScatteringAngleBinning", "0.05",
56 std::make_shared<RebinParamsValidator>(), Direction::Input),
57 "A comma separated list of the first scattering angle, the scattering "
58 "angle step size and the final scattering angle. Optionally this can "
59 "also be a single number, which is the angle step size. In this case, "
60 "the boundary of binning will be determined by minimum and maximum "
61 "scattering angle present in the workspaces.");
62 declareProperty(std::make_unique<PropertyWithValue<bool>>("CropNegativeScatteringAngles", false, Direction::Input),
63 "If true the negative scattering angles are cropped (ignored).");
65 std::make_unique<ArrayProperty<double>>("HeightAxis", std::make_shared<RebinParamsValidator>(true, true)),
66 "A comma separated list of the first y value, the y value step size and "
67 "the final y value. This can also be a single number, which "
68 "is the y value step size. In this case, the boundary of binning will "
69 "be determined by minimum and maximum y values present in the "
70 "workspaces. This can also be two numbers to give the range desired.");
71 declareProperty(std::make_unique<PropertyWithValue<bool>>("Normalise", true, Direction::Input),
72 "If true normalise to the number of entries added for a particular "
73 "scattering angle. ");
74 declareProperty(std::make_unique<PropertyWithValue<bool>>("MirrorScatteringAngles", false, Direction::Input),
75 "A flag to mirror the signed 2thetas. ");
76 declareProperty(std::make_unique<PropertyWithValue<bool>>("SplitCounts", false, Direction::Input),
77 "A flag to split the counts between adjacent bins");
78 auto toleranceValidator = std::make_shared<BoundedValidator<double>>(0.0, 0.0);
79 toleranceValidator->clearUpper();
80 declareProperty("ScatteringAngleTolerance", 0.0, toleranceValidator,
81 "The relative tolerance for the scattering angles before the "
82 "counts are split.");
83 setPropertySettings("ScatteringAngleTolerance",
84 std::make_unique<Kernel::EnabledWhenProperty>("SplitCounts", IS_NOT_DEFAULT));
85}
86
89
90 m_progress = std::make_unique<Progress>(this, 0.0, 1.0, m_workspaceList.size());
91
92 // we need histogram data with m_numPoints bins
93 HistogramData::BinEdges x(m_numPoints + 1, LinearGenerator(m_startScatteringAngle, m_stepScatteringAngle));
94
95 MatrixWorkspace_sptr outputWS = create<Workspace2D>(m_numHistograms, x);
96 outputWS->setDistribution(false);
97 outputWS->setSharedRun(m_workspaceList.front()->sharedRun());
98
99 auto newAxis = std::make_unique<NumericAxis>(m_heightAxis);
100 newAxis->setUnit("Label");
101 auto yLabelUnit = std::dynamic_pointer_cast<Kernel::Units::Label>(newAxis->unit());
102 yLabelUnit->setLabel("Height", "m");
103 newAxis->unit() = yLabelUnit;
104 outputWS->replaceAxis(1, std::move(newAxis));
105
106 outputWS->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("Label");
107 Unit_sptr xUnit = outputWS->getAxis(0)->unit();
108 std::shared_ptr<Units::Label> xLabel = std::dynamic_pointer_cast<Units::Label>(xUnit);
109 xLabel->setLabel("Scattering Angle", "degrees");
110
111 const auto normalisation = performBinning(outputWS);
112
113 if (getProperty("Normalise")) {
115 for (int64_t ii = 0; ii < static_cast<int64_t>(m_numPoints); ++ii)
116 for (size_t j = 0; j < m_numHistograms; ++j) {
117 // Avoid spurious normalisation for low counting cells
118 const auto i = static_cast<size_t>(ii);
119 if (normalisation[j][i] < 1e-15)
120 continue;
121 outputWS->mutableY(j)[i] /= normalisation[j][i];
122 outputWS->mutableE(j)[i] /= normalisation[j][i];
123 }
124 }
125
126 setProperty("OutputWorkspace", outputWS);
127}
128
130
131 // This is flag for flipping the sign of 2theta
132 m_mirrorDetectors = getProperty("MirrorScatteringAngles") ? -1 : 1;
133
134 const std::vector<std::string> inputWorkspaces = getProperty("InputWorkspaces");
135 auto workspaces = RunCombinationHelper::unWrapGroups(inputWorkspaces);
136 RunCombinationHelper combHelper;
137 m_workspaceList = combHelper.validateInputWorkspaces(workspaces, g_log);
138 m_outputType = getPropertyValue("OutputType");
139 const auto &instrument = m_workspaceList.front()->getInstrument();
140 std::string componentName = "";
141 auto componentNameParam = instrument->getStringParameter("detector_for_height_axis");
142 if (!componentNameParam.empty())
143 componentName = componentNameParam[0];
145 getHeightAxis(componentName);
146}
147
150 m_endScatteringAngle = -180.0;
151
152 // Loop to check minimum and maximum extents for workspace
153 for (auto &ws : m_workspaceList) {
154 const auto &specInfo = ws->spectrumInfo();
155 for (size_t i = 0; i < specInfo.size(); ++i) {
156 if (specInfo.isMonitor(i) || specInfo.isMasked(i))
157 continue;
158 const auto &pos = specInfo.position(i);
159 const double theta = atan2(pos.X(), pos.Z()) * m_mirrorDetectors * 180 / M_PI;
162 }
163 }
164
165 const std::vector<double> scatteringBinning = getProperty("ScatteringAngleBinning");
166 if (scatteringBinning.size() == 1) {
167 m_stepScatteringAngle = scatteringBinning[0];
168 // Extend the boundaries by half of the step size
171 } else if (scatteringBinning.size() == 3) {
172 if (scatteringBinning[0] > m_startScatteringAngle || scatteringBinning[2] < m_endScatteringAngle)
173 g_log.warning() << "Some detectors outside of scattering angle range.\n";
174 m_startScatteringAngle = scatteringBinning[0];
175 m_stepScatteringAngle = scatteringBinning[1];
176 m_endScatteringAngle = scatteringBinning[2];
177 }
178
179 if (getProperty("CropNegativeScatteringAngles")) {
180 if (m_endScatteringAngle < 0) {
181 throw std::runtime_error("No positive scattering angle range");
182 }
183 if (m_startScatteringAngle < 0) {
185 }
186 }
187
189 static_cast<size_t>(std::floor((m_endScatteringAngle - m_startScatteringAngle) / m_stepScatteringAngle));
190 g_log.information() << "Number of bins:" << m_numPoints << std::endl;
191 g_log.information() << "Scattering angle binning:" << m_startScatteringAngle << ", " << m_stepScatteringAngle << ", "
192 << m_endScatteringAngle << "\n";
193
195 throw std::runtime_error("Wrong scattering angle range, check your binning/data");
196 }
197}
198
199void SumOverlappingTubes::getHeightAxis(const std::string &componentName) {
200 std::vector<double> heightBinning = getProperty("HeightAxis");
201 m_heightAxis.clear();
202 if (componentName.length() == 0 && heightBinning.empty())
203 throw std::runtime_error("No detector_for_height_axis parameter for this "
204 "instrument. Please enter a value for the "
205 "HeightAxis parameter.");
206 if ((componentName.length() > 0 && heightBinning.empty()) || (m_outputType != "1D" && heightBinning.size() == 2)) {
207 // Try to get the component. It should be a tube with pixels in the
208 // y-direction, the height bins are then taken as the detector positions.
209 const auto &componentInfo = m_workspaceList.front()->componentInfo();
210 const auto componentIndex = componentInfo.indexOfAny(componentName);
211 const auto &detsInSubtree = componentInfo.detectorsInSubtree(componentIndex);
212 for (const auto detIndex : detsInSubtree) {
213 const auto posY = componentInfo.position({detIndex, 0}).Y();
214 if (heightBinning.size() == 2 && (posY < heightBinning[0] || posY > heightBinning[1]))
215 continue;
216 m_heightAxis.emplace_back(posY);
217 }
218 } else {
219 if (heightBinning.size() != 3) {
220 if (heightBinning.size() == 2 && m_outputType == "1D") {
221 m_heightAxis.emplace_back(heightBinning[0]);
222 m_heightAxis.emplace_back(heightBinning[1]);
223 } else
224 throw std::runtime_error("Height binning must have start, step and end "
225 "values (except for 1D option).");
226 } else if (m_outputType == "1D") {
227 m_heightAxis.emplace_back(heightBinning[0]);
228 m_heightAxis.emplace_back(heightBinning[2]);
229 } else {
230 double height = heightBinning[0];
231 while (height < heightBinning[2]) {
232 m_heightAxis.emplace_back(height);
233 height += heightBinning[1];
234 }
235 }
236 }
237
238 m_startHeight = *min_element(m_heightAxis.begin(), m_heightAxis.end());
239 m_endHeight = *max_element(m_heightAxis.begin(), m_heightAxis.end());
240
241 if (m_outputType == "1D")
242 m_heightAxis = {(m_heightAxis.front() + m_heightAxis.back()) * 0.5};
243
245
246 g_log.information() << "Number of histograms in output workspace:" << m_numHistograms << ".\n";
247 g_log.information() << "Height axis:" << m_heightAxis[0] << " to " << m_heightAxis[m_numHistograms - 1] << " with "
248 << m_heightAxis.size() << " entries.\n";
249}
250
251std::vector<std::vector<double>> SumOverlappingTubes::performBinning(MatrixWorkspace_sptr &outputWS) {
252 const double scatteringAngleTolerance = getProperty("ScatteringAngleTolerance");
253 const bool splitCounts = getProperty("SplitCounts");
254
255 std::vector<std::vector<double>> normalisation(m_numHistograms, std::vector<double>(m_numPoints, 0.0));
256
257 // loop over all workspaces
258 for (auto &ws : m_workspaceList) {
259 m_progress->report("Processing workspace " + std::string(ws->getName()));
260 // loop over spectra
261 const auto &specInfo = ws->spectrumInfo();
262 PARALLEL_FOR_IF(Kernel::threadSafe(*ws, *outputWS))
263 for (int i = 0; i < static_cast<int>(specInfo.size()); ++i) {
265 if (specInfo.isMonitor(i) || specInfo.isMasked(i))
266 continue;
267
268 const auto &pos = specInfo.position(i);
269 const auto height = pos.Y();
270
271 const double tolerance = 1e-6;
272 if (height < m_startHeight - tolerance || height > m_endHeight + tolerance)
273 continue;
274
275 size_t heightIndex;
276 try {
278 } catch (std::out_of_range &) {
279 continue;
280 }
281
282 double angle;
283 if (m_outputType == "2DTubes")
284 angle = atan2(pos.X(), pos.Z());
285 else
286 angle = specInfo.signedTwoTheta(i);
287 angle *= m_mirrorDetectors * 180.0 / M_PI;
288
289 const auto angleIndex = static_cast<int>(std::floor((angle - m_startScatteringAngle) / m_stepScatteringAngle));
290
291 // point is out of range, a warning should have been generated already for
292 // the theta index
293 if (angleIndex < 0 || angleIndex >= int(m_numPoints))
294 continue;
295
296 const double deltaAngle = distanceFromAngle(angleIndex, angle);
297 const auto counts = ws->histogram(i).y()[0];
298 const auto error = ws->histogram(i).e()[0];
299
300 PARALLEL_CRITICAL(Histogramming2ThetaVsHeight) {
301 auto &yData = outputWS->mutableY(heightIndex);
302 auto &eData = outputWS->mutableE(heightIndex);
303 // counts are split between bins if outside this tolerance
304 if (splitCounts && deltaAngle > m_stepScatteringAngle * scatteringAngleTolerance) {
305 int angleIndexNeighbor;
306 if (distanceFromAngle(angleIndex - 1, angle) < distanceFromAngle(angleIndex + 1, angle))
307 angleIndexNeighbor = angleIndex - 1;
308 else
309 angleIndexNeighbor = angleIndex + 1;
310
311 double deltaAngleNeighbor = distanceFromAngle(angleIndexNeighbor, angle);
312
313 const auto scalingFactor = deltaAngleNeighbor / m_stepScatteringAngle;
314 const auto newError = error * scalingFactor;
315 yData[angleIndex] += counts * scalingFactor;
316 eData[angleIndex] = sqrt(eData[angleIndex] * eData[angleIndex] + newError * newError);
317
318 normalisation[heightIndex][angleIndex] += (deltaAngleNeighbor / m_stepScatteringAngle);
319
320 if (angleIndexNeighbor >= 0 && angleIndexNeighbor < int(m_numPoints)) {
321 const auto scalingFactorNeighbor = deltaAngle / m_stepScatteringAngle;
322 const auto newErrorNeighbor = error * scalingFactorNeighbor;
323 yData[angleIndexNeighbor] += counts * scalingFactorNeighbor;
324 eData[angleIndexNeighbor] =
325 sqrt(eData[angleIndexNeighbor] * eData[angleIndexNeighbor] + newErrorNeighbor * newErrorNeighbor);
326 normalisation[heightIndex][angleIndexNeighbor] += (deltaAngle / m_stepScatteringAngle);
327 }
328 } else {
329 yData[angleIndex] += counts;
330 eData[angleIndex] = sqrt(eData[angleIndex] * eData[angleIndex] + error * error);
331 normalisation[heightIndex][angleIndex]++;
332 }
333 }
335 }
337 }
338
339 return normalisation;
340}
341
342double SumOverlappingTubes::distanceFromAngle(const int angleIndex, const double angle) const {
343 return fabs(m_startScatteringAngle + double(angleIndex) * m_stepScatteringAngle - angle);
344}
345
346} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double height
Definition: GetAllEi.cpp:155
double error
Definition: IndexPeaks.cpp:133
#define fabs(x)
Definition: Matrix.cpp:22
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_FOR_NO_WSP_CHECK()
#define PARALLEL_CRITICAL(name)
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
double tolerance
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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
A property class for workspaces.
static std::vector< std::string > unWrapGroups(const std::vector< std::string > &)
Flattens the list of group workspaces (if any) into list of workspaces.
std::list< API::MatrixWorkspace_sptr > validateInputWorkspaces(const std::vector< std::string > &inputWorkspaces, Kernel::Logger &g_log)
Checks that the input workspace all exist, that they are the same size, have the same units and the s...
void exec() override
Virtual method - must be overridden by concrete algorithm.
double distanceFromAngle(const int angleIndex, const double angle) const
std::list< API::MatrixWorkspace_sptr > m_workspaceList
void init() override
Virtual method - must be overridden by concrete algorithm.
std::vector< std::vector< double > > performBinning(API::MatrixWorkspace_sptr &outputWS)
std::unique_ptr< API::Progress > m_progress
holds the sign flipper for 2theta
void getHeightAxis(const std::string &componentName)
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
size_t MANTID_KERNEL_DLL indexOfValueFromCenters(const std::vector< double > &bin_centers, const double value)
Gets the bin of a value from a vector of bin centers and throws exception if out of range.
std::shared_ptr< Unit > Unit_sptr
Shared pointer to the Unit base class.
Definition: Unit.h:229
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54