SofQWNormalisedPolygon.cpp 15.4 KB
Newer Older
1
2
3
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
4
#include "MantidAlgorithms/SofQWNormalisedPolygon.h"
5
#include "MantidAlgorithms/SofQW.h"
6
#include "MantidAPI/BinEdgeAxis.h"
Dan Nixon's avatar
Dan Nixon committed
7
#include "MantidAPI/SpectrumDetectorMapping.h"
8
#include "MantidAPI/WorkspaceFactory.h"
9
#include "MantidDataObjects/FractionalRebinning.h"
10
#include "MantidGeometry/Instrument.h"
11
#include "MantidGeometry/Instrument/DetectorGroup.h"
12
#include "MantidGeometry/Instrument/ReferenceFrame.h"
13
14
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VectorHelper.h"
15

16
17
18
namespace Mantid {
namespace Algorithms {
// Setup typedef for later use
19
typedef std::map<specnum_t, Mantid::Kernel::V3D> SpectraDistanceMap;
20
21
22
typedef Geometry::IDetector_const_sptr DetConstPtr;

// Register the algorithm into the AlgorithmFactory
23
DECLARE_ALGORITHM(SofQWNormalisedPolygon)
24
25
26
27
28
29
30

using namespace Mantid::Kernel;
using namespace Mantid::Geometry;
using namespace Mantid::API;
using namespace Mantid::DataObjects;

/// Default constructor
31
SofQWNormalisedPolygon::SofQWNormalisedPolygon()
32
33
34
35
36
37
38
    : Rebin2D(), m_Qout(), m_thetaWidth(0.0), m_detNeighbourOffset(-1) {}

//----------------------------------------------------------------------------------------------

/**
 * @return the name of the Algorithm
 */
39
40
41
const std::string SofQWNormalisedPolygon::name() const {
  return "SofQWNormalisedPolygon";
}
42
43
44
45

/**
 * @return the version number of the Algorithm
 */
46
int SofQWNormalisedPolygon::version() const { return 1; }
47
48
49
50

/**
 * @return the category list for the Algorithm
 */
51
const std::string SofQWNormalisedPolygon::category() const {
Nick Draper's avatar
Nick Draper committed
52
  return "Inelastic\\SofQW";
53
}
54
55
56
57

/**
 * Initialize the algorithm
 */
58
59
60
void SofQWNormalisedPolygon::init() {
  SofQW::createCommonInputProperties(*this);
}
61
62
63
64

/**
 * Execute the algorithm.
 */
65
void SofQWNormalisedPolygon::exec() {
66
67
68
69
70
  MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
  // Do the full check for common binning
  if (!WorkspaceHelpers::commonBoundaries(inputWS)) {
    throw std::invalid_argument(
        "The input workspace must have common binning across all spectra");
71
72
  }

73
74
75
76
77
78
79
80
  RebinnedOutput_sptr outputWS =
      this->setUpOutputWorkspace(inputWS, getProperty("QAxisBinning"), m_Qout);
  g_log.debug() << "Workspace type: " << outputWS->id() << std::endl;
  setProperty("OutputWorkspace", outputWS);
  const size_t nEnergyBins = inputWS->blocksize();
  const size_t nHistos = inputWS->getNumberHistograms();

  // Holds the spectrum-detector mapping
81
  std::vector<specnum_t> specNumberMapping;
82
83
84
85
86
87
88
89
  std::vector<detid_t> detIDMapping;

  // Progress reports & cancellation
  const size_t nreports(nHistos * nEnergyBins);
  m_progress = boost::shared_ptr<API::Progress>(
      new API::Progress(this, 0.0, 1.0, nreports));

  // Compute input caches
90
  m_EmodeProperties.initCachedValues(*inputWS, this);
91
92
93
94
95
96
97
98
99
100

  std::vector<double> par =
      inputWS->getInstrument()->getNumberParameter("detector-neighbour-offset");
  if (par.empty()) {
    // Index theta cache
    this->initAngularCachesNonPSD(inputWS);
  } else {
    g_log.debug() << "Offset: " << par[0] << std::endl;
    this->m_detNeighbourOffset = static_cast<int>(par[0]);
    this->initAngularCachesPSD(inputWS);
101
102
  }

103
104
  const MantidVec &X = inputWS->readX(0);
  int emode = m_EmodeProperties.m_emode;
105
106

  PARALLEL_FOR2(inputWS, outputWS)
107
108
  for (int64_t i = 0; i < static_cast<int64_t>(nHistos);
       ++i) // signed for openmp
109
  {
110
    PARALLEL_START_INTERUPT_REGION
111

112
113
114
    DetConstPtr detector = inputWS->getDetector(i);
    if (detector->isMasked() || detector->isMonitor()) {
      continue;
115
116
    }

117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
    double theta = this->m_theta[i];
    double phi = this->m_phi[i];
    double thetaWidth = this->m_thetaWidths[i];
    double phiWidth = this->m_phiWidths[i];

    // Compute polygon points
    double thetaHalfWidth = 0.5 * thetaWidth;
    double phiHalfWidth = 0.5 * phiWidth;

    const double thetaLower = theta - thetaHalfWidth;
    const double thetaUpper = theta + thetaHalfWidth;

    const double phiLower = phi - phiHalfWidth;
    const double phiUpper = phi + phiHalfWidth;

132
    const double efixed = m_EmodeProperties.getEFixed(*detector);
133
    const specnum_t specNo = inputWS->getSpectrum(i)->getSpectrumNo();
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
    std::stringstream logStream;
    for (size_t j = 0; j < nEnergyBins; ++j) {
      m_progress->report("Computing polygon intersections");
      // For each input polygon test where it intersects with
      // the output grid and assign the appropriate weights of Y/E
      const double dE_j = X[j];
      const double dE_jp1 = X[j + 1];

      const double lrQ =
          this->calculateQ(efixed, emode, dE_jp1, thetaLower, phiLower);

      const V2D ll(dE_j,
                   this->calculateQ(efixed, emode, dE_j, thetaLower, phiLower));
      const V2D lr(dE_jp1, lrQ);
      const V2D ur(dE_jp1, this->calculateQ(efixed, emode, dE_jp1, thetaUpper,
                                            phiUpper));
      const V2D ul(dE_j,
                   this->calculateQ(efixed, emode, dE_j, thetaUpper, phiUpper));
      if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
        logStream << "Spectrum=" << specNo << ", theta=" << theta
                  << ",thetaWidth=" << thetaWidth << ", phi=" << phi
                  << ", phiWidth=" << phiWidth << ". QE polygon: ll=" << ll
                  << ", lr=" << lr << ", ur=" << ur << ", ul=" << ul << "\n";
157
158
      }

159
      Quadrilateral inputQ = Quadrilateral(ll, lr, ur, ul);
160

161
162
      FractionalRebinning::rebinToFractionalOutput(inputQ, inputWS, i, j,
                                                   outputWS, m_Qout);
Dan Nixon's avatar
Dan Nixon committed
163

164
165
166
167
168
      // Find which q bin this point lies in
      const MantidVec::difference_type qIndex =
          std::upper_bound(m_Qout.begin(), m_Qout.end(), lrQ) - m_Qout.begin();
      if (qIndex != 0 && qIndex < static_cast<int>(m_Qout.size())) {
        // Add this spectra-detector pair to the mapping
169
        PARALLEL_CRITICAL(SofQWNormalisedPolygon_spectramap) {
170
171
172
          specNumberMapping.push_back(
              outputWS->getSpectrum(qIndex - 1)->getSpectrumNo());
          detIDMapping.push_back(detector->getID());
173
        }
174
175
      }
    }
176
177
178
    if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
      g_log.debug(logStream.str());
    }
Dan Nixon's avatar
Dan Nixon committed
179

180
    PARALLEL_END_INTERUPT_REGION
181
  }
182
  PARALLEL_CHECK_INTERUPT_REGION
183
184

  outputWS->finalize();
185
  FractionalRebinning::normaliseOutput(outputWS, inputWS, m_progress);
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201

  // Set the output spectrum-detector mapping
  SpectrumDetectorMapping outputDetectorMap(specNumberMapping, detIDMapping);
  outputWS->updateSpectraUsing(outputDetectorMap);
}

/**
 * Calculate the Q value for a given set of energy transfer, scattering
 * and azimuthal angle.
 * @param efixed :: An fixed energy value
 * @param emode  :: the energy evaluation mode
 * @param deltaE :: The energy change
 * @param twoTheta :: The value of the scattering angle
 * @param azimuthal :: The value of the azimuthual angle
 * @return The value of Q
 */
202
203
204
205
double SofQWNormalisedPolygon::calculateQ(const double efixed, int emode,
                                          const double deltaE,
                                          const double twoTheta,
                                          const double azimuthal) const {
206
207
208
209
210
211
212
213
  double ki = 0.0;
  double kf = 0.0;
  if (emode == 1) {
    ki = std::sqrt(efixed * SofQW::energyToK());
    kf = std::sqrt((efixed - deltaE) * SofQW::energyToK());
  } else if (emode == 2) {
    ki = std::sqrt((deltaE + efixed) * SofQW::energyToK());
    kf = std::sqrt(efixed * SofQW::energyToK());
214
  }
215
216
217
218
219
220
221
222
223
224
225
226
227
  const double Qx = ki - kf * std::cos(twoTheta);
  const double Qy = -kf * std::sin(twoTheta) * std::cos(azimuthal);
  const double Qz = -kf * std::sin(twoTheta) * std::sin(azimuthal);
  return std::sqrt(Qx * Qx + Qy * Qy + Qz * Qz);
}
/**
 * A map detector ID and Q ranges
 * This method looks unnecessary as it could be calculated on the fly but
 * the parallelization means that lazy instantation slows it down due to the
 * necessary CRITICAL sections required to update the cache. The Q range
 * values are required very frequently so the total time is more than
 * offset by this precaching step
 */
228
void SofQWNormalisedPolygon::initAngularCachesNonPSD(
229
230
231
232
233
234
235
236
237
238
239
240
    const API::MatrixWorkspace_const_sptr &workspace) {
  const size_t nhist = workspace->getNumberHistograms();
  this->m_theta = std::vector<double>(nhist);
  this->m_thetaWidths = std::vector<double>(nhist);
  // Force phi widths to zero
  this->m_phi = std::vector<double>(nhist, 0.0);
  this->m_phiWidths = std::vector<double>(nhist, 0.0);

  auto inst = workspace->getInstrument();
  const PointingAlong upDir = inst->getReferenceFrame()->pointingUp();

  for (size_t i = 0; i < nhist; ++i) // signed for OpenMP
241
  {
242
243
244
245
246
247
    m_progress->report("Calculating detector angles");
    IDetector_const_sptr det;
    try {
      det = workspace->getDetector(i);
      // Check to see if there is an EFixed, if not skip it
      try {
248
        m_EmodeProperties.getEFixed(*det);
249
250
      } catch (std::runtime_error &) {
        det.reset();
251
      }
252
253
254
255
256
257
258
259
260
261
262
263
264
    } catch (Kernel::Exception::NotFoundError &) {
      // Catch if no detector. Next line tests whether this happened - test
      // placed
      // outside here because Mac Intel compiler doesn't like 'continue' in a
      // catch
      // in an openmp block.
    }
    // If no detector found, skip onto the next spectrum
    if (!det || det->isMonitor()) {
      this->m_theta[i] = -1.0; // Indicates a detector to skip
      this->m_thetaWidths[i] = -1.0;
      continue;
    }
265
    this->m_theta[i] = workspace->detectorTwoTheta(*det);
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

    /**
     * Determine width from shape geometry. A group is assumed to contain
     * detectors with the same shape & r, theta value, i.e. a ring mapped-group
     * The shape is retrieved and rotated to match the rotation of the detector.
     * The angular width is computed using the l2 distance from the sample
     */
    if (auto group = boost::dynamic_pointer_cast<const DetectorGroup>(det)) {
      // assume they all have same shape and same r,theta
      auto dets = group->getDetectors();
      det = dets[0];
    }
    const auto pos = det->getPos();
    double l2(0.0), t(0.0), p(0.0);
    pos.getSpherical(l2, t, p);
    // Get the shape
    auto shape =
        det->shape(); // Defined in its own reference frame with centre at 0,0,0
    auto rot = det->getRotation();
    BoundingBox bbox = shape->getBoundingBox();
    auto maxPoint(bbox.maxPoint());
    rot.rotate(maxPoint);
    double boxWidth = maxPoint[upDir];

    m_thetaWidths[i] = std::fabs(2.0 * std::atan(boxWidth / l2));
    if (g_log.is(Logger::Priority::PRIO_DEBUG)) {
      g_log.debug() << "Detector at spectrum ="
                    << workspace->getSpectrum(i)->getSpectrumNo()
                    << ", width=" << m_thetaWidths[i] * 180.0 / M_PI
                    << " degrees\n";
296
297
    }
  }
298
299
300
301
302
303
304
305
}

/**
 * Function that retrieves the two-theta and azimuthal angles from a given
 * detector. It then looks up the nearest neighbours. Using those detectors,
 * it calculates the two-theta and azimuthal angle widths.
 * @param workspace : the workspace containing the needed detector information
 */
306
307
void SofQWNormalisedPolygon::initAngularCachesPSD(
    const API::MatrixWorkspace_const_sptr &workspace) {
308
309
310
311
312
313
314
315
316
317
318
319
320
321
  // Trigger a build of the nearst neighbors outside the OpenMP loop
  const int numNeighbours = 4;
  const size_t nHistos = workspace->getNumberHistograms();
  g_log.debug() << "Number of Histograms: " << nHistos << std::endl;

  this->m_theta = std::vector<double>(nHistos);
  this->m_thetaWidths = std::vector<double>(nHistos);
  this->m_phi = std::vector<double>(nHistos);
  this->m_phiWidths = std::vector<double>(nHistos);

  for (size_t i = 0; i < nHistos; ++i) {
    m_progress->report("Calculating detector angular widths");
    DetConstPtr detector = workspace->getDetector(i);
    g_log.debug() << "Current histogram: " << i << std::endl;
322
    specnum_t inSpec = workspace->getSpectrum(i)->getSpectrumNo();
323
324
325
326
327
328
329
330
331
    SpectraDistanceMap neighbours =
        workspace->getNeighboursExact(inSpec, numNeighbours, true);

    g_log.debug() << "Current ID: " << inSpec << std::endl;
    // Convert from spectrum numbers to workspace indices
    double thetaWidth = -DBL_MAX;
    double phiWidth = -DBL_MAX;

    // Find theta and phi widths
332
    double theta = workspace->detectorTwoTheta(*detector);
333
334
    double phi = detector->getPhi();

335
336
337
338
    specnum_t deltaPlus1 = inSpec + 1;
    specnum_t deltaMinus1 = inSpec - 1;
    specnum_t deltaPlusT = inSpec + this->m_detNeighbourOffset;
    specnum_t deltaMinusT = inSpec - this->m_detNeighbourOffset;
339

340
    for (auto &neighbour : neighbours) {
341
      specnum_t spec = neighbour.first;
342
343
344
345
      g_log.debug() << "Neighbor ID: " << spec << std::endl;
      if (spec == deltaPlus1 || spec == deltaMinus1 || spec == deltaPlusT ||
          spec == deltaMinusT) {
        DetConstPtr detector_n = workspace->getDetector(spec - 1);
346
        double theta_n = workspace->detectorTwoTheta(*detector_n) * 0.5;
347
348
349
350
351
352
353
354
355
356
357
358
359
        double phi_n = detector_n->getPhi();

        double dTheta = std::fabs(theta - theta_n);
        double dPhi = std::fabs(phi - phi_n);
        if (dTheta > thetaWidth) {
          thetaWidth = dTheta;
          g_log.information()
              << "Current ThetaWidth: " << thetaWidth * 180 / M_PI << std::endl;
        }
        if (dPhi > phiWidth) {
          phiWidth = dPhi;
          g_log.information() << "Current PhiWidth: " << phiWidth * 180 / M_PI
                              << std::endl;
360
361
362
        }
      }
    }
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
    this->m_theta[i] = theta;
    this->m_phi[i] = phi;
    this->m_thetaWidths[i] = thetaWidth;
    this->m_phiWidths[i] = phiWidth;
  }
}

/** Creates the output workspace, setting the axes according to the input
 * binning parameters
 *  @param[in]  inputWorkspace The input workspace
 *  @param[in]  binParams The bin parameters from the user
 *  @param[out] newAxis        The 'vertical' axis defined by the given
 * parameters
 *  @return A pointer to the newly-created workspace
 */
378
379
380
RebinnedOutput_sptr SofQWNormalisedPolygon::setUpOutputWorkspace(
    API::MatrixWorkspace_const_sptr inputWorkspace,
    const std::vector<double> &binParams, std::vector<double> &newAxis) {
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
  // Create vector to hold the new X axis values
  MantidVecPtr xAxis;
  xAxis.access() = inputWorkspace->readX(0);
  const int xLength = static_cast<int>(xAxis->size());
  // Create a vector to temporarily hold the vertical ('y') axis and populate
  // that
  const int yLength = static_cast<int>(
      VectorHelper::createAxisFromRebinParams(binParams, newAxis));

  // Create the output workspace
  MatrixWorkspace_sptr temp = WorkspaceFactory::Instance().create(
      "RebinnedOutput", yLength - 1, xLength, xLength - 1);
  RebinnedOutput_sptr outputWorkspace =
      boost::static_pointer_cast<RebinnedOutput>(temp);
  WorkspaceFactory::Instance().initializeFromParent(inputWorkspace,
                                                    outputWorkspace, true);

  // Create a binned numeric axis to replace the default vertical one
  Axis *const verticalAxis = new BinEdgeAxis(newAxis);
  outputWorkspace->replaceAxis(1, verticalAxis);

  // Now set the axis values
  for (int i = 0; i < yLength - 1; ++i) {
    outputWorkspace->setX(i, xAxis);
405
406
  }

407
408
409
  // Set the axis units
  verticalAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
  verticalAxis->title() = "|Q|";
410

411
412
  // Set the X axis title (for conversion to MD)
  outputWorkspace->getAxis(0)->title() = "Energy transfer";
413

414
415
416
  outputWorkspace->setYUnit("");
  outputWorkspace->setYUnitLabel("Intensity");

417
418
  return outputWorkspace;
}
419

420
421
} // namespace Mantid
} // namespace Algorithms