Mantid
Loading...
Searching...
No Matches
TOFSANSResolutionByPixel.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 "boost/lexical_cast.hpp"
22
23namespace Mantid::Algorithms {
24
25// Register the algorithm into the AlgorithmFactory
26DECLARE_ALGORITHM(TOFSANSResolutionByPixel)
27
28using namespace Kernel;
29using namespace API;
30using namespace Geometry;
31using namespace DataObjects;
32
34
36 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
37 std::make_shared<WorkspaceUnitValidator>("Wavelength")),
38 "Name the workspace to calculate the resolution for, for "
39 "each pixel and wavelength");
40 declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output),
41 "Name of the newly created workspace which contains the Q resolution.");
42 auto positiveDouble = std::make_shared<BoundedValidator<double>>();
43 positiveDouble->setLower(0);
44 declareProperty("DeltaR", 0.0, positiveDouble, "Virtual ring width on the detector (mm).");
45 declareProperty("SampleApertureRadius", 0.0, positiveDouble, "Sample aperture radius, R2 (mm).");
46 declareProperty("SourceApertureRadius", 0.0, positiveDouble, "Source aperture radius, R1 (mm).");
47 declareProperty(std::make_unique<WorkspaceProperty<>>("SigmaModerator", "", Direction::Input,
48 std::make_shared<WorkspaceUnitValidator>("Wavelength")),
49 "Moderator time spread (microseconds) as a"
50 "function of wavelength (Angstroms).");
51 declareProperty("CollimationLength", 0.0, positiveDouble, "Collimation length (m)");
52 declareProperty("AccountForGravity", false, "Whether to correct for the effects of gravity");
53 declareProperty("ExtraLength", 0.0, positiveDouble, "Additional length for gravity correction.");
54}
55
56/*
57 * Double Boltzmann fit to the TOF resolution as a function of wavelength
58 */
60 UNUSED_ARG(wl);
61 return m_wl_resolution;
62}
63
65 MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
66 double deltaR = getProperty("DeltaR");
67 double R1 = getProperty("SourceApertureRadius");
68 double R2 = getProperty("SampleApertureRadius");
69 const bool doGravity = getProperty("AccountForGravity");
70
71 // Check the input
72 checkInput(inWS);
73
74 // Setup outputworkspace
75 auto outWS = setupOutputWorkspace(inWS);
76
77 // Convert to meters
78 deltaR /= 1000.0;
79 R1 /= 1000.0;
80 R2 /= 1000.0;
81
82 // The moderator workspace needs to match the data workspace
83 // in terms of wavelength binning
84 const MatrixWorkspace_sptr sigmaModeratorVSwavelength = getModeratorWorkspace(inWS);
85
86 // create interpolation table from sigmaModeratorVSwavelength
87 Kernel::Interpolation lookUpTable;
88
89 const auto &xInterpolate = sigmaModeratorVSwavelength->points(0);
90 const auto &yInterpolate = sigmaModeratorVSwavelength->y(0);
91
92 // prefer the input to be a pointworkspace and create interpolation function
93 if (sigmaModeratorVSwavelength->isHistogramData()) {
94 g_log.notice() << "mid-points of SigmaModerator histogram bins will be "
95 "used for interpolation.";
96 }
97
98 for (size_t i = 0; i < xInterpolate.size(); ++i) {
99 lookUpTable.addPoint(xInterpolate[i], yInterpolate[i]);
100 }
101
102 // Calculate the L1 distance
103 const V3D samplePos = inWS->getInstrument()->getSample()->getPos();
104 const V3D sourcePos = inWS->getInstrument()->getSource()->getPos();
105 const V3D SSD = samplePos - sourcePos;
106 const double L1 = SSD.norm();
107
108 // Get the collimation length
109 double LCollim = getProperty("CollimationLength");
110
111 if (LCollim == 0.0) {
112 auto collimationLengthEstimator = SANSCollimationLengthEstimator();
113 LCollim = collimationLengthEstimator.provideCollimationLength(inWS);
114 g_log.information() << "No collimation length was specified. A default "
115 "collimation length was estimated to be "
116 << LCollim << '\n';
117 } else {
118 g_log.information() << "The collimation length is " << LCollim << '\n';
119 }
120
121 const auto numberOfSpectra = static_cast<int>(inWS->getNumberHistograms());
122 Progress progress(this, 0.0, 1.0, numberOfSpectra);
123
124 const auto &spectrumInfo = inWS->spectrumInfo();
125 for (int i = 0; i < numberOfSpectra; i++) {
127 if (!spectrumInfo.hasDetectors(i)) {
128 g_log.information() << "Workspace index " << i << " has no detector assigned to it - discarding\n";
129 continue;
130 }
131 // If no detector found or if it's masked or a monitor, skip onto the next
132 // spectrum
133 if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
134 continue;
135
136 const double L2 = spectrumInfo.l2(i);
138 const double waveLengthIndependentFactor = calculator.getWavelengthIndependentFactor(R1, R2, deltaR, LCollim, L2);
139
140 // Multiplicative factor to go from lambda to Q
141 // Don't get fooled by the function name...
142 const double theta = spectrumInfo.twoTheta(i);
143 double sinTheta = sin(0.5 * theta);
144 double factor = 4.0 * M_PI * sinTheta;
145
146 const auto &xIn = inWS->x(i);
147 const size_t xLength = xIn.size();
148
149 // Gravity correction
150 std::unique_ptr<GravitySANSHelper> grav;
151 if (doGravity) {
152 grav = std::make_unique<GravitySANSHelper>(spectrumInfo, i, getProperty("ExtraLength"));
153 }
154
155 // Get handles on the outputWorkspace
156 auto &yOut = outWS->mutableY(i);
157 // for each wavelenght bin of each pixel calculate a q-resolution
158 for (size_t j = 0; j < xLength - 1; j++) {
159 // use the midpoint of each bin
160 const double wl = (xIn[j + 1] + xIn[j]) / 2.0;
161 // Calculate q. Alternatively q could be calculated using ConvertUnit
162 // If we include a gravity correction we need to adjust sinTheta
163 // for each wavelength (in Angstrom)
164 if (doGravity) {
165 double sinThetaGrav = grav->calcSinTheta(wl);
166 factor = 4.0 * M_PI * sinThetaGrav;
167 }
168 const double q = factor / wl;
169
170 // wavelenght spread from bin assumed to be
171 const double sigmaSpreadFromBin = xIn[j + 1] - xIn[j];
172
173 // Get the uncertainty in Q
174 auto sigmaQ = calculator.getSigmaQValue(lookUpTable.value(wl), waveLengthIndependentFactor, q, wl,
175 sigmaSpreadFromBin, L1, L2);
176
177 // Insert the Q value and the Q resolution into the outputworkspace
178 yOut[j] = sigmaQ;
179 }
180 progress.report("Computing Q resolution");
181 }
182
183 // Set the y axis label
184 outWS->setYUnitLabel("QResolution");
185
186 setProperty("OutputWorkspace", outWS);
187}
188
195 auto duplicate = createChildAlgorithm("CloneWorkspace");
196 duplicate->initialize();
197 duplicate->setProperty<Workspace_sptr>("InputWorkspace", inputWorkspace);
198 duplicate->execute();
199 Workspace_sptr temp = duplicate->getProperty("OutputWorkspace");
200 return std::dynamic_pointer_cast<MatrixWorkspace>(temp);
201}
202
210
211 MatrixWorkspace_sptr sigmaModerator = getProperty("SigmaModerator");
212 auto rebinned = createChildAlgorithm("RebinToWorkspace");
213 rebinned->initialize();
214 rebinned->setProperty("WorkspaceToRebin", sigmaModerator);
215 rebinned->setProperty("WorkspaceToMatch", inputWorkspace);
216 rebinned->setPropertyValue("OutputWorkspace", "SigmaModeratorRebinned");
217 rebinned->execute();
218 MatrixWorkspace_sptr outWS = rebinned->getProperty("OutputWorkspace");
219 return outWS;
220}
221
227 // Make sure that input workspace has an instrument as we rely heavily on
228 // thisa
229 auto inst = inWS->getInstrument();
230 if (inst->getName().empty()) {
231 throw std::invalid_argument("TOFSANSResolutionByPixel: The input workspace "
232 "does not contain an instrument");
233 }
234}
235
236} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
Helper class which provides the uncertainty calculations for the TOFSANSResolutionByPixel class.
double getSigmaQValue(double moderatorValue, double wavlengthIndependentFactor, double q, double wavelength, double deltaWavelength, double l1, double l2) const
double getWavelengthIndependentFactor(double r1, double r2, double deltaR, double lCollim, double l2) const
Calculates the wavelength-independent prefactor for the resolutin calculation.
Mantid::API::MatrixWorkspace_sptr setupOutputWorkspace(const Mantid::API::MatrixWorkspace_sptr &inputWorkspace)
Create an output workspace.
Mantid::API::MatrixWorkspace_sptr getModeratorWorkspace(const Mantid::API::MatrixWorkspace_sptr &inputWorkspace)
Get the moderator workspace.
void checkInput(const Mantid::API::MatrixWorkspace_sptr &inWS)
Check input.
virtual double getTOFResolution(double wl)
Return the TOF resolution for a particular wavelength.
double m_wl_resolution
Wavelength resolution (constant for all wavelengths)
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
Provide interpolation over a series of points.
Definition: Interpolation.h:29
double value(const double &at) const
get interpolated value at location at
void addPoint(const double &xx, const double &yy)
add data point
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Class for 3D vectors.
Definition: V3D.h:34
double norm() const noexcept
Definition: V3D.h:263
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< const Mantid::Geometry::IDetector > IDetector_const_sptr
Shared pointer to IDetector (const version)
Definition: IDetector.h:102
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54