Mantid
Loading...
Searching...
No Matches
EstimateResolutionDiffraction.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 +
9#include "MantidAPI/Run.h"
22#include "MantidKernel/V3D.h"
23#include "MantidTypes/SpectrumDefinition.h"
24
25#include <cmath>
26
27using namespace Mantid;
28using namespace Mantid::API;
29using namespace Mantid::Geometry;
30using namespace Mantid::Kernel;
31using namespace std;
32
33namespace Mantid::Algorithms {
34
35DECLARE_ALGORITHM(EstimateResolutionDiffraction)
36
37namespace { // hide these constants
39const double MICROSEC_TO_SEC = 1.0E-6;
41const double WAVELENGTH_TO_VELOCITY = 1.0E10 * PhysicalConstants::h / PhysicalConstants::NeutronMass;
43const double WAVELENGTH_MAX = 1000.;
44} // namespace
45
46const std::string EstimateResolutionDiffraction::name() const { return "EstimateResolutionDiffraction"; }
47
48const std::string EstimateResolutionDiffraction::alias() const { return "EstimatePDDetectorResolution"; }
49
50const std::string EstimateResolutionDiffraction::summary() const {
51 return "Estimate the resolution of each detector pixel for a powder "
52 "diffractometer";
53}
54
56
57const std::string EstimateResolutionDiffraction::category() const { return "Diffraction\\Utility"; }
58
60 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
61 "Name of the workspace to have detector resolution calculated");
62 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("DivergenceWorkspace", "", Direction::Input,
64 "Workspace containing the divergence");
65 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
66 "Name of the output workspace containing delta(d)/d of each "
67 "detector/spectrum");
68
69 auto positiveDeltaTOF = std::make_shared<BoundedValidator<double>>();
70 positiveDeltaTOF->setLower(0.);
71 positiveDeltaTOF->setLowerExclusive(true);
72 declareProperty("DeltaTOF", 0., positiveDeltaTOF, "DeltaT as the resolution of TOF with unit microsecond");
73
74 auto positiveWavelength = std::make_shared<BoundedValidator<double>>();
75 positiveWavelength->setLower(0.);
76 positiveWavelength->setLowerExclusive(true);
77 declareProperty("Wavelength", EMPTY_DBL(), positiveWavelength,
78 "Wavelength setting in Angstroms. This overrides what is in "
79 "the dataset.");
80
82 std::make_unique<WorkspaceProperty<WorkspaceGroup>>("PartialResolutionWorkspaces", "", Direction::Output),
83 "Workspaces created showing the various resolution terms");
84}
85
88
90
91 // create all of the output workspaces
92 std::string partials_prefix = getPropertyValue("PartialResolutionWorkspaces");
93 m_resTof = DataObjects::create<DataObjects::Workspace2D>(*m_inputWS, HistogramData::Points(1));
94 m_resPathLength = DataObjects::create<DataObjects::Workspace2D>(*m_inputWS, HistogramData::Points(1));
95 m_resAngle = DataObjects::create<DataObjects::Workspace2D>(*m_inputWS, HistogramData::Points(1));
96 m_outputWS = DataObjects::create<DataObjects::Workspace2D>(*m_inputWS, HistogramData::Points(1));
97
99
100 setProperty("OutputWorkspace", m_outputWS);
101
102 // put together the output group
103 auto partialsGroup = std::make_shared<WorkspaceGroup>();
104 API::AnalysisDataService::Instance().addOrReplace(partials_prefix + "_tof", m_resTof);
105 API::AnalysisDataService::Instance().addOrReplace(partials_prefix + "_length", m_resPathLength);
106 API::AnalysisDataService::Instance().addOrReplace(partials_prefix + "_angle", m_resAngle);
107 partialsGroup->addWorkspace(m_resTof);
108 partialsGroup->addWorkspace(m_resPathLength);
109 partialsGroup->addWorkspace(m_resAngle);
110 setProperty("PartialResolutionWorkspaces", partialsGroup);
111}
112
114 m_inputWS = getProperty("InputWorkspace");
115 m_divergenceWS = getProperty("DivergenceWorkspace");
116
117 m_deltaT = getProperty("DeltaTOF");
118 m_deltaT *= MICROSEC_TO_SEC; // convert to meter
119}
120
122 double wavelength = getProperty("Wavelength");
123 if (!isEmpty(wavelength)) {
124 return wavelength;
125 }
126
127 Property *cwlproperty = m_inputWS->run().getProperty("LambdaRequest");
128 if (!cwlproperty)
129 throw runtime_error("Unable to locate property LambdaRequest as central wavelength. ");
130
131 auto *cwltimeseries = dynamic_cast<TimeSeriesProperty<double> *>(cwlproperty);
132
133 if (!cwltimeseries)
134 throw runtime_error("LambdaReqeust is not a TimeSeriesProperty in double. ");
135
136 string unit = cwltimeseries->units();
137 if ((unit != "Angstrom") && (unit != "A")) {
138 throw runtime_error("Unit is not recognized: " + unit);
139 }
140
141 return cwltimeseries->timeAverageValue();
142}
143
145 double centrewavelength = getWavelength();
146 g_log.notice() << "Centre wavelength = " << centrewavelength << " Angstrom\n";
147 if (centrewavelength > WAVELENGTH_MAX) {
148 throw runtime_error("unphysical wavelength used");
149 }
150
151 // Calculate centre neutron velocity
152 m_centreVelocity = WAVELENGTH_TO_VELOCITY / centrewavelength;
153 g_log.notice() << "Centre neutron velocity = " << m_centreVelocity << "\n";
154}
155
157 const auto &spectrumInfo = m_inputWS->spectrumInfo();
158 const auto l1 = spectrumInfo.l1();
159 const auto &componentInfo = m_inputWS->componentInfo();
160 const auto &detectorInfo = m_inputWS->detectorInfo();
161 g_log.notice() << "L1 = " << l1 << "\n";
162 const auto samplepos = spectrumInfo.samplePosition();
163
164 const size_t numspec = m_inputWS->getNumberHistograms();
165
166 double mintwotheta = 2. * M_PI; // a bit more than 2*pi
167 double maxtwotheta = 0.;
168
169 double mint3 = 1.;
170 double maxt3 = 0.;
171
172 size_t count_nodetsize = 0;
173
174 for (size_t i = 0; i < numspec; ++i) {
175 const auto &det = spectrumInfo.detector(i);
176 double detdim;
177 const auto realdet = dynamic_cast<const Detector *>(&det);
178 if (realdet) {
179 const double dy = realdet->getHeight();
180 const double dx = realdet->getWidth();
181 detdim = sqrt(dx * dx + dy * dy) * 0.5;
182 } else {
183 // Use detector dimension as 0 as no-information
184 detdim = 0;
185 ++count_nodetsize;
186 }
187
188 // Get the distance from detector to source
189 const double l2 = spectrumInfo.l2(i);
190
191 // Calculate T
192 const double centraltof = (l1 + l2) / m_centreVelocity;
193
194 // Angle
195 const double twotheta = spectrumInfo.isMonitor(i) ? 0.0 : spectrumInfo.twoTheta(i);
196 const double theta = 0.5 * twotheta;
197
198 double deltatheta = 0.;
199 if (m_divergenceWS) {
200 deltatheta = m_divergenceWS->y(i)[0];
201 } else {
202 auto &spectrumDefinition = spectrumInfo.spectrumDefinition(i);
203 const double solidangle =
204 std::accumulate(spectrumDefinition.cbegin(), spectrumDefinition.cend(), 0.,
205 [&componentInfo, &detectorInfo, &samplepos](const auto sum, const auto &index) {
206 if (!detectorInfo.isMasked(index.first)) {
207 return sum + componentInfo.solidAngle(index.first, samplepos);
208 } else {
209 return sum;
210 }
211 });
212 deltatheta = sqrt(solidangle);
213 }
214
215 // Resolution
216 const double t1 = m_deltaT / centraltof;
217 const double t2 = detdim / (l1 + l2);
218 const double t3 = deltatheta / tan(theta);
219
220 if (spectrumInfo.isMonitor(i)) {
221 m_resTof->mutableY(i) = 0.;
222 m_resPathLength->mutableY(i) = 0.;
223 m_resAngle->mutableY(i) = 0.;
224 m_outputWS->mutableY(i) = 0.;
225 } else { // not a monitor
226 const double resolution = sqrt(t1 * t1 + t2 * t2 + t3 * t3);
227 m_resTof->mutableY(i) = t1;
228 m_resPathLength->mutableY(i) = t2;
229 m_resAngle->mutableY(i) = t3;
230 m_outputWS->mutableY(i) = resolution;
231 }
232
233 m_resTof->mutableX(i) = static_cast<double>(i);
234 m_resPathLength->mutableX(i) = static_cast<double>(i);
235 m_resAngle->mutableX(i) = static_cast<double>(i);
236 m_outputWS->mutableX(i) = static_cast<double>(i);
237
238 maxtwotheta = std::max(twotheta, maxtwotheta);
239 mintwotheta = std::min(twotheta, mintwotheta);
240
241 if (fabs(t3) < mint3)
242 mint3 = fabs(t3);
243 else if (fabs(t3) > maxt3)
244 maxt3 = fabs(t3);
245
246 g_log.debug() << det.type() << " " << i << "\t\t" << twotheta << "\t\tdT/T = " << t1 * t1 << "\t\tdL/L = " << t2
247 << "\t\tdTheta*cotTheta = " << t3 << "\n";
248 }
249
250 g_log.notice() << "2theta range: " << mintwotheta << ", " << maxtwotheta << "\n";
251 g_log.notice() << "t3 range: " << mint3 << ", " << maxt3 << "\n";
252 g_log.notice() << "Number of detector having NO size information = " << count_nodetsize << "\n";
253}
254
255} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define fabs(x)
Definition: Matrix.cpp:22
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
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
A property class for workspaces.
API::MatrixWorkspace_sptr m_resTof
workspace holding the term for just the time-of-flight portion of the resolution
API::MatrixWorkspace_sptr m_resAngle
workspace holding the term for just the angular/solid angle portion of the resolution
double getWavelength()
Returns the wavelength from either the property or the input workspace.
API::MatrixWorkspace_sptr m_outputWS
Output workspace.
void processAlgProperties()
Process input properties for algorithm.
const std::string category() const override
Algorithm's category for identification overriding a virtual method.
void exec() override
Implement abstract Algorithm methods.
API::MatrixWorkspace_sptr m_inputWS
Input workspace.
void init() override
Implement abstract Algorithm methods.
const std::string alias() const override
function to return any aliases to the algorithm
API::MatrixWorkspace_sptr m_divergenceWS
Workspace with custom divergence term.
API::MatrixWorkspace_sptr m_resPathLength
workspace holding the term for just the flight path portion of the resolution
const std::string summary() const override
Summary of algorithms purpose.
int version() const override
Algorithm's version for identification overriding a virtual method.
const std::string name() const override
Algorithm's name for identification overriding a virtual method.
This class represents a detector - i.e.
Definition: Detector.h:30
virtual double getHeight() const
get Height (Y-dimension) value for component
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
Base class for properties.
Definition: Property.h:94
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A specialised Property class for holding a series of time-value pairs.
Kernel::Logger g_log("ExperimentInfo")
static logger object
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
Helper class which provides the Collimation Length for SANS instruments.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54