Mantid
Loading...
Searching...
No Matches
SANSSolidAngleCorrection.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 +
18
20
21using namespace Kernel;
22using namespace API;
23using namespace Geometry;
24using namespace DataObjects;
25
26// Register the algorithm into the AlgorithmFactory
27DECLARE_ALGORITHM(SANSSolidAngleCorrection)
28
29namespace {
30
33static double getYTubeAngle(const SpectrumInfo &spectrumInfo, size_t i) {
34 const V3D samplePos = spectrumInfo.samplePosition();
35
36 // Get the vector from the sample position to the detector pixel
37 V3D sampleDetVec = spectrumInfo.position(i) - samplePos;
38
39 // Get the projection of that vector on the X-Z plane
40 V3D inPlane = V3D(sampleDetVec);
41 inPlane.setY(0.0);
42
43 // This is the angle between the sample-to-detector vector
44 // and its project on the X-Z plane.
45 return sampleDetVec.angle(inPlane);
46}
47
48} // namespace
49
50//----------------------------------------------------------------------------------------------
52 auto wsValidator = std::make_shared<CompositeValidator>();
53 wsValidator->add<WorkspaceUnitValidator>("Wavelength");
54 wsValidator->add<HistogramValidator>();
55 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator));
56 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output));
57 declareProperty("DetectorTubes", false,
58 "If true, the algorithm will assume "
59 "that the detectors are tubes in the "
60 "Y direction.");
61 declareProperty("DetectorWing", false,
62 "If true, the algorithm will assume "
63 "that the detector is curved around the sample. E.g. BIOSANS "
64 "Wing detector.");
65 declareProperty("OutputMessage", "", Direction::Output);
66 declareProperty("ReductionProperties", "__sans_reduction_properties", Direction::Input);
67}
68
70 // Reduction property manager
71 const std::string reductionManagerName = getProperty("ReductionProperties");
72 std::shared_ptr<PropertyManager> reductionManager;
73 if (PropertyManagerDataService::Instance().doesExist(reductionManagerName)) {
74 reductionManager = PropertyManagerDataService::Instance().retrieve(reductionManagerName);
75 } else {
76 reductionManager = std::make_shared<PropertyManager>();
77 PropertyManagerDataService::Instance().addOrReplace(reductionManagerName, reductionManager);
78 }
79
80 // If the solid angle algorithm isn't in the reduction properties, add it
81 if (!reductionManager->existsProperty("SANSSolidAngleCorrection")) {
82 auto algProp = std::make_unique<AlgorithmProperty>("SANSSolidAngleCorrection");
83 algProp->setValue(toString());
84 reductionManager->declareProperty(std::move(algProp));
85 }
86
87 MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
88 DataObjects::EventWorkspace_const_sptr inputEventWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWS);
89 if (inputEventWS)
90 return execEvent();
91
92 // Now create the output workspace
93 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
94 if (outputWS != inputWS) {
95 outputWS = WorkspaceFactory::Instance().create(inputWS);
96 outputWS->setDistribution(true);
97 outputWS->setYUnit("");
98 outputWS->setYUnitLabel("Steradian");
99 setProperty("OutputWorkspace", outputWS);
100 }
101
102 const auto numHists = static_cast<int>(inputWS->getNumberHistograms());
103 Progress progress(this, 0.0, 1.0, numHists);
104
105 // Number of X bins
106 const auto xLength = static_cast<int>(inputWS->y(0).size());
107
108 const auto &spectrumInfo = inputWS->spectrumInfo();
109 PARALLEL_FOR_IF(Kernel::threadSafe(*outputWS, *inputWS))
110 for (int i = 0; i < numHists; ++i) {
112 outputWS->setSharedX(i, inputWS->sharedX(i));
113
114 if (!spectrumInfo.hasDetectors(i)) {
115 g_log.warning() << "Workspace index " << i << " has no detector assigned to it - discarding\n";
116 continue;
117 }
118
119 // Skip if we have a monitor or if the detector is masked.
120 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
121 continue;
122
123 const auto &YIn = inputWS->y(i);
124 const auto &EIn = inputWS->e(i);
125
126 auto &YOut = outputWS->mutableY(i);
127 auto &EOut = outputWS->mutableE(i);
128
129 // Compute solid angle correction factor
130 const bool is_tube = getProperty("DetectorTubes");
131 const bool is_wing = getProperty("DetectorWing");
132
133 const double tanTheta = tan(spectrumInfo.twoTheta(i));
134 const double theta_term = sqrt(tanTheta * tanTheta + 1.0);
135 double corr;
136 if (is_tube || is_wing) {
137 const double tanAlpha = tan(getYTubeAngle(spectrumInfo, i));
138 const double alpha_term = sqrt(tanAlpha * tanAlpha + 1.0);
139 if (is_tube)
140 corr = alpha_term * theta_term * theta_term;
141 else { // is_wing
142 corr = alpha_term * alpha_term * alpha_term;
143 }
144 } else {
145 corr = theta_term * theta_term * theta_term;
146 }
147
148 // Correct data for all X bins
149 for (int j = 0; j < xLength; j++) {
150 YOut[j] = YIn[j] * corr;
151 EOut[j] = fabs(EIn[j] * corr);
152 }
153 progress.report("Solid Angle Correction");
155 }
157 setProperty("OutputMessage", "Solid angle correction applied");
158}
159
161 MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
162
163 // generate the output workspace pointer
164 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
165 if (outputWS != inputWS) {
166 outputWS = inputWS->clone();
167 setProperty("OutputWorkspace", outputWS);
168 }
169 auto outputEventWS = std::dynamic_pointer_cast<EventWorkspace>(outputWS);
170
171 const auto numberOfSpectra = static_cast<int>(outputEventWS->getNumberHistograms());
172 Progress progress(this, 0.0, 1.0, numberOfSpectra);
173 progress.report("Solid Angle Correction");
174
175 const auto &spectrumInfo = outputEventWS->spectrumInfo();
176 PARALLEL_FOR_IF(Kernel::threadSafe(*outputEventWS))
177 for (int i = 0; i < numberOfSpectra; i++) {
179
180 if (!spectrumInfo.hasDetectors(i)) {
181 g_log.warning() << "Workspace index " << i << " has no detector assigned to it - discarding\n";
182 continue;
183 }
184
185 // Skip if we have a monitor or if the detector is masked.
186 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
187 continue;
188
189 // Compute solid angle correction factor
190 const bool is_tube = getProperty("DetectorTubes");
191 const double tanTheta = tan(spectrumInfo.twoTheta(i));
192 const double theta_term = sqrt(tanTheta * tanTheta + 1.0);
193 double corr;
194 if (is_tube) {
195 const double tanAlpha = tan(getYTubeAngle(spectrumInfo, i));
196 const double alpha_term = sqrt(tanAlpha * tanAlpha + 1.0);
197 corr = alpha_term * theta_term * theta_term;
198 } else {
199 corr = theta_term * theta_term * theta_term;
200 }
201 EventList &el = outputEventWS->getSpectrum(i);
202 el *= corr;
203 progress.report("Solid Angle Correction");
205 }
207
208 setProperty("OutputMessage", "Solid angle correction applied");
209}
210
211} // namespace Mantid::WorkflowAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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_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
std::string toString() const override
Serialize an object to a string.
Definition: Algorithm.cpp:905
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
A class for holding :
Definition: EventList.h:56
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
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
double angle(const V3D &) const
Angle between this and another vector.
Definition: V3D.cpp:165
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
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
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