Mantid
Loading...
Searching...
No Matches
Rebin2D.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 +
21
22namespace Mantid::Algorithms {
23
24// Register the algorithm into the AlgorithmFactory
25DECLARE_ALGORITHM(Rebin2D)
26
27using namespace API;
28using namespace DataObjects;
29using namespace Geometry;
30using namespace Mantid::HistogramData;
31
32//--------------------------------------------------------------------------
33// Private methods
34//--------------------------------------------------------------------------
44 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input), "An input workspace.");
45 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
46 "An output workspace.");
47 const std::string docString = "A comma separated list of first bin boundary, width, last bin boundary. "
48 "Optionally "
49 "this can be followed by a comma and more widths and last boundary "
50 "pairs. "
51 "Negative width values indicate logarithmic binning.";
52 auto rebinValidator = std::make_shared<RebinParamsValidator>();
53 declareProperty(std::make_unique<ArrayProperty<double>>("Axis1Binning", rebinValidator), docString);
54 declareProperty(std::make_unique<ArrayProperty<double>>("Axis2Binning", rebinValidator), docString);
55 declareProperty(std::make_unique<PropertyWithValue<bool>>("UseFractionalArea", false, Direction::Input),
56 "Flag to turn on the using the fractional area tracking RebinnedOutput "
57 "workspace\n."
58 "Default is false.");
59 declareProperty(std::make_unique<PropertyWithValue<bool>>("Transpose", false),
60 "Run the Transpose algorithm on the resulting matrix.");
61}
62
67 // Information to form input grid
68 MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
69 const NumericAxis *oldAxis2 = dynamic_cast<API::NumericAxis *>(inputWS->getAxis(1));
70 if (!oldAxis2) {
71 throw std::invalid_argument("Vertical axis is not a numeric axis, cannot rebin. "
72 "If it is a spectra axis try running ConvertSpectrumAxis first.");
73 }
74
75 const auto &oldXEdges = inputWS->binEdges(0);
76 const size_t numXBins = inputWS->blocksize();
77 const size_t numYBins = inputWS->getNumberHistograms();
78 // This will convert plain NumericAxis to bin edges while
79 // BinEdgeAxis will just return its edges.
80 const std::vector<double> oldYEdges = oldAxis2->createBinBoundaries();
81
82 // Output grid and workspace. Fills in the new X and Y bin vectors
83 // MantidVecPtr newXBins;
84 BinEdges newXBins(oldXEdges.size());
85 BinEdges newYBins(oldXEdges.size());
86
87 // Flag for using a RebinnedOutput workspace
88 // NB. This is now redundant because if the input is a MatrixWorkspace,
89 // useFractionArea=false is forced since there is no fractional area info.
90 // But if the input is RebinnedOutput, useFractionalArea=true is forced to
91 // give correct signal/errors. It is kept for compatibility with old scripts.
92 bool useFractionalArea = getProperty("UseFractionalArea");
93 auto inputHasFA = std::dynamic_pointer_cast<const RebinnedOutput>(inputWS);
94 // For MatrixWorkspace, only UseFractionalArea=False makes sense.
95 if (useFractionalArea && !inputHasFA) {
96 g_log.warning("Fractional area tracking was requested but input workspace does "
97 "not have calculated bin fractions. Assuming bins are exact "
98 "(fractions are unity). The results may not be accurate if this "
99 "workspace was previously rebinned.");
100 }
101 // For RebinnedOutput, should always use useFractionalArea to get the
102 // correct signal and errors (so that weights of input ws is accounted for).
103 if (inputHasFA && !useFractionalArea) {
104 g_log.warning("Input workspace has bin fractions (e.g. from a "
105 "parallelpiped rebin like SofQW3). To give accurate results, "
106 "fractional area tracking has been turn on.");
107 useFractionalArea = true;
108 }
109
110 auto outputWS = createOutputWorkspace(inputWS, newXBins, newYBins, useFractionalArea);
111 auto outputRB = std::dynamic_pointer_cast<RebinnedOutput>(outputWS);
112
113 // Progress reports & cancellation
114 const auto nreports(static_cast<size_t>(numYBins));
115 m_progress = std::make_unique<API::Progress>(this, 0.0, 1.0, nreports);
116
117 PARALLEL_FOR_IF(Kernel::threadSafe(*inputWS, *outputWS))
118 for (int64_t i = 0; i < static_cast<int64_t>(numYBins); ++i) {
120
121 m_progress->report("Computing polygon intersections");
122 const double vlo = oldYEdges[i];
123 const double vhi = oldYEdges[i + 1];
124
125 for (size_t j = 0; j < numXBins; ++j) {
126 // For each input polygon test where it intersects with
127 // the output grid and assign the appropriate weights of Y/E
128 const double x_j = oldXEdges[j];
129 const double x_jp1 = oldXEdges[j + 1];
130 Quadrilateral inputQ(x_j, x_jp1, vlo, vhi);
131 if (!useFractionalArea) {
132 FractionalRebinning::rebinToOutput(inputQ, inputWS, i, j, *outputWS, newYBins.rawData());
133 } else {
134 FractionalRebinning::rebinToFractionalOutput(inputQ, inputWS, i, j, *outputRB, newYBins.rawData(), inputHasFA);
135 }
136 }
137
139 }
141 if (useFractionalArea) {
143 outputRB->finalize(true);
144 }
145
146 FractionalRebinning::normaliseOutput(outputWS, inputWS, m_progress.get());
147
148 bool Transpose = this->getProperty("Transpose");
149 if (Transpose) {
150 auto alg = createChildAlgorithm("Transpose", 0.9, 1.0);
151 alg->setProperty("InputWorkspace", outputWS);
152 alg->setPropertyValue("OutputWorkspace", "__anonymous");
153 alg->execute();
154 outputWS = alg->getProperty("OutputWorkspace");
155 }
156
157 setProperty("OutputWorkspace", outputWS);
158}
159
171 BinEdges &newYBins, const bool useFractionalArea) const {
173
174 auto &newY = newYBins.mutableRawData();
175 // First create the two sets of bin boundaries
176 static_cast<void>(createAxisFromRebinParams(getProperty("Axis1Binning"), newXBins.mutableRawData()));
177 const int newYSize = createAxisFromRebinParams(getProperty("Axis2Binning"), newY);
178 // and now the workspace
179 HistogramData::BinEdges binEdges(newXBins);
180 Workspace2D_sptr outputWS;
181 if (!useFractionalArea) {
182 outputWS = create<Workspace2D>(*parent, newYSize - 1, binEdges);
183 } else {
184 outputWS = create<RebinnedOutput>(*parent, newYSize - 1, binEdges);
185 }
186 auto verticalAxis = std::make_unique<BinEdgeAxis>(newY);
187 // Meta data
188 verticalAxis->unit() = parent->getAxis(1)->unit();
189 verticalAxis->title() = parent->getAxis(1)->title();
190 outputWS->replaceAxis(1, std::move(verticalAxis));
191
192 return outputWS;
193}
194
195} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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
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
Class to represent a numeric axis of a workspace.
Definition: NumericAxis.h:29
virtual std::vector< double > createBinBoundaries() const
Create bin boundaries from the point values.
A property class for workspaces.
void exec() override
Run the algorithm.
Definition: Rebin2D.cpp:66
std::unique_ptr< API::Progress > m_progress
Progress reporter.
Definition: Rebin2D.h:41
void init() override
Initialise the properties.
Definition: Rebin2D.cpp:38
API::MatrixWorkspace_sptr createOutputWorkspace(const API::MatrixWorkspace_const_sptr &parent, HistogramData::BinEdges &newXBins, HistogramData::BinEdges &newYBins, const bool useFractionalArea) const
Setup the output workspace.
Definition: Rebin2D.cpp:170
This algorithm "transposes" the bins of the input workspace into a single spectra.
Definition: Transpose.h:33
A ConvexPolygon with only 4 vertices.
Definition: Quadrilateral.h:24
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
The concrete, templated class for properties.
Validator to check the format of a vector providing the rebin parameters to an algorithm.
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
MANTID_DATAOBJECTS_DLL void normaliseOutput(const API::MatrixWorkspace_sptr &outputWS, const API::MatrixWorkspace_const_sptr &inputWS, API::Progress *progress=nullptr)
Compute sqrt of errors and put back in bin width division if necessary.
MANTID_DATAOBJECTS_DLL void rebinToFractionalOutput(const Geometry::Quadrilateral &inputQ, const API::MatrixWorkspace_const_sptr &inputWS, const size_t i, const size_t j, DataObjects::RebinnedOutput &outputWS, const std::vector< double > &verticalAxis, const DataObjects::RebinnedOutput_const_sptr &inputRB=nullptr)
Rebin the input quadrilateral to to output grid.
MANTID_DATAOBJECTS_DLL void rebinToOutput(const Geometry::Quadrilateral &inputQ, const API::MatrixWorkspace_const_sptr &inputWS, const size_t i, const size_t j, API::MatrixWorkspace &outputWS, const std::vector< double > &verticalAxis)
Rebin the input quadrilateral to to output grid.
MANTID_DATAOBJECTS_DLL void finalizeFractionalRebin(DataObjects::RebinnedOutput &outputWS)
Set finalize flag after fractional rebinning loop.
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
int MANTID_KERNEL_DLL createAxisFromRebinParams(const std::vector< double > &params, std::vector< double > &xnew, const bool resize_xnew=true, const bool full_bins_only=false, const double xMinHint=std::nan(""), const double xMaxHint=std::nan(""), const bool useReverseLogarithmic=false, const double power=-1)
Creates a new output X array given a 'standard' set of rebinning parameters.
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
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54