Mantid
Loading...
Searching...
No Matches
IntegrateByComponent.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 +
13#include "MantidTypes/SpectrumDefinition.h"
14
15#include <gsl/gsl_statistics.h>
16#include <unordered_map>
17
18namespace Mantid::Algorithms {
19
20// Register the algorithm into the AlgorithmFactory
21DECLARE_ALGORITHM(IntegrateByComponent)
22
23using namespace Mantid::API;
24using namespace Mantid::Kernel;
25
26//----------------------------------------------------------------------------------------------
28const std::string IntegrateByComponent::name() const { return "IntegrateByComponent"; }
29
31int IntegrateByComponent::version() const { return 1; }
32
34const std::string IntegrateByComponent::category() const { return "Utility\\Workspaces"; }
35
36//----------------------------------------------------------------------------------------------
37
38//----------------------------------------------------------------------------------------------
42 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input,
43 std::make_shared<HistogramValidator>()),
44 "The input workspace.");
45 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
46 "The output workspace.");
47 auto mustBePosInt = std::make_shared<BoundedValidator<int>>();
48 mustBePosInt->setLower(0);
49 declareProperty("LevelsUp", 0, mustBePosInt,
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");
53}
54
55//----------------------------------------------------------------------------------------------
59 MatrixWorkspace_sptr inputWS = this->getProperty("InputWorkspace");
60 int parents = getProperty("LevelsUp");
61 // Make sure it's integrated
62 auto childAlg = createChildAlgorithm("Integration", 0, 0.2);
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();
70 MatrixWorkspace_sptr integratedWS = childAlg->getProperty("OutputWorkspace");
71
72 if (parents > 0) {
73 std::vector<std::vector<size_t>> specmap = makeMap(integratedWS, parents);
74 API::Progress prog(this, 0.3, 1.0, specmap.size());
75 // calculate averages
76 const auto &spectrumInfo = integratedWS->spectrumInfo();
77 for (auto hists : specmap) {
78 prog.report();
79 std::vector<double> averageYInput, averageEInput;
80
82 for (int i = 0; i < static_cast<int>(hists.size()); ++i) { // NOLINT
84
85 if (spectrumInfo.isMonitor(hists[i]))
86 continue;
87 if (spectrumInfo.isMasked(hists[i]))
88 continue;
89
90 const double yValue = integratedWS->y(hists[i])[0];
91 const double eValue = integratedWS->e(hists[i])[0];
92
93 if (!std::isfinite(yValue) || !std::isfinite(eValue)) // NaNs/Infs
94 continue;
95
96 // Now we have a good value
97 PARALLEL_CRITICAL(IntegrateByComponent_good) {
98 averageYInput.emplace_back(yValue);
99 averageEInput.emplace_back(eValue * eValue);
100 }
101
103 }
105
106 double averageY, averageE;
107 if (averageYInput.empty()) {
108 g_log.information("some group has no valid histograms. Will use 0 for average.");
109 averageY = 0.;
110 averageE = 0.;
111 } else {
112 averageY = gsl_stats_mean(&averageYInput[0], 1, averageEInput.size());
113 averageE = std::sqrt(gsl_stats_mean(&averageEInput[0], 1, averageYInput.size()));
114 }
115
116 PARALLEL_FOR_IF(Kernel::threadSafe(*integratedWS))
117 for (int i = 0; i < static_cast<int>(hists.size()); ++i) { // NOLINT
119 if (spectrumInfo.isMonitor(hists[i]))
120 continue;
121 if (spectrumInfo.isMasked(hists[i]))
122 continue;
123
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)) // NaNs/Infs
127 continue;
128
129 // Now we have a good value
130 PARALLEL_CRITICAL(IntegrateByComponent_setaverage) {
131 integratedWS->dataY(hists[i])[0] = averageY;
132 integratedWS->dataE(hists[i])[0] = averageE;
133 }
134
136 }
138 }
139 }
140 // Assign it to the output workspace property
141 setProperty("OutputWorkspace", integratedWS);
142}
143
151std::vector<std::vector<size_t>> IntegrateByComponent::makeInstrumentMap(const API::MatrixWorkspace_sptr &countsWS) {
152 std::vector<std::vector<size_t>> mymap;
153 std::vector<size_t> single;
154
155 for (size_t i = 0; i < countsWS->getNumberHistograms(); i++) {
156 single.emplace_back(i);
157 }
158 mymap.emplace_back(single);
159 return mymap;
160}
161
169std::vector<std::vector<size_t>> IntegrateByComponent::makeMap(const API::MatrixWorkspace_sptr &countsWS, int parents) {
170 std::unordered_multimap<Mantid::Geometry::ComponentID, size_t> mymap;
171
172 if (parents == 0) // this should not happen in this file, but if one reuses
173 // the function and parents==0, the program has a sudden end
174 // without this check.
175 {
176 return makeInstrumentMap(countsWS);
177 }
178
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");
184 continue;
185 }
186
187 const auto detIdx = spectrumInfo.spectrumDefinition(i)[0].first;
188 std::vector<std::shared_ptr<const Mantid::Geometry::IComponent>> anc = detectorInfo.detector(detIdx).getAncestors();
189
190 if (anc.size() < static_cast<size_t>(parents)) {
191 g_log.warning("Too many levels up. Will ignore LevelsUp");
192 return makeInstrumentMap(countsWS);
193 }
194 mymap.emplace(anc[parents - 1]->getComponentID(), i);
195 }
196
197 std::vector<std::vector<size_t>> speclist;
198 std::vector<size_t> speclistsingle;
199
200 auto s_it = mymap.begin();
201 for (auto m_it = mymap.begin(); m_it != mymap.end(); m_it = s_it) {
202 Mantid::Geometry::ComponentID theKey = (*m_it).first;
203 auto keyRange = mymap.equal_range(theKey);
204
205 // Iterate over all map elements with key == theKey
206 speclistsingle.clear();
207 for (s_it = keyRange.first; s_it != keyRange.second; ++s_it) {
208 speclistsingle.emplace_back((*s_it).second);
209 }
210 speclist.emplace_back(speclistsingle);
211 }
212
213 return speclist;
214}
215
216} // 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_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.
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
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
Definition: IComponent.h:51
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.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
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.
Definition: MultiThreaded.h:22
constexpr int EMPTY_INT() noexcept
Returns what we consider an "empty" integer within a property.
Definition: EmptyValues.h:25
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54