Mantid
Loading...
Searching...
No Matches
SolidAngle.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 +
20
21#include <atomic>
22
23namespace Mantid::Algorithms {
24
25// Register with the algorithm factory
26DECLARE_ALGORITHM(SolidAngle)
27
28namespace SolidAngleMethods {
29static const std::string GENERIC_SHAPE = "GenericShape";
30static const std::string RECTANGLE = "Rectangle";
31static const std::string VERTICAL_TUBE = "VerticalTube";
32static const std::string HORIZONTAL_TUBE = "HorizontalTube";
33static const std::string VERTICAL_WING = "VerticalWing";
34static const std::string HORIZONTAL_WING = "HorizontalWing";
35} // namespace SolidAngleMethods
36
37using namespace Kernel;
38using namespace API;
39using namespace Geometry;
40using namespace SolidAngleMethods;
41
42namespace SolidAngleHelpers {
43
44constexpr double MM_TO_METERS = 1. / 1000.;
45
54 : m_detectorInfo(detectorInfo), m_samplePos(detectorInfo.samplePosition()) {}
55 double getAlpha(size_t index) const {
56 const auto sampleDetVec = m_detectorInfo.position(index) - m_samplePos;
57 auto inPlane = sampleDetVec;
58 project(inPlane);
59 return sampleDetVec.cosAngle(inPlane);
60 }
61 virtual void project(V3D &v) const = 0;
62 virtual ~AlphaAngleCalculator() = default;
63
64private:
67};
68
71 void project(V3D &v) const override { v.setY(0.0); }
72};
73
76 void project(V3D &v) const override { v.setX(0.0); }
77};
78
83 SolidAngleCalculator(const ComponentInfo &componentInfo, const DetectorInfo &detectorInfo, const std::string &method,
84 const double pixelArea)
85 : m_componentInfo(componentInfo), m_detectorInfo(detectorInfo), m_pixelArea(pixelArea),
86 m_samplePos(detectorInfo.samplePosition()), m_beamLine(m_samplePos - detectorInfo.sourcePosition()) {
87 if (method.find("Vertical") != std::string::npos) {
88 m_alphaAngleCalculator = std::make_unique<AlphaAngleVertical>(detectorInfo);
89 } else if (method.find("Horizontal") != std::string::npos) {
90 m_alphaAngleCalculator = std::make_unique<AlphaAngleHorizontal>(detectorInfo);
91 }
92 }
93 virtual double solidAngle(size_t index) const = 0;
94 virtual ~SolidAngleCalculator() = default;
95
96protected:
99 const double m_pixelArea;
102 std::unique_ptr<const AlphaAngleCalculator> m_alphaAngleCalculator;
103};
104
107 double solidAngle(size_t index) const override { return m_detectorInfo.detector(index).solidAngle(m_samplePos); }
108};
109
112 double solidAngle(size_t index) const override {
113 const V3D sampleDetVec = m_detectorInfo.position(index) - m_samplePos;
114 const double cosTheta = sampleDetVec.cosAngle(m_beamLine);
115 const double l2 = m_detectorInfo.l2(index);
116 const V3D scaleFactor = m_componentInfo.scaleFactor(index);
117 const double scaledPixelArea = m_pixelArea * scaleFactor[0] * scaleFactor[1];
118 return scaledPixelArea * cosTheta / (l2 * l2);
119 }
120};
121
122struct Tube : public SolidAngleCalculator {
124 double solidAngle(size_t index) const override {
125 const double cosAlpha = m_alphaAngleCalculator->getAlpha(index);
126 const double l2 = m_detectorInfo.l2(index);
127 const V3D scaleFactor = m_componentInfo.scaleFactor(index);
128 const double scaledPixelArea = m_pixelArea * scaleFactor[0] * scaleFactor[1];
129 return scaledPixelArea * cosAlpha / (l2 * l2);
130 }
131};
132
133struct Wing : public SolidAngleCalculator {
135 double solidAngle(size_t index) const override {
136 const V3D sampleDetVec = m_detectorInfo.position(index) - m_samplePos;
137 const double cosTheta = sampleDetVec.cosAngle(m_beamLine);
138 const double cosAlpha = m_alphaAngleCalculator->getAlpha(index);
139 const double l2 = m_detectorInfo.l2(index);
140 const V3D scaleFactor = m_componentInfo.scaleFactor(index);
141 const double scaledPixelArea = m_pixelArea * scaleFactor[0] * scaleFactor[1];
142 return scaledPixelArea * cosAlpha * cosAlpha * cosAlpha / (l2 * l2 * cosTheta * cosTheta);
143 }
144};
145
146} // namespace SolidAngleHelpers
147
150 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input,
151 std::make_shared<InstrumentValidator>()),
152 "This workspace is used to identify the instrument to use "
153 "and also which\n"
154 "spectra to create a solid angle for. If the Max and Min "
155 "spectra values are\n"
156 "not provided one solid angle will be created for each "
157 "spectra in the input\n"
158 "workspace");
159 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
160 "The name of the workspace to be created as the output of "
161 "the algorithm. A workspace of this name will be created "
162 "and stored in the Analysis Data Service.");
163
164 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
165 mustBePositive->setLower(0);
166 declareProperty("StartWorkspaceIndex", 0, mustBePositive,
167 "The index number of the first spectrum for which to find "
168 "the solid angle\n"
169 "(default: 0)");
170 declareProperty("EndWorkspaceIndex", EMPTY_INT(), mustBePositive,
171 "The index of the last spectrum whose solid angle is to be "
172 "found (default: the\n"
173 "last spectrum in the workspace)");
174
175 const std::vector<std::string> methods{GENERIC_SHAPE, RECTANGLE, VERTICAL_TUBE,
177 declareProperty("Method", GENERIC_SHAPE, std::make_shared<StringListValidator>(methods),
178 "Select the method to calculate the Solid Angle.");
179}
180
184 // Get the workspaces
185 MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
186 int m_MinSpec = getProperty("StartWorkspaceIndex");
187 int m_MaxSpec = getProperty("EndWorkspaceIndex");
188
189 const auto numberOfSpectra = static_cast<int>(inputWS->getNumberHistograms());
190
191 // Check 'StartSpectrum' is in range 0-numberOfSpectra
192 if (m_MinSpec > numberOfSpectra) {
193 g_log.warning("StartWorkspaceIndex out of range! Set to 0.");
194 m_MinSpec = 0;
195 }
196 if (isEmpty(m_MaxSpec))
197 m_MaxSpec = numberOfSpectra - 1;
198 if (m_MaxSpec > numberOfSpectra - 1 || m_MaxSpec < m_MinSpec) {
199 g_log.warning("EndWorkspaceIndex out of range! Set to max detector number");
200 m_MaxSpec = numberOfSpectra - 1;
201 }
202
203 MatrixWorkspace_sptr outputWS = WorkspaceFactory::Instance().create(inputWS, numberOfSpectra, 2, 1);
204 // The result of this will be a distribution
205 outputWS->setDistribution(true);
206 outputWS->setYUnit("");
207 outputWS->setYUnitLabel("Steradian");
208 setProperty("OutputWorkspace", outputWS);
209
210 const auto &spectrumInfo = inputWS->spectrumInfo();
211 const auto &detectorInfo = inputWS->detectorInfo();
212 const auto &componentInfo = inputWS->componentInfo();
213
214 // this is the pixel area that is supposed to be constant for the whole
215 // instrument this is used only if Method != GenericShape
216 double pixelArea = 0.;
217 const std::string method = getProperty("Method");
218
219 using namespace SolidAngleHelpers;
220 if (method != GENERIC_SHAPE) {
221 const auto instrument = inputWS->getInstrument();
222 if (instrument->hasParameter("x-pixel-size") && instrument->hasParameter("y-pixel-size")) {
223 const double pixelSizeX = instrument->getNumberParameter("x-pixel-size")[0] * MM_TO_METERS;
224 const double pixelSizeY = instrument->getNumberParameter("y-pixel-size")[0] * MM_TO_METERS;
225 pixelArea = pixelSizeX * pixelSizeY; // l2 is retrieved per pixel
226 } else {
227 // TODO: try to get the pixel sizes from bounding box
228 throw std::runtime_error("Missing necessary instrument parameters for non generic shape: "
229 "x-pixel-size and y-pixel-size [in mm].");
230 }
231 }
232
233 std::unique_ptr<SolidAngleCalculator> solidAngleCalculator;
234 if (method == GENERIC_SHAPE) {
235 solidAngleCalculator = std::make_unique<GenericShape>(componentInfo, detectorInfo, method, pixelArea);
236 } else if (method == RECTANGLE) {
237 solidAngleCalculator = std::make_unique<Rectangle>(componentInfo, detectorInfo, method, pixelArea);
238 } else if (method == VERTICAL_TUBE || method == HORIZONTAL_TUBE) {
239 solidAngleCalculator = std::make_unique<Tube>(componentInfo, detectorInfo, method, pixelArea);
240 } else if (method == VERTICAL_WING || method == HORIZONTAL_WING) {
241 solidAngleCalculator = std::make_unique<Wing>(componentInfo, detectorInfo, method, pixelArea);
242 }
243
244 std::atomic<size_t> failCount{0};
245 Progress prog(this, 0.0, 1.0, numberOfSpectra);
246 // Loop over the histograms (detector spectra)
247 PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS, *inputWS))
248 for (int j = m_MinSpec; j <= m_MaxSpec; ++j) {
250 initSpectrum(*inputWS, *outputWS, j);
251 if (spectrumInfo.hasDetectors(j)) {
252 double solidAngle = 0.0;
253 for (const auto detID : inputWS->getSpectrum(j).getDetectorIDs()) {
254 const auto index = detectorInfo.indexOf(detID);
255 if (!detectorInfo.isMasked(index) && !detectorInfo.isMonitor(index)) {
256 solidAngle += solidAngleCalculator->solidAngle(index);
257 }
258 }
259 outputWS->mutableY(j)[0] = solidAngle;
260 } else {
261 ++failCount;
262 }
263 prog.report();
265 } // loop over spectra
267 g_log.warning() << "min max total" << m_MinSpec << " " << m_MaxSpec << " " << numberOfSpectra;
268
269 auto &outputSpectrumInfo = outputWS->mutableSpectrumInfo();
270 // Loop over the histograms (detector spectra)
271 for (int j = 0; j < m_MinSpec; ++j) {
272 initSpectrum(*inputWS, *outputWS, j);
273 // SpectrumInfo::setMasked is NOT threadsafe.
274 outputSpectrumInfo.setMasked(j, true);
275 prog.report();
276 } // loop over spectra
277
278 // Loop over the histograms (detector spectra)
279 for (int j = m_MaxSpec + 1; j < numberOfSpectra; ++j) {
280 initSpectrum(*inputWS, *outputWS, j);
281 // SpectrumInfo::setMasked is NOT threadsafe.
282 outputSpectrumInfo.setMasked(j, true);
283 prog.report();
284 } // loop over spectra
285
286 if (failCount != 0) {
287 g_log.information() << "Unable to calculate solid angle for " << failCount
288 << " spectra. The solid angle will be set to zero for "
289 "those detectors.\n";
290 }
291}
292
297void SolidAngle::initSpectrum(const MatrixWorkspace &inputWS, MatrixWorkspace &outputWS, const size_t wsIndex) {
298 outputWS.mutableX(wsIndex)[0] = inputWS.x(wsIndex).front();
299 outputWS.mutableX(wsIndex)[1] = inputWS.x(wsIndex).back();
300 outputWS.mutableE(wsIndex) = 0.;
301 outputWS.mutableY(wsIndex) = 0.; // default value for not calculated
302}
303
304} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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_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.
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
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
Base MatrixWorkspace Abstract Class.
HistogramData::HistogramX & mutableX(const size_t index) &
const HistogramData::HistogramX & x(const size_t index) const
HistogramData::HistogramE & mutableE(const size_t index) &
HistogramData::HistogramY & mutableY(const size_t index) &
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void initSpectrum(const API::MatrixWorkspace &input, API::MatrixWorkspace &output, const size_t j)
SolidAngle::initSpectrum Sets the default value for the spectra for which solid angle is not calculat...
Definition: SolidAngle.cpp:297
void exec() override
Executes the algorithm.
Definition: SolidAngle.cpp:183
void init() override
Initialisation method.
Definition: SolidAngle.cpp:149
ComponentInfo : Provides a component centric view on to the instrument.
Definition: ComponentInfo.h:40
Kernel::V3D scaleFactor(const size_t componentIndex) const
Geometry::DetectorInfo is an intermediate step towards a DetectorInfo that is part of Instrument-2....
Definition: DetectorInfo.h:49
double l2(const size_t index) const
Returns L2 (distance from sample to spectrum).
Kernel::V3D position(const size_t index) const
Returns the position of the detector with given index.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector with given index.
virtual double solidAngle(const Kernel::V3D &observer) const =0
Finds the approximate solid angle covered by the component when viewed from the point given.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
void setX(const double xx) noexcept
Set is x position.
Definition: V3D.h:218
double cosAngle(const V3D &) const
cos(Angle) between this and another vector
Definition: V3D.cpp:180
void setY(const double yy) noexcept
Set is y position.
Definition: V3D.h:224
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
static const std::string HORIZONTAL_WING
Definition: SolidAngle.cpp:34
static const std::string VERTICAL_TUBE
Definition: SolidAngle.cpp:31
static const std::string GENERIC_SHAPE
Definition: SolidAngle.cpp:29
static const std::string HORIZONTAL_TUBE
Definition: SolidAngle.cpp:32
static const std::string RECTANGLE
Definition: SolidAngle.cpp:30
static const std::string VERTICAL_WING
Definition: SolidAngle.cpp:33
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
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
Returns the angle between the sample-to-pixel vector and its projection on the X-Z (vertical tube) or...
Definition: SolidAngle.cpp:52
AlphaAngleCalculator(const DetectorInfo &detectorInfo)
Definition: SolidAngle.cpp:53
double solidAngle(size_t index) const override
Definition: SolidAngle.cpp:107
double solidAngle(size_t index) const override
Definition: SolidAngle.cpp:112
Creates the solid angle calculator based on the selected method.
Definition: SolidAngle.cpp:82
std::unique_ptr< const AlphaAngleCalculator > m_alphaAngleCalculator
Definition: SolidAngle.cpp:102
virtual double solidAngle(size_t index) const =0
SolidAngleCalculator(const ComponentInfo &componentInfo, const DetectorInfo &detectorInfo, const std::string &method, const double pixelArea)
Definition: SolidAngle.cpp:83
double solidAngle(size_t index) const override
Definition: SolidAngle.cpp:124
double solidAngle(size_t index) const override
Definition: SolidAngle.cpp:135
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54