PredictSatellitePeaks.cpp 15.3 KB
Newer Older
1
/*
Lynch, Vickie's avatar
Lynch, Vickie committed
2
3
4
5
6
 * PredictSatellitePeaks.cpp
 *
 *  Created on: July 15, 2018
 *      Author: Vickie Lynch
 */
7
#include "MantidCrystal/PredictSatellitePeaks.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
8
#include "MantidAPI/Run.h"
9
10
11
#include "MantidAPI/Sample.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/PeaksWorkspace.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
12
13
#include "MantidGeometry/Crystal/BasicHKLFilters.h"
#include "MantidGeometry/Crystal/HKLGenerator.h"
14
15
16
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Objects/InstrumentRayTracer.h"
#include "MantidKernel/ArrayLengthValidator.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
17
#include "MantidKernel/ArrayProperty.h"
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include "MantidKernel/EnabledWhenProperty.h"
#include <boost/math/special_functions/round.hpp>

namespace Mantid {
using namespace Mantid::DataObjects;
using namespace Mantid::API;
using namespace std;
using namespace Mantid::Geometry;
using namespace Mantid::Kernel;
namespace Crystal {

DECLARE_ALGORITHM(PredictSatellitePeaks)

/// Initialise the properties
void PredictSatellitePeaks::init() {
  declareProperty(
      make_unique<WorkspaceProperty<PeaksWorkspace>>("Peaks", "",
                                                     Direction::Input),
      "Workspace of Peaks with orientation matrix that indexed the peaks and "
      "instrument loaded");

  declareProperty(
      make_unique<WorkspaceProperty<PeaksWorkspace>>("SatellitePeaks", "",
                                                     Direction::Output),
      "Workspace of Peaks with peaks with fractional h,k, and/or l values");
  declareProperty(Kernel::make_unique<Kernel::ArrayProperty<double>>(
Lynch, Vickie's avatar
Lynch, Vickie committed
44
45
                      string("ModVector1"), "0.0,0.0,0.0"),
                  "Offsets for h, k, l directions ");
46
  declareProperty(Kernel::make_unique<Kernel::ArrayProperty<double>>(
Lynch, Vickie's avatar
Lynch, Vickie committed
47
48
                      string("ModVector2"), "0.0,0.0,0.0"),
                  "Offsets for h, k, l directions ");
49
  declareProperty(Kernel::make_unique<Kernel::ArrayProperty<double>>(
Lynch, Vickie's avatar
Lynch, Vickie committed
50
51
                      string("ModVector3"), "0.0,0.0,0.0"),
                  "Offsets for h, k, l directions ");
Lynch, Vickie's avatar
Lynch, Vickie committed
52
53
54
  declareProperty(
      make_unique<PropertyWithValue<int>>("MaxOrder", 0, Direction::Input),
      "Maximum order to apply ModVectors. Default = 0");
55

56
57
58
59
  declareProperty(make_unique<PropertyWithValue<bool>>("CrossTerms", false,
                                                       Direction::Input),
                  "Include cross terms (false)");

60
  declareProperty(
Lynch, Vickie's avatar
Lynch, Vickie committed
61
62
63
      "IncludeIntegerHKL", true,
      "If false order 0 peaks are not included in workspace (integer HKL)");

Lynch, Vickie's avatar
Lynch, Vickie committed
64
65
66
67
  declareProperty("IncludeAllPeaksInRange", false,
                  "If false only offsets from "
                  "peaks from Peaks workspace "
                  "in input are used");
Lynch, Vickie's avatar
Lynch, Vickie committed
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

  declareProperty(make_unique<PropertyWithValue<double>>("WavelengthMin", 0.1,
                                                         Direction::Input),
                  "Minimum wavelength limit at which to start looking for "
                  "single-crystal peaks.");
  declareProperty(make_unique<PropertyWithValue<double>>("WavelengthMax", 100.0,
                                                         Direction::Input),
                  "Maximum wavelength limit at which to start looking for "
                  "single-crystal peaks.");
  declareProperty(make_unique<PropertyWithValue<double>>("MinDSpacing", 0.1,
                                                         Direction::Input),
                  "Minimum d-spacing of peaks to consider. Default = 1.0");
  declareProperty(make_unique<PropertyWithValue<double>>("MaxDSpacing", 100.0,
                                                         Direction::Input),
                  "Maximum d-spacing of peaks to consider");
83
84

  setPropertySettings(
Lynch, Vickie's avatar
Lynch, Vickie committed
85
86
87
      "WavelengthMin",
      Kernel::make_unique<Kernel::EnabledWhenProperty>(
          string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
88
89

  setPropertySettings(
Lynch, Vickie's avatar
Lynch, Vickie committed
90
91
92
      "WavelengthMax",
      Kernel::make_unique<Kernel::EnabledWhenProperty>(
          string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
93
  setPropertySettings(
Lynch, Vickie's avatar
Lynch, Vickie committed
94
95
96
      "MinDSpacing",
      Kernel::make_unique<Kernel::EnabledWhenProperty>(
          string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
97
98

  setPropertySettings(
Lynch, Vickie's avatar
Lynch, Vickie committed
99
100
101
      "MaxDSpacing",
      Kernel::make_unique<Kernel::EnabledWhenProperty>(
          string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
102
103
104
105
}

/// Run the algorithm
void PredictSatellitePeaks::exec() {
106
107
108
109
110
  bool includePeaksInRange = getProperty("IncludeAllPeaksInRange");
  if (!includePeaksInRange) {
    exec_peaks();
    return;
  }
111
112
113
114
115
  PeaksWorkspace_sptr Peaks = getProperty("Peaks");
  if (!Peaks)
    throw std::invalid_argument(
        "Input workspace is not a PeaksWorkspace. Type=" + Peaks->id());

116
117
118
  V3D offsets1 = getOffsetVector("ModVector1");
  V3D offsets2 = getOffsetVector("ModVector2");
  V3D offsets3 = getOffsetVector("ModVector3");
Lynch, Vickie's avatar
Lynch, Vickie committed
119
  int maxOrder = getProperty("MaxOrder");
120
  bool crossTerms = getProperty("CrossTerms");
121

122
  bool includeOrderZero = getProperty("IncludeIntegerHKL");
123
124
125
126
127
128

  if (Peaks->getNumberPeaks() <= 0) {
    g_log.error() << "There are No peaks in the input PeaksWorkspace\n";
    return;
  }

129
  API::Sample sample = Peaks->mutableSample();
130

131
  OrientedLattice lattice = sample.getOrientedLattice();
132

133
  const auto instrument = Peaks->getInstrument();
134
135
136

  auto OutPeaks = boost::dynamic_pointer_cast<IPeaksWorkspace>(
      WorkspaceFactory::Instance().createPeaks());
137
  OutPeaks->setInstrument(instrument);
138
  OutPeaks->mutableSample().setOrientedLattice(&lattice);
139

140
141
  Kernel::Matrix<double> goniometer;
  goniometer.identityMatrix();
142

143
144
  const double lambdaMin = getProperty("WavelengthMin");
  const double lambdaMax = getProperty("WavelengthMax");
145
146
  IPeak &peak0 = Peaks->getPeak(0);

147
  std::vector<V3D> possibleHKLs;
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
  const double dMin = getProperty("MinDSpacing");
  const double dMax = getProperty("MaxDSpacing");
  Geometry::HKLGenerator gen(lattice, dMin);
  auto dSpacingFilter =
      boost::make_shared<HKLFilterDRange>(lattice, dMin, dMax);

  V3D hkl = *(gen.begin());
  g_log.information() << "HKL range for d_min of " << dMin << " to d_max of "
                      << dMax << " is from " << hkl << " to " << hkl * -1.0
                      << ", a total of " << gen.size() << " possible HKL's\n";
  if (gen.size() > MAX_NUMBER_HKLS)
    throw std::invalid_argument("More than 10 billion HKLs to search. Is "
                                "your d_min value too small?");

  possibleHKLs.clear();
  possibleHKLs.reserve(gen.size());
Lynch, Vickie's avatar
Lynch, Vickie committed
164
  std::remove_copy_if(gen.begin(), gen.end(), std::back_inserter(possibleHKLs),
165
                      (~dSpacingFilter)->fn());
166

167
168
  size_t N = possibleHKLs.size();
  N = max<size_t>(100, N);
169
170
  auto &UB = lattice.getUB();
  goniometer = peak0.getGoniometerMatrix();
171
  Progress prog(this, 0.0, 1.0, N);
172
  vector<vector<int>> AlreadyDonePeaks;
173
  auto orientedUB = goniometer * UB;
174
  HKLFilterWavelength lambdaFilter(orientedUB, lambdaMin, lambdaMax);
Lynch, Vickie's avatar
Lynch, Vickie committed
175
  int seqNum = 0;
Lynch, Vickie's avatar
Lynch, Vickie committed
176
177
178
179
  Sample &sampleOut = OutPeaks->mutableSample();
  sampleOut.getOrientedLattice().setModVec1(offsets1);
  sampleOut.getOrientedLattice().setModVec2(offsets2);
  sampleOut.getOrientedLattice().setModVec3(offsets3);
Lynch, Vickie's avatar
Lynch, Vickie committed
180
  for (auto it = possibleHKLs.begin(); it != possibleHKLs.end(); ++it) {
Lynch, Vickie's avatar
Lynch, Vickie committed
181
    V3D hkl = *it;
182
    if (crossTerms) {
Lynch, Vickie's avatar
Lynch, Vickie committed
183
184
185
186
      predictOffsetsWithCrossTerms(Peaks, OutPeaks, offsets1, offsets2,
                                   offsets3, maxOrder, hkl, lambdaFilter,
                                   includePeaksInRange, includeOrderZero,
                                   seqNum, AlreadyDonePeaks);
187
188
189
190
191
192
193
194
195
196
197
    } else {
      predictOffsets(Peaks, OutPeaks, 0, offsets1, maxOrder, hkl, lambdaFilter,
                     includePeaksInRange, includeOrderZero, seqNum,
                     AlreadyDonePeaks);
      predictOffsets(Peaks, OutPeaks, 1, offsets2, maxOrder, hkl, lambdaFilter,
                     includePeaksInRange, includeOrderZero, seqNum,
                     AlreadyDonePeaks);
      predictOffsets(Peaks, OutPeaks, 2, offsets3, maxOrder, hkl, lambdaFilter,
                     includePeaksInRange, includeOrderZero, seqNum,
                     AlreadyDonePeaks);
    }
Lynch, Vickie's avatar
Lynch, Vickie committed
198
  }
199
200
201
202
203
204
205
206
207
208
209
210
211
  setProperty("SatellitePeaks", OutPeaks);
}

void PredictSatellitePeaks::exec_peaks() {
  PeaksWorkspace_sptr Peaks = getProperty("Peaks");
  if (!Peaks)
    throw std::invalid_argument(
        "Input workspace is not a PeaksWorkspace. Type=" + Peaks->id());

  V3D offsets1 = getOffsetVector("ModVector1");
  V3D offsets2 = getOffsetVector("ModVector2");
  V3D offsets3 = getOffsetVector("ModVector3");
  int maxOrder = getProperty("MaxOrder");
212
  bool crossTerms = getProperty("CrossTerms");
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235

  bool includePeaksInRange = false;
  bool includeOrderZero = getProperty("IncludeIntegerHKL");

  if (Peaks->getNumberPeaks() <= 0) {
    g_log.error() << "There are No peaks in the input PeaksWorkspace\n";
    return;
  }

  API::Sample sample = Peaks->mutableSample();

  OrientedLattice lattice = sample.getOrientedLattice();

  const auto instrument = Peaks->getInstrument();

  auto OutPeaks = boost::dynamic_pointer_cast<IPeaksWorkspace>(
      WorkspaceFactory::Instance().createPeaks());
  OutPeaks->setInstrument(instrument);
  OutPeaks->mutableSample().setOrientedLattice(&lattice);

  vector<vector<int>> AlreadyDonePeaks;
  int seqNum = 0;
  HKLFilterWavelength lambdaFilter(DblMatrix(3, 3, true), 0.1, 100.);
Lynch, Vickie's avatar
Lynch, Vickie committed
236
237
238
239
  Sample &sampleOut = OutPeaks->mutableSample();
  sampleOut.getOrientedLattice().setModVec1(offsets1);
  sampleOut.getOrientedLattice().setModVec2(offsets2);
  sampleOut.getOrientedLattice().setModVec3(offsets3);
Lynch, Vickie's avatar
Lynch, Vickie committed
240
241
  std::vector<Peak> peaks = Peaks->getPeaks();
  for (auto it = peaks.begin(); it != peaks.end(); ++it) {
Lynch, Vickie's avatar
Lynch, Vickie committed
242
243
    auto peak = *it;
    V3D hkl = peak.getHKL();
244
    if (crossTerms) {
Lynch, Vickie's avatar
Lynch, Vickie committed
245
246
247
248
      predictOffsetsWithCrossTerms(Peaks, OutPeaks, offsets1, offsets2,
                                   offsets3, maxOrder, hkl, lambdaFilter,
                                   includePeaksInRange, includeOrderZero,
                                   seqNum, AlreadyDonePeaks);
249
250
251
252
253
254
255
256
257
258
259
    } else {
      predictOffsets(Peaks, OutPeaks, 0, offsets1, maxOrder, hkl, lambdaFilter,
                     includePeaksInRange, includeOrderZero, seqNum,
                     AlreadyDonePeaks);
      predictOffsets(Peaks, OutPeaks, 1, offsets2, maxOrder, hkl, lambdaFilter,
                     includePeaksInRange, includeOrderZero, seqNum,
                     AlreadyDonePeaks);
      predictOffsets(Peaks, OutPeaks, 2, offsets3, maxOrder, hkl, lambdaFilter,
                     includePeaksInRange, includeOrderZero, seqNum,
                     AlreadyDonePeaks);
    }
260
261
262
263
  }
  setProperty("SatellitePeaks", OutPeaks);
}

Lynch, Vickie's avatar
Lynch, Vickie committed
264
265
void PredictSatellitePeaks::predictOffsets(
    DataObjects::PeaksWorkspace_sptr Peaks,
Lynch, Vickie's avatar
Lynch, Vickie committed
266
267
    boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks, int iVector,
    V3D offsets, int &maxOrder, V3D &hkl, HKLFilterWavelength &lambdaFilter,
Lynch, Vickie's avatar
Lynch, Vickie committed
268
    bool &includePeaksInRange, bool &includeOrderZero, int &seqNum,
Lynch, Vickie's avatar
Lynch, Vickie committed
269
    vector<vector<int>> &AlreadyDonePeaks) {
Lynch, Vickie's avatar
Lynch, Vickie committed
270
271
  if (offsets == V3D(0, 0, 0))
    return;
272
  const Kernel::DblMatrix &UB = Peaks->sample().getOrientedLattice().getUB();
Lynch, Vickie's avatar
Lynch, Vickie committed
273
  IPeak &peak1 = Peaks->getPeak(0);
274
275
  Kernel::Matrix<double> goniometer = peak1.getGoniometerMatrix();
  auto RunNumber = peak1.getRunNumber();
276
  Geometry::InstrumentRayTracer tracer(Peaks->getInstrument());
Lynch, Vickie's avatar
Lynch, Vickie committed
277
278
279
  for (int order = -maxOrder; order <= maxOrder; order++) {
    if (order == 0 && !includeOrderZero)
      continue; // exclude order 0
280
    V3D satelliteHKL(hkl);
Lynch, Vickie's avatar
Lynch, Vickie committed
281

282
283
    satelliteHKL += offsets * order;
    if (!lambdaFilter.isAllowed(satelliteHKL) && includePeaksInRange)
Lynch, Vickie's avatar
Lynch, Vickie committed
284
      continue;
285

286
    Kernel::V3D Qs = goniometer * UB * satelliteHKL * 2.0 * M_PI;
Lynch, Vickie's avatar
Lynch, Vickie committed
287

288
    // Check if Q is non-physical
Lynch, Vickie's avatar
Lynch, Vickie committed
289
290
    if (Qs[2] <= 0)
      continue;
Lynch, Vickie's avatar
Lynch, Vickie committed
291

292
    auto peak(Peaks->createPeak(Qs, 1));
Lynch, Vickie's avatar
Lynch, Vickie committed
293

294
    peak->setGoniometerMatrix(goniometer);
Lynch, Vickie's avatar
Lynch, Vickie committed
295

Lynch, Vickie's avatar
Lynch, Vickie committed
296
297
    if (!peak->findDetector(tracer))
      continue;
Lynch, Vickie's avatar
Lynch, Vickie committed
298
    vector<int> SavPk{RunNumber, boost::math::iround(1000.0 * satelliteHKL[0]),
Lynch, Vickie's avatar
Lynch, Vickie committed
299
300
301
302
303
                      boost::math::iround(1000.0 * satelliteHKL[1]),
                      boost::math::iround(1000.0 * satelliteHKL[2])};

    bool foundPeak =
        binary_search(AlreadyDonePeaks.begin(), AlreadyDonePeaks.end(), SavPk);
Lynch, Vickie's avatar
Lynch, Vickie committed
304

305
    if (!foundPeak) {
Lynch, Vickie's avatar
Lynch, Vickie committed
306
307
308
      AlreadyDonePeaks.push_back(SavPk);
    } else {
      continue;
Lynch, Vickie's avatar
Lynch, Vickie committed
309
    }
Lynch, Vickie's avatar
Lynch, Vickie committed
310
311
312
313
314
315

    peak->setHKL(satelliteHKL);
    peak->setIntHKL(hkl);
    peak->setPeakNumber(seqNum);
    seqNum++;
    peak->setRunNumber(RunNumber);
316
317
318
319
320
321
322
323
324
325
    V3D mnp;
    mnp[iVector] = order;
    peak->setIntMNP(mnp);
    OutPeaks->addPeak(*peak);
  }
}

void PredictSatellitePeaks::predictOffsetsWithCrossTerms(
    DataObjects::PeaksWorkspace_sptr Peaks,
    boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks, V3D offsets1,
Lynch, Vickie's avatar
Lynch, Vickie committed
326
327
328
    V3D offsets2, V3D offsets3, int &maxOrder, V3D &hkl,
    HKLFilterWavelength &lambdaFilter, bool &includePeaksInRange,
    bool &includeOrderZero, int &seqNum,
329
    vector<vector<int>> &AlreadyDonePeaks) {
Lynch, Vickie's avatar
Lynch, Vickie committed
330
331
  if (offsets1 == V3D(0, 0, 0) && offsets2 == V3D(0, 0, 0) &&
      offsets3 == V3D(0, 0, 0))
332
333
334
335
336
337
338
    return;
  const Kernel::DblMatrix &UB = Peaks->sample().getOrientedLattice().getUB();
  IPeak &peak1 = Peaks->getPeak(0);
  Kernel::Matrix<double> goniometer = peak1.getGoniometerMatrix();
  auto RunNumber = peak1.getRunNumber();
  Geometry::InstrumentRayTracer tracer(Peaks->getInstrument());
  DblMatrix offsetsMat(3, 3);
339
340
341
  offsetsMat.setColumn(0, offsets1);
  offsetsMat.setColumn(1, offsets2);
  offsetsMat.setColumn(2, offsets3);
Lynch, Vickie's avatar
Lynch, Vickie committed
342
  int maxOrder1 = maxOrder;
Lynch, Vickie's avatar
Lynch, Vickie committed
343
344
  if (offsets1 == V3D(0, 0, 0))
    maxOrder1 = 0;
Lynch, Vickie's avatar
Lynch, Vickie committed
345
  int maxOrder2 = maxOrder;
Lynch, Vickie's avatar
Lynch, Vickie committed
346
347
  if (offsets2 == V3D(0, 0, 0))
    maxOrder2 = 0;
Lynch, Vickie's avatar
Lynch, Vickie committed
348
  int maxOrder3 = maxOrder;
Lynch, Vickie's avatar
Lynch, Vickie committed
349
350
  if (offsets3 == V3D(0, 0, 0))
    maxOrder3 = 0;
Lynch, Vickie's avatar
Lynch, Vickie committed
351
352
353
  for (int m = -maxOrder1; m <= maxOrder1; m++)
    for (int n = -maxOrder2; n <= maxOrder2; n++)
      for (int p = -maxOrder3; p <= maxOrder3; p++) {
354
355
356
357
358
359
360
361
        if (m == 0 && n == 0 && p == 0 && !includeOrderZero)
          continue; // exclude 0,0,0
        V3D satelliteHKL(hkl);
        V3D mnp = V3D(m, n, p);
        satelliteHKL -= offsetsMat * mnp;
        if (!lambdaFilter.isAllowed(satelliteHKL) && includePeaksInRange)
          continue;

Lynch, Vickie's avatar
Lynch, Vickie committed
362
        Kernel::V3D Qs = goniometer * UB * satelliteHKL * 2.0 * M_PI;
363

Lynch, Vickie's avatar
Lynch, Vickie committed
364
365
366
        // Check if Q is non-physical
        if (Qs[2] <= 0)
          continue;
367

Lynch, Vickie's avatar
Lynch, Vickie committed
368
        auto peak(Peaks->createPeak(Qs, 1));
369

Lynch, Vickie's avatar
Lynch, Vickie committed
370
        peak->setGoniometerMatrix(goniometer);
371

Lynch, Vickie's avatar
Lynch, Vickie committed
372
373
374
375
376
377
        if (!peak->findDetector(tracer))
          continue;
        vector<int> SavPk{RunNumber,
                          boost::math::iround(1000.0 * satelliteHKL[0]),
                          boost::math::iround(1000.0 * satelliteHKL[1]),
                          boost::math::iround(1000.0 * satelliteHKL[2])};
378

Lynch, Vickie's avatar
Lynch, Vickie committed
379
380
        bool foundPeak = binary_search(AlreadyDonePeaks.begin(),
                                       AlreadyDonePeaks.end(), SavPk);
381

Lynch, Vickie's avatar
Lynch, Vickie committed
382
383
384
385
386
387
388
389
390
391
392
393
394
395
        if (!foundPeak) {
          AlreadyDonePeaks.push_back(SavPk);
        } else {
          continue;
        }

        peak->setHKL(satelliteHKL);
        peak->setIntHKL(hkl);
        peak->setPeakNumber(seqNum);
        seqNum++;
        peak->setRunNumber(RunNumber);
        peak->setIntMNP(V3D(m, n, p));
        OutPeaks->addPeak(*peak);
      }
396
397
}

398
399
400
401
402
403
404
V3D PredictSatellitePeaks::getOffsetVector(const std::string &label) {
  vector<double> offsets = getProperty(label);
  if (offsets.empty()) {
    offsets.push_back(0.0);
    offsets.push_back(0.0);
    offsets.push_back(0.0);
  }
Lynch, Vickie's avatar
Lynch, Vickie committed
405
  V3D offsets1 = V3D(offsets[0], offsets[1], offsets[2]);
406
407
408
  return offsets1;
}

409
410
} // namespace Crystal
} // namespace Mantid