Mantid
Loading...
Searching...
No Matches
CorrectToFile.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"
14#include "MantidKernel/Unit.h"
16
17namespace Mantid::Algorithms {
18using namespace API;
19using namespace DataObjects;
20
21// Register the algorithm into the AlgorithmFactory
22DECLARE_ALGORITHM(CorrectToFile)
23// estimate that this algorithm will spend half it's time loading the file
24const double CorrectToFile::LOAD_TIME = 0.5;
25
27 declareProperty(std::make_unique<API::WorkspaceProperty<>>("WorkspaceToCorrect", "", Kernel::Direction::Input),
28 "Name of the input workspace");
29 declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Load),
30 "The file containing the correction factors");
31
32 std::vector<std::string> propOptions = Kernel::UnitFactory::Instance().getKeys();
33 propOptions.emplace_back("SpectrumNumber");
34 declareProperty("FirstColumnValue", "Wavelength", std::make_shared<Kernel::StringListValidator>(propOptions),
35 "The units of the first column of the correction file "
36 "(default wavelength)");
37
38 std::vector<std::string> operations{"Divide", "Multiply"};
39 declareProperty("WorkspaceOperation", "Divide", std::make_shared<Kernel::StringListValidator>(operations),
40 "Allowed values: Divide, Multiply (default is divide)");
41 declareProperty(std::make_unique<API::WorkspaceProperty<>>("OutputWorkspace", "", Kernel::Direction::Output),
42 "Name of the output workspace to store the results in");
43}
44
46 // The input workspace is the uncorrected data
47 MatrixWorkspace_sptr toCorrect = getProperty("WorkspaceToCorrect");
48 // This workspace is loaded from the RKH compatible file
49 MatrixWorkspace_sptr rkhInput = loadInFile(getProperty("Filename"));
50 // Only create the output workspace if it's not the same as the input one
51 MatrixWorkspace_sptr outputWS = create<HistoWorkspace>(*toCorrect);
52 const std::string operation = getProperty("WorkspaceOperation");
53
54 if (getPropertyValue("FirstColumnValue") == "SpectrumNumber") {
55 // the workspace (probably) contains many spectra but each with only 1 bin
56 doWkspAlgebra(toCorrect, rkhInput, operation, outputWS);
57 } else // interpolation the correction values and divide or multiply the input
58 // by these values
59 { // the correction values should be all contained in 1 spectrum
60 // Check that the workspace to rebin has the same units as the one that we
61 // are matching to
62 // However, just print a warning if it isn't, don't abort (since user
63 // provides the file's unit)
64 if (toCorrect->getAxis(0)->unit()->unitID() != rkhInput->getAxis(0)->unit()->unitID()) {
65 g_log.warning("Unit on input workspace is different to that specified in "
66 "'FirstColumnValue' property");
67 }
68
69 // Get references to the correction factors
70 auto &Xcor = rkhInput->x(0);
71 auto &Ycor = rkhInput->y(0);
72 auto &Ecor = rkhInput->e(0);
73
74 const bool divide = operation == "Divide";
75 double Yfactor, correctError;
76
77 const auto nOutSpec = static_cast<int64_t>(outputWS->getNumberHistograms());
78 const size_t nbins = outputWS->blocksize();
79 // Set the progress bar
80 Progress prg(this, 0 /*LOAD_TIME*/, 1.0, nOutSpec);
81
82 for (int64_t i = 0; i < nOutSpec; ++i) {
83 const auto xIn = toCorrect->points(i);
84 outputWS->setSharedX(i, toCorrect->sharedX(i));
85
86 auto &yOut = outputWS->mutableY(i);
87 auto &eOut = outputWS->mutableE(i);
88
89 auto &yIn = toCorrect->y(i);
90 auto &eIn = toCorrect->e(i);
91
92 for (size_t j = 0; j < nbins; ++j) {
93 const double currentX = xIn[j];
94 // Find out the index of the first correction point after this value
95 auto pos = std::lower_bound(Xcor.cbegin(), Xcor.cend(), currentX);
96 const size_t index = pos - Xcor.begin();
97 if (index == Xcor.size()) {
98 // If we're past the end of the correction factors vector, use the
99 // last point
100 Yfactor = Ycor[index - 1];
101 correctError = Ecor[index - 1];
102 } else if (index) {
103 // Calculate where between the two closest points our current X value
104 // is
105 const double fraction = (currentX - Xcor[index - 1]) / (Xcor[index] - Xcor[index - 1]);
106 // Now linearly interpolate to find the correction factors to use
107 Yfactor = Ycor[index - 1] + fraction * (Ycor[index] - Ycor[index - 1]);
108 correctError = Ecor[index - 1] + fraction * (Ecor[index] - Ecor[index - 1]);
109 } else {
110 // If we're before the start of the correction factors vector, use the
111 // first point
112 Yfactor = Ycor[0];
113 correctError = Ecor[0];
114 }
115
116 // Now do the correction on the current point
117 if (divide) {
118 yOut[j] = yIn[j] / Yfactor;
119 // the proportional error is equal to the sum of the proportional
120 // errors
121 // re-arrange so that you don't get infinity if leftY==0. Sa = error
122 // on a, etc.
123 // c = a/b
124 // (Sa/a)2 + (Sb/b)2 = (Sc/c)2
125 // (Sa c/a)2 + (Sb c/b)2 = (Sc)2
126 // = (Sa 1/b)2 + (Sb (a/b2))2
127 // (Sc)2 = (1/b)2( (Sa)2 + (Sb a/b)2 )
128 eOut[j] = sqrt(pow(eIn[j], 2) + pow(yIn[j] * correctError / Yfactor, 2)) / Yfactor;
129 } else {
130 yOut[j] = yIn[j] * Yfactor;
131 // error multiplying two uncorrelated numbers, re-arrange so that you
132 // don't get infinity if leftY or rightY == 0
133 // Sa = error on a, etc.
134 // c = a*b
135 // (Sa/a)2 + (Sb/b)2 = (Sc/c)2
136 // (Sc)2 = (Sa c/a)2 + (Sb c/b)2 = (Sa b)2 + (Sb a)2
137 eOut[j] = sqrt(pow(eIn[j] * Yfactor, 2) + pow(correctError * yIn[j], 2));
138 }
139 }
140 prg.report("CorrectToFile: applying " + operation);
141 }
142 }
143
144 // Set the resulting workspace
145 setProperty("Outputworkspace", outputWS);
146}
153 g_log.information() << "Loading file " << corrFile << '\n';
154 progress(0, "Loading file");
155 auto loadRKH = createChildAlgorithm("LoadRKH", 0, 1.0 /*LOAD_TIME*/);
156 std::string rkhfile = getProperty("Filename");
157 loadRKH->setPropertyValue("Filename", rkhfile);
158 loadRKH->setPropertyValue("OutputWorkspace", "rkhout");
159 std::string columnValue = getProperty("FirstColumnValue");
160 loadRKH->setPropertyValue("FirstColumnValue", columnValue);
161 loadRKH->executeAsChildAlg();
162
163 g_log.debug() << corrFile << " loaded\n";
164 return loadRKH->getProperty("OutputWorkspace");
165}
175 const std::string &algName, API::MatrixWorkspace_sptr &result) {
176 g_log.information() << "Initalising the algorithm " << algName << '\n';
177 progress(LOAD_TIME, "Applying correction");
178 auto algebra = createChildAlgorithm(algName, LOAD_TIME, 1.0);
179 algebra->setProperty("LHSWorkspace", lhs);
180 algebra->setProperty("RHSWorkspace", rhs);
181 algebra->setProperty("OutputWorkspace", result);
182
183 try {
184 algebra->execute();
185 } catch (std::runtime_error &) {
186 g_log.warning() << "Error encountered while running algorithm " << algName << '\n';
187 g_log.error() << "Correction file " << getPropertyValue("Filename") + " can't be used to correct workspace "
188 << getPropertyValue("WorkspaceToCorrect") << '\n';
189 g_log.error() << "Mismatched number of spectra?\n";
190 throw std::runtime_error("Correct to file failed, see log for details");
191 }
192
193 result = algebra->getProperty("OutputWorkspace");
194 g_log.debug() << algName << " complete\n";
195}
196} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
const std::vector< double > & rhs
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
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
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
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void init() override
Initialisation code.
void exec() override
Execution code.
static const double LOAD_TIME
used for the progress bar: the, approximate, fraction of processing time that taken up with loading t...
Definition: CorrectToFile.h:45
API::MatrixWorkspace_sptr loadInFile(const std::string &corrFile)
Load in the RKH file for that has the correction information.
void doWkspAlgebra(const API::MatrixWorkspace_sptr &lhs, const API::MatrixWorkspace_sptr &rhs, const std::string &algName, API::MatrixWorkspace_sptr &result)
Multiply or divide the input workspace as specified by the user.
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 error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54