Skip to content

Commit

Permalink
Re #11619. RemoveMaskedSpectra and CropWorkspace use ExtractSpectra.
Browse files Browse the repository at this point in the history
  • Loading branch information
mantid-roman committed Apr 24, 2015
1 parent cfd9484 commit 2e0326d
Show file tree
Hide file tree
Showing 6 changed files with 91 additions and 119 deletions.
Expand Up @@ -46,6 +46,7 @@ class DLLExport ExtractSpectra : public API::Algorithm {
private:
void init();
void exec();
void execHistogram();
void execEvent();

void checkProperties();
Expand All @@ -61,16 +62,14 @@ class DLLExport ExtractSpectra : public API::Algorithm {
std::size_t m_minX;
/// The bin index to end the cropped workspace at
std::size_t m_maxX;
/// The spectrum index to start the cropped workspace from
specid_t m_minSpec;
/// The spectrum index to end the cropped workspace at
specid_t m_maxSpec;
/// Flag indicating whether the input workspace has common boundaries
bool m_commonBoundaries;
/// Flag indicating whether we're dealing with histogram data
bool m_histogram;
/// Flag indicating whether XMin and/or XMax has been set
bool m_croppingInX;
/// The list of spectra to extract.
std::vector<specid_t> m_spectrumList;
};

} // namespace Algorithms
Expand Down
Expand Up @@ -43,7 +43,7 @@ class DLLExport RemoveMaskedSpectra : public API::Algorithm {
private:
void init();
void exec();
void makeIndexList(std::vector<size_t> &indices,
void makeIndexList(std::vector<specid_t> &indices,
const API::MatrixWorkspace *maskedWorkspace);
};

Expand Down
17 changes: 16 additions & 1 deletion Code/Mantid/Framework/Algorithms/src/CropWorkspace.cpp
Expand Up @@ -71,8 +71,23 @@ void CropWorkspace::exec() {

auto extract = createChildAlgorithm("ExtractSpectra");
extract->initialize();
extract->copyPropertiesFrom(*this);
extract->setRethrows(true);

MatrixWorkspace_sptr inputWorkspace = getProperty("InputWorkspace");
extract->setProperty("InputWorkspace", inputWorkspace);

double xmin = getProperty("XMin");
extract->setProperty("XMin", xmin);

double xmax = getProperty("XMax");
extract->setProperty("XMax", xmax);

int start = getProperty("StartWorkspaceIndex");
extract->setProperty("StartWorkspaceIndex", start);

int end = getProperty("EndWorkspaceIndex");
extract->setProperty("EndWorkspaceIndex", end);

extract->execute();

MatrixWorkspace_sptr outputWorkspace = extract->getProperty("OutputWorkspace");
Expand Down
113 changes: 53 additions & 60 deletions Code/Mantid/Framework/Algorithms/src/ExtractSpectra.cpp
Expand Up @@ -4,6 +4,7 @@
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/TextAxis.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/VectorHelper.h"

Expand All @@ -28,7 +29,7 @@ DECLARE_ALGORITHM(ExtractSpectra)
/** Constructor
*/
ExtractSpectra::ExtractSpectra()
: Algorithm(), m_minX(0), m_maxX(0), m_minSpec(-1), m_maxSpec(-1),
: Algorithm(), m_minX(0), m_maxX(0), //m_minSpec(-1), m_maxSpec(-1),
m_commonBoundaries(false), m_histogram(false), m_croppingInX(false) {}

//----------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -80,12 +81,15 @@ void ExtractSpectra::init() {
"will be loaded\n"
"(default: first entry in the Workspace)");
// As the property takes ownership of the validator pointer, have to take care
// to pass in a unique
// pointer to each property.
// to pass in a unique pointer to each property.
declareProperty(
"EndWorkspaceIndex", EMPTY_INT(), mustBePositive,
"The index number of the last entry in the Workspace to be loaded\n"
"(default: last entry in the Workspace)");
declareProperty(
new ArrayProperty<specid_t>("SpectrumList"),
"A comma-separated list of individual spectra to read. Only used if\n"
"explicitly set.");
}

//----------------------------------------------------------------------------------------------
Expand All @@ -101,9 +105,14 @@ void ExtractSpectra::exec() {
if (eventW != NULL) {
// Input workspace is an event workspace. Use the other exec method
this->execEvent();
return;
} else {
// Otherwise it's a Workspace2D
this->execHistogram();
}
}

/// Execute the algorithm in case of a histogrammed data.
void ExtractSpectra::execHistogram() {
m_histogram = m_inputWorkspace->isHistogramData();
// Check for common boundaries in input workspace
m_commonBoundaries = WorkspaceHelpers::commonBoundaries(m_inputWorkspace);
Expand All @@ -113,7 +122,7 @@ void ExtractSpectra::exec() {

// Create the output workspace
MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create(
m_inputWorkspace, m_maxSpec - m_minSpec + 1, m_maxX - m_minX,
m_inputWorkspace, m_spectrumList.size(), m_maxX - m_minX,
m_maxX - m_minX - m_histogram);

// If this is a Workspace2D, get the spectra axes for copying in the spectraNo
Expand All @@ -128,12 +137,13 @@ void ExtractSpectra::exec() {

cow_ptr<MantidVec> newX;
if (m_commonBoundaries) {
const MantidVec &oldX = m_inputWorkspace->readX(m_minSpec);
const MantidVec &oldX = m_inputWorkspace->readX(m_spectrumList.front());
newX.access().assign(oldX.begin() + m_minX, oldX.begin() + m_maxX);
}
Progress prog(this, 0.0, 1.0, (m_maxSpec - m_minSpec));
Progress prog(this, 0.0, 1.0, (m_spectrumList.size()));
// Loop over the required spectra, copying in the desired bins
for (int i = m_minSpec, j = 0; i <= m_maxSpec; ++i, ++j) {
for (int j = 0; j < m_spectrumList.size(); ++j) {
auto i = m_spectrumList[j];
// Preserve/restore sharing if X vectors are the same
if (m_commonBoundaries) {
outputWorkspace->setX(j, newX);
Expand Down Expand Up @@ -185,31 +195,6 @@ void ExtractSpectra::exec() {
setProperty("OutputWorkspace", outputWorkspace);
}

template <typename T>
std::size_t lowerBound(const std::vector<T> &events, const double tof) {
typename std::vector<T>::const_iterator first = events.begin();
typename std::vector<T>::const_iterator last = events.end();
typename std::vector<T>::const_iterator it;
typename std::vector<T>::difference_type count = distance(first, last);
typename std::vector<T>::difference_type step;

while (count > 0) {
it = first;
step = count / 2;
advance(it, step);
if (it->tof() < tof) // or: if (comp(*it,value)), for the comp version
{
first = ++it;
count -= step + 1;
} else
count = step;
}
if (first == events.end())
return events.size();
else
return distance(events.begin(), first);
}

/** Executes the algorithm
* @throw std::out_of_range If a property is set to an invalid value for the
* input workspace
Expand All @@ -230,7 +215,7 @@ void ExtractSpectra::execEvent() {
this->checkProperties();
cow_ptr<MantidVec> XValues_new;
if (m_commonBoundaries) {
const MantidVec &oldX = m_inputWorkspace->readX(m_minSpec);
const MantidVec &oldX = m_inputWorkspace->readX(m_spectrumList.front());
XValues_new.access().assign(oldX.begin() + m_minX, oldX.begin() + m_maxX);
}
size_t ntcnew = m_maxX - m_minX;
Expand All @@ -255,21 +240,21 @@ void ExtractSpectra::execEvent() {
EventWorkspace_sptr outputWorkspace =
boost::dynamic_pointer_cast<EventWorkspace>(
API::WorkspaceFactory::Instance().create(
"EventWorkspace", m_maxSpec - m_minSpec + 1, ntcnew,
"EventWorkspace", m_spectrumList.size(), ntcnew,
ntcnew - m_histogram));
eventW->sortAll(TOF_SORT, NULL);
outputWorkspace->sortAll(TOF_SORT, NULL);
// Copy required stuff from it
API::WorkspaceFactory::Instance().initializeFromParent(m_inputWorkspace,
outputWorkspace, true);

Progress prog(this, 0.0, 1.0, 2 * (m_maxSpec - m_minSpec));
Progress prog(this, 0.0, 1.0, 2 * m_spectrumList.size());
eventW->sortAll(Mantid::DataObjects::TOF_SORT, &prog);
// Loop over the required spectra, copying in the desired bins
PARALLEL_FOR2(m_inputWorkspace, outputWorkspace)
for (int i = m_minSpec; i <= m_maxSpec; ++i) {
for (int j = 0; j < m_spectrumList.size(); ++j) {
PARALLEL_START_INTERUPT_REGION
int j = i - m_minSpec;
auto i = m_spectrumList[j];
const EventList &el = eventW->getEventList(i);
// The output event list
EventList &outEL = outputWorkspace->getOrAddEventList(j);
Expand Down Expand Up @@ -390,27 +375,35 @@ void ExtractSpectra::checkProperties() {
if (!m_commonBoundaries)
m_maxX = static_cast<int>(m_inputWorkspace->readX(0).size());

m_minSpec = getProperty("StartWorkspaceIndex");
const int numberOfSpectra =
static_cast<int>(m_inputWorkspace->getNumberHistograms());
m_maxSpec = getProperty("EndWorkspaceIndex");
if (isEmpty(m_maxSpec))
m_maxSpec = numberOfSpectra - 1;

// Check 'StartSpectrum' is in range 0-numberOfSpectra
if (m_minSpec > numberOfSpectra - 1) {
g_log.error("StartWorkspaceIndex out of range!");
throw std::out_of_range("StartSpectrum out of range!");
}
if (m_maxSpec > numberOfSpectra - 1) {
g_log.error("EndWorkspaceIndex out of range!");
throw std::out_of_range("EndWorkspaceIndex out of range!");
}
if (m_maxSpec < m_minSpec) {
g_log.error(
"StartWorkspaceIndex must be less than or equal to EndWorkspaceIndex");
throw std::out_of_range(
"StartWorkspaceIndex must be less than or equal to EndWorkspaceIndex");
m_spectrumList = getProperty("SpectrumList");

if (m_spectrumList.empty()) {
int minSpec = getProperty("StartWorkspaceIndex");
const int numberOfSpectra =
static_cast<int>(m_inputWorkspace->getNumberHistograms());
int maxSpec = getProperty("EndWorkspaceIndex");
if (isEmpty(maxSpec))
maxSpec = numberOfSpectra - 1;

// Check 'StartSpectrum' is in range 0-numberOfSpectra
if (minSpec > numberOfSpectra - 1) {
g_log.error("StartWorkspaceIndex out of range!");
throw std::out_of_range("StartSpectrum out of range!");
}
if (maxSpec > numberOfSpectra - 1) {
g_log.error("EndWorkspaceIndex out of range!");
throw std::out_of_range("EndWorkspaceIndex out of range!");
}
if (maxSpec < minSpec) {
g_log.error("StartWorkspaceIndex must be less than or equal to "
"EndWorkspaceIndex");
throw std::out_of_range("StartWorkspaceIndex must be less than or equal "
"to EndWorkspaceIndex");
}
m_spectrumList.reserve(maxSpec - minSpec + 1);
for (specid_t i = minSpec; i <= maxSpec; ++i) {
m_spectrumList.push_back(i);
}
}
}

Expand Down Expand Up @@ -477,7 +470,7 @@ size_t ExtractSpectra::getXMax(const int wsIndex) {
* workspace
*/
void ExtractSpectra::cropRagged(API::MatrixWorkspace_sptr outputWorkspace,
int inIndex, int outIndex) {
int inIndex, int outIndex) {
MantidVec &Y = outputWorkspace->dataY(outIndex);
MantidVec &E = outputWorkspace->dataE(outIndex);
const size_t size = Y.size();
Expand Down
64 changes: 11 additions & 53 deletions Code/Mantid/Framework/Algorithms/src/RemoveMaskedSpectra.cpp
Expand Up @@ -77,61 +77,19 @@ void RemoveMaskedSpectra::exec() {
}

// Find indices of the unmasked spectra.
std::vector<size_t> indices;
std::vector<specid_t> indices;
makeIndexList(indices, maskedWorkspace.get());

// Number of spectra in the cropped workspace.
size_t nSpectra = indices.size();
// Number of bins/data points in the cropped workspace.
size_t nBins = inputWorkspace->blocksize();
size_t histogram = inputWorkspace->isHistogramData() ? 1 : 0;

// Create the output workspace.
MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create(
inputWorkspace, nSpectra, nBins + histogram, nBins);

// If this is a Workspace2D, get the spectra axes for copying in the spectraNo
// later.
Axis *inAxis1(NULL), *outAxis1(NULL);
TextAxis *outTxtAxis(NULL);
if (inputWorkspace->axes() > 1) {
inAxis1 = inputWorkspace->getAxis(1);
outAxis1 = outputWorkspace->getAxis(1);
outTxtAxis = dynamic_cast<TextAxis *>(outAxis1);
}
auto extract = createChildAlgorithm("ExtractSpectra");
extract->initialize();
extract->setRethrows(true);

// Check for common boundaries in input workspace.
bool commonBoundaries = WorkspaceHelpers::commonBoundaries(inputWorkspace);
extract->setProperty("InputWorkspace", inputWorkspace);
extract->setProperty("SpectrumList", indices);

MantidVecPtr newX;
if (commonBoundaries) newX = inputWorkspace->refX(0);
extract->execute();

// Add spectra to the output workspace.
for (size_t j = 0; j < nSpectra; ++j) {
auto i = indices[j];
// copy the x bins
if (commonBoundaries) {
outputWorkspace->setX(j, newX);
} else {
outputWorkspace->setX(j, inputWorkspace->refX(i));
}
// Copy the y values and errors.
outputWorkspace->getSpectrum(j)->setData(inputWorkspace->readY(i), inputWorkspace->readE(i));
// Copy spectrum number & detectors
outputWorkspace->getSpectrum(j)
->copyInfoFrom(*inputWorkspace->getSpectrum(i));
// Copy over the axis entry for each spectrum, regardless of the type of
// axes present.
if (inAxis1) {
if (outAxis1->isText()) {
outTxtAxis->setLabel(j, inAxis1->label(i));
} else if (!outAxis1->isSpectra()) // handled by copyInfoFrom line
{
dynamic_cast<NumericAxis *>(outAxis1)
->setValue(j, inAxis1->operator()(i));
}
}
}
MatrixWorkspace_sptr outputWorkspace = extract->getProperty("OutputWorkspace");
setProperty("OutputWorkspace", outputWorkspace);
}

Expand All @@ -140,12 +98,12 @@ void RemoveMaskedSpectra::exec() {
/// @param indices :: A reference to a vector to fill with the indices.
/// @param maskedWorkspace :: A workspace with masking information.
void RemoveMaskedSpectra::makeIndexList(
std::vector<size_t> &indices, const API::MatrixWorkspace *maskedWorkspace) {
std::vector<specid_t> &indices, const API::MatrixWorkspace *maskedWorkspace) {
auto mask = dynamic_cast<const DataObjects::MaskWorkspace *>(maskedWorkspace);
if (mask) {
for (size_t i = 0; i < mask->getNumberHistograms(); ++i) {
if (mask->readY(i)[0] == 0.0) {
indices.push_back(i);
indices.push_back(static_cast<specid_t>(i));
}
}
} else {
Expand All @@ -157,7 +115,7 @@ void RemoveMaskedSpectra::makeIndexList(
continue;
}
if (det->isMasked()) {
indices.push_back(i);
indices.push_back(static_cast<specid_t>(i));
}
}
}
Expand Down

0 comments on commit 2e0326d

Please sign in to comment.