Mantid
Loading...
Searching...
No Matches
SparseWorkspace.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 +
8
10#include "MantidAPI/Run.h"
19#include "MantidHistogramData/HistogramIterator.h"
20#include "MantidKernel/Logger.h"
21
22#include <Poco/DOM/AutoPtr.h>
23#include <Poco/DOM/Document.h>
24
25namespace {
31bool constantIndirectEFixed(const Mantid::API::ExperimentInfo &info, const std::vector<Mantid::detid_t> &detIDs) {
32 const auto e = info.getEFixed(detIDs[0]);
33 for (size_t i = 1; i < detIDs.size(); ++i) {
34 if (e != info.getEFixed(detIDs[i])) {
35 return false;
36 }
37 }
38 return true;
39}
40
41constexpr double R = 1.0; // This will be the default L2 distance.
42
44Mantid::Kernel::Logger g_log("SparseWorkspace");
45} // namespace
46
47namespace Mantid::Algorithms {
48
49SparseWorkspace::SparseWorkspace(const API::MatrixWorkspace &modelWS, const size_t wavelengthPoints, const size_t rows,
50 const size_t columns)
51 : Workspace2D() {
52 double minLat, maxLat, minLong, maxLong;
53 std::tie(minLat, maxLat, minLong, maxLong) = extremeAngles(modelWS);
54 m_gridDef = std::make_unique<Algorithms::DetectorGridDefinition>(minLat, maxLat, rows, minLong, maxLong, columns);
55 if ((rows < 3) || (columns < 3)) {
56 g_log.warning("Can't calculate errors on a sparse workspace with lat or "
57 "long dimension < 3");
58 }
59 const size_t numSpectra = rows * columns;
60 const auto h = modelHistogram(modelWS, wavelengthPoints);
61 initialize(numSpectra, h);
62
63 // Build a quite standard and somewhat complete instrument.
64 auto instrument = std::make_shared<Geometry::Instrument>("MC_simulation_instrument");
65 const auto refFrame = modelWS.getInstrument()->getReferenceFrame();
66
67 instrument->setReferenceFrame(std::make_shared<Geometry::ReferenceFrame>(*refFrame));
68 // The sparse instrument is build around origin.
69 constexpr Kernel::V3D samplePos{0.0, 0.0, 0.0};
70 auto sample = std::make_unique<Geometry::Component>("sample", instrument.get());
71 sample->setPos(samplePos);
72 instrument->add(sample.get());
73 instrument->markAsSamplePos(sample.release());
74 // Add source behind the sample.
75 const Kernel::V3D sourcePos = [&]() {
77 p[refFrame->pointingAlongBeam()] = -2.0 * R;
78 return p;
79 }();
80 auto source = std::make_unique<Geometry::ObjComponent>("source", nullptr, instrument.get());
81 source->setPos(sourcePos);
82 instrument->add(source.get());
83 instrument->markAsSource(source.release());
84
85 auto detShape = makeCubeShape();
86 for (size_t col = 0; col < columns; ++col) {
87 const auto lon = m_gridDef->longitudeAt(col);
88 for (size_t row = 0; row < rows; ++row) {
89 const auto lat = m_gridDef->latitudeAt(row);
90 const size_t index = col * rows + row;
91 const auto detID = static_cast<int>(index);
92 std::ostringstream detName;
93 detName << "det-" << detID;
94 auto det = std::make_unique<Geometry::Detector>(detName.str(), detID, detShape, instrument.get());
95 const Kernel::V3D pos = [&]() {
97 p[refFrame->pointingHorizontal()] = R * std::sin(lon) * std::cos(lat);
98 p[refFrame->pointingUp()] = R * std::sin(lat);
99 p[refFrame->pointingAlongBeam()] = R * std::cos(lon) * std::cos(lat);
100 return p;
101 }();
102 det->setPos(pos);
104 instrument->add(det.get());
105 instrument->markAsDetector(det.release());
106 }
107 }
108 setInstrument(instrument);
109
110 // Copy things needed for the simulation from the model workspace.
111 auto &paramMap = instrumentParameters();
112 auto parametrizedInstrument = getInstrument();
113 // Copy beam parameters.
114 const auto modelSource = modelWS.getInstrument()->getSource();
115 const auto parametrizedSource = parametrizedInstrument->getSource();
116 const auto beamShapeParam = modelSource->getParameterAsString("beam-shape");
117 if (!beamShapeParam.empty())
118 paramMap.add("string", parametrizedSource.get(), "beam-shape", beamShapeParam);
119 const auto beamWidthParam = modelSource->getNumberParameter("beam-width");
120 const auto beamHeightParam = modelSource->getNumberParameter("beam-height");
121 if (beamWidthParam.size() == 1 && beamHeightParam.size() == 1) {
122 paramMap.add("double", parametrizedSource.get(), "beam-width", beamWidthParam[0]);
123 paramMap.add("double", parametrizedSource.get(), "beam-height", beamHeightParam[0]);
124 }
125 const auto beamRadiusParam = modelSource->getNumberParameter("beam-radius");
126 if (beamRadiusParam.size() == 1) {
127 paramMap.add("double", parametrizedSource.get(), "beam-radius", beamRadiusParam[0]);
128 }
129 // Add information about EFixed in a proper place.
130 const auto eMode = modelWS.getEMode();
131 mutableRun().addProperty("deltaE-mode", Kernel::DeltaEMode::asString(eMode));
132 if (eMode == Kernel::DeltaEMode::Direct) {
133 mutableRun().addProperty("Ei", modelWS.getEFixed());
134 } else if (eMode == Kernel::DeltaEMode::Indirect) {
135 const auto &detIDs = modelWS.detectorInfo().detectorIDs();
136 if (!constantIndirectEFixed(modelWS, detIDs)) {
137 throw std::runtime_error("Sparse instrument with variable EFixed not supported.");
138 }
139 const auto e = modelWS.getEFixed(detIDs[0]);
140 const auto &sparseDetIDs = detectorInfo().detectorIDs();
141 for (int sparseDetID : sparseDetIDs) {
142 setEFixed(sparseDetID, e);
143 }
144 }
145}
146
148 : Workspace2D(other), m_gridDef(std::make_unique<Algorithms::DetectorGridDefinition>(*other.m_gridDef)) {}
149
156std::tuple<double, double, double, double> SparseWorkspace::extremeAngles(const API::MatrixWorkspace &ws) {
157 const auto &spectrumInfo = ws.spectrumInfo();
158 double minLat = std::numeric_limits<double>::max();
159 double maxLat = std::numeric_limits<double>::lowest();
160 double minLong = std::numeric_limits<double>::max();
161 double maxLong = std::numeric_limits<double>::lowest();
162 for (size_t i = 0; i < ws.getNumberHistograms(); ++i) {
163 if (spectrumInfo.hasDetectors(i)) {
164 double lat, lon;
165 std::tie(lat, lon) = spectrumInfo.geographicalAngles(i);
166 if (lat < minLat) {
167 minLat = lat;
168 }
169 if (lat > maxLat) {
170 maxLat = lat;
171 }
172 if (lon < minLong) {
173 minLong = lon;
174 }
175 if (lon > maxLong) {
176 maxLong = lon;
177 }
178 }
179 }
180 return std::make_tuple(minLat, maxLat, minLong, maxLong);
181}
182
187std::tuple<double, double> SparseWorkspace::extremeWavelengths(const API::MatrixWorkspace &ws) {
188 double currentMin = std::numeric_limits<double>::max();
189 double currentMax = std::numeric_limits<double>::lowest();
190 for (size_t i = 0; i < ws.getNumberHistograms(); ++i) {
191 const auto h = ws.histogram(i);
192 const auto first = h.begin();
193 currentMin = std::min(first->center(), currentMin);
194 const auto last = std::prev(h.end());
195 currentMax = std::max(last->center(), currentMax);
196 }
197 return std::make_tuple(currentMin, currentMax);
198}
199
205Mantid::HistogramData::Histogram SparseWorkspace::modelHistogram(const API::MatrixWorkspace &modelWS,
206 const size_t wavelengthPoints) {
207 double minWavelength, maxWavelength;
208 std::tie(minWavelength, maxWavelength) = extremeWavelengths(modelWS);
209 HistogramData::Frequencies ys(wavelengthPoints, 0.0);
210 HistogramData::FrequencyVariances es(wavelengthPoints, 0.0);
211 HistogramData::Points ps(wavelengthPoints, 0.0);
212 HistogramData::Histogram h(ps, ys, es);
213 auto &xs = h.mutableX();
214 if (wavelengthPoints > 1) {
215 const double step = (maxWavelength - minWavelength) / static_cast<double>(wavelengthPoints - 1);
216 for (size_t i = 0; i < xs.size() - 1; ++i) {
217 xs[i] = minWavelength + step * static_cast<double>(i);
218 }
219 // Force last point as otherwise it might be slightly off due to
220 // small rounding errors in the calculation above.
221 xs.back() = maxWavelength;
222 } else {
223 xs.front() = (minWavelength + maxWavelength) / 2.0;
224 }
225 return h;
226}
227
232 using namespace Poco::XML;
233 const double dimension = 0.05;
234 AutoPtr<Document> shapeDescription = new Document;
235 AutoPtr<Element> typeElement = shapeDescription->createElement("type");
236 typeElement->setAttribute("name", "detector");
237 AutoPtr<Element> shapeElement = shapeDescription->createElement("cuboid");
238 shapeElement->setAttribute("id", "cube");
239 const std::string posCoord = std::to_string(dimension / 2);
240 const std::string negCoord = std::to_string(-dimension / 2);
241 AutoPtr<Element> element = shapeDescription->createElement("left-front-bottom-point");
242 element->setAttribute("x", negCoord);
243 element->setAttribute("y", negCoord);
244 element->setAttribute("z", posCoord);
245 shapeElement->appendChild(element);
246 element = shapeDescription->createElement("left-front-top-point");
247 element->setAttribute("x", negCoord);
248 element->setAttribute("y", posCoord);
249 element->setAttribute("z", posCoord);
250 shapeElement->appendChild(element);
251 element = shapeDescription->createElement("left-back-bottom-point");
252 element->setAttribute("x", negCoord);
253 element->setAttribute("y", negCoord);
254 element->setAttribute("z", negCoord);
255 shapeElement->appendChild(element);
256 element = shapeDescription->createElement("right-front-bottom-point");
257 element->setAttribute("x", posCoord);
258 element->setAttribute("y", negCoord);
259 element->setAttribute("z", posCoord);
260 shapeElement->appendChild(element);
261 typeElement->appendChild(shapeElement);
262 AutoPtr<Element> algebraElement = shapeDescription->createElement("algebra");
263 algebraElement->setAttribute("val", "cube");
264 typeElement->appendChild(algebraElement);
265 Geometry::ShapeFactory shapeFactory;
266 return shapeFactory.createShape(typeElement);
267}
268
276double SparseWorkspace::greatCircleDistance(const double lat1, const double long1, const double lat2,
277 const double long2) {
278 const double latD = std::sin((lat2 - lat1) / 2.0);
279 const double longD = std::sin((long2 - long1) / 2.0);
280 const double S = latD * latD + std::cos(lat1) * std::cos(lat2) * longD * longD;
281 return 2.0 * std::asin(std::sqrt(S));
282}
283
288std::array<double, 4> SparseWorkspace::inverseDistanceWeights(const std::array<double, 4> &distances) {
289 std::array<double, 4> weights;
290 for (size_t i = 0; i < weights.size(); ++i) {
291 if (distances[i] == 0.0) {
292 weights.fill(0.0);
293 weights[i] = 1.0;
294 return weights;
295 }
296 weights[i] = 1.0 / distances[i] / distances[i];
297 }
298 return weights;
299}
300
306HistogramData::HistogramY SparseWorkspace::secondDerivative(const std::array<size_t, 3> indices,
307 const double distanceStep) const {
308 HistogramData::HistogramY avgSecondDeriv(blocksize());
309
310 HistogramData::HistogramY sumSecondDeriv(blocksize());
311 auto firstDerivB = (y(indices[2]) - y(indices[1])) / distanceStep;
312 auto firstDerivA = (y(indices[1]) - y(indices[0])) / distanceStep;
313 auto secondDerivA = (firstDerivB - firstDerivA) / distanceStep;
314
315 return secondDerivA;
316}
317
323HistogramData::Histogram SparseWorkspace::interpolateFromDetectorGrid(const double lat, const double lon) const {
324 const auto indices = m_gridDef->nearestNeighbourIndices(lat, lon);
325
326 auto h = histogram(0);
327
328 std::array<double, 4> distances;
329 for (size_t i = 0; i < 4; ++i) {
330 double detLat, detLong;
331 std::tie(detLat, detLong) = spectrumInfo().geographicalAngles(indices[i]);
332 distances[i] = greatCircleDistance(lat, lon, detLat, detLong);
333 }
334 const auto weights = inverseDistanceWeights(distances);
335 auto weightSum = weights[0];
336 h.mutableY() = weights[0] * y(indices[0]);
337 for (size_t i = 1; i < 4; ++i) {
338 weightSum += weights[i];
339 h.mutableY() += weights[i] * y(indices[i]);
340 }
341 h.mutableY() /= weightSum;
342 return h;
343}
344
349HistogramData::HistogramE SparseWorkspace::esq(const HistogramData::HistogramE &e) const { return e * e; }
350
355HistogramData::HistogramE SparseWorkspace::esqrt(HistogramData::HistogramE e) const {
356 std::transform(e.cbegin(), e.cend(), e.begin(), [](double f) -> double { return sqrt(f); });
357 return e;
358}
359
369HistogramData::Histogram SparseWorkspace::bilinearInterpolateFromDetectorGrid(const double lat,
370 const double lon) const {
371
372 size_t nearestLatIndex, nearestLonIndex;
373 std::tie(nearestLatIndex, nearestLonIndex) = m_gridDef->getNearestVertex(lat, lon);
374
375 std::array<std::array<size_t, 2>, 2> detIndices;
376 for (size_t i = 0; i < 2; i++) {
377 for (size_t j = 0; j < 2; j++) {
378 detIndices[i][j] = m_gridDef->getDetectorIndex(nearestLatIndex + j, nearestLonIndex + i);
379 }
380 }
381
382 double latLow, longLow, latHigh, longHigh;
383 std::tie(latLow, longLow) = spectrumInfo().geographicalAngles(detIndices[0][0]);
384 std::tie(latHigh, longHigh) = spectrumInfo().geographicalAngles(detIndices[1][1]);
385
386 // interpolate across different longitudes
387 auto ylat1 = ((longHigh - lon) * y(detIndices[0][0]) + (lon - longLow) * y(detIndices[1][0])) / (longHigh - longLow);
388 auto errLat1 = esqrt(esq((longHigh - lon) * e(detIndices[0][0])) + esq((lon - longLow) * e(detIndices[1][0]))) /
389 (longHigh - longLow);
390
391 auto ylat2 = ((longHigh - lon) * y(detIndices[0][1]) + (lon - longLow) * y(detIndices[1][1])) / (longHigh - longLow);
392 auto errLat2 = esqrt(esq((longHigh - lon) * e(detIndices[0][1])) + esq((lon - longLow) * e(detIndices[1][1]))) /
393 (longHigh - longLow);
394
395 // interpolate across different latitudes
396 auto interpY = ((latHigh - lat) * ylat1 + (lat - latLow) * ylat2) / (latHigh - latLow);
397 auto errFromSourcePoints = esqrt(esq((latHigh - lat) * errLat1) + esq((lat - latLow) * errLat2)) / (latHigh - latLow);
398
399 // calculate interpolation errors if possible
400 HistogramData::HistogramE interpolationError(e(0).size(), 0);
401 if ((m_gridDef->numberColumns() > 2) && (m_gridDef->numberRows() > 2)) {
402 // step back a further row\col (if possible) to estimate the second
403 // derivative
404 auto nearestLatIndexSec = nearestLatIndex > 0 ? nearestLatIndex - 1 : 0;
405 auto nearestLonIndexSec = nearestLonIndex > 0 ? nearestLonIndex - 1 : 0;
406 std::array<size_t, 3> threeIndices;
407
408 // 2nd derivative in longitude
409 for (int i = 0; i < 3; i++) {
410 threeIndices[i] = m_gridDef->getDetectorIndex(nearestLatIndex, nearestLonIndexSec + i);
411 }
412 auto avgSecondDerivLong = secondDerivative(threeIndices, longHigh - longLow);
413
414 // 2nd derivative in latitude
415 for (int i = 0; i < 3; i++) {
416 threeIndices[i] = m_gridDef->getDetectorIndex(nearestLatIndexSec + i, nearestLonIndex);
417 }
418 auto avgSecondDerivLat = secondDerivative(threeIndices, latHigh - latLow);
419
420 // calculate the interpolation error according to a Taylor expansion from
421 // the low lat and low long points. Doesn't make any difference if do this
422 // in angle or distance because scaling factor (~ radius) cancels out
423 HistogramData::HistogramY interpolationErrorAsYData =
424 0.5 * (lon - longLow) * (longHigh - lon) * avgSecondDerivLong +
425 0.5 * (lat - latLow) * (latHigh - lat) * avgSecondDerivLat;
426 // convert from HistogramY to HistogramE
427 interpolationError = interpolationErrorAsYData.rawData();
428 }
429
430 auto h = histogram(0);
431 h.mutableY() = interpY;
432 // combine the (independent) errors
433 // Note - interpolationError may be negative if the 2nd derivative is negative
434 // so the following line also has useful side effect of taking abs value
435 h.mutableE() = esqrt(esq(interpolationError) + esq(errFromSourcePoints));
436 return h;
437}
438
440
441} // namespace Mantid::Algorithms
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
This class is shared by a few Workspace types and holds information related to a particular experimen...
Run & mutableRun()
Writable version of the run object.
const SpectrumInfo & spectrumInfo() const
Return a reference to the SpectrumInfo object.
const Geometry::DetectorInfo & detectorInfo() const
Return a const reference to the DetectorInfo object.
double getEFixed(const detid_t detID) const
Easy access to the efixed value for this run & detector ID.
Geometry::Instrument_const_sptr getInstrument() const
Returns the parameterized instrument.
Kernel::DeltaEMode::Type getEMode() const
Returns the emode for this run.
const Sample & sample() const
Sample accessors.
void setEFixed(const detid_t detID, const double value)
Set the efixed value for a given detector ID.
const Geometry::ParameterMap & instrumentParameters() const
Returns the set of parameters modifying the base instrument (const-version)
void setInstrument(const Geometry::Instrument_const_sptr &instr)
Instrument accessors.
void setDetectorID(const detid_t detID)
Clear the list of detector IDs, then add one.
Definition: ISpectrum.cpp:84
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
Definition: LogManager.h:79
Base MatrixWorkspace Abstract Class.
const HistogramData::HistogramE & e(const size_t index) const
void initialize(const std::size_t &NVectors, const std::size_t &XLength, const std::size_t &YLength)
Initialize the workspace.
virtual std::size_t getNumberHistograms() const =0
Returns the number of histograms in the workspace.
HistogramData::Histogram histogram(const size_t index) const
Returns the Histogram at the given workspace index.
bool hasDetectors(const size_t index) const
Returns true if the spectrum is associated with detectors in the instrument.
std::pair< double, double > geographicalAngles(const size_t index) const
Calculate latitude and longitude for given spectrum index.
DetectorGridDefinition is a helper class for building the sparse instrument in MonteCarloAbsorption.
Defines functions and utilities to create and deal with sparse instruments.
static std::array< double, 4 > inverseDistanceWeights(const std::array< double, 4 > &distances)
Calculate the inverse distance weights for the given distances.
HistogramData::HistogramE esq(const HistogramData::HistogramE &e) const
Square the error values in a histogram.
static std::tuple< double, double > extremeWavelengths(const API::MatrixWorkspace &ws)
Find the maximum and minimum wavelength points over the entire workpace.
static HistogramData::Histogram modelHistogram(const API::MatrixWorkspace &modelWS, const size_t wavelengthPoints)
Create a template histogram for the sparse instrument workspace.
HistogramData::HistogramY secondDerivative(const std::array< size_t, 3 > indices, const double distanceStep) const
Calculate the second derivative of a histogram along a row of indices.
std::unique_ptr< Algorithms::DetectorGridDefinition > m_gridDef
static double greatCircleDistance(const double lat1, const double long1, const double lat2, const double long2)
Calculate the distance between two points on a unit sphere.
HistogramData::HistogramE esqrt(HistogramData::HistogramE e) const
Square root the error values in a histogram.
static std::tuple< double, double, double, double > extremeAngles(const API::MatrixWorkspace &ws)
Find the latitude and longitude intervals the detectors of the given workspace span as seen from the ...
Mantid::Geometry::IObject_sptr makeCubeShape()
Creates a rectangular cuboid shape.
virtual HistogramData::Histogram interpolateFromDetectorGrid(const double lat, const double lon) const
Spatially interpolate a single histogram from nearby detectors.
SparseWorkspace * doClone() const override
Virtual clone method. Not implemented to force implementation in children.
virtual HistogramData::Histogram bilinearInterpolateFromDetectorGrid(const double lat, const double lon) const
Spatially interpolate a single histogram from nearby detectors using bilinear interpolation method.
SparseWorkspace(const API::MatrixWorkspace &modelWS, const size_t wavelengthPoints, const size_t rows, const size_t columns)
Concrete workspace implementation.
Definition: Workspace2D.h:29
std::size_t size() const override
get pseudo size
Histogram1D & getSpectrum(const size_t index) override
Return the underlying ISpectrum ptr at the given workspace index.
Definition: Workspace2D.h:62
std::size_t blocksize() const override
get the size of each vector
const std::vector< detid_t > & detectorIDs() const
Returns a sorted vector of all detector IDs.
Class originally intended to be used with the DataHandling 'LoadInstrument' algorithm.
Definition: ShapeFactory.h:89
std::shared_ptr< CSGObject > createShape(Poco::XML::Element *pElem)
Creates a geometric object from a DOM-element-node pointing to an element whose child nodes contain t...
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
Class for 3D vectors.
Definition: V3D.h:34
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< IObject > IObject_sptr
Typdef for a shared pointer.
Definition: IObject.h:92
static constexpr double h
Planck constant in J*s.
STL namespace.
std::string to_string(const wide_integer< Bits, Signed > &n)
static std::string asString(const Type mode)
Return a string representation of the given mode.
Definition: DeltaEMode.cpp:52