13#include "MantidTypes/SpectrumDefinition.h"
15#include <gsl/gsl_statistics.h>
16#include <unordered_map>
43 std::make_shared<HistogramValidator>()),
44 "The input workspace.");
46 "The output workspace.");
47 auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
48 mustBePosInt->setLower(0);
50 "Levels above pixel that will be used to compute the average.\n"
51 "If no level is specified, the median is over the whole instrument.\n If "
52 "0, it will just return the integrated values in each pixel");
63 childAlg->setProperty(
"InputWorkspace", inputWS);
64 childAlg->setProperty(
"StartWorkspaceIndex", 0);
65 childAlg->setProperty(
"EndWorkspaceIndex",
EMPTY_INT());
66 childAlg->setProperty(
"RangeLower", inputWS->getXMin());
67 childAlg->setProperty(
"RangeUpper", inputWS->getXMax());
68 childAlg->setPropertyValue(
"IncludePartialBins",
"1");
69 childAlg->executeAsChildAlg();
73 std::vector<std::vector<size_t>> specmap =
makeMap(integratedWS, parents);
76 const auto &spectrumInfo = integratedWS->spectrumInfo();
77 for (
auto hists : specmap) {
79 std::vector<double> averageYInput, averageEInput;
82 for (
int i = 0; i < static_cast<int>(hists.size()); ++i) {
85 if (spectrumInfo.isMonitor(hists[i]))
87 if (spectrumInfo.isMasked(hists[i]))
90 const double yValue = integratedWS->y(hists[i])[0];
91 const double eValue = integratedWS->e(hists[i])[0];
93 if (!std::isfinite(yValue) || !std::isfinite(eValue))
98 averageYInput.emplace_back(yValue);
99 averageEInput.emplace_back(eValue * eValue);
106 double averageY, averageE;
107 if (averageYInput.empty()) {
108 g_log.
information(
"some group has no valid histograms. Will use 0 for average.");
112 averageY = gsl_stats_mean(&averageYInput[0], 1, averageEInput.size());
113 averageE = std::sqrt(gsl_stats_mean(&averageEInput[0], 1, averageYInput.size()));
117 for (
int i = 0; i < static_cast<int>(hists.size()); ++i) {
119 if (spectrumInfo.isMonitor(hists[i]))
121 if (spectrumInfo.isMasked(hists[i]))
124 const double yValue = integratedWS->y(hists[i])[0];
125 const double eValue = integratedWS->e(hists[i])[0];
126 if (!std::isfinite(yValue) || !std::isfinite(eValue))
131 integratedWS->dataY(hists[i])[0] = averageY;
132 integratedWS->dataE(hists[i])[0] = averageE;
152 std::vector<std::vector<size_t>> mymap;
153 std::vector<size_t> single;
155 for (
size_t i = 0; i < countsWS->getNumberHistograms(); i++) {
156 single.emplace_back(i);
158 mymap.emplace_back(single);
170 std::unordered_multimap<Mantid::Geometry::ComponentID, size_t> mymap;
179 const auto spectrumInfo = countsWS->spectrumInfo();
180 const auto &detectorInfo = countsWS->detectorInfo();
181 for (
size_t i = 0; i < countsWS->getNumberHistograms(); i++) {
182 if (!spectrumInfo.hasDetectors(i)) {
183 g_log.
debug(
"Spectrum has no detector, skipping");
187 const auto detIdx = spectrumInfo.spectrumDefinition(i)[0].first;
188 std::vector<std::shared_ptr<const Mantid::Geometry::IComponent>> anc = detectorInfo.detector(detIdx).getAncestors();
190 if (anc.size() <
static_cast<size_t>(parents)) {
191 g_log.
warning(
"Too many levels up. Will ignore LevelsUp");
194 mymap.emplace(anc[parents - 1]->getComponentID(), i);
197 std::vector<std::vector<size_t>> speclist;
198 std::vector<size_t> speclistsingle;
200 auto s_it = mymap.begin();
201 for (
auto m_it = mymap.begin(); m_it != mymap.end(); m_it = s_it) {
203 auto keyRange = mymap.equal_range(theKey);
206 speclistsingle.clear();
207 for (s_it = keyRange.first; s_it != keyRange.second; ++s_it) {
208 speclistsingle.emplace_back((*s_it).second);
210 speclist.emplace_back(speclistsingle);
#define DECLARE_ALGORITHM(classname)
#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...
#define PARALLEL_CRITICAL(name)
#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.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
Helper class for reporting progress from algorithms.
A property class for workspaces.
void init() override
Initialize the algorithm's properties.
std::vector< std::vector< size_t > > makeMap(const API::MatrixWorkspace_sptr &countsWS, int parents)
method to check which spectra should be averaged
std::vector< std::vector< size_t > > makeInstrumentMap(const API::MatrixWorkspace_sptr &countsWS)
method to create the map with all spectra
int version() const override
Algorithm's version for identification.
const std::string name() const override
Algorithm's name for identification.
const std::string category() const override
Algorithm's category for identification.
void exec() override
Execute the algorithm.
base class for Geometric IComponent
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
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.
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
@ Input
An input workspace.
@ Output
An output workspace.