Mantid
Loading...
Searching...
No Matches
LorentzCorrection.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 +
8#include "MantidAPI/Axis.h"
10#include "MantidAPI/Progress.h"
16#include "MantidKernel/Unit.h"
17#include <cmath>
18#include <memory>
19
20namespace Mantid::Algorithms {
21using namespace Mantid::Kernel;
22using namespace Mantid::API;
23using namespace Mantid::Geometry;
24using namespace Mantid::Kernel;
28
29// Register the algorithm into the AlgorithmFactory
30DECLARE_ALGORITHM(LorentzCorrection)
31
32namespace { // anonymous
33namespace PropertyNames {
34const std::string INPUT_WKSP("InputWorkspace");
35const std::string OUTPUT_WKSP("OutputWorkspace");
36const std::string TYPE("Type");
37} // namespace PropertyNames
38
39const std::string TOF_SCD("SingleCrystalTOF");
40const std::string TOF_PD("PowderTOF");
41
42inline double sinTheta(const API::SpectrumInfo &spectrumInfo, int64_t index) {
43 const double twoTheta = spectrumInfo.isMonitor(index) ? 0.0 : spectrumInfo.twoTheta(index);
44 return std::sin(0.5 * twoTheta);
45}
46} // namespace
47
49int LorentzCorrection::version() const { return 1; }
50
52const std::string LorentzCorrection::category() const { return "Crystal\\Corrections"; }
53
55const std::string LorentzCorrection::summary() const { return "Performs a white beam Lorentz Correction"; }
56
57const std::string LorentzCorrection::name() const { return "LorentzCorrection"; }
58
62 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>(PropertyNames::INPUT_WKSP, "", Direction::Input,
64 // std::make_shared<WorkspaceUnitValidator>("Wavelength")),
65 "Input workspace to correct in Wavelength.");
67 std::make_unique<WorkspaceProperty<MatrixWorkspace>>(PropertyNames::OUTPUT_WKSP, "", Direction::Output),
68 "An output workspace.");
69 const std::vector<std::string> correction_types{TOF_SCD, TOF_PD};
70 declareProperty(PropertyNames::TYPE, correction_types[0], std::make_shared<StringListValidator>(correction_types),
71 "Type of Lorentz correction to do");
72}
73
74std::map<std::string, std::string> LorentzCorrection::validateInputs() {
75 std::map<std::string, std::string> result;
76
77 const auto processingType = this->getPropertyValue(PropertyNames::TYPE);
78 // check units if the SCD option is selected
79 if (processingType == TOF_SCD) {
80 MatrixWorkspace_const_sptr wksp = this->getProperty(PropertyNames::INPUT_WKSP);
81 if (!wksp) {
82 result[PropertyNames::INPUT_WKSP] = "The workspace must be a MatrixWorkspace.";
83 return result;
84 }
85 // code is a variant of private method from WorkspaceUnitValidator
86 const auto unit = wksp->getAxis(0)->unit();
87 if ((!unit) || (unit->unitID().compare("Wavelength"))) {
88 result[PropertyNames::INPUT_WKSP] = "The workspace must have units of Wavelength";
89 }
90 }
91
92 return result;
93}
94
98 MatrixWorkspace_sptr inWS = this->getProperty(PropertyNames::INPUT_WKSP);
99
100 // clone the workspace - does nothing for inplace operation
101 auto cloneAlg = this->createChildAlgorithm("CloneWorkspace", 0, 0.1);
102 cloneAlg->initialize();
103 cloneAlg->setProperty("InputWorkspace", inWS);
104 cloneAlg->setPropertyValue("OutputWorkspace", this->getPropertyValue(PropertyNames::OUTPUT_WKSP));
105 cloneAlg->execute();
106 Workspace_sptr temp = cloneAlg->getProperty("OutputWorkspace");
107 MatrixWorkspace_sptr outWS = std::dynamic_pointer_cast<MatrixWorkspace>(temp);
108
109 Progress prog(this, 0.1, 1.0, inWS->getNumberHistograms());
110
111 const auto processingType = this->getPropertyValue(PropertyNames::TYPE);
112 if (processingType == TOF_SCD)
113 processTOF_SCD(outWS, prog);
114 else if (processingType == TOF_PD)
115 processTOF_PD(outWS, prog);
116 else {
117 std::stringstream msg;
118 msg << "Do not understand know how to process Type=\"" << processingType
119 << "\" - developer forgot to fill in if/else tree";
120 throw std::runtime_error(msg.str());
121 }
122
123 this->setProperty(PropertyNames::OUTPUT_WKSP, outWS);
124}
125
127 const int64_t numHistos = static_cast<int64_t>(wksp->getNumberHistograms());
128 const auto &spectrumInfo = wksp->spectrumInfo();
129
131 for (int64_t i = 0; i < numHistos; ++i) {
133
134 if (!spectrumInfo.hasDetectors(i))
135 continue;
136
137 const double sinTheta_v = sinTheta(spectrumInfo, i);
138 const double sinThetaSq = sinTheta_v * sinTheta_v;
139
140 auto &outY = wksp->mutableY(i);
141 auto &outE = wksp->mutableE(i);
142 const auto points = wksp->points(i);
143 const auto num_points = points.size();
144 const auto pos = std::find(cbegin(points), cend(points), 0.0);
145 if (pos != cend(points)) {
146 std::stringstream buffer;
147 buffer << "Cannot have zero values Wavelength. At workspace index: " << pos - cbegin(points);
148 throw std::runtime_error(buffer.str());
149 }
150 for (size_t j = 0; j < num_points; ++j) {
151 double weight = sinThetaSq / (points[j] * points[j] * points[j] * points[j]);
152 outY[j] *= weight;
153 outE[j] *= weight;
154 }
155
156 prog.report();
157
159 }
161}
162
164 const int64_t numHistos = static_cast<int64_t>(wksp->getNumberHistograms());
165 const auto &spectrumInfo = wksp->spectrumInfo();
166
167 EventWorkspace_sptr wkspEvent = std::dynamic_pointer_cast<EventWorkspace>(wksp);
168 bool isEvent = bool(wkspEvent);
169
171 for (int64_t i = 0; i < numHistos; ++i) {
173
174 if (!spectrumInfo.hasDetectors(i))
175 continue;
176
177 const double sinTheta_v = sinTheta(spectrumInfo, i);
178
179 const auto points = wksp->points(i);
180 const auto pos = std::find(cbegin(points), cend(points), 0.0);
181 if (pos != cend(points)) {
182 std::stringstream buffer;
183 buffer << "Cannot have zero values Wavelength. At workspace index: " << pos - cbegin(points);
184 throw std::runtime_error(buffer.str());
185 }
186
187 if (isEvent) {
188 wkspEvent->getSpectrum(i) *= sinTheta_v;
189 } else {
190 wksp->mutableY(i) *= sinTheta_v;
191 wksp->mutableE(i) *= sinTheta_v;
192 }
193
194 prog.report();
195
197 }
199}
200
201} // 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
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
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
Definition: SpectrumInfo.h:53
bool isMonitor(const size_t index) const
Returns true if the detector(s) associated with the spectrum are monitors.
bool hasDetectors(const size_t index) const
Returns true if the spectrum is associated with detectors in the instrument.
double twoTheta(const size_t index) const
Returns the scattering angle 2 theta in radians (angle w.r.t.
A property class for workspaces.
void processTOF_PD(Mantid::API::MatrixWorkspace_sptr &wksp, Mantid::API::Progress &prog)
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
void processTOF_SCD(Mantid::API::MatrixWorkspace_sptr &wksp, Mantid::API::Progress &prog)
int version() const override
Algorithm's version for identification.
void exec() override
Execute the algorithm.
const std::string name() const override
function to return a name of the algorithm, must be overridden in all algorithms
const std::string category() const override
Algorithm's category for identification.
void init() override
Initialize the algorithm's properties.
std::map< std::string, std::string > validateInputs() override
Method checking errors on ALL the inputs, before execution.
This class is intended to fulfill the design specified in <https://github.com/mantidproject/documents...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
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< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
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