Mantid
Loading...
Searching...
No Matches
SpectrumInfo.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 +
10#include "MantidBeamline/SpectrumInfo.h"
15#include "MantidKernel/Logger.h"
17#include "MantidTypes/SpectrumDefinition.h"
18
19#include <algorithm>
20#include <limits>
21#include <memory>
22
23using namespace Mantid::Kernel;
24
25namespace Mantid::API {
27Kernel::Logger g_log("ExperimentInfo");
28
29SpectrumInfo::SpectrumInfo(const Beamline::SpectrumInfo &spectrumInfo, const ExperimentInfo &experimentInfo,
30 Geometry::DetectorInfo &detectorInfo)
31 : m_experimentInfo(experimentInfo), m_detectorInfo(detectorInfo), m_spectrumInfo(spectrumInfo),
32 m_lastDetector(PARALLEL_GET_MAX_THREADS), m_lastIndex(PARALLEL_GET_MAX_THREADS, -1) {}
33
34// Defined as default in source for forward declaration with std::unique_ptr.
36
38size_t SpectrumInfo::size() const { return m_spectrumInfo.size(); }
39
40size_t SpectrumInfo::detectorCount() const { return m_spectrumInfo.detectorCount(); }
41
43const SpectrumDefinition &SpectrumInfo::spectrumDefinition(const size_t index) const {
45 return m_spectrumInfo.spectrumDefinition(index);
46}
47
49 for (size_t i = 0; i < size(); ++i)
51 return m_spectrumInfo.sharedSpectrumDefinitions();
52}
53
55bool SpectrumInfo::isMonitor(const size_t index) const {
57 return std::all_of(spectrumDefinition.cbegin(), spectrumDefinition.cend(),
58 [this](const std::pair<size_t, size_t> &detIndex) { return m_detectorInfo.isMonitor(detIndex); });
59}
60
62bool SpectrumInfo::isMasked(const size_t index) const {
64 return std::all_of(spectrumDefinition.cbegin(), spectrumDefinition.cend(),
65 [this](const std::pair<size_t, size_t> &detIndex) { return m_detectorInfo.isMasked(detIndex); });
66}
67
73double SpectrumInfo::l2(const size_t index) const {
75 auto l2 = std::accumulate(
76 spectrumDefinition.cbegin(), spectrumDefinition.cend(), 0.0,
77 [this](double x, const std::pair<size_t, size_t> &detIndex) { return x + m_detectorInfo.l2(detIndex); });
78 return l2 / static_cast<double>(spectrumDefinition.size());
79}
80
86double SpectrumInfo::twoTheta(const size_t index) const {
88 auto twoTheta = std::accumulate(
89 spectrumDefinition.cbegin(), spectrumDefinition.cend(), 0.0,
90 [this](double x, const std::pair<size_t, size_t> &detIndex) { return x + m_detectorInfo.twoTheta(detIndex); });
91 return twoTheta / static_cast<double>(spectrumDefinition.size());
92}
93
99double SpectrumInfo::signedTwoTheta(const size_t index) const {
101 auto signedTwoTheta = std::accumulate(spectrumDefinition.cbegin(), spectrumDefinition.cend(), 0.0,
102 [this](double x, const std::pair<size_t, size_t> &detIndex) {
103 return x + m_detectorInfo.signedTwoTheta(detIndex);
104 });
105 return signedTwoTheta / static_cast<double>(spectrumDefinition.size());
106}
107
113double SpectrumInfo::azimuthal(const size_t index) const {
115 auto phi = std::accumulate(
116 spectrumDefinition.cbegin(), spectrumDefinition.cend(), 0.0,
117 [this](double x, const std::pair<size_t, size_t> &detIndex) { return x + m_detectorInfo.azimuthal(detIndex); });
118 return phi / static_cast<double>(spectrumDefinition.size());
119}
120
125std::pair<double, double> SpectrumInfo::geographicalAngles(const size_t index) const {
126 double lat{0.0}, lon{0.0};
127 for (const auto &detIndex : checkAndGetSpectrumDefinition(index)) {
128 auto latlong = m_detectorInfo.geographicalAngles(detIndex);
129 lat += latlong.first;
130 lon += latlong.second;
131 }
132 return std::pair<double, double>(lat / static_cast<double>(spectrumDefinition(index).size()),
133 lon / static_cast<double>(spectrumDefinition(index).size()));
134}
135
139 auto newPos = std::accumulate(spectrumDefinition.cbegin(), spectrumDefinition.cend(), Kernel::V3D(),
140 [this](const auto &x, const std::pair<size_t, size_t> &detIndex) {
141 return x + m_detectorInfo.position(detIndex);
142 });
143 return newPos / static_cast<double>(spectrumDefinition.size());
144}
145
155UnitParametersMap SpectrumInfo::diffractometerConstants(const size_t index, std::vector<detid_t> &warningDets) const {
157 throw std::runtime_error("Retrieval of diffractometer constants not "
158 "implemented for scanning instrument");
159 }
160 auto spectrumDef = checkAndGetSpectrumDefinition(index);
161 std::vector<size_t> detectorIndicesOnly;
162 std::vector<detid_t> calibratedDets;
163 std::vector<detid_t> uncalibratedDets;
164 std::transform(spectrumDef.begin(), spectrumDef.end(), std::back_inserter(detectorIndicesOnly),
165 [](auto const &pair) { return pair.first; });
166 double difa{0.}, difc{0.}, tzero{0.};
167 for (const auto &detIndex : detectorIndicesOnly) {
168 auto newDiffConstants = m_detectorInfo.diffractometerConstants(detIndex, calibratedDets, uncalibratedDets);
169 difa += std::get<0>(newDiffConstants);
170 difc += std::get<1>(newDiffConstants);
171 tzero += std::get<2>(newDiffConstants);
172 }
173
174 if (calibratedDets.size() > 0 && uncalibratedDets.size() > 0) {
175 warningDets.insert(warningDets.end(), uncalibratedDets.begin(), uncalibratedDets.end());
176 };
177 // if no calibration is found then return difc only based on the average
178 // of the detector L2 and twoThetas.
179 if (calibratedDets.size() == 0) {
180 return {{UnitParams::difc, difcUncalibrated(index)}};
181 }
182 return {{UnitParams::difa, difa / static_cast<double>(spectrumDefinition(index).size())},
183 {UnitParams::difc, difc / static_cast<double>(spectrumDefinition(index).size())},
184 {UnitParams::tzero, tzero / static_cast<double>(spectrumDefinition(index).size())}};
185}
186
194 std::vector<int> warningDets;
195 return diffractometerConstants(index, warningDets);
196}
197
203double SpectrumInfo::difcUncalibrated(const size_t index) const {
204 // calculate difc based on the average of the detector L2 and twoThetas.
205 // This will be different to the average of the per detector difcs. This is
206 // for backwards compatibility because Mantid always used to calculate
207 // spectrum level difc's this way
209}
210
222void SpectrumInfo::getDetectorValues(const Kernel::Unit &inputUnit, const Kernel::Unit &outputUnit,
223 const Kernel::DeltaEMode::Type emode, const bool signedTheta, int64_t wsIndex,
224 UnitParametersMap &pmap) const {
225 if (!hasDetectors(wsIndex))
226 return;
227 pmap[UnitParams::l2] = l2(wsIndex);
228
229 if (!isMonitor(wsIndex)) {
230 // The scattering angle for this detector (in radians).
231 try {
232 if (signedTheta)
233 pmap[UnitParams::twoTheta] = signedTwoTheta(wsIndex);
234 else
235 pmap[UnitParams::twoTheta] = twoTheta(wsIndex);
236 } catch (const std::runtime_error &e) {
237 g_log.warning(e.what());
238 }
239 if (emode != Kernel::DeltaEMode::Elastic && pmap.find(UnitParams::efixed) == pmap.end()) {
240 std::shared_ptr<const Geometry::IDetector> det(&detector(wsIndex), Mantid::NoDeleting());
241 try {
242 pmap[UnitParams::efixed] = m_experimentInfo.getEFixedGivenEMode(det, emode);
243 g_log.debug() << "Detector: " << det->getID() << " EFixed: " << pmap[UnitParams::efixed] << "\n";
244 } catch (std::runtime_error &) {
245 // let the unit classes work out if this is a problem
246 }
247 }
248
249 std::vector<detid_t> warnDetIds;
250 try {
251 std::set<std::string> diffConstUnits = {"dSpacing", "MomentumTransfer", "Empty"};
252 if ((emode == Kernel::DeltaEMode::Elastic) &&
253 (diffConstUnits.count(inputUnit.unitID()) || diffConstUnits.count(outputUnit.unitID()))) {
254 auto diffConstsMap = diffractometerConstants(wsIndex, warnDetIds);
255 pmap.insert(diffConstsMap.begin(), diffConstsMap.end());
256 if (warnDetIds.size() > 0) {
257 createDetectorIdLogMessages(warnDetIds, wsIndex);
258 }
259 } else {
260 pmap[UnitParams::difc] = difcUncalibrated(wsIndex);
261 }
262 } catch (const std::runtime_error &e) {
263 g_log.warning(e.what());
264 }
265 } else {
266 pmap[UnitParams::twoTheta] = 0.0;
267 pmap[UnitParams::efixed] = std::numeric_limits<double>::min();
268 // Energy transfer is meaningless for a monitor, so set l2 to 0.
269 if (outputUnit.unitID().find("DeltaE") != std::string::npos) {
270 pmap[UnitParams::l2] = 0.0;
271 }
272 pmap[UnitParams::difc] = 0;
273 }
274}
275
276void SpectrumInfo::createDetectorIdLogMessages(const std::vector<detid_t> &detids, int64_t wsIndex) const {
277 std::string detIDstring;
278 auto iter = detids.begin();
279 auto itEnd = detids.end();
280 for (; iter != itEnd; ++iter) {
281 detIDstring += std::to_string(*iter) + ",";
282 }
283
284 if (!detIDstring.empty()) {
285 detIDstring.pop_back(); // Drop last comma
286 }
287 g_log.warning("Incomplete set of calibrated diffractometer constants found for "
288 "workspace index" +
289 std::to_string(wsIndex) + ". Using uncalibrated values for detectors " + detIDstring);
290}
291
294bool SpectrumInfo::hasDetectors(const size_t index) const {
295 // Workspaces can contain invalid detector IDs. Those IDs will be silently
296 // ignored here until this is fixed.
297 return spectrumDefinition(index).size() > 0;
298}
299
301bool SpectrumInfo::hasUniqueDetector(const size_t index) const {
302 // Workspaces can contain invalid detector IDs. Those IDs will be silently
303 // ignored here until this is fixed.
304 return spectrumDefinition(index).size() == 1;
305}
306
310void SpectrumInfo::setMasked(const size_t index, bool masked) {
311 for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
312 m_detectorInfo.setMasked(detIndex, masked);
313}
314
317const Geometry::IDetector &SpectrumInfo::detector(const size_t index) const { return getDetector(index); }
318
321
324
326double SpectrumInfo::l1() const { return m_detectorInfo.l1(); }
327
329 auto thread = static_cast<size_t>(PARALLEL_THREAD_NUMBER);
330 if (m_lastIndex[thread] == index)
331 return *m_lastDetector[thread];
332
333 // Note: This function body has big overlap with the method
334 // MatrixWorkspace::getDetector(). The plan is to eventually remove the
335 // latter, once SpectrumInfo is in widespread use.
336 const auto &specDef = spectrumDefinition(index);
337 const size_t ndets = specDef.size();
338 if (ndets == 1) {
339 // If only 1 detector for the spectrum number, just return it
340 const auto detIndex = specDef[0].first;
341 m_lastDetector[thread] = m_detectorInfo.getDetectorPtr(detIndex);
342 } else if (ndets == 0) {
343 throw Kernel::Exception::NotFoundError("MatrixWorkspace::getDetector(): No "
344 "detectors for this workspace "
345 "index.",
347 } else {
348 // Else need to construct a DetectorGroup and use that
349 std::vector<std::shared_ptr<const Geometry::IDetector>> det_ptrs;
350 for (const auto &specDefIndex : specDef) {
351 const auto detIndex = specDefIndex.first;
352 det_ptrs.emplace_back(m_detectorInfo.getDetectorPtr(detIndex));
353 }
354 m_lastDetector[thread] = std::make_shared<Geometry::DetectorGroup>(det_ptrs);
355 }
356 m_lastIndex[thread] = index;
357 return *m_lastDetector[thread];
358}
359
360const SpectrumDefinition &SpectrumInfo::checkAndGetSpectrumDefinition(const size_t index) const {
361 if (spectrumDefinition(index).size() == 0)
362 throw Kernel::Exception::NotFoundError("SpectrumInfo: No detectors for this workspace index.",
365}
366
367// Begin method for iterator
369
370// End method for iterator
372
374
375// End method for iterator
377
378} // namespace Mantid::API
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define PARALLEL_THREAD_NUMBER
#define PARALLEL_GET_MAX_THREADS
This class is shared by a few Workspace types and holds information related to a particular experimen...
void updateSpectrumDefinitionIfNecessary(const size_t index) const
double getEFixedGivenEMode(const std::shared_ptr< const Geometry::IDetector > &detector, const Kernel::DeltaEMode::Type emode) const
Easy access to the efixed value for this run & detector.
double signedTwoTheta(const size_t index) const
Returns the signed scattering angle 2 theta in radians (angle w.r.t.
SpectrumInfoIterator< SpectrumInfo > end()
bool isMonitor(const size_t index) const
Returns true if the detector(s) associated with the spectrum are monitors.
const Kernel::cow_ptr< std::vector< SpectrumDefinition > > & sharedSpectrumDefinitions() const
Kernel::V3D sourcePosition() const
Returns the source position.
const SpectrumInfoIterator< const SpectrumInfo > cbegin() const
Kernel::UnitParametersMap diffractometerConstants(const size_t index, std::vector< detid_t > &uncalibratedDets) const
Calculate average diffractometer constants (DIFA, DIFC, TZERO) of detectors associated with this spec...
void setMasked(const size_t index, bool masked)
Set the mask flag of the spectrum with given index.
void createDetectorIdLogMessages(const std::vector< detid_t > &detids, int64_t wsIndex) const
Kernel::V3D samplePosition() const
Returns the sample position.
bool hasDetectors(const size_t index) const
Returns true if the spectrum is associated with detectors in the instrument.
double twoTheta(const size_t index) const
Returns the scattering angle 2 theta in radians (angle w.r.t.
double difcUncalibrated(const size_t index) const
Calculate average uncalibrated DIFC value of detectors associated with this spectrum.
Kernel::V3D position(const size_t index) const
Returns the position of the spectrum with given index.
double l2(const size_t index) const
Returns L2 (distance from sample to spectrum).
double azimuthal(const size_t index) const
Returns the out-of-plane angle in radians (angle w.r.t.
bool isMasked(const size_t index) const
Returns true if the detector(s) associated with the spectrum are masked.
const SpectrumDefinition & checkAndGetSpectrumDefinition(const size_t index) const
SpectrumInfo(const Beamline::SpectrumInfo &spectrumInfo, const ExperimentInfo &experimentInfo, Geometry::DetectorInfo &detectorInfo)
bool hasUniqueDetector(const size_t index) const
Returns true if the spectrum is associated with exactly one detector.
const ExperimentInfo & m_experimentInfo
Definition: SpectrumInfo.h:108
const SpectrumDefinition & spectrumDefinition(const size_t index) const
Returns a const reference to the SpectrumDefinition of the spectrum.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector or detector group of the spectrum with given index.
size_t size() const
Returns the size of the SpectrumInfo, i.e., the number of spectra.
double l1() const
Returns L1 (distance from source to sample).
const SpectrumInfoIterator< const SpectrumInfo > cend() const
void getDetectorValues(const Kernel::Unit &inputUnit, const Kernel::Unit &outputUnit, const Kernel::DeltaEMode::Type emode, const bool signedTheta, int64_t wsIndex, Kernel::UnitParametersMap &pmap) const
Get the detector values relevant to unit conversion for a workspace index.
Geometry::DetectorInfo & m_detectorInfo
Definition: SpectrumInfo.h:109
size_t detectorCount() const
SpectrumInfoIterator< SpectrumInfo > begin()
std::pair< double, double > geographicalAngles(const size_t index) const
Calculate latitude and longitude for given spectrum index.
std::vector< std::shared_ptr< const Geometry::IDetector > > m_lastDetector
Definition: SpectrumInfo.h:111
std::vector< size_t > m_lastIndex
Definition: SpectrumInfo.h:112
const Geometry::IDetector & getDetector(const size_t index) const
const Beamline::SpectrumInfo & m_spectrumInfo
Definition: SpectrumInfo.h:110
Geometry::DetectorInfo is an intermediate step towards a DetectorInfo that is part of Instrument-2....
Definition: DetectorInfo.h:49
Kernel::V3D samplePosition() const
Returns the sample position.
std::pair< double, double > geographicalAngles(const size_t index) const
void setMasked(const size_t index, bool masked)
Set the mask flag of the detector with given index. Not thread safe.
double l1() const
Returns L1 (distance from source to sample).
Kernel::V3D sourcePosition() const
Returns the source position.
std::tuple< double, double, double > diffractometerConstants(const size_t index, std::vector< detid_t > &calibratedDets, std::vector< detid_t > &uncalibratedDets) const
bool isScanning() const
Returns true if the beamline has scanning detectors.
std::shared_ptr< const Geometry::IDetector > getDetectorPtr(const size_t index) const
Helper used by SpectrumInfo.
Interface class for detector objects.
Definition: IDetector.h:43
Exception for when an item is not found in a collection.
Definition: Exception.h:145
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
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
The base units (abstract) class.
Definition: Unit.h:41
virtual const std::string unitID() const =0
The name of the unit.
Class for 3D vectors.
Definition: V3D.h:34
Implements a copy on write data template.
Definition: cow_ptr.h:41
An object for constructing a shared_ptr that won't ever delete its pointee.
Definition: IComponent.h:166
Kernel::Logger g_log("ExperimentInfo")
static logger object
SpectrumInfoIterator< const SpectrumInfo > SpectrumInfoConstIt
Definition: SpectrumInfo.h:116
SpectrumInfoIterator< SpectrumInfo > SpectrumInfoIt
Definition: SpectrumInfo.h:115
MANTID_KERNEL_DLL double tofToDSpacingFactor(const double l1, const double l2, const double twoTheta, const double offset)
Calculate and return conversion factor from tof to d-spacing.
Definition: Unit.cpp:568
std::unordered_map< UnitParams, double > UnitParametersMap
Definition: Unit.h:30
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
Type
Define the available energy transfer modes It is important to assign enums proper numbers,...
Definition: DeltaEMode.h:29