IntegratePeaksMD2.cpp 41.6 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 +
Lynch, Vickie's avatar
Lynch, Vickie committed
7
8
#include "MantidMDAlgorithms/IntegratePeaksMD2.h"
#include "MantidAPI/AnalysisDataService.h"
9
#include "MantidAPI/Column.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
10
11
12
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionFactory.h"
13
14
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/IMDEventWorkspace.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
15
#include "MantidAPI/IPeakFunction.h"
16
#include "MantidAPI/Progress.h"
17
#include "MantidAPI/Run.h"
18
19
20
21
#include "MantidAPI/TableRow.h"
#include "MantidAPI/TextAxis.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/CoordTransformDistance.h"
22
#include "MantidDataObjects/MDBoxIterator.h"
23
24
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidDataObjects/Peak.h"
25
#include "MantidDataObjects/PeakShapeEllipsoid.h"
26
27
28
#include "MantidDataObjects/PeakShapeSpherical.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/Workspace2D.h"
LamarMoore's avatar
LamarMoore committed
29
#include "MantidGeometry/Instrument/DetectorInfo.h"
30
31
32
33
34
#include "MantidHistogramData/LinearGenerator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/System.h"
#include "MantidKernel/Utils.h"
#include "MantidMDAlgorithms/GSLFunctions.h"
35
#include "MantidMDAlgorithms/MDBoxMaskFunction.h"
36

37
#include <cmath>
Lynch, Vickie's avatar
Lynch, Vickie committed
38
#include <fstream>
LamarMoore's avatar
LamarMoore committed
39
#include <gsl/gsl_integration.h>
Lynch, Vickie's avatar
Lynch, Vickie committed
40

41
42
43
44
45
46
namespace Mantid {
namespace MDAlgorithms {

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(IntegratePeaksMD2)

47
using namespace Mantid::HistogramData;
48
49
using namespace Mantid::Kernel;
using namespace Mantid::API;
50
using namespace Mantid::DataObjects;
51
52
53
54
55
56
using namespace Mantid::DataObjects;
using namespace Mantid::Geometry;

/** Initialize the algorithm's properties.
 */
void IntegratePeaksMD2::init() {
57
  declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>(
58
                      "InputWorkspace", "", Direction::Input),
59
60
61
                  "An input MDEventWorkspace.");

  declareProperty(
62
      std::make_unique<PropertyWithValue<double>>("PeakRadius", 1.0,
Sam Jenkins's avatar
Sam Jenkins committed
63
                                                  Direction::Input),
64
65
66
67
      "Fixed radius around each peak position in which to integrate (in the "
      "same units as the workspace).");

  declareProperty(
68
      std::make_unique<PropertyWithValue<double>>("BackgroundInnerRadius", 0.0,
Sam Jenkins's avatar
Sam Jenkins committed
69
                                                  Direction::Input),
70
71
72
73
74
      "Inner radius to use to evaluate the background of the peak.\n"
      "If smaller than PeakRadius, then we assume BackgroundInnerRadius = "
      "PeakRadius.");

  declareProperty(
75
      std::make_unique<PropertyWithValue<double>>("BackgroundOuterRadius", 0.0,
Sam Jenkins's avatar
Sam Jenkins committed
76
                                                  Direction::Input),
77
78
79
80
81
82
      "Outer radius to use to evaluate the background of the peak.\n"
      "The signal density around the peak (BackgroundInnerRadius < r < "
      "BackgroundOuterRadius) is used to estimate the background under the "
      "peak.\n"
      "If smaller than PeakRadius, no background measurement is done.");

83
  declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>(
84
                      "PeaksWorkspace", "", Direction::Input),
85
86
87
                  "A PeaksWorkspace containing the peaks to integrate.");

  declareProperty(
88
      std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "",
Sam Jenkins's avatar
Sam Jenkins committed
89
                                                          Direction::Output),
90
91
92
93
94
95
96
97
98
99
100
101
102
      "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
      "with the peaks' integrated intensities.");

  declareProperty("ReplaceIntensity", true,
                  "Always replace intensity in PeaksWorkspacem (default).\n"
                  "If false, then do not replace intensity if calculated value "
                  "is 0 (used for SNSSingleCrystalReduction)");

  declareProperty(
      "IntegrateIfOnEdge", true,
      "Only warning if all of peak outer radius is not on detector (default).\n"
      "If false, do not integrate if the outer radius is not on a detector.");

Lynch, Vickie's avatar
Lynch, Vickie committed
103
  declareProperty("AdaptiveQBackground", false,
104
                  "Default is false.   If true, "
105
                  "BackgroundOuterRadius + AdaptiveQMultiplier * **|Q|** and "
106
                  "BackgroundInnerRadius + AdaptiveQMultiplier * **|Q|**");
107

108
109
110
  declareProperty("Ellipsoid", false, "Default is sphere.");
  declareProperty("FixQAxis", false, "Default is sphere.");

111
112
113
114
  declareProperty("Cylinder", false,
                  "Default is sphere.  Use next five parameters for cylinder.");

  declareProperty(
115
      std::make_unique<PropertyWithValue<double>>("CylinderLength", 0.0,
Sam Jenkins's avatar
Sam Jenkins committed
116
                                                  Direction::Input),
117
118
119
      "Length of cylinder in which to integrate (in the same units as the "
      "workspace).");

Sam Jenkins's avatar
Sam Jenkins committed
120
121
  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "PercentBackground", 0.0, Direction::Input),
122
                  "Percent of CylinderLength that is background (20 is 20%)");
123
124
125

  std::vector<std::string> peakNames =
      FunctionFactory::Instance().getFunctionNames<IPeakFunction>();
126
  peakNames.emplace_back("NoFit");
127
  declareProperty("ProfileFunction", "Gaussian",
128
                  std::make_shared<StringListValidator>(peakNames),
129
130
131
132
133
134
135
                  "Fitting function for profile that is used only with "
                  "Cylinder integration.");

  std::vector<std::string> integrationOptions(2);
  integrationOptions[0] = "Sum";
  integrationOptions[1] = "GaussianQuadrature";
  auto integrationvalidator =
136
      std::make_shared<StringListValidator>(integrationOptions);
137
138
139
140
141
142
  declareProperty("IntegrationOption", "GaussianQuadrature",
                  integrationvalidator,
                  "Integration method for calculating intensity "
                  "used only with Cylinder integration.");

  declareProperty(
Sam Jenkins's avatar
Sam Jenkins committed
143
144
145
      std::make_unique<FileProperty>("ProfilesFile", "",
                                     FileProperty::OptionalSave,
                                     std::vector<std::string>(1, "profiles")),
146
      "Save (Optionally) as Isaw peaks file with profiles included");
Lynch, Vickie's avatar
Lynch, Vickie committed
147
148

  declareProperty("AdaptiveQMultiplier", 0.0,
149
                  "PeakRadius + AdaptiveQMultiplier * **|Q|** "
150
                  "so each peak has a "
Lynch, Vickie's avatar
Lynch, Vickie committed
151
                  "different integration radius.  Q includes the 2*pi factor.");
152
153
154
155

  declareProperty(
      "CorrectIfOnEdge", false,
      "Only warning if all of peak outer radius is not on detector (default).\n"
Lynch, Vickie's avatar
Lynch, Vickie committed
156
157
      "If false, correct for volume off edge for both background and "
      "intensity.");
158

159
  declareProperty("UseOnePercentBackgroundCorrection", true,
160
161
162
                  "If this options is enabled, then the the top 1% of the "
                  "background will be removed"
                  "before the background subtraction.");
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
}

//----------------------------------------------------------------------------------------------
/** Integrate the peaks of the workspace using parameters saved in the algorithm
 * class
 * @param ws ::  MDEventWorkspace to integrate
 */
template <typename MDE, size_t nd>
void IntegratePeaksMD2::integrate(typename MDEventWorkspace<MDE, nd>::sptr ws) {
  if (nd != 3)
    throw std::invalid_argument("For now, we expect the input MDEventWorkspace "
                                "to have 3 dimensions only.");

  /// Peak workspace to integrate
  Mantid::DataObjects::PeaksWorkspace_sptr inPeakWS =
      getProperty("PeaksWorkspace");

  /// Output peaks workspace, create if needed
  Mantid::DataObjects::PeaksWorkspace_sptr peakWS =
      getProperty("OutputWorkspace");
  if (peakWS != inPeakWS)
184
    peakWS = inPeakWS->clone();
185
186
  // This only fails in the unit tests which say that MaskBTP is not registered
  try {
187
188
    runMaskDetectors(inPeakWS, "Tube", "edges");
    runMaskDetectors(inPeakWS, "Pixel", "edges");
189
  } catch (...) {
190
191
    g_log.error("Can't execute MaskBTP algorithm for this instrument to set "
                "edge for IntegrateIfOnEdge option");
Lynch, Vickie's avatar
Lynch, Vickie committed
192
193
  }

194
  calculateE1(inPeakWS->detectorInfo()); // fill E1Vec for use in detectorQ
195
196
  Mantid::Kernel::SpecialCoordinateSystem CoordinatesToUse =
      ws->getSpecialCoordinateSystem();
197
198
199
200
201
202
203

  /// Radius to use around peaks
  double PeakRadius = getProperty("PeakRadius");
  /// Background (end) radius
  double BackgroundOuterRadius = getProperty("BackgroundOuterRadius");
  /// Start radius of the background
  double BackgroundInnerRadius = getProperty("BackgroundInnerRadius");
204
  /// One percent background correction
205
206
  bool useOnePercentBackgroundCorrection =
      getProperty("UseOnePercentBackgroundCorrection");
207

208
209
  if (BackgroundInnerRadius < PeakRadius)
    BackgroundInnerRadius = PeakRadius;
210
  // Ellipsoid
211
212
  bool isEllipse = getProperty("Ellipsoid");
  bool qAxisIsFixed = getProperty("FixQAxis");
213
214
215
216
217
  /// Cylinder Length to use around peaks for cylinder
  double cylinderLength = getProperty("CylinderLength");
  Workspace2D_sptr wsProfile2D, wsFit2D, wsDiff2D;
  size_t numSteps = 0;
  bool cylinderBool = getProperty("Cylinder");
Lynch, Vickie's avatar
Lynch, Vickie committed
218
219
220
  bool adaptiveQBackground = getProperty("AdaptiveQBackground");
  double adaptiveQMultiplier = getProperty("AdaptiveQMultiplier");
  double adaptiveQBackgroundMultiplier = 0.0;
221
222
  if (adaptiveQBackground)
    adaptiveQBackgroundMultiplier = adaptiveQMultiplier;
223
224
225
226
227
228
229
230
231
232
  std::vector<double> PeakRadiusVector(peakWS->getNumberPeaks(), PeakRadius);
  std::vector<double> BackgroundInnerRadiusVector(peakWS->getNumberPeaks(),
                                                  BackgroundInnerRadius);
  std::vector<double> BackgroundOuterRadiusVector(peakWS->getNumberPeaks(),
                                                  BackgroundOuterRadius);
  if (cylinderBool) {
    numSteps = 100;
    size_t histogramNumber = peakWS->getNumberPeaks();
    Workspace_sptr wsProfile = WorkspaceFactory::Instance().create(
        "Workspace2D", histogramNumber, numSteps, numSteps);
233
    wsProfile2D = std::dynamic_pointer_cast<Workspace2D>(wsProfile);
234
235
236
    AnalysisDataService::Instance().addOrReplace("ProfilesData", wsProfile2D);
    Workspace_sptr wsFit = WorkspaceFactory::Instance().create(
        "Workspace2D", histogramNumber, numSteps, numSteps);
237
    wsFit2D = std::dynamic_pointer_cast<Workspace2D>(wsFit);
238
239
240
    AnalysisDataService::Instance().addOrReplace("ProfilesFit", wsFit2D);
    Workspace_sptr wsDiff = WorkspaceFactory::Instance().create(
        "Workspace2D", histogramNumber, numSteps, numSteps);
241
    wsDiff2D = std::dynamic_pointer_cast<Workspace2D>(wsDiff);
242
    AnalysisDataService::Instance().addOrReplace("ProfilesFitDiff", wsDiff2D);
243
244
245
    auto newAxis1 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
    auto newAxis2 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
    auto newAxis3 = std::make_unique<TextAxis>(peakWS->getNumberPeaks());
246
247
248
    auto newAxis1Raw = newAxis1.get();
    auto newAxis2Raw = newAxis2.get();
    auto newAxis3Raw = newAxis3.get();
249
250
251
    wsProfile2D->replaceAxis(1, std::move(newAxis1));
    wsFit2D->replaceAxis(1, std::move(newAxis2));
    wsDiff2D->replaceAxis(1, std::move(newAxis3));
252
    for (int i = 0; i < peakWS->getNumberPeaks(); ++i) {
253
      // Get a direct ref to that peak.
254
      IPeak &p = peakWS->getPeak(i);
255
      std::ostringstream label;
256
257
      label << Utils::round(p.getH()) << "_" << Utils::round(p.getK()) << "_"
            << Utils::round(p.getL()) << "_" << p.getRunNumber();
258
259
260
      newAxis1Raw->setLabel(i, label.str());
      newAxis2Raw->setLabel(i, label.str());
      newAxis3Raw->setLabel(i, label.str());
Lynch, Vickie's avatar
Lynch, Vickie committed
261
    }
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
  }
  double percentBackground = getProperty("PercentBackground");
  size_t peakMin = 0;
  size_t peakMax = numSteps;
  double ratio = 0.0;
  if (cylinderBool) {
    peakMin = static_cast<size_t>(static_cast<double>(numSteps) *
                                  percentBackground / 100.);
    peakMax = numSteps - peakMin - 1;
    size_t numPeakCh = peakMax - peakMin + 1; // number of peak channels
    size_t numBkgCh = numSteps - numPeakCh;   // number of background channels
    ratio = static_cast<double>(numPeakCh) / static_cast<double>(numBkgCh);
  }
  /// Replace intensity with 0
  bool replaceIntensity = getProperty("ReplaceIntensity");
  bool integrateEdge = getProperty("IntegrateIfOnEdge");
278
  bool correctEdge = getProperty("CorrectIfOnEdge");
279

280
281
  std::string profileFunction = getProperty("ProfileFunction");
  std::string integrationOption = getProperty("IntegrationOption");
282
  std::ofstream out;
283
  if (cylinderBool && profileFunction != "NoFit") {
284
285
286
    std::string outFile = getProperty("InputWorkspace");
    outFile.append(profileFunction);
    outFile.append(".dat");
287
288
    std::string save_path =
        ConfigService::Instance().getString("defaultsave.directory");
289
290
    outFile = save_path + outFile;
    out.open(outFile.c_str(), std::ofstream::out);
291
  }
292
  // volume of Background sphere with inner volume subtracted
LamarMoore's avatar
LamarMoore committed
293
294
295
  double volumeBkg =
      4.0 / 3.0 * M_PI *
      (std::pow(BackgroundOuterRadius, 3) - std::pow(BackgroundOuterRadius, 3));
296
  // volume of PeakRadius sphere
Lynch, Vickie's avatar
Lynch, Vickie committed
297
  double volumeRadius = 4.0 / 3.0 * M_PI * std::pow(PeakRadius, 3);
298
299
300
301
302
303
304
305
306
307
  //
  // If the following OMP pragma is included, this algorithm seg faults
  // sporadically when processing multiple TOPAZ runs in a script, on
  // Scientific Linux 6.2.  Typically, it seg faults after 2 to 6 runs are
  // processed, though occasionally it will process all 8 requested in the
  // script without crashing.  Since the lower level codes already use OpenMP,
  // parallelizing at this level is only marginally useful, giving about a
  // 5-10% speedup.  Perhaps is should just be removed permanantly, but for
  // now it is commented out to avoid the seg faults.  Refs #5533
  // PRAGMA_OMP(parallel for schedule(dynamic, 10) )
308
  // Initialize progress reporting
309
  int nPeaks = peakWS->getNumberPeaks();
310
  Progress progress(this, 0., 1., nPeaks);
311
312
313
  for (int i = 0; i < nPeaks; ++i) {
    if (this->getCancel())
      break; // User cancellation
314
    progress.report();
315

316
317
318
319
320
    // Get a direct ref to that peak.
    IPeak &p = peakWS->getPeak(i);

    // Get the peak center as a position in the dimensions of the workspace
    V3D pos;
321
    if (CoordinatesToUse == Mantid::Kernel::QLab) //"Q (lab frame)"
322
      pos = p.getQLabFrame();
323
    else if (CoordinatesToUse == Mantid::Kernel::QSample) //"Q (sample frame)"
324
      pos = p.getQSampleFrame();
325
    else if (CoordinatesToUse == Mantid::Kernel::HKL) //"HKL"
326
327
328
329
      pos = p.getHKL();

    // Do not integrate if sphere is off edge of detector

Lynch, Vickie's avatar
Lynch, Vickie committed
330
331
    double edge = detectorQ(p.getQLabFrame(),
                            std::max(BackgroundOuterRadius, PeakRadius));
332
    if (edge < std::max(BackgroundOuterRadius, PeakRadius)) {
333
      g_log.warning() << "Warning: sphere/cylinder for integration is off edge "
LamarMoore's avatar
LamarMoore committed
334
335
                         "of detector for peak "
                      << i << "; radius of edge =  " << edge << '\n';
336
337
338
339
      if (!integrateEdge) {
        if (replaceIntensity) {
          p.setIntensity(0.0);
          p.setSigmaIntensity(0.0);
Lynch, Vickie's avatar
Lynch, Vickie committed
340
        }
341
        continue;
Lynch, Vickie's avatar
Lynch, Vickie committed
342
      }
343
344
345
346
347
348
349
350
351
    }

    // Build the sphere transformation
    bool dimensionsUsed[nd];
    coord_t center[nd];
    for (size_t d = 0; d < nd; ++d) {
      dimensionsUsed[d] = true; // Use all dimensions
      center[d] = static_cast<coord_t>(pos[d]);
    }
352
353
354
355
    signal_t signal = 0;
    signal_t errorSquared = 0;
    signal_t bgSignal = 0;
    signal_t bgErrorSquared = 0;
356
    double background_total = 0.0;
357
    if (!cylinderBool) {
358
      // modulus of Q
Lynch, Vickie's avatar
Lynch, Vickie committed
359
      coord_t lenQpeak = 0.0;
360
      if (adaptiveQMultiplier != 0.0) {
361
        lenQpeak = 0.0;
362
        for (size_t d = 0; d < nd; ++d) {
363
364
365
366
          lenQpeak += center[d] * center[d];
        }
        lenQpeak = std::sqrt(lenQpeak);
      }
367
      double adaptiveRadius = adaptiveQMultiplier * lenQpeak + PeakRadius;
368
      if (adaptiveRadius <= 0.0) {
369
370
        g_log.error() << "Error: Radius for integration sphere of peak " << i
                      << " is negative =  " << adaptiveRadius << '\n';
371
372
373
374
375
376
377
378
        adaptiveRadius = 0.;
        p.setIntensity(0.0);
        p.setSigmaIntensity(0.0);
        PeakRadiusVector[i] = 0.0;
        BackgroundInnerRadiusVector[i] = 0.0;
        BackgroundOuterRadiusVector[i] = 0.0;
        continue;
      }
379
      PeakRadiusVector[i] = adaptiveRadius;
380
381
382
383
      BackgroundInnerRadiusVector[i] =
          adaptiveQBackgroundMultiplier * lenQpeak + BackgroundInnerRadius;
      BackgroundOuterRadiusVector[i] =
          adaptiveQBackgroundMultiplier * lenQpeak + BackgroundOuterRadius;
384
385
386
      // define the radius squared for a sphere intially
      CoordTransformDistance getRadiusSq(nd, center, dimensionsUsed);
      // set spherical shape
387
      if (auto *shapeablePeak = dynamic_cast<Peak *>(&p)) {
388
        PeakShape *sphereShape = new PeakShapeSpherical(
389
390
391
            PeakRadiusVector[i], BackgroundInnerRadiusVector[i],
            BackgroundOuterRadiusVector[i], CoordinatesToUse, this->name(),
            this->version());
392
        shapeablePeak->setPeakShape(sphereShape);
393
      }
394
395
396
397
      const double scaleFactor = pow(PeakRadiusVector[i], 3) /
                                 (pow(BackgroundOuterRadiusVector[i], 3) -
                                  pow(BackgroundInnerRadiusVector[i], 3));
      // Integrate spherical background shell if specified
398
      if (BackgroundOuterRadius > PeakRadius) {
399
        // Get the total signal inside background shell
400
        ws->getBox()->integrateSphere(
401
402
            getRadiusSq,
            static_cast<coord_t>(pow(BackgroundOuterRadiusVector[i], 2)),
403
            bgSignal, bgErrorSquared,
404
            static_cast<coord_t>(pow(BackgroundInnerRadiusVector[i], 2)),
405
            useOnePercentBackgroundCorrection);
406
        // correct bg signal by Vpeak/Vshell (same for sphere and ellipse)
407
408
        bgSignal *= scaleFactor;
        bgErrorSquared *= scaleFactor * scaleFactor;
409
      }
410
411
      // if ellipsoid find covariance and centroid in spherical region
      // using one-pass algorithm from https://doi.org/10.1145/359146.359153
412
      if (isEllipse) {
RichardWaiteSTFC's avatar
RichardWaiteSTFC committed
413
        // flat bg to subtract
414
415
416
417
418
419
        const auto bgDensity =
            bgSignal / (4 * M_PI * pow(PeakRadiusVector[i], 3) / 3);
        std::vector<V3D> eigenvects;
        std::vector<double> eigenvals;
        findEllipsoid<MDE, nd>(
            ws, getRadiusSq, pos,
420
            static_cast<coord_t>(pow(PeakRadiusVector[i], 2)), qAxisIsFixed,
421
            bgDensity, eigenvects, eigenvals);
422

423
424
425
426
427
428
429
430
431
432
        // transform ellispoid onto sphere of radius = R
        getRadiusSq =
            CoordTransformDistance(nd, center, dimensionsUsed, 1, /* outD */
                                   eigenvects, eigenvals);
        // Integrate ellipsoid background shell if specified
        if (BackgroundOuterRadius > PeakRadius) {
          // Get the total signal inside "BackgroundOuterRadius"
          bgSignal = 0;
          bgErrorSquared = 0;
          ws->getBox()->integrateSphere(
RichardWaiteSTFC's avatar
RichardWaiteSTFC committed
433
434
435
436
              getRadiusSq,
              static_cast<coord_t>(pow(BackgroundOuterRadiusVector[i], 2)),
              bgSignal, bgErrorSquared,
              static_cast<coord_t>(pow(BackgroundInnerRadiusVector[i], 2)),
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
              useOnePercentBackgroundCorrection);
          // correct bg signal by Vpeak/Vshell (same as previously
          // calculated for sphere)
          bgSignal *= scaleFactor;
          bgErrorSquared *= scaleFactor * scaleFactor;
        }
        // set peak shape
        if (auto *shapeablePeak = dynamic_cast<Peak *>(&p)) {
          // get radii in same proprtion as eigenvalues
          auto max_stdev =
              pow(*std::max_element(eigenvals.begin(), eigenvals.end()), 0.5);
          std::vector<double> peakRadii(3, 0.0);
          std::vector<double> backgroundInnerRadii(3, 0.0);
          std::vector<double> backgroundOuterRadii(3, 0.0);
          for (size_t irad = 0; irad < peakRadii.size(); irad++) {
            auto scale = pow(eigenvals[irad], 0.5) / max_stdev;
            peakRadii[irad] = PeakRadiusVector[i] * scale;
            backgroundInnerRadii[irad] = BackgroundInnerRadiusVector[i] * scale;
            backgroundOuterRadii[irad] = BackgroundOuterRadiusVector[i] * scale;
          }
          PeakShape *ellipsoidShape = new PeakShapeEllipsoid(
              eigenvects, peakRadii, backgroundInnerRadii, backgroundOuterRadii,
              CoordinatesToUse, this->name(), this->version());
          shapeablePeak->setPeakShape(ellipsoidShape);
        }
462
      }
463
464
465
466
467
      // spherical integration of signal
      ws->getBox()->integrateSphere(
          getRadiusSq, static_cast<coord_t>(adaptiveRadius * adaptiveRadius),
          signal, errorSquared, 0.0 /* innerRadiusSquared */,
          useOnePercentBackgroundCorrection);
468
    } else {
469
470
471
      CoordTransformDistance cylinder(nd, center, dimensionsUsed, 2);

      // Perform the integration into whatever box is contained within.
472
473
      Counts signal_fit(numSteps);
      signal_fit = 0;
474

475
476
477
478
      ws->getBox()->integrateCylinder(
          cylinder, static_cast<coord_t>(PeakRadius),
          static_cast<coord_t>(cylinderLength), signal, errorSquared,
          signal_fit.mutableRawData());
479
480

      // Integrate around the background radius
481
      if (BackgroundOuterRadius > PeakRadius) {
482
        // Get the total signal inside "BackgroundOuterRadius"
483
        signal_fit = 0;
484

485
486
487
        ws->getBox()->integrateCylinder(
            cylinder, static_cast<coord_t>(BackgroundOuterRadius),
            static_cast<coord_t>(cylinderLength), bgSignal, bgErrorSquared,
488
            signal_fit.mutableRawData());
489

490
        Points points(signal_fit.size(), LinearGenerator(0, 1));
491
        wsProfile2D->setHistogram(i, points, signal_fit);
492
493
494
495
496
497

        // Evaluate the signal inside "BackgroundInnerRadius"
        signal_t interiorSignal = 0;
        signal_t interiorErrorSquared = 0;

        // Integrate this 3rd radius, if needed
498
499
500
501
        if (BackgroundInnerRadius != PeakRadius) {
          ws->getBox()->integrateCylinder(
              cylinder, static_cast<coord_t>(BackgroundInnerRadius),
              static_cast<coord_t>(cylinderLength), interiorSignal,
502
              interiorErrorSquared, signal_fit.mutableRawData());
503
        } else {
504
505
          // PeakRadius == BackgroundInnerRadius, so use the previous
          // value
506
507
508
          interiorSignal = signal;
          interiorErrorSquared = errorSquared;
        }
509
510
511
        // Subtract the peak part to get the intensity in the shell
        // (BackgroundInnerRadius < r < BackgroundOuterRadius)
        bgSignal -= interiorSignal;
512
513
514
        // We can subtract the error (instead of adding) because the two
        // values are 100% dependent; this is the same as integrating a
        // shell.
515
516
        bgErrorSquared -= interiorErrorSquared;
        // Relative volume of peak vs the BackgroundOuterRadius cylinder
517
518
        const double radiusRatio = (PeakRadius / BackgroundOuterRadius);
        const double peakVolume = radiusRatio * radiusRatio * cylinderLength;
519

520
521
        // Relative volume of the interior of the shell vs overall
        // background
522
523
        const double interiorRatio =
            (BackgroundInnerRadius / BackgroundOuterRadius);
524
525
        // Volume of the bg shell, relative to the volume of the
        // BackgroundOuterRadius cylinder
526
527
        const double bgVolume =
            1.0 - interiorRatio * interiorRatio * cylinderLength;
528
529
530

        // Finally, you will multiply the bg intensity by this to get the
        // estimated background under the peak volume
531
        const double scaleFactor = peakVolume / bgVolume;
532
533
534
        bgSignal *= scaleFactor;
        bgErrorSquared *= scaleFactor * scaleFactor;
      } else {
535
        Points points(signal_fit.size(), LinearGenerator(0, 1));
536
        wsProfile2D->setHistogram(i, points, signal_fit);
537
538
      }

539
      if (profileFunction == "NoFit") {
540
        signal = 0.;
541
542
        for (size_t j = 0; j < numSteps; j++) {
          if (j < peakMin || j > peakMax)
543
            background_total = background_total + wsProfile2D->mutableY(i)[j];
544
          else
545
            signal = signal + wsProfile2D->mutableY(i)[j];
546
        }
547
        errorSquared = std::fabs(signal);
548
      } else {
Lynch, Vickie's avatar
Lynch, Vickie committed
549

Lynch, Vickie's avatar
Lynch, Vickie committed
550
551
552
553
554
555
556
557
558
559
560
561
562
        IAlgorithm_sptr fitAlgorithm =
            createChildAlgorithm("Fit", -1, -1, false);
        // fitAlgorithm->setProperty("CreateOutput", true);
        // fitAlgorithm->setProperty("Output", "FitPeaks1D");
        std::string myFunc =
            std::string("name=LinearBackground;name=") + profileFunction;
        auto maxPeak = std::max_element(signal_fit.begin(), signal_fit.end());

        std::ostringstream strs;
        strs << maxPeak[0];
        std::string strMax = strs.str();
        if (profileFunction == "Gaussian") {
          myFunc += ", PeakCentre=50, Height=" + strMax;
563
          fitAlgorithm->setProperty("Constraints", "40<f1.PeakCentre<60");
Lynch, Vickie's avatar
Lynch, Vickie committed
564
565
566
        } else if (profileFunction == "BackToBackExponential" ||
                   profileFunction == "IkedaCarpenterPV") {
          myFunc += ", X0=50, I=" + strMax;
567
          fitAlgorithm->setProperty("Constraints", "40<f1.X0<60");
Lynch, Vickie's avatar
Lynch, Vickie committed
568
569
570
571
572
        }
        fitAlgorithm->setProperty("CalcErrors", true);
        fitAlgorithm->setProperty("Function", myFunc);
        fitAlgorithm->setProperty("InputWorkspace", wsProfile2D);
        fitAlgorithm->setProperty("WorkspaceIndex", static_cast<int>(i));
573
        try {
Lynch, Vickie's avatar
Lynch, Vickie committed
574
          fitAlgorithm->executeAsChildAlg();
575
        } catch (...) {
Lynch, Vickie's avatar
Lynch, Vickie committed
576
          g_log.error("Can't execute Fit algorithm");
577
578
          continue;
        }
579

Lynch, Vickie's avatar
Lynch, Vickie committed
580
        IFunction_sptr ifun = fitAlgorithm->getProperty("Function");
Lynch, Vickie's avatar
Lynch, Vickie committed
581
        if (i == 0) {
Lynch, Vickie's avatar
Lynch, Vickie committed
582
583
          out << std::setw(20) << "spectrum"
              << " ";
Lynch, Vickie's avatar
Lynch, Vickie committed
584
585
          for (size_t j = 0; j < ifun->nParams(); ++j)
            out << std::setw(20) << ifun->parameterName(j) << " ";
Lynch, Vickie's avatar
Lynch, Vickie committed
586
587
          out << std::setw(20) << "chi2"
              << " ";
Lynch, Vickie's avatar
Lynch, Vickie committed
588
589
          out << "\n";
        }
590
        out << std::setw(20) << i << " ";
Lynch, Vickie's avatar
Lynch, Vickie committed
591
592
593
594
        for (size_t j = 0; j < ifun->nParams(); ++j) {
          out << std::setw(20) << std::fixed << std::setprecision(10)
              << ifun->getParameter(j) << " ";
        }
595
        double chi2 = fitAlgorithm->getProperty("OutputChi2overDoF");
Lynch, Vickie's avatar
Lynch, Vickie committed
596
597
        out << std::setw(20) << std::fixed << std::setprecision(10) << chi2
            << "\n";
598

599
600
        std::shared_ptr<const CompositeFunction> fun =
            std::dynamic_pointer_cast<const CompositeFunction>(ifun);
601

602
        const auto &x = wsProfile2D->x(i);
603
604
605
        wsFit2D->setSharedX(i, wsProfile2D->sharedX(i));
        wsDiff2D->setSharedX(i, wsProfile2D->sharedX(i));

606
        FunctionDomain1DVector domain(x.rawData());
607
608
        FunctionValues yy(domain);
        fun->function(domain, yy);
609
        auto funcValues = yy.toVector();
610

611
        wsFit2D->mutableY(i) = std::move(funcValues);
612
        wsDiff2D->setSharedY(i, wsProfile2D->sharedY(i));
613
        wsDiff2D->mutableY(i) -= wsFit2D->y(i);
614

615
        // Calculate intensity
616
        signal = 0.0;
617
        if (integrationOption == "Sum") {
618
          for (size_t j = peakMin; j <= peakMax; j++)
619
            if (std::isfinite(yy[j]))
620
621
622
              signal += yy[j];
        } else {
          gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
623
624
625
626
627
628
629

          double error;

          gsl_function F;
          F.function = &Mantid::MDAlgorithms::f_eval2;
          F.params = &fun;

630
631
          gsl_integration_qags(&F, x[peakMin], x[peakMax], 0, 1e-7, 1000, w,
                               &signal, &error);
632

633
          gsl_integration_workspace_free(w);
634
635
636
        }
        errorSquared = std::fabs(signal);
        // Get background counts
637
        for (size_t j = 0; j < numSteps; j++) {
Lynch, Vickie's avatar
Lynch, Vickie committed
638
639
          double background =
              ifun->getParameter(0) + ifun->getParameter(1) * x[j];
640
641
642
643
          if (j < peakMin || j > peakMax)
            background_total = background_total + background;
        }
      }
Lynch, Vickie's avatar
Lynch, Vickie committed
644
    }
645
646
647
648
649
    checkOverlap(
        i, peakWS, CoordinatesToUse,
        2.0 * std::max(PeakRadiusVector[i], BackgroundOuterRadiusVector[i]));
    // Save it back in the peak object.
    if (signal != 0. || replaceIntensity) {
650
      double edgeMultiplier = 1.0;
651
      double peakMultiplier = 1.0;
Lynch, Vickie's avatar
Lynch, Vickie committed
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
      if (correctEdge) {
        if (edge < BackgroundOuterRadius) {
          double e1 = BackgroundOuterRadius - edge;
          // volume of cap of sphere with h = edge
          double f1 =
              M_PI * std::pow(e1, 2) / 3 * (3 * BackgroundOuterRadius - e1);
          edgeMultiplier = volumeBkg / (volumeBkg - f1);
        }
        if (edge < PeakRadius) {
          double sigma = PeakRadius / 3.0;
          // assume gaussian peak
          double e1 =
              std::exp(-std::pow(edge, 2) / (2 * sigma * sigma)) * PeakRadius;
          // volume of cap of sphere with h = edge
          double f1 = M_PI * std::pow(e1, 2) / 3 * (3 * PeakRadius - e1);
          peakMultiplier = volumeRadius / (volumeRadius - f1);
        }
669
      }
670

Lynch, Vickie's avatar
Lynch, Vickie committed
671
672
673
674
675
676
      p.setIntensity(peakMultiplier * signal -
                     edgeMultiplier * (ratio * background_total + bgSignal));
      p.setSigmaIntensity(
          sqrt(peakMultiplier * errorSquared +
               edgeMultiplier * (ratio * ratio * std::fabs(background_total) +
                                 bgErrorSquared)));
677
678
679
680
681
682
683
684
    }

    g_log.information() << "Peak " << i << " at " << pos << ": signal "
                        << signal << " (sig^2 " << errorSquared
                        << "), with background "
                        << bgSignal + ratio * background_total << " (sig^2 "
                        << bgErrorSquared +
                               ratio * ratio * std::fabs(background_total)
685
                        << ") subtracted.\n";
686
  }
687
688
  // This flag is used by the PeaksWorkspace to evaluate whether it has
  // been integrated.
689
690
691
692
693
694
695
696
697
698
699
  peakWS->mutableRun().addProperty("PeaksIntegrated", 1, true);
  // These flags are specific to the algorithm.
  peakWS->mutableRun().addProperty("PeakRadius", PeakRadiusVector, true);
  peakWS->mutableRun().addProperty("BackgroundInnerRadius",
                                   BackgroundInnerRadiusVector, true);
  peakWS->mutableRun().addProperty("BackgroundOuterRadius",
                                   BackgroundOuterRadiusVector, true);

  // save profiles in peaks file
  const std::string outfile = getProperty("ProfilesFile");
  if (outfile.length() > 0) {
700
    IAlgorithm_sptr alg;
701
702
    try {
      alg = createChildAlgorithm("SaveIsawPeaks", -1, -1, false);
703
    } catch (Exception::NotFoundError &) {
704
705
      g_log.error("Can't locate SaveIsawPeaks algorithm");
      throw;
706
707
708
709
710
    }
    alg->setProperty("InputWorkspace", peakWS);
    alg->setProperty("ProfileWorkspace", wsProfile2D);
    alg->setPropertyValue("Filename", outfile);
    alg->execute();
Lynch, Vickie's avatar
Lynch, Vickie committed
711
  }
712
713
714
715
  // Save the output
  setProperty("OutputWorkspace", peakWS);
}

716
717
718
719
720
/**
 * Calculate the covariance matrix of a spherical region and store the
 * eigenvectors and eigenvalues that diagonalise the covariance matrix in the
 * vectors provided
 *
RichardWaiteSTFC's avatar
RichardWaiteSTFC committed
721
722
723
724
 *  @param ws             input workspace
 *  @param getRadiusSq    Coord transfrom for sphere
 *  @param pos            V3D of peak centre
 *  @param radiusSquared  radius that defines spherical region for covarariance
725
 *  @param qAxisIsFixed      bool to fix an eigenvector along direction pos
RichardWaiteSTFC's avatar
RichardWaiteSTFC committed
726
727
728
 *  @param bgDensity      background counts per unit volume
 *  @param eigenvects     eigenvectors of covariance matrix of spherical region
 *  @param eigenvals      eigenvectors of covariance matrix of spherical region
729
730
731
732
733
 */
template <typename MDE, size_t nd>
void IntegratePeaksMD2::findEllipsoid(
    typename MDEventWorkspace<MDE, nd>::sptr ws,
    const CoordTransform &getRadiusSq, const V3D &pos,
734
    const coord_t &radiusSquared, const bool &qAxisIsFixed,
735
736
737
738
739
740
    const double &bgDensity, std::vector<V3D> &eigenvects,
    std::vector<double> &eigenvals) {
  double w_sum = 0.0;  // sum of weights
  Matrix<double> Evec; // hold eigenvectors
  Matrix<double> Eval; // hold eigenvals in diag
  // get leaf-only iterators over all boxes in ws
RichardWaiteSTFC's avatar
RichardWaiteSTFC committed
741
742
  auto function = std::make_unique<Geometry::MDAlgorithms::MDBoxMaskFunction>(
      pos, radiusSquared);
743
744
  MDBoxBase<MDE, nd> *baseBox = ws->getBox();
  MDBoxIterator<MDE, nd> MDiter(baseBox, 1000, true, function.get());
745
746
747
748
749
750
751

  Matrix<double> cov_mat;
  Matrix<double> Pinv(3, 3);
  if (qAxisIsFixed) {
    // 2D covar in plane perp to Q (uhat,vhat basis)
    cov_mat = Matrix<double>(2, 2);
    // transformation from Qlab to Qhat, vhat and uhat,
752
    getPinv(pos, Pinv);
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
  } else {
    cov_mat = Matrix<double>(3, 3);
  }
  std::vector<double> mean(3, 0.0);
  double var_Qhat = 0.0; //  variance parallel to Q (used if fix Q axis)

  // loop over all boxes inside radius
  do {
    auto *box = dynamic_cast<MDBox<MDE, nd> *>(MDiter.getBox());
    if (box && !box->getIsMasked()) {
      const std::vector<MDE> &events = box->getConstEvents();
      auto bg = bgDensity / (static_cast<double>(events.size()) *
                             (box->getInverseVolume()));
      // For each event
      for (const auto &evnt : events) {
        std::vector<coord_t> center(nd);
        for (size_t d = 0; d < nd; ++d) {
          center[d] = evnt.getCenter(d);
        }
        coord_t out[nd];
        getRadiusSq.apply(center.data(), out);

        if (evnt.getSignal() > bg && out[0] < radiusSquared) {
          auto signal = (evnt.getSignal() - bg);
          w_sum += signal;

          if (qAxisIsFixed) {
            // transform coords to Q, uhat, vhat basis
            // use V3D for matrix algebra
            auto tmp = Pinv * V3D(static_cast<double>(center[0]),
RichardWaiteSTFC's avatar
RichardWaiteSTFC committed
783
784
                                  static_cast<double>(center[1]),
                                  static_cast<double>(center[2]));
785
786
            for (size_t d = 0; d < center.size(); ++d) {
              center[d] = tmp[d];
787
            }
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
          }

          // update mean
          for (size_t d = 0; d < mean.size(); ++d) {
            mean[d] += (signal / w_sum) * (center[d] - mean[d]);
          }

          if (qAxisIsFixed) {
            // get variance along Q
            var_Qhat += signal * pow((center[0] - mean[0]), 2);
          }
          for (size_t row = 0; row < cov_mat.numRows(); ++row) {
            for (size_t col = 0; col < cov_mat.numRows(); ++col) {
              // symmeteric matrix
              if (row <= col) {
                double cov = 0.0;
                if (!qAxisIsFixed) {
                  cov = signal * (center[row] - mean[row]) *
                        (center[col] - mean[col]);
                } else {
                  cov = signal * (center[row + 1] - mean[row + 1]) *
                        (center[col + 1] - mean[col + 1]);
                }
                if (row == col) {
                  cov_mat[row][col] += cov;
                } else {
                  cov_mat[row][col] += cov;
                  cov_mat[col][row] += cov;
816
817
818
819
820
821
                }
              }
            }
          }
        }
      }
822
823
824
825
826
827
828
829
830
831
    }
    box->releaseEvents();
  } while (MDiter.next());
  // normalise the covariance matrix
  cov_mat /= w_sum;                // normalise by sum of weights
  cov_mat.Diagonalise(Evec, Eval); // 3 x 3 matrices
  eigenvals = Eval.Diagonal();
  if (qAxisIsFixed) {
    // insert variance along Q (first eigenvector)
    eigenvals.insert(eigenvals.begin(), var_Qhat / w_sum);
832
    eigenvects.push_back(pos / pos.norm());
833
834
835
836
837
838
839
840
841
842
843
  }
  // set min eigenval to be small but non-zero (1e-6)
  // when no discernible peak above background
  std::replace_if(eigenvals.begin(), eigenvals.end(),
                  [&](auto x) { return x < 1e-6; }, 1e-6);
  // populate rest of eigenvect/vals
  for (size_t ivect = 0; ivect < cov_mat.numRows(); ++ivect) {
    if (!qAxisIsFixed) {
      eigenvects.push_back(V3D(Evec[0][ivect], Evec[1][ivect], Evec[2][ivect]));
    } else {
      // transform back to Qlab basis
844
      eigenvects.push_back(V3D(0.0, 0.0, 0.0));
845
      for (size_t ibasis = 0; ibasis < Evec.numRows(); ++ibasis) {
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
        eigenvects.back() +=
            V3D(Pinv[ibasis + 1][0], Pinv[ibasis + 1][1], Pinv[ibasis + 1][2]) *
            Evec[ibasis][ivect];
      }
    }
  }
}

/**
 * Get the inverse of the matrix P which transforms from Qlab to basis Qhat,
 * and uhat,vhat in plane perpendicular to Q. P is a matrix with columns
 * corresponding to new basis vectors. The inverse of P is equivilent to the
 * transpose (as for any roation matrix)
 *
 *  @param q     Qlab of peak center.
 *  @param Pinv  3 x 3 matrix with rows correpsonding to new basis vectors
 */
void IntegratePeaksMD2::getPinv(const V3D &q, Kernel::Matrix<double> &Pinv) {
  // loop over 3 mutually-orthogonal vectors  until get one with
  // a component perp to Q (within tolerance)
  double dotprod = 1;
  size_t ii = 0;
  V3D qhat = q / q.norm();
  V3D tmp;
  do {
    tmp = V3D(0, 0, 0); // reset u
    tmp[ii] = 1.0;
    dotprod = qhat.scalar_prod(tmp);
    ii++;
  } while (abs(dotprod) > 1.0 - 1e-6);
  // populate Pinv with basis vector rows
  Pinv.setRow(0, qhat);
  tmp = qhat.cross_prod(tmp);
  Pinv.setRow(1, tmp / tmp.norm());
  tmp = qhat.cross_prod(tmp);
  Pinv.setRow(2, tmp / tmp.norm());
}

884
/*
885
886
887
888
889
 * Define edges for each instrument by masking. For CORELLI, tubes 1 and
 *16, and pixels 0 and 255. Get Q in the lab frame for every peak, call it
 *C For every point on the edge, the trajectory in reciprocal space is a
 *straight line, going through O=V3D(0,0,0). Calculate a point at a fixed
 *momentum, say k=1. Q in the lab frame
890
891
892
 *E=V3D(-k*sin(tt)*cos(ph),-k*sin(tt)*sin(ph),k-k*cos(ph)).
 * Normalize E to 1: E=E*(1./E.norm())
 *
893
 * @param inst: instrument
894
 */
895
896
void IntegratePeaksMD2::calculateE1(
    const Geometry::DetectorInfo &detectorInfo) {
897
898
  for (size_t i = 0; i < detectorInfo.size(); ++i) {
    if (detectorInfo.isMonitor(i))
899
      continue; // skip monitor
900
    if (!detectorInfo.isMasked(i))
901
      continue; // edge is masked so don't check if not masked
902
903
904
    const auto &det = detectorInfo.detector(i);
    double tt1 = det.getTwoTheta(V3D(0, 0, 0), V3D(0, 0, 1)); // two theta
    double ph1 = det.getPhi();                                // phi
905
906
907
    V3D E1 = V3D(-std::sin(tt1) * std::cos(ph1), -std::sin(tt1) * std::sin(ph1),
                 1. - std::cos(tt1)); // end of trajectory
    E1 = E1 * (1. / E1.norm());       // normalize
908
    E1Vec.emplace_back(E1);
909
  }
910
}
911

912
913
/** Calculate if this Q is on a detector
 * The distance from C to OE is given by dv=C-E*(C.scalar_prod(E))
914
915
916
917
 * If dv.norm<integration_radius, one of the detector trajectories on the
 *edge is too close to the peak This method is applied to all masked
 *pixels. If there are masked pixels trajectories inside an integration
 *volume, the peak must be rejected.
918
919
920
921
 *
 * @param QLabFrame: The Peak center.
 * @param r: Peak radius.
 */
922
923
double IntegratePeaksMD2::detectorQ(Mantid::Kernel::V3D QLabFrame, double r) {
  double edge = r;
924
  for (auto &E1 : E1Vec) {
925
926
927
    V3D distv = QLabFrame - E1 * (QLabFrame.scalar_prod(E1)); // distance to the
                                                              // trajectory as a
                                                              // vector
928
    if (distv.norm() < r) {
929
      edge = distv.norm();
930
    }
931
  }
932
  return edge;
933
}
Hahn, Steven's avatar
Hahn, Steven committed
934

935
void IntegratePeaksMD2::runMaskDetectors(
David Fairbrother's avatar
David Fairbrother committed
936
937
    const Mantid::DataObjects::PeaksWorkspace_sptr &peakWS,
    const std::string &property, const std::string &values) {
Lynch, Vickie's avatar
Lynch, Vickie committed
938
  // For CORELLI do not count as edge if next to another detector bank
Lynch, Vickie's avatar
Lynch, Vickie committed
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
  if (property == "Tube" && peakWS->getInstrument()->getName() == "CORELLI") {
    IAlgorithm_sptr alg = createChildAlgorithm("MaskBTP");
    alg->setProperty<Workspace_sptr>("Workspace", peakWS);
    alg->setProperty("Bank", "1,7,12,17,22,27,30,59,63,69,74,79,84,89");
    alg->setProperty(property, "1");
    if (!alg->execute())
      throw std::runtime_error(
          "MaskDetectors Child Algorithm has not executed successfully");
    IAlgorithm_sptr alg2 = createChildAlgorithm("MaskBTP");
    alg2->setProperty<Workspace_sptr>("Workspace", peakWS);
    alg2->setProperty("Bank", "6,11,16,21,26,29,58,62,68,73,78,83,88,91");
    alg2->setProperty(property, "16");
    if (!alg2->execute())
      throw std::runtime_error(
          "MaskDetectors Child Algorithm has not executed successfully");
  } else {
    IAlgorithm_sptr alg = createChildAlgorithm("MaskBTP");
    alg->setProperty<Workspace_sptr>("Workspace", peakWS);
    alg->setProperty(property, values);
    if (!alg->execute())
      throw std::runtime_error(
          "MaskDetectors Child Algorithm has not executed successfully");
961
  }
962
963
}

964
void IntegratePeaksMD2::checkOverlap(
David Fairbrother's avatar
David Fairbrother committed
965
    int i, const Mantid::DataObjects::PeaksWorkspace_sptr &peakWS,
966
    Mantid::Kernel::SpecialCoordinateSystem CoordinatesToUse, double radius) {
967
968
969
  // Get a direct ref to that peak.
  IPeak &p1 = peakWS->getPeak(i);
  V3D pos1;
970
  if (CoordinatesToUse == Kernel::QLab) //"Q (lab frame)"
971
    pos1 = p1.getQLabFrame();
972
  else if (CoordinatesToUse == Kernel::QSample) //"Q (sample frame)"
973
    pos1 = p1.getQSampleFrame();
974
  else if (CoordinatesToUse == Kernel::HKL) //"HKL"
975
976
977
978
979
    pos1 = p1.getHKL();
  for (int j = i + 1; j < peakWS->getNumberPeaks(); ++j) {
    // Get a direct ref to rest of peaks peak.
    IPeak &p2 = peakWS->getPeak(j);
    V3D pos2;
980
    if (CoordinatesToUse == Kernel::QLab) //"Q (lab frame)"
981
      pos2 = p2.getQLabFrame();
982
    else if (CoordinatesToUse == Kernel::QSample) //"Q (sample frame)"
983
      pos2 = p2.getQSampleFrame();
984
    else if (CoordinatesToUse == Kernel::HKL) //"HKL"
985
986
987
988
      pos2 = p2.getHKL();
    if (pos1.distance(pos2) < radius) {
      g_log.warning() << " Warning:  Peak integration spheres for peaks " << i
                      << " and " << j << " overlap.  Distance between peaks is "
989
                      << pos1.distance(pos2) << '\n';
990
    }
Lynch, Vickie's avatar
Lynch, Vickie committed
991
  }
992
}
993

994
995
996
997
998
999
1000
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void IntegratePeaksMD2::exec() {
  inWS = getProperty("InputWorkspace");

  CALL_MDEVENT_FUNCTION(this->integrate, inWS);
For faster browsing, not all history is shown. View entire blame