Commit a669155d authored by Peterson, Peter's avatar Peterson, Peter
Browse files

Added more const and some minor renames

parent 7e787421
......@@ -68,8 +68,6 @@ private:
/// Helper functions
void findCenterOfMass(API::MatrixWorkspace_sptr inputWS, double &centerX, double &centerY, const int numSpec,
API::Progress &progress);
API::MatrixWorkspace_sptr sumUsingSpectra(DataObjects::EventWorkspace_const_sptr inputEventWS, const int numSpec,
API::Progress &progress);
void storeOutputWorkspace(double centerX, double centerY);
// Iteration cutoff
const int m_maxIteration = 200;
......
......@@ -24,54 +24,54 @@ Kernel::Logger g_log("WorkspaceBoundingBox");
*/
class WorkspaceBoundingBox {
public:
WorkspaceBoundingBox(API::MatrixWorkspace_sptr workspace);
WorkspaceBoundingBox(API::MatrixWorkspace_const_sptr workspace);
WorkspaceBoundingBox();
~WorkspaceBoundingBox();
API::MatrixWorkspace_sptr getWorkspace() { return workspace; }
double getX() { return x; }
double getY() { return y; }
double getCenterX() { return centerX; }
double getCenterY() { return centerY; }
double getXMin() { return xMin; }
double getXMax() { return xMax; }
double getYMin() { return yMin; }
double getYMax() { return yMax; }
API::MatrixWorkspace_const_sptr getWorkspace() { return m_workspace; }
double getX() const { return x; }
double getY() const { return y; }
double getCenterX() const { return centerX; }
double getCenterY() const { return centerY; }
double getXMin() const { return xMin; }
double getXMax() const { return xMax; }
double getYMin() const { return yMin; }
double getYMax() const { return yMax; }
void setPosition(double x, double y);
void setCenter(double x, double y);
void setBounds(double xMin, double xMax, double yMin, double yMax);
double calculateDistance();
double calculateRadiusX();
double calculateRadiusY();
double calculateDistance() const;
double calculateRadiusX() const;
double calculateRadiusY() const;
double updatePositionAndReturnCount(int index);
int findFirstValidWs(const int numSpec);
bool isValidWs(int index);
int findFirstValidWs(const int numSpec) const;
bool isValidWs(int index) const;
bool isOutOfBoundsOfNonDirectBeam(const double beamRadius, int index, const bool directBeam);
bool containsPoint(double x, double y);
void normalizePosition(double x, double y);
void updateMinMax(int index);
private:
Kernel::V3D &position(int index);
const HistogramData::HistogramY &histogramY(int index);
API::MatrixWorkspace_sptr workspace;
const API::SpectrumInfo *spectrumInfo;
double x = 0;
double y = 0;
double centerX;
double centerY;
double xMin = 0;
double xMax = 0;
double yMin = 0;
double yMax = 0;
const int m_specID = 0;
int m_cachedPositionIndex = -1;
int m_cachedHistogramYIndex = -1;
Kernel::V3D m_cachedPosition;
const HistogramData::HistogramY *m_cachedHistogramY;
Kernel::V3D &position(int index) const;
double yValue(const int index) const;
API::MatrixWorkspace_const_sptr m_workspace;
const API::SpectrumInfo *m_spectrumInfo;
double x{0};
double y{0};
double centerX{0};
double centerY{0};
double xMin{0};
double xMax{0};
double yMin{0};
double yMax{0};
// cache information
mutable int m_cachedPositionIndex{-1};
mutable Kernel::V3D m_cachedPosition;
mutable int m_cachedHistogramYIndex{-1};
mutable double m_cachedYValue;
};
} // namespace Algorithms
......
......@@ -31,6 +31,12 @@ using namespace API;
using namespace Geometry;
using namespace DataObjects;
namespace {
bool equals(const double a, const double b) { // anonymous namespace
return fabs(a - b) < std::numeric_limits<double>::min();
} // anonymous namespace
} // namespace
void FindCenterOfMassPosition2::init() {
const auto wsValidator = std::make_shared<CompositeValidator>();
const auto positiveDouble = std::make_shared<BoundedValidator<double>>();
......@@ -60,40 +66,6 @@ void FindCenterOfMassPosition2::init() {
"the center of mass of the scattering pattern.");
}
/** Integrate events in the inputEventWS to determine thier sum and error values
*
* @param inputEventWS :: event workspace to ingegrate
* @param numSpec :: number of events in the workspace to iterate through
* @param progress :: object for reporting progress of the operation
* @return workspace containing the calculated x, y, and e data over numSpec
*/
API::MatrixWorkspace_sptr
FindCenterOfMassPosition2::sumUsingSpectra(DataObjects::EventWorkspace_const_sptr inputEventWS, const int numSpec,
Progress &progress) {
std::vector<double> yValues(numSpec);
std::vector<double> eValues(numSpec);
PARALLEL_FOR_NO_WSP_CHECK()
for (int i = 0; i < numSpec; i++) {
double sum_i(0), err_i(0);
progress.report("Integrating events");
const EventList &el = inputEventWS->getSpectrum(i);
el.integrate(0, 0, true, sum_i, err_i);
yValues[i] = sum_i;
eValues[i] = err_i;
}
IAlgorithm_sptr algo = createChildAlgorithm("CreateWorkspace", 0.7, 1.0);
algo->setProperty<std::vector<double>>("DataX", std::vector<double>(2, 0.0));
algo->setProperty<std::vector<double>>("DataY", yValues);
algo->setProperty<std::vector<double>>("DataE", eValues);
algo->setProperty<int>("NSpec", numSpec);
algo->execute();
MatrixWorkspace_sptr inputWS = algo->getProperty("OutputWorkspace");
return inputWS;
}
// find the min/max for x/y coords in set of valid spectra, update position of bounding box
double initBoundingBox(WorkspaceBoundingBox &boundingBox, const int numSpec, const double beamRadius,
const bool directBeam) {
......@@ -134,8 +106,6 @@ double updateBoundingBox(WorkspaceBoundingBox &boundingBox, WorkspaceBoundingBox
return totalCount;
}
bool equals(double a, double b) { return fabs(a - b) < std::numeric_limits<double>::min(); }
/** Iterates through spectrum in the input workspace finding the center of mass until
* we converge to within the tolerance specified in meters
*
......@@ -188,8 +158,10 @@ void FindCenterOfMassPosition2::findCenterOfMass(API::MatrixWorkspace_sptr input
}
boundingBox.setCenter(boundingBox.getX(), boundingBox.getY());
previousBoundingBox.setBounds(boundingBox.getCenterX() - radiusX, boundingBox.getCenterX() + radiusX,
boundingBox.getCenterY() - radiusY, boundingBox.getCenterY() + radiusY);
const auto oldCenterX = boundingBox.getCenterX();
const auto oldCenterY = boundingBox.getCenterY();
previousBoundingBox.setBounds(oldCenterX - radiusX, oldCenterX + radiusX, oldCenterY - radiusY,
oldCenterY + radiusY);
// Check to see if we have the same result
// as the previous iteration
......@@ -275,23 +247,18 @@ void FindCenterOfMassPosition2::exec() {
MatrixWorkspace_sptr inputWS;
// Sum up all the wavelength bins
IAlgorithm_sptr childAlg = createChildAlgorithm("Integration", 0., 0.5); // first half is integrating
childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWSWvl);
childAlg->executeAsChildAlg();
inputWS = childAlg->getProperty("OutputWorkspace");
// Get the number of monitors. We assume that all monitors are stored in the
// first spectra
const auto numSpec = static_cast<int>(inputWSWvl->getNumberHistograms());
// Set up the progress reporting object
Progress progress(this, 0.0, 1.0, m_maxIteration);
EventWorkspace_const_sptr inputEventWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWSWvl);
if (inputEventWS) {
inputWS = sumUsingSpectra(inputEventWS, numSpec, progress);
WorkspaceFactory::Instance().initializeFromParent(*inputWSWvl, *inputWS, false);
} else {
// Sum up all the wavelength bins
IAlgorithm_sptr childAlg = createChildAlgorithm("Integration");
childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWSWvl);
childAlg->executeAsChildAlg();
inputWS = childAlg->getProperty("OutputWorkspace");
}
// Set up the progress reporting object
Progress progress(this, 0.5, 1.0, m_maxIteration);
findCenterOfMass(inputWS, centerX, centerY, numSpec, progress);
storeOutputWorkspace(centerX, centerY);
......
......@@ -4,30 +4,40 @@
namespace Mantid {
namespace Algorithms {
WorkspaceBoundingBox::WorkspaceBoundingBox(API::MatrixWorkspace_sptr workspace) {
this->workspace = workspace;
this->spectrumInfo = &workspace->spectrumInfo();
namespace {
constexpr int HISTOGRAM_INDEX{0};
}
WorkspaceBoundingBox::WorkspaceBoundingBox() {}
WorkspaceBoundingBox::WorkspaceBoundingBox(API::MatrixWorkspace_const_sptr workspace) : m_workspace(workspace) {
if (m_workspace->y(0).size() != 1)
throw std::runtime_error("This object only works with integrated workspaces");
m_spectrumInfo = &workspace->spectrumInfo();
}
WorkspaceBoundingBox::WorkspaceBoundingBox() {
m_spectrumInfo = nullptr; // certain functionality is not available
}
WorkspaceBoundingBox::~WorkspaceBoundingBox() {}
Kernel::V3D &WorkspaceBoundingBox::position(int index) {
if (m_cachedPositionIndex == -1 || m_cachedPositionIndex != index) {
m_cachedPosition = this->spectrumInfo->position(index);
Kernel::V3D &WorkspaceBoundingBox::position(int index) const {
if (m_cachedPositionIndex != index) {
if (!m_spectrumInfo)
throw std::runtime_error("SpectrumInfo object is not initialized");
m_cachedPosition = m_spectrumInfo->position(index);
m_cachedPositionIndex = index;
}
return m_cachedPosition;
}
const HistogramData::HistogramY &WorkspaceBoundingBox::histogramY(int index) {
if (m_cachedHistogramYIndex == -1 || m_cachedHistogramYIndex != index) {
auto &YIn = this->workspace->y(index);
m_cachedHistogramY = &YIn;
double WorkspaceBoundingBox::yValue(const int index) const {
if (m_cachedHistogramYIndex != index) {
m_cachedYValue = m_workspace->y(index)[HISTOGRAM_INDEX];
m_cachedHistogramYIndex = index;
}
return *m_cachedHistogramY;
return m_cachedYValue;
}
void WorkspaceBoundingBox::setPosition(double x, double y) {
......@@ -53,19 +63,21 @@ void WorkspaceBoundingBox::setBounds(double xMin, double xMax, double yMin, doub
* @param index :: index of spectrum data
* @return true/false if its valid
*/
bool WorkspaceBoundingBox::isValidWs(int index) {
if (!this->spectrumInfo->hasDetectors(index)) {
bool WorkspaceBoundingBox::isValidWs(int index) const {
if (!m_spectrumInfo)
throw std::runtime_error("SpectrumInfo object is not initialized");
if (!m_spectrumInfo->hasDetectors(index)) {
g_log.warning() << "Workspace index " << index << " has no detector assigned to it - discarding\n";
return false;
}
// Skip if we have a monitor or if the detector is masked.
if (this->spectrumInfo->isMonitor(index) || this->spectrumInfo->isMasked(index))
if (this->m_spectrumInfo->isMonitor(index) || this->m_spectrumInfo->isMasked(index))
return false;
// Get the current spectrum
auto &YIn = this->histogramY(index);
const auto YIn = this->yValue(index);
// Skip if NaN of inf
if (std::isnan(YIn[m_specID]) || std::isinf(YIn[m_specID]))
if (std::isnan(YIn) || std::isinf(YIn))
return false;
return true;
}
......@@ -75,7 +87,7 @@ bool WorkspaceBoundingBox::isValidWs(int index) {
* @param numSpec :: the number of spectrum in the workspace to search through
* @return index of first valid spectrum
*/
int WorkspaceBoundingBox::findFirstValidWs(const int numSpec) {
int WorkspaceBoundingBox::findFirstValidWs(const int numSpec) const {
int i;
for (i = 0; i < numSpec; ++i) {
if (isValidWs(i))
......@@ -91,12 +103,13 @@ int WorkspaceBoundingBox::findFirstValidWs(const int numSpec) {
* @return number of points of histogram data at index
*/
double WorkspaceBoundingBox::updatePositionAndReturnCount(int index) {
auto &YIn = this->histogramY(index);
double x = this->position(index).X();
double y = this->position(index).Y();
this->x += YIn[m_specID] * x;
this->y += YIn[m_specID] * y;
return YIn[m_specID];
const auto YIn = this->yValue(index);
const auto &position = this->position(index);
this->x += YIn * position.X();
this->y += YIn * position.Y();
return YIn;
}
/** Compare current mins and maxs to the coordinates of the spectrum at index
......@@ -105,8 +118,10 @@ double WorkspaceBoundingBox::updatePositionAndReturnCount(int index) {
* @param index :: index of spectrum data
*/
void WorkspaceBoundingBox::updateMinMax(int index) {
double x = this->position(index).X();
double y = this->position(index).Y();
const auto &position = this->position(index);
const double x = position.X();
const double y = position.Y();
this->xMin = std::min(x, this->xMin);
this->xMax = std::max(x, this->xMax);
this->yMin = std::min(y, this->yMin);
......@@ -121,22 +136,25 @@ void WorkspaceBoundingBox::updateMinMax(int index) {
* @return number of points of histogram data at index
*/
bool WorkspaceBoundingBox::isOutOfBoundsOfNonDirectBeam(const double beamRadius, int index, const bool directBeam) {
double x = this->position(index).X();
double y = this->position(index).Y();
if (!directBeam) {
double dx = x - this->centerX;
double dy = y - this->centerY;
const auto &position = this->position(index);
const double dx = position.X() - this->centerX;
const double dy = position.Y() - this->centerY;
if (dx * dx + dy * dy < beamRadius * beamRadius)
return false;
}
return true;
}
double WorkspaceBoundingBox::calculateDistance() {
return sqrt((centerX - x) * (centerX - x) + (centerY - y) * (centerY - y));
double WorkspaceBoundingBox::calculateDistance() const {
const auto xExtent = (centerX - x);
const auto yExtent = (centerY - y);
return sqrt(xExtent * xExtent + yExtent * yExtent);
}
double WorkspaceBoundingBox::calculateRadiusX() { return std::min((x - xMin), (xMax - x)); }
double WorkspaceBoundingBox::calculateRadiusY() { return std::min((y - yMin), (yMax - y)); }
double WorkspaceBoundingBox::calculateRadiusX() const { return std::min((x - xMin), (xMax - x)); }
double WorkspaceBoundingBox::calculateRadiusY() const { return std::min((y - yMin), (yMax - y)); }
/** Perform normalization on x/y coords over given values
*
......@@ -144,8 +162,8 @@ double WorkspaceBoundingBox::calculateRadiusY() { return std::min((y - yMin), (y
* @param y :: value to normalize member y over
*/
void WorkspaceBoundingBox::normalizePosition(double x, double y) {
this->x /= x;
this->y /= y;
this->x /= std::fabs(x);
this->y /= std::fabs(y);
}
/** Checks if a given x/y coord is within the bounding box
......@@ -155,10 +173,8 @@ void WorkspaceBoundingBox::normalizePosition(double x, double y) {
* @return true/false if it is within the mins/maxs of the box
*/
bool WorkspaceBoundingBox::containsPoint(double x, double y) {
if (x > this->xMax || x < this->xMin || y > yMax || y < yMin)
return false;
return true;
return (x <= this->xMax && x >= this->xMin && y <= yMax && y >= yMin);
}
} // namespace Algorithms
} // namespace Mantid
\ No newline at end of file
} // namespace Mantid
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment