SpectrumInfo.cpp 16.3 KB
Newer Older
1
2
3
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4
5
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6
// SPDX - License - Identifier: GPL - 3.0 +
7
#include "MantidAPI/SpectrumInfo.h"
LamarMoore's avatar
LamarMoore committed
8
#include "MantidAPI/ExperimentInfo.h"
9
#include "MantidAPI/SpectrumInfoIterator.h"
LamarMoore's avatar
LamarMoore committed
10
#include "MantidBeamline/SpectrumInfo.h"
11
#include "MantidGeometry/Instrument/DetectorGroup.h"
12
#include "MantidGeometry/Instrument/DetectorInfo.h"
13
#include "MantidKernel/Exception.h"
Danny Hindson's avatar
Danny Hindson committed
14
#include "MantidKernel/Logger.h"
15
#include "MantidKernel/MultiThreaded.h"
16
#include "MantidTypes/SpectrumDefinition.h"
17

18
#include <algorithm>
19
#include <memory>
20

Danny Hindson's avatar
Danny Hindson committed
21
22
using namespace Mantid::Kernel;

23
24
namespace Mantid {
namespace API {
Danny Hindson's avatar
Danny Hindson committed
25
26
/// static logger object
Kernel::Logger g_log("ExperimentInfo");
27

28
SpectrumInfo::SpectrumInfo(const Beamline::SpectrumInfo &spectrumInfo, const ExperimentInfo &experimentInfo,
29
                           Geometry::DetectorInfo &detectorInfo)
30
31
    : m_experimentInfo(experimentInfo), m_detectorInfo(detectorInfo), m_spectrumInfo(spectrumInfo),
      m_lastDetector(PARALLEL_GET_MAX_THREADS), m_lastIndex(PARALLEL_GET_MAX_THREADS, -1) {}
32
33
34
35

// Defined as default in source for forward declaration with std::unique_ptr.
SpectrumInfo::~SpectrumInfo() = default;

36
/// Returns the size of the SpectrumInfo, i.e., the number of spectra.
37
size_t SpectrumInfo::size() const { return m_spectrumInfo.size(); }
38

39
size_t SpectrumInfo::detectorCount() const { return m_spectrumInfo.detectorCount(); }
Owen Arnold's avatar
Owen Arnold committed
40

41
/// Returns a const reference to the SpectrumDefinition of the spectrum.
42
const SpectrumDefinition &SpectrumInfo::spectrumDefinition(const size_t index) const {
43
44
45
46
  m_experimentInfo.updateSpectrumDefinitionIfNecessary(index);
  return m_spectrumInfo.spectrumDefinition(index);
}

47
const Kernel::cow_ptr<std::vector<SpectrumDefinition>> &SpectrumInfo::sharedSpectrumDefinitions() const {
48
49
50
51
52
  for (size_t i = 0; i < size(); ++i)
    m_experimentInfo.updateSpectrumDefinitionIfNecessary(i);
  return m_spectrumInfo.sharedSpectrumDefinitions();
}

53
/// Returns true if the detector(s) associated with the spectrum are monitors.
54
bool SpectrumInfo::isMonitor(const size_t index) const {
55
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
56
57
58
    if (!m_detectorInfo.isMonitor(detIndex))
      return false;
  return true;
59
60
}

61
/// Returns true if the detector(s) associated with the spectrum are masked.
62
bool SpectrumInfo::isMasked(const size_t index) const {
63
  bool masked = true;
64
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
65
66
    masked &= m_detectorInfo.isMasked(detIndex);
  return masked;
67
68
69
70
71
72
73
74
}

/** Returns L2 (distance from sample to spectrum).
 *
 * For monitors this is defined such that L1+L2 = source-detector distance,
 * i.e., for a monitor in the beamline between source and sample L2 is negative.
 */
double SpectrumInfo::l2(const size_t index) const {
75
  double l2{0.0};
76
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
77
    l2 += m_detectorInfo.l2(detIndex);
78
  return l2 / static_cast<double>(spectrumDefinition(index).size());
79
80
}

David Fairbrother's avatar
David Fairbrother committed
81
82
/** Returns the scattering angle 2 theta in radians (angle w.r.t. to beam
 *direction).
83
 *
84
 * Throws an exception if the spectrum is a monitor.
85
 */
86
double SpectrumInfo::twoTheta(const size_t index) const {
87
  double twoTheta{0.0};
88
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
89
    twoTheta += m_detectorInfo.twoTheta(detIndex);
90
  return twoTheta / static_cast<double>(spectrumDefinition(index).size());
91
92
}

93
/** Returns the signed scattering angle 2 theta in radians (angle w.r.t. to beam
94
95
 * direction).
 *
96
 * Throws an exception if the spectrum is a monitor.
97
 */
98
double SpectrumInfo::signedTwoTheta(const size_t index) const {
99
  double signedTwoTheta{0.0};
100
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
101
    signedTwoTheta += m_detectorInfo.signedTwoTheta(detIndex);
102
  return signedTwoTheta / static_cast<double>(spectrumDefinition(index).size());
103
104
}

105
106
/** Returns the out-of-plane angle in radians (angle w.r.t. to
 * vecPointingHorizontal direction).
107
108
 *
 * Throws an exception if the spectrum is a monitor.
109
110
 */
double SpectrumInfo::azimuthal(const size_t index) const {
111
112
113
114
115
116
  double phi{0.0};
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
    phi += m_detectorInfo.azimuthal(detIndex);
  return phi / static_cast<double>(spectrumDefinition(index).size());
}

117
118
119
120
/** Calculate latitude and longitude for given spectrum index.
 *  @param index Index of the spectrum that lat/long are required for
 *  @return A pair containing the latitude and longitude values.
 */
121
std::pair<double, double> SpectrumInfo::geographicalAngles(const size_t index) const {
122
123
124
125
126
127
  double lat{0.0}, lon{0.0};
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index)) {
    auto latlong = m_detectorInfo.geographicalAngles(detIndex);
    lat += latlong.first;
    lon += latlong.second;
  }
128
129
  return std::pair<double, double>(lat / static_cast<double>(spectrumDefinition(index).size()),
                                   lon / static_cast<double>(spectrumDefinition(index).size()));
130
131
}

132
/// Returns the position of the spectrum with given index.
133
Kernel::V3D SpectrumInfo::position(const size_t index) const {
134
  Kernel::V3D newPos;
135
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
136
    newPos += m_detectorInfo.position(detIndex);
137
  return newPos / static_cast<double>(spectrumDefinition(index).size());
138
139
}

140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
/** Calculate average diffractometer constants (DIFA, DIFC, TZERO) of detectors
 * associated with this spectrum. Use calibrated values where possible, filling
 * in with uncalibrated values where they're missing
 *  @param index Index of the spectrum that constants are required for
 *  @param warningDets A vector containing the det ids where an uncalibrated
 * value was used in the situation where some dets have calibrated values and
 * some don't
 *  @return tuple containing the average constants
 */
std::tuple<double, double, double>
SpectrumInfo::diffractometerConstants(const size_t index,
                                      std::vector<detid_t> &warningDets) const {
  if (m_detectorInfo.isScanning()) {
    throw std::runtime_error("Retrieval of diffractometer constants not "
                             "implemented for scanning instrument");
  }
  auto spectrumDef = checkAndGetSpectrumDefinition(index);
  std::vector<size_t> detectorIndicesOnly;
  std::vector<detid_t> calibratedDets;
  std::vector<detid_t> uncalibratedDets;
  std::transform(spectrumDef.begin(), spectrumDef.end(),
                 std::back_inserter(detectorIndicesOnly),
                 [](auto const &pair) { return pair.first; });
  double difa{0.}, difc{0.}, tzero{0.};
  for (const auto &detIndex : detectorIndicesOnly) {
    auto newDiffConstants = m_detectorInfo.diffractometerConstants(
        detIndex, calibratedDets, uncalibratedDets);
    difa += std::get<0>(newDiffConstants);
    difc += std::get<1>(newDiffConstants);
    tzero += std::get<2>(newDiffConstants);
  }

  if (calibratedDets.size() > 0 && uncalibratedDets.size() > 0) {
    warningDets.insert(warningDets.end(), uncalibratedDets.begin(),
                       uncalibratedDets.end());
  };
Danny Hindson's avatar
Danny Hindson committed
176
  // if no calibration is found then return difc only based on the average
177
  // of the detector L2 and twoThetas.
Danny Hindson's avatar
Danny Hindson committed
178
  if (calibratedDets.size() == 0) {
179
    return {0., difcUncalibrated(index), 0.};
Danny Hindson's avatar
Danny Hindson committed
180
  }
181
182
183
184
185
  return {difa / static_cast<double>(spectrumDefinition(index).size()),
          difc / static_cast<double>(spectrumDefinition(index).size()),
          tzero / static_cast<double>(spectrumDefinition(index).size())};
}

186
187
188
189
190
191
192
193
194
195
196
197
/** Calculate average diffractometer constants (DIFA, DIFC, TZERO) of detectors
 * associated with this spectrum. Use calibrated values where possible, filling
 * in with uncalibrated values where they're missing
 *  @param index Index of the spectrum that constants are required for
 *  @return tuple containing the average constants
 */
std::tuple<double, double, double>
SpectrumInfo::diffractometerConstants(const size_t index) const {
  std::vector<int> warningDets;
  return diffractometerConstants(index, warningDets);
}

198
199
200
201
202
203
/** Calculate average uncalibrated DIFC value of detectors associated with this
 * spectrum
 *  @param index Index of the spectrum that DIFC is required for
 *  @return The average DIFC
 */
double SpectrumInfo::difcUncalibrated(const size_t index) const {
204
205
206
207
208
209
  // calculate difc based on the average of the detector L2 and twoThetas. This
  // will be different to the average of the per detector difcs. This is for
  // backwards compatibility because Mantid always used to calculate spectrum
  // level difc's this way
  return 1. / Mantid::Geometry::Conversion::tofToDSpacingFactor(
                  l1(), l2(index), twoTheta(index), 0.);
210
211
}

Danny Hindson's avatar
Danny Hindson committed
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
/** Get the detector values relevant to unit conversion for a workspace index
 * @param inputUnit :: The input unit (Empty implies "all")
 * @param outputUnit :: The output unit (Empty implies "all")
 * @param emode :: The energy mode
 * @param signedTheta :: Return twotheta with sign or without
 * @param wsIndex :: The workspace index
 * @param pmap :: a map containing values for conversion parameters that are
required by unit classes to perform their conversions eg efixed. It can contain
values on the way in if a look up isn't desired here eg if value supplied in
parameters to the calling algorithm
 */
void SpectrumInfo::getDetectorValues(const Kernel::Unit &inputUnit,
                                     const Kernel::Unit &outputUnit,
                                     const Kernel::DeltaEMode::Type emode,
                                     const bool signedTheta, int64_t wsIndex,
                                     UnitParametersMap &pmap) const {
  if (!hasDetectors(wsIndex))
    return;
  pmap[UnitParams::l2] = l2(wsIndex);

  if (!isMonitor(wsIndex)) {
    // The scattering angle for this detector (in radians).
    try {
      if (signedTheta)
        pmap[UnitParams::twoTheta] = signedTwoTheta(wsIndex);
      else
        pmap[UnitParams::twoTheta] = twoTheta(wsIndex);
    } catch (const std::runtime_error &e) {
      g_log.warning(e.what());
      pmap[UnitParams::twoTheta] = std::numeric_limits<double>::quiet_NaN();
    }
    if (emode != Kernel::DeltaEMode::Elastic &&
        pmap.find(UnitParams::efixed) == pmap.end()) {
      std::shared_ptr<const Geometry::IDetector> det(&detector(wsIndex),
                                                     Mantid::NoDeleting());
      try {
        pmap[UnitParams::efixed] =
            m_experimentInfo.getEFixedGivenEMode(det, emode);
        g_log.debug() << "Detector: " << det->getID()
                      << " EFixed: " << pmap[UnitParams::efixed] << "\n";
      } catch (std::runtime_error) {
        // let the unit classes work out if this is a problem
      }
    }

    std::vector<detid_t> warnDetIds;
    try {
259
260
      std::vector<std::string> diffConstUnits = {"dSpacing", "MomentumTransfer",
                                                 "Empty"};
Danny Hindson's avatar
Danny Hindson committed
261
      if ((emode == Kernel::DeltaEMode::Elastic) &&
262
263
264
265
          (std::find(diffConstUnits.begin(), diffConstUnits.end(),
                     inputUnit.unitID()) != diffConstUnits.end()) &&
          (std::find(diffConstUnits.begin(), diffConstUnits.end(),
                     outputUnit.unitID()) != diffConstUnits.end())) {
Danny Hindson's avatar
Danny Hindson committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
        auto [difa, difc, tzero] = diffractometerConstants(wsIndex, warnDetIds);
        pmap[UnitParams::difa] = difa;
        pmap[UnitParams::difc] = difc;
        pmap[UnitParams::tzero] = tzero;
        if (warnDetIds.size() > 0) {
          createDetectorIdLogMessages(warnDetIds, wsIndex);
        }
      } else {
        pmap[UnitParams::difc] = difcUncalibrated(wsIndex);
      }
    } catch (const std::runtime_error &e) {
      g_log.warning(e.what());
    }
  } else {
    pmap[UnitParams::twoTheta] = 0.0;
    pmap[UnitParams::efixed] = DBL_MIN;
    // Energy transfer is meaningless for a monitor, so set l2 to 0.
    if (outputUnit.unitID().find("DeltaE") != std::string::npos) {
      pmap[UnitParams::l2] = 0.0;
    }
    pmap[UnitParams::difc] = 0;
  }
}

void SpectrumInfo::createDetectorIdLogMessages(
    const std::vector<detid_t> &detids, int64_t wsIndex) const {
  std::string detIDstring;
  auto iter = detids.begin();
  auto itEnd = detids.end();
  for (; iter != itEnd; ++iter) {
    detIDstring += std::to_string(*iter) + ",";
  }

  if (!detIDstring.empty()) {
    detIDstring.pop_back(); // Drop last comma
  }
  g_log.warning(
      "Incomplete set of calibrated diffractometer constants found for "
      "workspace index" +
      std::to_string(wsIndex) + ". Using uncalibrated values for detectors " +
      detIDstring);
}

309
/// Returns true if the spectrum is associated with detectors in the instrument.
310
bool SpectrumInfo::hasDetectors(const size_t index) const {
311
312
  // Workspaces can contain invalid detector IDs. Those IDs will be silently
  // ignored here until this is fixed.
313
  return spectrumDefinition(index).size() > 0;
314
315
}

316
/// Returns true if the spectrum is associated with exactly one detector.
317
bool SpectrumInfo::hasUniqueDetector(const size_t index) const {
318
319
  // Workspaces can contain invalid detector IDs. Those IDs will be silently
  // ignored here until this is fixed.
320
  return spectrumDefinition(index).size() == 1;
321
322
}

323
/** Set the mask flag of the spectrum with given index. Not thread safe.
324
325
326
 *
 * Currently this simply sets the mask flags for the underlying detectors. */
void SpectrumInfo::setMasked(const size_t index, bool masked) {
327
  for (const auto &detIndex : checkAndGetSpectrumDefinition(index))
328
    m_detectorInfo.setMasked(detIndex, masked);
329
330
}

331
332
/// Return a const reference to the detector or detector group of the spectrum
/// with given index.
333
const Geometry::IDetector &SpectrumInfo::detector(const size_t index) const { return getDetector(index); }
334

335
/// Returns the source position.
336
Kernel::V3D SpectrumInfo::sourcePosition() const { return m_detectorInfo.sourcePosition(); }
337
338

/// Returns the sample position.
339
Kernel::V3D SpectrumInfo::samplePosition() const { return m_detectorInfo.samplePosition(); }
340

341
/// Returns L1 (distance from source to sample).
342
double SpectrumInfo::l1() const { return m_detectorInfo.l1(); }
343

344
const Geometry::IDetector &SpectrumInfo::getDetector(const size_t index) const {
345
  auto thread = static_cast<size_t>(PARALLEL_THREAD_NUMBER);
346
  if (m_lastIndex[thread] == index)
347
    return *m_lastDetector[thread];
348

349
350
351
  // Note: This function body has big overlap with the method
  // MatrixWorkspace::getDetector(). The plan is to eventually remove the
  // latter, once SpectrumInfo is in widespread use.
352
  const auto &specDef = spectrumDefinition(index);
353
  const size_t ndets = specDef.size();
354
355
  if (ndets == 1) {
    // If only 1 detector for the spectrum number, just return it
356
    const auto detIndex = specDef[0].first;
357
    m_lastDetector[thread] = m_detectorInfo.getDetectorPtr(detIndex);
358
359
360
361
  } else if (ndets == 0) {
    throw Kernel::Exception::NotFoundError("MatrixWorkspace::getDetector(): No "
                                           "detectors for this workspace "
                                           "index.",
362
                                           std::to_string(index));
363
364
  } else {
    // Else need to construct a DetectorGroup and use that
365
    std::vector<std::shared_ptr<const Geometry::IDetector>> det_ptrs;
366
367
    for (const auto &index : specDef) {
      const auto detIndex = index.first;
368
      det_ptrs.emplace_back(m_detectorInfo.getDetectorPtr(detIndex));
369
    }
370
    m_lastDetector[thread] = std::make_shared<Geometry::DetectorGroup>(det_ptrs);
371
  }
372
  m_lastIndex[thread] = index;
373
  return *m_lastDetector[thread];
374
375
}

376
const SpectrumDefinition &SpectrumInfo::checkAndGetSpectrumDefinition(const size_t index) const {
377
  if (spectrumDefinition(index).size() == 0)
378
379
    throw Kernel::Exception::NotFoundError("SpectrumInfo: No detectors for this workspace index.",
                                           std::to_string(index));
380
  return spectrumDefinition(index);
381
382
}

Bhuvan Bezawada's avatar
Bhuvan Bezawada committed
383
// Begin method for iterator
384
385
386
387
388
SpectrumInfoIt SpectrumInfo::begin() { return SpectrumInfoIt(*this, 0); }

// End method for iterator
SpectrumInfoIt SpectrumInfo::end() { return SpectrumInfoIt(*this, size()); }

389
const SpectrumInfoConstIt SpectrumInfo::cbegin() const { return SpectrumInfoConstIt(*this, 0); }
Bhuvan Bezawada's avatar
Bhuvan Bezawada committed
390
391

// End method for iterator
392
const SpectrumInfoConstIt SpectrumInfo::cend() const { return SpectrumInfoConstIt(*this, size()); }
Bhuvan Bezawada's avatar
Bhuvan Bezawada committed
393

394
395
} // namespace API
} // namespace Mantid