PredictSatellitePeaks.cpp 10.7 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
101
102
103
104
105
}

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

Lynch, Vickie's avatar
Lynch, Vickie committed
106
107
108
109
  vector<double> offsets1 = getProperty("ModVector1");
  vector<double> offsets2 = getProperty("ModVector2");
  vector<double> offsets3 = getProperty("ModVector3");
  int maxOrder = getProperty("MaxOrder");
110
111
112
113
114
115
116
117
118
119
120
121
122
123
  if (offsets1.empty()) {
    offsets1.push_back(0.0);
    offsets1.push_back(0.0);
    offsets1.push_back(0.0);
  }
  if (offsets2.empty()) {
    offsets2.push_back(0.0);
    offsets2.push_back(0.0);
    offsets2.push_back(0.0);
  }
  if (offsets3.empty()) {
    offsets3.push_back(0.0);
    offsets3.push_back(0.0);
    offsets3.push_back(0.0);
Lynch, Vickie's avatar
Lynch, Vickie committed
124
  }
125
126

  bool includePeaksInRange = getProperty("IncludeAllPeaksInRange");
127
  bool includeOrderZero = getProperty("IncludeIntegerHKL");
128
129
130
131
132
133

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

134
  API::Sample sample = Peaks->sample();
135

136
  const auto &lattice = sample.getOrientedLattice();
137

138
  const auto instrument = Peaks->getInstrument();
139
140
141

  auto OutPeaks = boost::dynamic_pointer_cast<IPeaksWorkspace>(
      WorkspaceFactory::Instance().createPeaks());
142
  OutPeaks->setInstrument(instrument);
143
144
145
146

  V3D hkl;
  int peakNum = 0;
  const auto NPeaks = Peaks->getNumberPeaks();
147
148
  Kernel::Matrix<double> goniometer;
  goniometer.identityMatrix();
149

150
151
  const double lambdaMin = getProperty("WavelengthMin");
  const double lambdaMax = getProperty("WavelengthMax");
152
153
  IPeak &peak0 = Peaks->getPeak(0);

154
155
156
157
  std::vector<V3D> possibleHKLs;
  if (includePeaksInRange) {
    const double dMin = getProperty("MinDSpacing");
    const double dMax = getProperty("MaxDSpacing");
158
159
    Geometry::HKLGenerator gen(lattice, dMin);
    auto filter = boost::make_shared<HKLFilterDRange>(lattice, dMin, dMax);
Lynch, Vickie's avatar
Lynch, Vickie committed
160

161
162
    V3D hkl = *(gen.begin());
    g_log.information() << "HKL range for d_min of " << dMin << " to d_max of "
Lynch, Vickie's avatar
Lynch, Vickie committed
163
164
165
166
167
168
169
170
171
172
                        << dMax << " is from " << hkl << " to " << hkl * -1.0
                        << ", a total of " << gen.size() << " possible HKL's\n";
    if (gen.size() > 10000000000)
      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), (~filter)->fn());
173
174
175
176
177
178
  } else {
    hkl[0] = peak0.getH();
    hkl[1] = peak0.getK();
    hkl[2] = peak0.getL();
  }

Lynch, Vickie's avatar
Lynch, Vickie committed
179
  size_t N = NPeaks * (1 + 2 * maxOrder);
180
181
182
183
184
  if (includePeaksInRange) {
    N = possibleHKLs.size();
    N = max<size_t>(100, N);
  }
  auto RunNumber = peak0.getRunNumber();
185
186
  auto &UB = lattice.getUB();
  goniometer = peak0.getGoniometerMatrix();
187
  Progress prog(this, 0.0, 1.0, N);
188
189
  vector<vector<int>> AlreadyDonePeaks;
  bool done = false;
190
  DblMatrix orientedUB = goniometer * UB;
191
  HKLFilterWavelength lambdaFilter(orientedUB, lambdaMin, lambdaMax);
Lynch, Vickie's avatar
Lynch, Vickie committed
192
  int seqNum = 0;
193
  size_t next = 0;
194
  while (!done) {
195
    std::string offsetName = "Offset1";
Lynch, Vickie's avatar
Lynch, Vickie committed
196
197
198
    predictOffsets(Peaks, OutPeaks, offsets1, offsetName, maxOrder, hkl,
                   goniometer, UB, lambdaFilter, includePeaksInRange,
                   includeOrderZero, RunNumber, seqNum, AlreadyDonePeaks);
199
    offsetName = "Offset2";
Lynch, Vickie's avatar
Lynch, Vickie committed
200
201
202
    predictOffsets(Peaks, OutPeaks, offsets2, offsetName, maxOrder, hkl,
                   goniometer, UB, lambdaFilter, includePeaksInRange,
                   includeOrderZero, RunNumber, seqNum, AlreadyDonePeaks);
203
    offsetName = "Offset3";
Lynch, Vickie's avatar
Lynch, Vickie committed
204
205
206
    predictOffsets(Peaks, OutPeaks, offsets3, offsetName, maxOrder, hkl,
                   goniometer, UB, lambdaFilter, includePeaksInRange,
                   includeOrderZero, RunNumber, seqNum, AlreadyDonePeaks);
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
    if (includePeaksInRange) {
      next++;
      if (next == possibleHKLs.size())
        break;
      hkl = possibleHKLs[next];
    } else {
      peakNum++;
      if (peakNum >= NPeaks)
        done = true;
      else { // peak0= Peaks->getPeak(peakNum);
        IPeak &peak1 = Peaks->getPeak(peakNum);
        //??? could not assign to peak0 above. Did not work
        // the peak that peak0 was associated with did NOT change
        hkl[0] = peak1.getH();
        hkl[1] = peak1.getK();
        hkl[2] = peak1.getL();
        goniometer = peak1.getGoniometerMatrix();
        RunNumber = peak1.getRunNumber();
      }
    }
    prog.report();
  }

  setProperty("SatellitePeaks", OutPeaks);
}

Lynch, Vickie's avatar
Lynch, Vickie committed
233
234
235
236
237
238
239
240
241
void PredictSatellitePeaks::predictOffsets(
    DataObjects::PeaksWorkspace_sptr Peaks,
    boost::shared_ptr<Mantid::API::IPeaksWorkspace> &OutPeaks,
    std::vector<double> offsets, std::string &label, int &maxOrder, V3D &hkl,
    Kernel::Matrix<double> &goniometer, const Kernel::DblMatrix &UB,
    HKLFilterWavelength &lambdaFilter, bool &includePeaksInRange,
    bool &includeOrderZero, int &RunNumber, int &seqNum,
    vector<vector<int>> &AlreadyDonePeaks) {
  OutPeaks->mutableRun().addProperty<std::vector<double>>(label, offsets, true);
242
  Geometry::InstrumentRayTracer tracer(Peaks->getInstrument());
Lynch, Vickie's avatar
Lynch, Vickie committed
243
244
245
246
  for (int order = -maxOrder; order <= maxOrder; order++) {
    if (order == 0 && !includeOrderZero)
      continue; // exclude order 0
    V3D hkl1(hkl);
Lynch, Vickie's avatar
Lynch, Vickie committed
247

Lynch, Vickie's avatar
Lynch, Vickie committed
248
249
250
251
252
    hkl1[0] += order * offsets[0];
    hkl1[1] += order * offsets[1];
    hkl1[2] += order * offsets[2];
    if (!lambdaFilter.isAllowed(hkl1) && includePeaksInRange)
      continue;
253

Lynch, Vickie's avatar
Lynch, Vickie committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
    Kernel::V3D Qs = goniometer * UB * hkl1 * 2.0 * M_PI;

    if (Qs[2] <= 0)
      continue;

    try {
      auto peak(Peaks->createPeak(Qs, 1));

      peak->setGoniometerMatrix(goniometer);

      if (peak->findDetector(tracer)) {
        vector<int> SavPk{RunNumber, boost::math::iround(1000.0 * hkl1[0]),
                          boost::math::iround(1000.0 * hkl1[1]),
                          boost::math::iround(1000.0 * hkl1[2])};
Lynch, Vickie's avatar
Lynch, Vickie committed
268

Lynch, Vickie's avatar
Lynch, Vickie committed
269
270
271
272
273
274
275
        auto it = find(AlreadyDonePeaks.begin(), AlreadyDonePeaks.end(), SavPk);

        if (it == AlreadyDonePeaks.end()) {
          AlreadyDonePeaks.push_back(SavPk);
          std::sort(AlreadyDonePeaks.begin(), AlreadyDonePeaks.end());
        } else {
          continue;
Lynch, Vickie's avatar
Lynch, Vickie committed
276
        }
Lynch, Vickie's avatar
Lynch, Vickie committed
277
278
279
280
281
282
283
284

        peak->setHKL(hkl1);
        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
285
      }
Lynch, Vickie's avatar
Lynch, Vickie committed
286
287
288
    } catch (std::runtime_error &) {
      throw std::invalid_argument("Invalid Q vector");
    }
Lynch, Vickie's avatar
Lynch, Vickie committed
289
  }
290
291
292
293
}

} // namespace Crystal
} // namespace Mantid