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
158
159
160
  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());
  std::remove_copy_if(gen.begin(), gen.end(),
                      std::back_inserter(possibleHKLs),
                      (~dSpacingFilter)->fn());
161

162
163
  size_t N = possibleHKLs.size();
  N = max<size_t>(100, N);
164
165
  auto &UB = lattice.getUB();
  goniometer = peak0.getGoniometerMatrix();
166
  Progress prog(this, 0.0, 1.0, N);
167
  vector<vector<int>> AlreadyDonePeaks;
168
  auto orientedUB = goniometer * UB;
169
  HKLFilterWavelength lambdaFilter(orientedUB, lambdaMin, lambdaMax);
Lynch, Vickie's avatar
Lynch, Vickie committed
170
  int seqNum = 0;
Lynch, Vickie's avatar
Lynch, Vickie committed
171
172
173
174
175
176
177
  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
178
    V3D hkl = *it;
Lynch, Vickie's avatar
Lynch, Vickie committed
179
180
181
182
183
184
185
186
187
188
    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);
  }
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
224
  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
225
226
227
228
229
230
231
232
  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
233
234
    auto peak = *it;
    V3D hkl = peak.getHKL();
Lynch, Vickie's avatar
Lynch, Vickie committed
235
236
237
238
239
240
241
242
243
    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);
244
245
246
247
  }
  setProperty("SatellitePeaks", OutPeaks);
}

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

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

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

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

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

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

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

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

    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
302
  }
303
304
}

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

316
317
} // namespace Crystal
} // namespace Mantid