PredictSatellitePeaks.cpp 12.2 KB
Newer Older
1
2
3
/*
* PredictSatellitePeaks.cpp
*
4
5
*  Created on: July 15, 2018
*      Author: Vickie Lynch
6
7
8
9
10
11
12
13
14
15
*/
#include "MantidCrystal/PredictSatellitePeaks.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Objects/InstrumentRayTracer.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/ArrayLengthValidator.h"
#include "MantidKernel/EnabledWhenProperty.h"
16
#include "MantidAPI/Run.h"
17
18
#include "MantidGeometry/Crystal/BasicHKLFilters.h"
#include "MantidGeometry/Crystal/HKLGenerator.h"
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 <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

  declareProperty(
Lynch, Vickie's avatar
Lynch, Vickie committed
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
      "IncludeIntegerHKL", true,
      "If false order 0 peaks are not included in workspace (integer HKL)");

  declareProperty("IncludeAllPeaksInRange", false, "If false only offsets from "
                                                   "peaks from Peaks workspace "
                                                   "in input are used");

  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");
78
79

  setPropertySettings(
Lynch, Vickie's avatar
Lynch, Vickie committed
80
81
82
      "WavelengthMin",
      Kernel::make_unique<Kernel::EnabledWhenProperty>(
          string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
83
84

  setPropertySettings(
Lynch, Vickie's avatar
Lynch, Vickie committed
85
86
87
      "WavelengthMax",
      Kernel::make_unique<Kernel::EnabledWhenProperty>(
          string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
88
  setPropertySettings(
Lynch, Vickie's avatar
Lynch, Vickie committed
89
90
91
      "MinDSpacing",
      Kernel::make_unique<Kernel::EnabledWhenProperty>(
          string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
92
93

  setPropertySettings(
Lynch, Vickie's avatar
Lynch, Vickie committed
94
95
96
      "MaxDSpacing",
      Kernel::make_unique<Kernel::EnabledWhenProperty>(
          string("IncludeAllPeaksInRange"), Kernel::IS_EQUAL_TO, "1"));
97
98
99
100
}

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

111
112
113
  V3D offsets1 = getOffsetVector("ModVector1");
  V3D offsets2 = getOffsetVector("ModVector2");
  V3D offsets3 = getOffsetVector("ModVector3");
Lynch, Vickie's avatar
Lynch, Vickie committed
114
  int maxOrder = getProperty("MaxOrder");
115

116
  bool includeOrderZero = getProperty("IncludeIntegerHKL");
117
118
119
120
121
122

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

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

125
  OrientedLattice lattice = sample.getOrientedLattice();
126

127
  const auto instrument = Peaks->getInstrument();
128
129
130

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

134
135
  Kernel::Matrix<double> goniometer;
  goniometer.identityMatrix();
136

137
138
  const double lambdaMin = getProperty("WavelengthMin");
  const double lambdaMax = getProperty("WavelengthMax");
139
140
  IPeak &peak0 = Peaks->getPeak(0);

141
  std::vector<V3D> possibleHKLs;
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
  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
158
  std::remove_copy_if(gen.begin(), gen.end(), std::back_inserter(possibleHKLs),
159
                      (~dSpacingFilter)->fn());
160

161
162
  size_t N = possibleHKLs.size();
  N = max<size_t>(100, N);
163
164
  auto &UB = lattice.getUB();
  goniometer = peak0.getGoniometerMatrix();
165
  Progress prog(this, 0.0, 1.0, N);
166
  vector<vector<int>> AlreadyDonePeaks;
167
  auto orientedUB = goniometer * UB;
168
  HKLFilterWavelength lambdaFilter(orientedUB, lambdaMin, lambdaMax);
Lynch, Vickie's avatar
Lynch, Vickie committed
169
  int seqNum = 0;
Lynch, Vickie's avatar
Lynch, Vickie committed
170
171
172
173
174
175
176
  OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset1", offsets1,
                                                          true);
  OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset2", offsets2,
                                                          true);
  OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset3", offsets3,
                                                          true);
  for (auto it = possibleHKLs.begin(); it != possibleHKLs.end(); ++it) {
Lynch, Vickie's avatar
Lynch, Vickie committed
177
    V3D hkl = *it;
Lynch, Vickie's avatar
Lynch, Vickie committed
178
179
180
181
182
183
184
185
186
187
    predictOffsets(Peaks, OutPeaks, offsets1, maxOrder, hkl, lambdaFilter,
                   includePeaksInRange, includeOrderZero, seqNum,
                   AlreadyDonePeaks);
    predictOffsets(Peaks, OutPeaks, offsets2, maxOrder, hkl, lambdaFilter,
                   includePeaksInRange, includeOrderZero, seqNum,
                   AlreadyDonePeaks);
    predictOffsets(Peaks, OutPeaks, offsets3, maxOrder, hkl, lambdaFilter,
                   includePeaksInRange, includeOrderZero, seqNum,
                   AlreadyDonePeaks);
  }
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
  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");

  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
224
225
226
227
228
229
230
231
  OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset1", offsets1,
                                                          true);
  OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset2", offsets2,
                                                          true);
  OutPeaks->mutableRun().addProperty<std::vector<double>>("Offset3", offsets3,
                                                          true);
  std::vector<Peak> peaks = Peaks->getPeaks();
  for (auto it = peaks.begin(); it != peaks.end(); ++it) {
Lynch, Vickie's avatar
Lynch, Vickie committed
232
233
    auto peak = *it;
    V3D hkl = peak.getHKL();
Lynch, Vickie's avatar
Lynch, Vickie committed
234
235
236
237
238
239
240
241
242
    predictOffsets(Peaks, OutPeaks, offsets1, maxOrder, hkl, lambdaFilter,
                   includePeaksInRange, includeOrderZero, seqNum,
                   AlreadyDonePeaks);
    predictOffsets(Peaks, OutPeaks, offsets2, maxOrder, hkl, lambdaFilter,
                   includePeaksInRange, includeOrderZero, seqNum,
                   AlreadyDonePeaks);
    predictOffsets(Peaks, OutPeaks, offsets3, maxOrder, hkl, lambdaFilter,
                   includePeaksInRange, includeOrderZero, seqNum,
                   AlreadyDonePeaks);
243
244
245
246
  }
  setProperty("SatellitePeaks", OutPeaks);
}

Lynch, Vickie's avatar
Lynch, Vickie committed
247
248
void PredictSatellitePeaks::predictOffsets(
    DataObjects::PeaksWorkspace_sptr Peaks,
Lynch, Vickie's avatar
Lynch, Vickie committed
249
250
251
    boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks, V3D offsets,
    int &maxOrder, V3D &hkl, HKLFilterWavelength &lambdaFilter,
    bool &includePeaksInRange, bool &includeOrderZero, int &seqNum,
Lynch, Vickie's avatar
Lynch, Vickie committed
252
    vector<vector<int>> &AlreadyDonePeaks) {
Lynch, Vickie's avatar
Lynch, Vickie committed
253
254
  if (offsets == V3D(0, 0, 0))
    return;
255
  const Kernel::DblMatrix &UB = Peaks->sample().getOrientedLattice().getUB();
Lynch, Vickie's avatar
Lynch, Vickie committed
256
  IPeak &peak1 = Peaks->getPeak(0);
257
258
  Kernel::Matrix<double> goniometer = peak1.getGoniometerMatrix();
  auto RunNumber = peak1.getRunNumber();
259
  Geometry::InstrumentRayTracer tracer(Peaks->getInstrument());
Lynch, Vickie's avatar
Lynch, Vickie committed
260
261
262
  for (int order = -maxOrder; order <= maxOrder; order++) {
    if (order == 0 && !includeOrderZero)
      continue; // exclude order 0
263
    V3D satelliteHKL(hkl);
Lynch, Vickie's avatar
Lynch, Vickie committed
264

265
266
    satelliteHKL += offsets * order;
    if (!lambdaFilter.isAllowed(satelliteHKL) && includePeaksInRange)
Lynch, Vickie's avatar
Lynch, Vickie committed
267
      continue;
268

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

271
    // Check if Q is non-physical
Lynch, Vickie's avatar
Lynch, Vickie committed
272
273
    if (Qs[2] <= 0)
      continue;
Lynch, Vickie's avatar
Lynch, Vickie committed
274

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

277
    peak->setGoniometerMatrix(goniometer);
Lynch, Vickie's avatar
Lynch, Vickie committed
278

Lynch, Vickie's avatar
Lynch, Vickie committed
279
280
    if (!peak->findDetector(tracer))
      continue;
Lynch, Vickie's avatar
Lynch, Vickie committed
281
    vector<int> SavPk{RunNumber, boost::math::iround(1000.0 * satelliteHKL[0]),
Lynch, Vickie's avatar
Lynch, Vickie committed
282
283
284
285
286
                      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
287

288
    if (!foundPeak) {
Lynch, Vickie's avatar
Lynch, Vickie committed
289
290
291
      AlreadyDonePeaks.push_back(SavPk);
    } else {
      continue;
Lynch, Vickie's avatar
Lynch, Vickie committed
292
    }
Lynch, Vickie's avatar
Lynch, Vickie committed
293
294
295
296
297
298
299
300

    peak->setHKL(satelliteHKL);
    peak->setIntHKL(hkl);
    peak->setPeakNumber(seqNum);
    seqNum++;
    peak->setRunNumber(RunNumber);
    peak->setIntMNP(V3D(order, 0, 0));
    OutPeaks->addPeak(*peak);
Lynch, Vickie's avatar
Lynch, Vickie committed
301
  }
302
303
}

304
305
306
307
308
309
310
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
311
  V3D offsets1 = V3D(offsets[0], offsets[1], offsets[2]);
312
313
314
  return offsets1;
}

315
316
} // namespace Crystal
} // namespace Mantid