19#include "MantidHistogramData/HistogramIterator.h"
22#include <Poco/DOM/AutoPtr.h>
23#include <Poco/DOM/Document.h>
33 for (
size_t i = 1; i < detIDs.size(); ++i) {
41constexpr double R = 1.0;
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");
59 const size_t numSpectra = rows * columns;
64 auto instrument = std::make_shared<Geometry::Instrument>(
"MC_simulation_instrument");
65 const auto refFrame = modelWS.
getInstrument()->getReferenceFrame();
67 instrument->setReferenceFrame(std::make_shared<Geometry::ReferenceFrame>(*refFrame));
70 auto sample = std::make_unique<Geometry::Component>(
"sample", instrument.get());
72 instrument->add(
sample.get());
73 instrument->markAsSamplePos(
sample.release());
77 p[refFrame->pointingAlongBeam()] = -2.0 * R;
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());
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());
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);
104 instrument->add(det.get());
105 instrument->markAsDetector(det.release());
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]);
125 const auto beamRadiusParam = modelSource->getNumberParameter(
"beam-radius");
126 if (beamRadiusParam.size() == 1) {
127 paramMap.add(
"double", parametrizedSource.get(),
"beam-radius", beamRadiusParam[0]);
130 const auto eMode = modelWS.
getEMode();
136 if (!constantIndirectEFixed(modelWS, detIDs)) {
137 throw std::runtime_error(
"Sparse instrument with variable EFixed not supported.");
141 for (
int sparseDetID : sparseDetIDs) {
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();
180 return std::make_tuple(minLat, maxLat, minLong, maxLong);
188 double currentMin = std::numeric_limits<double>::max();
189 double currentMax = std::numeric_limits<double>::lowest();
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);
197 return std::make_tuple(currentMin, currentMax);
206 const size_t wavelengthPoints) {
207 double minWavelength, maxWavelength;
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);
221 xs.back() = maxWavelength;
223 xs.front() = (minWavelength + maxWavelength) / 2.0;
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");
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);
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));
289 std::array<double, 4> weights;
290 for (
size_t i = 0; i < weights.size(); ++i) {
291 if (distances[i] == 0.0) {
296 weights[i] = 1.0 / distances[i] / distances[i];
307 const double distanceStep)
const {
308 HistogramData::HistogramY avgSecondDeriv(
blocksize());
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;
324 const auto indices =
m_gridDef->nearestNeighbourIndices(lat, lon);
328 std::array<double, 4> distances;
329 for (
size_t i = 0; i < 4; ++i) {
330 double detLat, detLong;
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]);
341 h.mutableY() /= weightSum;
356 std::transform(
e.cbegin(),
e.cend(),
e.begin(), [](
double f) ->
double { return sqrt(f); });
370 const double lon)
const {
372 size_t nearestLatIndex, nearestLonIndex;
373 std::tie(nearestLatIndex, nearestLonIndex) =
m_gridDef->getNearestVertex(lat, lon);
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);
382 double latLow, longLow, latHigh, longHigh;
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);
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);
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);
400 HistogramData::HistogramE interpolationError(
e(0).
size(), 0);
404 auto nearestLatIndexSec = nearestLatIndex > 0 ? nearestLatIndex - 1 : 0;
405 auto nearestLonIndexSec = nearestLonIndex > 0 ? nearestLonIndex - 1 : 0;
406 std::array<size_t, 3> threeIndices;
409 for (
int i = 0; i < 3; i++) {
410 threeIndices[i] =
m_gridDef->getDetectorIndex(nearestLatIndex, nearestLonIndexSec + i);
412 auto avgSecondDerivLong =
secondDerivative(threeIndices, longHigh - longLow);
415 for (
int i = 0; i < 3; i++) {
416 threeIndices[i] =
m_gridDef->getDetectorIndex(nearestLatIndexSec + i, nearestLonIndex);
423 HistogramData::HistogramY interpolationErrorAsYData =
424 0.5 * (lon - longLow) * (longHigh - lon) * avgSecondDerivLong +
425 0.5 * (lat - latLow) * (latHigh - lat) * avgSecondDerivLat;
427 interpolationError = interpolationErrorAsYData.rawData();
431 h.mutableY() = interpY;
435 h.mutableE() =
esqrt(
esq(interpolationError) +
esq(errFromSourcePoints));
std::map< DeltaEMode::Type, std::string > index
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.
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
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.
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.
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.
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.
void warning(const std::string &msg)
Logs at warning level.
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< IObject > IObject_sptr
Typdef for a shared pointer.
static constexpr double h
Planck constant in J*s.
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.