Mantid
Loading...
Searching...
No Matches
PointByPointVCorrection.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 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
11#include "MantidAPI/Axis.h"
16
17#include <cfloat>
18#include <cmath>
19#include <numeric>
20
21namespace Mantid::Algorithms {
22
23// Register with the algorithm factory
24DECLARE_ALGORITHM(PointByPointVCorrection)
25
26using namespace Kernel;
27using namespace API;
28using namespace DataObjects;
29
32
33// Destructor
35
37 declareProperty(std::make_unique<WorkspaceProperty<>>("InputW1", "", Direction::Input),
38 "Name of the Sample workspace.");
39 declareProperty(std::make_unique<WorkspaceProperty<>>("InputW2", "", Direction::Input),
40 "Name of the Vanadium workspace.");
41 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
42 "Name of the output workspace.");
43}
44
46 // Get the input workspace and output workspace
47 MatrixWorkspace_const_sptr inputWS1 = getProperty("InputW1");
48 MatrixWorkspace_const_sptr inputWS2 = getProperty("InputW2");
49 MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
50
51 // Check that everything is OK.
52 check_validity(inputWS1, inputWS2, outputWS);
53
54 // Now do the normalisation
55 const auto size = static_cast<int>(inputWS1->x(0).size());
56 const auto nHist = static_cast<int>(inputWS1->getNumberHistograms());
57 Progress prog(this, 0.0, 1.0, nHist);
58
59 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS1, *inputWS2, *outputWS))
60 for (int i = 0; i < nHist; i++) // Looping on all histograms
61 {
63
64 outputWS->setSharedX(i, inputWS1->sharedX(i));
65
66 const auto &X = inputWS1->x(i);
67 const auto &Y1 = inputWS1->y(i);
68 const auto &Y2 = inputWS2->y(i);
69 const auto &E1 = inputWS1->e(i);
70 const auto &E2 = inputWS2->e(i);
71 auto &resultY = outputWS->mutableY(i);
72 auto &resultE = outputWS->mutableE(i);
73
74 // Work on the Y data
75 MantidVec binwidths(size); // MantidVec for bin widths
76 MantidVec errors(size - 1); // MantidVec for temporary errors
77 std::adjacent_difference(X.begin(), X.end(),
78 binwidths.begin()); // Calculate the binwidths
79 std::transform(binwidths.begin() + 1, binwidths.end(), Y2.begin(), resultY.begin(),
81 std::transform(Y1.begin(), Y1.end(), resultY.begin(), resultY.begin(),
82 std::multiplies<double>()); // Now resultY contains the
83 // A_i=s_i/v_i*Dlam_i
84
85 // Calculate the errors squared related to A_i at this point
86 for (int j = 0; j < size - 1; j++) {
87 double r = 0.0;
88 if (std::abs(Y1[j]) > 1e-7)
89 r += std::pow(E1[j] / Y1[j], 2);
90 if (std::abs(Y2[j]) > 1e-7)
91 r += std::pow(E2[j] / Y2[j], 2);
92 errors[j] = r; // This are the errors^2 of S_i/v_i*Dlam_i
93 if (errors[j] > DBL_MAX || errors[j] < -DBL_MAX)
94 errors[j] = 0;
95 }
96
97 // Calculate the normaliser
98 double factor1 = std::accumulate(Y1.begin(), Y1.end(), 0.0);
99 double factor2 = std::accumulate(resultY.begin(), resultY.end(), 0.0);
100 double factor = factor1 / factor2;
101
102 // Now propagate the error bars due to the normaliser
103 double error2_factor1 = std::inner_product(E1.begin(), E1.end(), E1.begin(), 0.0);
104 double error2_factor2 = 0;
105
106 for (int j = 0; j < size - 1; j++) {
107 double test = std::abs(std::pow(resultY[j], 2));
108 if (test > DBL_MAX)
109 test = 0;
110 error2_factor2 += errors[j] * test / factor2 / factor2;
111 }
112 double error2_factor = (error2_factor1 / factor1 / factor1 + error2_factor2);
113
114 // Calculate the normalized Y values
115 using std::placeholders::_1;
116 // NOTE: Previously, we had been using std::transform with
117 // std::bind(std::multiplies<double>(), _1,factor)
118 // here, but that seemed to have strange effects in Windows Debug
119 // builds which caused the unit tests
120 // to sometimes fail. Maybe this is some compiler bug to do with
121 // using bind2nd within the parrallel macros.
122
123 // Now result is s_i/v_i*Dlam_i*(sum_i s_i)/(sum_i S_i/v_i*Dlam_i)
124 std::transform(resultY.begin(), resultY.end(), resultY.begin(), [&factor](double rY) { return rY * factor; });
125
126 // Finally get the normalized errors
127 for (int j = 0; j < size - 1; j++)
128 resultE[j] = resultY[j] * sqrt(errors[j] + error2_factor);
129
130 // Check that any masking matches, print a warning if not
131 check_masks(inputWS1, inputWS2, i);
132
133 prog.report();
135 }
137
138 outputWS->setYUnitLabel("Counts normalised to a vanadium");
139 outputWS->setDistribution(false);
140}
141
150 // First check that the instrument matches for both input workspaces
151 if (w1->getInstrument()->getName() != w2->getInstrument()->getName()) {
152 g_log.error("The input workspaces have different instrument definitions");
153 throw std::runtime_error("The input workspaces have different instrument definitions");
154 }
155 // Check that the two workspaces are the same size
156 if (w1->size() != w2->size()) {
157 g_log.error("The input workspaces are not the same size");
158 throw std::runtime_error("The input workspaces are not the same size");
159 }
160 // Now check that the bins match
161 if (!WorkspaceHelpers::matchingBins(*w1, *w2)) {
162 g_log.error("The input workspaces have different binning");
163 throw std::runtime_error("The input workspaces have different binning");
164 }
165 const Mantid::API::Axis *const axis1 = w1->getAxis(1);
166 const Mantid::API::Axis *const axis2 = w2->getAxis(1);
167 if (!((*axis1) == (*axis2))) // Spectra axis are different, so division does
168 // not make any sense
169 {
170 g_log.error("The two workspaces InputW1 and InputW2 have different spectra list");
171 throw std::runtime_error("The two workspaces InputW1 and InputW2 have different spectra list");
172 }
173
174 if (out != w1 && out != w2) // Create a new workspace only if it is different
175 // from of the input ones.
176 {
177 out = create<MatrixWorkspace>(*w1);
178 setProperty("OutputWorkspace", out);
179 } else if (out == w2) {
180 g_log.warning("Any masking in the output workspaces will be taken from the "
181 "vanadium workspace (InputW2)");
182 }
183}
184
192 const API::MatrixWorkspace_const_sptr &w2, const int &index) const {
193 static bool warned = false;
194 if (!warned) {
195 const bool w1masked = w1->hasMaskedBins(index);
196 const bool w2masked = w2->hasMaskedBins(index);
197 if ((w1masked && w2masked && (w1->maskedBins(index) != w2->maskedBins(index))) || (w1masked && !w2masked) ||
198 (!w1masked && w2masked)) {
199 g_log.warning("The input workspaces do not have matching bin masking");
200 warned = true;
201 }
202 }
203}
204
205} // 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.
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
Class to represent the axis of a workspace.
Definition: Axis.h:30
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void check_masks(const API::MatrixWorkspace_const_sptr &w1, const API::MatrixWorkspace_const_sptr &w2, const int &index) const
Checks whether the two input workspaces have the same bin masking.
void check_validity(API::MatrixWorkspace_const_sptr &w1, API::MatrixWorkspace_const_sptr &w2, API::MatrixWorkspace_sptr &out)
Checks that the axes of the input workspaces match and creates the output workspace if necessary.
void exec() override
Virtual method - must be overridden by concrete algorithm.
void init() override
Virtual method - must be overridden by concrete algorithm.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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 report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
std::complex< double > MANTID_API_DLL E1(std::complex< double > z)
Integral for Gamma.
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::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
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
static bool matchingBins(const MatrixWorkspace &ws1, const MatrixWorkspace &ws2, const bool firstOnly=false)
Checks whether the bins (X values) of two workspace are the same.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54
Divide functor with result reset to 0 if the denominator is null.
Definition: VectorHelper.h:210