diff --git a/Framework/DataHandling/inc/MantidDataHandling/GroupDetectors2.h b/Framework/DataHandling/inc/MantidDataHandling/GroupDetectors2.h
index 40a3464743137b473018fbe01fda0fbb3d9114ec..9f5dcbb1e915bd6c2daf79434391ebc78f8e5854 100644
--- a/Framework/DataHandling/inc/MantidDataHandling/GroupDetectors2.h
+++ b/Framework/DataHandling/inc/MantidDataHandling/GroupDetectors2.h
@@ -204,10 +204,6 @@ private:
DataObjects::EventWorkspace_sptr outputWS,
const double prog4Copy);
- /// Returns true if detectors exists and is masked
- bool isMaskedDetector(const API::SpectrumInfo &detector,
- const size_t index) const;
-
/// Copy the ungrouped spectra from the input workspace to the output
template <class TIn, class TOut>
void moveOthers(const std::set<int64_t> &unGroupedSet, const TIn &inputWS,
diff --git a/Framework/DataHandling/src/GroupDetectors2.cpp b/Framework/DataHandling/src/GroupDetectors2.cpp
index d1ad725e416fe4074ce99053b0a5f4a6592d29c6..915af2ebb564e56c9184b87a0e7239cde74c1a2e 100644
--- a/Framework/DataHandling/src/GroupDetectors2.cpp
+++ b/Framework/DataHandling/src/GroupDetectors2.cpp
@@ -40,6 +40,7 @@ using namespace DataObjects;
using std::size_t;
namespace { // anonymous namespace
+enum class Behaviour { SUM, AVERAGE };
/**
* Translate the PerformIndexOperations processing instructions from a vector
@@ -958,70 +959,83 @@ size_t GroupDetectors2::formGroups(API::MatrixWorkspace_const_sptr inputWS,
const double prog4Copy, const bool keepAll,
const std::set<int64_t> &unGroupedSet,
Indexing::IndexInfo &indexInfo) {
- // get "Behaviour" string
- const std::string behaviour = getProperty("Behaviour");
- int bhv = 0;
- if (behaviour == "Average")
- bhv = 1;
-
- API::MatrixWorkspace_sptr beh = API::WorkspaceFactory::Instance().create(
- "Workspace2D", static_cast<int>(m_GroupWsInds.size()), 1, 1);
-
- g_log.debug() << name() << ": Preparing to group spectra into "
- << m_GroupWsInds.size() << " groups\n";
-
- // where we are copying spectra to, we start copying to the start of the
- // output workspace
+ const std::string behaviourChoice = getProperty("Behaviour");
+ const auto behaviour =
+ behaviourChoice == "Sum" ? Behaviour::SUM : Behaviour::AVERAGE;
size_t outIndex = 0;
- // Only used for averaging behaviour. We may have a 1:1 map where a Divide
- // would be waste as it would be just dividing by 1
- bool requireDivide(false);
const auto &spectrumInfo = inputWS->spectrumInfo();
-
+ const auto nFinalHistograms =
+ m_GroupWsInds.size() + (keepAll ? unGroupedSet.size() : 0);
auto spectrumGroups = std::vector<std::vector<size_t>>();
+ spectrumGroups.reserve(nFinalHistograms);
auto spectrumNumbers = std::vector<Indexing::SpectrumNumber>();
+ spectrumNumbers.reserve(nFinalHistograms);
- for (storage_map::const_iterator it = m_GroupWsInds.begin();
- it != m_GroupWsInds.end(); ++it) {
+ for (const auto &group : m_GroupWsInds) {
// This is the grouped spectrum
auto &outSpec = outputWS->getSpectrum(outIndex);
// The spectrum number of the group is the key
- spectrumNumbers.push_back(it->first);
+ spectrumNumbers.emplace_back(group.first);
// Start fresh with no detector IDs
outSpec.clearDetectorIDs();
// Copy over X data from first spectrum, the bin boundaries for all spectra
// are assumed to be the same here
outSpec.setSharedX(inputWS->sharedX(0));
- auto outputHistogram = outSpec.histogram();
// Keep track of number of detectors required for masking
- size_t nonMaskedSpectra(0);
- auto spectrumGroup = std::vector<size_t>();
-
- for (auto originalWI : it->second) {
- // detectors to add to firstSpecNum
- const auto &inputSpectrum = inputWS->getSpectrum(originalWI);
-
- outputHistogram += inputSpectrum.histogram();
- outSpec.addDetectorIDs(inputSpectrum.getDetectorIDs());
-
- if (!isMaskedDetector(spectrumInfo, originalWI))
- ++nonMaskedSpectra;
-
- spectrumGroup.push_back(originalWI);
+ std::vector<size_t> spectrumGroup;
+ spectrumGroup.reserve(group.second.size());
+
+ auto &Ys = outSpec.mutableY();
+ auto &Es = outSpec.mutableE();
+ std::vector<double> sum(Ys.size(), 0.);
+ std::vector<double> errorSum(Ys.size(), 0.);
+ std::vector<int> count(Ys.size(), 0);
+ for (auto originalWI : group.second) {
+ const auto &inSpec = inputWS->getSpectrum(originalWI);
+ outSpec.addDetectorIDs(inSpec.getDetectorIDs());
+ spectrumGroup.emplace_back(originalWI);
+ if (spectrumInfo.hasDetectors(originalWI) &&
+ spectrumInfo.isMasked(originalWI)) {
+ continue;
+ }
+ const auto &Ys = inputWS->y(originalWI);
+ const auto &Es = inputWS->e(originalWI);
+ if (inputWS->hasMaskedBins(originalWI)) {
+ const auto &maskedBins = inputWS->maskedBins(originalWI);
+ for (size_t binIndex = 0; binIndex < Ys.size(); ++binIndex) {
+ if (maskedBins.count(binIndex) == 0) {
+ sum[binIndex] += Ys[binIndex];
+ errorSum[binIndex] += Es[binIndex] * Es[binIndex];
+ count[binIndex] += 1;
+ }
+ }
+ } else {
+ for (size_t binIndex = 0; binIndex < Ys.size(); ++binIndex) {
+ sum[binIndex] += Ys[binIndex];
+ errorSum[binIndex] += Es[binIndex] * Es[binIndex];
+ count[binIndex] += 1;
+ }
+ }
+ }
+ spectrumGroups.emplace_back(std::move(spectrumGroup));
+ for (size_t binIndex = 0; binIndex < sum.size(); ++binIndex) {
+ errorSum[binIndex] = std::sqrt(errorSum[binIndex]);
+ if (behaviour == Behaviour::AVERAGE) {
+ const auto n = static_cast<double>(count[binIndex]);
+ if (n != 0) {
+ sum[binIndex] /= n;
+ errorSum[binIndex] /= n;
+ } else {
+ sum[binIndex] = 0;
+ errorSum[binIndex] = 0;
+ }
+ }
+ Ys[binIndex] = sum[binIndex];
+ Es[binIndex] = errorSum[binIndex];
}
-
- spectrumGroups.push_back(spectrumGroup);
-
- outSpec.setHistogram(outputHistogram);
-
- if (nonMaskedSpectra == 0)
- ++nonMaskedSpectra; // Avoid possible divide by zero
- if (!requireDivide)
- requireDivide = (nonMaskedSpectra > 1);
- beh->mutableY(outIndex)[0] = static_cast<double>(nonMaskedSpectra);
// make regular progress reports and check for cancelling the algorithm
if (outIndex % INTERVAL == 0) {
@@ -1042,28 +1056,15 @@ size_t GroupDetectors2::formGroups(API::MatrixWorkspace_const_sptr inputWS,
if (originalWI < 0)
continue;
- spectrumGroups.push_back(std::vector<size_t>(1, originalWI));
+ spectrumGroups.emplace_back(std::vector<size_t>(1, originalWI));
auto spectrumNumber = inputWS->getSpectrum(originalWI).getSpectrumNo();
- spectrumNumbers.push_back(spectrumNumber);
+ spectrumNumbers.emplace_back(spectrumNumber);
}
}
indexInfo = Indexing::group(inputWS->indexInfo(), std::move(spectrumNumbers),
spectrumGroups);
-
- if (bhv == 1 && requireDivide) {
- g_log.debug() << "Running Divide algorithm to perform averaging.\n";
- Mantid::API::IAlgorithm_sptr divide = createChildAlgorithm("Divide");
- divide->initialize();
- divide->setProperty<API::MatrixWorkspace_sptr>("LHSWorkspace", outputWS);
- divide->setProperty<API::MatrixWorkspace_sptr>("RHSWorkspace", beh);
- divide->setProperty<API::MatrixWorkspace_sptr>("OutputWorkspace", outputWS);
- divide->execute();
- }
-
- g_log.debug() << name() << " created " << outIndex
- << " new grouped spectra\n";
return outIndex;
}
@@ -1126,7 +1127,8 @@ GroupDetectors2::formGroupsEvent(DataObjects::EventWorkspace_const_sptr inputWS,
// detectors to add to the output spectrum
outEL.addDetectorIDs(fromEL.getDetectorIDs());
- if (!isMaskedDetector(spectrumInfo, originalWI)) {
+ if (!spectrumInfo.hasDetectors(originalWI) ||
+ !spectrumInfo.isMasked(originalWI)) {
++nonMaskedSpectra;
}
}
@@ -1162,16 +1164,6 @@ GroupDetectors2::formGroupsEvent(DataObjects::EventWorkspace_const_sptr inputWS,
return outIndex;
}
-bool GroupDetectors2::isMaskedDetector(const API::SpectrumInfo &spectrum,
- const size_t index) const {
- if (spectrum.hasDetectors(index)) {
- return spectrum.isMasked(index);
- } else {
- // Can't be masked if it doesn't exist
- return false;
- }
-}
-
// RangeHelper
/** Expands any ranges in the input string of non-negative integers, eg. "1 3-5
* 4" -> "1 3 4 5 4"
diff --git a/Framework/DataHandling/test/GroupDetectors2Test.h b/Framework/DataHandling/test/GroupDetectors2Test.h
index fb8c68270cd44a9d6b1b089772ecb8e8242fa54a..2a12cec0ec85c3f16499b2e756d6db75791ec163 100644
--- a/Framework/DataHandling/test/GroupDetectors2Test.h
+++ b/Framework/DataHandling/test/GroupDetectors2Test.h
@@ -608,6 +608,67 @@ public:
TS_ASSERT_EQUALS(output->y(0)[1], 1.5);
}
+ void testAverageBehaviourWithMaskedBins() {
+ createTestWorkspace(inputWSName, 0);
+ MatrixWorkspace_sptr input =
+ AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(
+ inputWSName);
+ input->flagMasked(0, 0);
+ GroupDetectors2 gd2;
+ gd2.initialize();
+ gd2.setChild(true);
+ gd2.setRethrows(true);
+ gd2.setPropertyValue("InputWorkspace", inputWSName);
+ gd2.setPropertyValue("OutputWorkspace", "_unused_for_child");
+ gd2.setPropertyValue("WorkspaceIndexList", "0,1");
+ gd2.setPropertyValue("Behaviour", "Average");
+ TS_ASSERT_THROWS_NOTHING(gd2.execute());
+ TS_ASSERT(gd2.isExecuted())
+ MatrixWorkspace_sptr output = gd2.getProperty("OutputWorkspace");
+ TS_ASSERT_EQUALS(output->getNumberHistograms(), 1)
+ const auto &spectrum = output->getSpectrum(0);
+ const auto &detIds = spectrum.getDetectorIDs();
+ TS_ASSERT_EQUALS(detIds.size(), 2)
+ TS_ASSERT_DIFFERS(detIds.find(0), detIds.end())
+ TS_ASSERT_DIFFERS(detIds.find(1), detIds.end())
+ const auto &y = output->y(0);
+ const auto &e = output->e(0);
+ for (size_t i = 0; i < y.size(); ++i) {
+ const double expectedSignal = i == 0 ? 2. : (1. + 2.) / 2.;
+ TS_ASSERT_EQUALS(y[i], expectedSignal)
+ const double expectedError = i == 0 ? 1. : std::sqrt(2.) / 2.;
+ TS_ASSERT_EQUALS(e[i], expectedError)
+ }
+ }
+
+ void testSumBehaviourWithMaskedBins() {
+ createTestWorkspace(inputWSName, 0);
+ MatrixWorkspace_sptr input =
+ AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(
+ inputWSName);
+ input->flagMasked(0, 0);
+ GroupDetectors2 gd2;
+ gd2.initialize();
+ gd2.setChild(true);
+ gd2.setRethrows(true);
+ gd2.setPropertyValue("InputWorkspace", inputWSName);
+ gd2.setPropertyValue("OutputWorkspace", "_unused_for_child");
+ gd2.setPropertyValue("WorkspaceIndexList", "0,1");
+ gd2.setPropertyValue("Behaviour", "Sum");
+ TS_ASSERT_THROWS_NOTHING(gd2.execute());
+ TS_ASSERT(gd2.isExecuted())
+ MatrixWorkspace_sptr output = gd2.getProperty("OutputWorkspace");
+ TS_ASSERT_EQUALS(output->getNumberHistograms(), 1)
+ const auto &y = output->y(0);
+ const auto &e = output->e(0);
+ for (size_t i = 0; i < y.size(); ++i) {
+ const double expectedSignal = i == 0 ? 2. : 1. + 2.;
+ TS_ASSERT_EQUALS(y[i], expectedSignal)
+ const double expectedError = i == 0 ? 1. : std::sqrt(2.);
+ TS_ASSERT_EQUALS(e[i], expectedError)
+ }
+ }
+
void testEvents() {
int numPixels = 5;
int numBins = 5;
@@ -1018,11 +1079,36 @@ public:
output->y(0)[0], 0.00001);
}
+ void test_masked_detids_get_propagated() {
+ createTestWorkspace(inputWSName, 0);
+ MatrixWorkspace_sptr input =
+ AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(
+ inputWSName);
+ input->mutableDetectorInfo().setMasked(0, true);
+ GroupDetectors2 gd2;
+ gd2.initialize();
+ gd2.setChild(true);
+ gd2.setRethrows(true);
+ gd2.setPropertyValue("InputWorkspace", inputWSName);
+ gd2.setPropertyValue("OutputWorkspace", "_unused_for_child");
+ gd2.setPropertyValue("WorkspaceIndexList", "0,1");
+ gd2.setPropertyValue("Behaviour", "Sum");
+ TS_ASSERT_THROWS_NOTHING(gd2.execute());
+ TS_ASSERT(gd2.isExecuted())
+ MatrixWorkspace_sptr output = gd2.getProperty("OutputWorkspace");
+ TS_ASSERT_EQUALS(output->getNumberHistograms(), 1)
+ const auto &spectrum = output->getSpectrum(0);
+ const auto &ids = spectrum.getDetectorIDs();
+ TS_ASSERT(ids.size() == 2)
+ TS_ASSERT_DIFFERS(ids.find(0), ids.end())
+ TS_ASSERT_DIFFERS(ids.find(1), ids.end())
+ }
+
private:
const std::string inputWSName, offsetWSName, outputWSNameBase, inputFile;
- enum constants { NHIST = 6, NBINS = 4 };
+ enum { NHIST = 6, NBINS = 4 };
- void createTestWorkspace(const std::string &name, const int offset) {
+ static void createTestWorkspace(const std::string &name, const int offset) {
// Set up a small workspace for testing
auto space2D = createWorkspace<Workspace2D>(NHIST, NBINS + 1, NBINS);
space2D->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
@@ -1101,48 +1187,63 @@ private:
class GroupDetectors2TestPerformance : public CxxTest::TestSuite {
public:
- void setUp() override {
- constexpr int numGroups = 2;
+ static GroupDetectors2TestPerformance *createSuite() {
+ return new GroupDetectors2TestPerformance();
+ }
+ static void destroySuite(GroupDetectors2TestPerformance *suite) {
+ delete suite;
+ }
+
+ GroupDetectors2TestPerformance()
+ : inputEventWs(nullptr), inputMatrixWs(nullptr), groupWs(nullptr), alg() {
+ constexpr int numGroups = 40;
// This controls speed of test
constexpr int bankPixelWidth = 30;
+ constexpr int numBins = 1000;
- inputWs = WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(
- numGroups, bankPixelWidth);
- AnalysisDataService::Instance().addOrReplace(nxsWSname, inputWs);
-
+ inputEventWs =
+ WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(
+ numGroups, bankPixelWidth);
+ inputMatrixWs =
+ WorkspaceCreationHelper::create2DWorkspaceWithRectangularInstrument(
+ numGroups, bankPixelWidth, numBins);
// Create an axis for each pixel.
- for (size_t pix = 0; pix < inputWs->getNumberHistograms(); pix++) {
- size_t xAxisSize = inputWs->x(pix).size();
+ for (size_t pix = 0; pix < inputEventWs->getNumberHistograms(); pix++) {
+ size_t xAxisSize = inputEventWs->x(pix).size();
Mantid::HistogramData::HistogramX axisVals(xAxisSize, 1.0);
- inputWs->mutableX(pix) = axisVals;
- inputWs->getSpectrum(pix).addEventQuickly(TofEvent(1000.0));
+ inputEventWs->mutableX(pix) = axisVals;
+ inputEventWs->getSpectrum(pix).addEventQuickly(TofEvent(1000.0));
}
-
setupGroupWS(numGroups);
alg.initialize();
- alg.setPropertyValue("InputWorkspace", nxsWSname);
- alg.setPropertyValue("OutputWorkspace", outputws);
- alg.setPropertyValue("CopyGroupingFromWorkspace", groupWSName);
-
+ alg.setProperty("OutputWorkspace", "_unused_for_child");
+ alg.setProperty("CopyGroupingFromWorkspace", groupWs);
+ alg.setChild(true);
alg.setRethrows(true);
}
- void testGroupDetectors2Performance() {
- TS_ASSERT_THROWS_NOTHING(alg.execute());
+ void testGroupDetectors2EventPerformance() {
+ alg.setProperty("InputWorkspace", inputEventWs);
+ for (size_t i = 0; i < 100; ++i) {
+ TS_ASSERT_THROWS_NOTHING(alg.execute());
+ }
}
- void tearDown() override {
- AnalysisDataService::Instance().remove(groupWSName);
- AnalysisDataService::Instance().remove(nxsWSname);
- AnalysisDataService::Instance().remove(outputws);
+ void testGroupDetectors2HistogramPerformance() {
+ alg.setProperty("InputWorkspace", inputMatrixWs);
+ for (size_t i = 0; i < 50; ++i) {
+ TS_ASSERT_THROWS_NOTHING(alg.execute());
+ }
}
+ void tearDown() override {}
+
void setupGroupWS(const size_t numGroups) {
// ------------ Create a grouping workspace to match -------------
- groupWs = boost::make_shared<GroupingWorkspace>(inputWs->getInstrument());
- AnalysisDataService::Instance().addOrReplace(groupWSName, groupWs);
+ groupWs =
+ boost::make_shared<GroupingWorkspace>(inputEventWs->getInstrument());
// fill in some groups
constexpr size_t startingGroupNo = 1;
@@ -1157,11 +1258,8 @@ public:
}
private:
- const std::string nxsWSname = "GroupDetectors2TestTarget_ws";
- const std::string groupWSName = nxsWSname + "_GROUP";
- const std::string outputws = nxsWSname + "_grouped";
-
- EventWorkspace_sptr inputWs;
+ EventWorkspace_sptr inputEventWs;
+ MatrixWorkspace_sptr inputMatrixWs;
GroupingWorkspace_sptr groupWs;
GroupDetectors2 alg;
diff --git a/docs/source/release/v3.14.0/framework.rst b/docs/source/release/v3.14.0/framework.rst
index b6bd62c5ee12de14ca469a4eef54599585e5cb95..7e8f0656a5ace58a655ecb023d38f7c9e303afb0 100644
--- a/docs/source/release/v3.14.0/framework.rst
+++ b/docs/source/release/v3.14.0/framework.rst
@@ -55,6 +55,7 @@ Improvements
- :ref:`LoadSampleShape <algm-LoadSampleShape-v1>` now supports loading from binary .stl files.
- :ref:`MaskDetectorsIf <algm-MaskDetectorsIf>` now supports masking a workspace in addition to writing the masking information to a calfile.
- :ref:`LoadSampleShape <algm-LoadSampleShape-v1>` now supports loading from binary .stl files.
+- :ref:`GroupDetectors <algm-GroupDetectors>` now takes masked bins correctly into account when processing histogram workspaces.
Bugfixes
########