VesuvioCalculateGammaBackground.cpp 24.5 KB
Newer Older
1
2
3
4
5
6
// 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
// SPDX - License - Identifier: GPL - 3.0 +
7
#include "MantidCurveFitting/Algorithms/VesuvioCalculateGammaBackground.h"
8
#include "MantidCurveFitting/Algorithms/ConvertToYSpace.h"
9
10
#include "MantidCurveFitting/Functions/ComptonProfile.h"
#include "MantidCurveFitting/Functions/VesuvioResolution.h"
11

12
13
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/FunctionFactory.h"
14
#include "MantidAPI/FunctionProperty.h"
15
#include "MantidAPI/HistogramValidator.h"
16
#include "MantidAPI/WorkspaceFactory.h"
17
#include "MantidAPI/WorkspaceUnitValidator.h"
18

19
#include "MantidKernel/ArrayProperty.h"
20
#include "MantidKernel/CompositeValidator.h"
21

22
#include "MantidGeometry/Instrument.h"
23
#include "MantidGeometry/Instrument/ReferenceFrame.h"
Simon Heybrock's avatar
Simon Heybrock committed
24
#include "MantidGeometry/Objects/BoundingBox.h"
25
#include "MantidGeometry/Objects/CSGObject.h"
26

27
28
namespace Mantid {
namespace CurveFitting {
29
namespace Algorithms {
30
31
using namespace API;
using namespace Kernel;
32
33
using namespace CurveFitting;
using namespace CurveFitting::Functions;
34
using namespace std;
35
using std::placeholders::_1;
36
37

// Subscription
38
DECLARE_ALGORITHM(VesuvioCalculateGammaBackground)
39
40
41
42
43
44
45
46
47
48
49
50

namespace {
/// Number of elements in theta direction integration
size_t NTHETA = 5;
/// Number of elements in up direction integration
size_t NUP = 5;

/// Degrees to radians
double DEG2RAD = M_PI / 180.0;
/// Wavelength of absorption (30603e-24 * 6e19). Constant came from VMS
double ABSORB_WAVELENGTH = 1.83618;
/// Start of forward scattering spectrum numbers (inclusive)
51
specnum_t FORWARD_SCATTER_SPECMIN = 135;
52
/// End of forward scattering spectrum numbers (inclusive)
53
specnum_t FORWARD_SCATTER_SPECMAX = 198;
LamarMoore's avatar
LamarMoore committed
54
} // namespace
55
56
57
58
59
60

//--------------------------------------------------------------------------------------------------------
// Public members
//--------------------------------------------------------------------------------------------------------

/// Default constructor
61
VesuvioCalculateGammaBackground::VesuvioCalculateGammaBackground()
62
63
64
    : Algorithm(), m_inputWS(), m_indices(), m_profileFunction(), m_npeaks(0),
      m_reversed(), m_samplePos(), m_l1(0.0), m_foilRadius(0.0),
      m_foilUpMin(0.0), m_foilUpMax(0.0), m_foils0(), m_foils1(),
65
      m_backgroundWS(), m_correctedWS() {}
66
67

/// Destructor
68
VesuvioCalculateGammaBackground::~VesuvioCalculateGammaBackground() {
69
70
71
72
73
74
75
  m_indices.clear();
}

//--------------------------------------------------------------------------------------------------------
// Private members
//--------------------------------------------------------------------------------------------------------

76
77
const std::string VesuvioCalculateGammaBackground::name() const {
  return "VesuvioCalculateGammaBackground";
78
79
}

80
int VesuvioCalculateGammaBackground::version() const { return 1; }
81

82
const std::string VesuvioCalculateGammaBackground::category() const {
83
  return "CorrectionFunctions\\BackgroundCorrections";
84
85
}

86
void VesuvioCalculateGammaBackground::init() {
87
88
89
90

  auto wsValidator = boost::make_shared<CompositeValidator>();
  wsValidator->add<WorkspaceUnitValidator>("TOF");
  wsValidator->add<HistogramValidator>(false); // point data
91
  declareProperty(std::make_unique<WorkspaceProperty<>>(
92
                      "InputWorkspace", "", Direction::Input, wsValidator),
93
94
95
                  "An input workspace containing TOF data");

  declareProperty(
Sam Jenkins's avatar
Sam Jenkins committed
96
97
      std::make_unique<API::FunctionProperty>("ComptonFunction",
                                              Direction::InOut),
98
99
100
      "Function that is able to compute the mass spectrum for the input data"
      "This will usually be the output from the Fitting");

101
  declareProperty(std::make_unique<ArrayProperty<int>>("WorkspaceIndexList"),
102
103
104
105
                  "Indices of the spectra to include in the correction. If "
                  "provided, the output only include these spectra\n"
                  "(Default: all spectra from input)");

Sam Jenkins's avatar
Sam Jenkins committed
106
107
  declareProperty(std::make_unique<WorkspaceProperty<>>("BackgroundWorkspace",
                                                        "", Direction::Output),
108
                  "A new workspace containing the calculated background.");
109
  declareProperty(
110
      std::make_unique<WorkspaceProperty<>>("CorrectedWorkspace", "",
Sam Jenkins's avatar
Sam Jenkins committed
111
                                            Direction::Output),
112
113
114
115
      "A new workspace containing the calculated background subtracted from "
      "the input.");
}

116
void VesuvioCalculateGammaBackground::exec() {
117
118
119
  retrieveInputs();
  createOutputWorkspaces();

120
  const auto nhist = static_cast<int64_t>(m_indices.size());
121
122
  const int64_t nreports =
      10 + nhist * (m_npeaks + 2 * m_foils0.size() * NTHETA * NUP * m_npeaks);
123
  m_progress = std::make_unique<Progress>(this, 0.0, 1.0, nreports);
124

125
126
  PARALLEL_FOR_IF(
      Kernel::threadSafe(*m_inputWS, *m_correctedWS, *m_backgroundWS))
127
128
129
130
131
132
133
134
135
  for (int64_t i = 0; i < nhist; ++i) {
    PARALLEL_START_INTERUPT_REGION
    const size_t outputIndex = i;
    auto indexIter = m_indices.cbegin();
    std::advance(indexIter, i);
    const size_t inputIndex = indexIter->second;

    if (!calculateBackground(inputIndex, outputIndex)) {
      g_log.information("No detector defined for index=" +
136
                        std::to_string(inputIndex) + ". Skipping correction.");
137
138
    }

139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  setProperty("BackgroundWorkspace", m_backgroundWS);
  setProperty("CorrectedWorkspace", m_correctedWS);
}

/**
 * Calculate the background from the input spectrum and assign the value to the
 * output one
 * @param inputIndex The index on the input workspace on which to operate
 * @param outputIndex The index on the output workspace where the results are
 * stored
 * @return True if the background was subtracted, false otherwise
 */
Elliot Oram's avatar
Elliot Oram committed
155
156
bool VesuvioCalculateGammaBackground::calculateBackground(
    const size_t inputIndex, const size_t outputIndex) {
157
  // Copy X values
158
159
  m_backgroundWS->setSharedX(outputIndex, m_inputWS->sharedX(inputIndex));
  m_correctedWS->setSharedX(outputIndex, m_inputWS->sharedX(inputIndex));
160
  // Copy errors to corrected
161
  m_correctedWS->setSharedE(outputIndex, m_inputWS->sharedE(inputIndex));
162
163

  try {
164
165
166
167
    const auto &inSpec = m_inputWS->getSpectrum(inputIndex);
    const specnum_t spectrumNo(inSpec.getSpectrumNo());
    m_backgroundWS->getSpectrum(outputIndex).copyInfoFrom(inSpec);
    m_correctedWS->getSpectrum(outputIndex).copyInfoFrom(inSpec);
168
169
170
171
172

    if (spectrumNo >= FORWARD_SCATTER_SPECMIN &&
        spectrumNo <= FORWARD_SCATTER_SPECMAX) {
      applyCorrection(inputIndex, outputIndex);
    } else {
173
      g_log.information("Spectrum " + std::to_string(spectrumNo) +
174
175
                        " not in forward scatter range. Skipping correction.");
      // Leave background at 0 and just copy data to corrected
176
      m_correctedWS->setSharedY(outputIndex, m_inputWS->sharedY(inputIndex));
177
    }
178
179
180
181
182
183
184
    return true;
  } catch (Exception::NotFoundError &) {
    return false;
  }
}

/**
LamarMoore's avatar
LamarMoore committed
185
186
187
188
189
190
191
 * Calculate & apply gamma correction for the given index of the
 * input workspace
 * @param inputIndex A workspace index that defines the input spectrum to
 * correct
 * @param outputIndex A workspace index that defines the output to hold the
 * results
 */
Elliot Oram's avatar
Elliot Oram committed
192
193
void VesuvioCalculateGammaBackground::applyCorrection(
    const size_t inputIndex, const size_t outputIndex) {
194
195
196
197
198
199
200
201
202
203
204
205
206
  m_progress->report("Computing TOF from detector");

  // results go straight in m_correctedWS to save memory allocations
  calculateSpectrumFromDetector(inputIndex, outputIndex);

  m_progress->report("Computing TOF foils");
  // Output goes to m_background to save memory allocations
  calculateBackgroundFromFoils(inputIndex, outputIndex);

  m_progress->report("Computing correction to input");
  // Compute total counts from input data, (detector-foil) contributions
  // assume constant binning
  const size_t nbins = m_correctedWS->blocksize();
207
208
209
210
211
  const auto &inY = m_inputWS->y(inputIndex);
  auto &detY = m_correctedWS->mutableY(outputIndex);
  auto &foilY = m_backgroundWS->mutableY(outputIndex);
  const double deltaT =
      m_correctedWS->x(outputIndex)[1] - m_correctedWS->x(outputIndex)[0];
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234

  double dataCounts(0.0), simulCounts(0.0);
  for (size_t j = 0; j < nbins; ++j) {
    dataCounts += inY[j] * deltaT;
    simulCounts += (detY[j] - foilY[j]) * deltaT;
  }

  // Now corrected for the calculated background
  const double corrFactor = dataCounts / simulCounts;
  if (g_log.is(Logger::Priority::PRIO_INFORMATION))
    g_log.information() << "Correction factor for background=" << corrFactor
                        << "\n";

  for (size_t j = 0; j < nbins; ++j) {
    // m_backgroundWS already contains the foil values, careful not to overwrite
    // them
    double &foilValue = foilY[j]; // non-const reference
    foilValue *= corrFactor;
    detY[j] = (inY[j] - foilValue);
  }
}

/**
LamarMoore's avatar
LamarMoore committed
235
236
237
238
239
 * Results are placed in the mapped index on the output corrected workspace
 * @param inputIndex Workspace index that defines the input spectrum to correct
 * @param outputIndex Workspace index that defines the spectrum to hold the
 * results
 */
240
void VesuvioCalculateGammaBackground::calculateSpectrumFromDetector(
241
242
243
244
    const size_t inputIndex, const size_t outputIndex) {
  // -- Setup detector & resolution parameters --
  DetectorParams detPar =
      ConvertToYSpace::getDetectorParameters(m_inputWS, inputIndex);
245
  CurveFitting::Functions::ResolutionParams detRes =
Nick Draper's avatar
Nick Draper committed
246
247
      CurveFitting::Functions::VesuvioResolution::getResolutionParameters(
          m_inputWS, inputIndex);
248
249
250
251

  // Compute a time of flight spectrum convolved with a Voigt resolution
  // function for each mass
  // at the detector point & sum to a single spectrum
252
  auto &ctdet = m_correctedWS->mutableY(outputIndex);
253
  std::vector<double> tmpWork(ctdet.size());
254
255
  ctdet = calculateTofSpectrum(ctdet.rawData(), tmpWork, outputIndex, detPar,
                               detRes);
256
257
258
  // Correct for distance to the detector: 0.5/l2^2
  const double detDistCorr = 0.5 / detPar.l2 / detPar.l2;
  std::transform(ctdet.begin(), ctdet.end(), ctdet.begin(),
259
                 std::bind(std::multiplies<double>(), _1, detDistCorr));
260
261
262
}

/**
LamarMoore's avatar
LamarMoore committed
263
264
265
266
267
268
 * Calculate & apply gamma correction for the given index of the
 * input workspace
 * @param inputIndex Workspace index that defines the input spectrum to correct
 * @param outputIndex Workspace index that defines the spectrum to hold the
 * results
 */
269
void VesuvioCalculateGammaBackground::calculateBackgroundFromFoils(
270
271
272
273
    const size_t inputIndex, const size_t outputIndex) {
  // -- Setup detector & resolution parameters --
  DetectorParams detPar =
      ConvertToYSpace::getDetectorParameters(m_inputWS, inputIndex);
274
  CurveFitting::Functions::ResolutionParams detRes =
Nick Draper's avatar
Nick Draper committed
275
276
      CurveFitting::Functions::VesuvioResolution::getResolutionParameters(
          m_inputWS, inputIndex);
277
278
279

  const size_t nxvalues = m_backgroundWS->blocksize();
  std::vector<double> foilSpectrum(nxvalues);
280
  auto &ctfoil = m_backgroundWS->mutableY(outputIndex);
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298

  // Compute (C1 - C0) where C1 is counts in pos 1 and C0 counts in pos 0
  assert(m_foils0.size() == m_foils1.size());
  for (size_t i = 0; i < m_foils0.size(); ++i) {
    foilSpectrum.assign(nxvalues, 0.0);
    calculateBackgroundSingleFoil(foilSpectrum, outputIndex, m_foils1[i],
                                  detPar.pos, detPar, detRes);
    // sum spectrum values from first position
    std::transform(ctfoil.begin(), ctfoil.end(), foilSpectrum.begin(),
                   ctfoil.begin(), std::plus<double>());

    foilSpectrum.assign(nxvalues, 0.0);
    calculateBackgroundSingleFoil(foilSpectrum, outputIndex, m_foils0[i],
                                  detPar.pos, detPar, detRes);
    // subtract spectrum values from zeroth position
    std::transform(ctfoil.begin(), ctfoil.end(), foilSpectrum.begin(),
                   ctfoil.begin(), std::minus<double>());
  }
LamarMoore's avatar
LamarMoore committed
299
300
301
  bool reversed =
      (m_reversed.count(m_inputWS->getSpectrum(inputIndex).getSpectrumNo()) !=
       0);
302
303
304
305
  // This is quicker than the if within the loop
  if (reversed) {
    // The reversed ones should be (C0 - C1)
    std::transform(ctfoil.begin(), ctfoil.end(), ctfoil.begin(),
306
                   std::bind(std::multiplies<double>(), _1, -1.0));
307
308
309
310
  }
}

/**
LamarMoore's avatar
LamarMoore committed
311
312
313
314
315
316
317
318
319
320
321
322
 * Integrates over the foil area defined by the foil radius to accumulate an
 * estimate of the counts
 * resulting from this region
 * @param ctfoil Output vector to hold results
 * @param wsIndex Index on output background workspaces currently operating
 * @param foilInfo Foil description object
 * @param detPos The pre-calculated detector V3D
 * @param detPar DetectorParams object that defines information on the detector
 * associated with spectrum at wsIndex
 * @param detRes ResolutionParams object that defines information on the
 * resolution associated with spectrum at wsIndex
 */
323
void VesuvioCalculateGammaBackground::calculateBackgroundSingleFoil(
324
325
326
327
    std::vector<double> &ctfoil, const size_t wsIndex, const FoilInfo &foilInfo,
    const V3D &detPos, const DetectorParams &detPar,
    const ResolutionParams &detRes) {
  /** Integrates over the foils
LamarMoore's avatar
LamarMoore committed
328
329
330
   *  by dividing into 2cm^2 elements
   *  The integration is performed in cylindrical coordinates
   */
331
332
333
334
335
336
337
338
339
340
341

  const double thetaStep =
      (foilInfo.thetaMax - foilInfo.thetaMin) / static_cast<double>(NTHETA);
  const double thetaStepRad = thetaStep * DEG2RAD;
  const double upStep = (m_foilUpMax - m_foilUpMin) / static_cast<double>(NUP);
  const double elementArea = abs(m_foilRadius * upStep * thetaStepRad);
  V3D elementPos; // reusable vector for computing distances

  // Structs to hold geometry & resolution information
  DetectorParams foilPar = detPar; // copy
  foilPar.t0 = 0.0;
342
  CurveFitting::Functions::ResolutionParams foilRes = detRes; // copy
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
  foilRes.dEnGauss = foilInfo.gaussWidth;
  foilRes.dEnLorentz = foilInfo.lorentzWidth;

  const size_t nvalues = ctfoil.size();
  std::vector<double> singleElement(nvalues), tmpWork(nvalues);

  for (size_t i = 0; i < NTHETA; ++i) {
    double thetaZeroRad =
        (foilInfo.thetaMin + (static_cast<double>(i) + 0.5) * thetaStep) *
        DEG2RAD;
    elementPos.setZ(m_foilRadius * cos(thetaZeroRad));
    elementPos.setX(m_foilRadius * sin(thetaZeroRad));
    for (size_t j = 0; j < NUP; ++j) {
      double ypos = m_foilUpMin + (static_cast<double>(j) + 0.5) * upStep;
      elementPos.setY(ypos);
      foilPar.l2 = m_samplePos.distance(elementPos);
      foilPar.theta = acos(m_foilRadius * cos(thetaZeroRad) /
                           foilPar.l2); // scattering angle in radians

      // Spectrum for this position
      singleElement.assign(nvalues, 0.0);
364
365
      singleElement = calculateTofSpectrum(singleElement, tmpWork, wsIndex,
                                           foilPar, foilRes);
366
367
368
369
370
371
372
373
374
375
376
377
378

      // Correct for absorption & distance
      const double den = (elementPos.Z() * cos(thetaZeroRad) +
                          elementPos.X() * sin(thetaZeroRad));
      const double absorbFactor =
          1.0 / (1.0 - exp(-ABSORB_WAVELENGTH * foilPar.l2 / den));
      const double elemDetDist = elementPos.distance(detPos);
      const double distFactor =
          elementArea /
          (4.0 * M_PI * foilPar.l2 * foilPar.l2 * elemDetDist * elemDetDist);
      // Add on to other contributions
      for (size_t k = 0; k < nvalues; ++k) {
        ctfoil[k] += singleElement[k] * absorbFactor * distFactor;
379
380
      }
    }
381
382
383
384
  }
}

/**
LamarMoore's avatar
LamarMoore committed
385
386
387
388
389
390
391
392
 * Uses the compton profile functions to compute a particular mass spectrum
 * @param inSpectrum The value of the computed spectrum
 * @param tmpWork Pre-allocated working area that will be overwritten
 * @param wsIndex Index on the output background workspace that gives the X
 * values to use
 * @param detpar Struct containing parameters about the detector
 * @param respar Struct containing parameters about the resolution
 */
393
std::vector<double> VesuvioCalculateGammaBackground::calculateTofSpectrum(
394
    const std::vector<double> &inSpectrum, std::vector<double> &tmpWork,
395
396
    const size_t wsIndex, const DetectorParams &detpar,
    const ResolutionParams &respar) {
397
  assert(inSpectrum.size() == tmpWork.size());
398
399

  // Assumes the input is in seconds, transform it temporarily
400
  auto &tseconds = m_backgroundWS->mutableX(wsIndex);
401
  std::transform(tseconds.begin(), tseconds.end(), tseconds.begin(),
402
                 std::bind(std::multiplies<double>(), _1, 1e-6));
403
404
405
406
407
408
409

  // retrieveInputs ensures we will get a composite function and that each
  // member is a ComptonProfile
  // we can't static_cast though due to the virtual inheritance with IFunction
  auto profileFunction = boost::dynamic_pointer_cast<CompositeFunction>(
      FunctionFactory::Instance().createInitialized(m_profileFunction));

410
  std::vector<double> correctedVals(inSpectrum);
411

412
  for (size_t i = 0; i < m_npeaks; ++i) {
Nick Draper's avatar
Nick Draper committed
413
414
415
    auto profile =
        boost::dynamic_pointer_cast<CurveFitting::Functions::ComptonProfile>(
            profileFunction->getFunction(i));
416
417
    profile->disableLogging();
    profile->setUpForFit();
418
419

    // Fix the Mass parameter
420
    profile->fix(0);
421

422
    profile->cacheYSpaceValues(m_backgroundWS->points(wsIndex), detpar, respar);
423
424
425

    profile->massProfile(tmpWork.data(), tmpWork.size());
    // Add to final result
426
427
428

    std::transform(correctedVals.begin(), correctedVals.end(), tmpWork.begin(),
                   correctedVals.begin(), std::plus<double>());
429
430
431
432
    m_progress->report();
  }
  // Put X back microseconds
  std::transform(tseconds.begin(), tseconds.end(), tseconds.begin(),
433
                 std::bind(std::multiplies<double>(), _1, 1e6));
434
  return correctedVals;
435
436
437
}

/**
LamarMoore's avatar
LamarMoore committed
438
439
 * Caches input details for the peak information
 */
440
void VesuvioCalculateGammaBackground::retrieveInputs() {
441
442
  m_inputWS = getProperty("InputWorkspace");
  m_profileFunction = getPropertyValue("ComptonFunction");
443
  if (m_profileFunction.find(';') == std::string::npos) // not composite
444
445
446
447
448
449
450
451
452
453
  {
    m_profileFunction = "composite=CompositeFunction;" + m_profileFunction;
  }

  IFunction_sptr profileFunction =
      FunctionFactory::Instance().createInitialized(m_profileFunction);
  if (auto composite =
          boost::dynamic_pointer_cast<CompositeFunction>(profileFunction)) {
    m_npeaks = composite->nFunctions();
    for (size_t i = 0; i < m_npeaks; ++i) {
Nick Draper's avatar
Nick Draper committed
454
455
456
      auto single =
          boost::dynamic_pointer_cast<CurveFitting::Functions::ComptonProfile>(
              composite->getFunction(i));
457
458
459
      if (!single) {
        throw std::invalid_argument("Invalid function. Composite must contain "
                                    "only ComptonProfile functions");
460
      }
461
    }
462
463
464
465
466
467
468
469
  } else {
    throw std::invalid_argument(
        "Invalid function found. Expected ComptonFunction to contain a "
        "composite of ComptonProfiles or a single ComptonProfile.");
  }

  // Spectrum numbers whose calculation of background from foils is reversed
  m_reversed.clear();
470
  for (specnum_t i = 143; i < 199; ++i) {
471
472
473
474
475
476
477
478
479
480
    if ((i >= 143 && i <= 150) || (i >= 159 && i <= 166) ||
        (i >= 175 && i <= 182) || (i >= 191 && i <= 198))
      m_reversed.insert(i);
  }

  // Workspace indices mapping input->output
  m_indices.clear();
  std::vector<int> requestedIndices = getProperty("WorkspaceIndexList");
  if (requestedIndices.empty()) {
    for (size_t i = 0; i < m_inputWS->getNumberHistograms(); ++i) {
481
      m_indices.emplace(i, i); // 1-to-1
482
    }
483
484
  } else {
    for (size_t i = 0; i < requestedIndices.size(); ++i) {
485
      m_indices.emplace(
486
          i, static_cast<size_t>(
487
                 requestedIndices[i])); // user-requested->increasing on output
488
    }
489
490
491
492
493
494
  }

  cacheInstrumentGeometry();
}

/**
LamarMoore's avatar
LamarMoore committed
495
496
 * Create & cache output workspaces
 */
497
void VesuvioCalculateGammaBackground::createOutputWorkspaces() {
498
499
500
501
502
503
  const size_t nhist = m_indices.size();
  m_backgroundWS = WorkspaceFactory::Instance().create(m_inputWS, nhist);
  m_correctedWS = WorkspaceFactory::Instance().create(m_inputWS, nhist);
}

/**
LamarMoore's avatar
LamarMoore committed
504
 */
505
void VesuvioCalculateGammaBackground::cacheInstrumentGeometry() {
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
  auto inst = m_inputWS->getInstrument();
  auto refFrame = inst->getReferenceFrame();
  auto upAxis = refFrame->pointingUp();
  auto source = inst->getSource();
  auto sample = inst->getSample();
  m_samplePos = sample->getPos();
  m_l1 = m_samplePos.distance(source->getPos());

  // foils
  auto changer = boost::dynamic_pointer_cast<const Geometry::IObjComponent>(
      inst->getComponentByName("foil-changer"));
  if (!changer) {
    throw std::invalid_argument(
        "Input workspace has no component named foil-changer. "
        "One is required to define integration area.");
  }

  // 'height' of box sets limits in beam direction
  Geometry::BoundingBox boundBox;
  changer->getBoundingBox(boundBox);
  m_foilUpMin = boundBox.minPoint()[upAxis];
  m_foilUpMax = boundBox.maxPoint()[upAxis];

  // foil geometry
  // there should be the same number in each position
  const auto &pmap = m_inputWS->constInstrumentParameters();
  auto foils0 = inst->getAllComponentsWithName("foil-pos0");
  auto foils1 = inst->getAllComponentsWithName("foil-pos1");
  const size_t nfoils = foils0.size();
  if (nfoils != foils1.size()) {
    std::ostringstream os;
    os << "Mismatch in number of foils between pos 0 & 1: pos0=" << nfoils
       << ", pos1=" << foils1.size();
    throw std::runtime_error(os.str());
  }
  // It is assumed that the foils all lie on a circle of the same radius from
  // the sample position
  auto firstFoilPos = foils0[0]->getPos();
  double dummy(0.0);
  firstFoilPos.getSpherical(m_foilRadius, dummy, dummy);

  // Cache min/max theta values
  m_foils0.resize(nfoils);
  m_foils1.resize(nfoils);
  for (size_t i = 0; i < nfoils; ++i) {
    const auto &foil0 = foils0[i];
    auto thetaRng0 = calculateThetaRange(foil0, m_foilRadius,
                                         refFrame->pointingHorizontal());
    FoilInfo descr;
    descr.thetaMin = thetaRng0.first;
    descr.thetaMax = thetaRng0.second;
    descr.lorentzWidth =
558
        ConvertToYSpace::getComponentParameter(*foil0, pmap, "hwhm_lorentz");
559
    descr.gaussWidth =
560
        ConvertToYSpace::getComponentParameter(*foil0, pmap, "sigma_gauss");
561
562
563
564
565
566
567
568
    m_foils0[i] = descr; // copy

    const auto &foil1 = foils1[i];
    auto thetaRng1 = calculateThetaRange(foil1, m_foilRadius,
                                         refFrame->pointingHorizontal());
    descr.thetaMin = thetaRng1.first;
    descr.thetaMax = thetaRng1.second;
    descr.lorentzWidth =
569
        ConvertToYSpace::getComponentParameter(*foil1, pmap, "hwhm_lorentz");
570
    descr.gaussWidth =
571
        ConvertToYSpace::getComponentParameter(*foil1, pmap, "sigma_gauss");
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
    m_foils1[i] = descr; // copy
  }

  if (g_log.is(Kernel::Logger::Priority::PRIO_INFORMATION)) {
    std::ostringstream os;
    os << "Instrument geometry:\n"
       << "  l1 = " << m_l1 << "m\n"
       << "  foil radius = " << m_foilRadius << "\n"
       << "  foil integration min = " << m_foilUpMin << "\n"
       << "  foil integration max = " << m_foilUpMax << "\n";
    std::ostringstream secondos;
    for (size_t i = 0; i < nfoils; ++i) {
      const auto &descr0 = m_foils0[i];
      os << "  foil theta range in position 0: theta_min=" << descr0.thetaMin
         << ", theta_max=" << descr0.thetaMax << "\n";
      const auto &descr1 = m_foils1[i];
      secondos << "  foil theta range in position 1: theta_min="
               << descr1.thetaMin << ", theta_max=" << descr1.thetaMax << "\n";
590
    }
591
592
593
594
595
    g_log.information() << os.str() << secondos.str();
  }
}

/**
LamarMoore's avatar
LamarMoore committed
596
597
598
599
600
601
602
 * @param foilComp A pointer to the foil component
 * @param radius The radius that gives the distance to the centre of the
 * bounding box
 * @param horizDir An enumeration defining which direction is horizontal
 * @return The min/max angle in theta(degrees) (horizontal direction if you
 * assume mid-point theta = 0)
 */
603
std::pair<double, double> VesuvioCalculateGammaBackground::calculateThetaRange(
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
    const Geometry::IComponent_const_sptr &foilComp, const double radius,
    const unsigned int horizDir) const {
  auto shapedObject =
      boost::dynamic_pointer_cast<const Geometry::IObjComponent>(foilComp);
  if (!shapedObject) {
    throw std::invalid_argument("A foil has been defined without a shape. "
                                "Please check instrument definition.");
  }

  // First get current theta position
  auto pos = foilComp->getPos();
  double theta(0.0), dummy(0.0);
  pos.getSpherical(dummy, theta, dummy); // absolute angle values
  if (pos[horizDir] < 0.0)
    theta *= -1.0; // negative quadrant for theta

  // Compute dtheta from bounding box & radius
  const auto &box = shapedObject->shape()->getBoundingBox();
  // box has center at 0,0,0
  double xmax = box.maxPoint()[0];
  double dtheta = std::asin(xmax / radius) * 180.0 / M_PI; // degrees
  return std::make_pair(theta - dtheta, theta + dtheta);
}

628
} // namespace Algorithms
629
} // namespace CurveFitting
630
} // namespace Mantid