21using namespace Kernel;
23using namespace Geometry;
24using namespace DataObjects;
33static double getYTubeAngle(
const SpectrumInfo &spectrumInfo,
size_t i) {
34 const V3D samplePos = spectrumInfo.samplePosition();
37 V3D sampleDetVec = spectrumInfo.position(i) - samplePos;
40 V3D inPlane =
V3D(sampleDetVec);
45 return sampleDetVec.
angle(inPlane);
52 auto wsValidator = std::make_shared<CompositeValidator>();
58 "If true, the algorithm will assume "
59 "that the detectors are tubes in the "
62 "If true, the algorithm will assume "
63 "that the detector is curved around the sample. E.g. BIOSANS "
71 const std::string reductionManagerName =
getProperty(
"ReductionProperties");
72 std::shared_ptr<PropertyManager> reductionManager;
76 reductionManager = std::make_shared<PropertyManager>();
81 if (!reductionManager->existsProperty(
"SANSSolidAngleCorrection")) {
82 auto algProp = std::make_unique<AlgorithmProperty>(
"SANSSolidAngleCorrection");
84 reductionManager->declareProperty(std::move(algProp));
94 if (outputWS != inputWS) {
96 outputWS->setDistribution(
true);
97 outputWS->setYUnit(
"");
98 outputWS->setYUnitLabel(
"Steradian");
102 const auto numHists =
static_cast<int>(inputWS->getNumberHistograms());
106 const auto xLength =
static_cast<int>(inputWS->y(0).size());
108 const auto &spectrumInfo = inputWS->spectrumInfo();
110 for (
int i = 0; i < numHists; ++i) {
112 outputWS->setSharedX(i, inputWS->sharedX(i));
114 if (!spectrumInfo.hasDetectors(i)) {
115 g_log.
warning() <<
"Workspace index " << i <<
" has no detector assigned to it - discarding\n";
120 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
123 const auto &YIn = inputWS->y(i);
124 const auto &EIn = inputWS->e(i);
126 auto &YOut = outputWS->mutableY(i);
127 auto &EOut = outputWS->mutableE(i);
133 const double tanTheta = tan(spectrumInfo.twoTheta(i));
134 const double theta_term = sqrt(tanTheta * tanTheta + 1.0);
136 if (is_tube || is_wing) {
137 const double tanAlpha = tan(getYTubeAngle(spectrumInfo, i));
138 const double alpha_term = sqrt(tanAlpha * tanAlpha + 1.0);
140 corr = alpha_term * theta_term * theta_term;
142 corr = alpha_term * alpha_term * alpha_term;
145 corr = theta_term * theta_term * theta_term;
149 for (
int j = 0; j < xLength; j++) {
150 YOut[j] = YIn[j] * corr;
151 EOut[j] =
fabs(EIn[j] * corr);
153 progress.report(
"Solid Angle Correction");
157 setProperty(
"OutputMessage",
"Solid angle correction applied");
165 if (outputWS != inputWS) {
166 outputWS = inputWS->clone();
169 auto outputEventWS = std::dynamic_pointer_cast<EventWorkspace>(outputWS);
171 const auto numberOfSpectra =
static_cast<int>(outputEventWS->getNumberHistograms());
173 progress.report(
"Solid Angle Correction");
175 const auto &spectrumInfo = outputEventWS->spectrumInfo();
177 for (
int i = 0; i < numberOfSpectra; i++) {
180 if (!spectrumInfo.hasDetectors(i)) {
181 g_log.
warning() <<
"Workspace index " << i <<
" has no detector assigned to it - discarding\n";
186 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
191 const double tanTheta = tan(spectrumInfo.twoTheta(i));
192 const double theta_term = sqrt(tanTheta * tanTheta + 1.0);
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;
199 corr = theta_term * theta_term * theta_term;
201 EventList &el = outputEventWS->getSpectrum(i);
203 progress.report(
"Solid Angle Correction");
208 setProperty(
"OutputMessage",
"Solid angle correction applied");
#define DECLARE_ALGORITHM(classname)
#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...
#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.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
std::string toString() const override
Serialize an object to a string.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
Helper class for reporting progress from algorithms.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
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.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
double angle(const V3D &) const
Angle between this and another vector.
void setY(const double yy) noexcept
Set is y position.
void exec() override
Execution code.
void init() override
Initialisation code.
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.
@ Input
An input workspace.
@ Output
An output workspace.