-
Notifications
You must be signed in to change notification settings - Fork 122
/
MDHistoToWorkspace2D.cpp
149 lines (127 loc) · 5.48 KB
/
MDHistoToWorkspace2D.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
/**
* This algorithm flattens a MDHistoWorkspace to a Workspace2D. Mantid has far
*more tools
* to deal with W2D then for MD ones.
*
* copyright: do not bother me, see Mantid copyright
*
* Mark Koennecke, November 2012
*
* Added copying of meta data. Mark Koennecke, July 2013
*/
#include "MantidSINQ/MDHistoToWorkspace2D.h"
#include "MantidAPI/IMDHistoWorkspace.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include "MantidKernel/Logger.h"
#include <cmath>
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MDHistoToWorkspace2D)
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;
using namespace Mantid;
using Mantid::HistogramData::Counts;
using Mantid::HistogramData::Points;
namespace {
Logger logger("MDHistoToWorkspace2D");
}
// A reference to the logger is provided by the base class, it is called g_log.
// It is used to print out information, warning and error messages
MDHistoToWorkspace2D::MDHistoToWorkspace2D() : Mantid::API::Algorithm(), m_rank(0), m_currentSpectra(0) {}
void MDHistoToWorkspace2D::init() {
declareProperty(std::make_unique<WorkspaceProperty<IMDHistoWorkspace>>("InputWorkspace", "", Direction::Input));
declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output));
}
void MDHistoToWorkspace2D::exec() {
IMDHistoWorkspace_sptr inWS = IMDHistoWorkspace_sptr(getProperty("InputWorkspace"));
m_rank = inWS->getNumDims();
size_t nSpectra = calculateNSpectra(inWS);
logger.debug() << "nSpectra = " << nSpectra << '\n';
std::shared_ptr<const IMDDimension> lastDim = inWS->getDimension(m_rank - 1);
logger.debug() << "spectraLength = " << lastDim->getNBins() << '\n';
Mantid::DataObjects::Workspace2D_sptr outWS;
outWS = std::dynamic_pointer_cast<Mantid::DataObjects::Workspace2D>(
WorkspaceFactory::Instance().create("Workspace2D", nSpectra, lastDim->getNBins(), lastDim->getNBins()));
for (size_t i = 0; i < nSpectra; ++i)
outWS->getSpectrum(i).setDetectorID(static_cast<detid_t>(i + 1));
outWS->setYUnit("Counts");
struct FreeDeleter {
void operator()(void *x) { free(x); }
};
auto pos = std::unique_ptr<coord_t, FreeDeleter>(reinterpret_cast<coord_t *>(malloc(m_rank * sizeof(coord_t))));
memset(pos.get(), 0, m_rank * sizeof(coord_t));
m_currentSpectra = 0;
recurseData(inWS, outWS, 0, pos.get());
copyMetaData(inWS, outWS);
setProperty("OutputWorkspace", std::dynamic_pointer_cast<Workspace>(outWS));
}
size_t MDHistoToWorkspace2D::calculateNSpectra(const IMDHistoWorkspace_sptr &inWS) {
size_t nSpectra = 1;
for (size_t i = 0; i < m_rank - 1; i++) {
std::shared_ptr<const IMDDimension> dim = inWS->getDimension(i);
nSpectra *= dim->getNBins();
}
return nSpectra;
}
void MDHistoToWorkspace2D::recurseData(const IMDHistoWorkspace_sptr &inWS, const Workspace2D_sptr &outWS,
size_t currentDim, coord_t *pos) {
std::shared_ptr<const IMDDimension> dim = inWS->getDimension(currentDim);
if (currentDim == m_rank - 1) {
Counts counts(dim->getNBins());
auto &Y = counts.mutableData();
for (unsigned int j = 0; j < dim->getNBins(); j++) {
pos[currentDim] = dim->getX(j);
Y[j] = inWS->getSignalAtCoord(pos, static_cast<Mantid::API::MDNormalization>(0));
}
Points points(dim->getNBins());
auto &xData = points.mutableData();
for (unsigned int i = 0; i < dim->getNBins(); i++) {
xData[i] = dim->getX(i);
}
outWS->setHistogram(m_currentSpectra, std::move(points), std::move(counts));
outWS->getSpectrum(m_currentSpectra).setSpectrumNo(static_cast<specnum_t>(m_currentSpectra));
m_currentSpectra++;
} else {
// recurse deeper
for (int i = 0; i < static_cast<int>(dim->getNBins()); i++) {
pos[currentDim] = dim->getX(i);
recurseData(inWS, outWS, currentDim + 1, pos);
}
}
}
void MDHistoToWorkspace2D::checkW2D(const Mantid::DataObjects::Workspace2D_sptr &outWS) {
size_t nSpectra = outWS->getNumberHistograms();
size_t length = outWS->blocksize();
g_log.information() << "W2D has " << nSpectra << " histograms of length " << length;
for (size_t i = 0; i < nSpectra; i++) {
auto &spec = outWS->getSpectrum(i);
auto &x = spec.x();
auto &y = spec.y();
auto &e = spec.e();
if (x.size() != length) {
g_log.information() << "Spectrum " << i << " x-size mismatch, is " << x.size() << " should be " << length << "\n";
}
if (y.size() != length) {
g_log.information() << "Spectrum " << i << " y-size mismatch, is " << y.size() << " should be " << length << "\n";
}
if (e.size() != length) {
g_log.information() << "Spectrum " << i << " e-size mismatch, is " << e.size() << " should be " << length << "\n";
}
}
}
void MDHistoToWorkspace2D::copyMetaData(const Mantid::API::IMDHistoWorkspace_sptr &inWS,
const Mantid::DataObjects::Workspace2D_sptr &outWS) {
if (inWS->getNumExperimentInfo() > 0) {
ExperimentInfo_sptr info = inWS->getExperimentInfo(0);
outWS->copyExperimentInfoFrom(info.get());
}
outWS->setTitle(inWS->getTitle());
}