LoadILLReflectometry.cpp 39 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 +
7
#include "MantidDataHandling/LoadILLReflectometry.h"
8

9
#include "MantidAPI/Axis.h"
10
#include "MantidAPI/CompositeFunction.h"
Yannick Raoul's avatar
Yannick Raoul committed
11
#include "MantidAPI/FileProperty.h"
12
13
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IPeakFunction.h"
14
#include "MantidAPI/MatrixWorkspace.h"
Yannick Raoul's avatar
Yannick Raoul committed
15
#include "MantidAPI/RegisterFileLoader.h"
16
#include "MantidAPI/Run.h"
17
18
#include "MantidAPI/SpectrumInfo.h"
#include "MantidDataObjects/TableWorkspace.h"
19
20
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
21
22
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
23
#include "MantidKernel/BoundedValidator.h"
24
#include "MantidKernel/ListValidator.h"
Hahn, Steven's avatar
Hahn, Steven committed
25
#include "MantidKernel/OptionalBool.h"
26
#include "MantidKernel/Quat.h"
27
#include "MantidKernel/UnitFactory.h"
28
#include "MantidKernel/V3D.h"
29

30
namespace {
31
32
/// Component coordinates for FIGARO, in meter.
namespace FIGARO {
33
34
constexpr double DH1Z{1.135}; // Motor DH1 horizontal position
constexpr double DH2Z{2.077}; // Motor DH2 horizontal position
LamarMoore's avatar
LamarMoore committed
35
} // namespace FIGARO
36

37
/// A struct for information needed for detector angle calibration.
38
struct PeakInfo {
39
40
  double detectorAngle;
  double detectorDistance;
41
  double peakCentre;
42
43
};

44
45
46
47
/** Convert degrees to radians.
 *  @param x an angle in degrees
 *  @return the angle in radians
 */
48
constexpr double degToRad(const double x) { return x * M_PI / 180.; }
49

50
51
52
53
/** Convert radians to degrees.
 *  @param x an angle in radians
 *  @return the angle in degrees
 */
54
constexpr double radToDeg(const double x) { return x * 180. / M_PI; }
55

56
57
58
59
/** Convert millimeters to meters.
 *  @param x a distance in millimeters
 *  @return the distance in meters
 */
60
constexpr double mmToMeter(const double x) { return x * 1.e-3; }
61

62
63
64
65
/** Create a table with data needed for detector angle calibration.
 * @param info data to be written to the table
 * @return a TableWorkspace containing the beam position info
 */
Antti Soininen's avatar
Antti Soininen committed
66
Mantid::API::ITableWorkspace_sptr
67
createPeakPositionTable(const PeakInfo &info) {
68
  auto table = std::make_shared<Mantid::DataObjects::TableWorkspace>();
69
70
  table->addColumn("double", "DetectorAngle");
  table->addColumn("double", "DetectorDistance");
71
  table->addColumn("double", "PeakCentre");
72
73
74
75
76
  table->appendRow();
  auto col = table->getColumn("DetectorAngle");
  col->cell<double>(0) = info.detectorAngle;
  col = table->getColumn("DetectorDistance");
  col->cell<double>(0) = info.detectorDistance;
77
78
  col = table->getColumn("PeakCentre");
  col->cell<double>(0) = info.peakCentre;
79
80
81
  return table;
}

82
83
84
85
/** Strip monitors from the beginning and end of a workspace.
 *  @param ws a workspace to work on
 *  @return begin and end ws indices for non-monitor histograms
 */
Antti Soininen's avatar
Antti Soininen committed
86
87
std::pair<size_t, size_t>
fitIntegrationWSIndexRange(const Mantid::API::MatrixWorkspace &ws) {
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
  const size_t nHisto = ws.getNumberHistograms();
  size_t begin = 0;
  const auto &spectrumInfo = ws.spectrumInfo();
  for (size_t i = 0; i < nHisto; ++i) {
    if (!spectrumInfo.isMonitor(i)) {
      break;
    }
    ++begin;
  }
  size_t end = nHisto - 1;
  for (ptrdiff_t i = static_cast<ptrdiff_t>(nHisto) - 1; i != 0; --i) {
    if (!spectrumInfo.isMonitor(i)) {
      break;
    }
    --end;
  }
  return std::pair<size_t, size_t>{begin, end};
}

107
108
109
110
/** Construct a DirectBeamMeasurement object from a beam position table.
 *  @param table a beam position TableWorkspace
 *  @return a DirectBeamMeasurement object corresonding to the table parameter.
 */
111
PeakInfo parseBeamPositionTable(const Mantid::API::ITableWorkspace &table) {
112
  if (table.rowCount() != 1) {
113
114
    throw std::runtime_error(
        "DirectBeamPosition table should have a single row.");
115
  }
116
  PeakInfo p;
117
  auto col = table.getColumn("DetectorAngle");
118
  p.detectorAngle = col->cell<double>(0);
119
  col = table.getColumn("DetectorDistance");
120
121
122
123
124
125
  p.detectorDistance = col->cell<double>(0);
  col = table.getColumn("PeakCentre");
  p.peakCentre = col->cell<double>(0);
  return p;
}

126
127
128
/** Fill the X values of the first histogram of ws with values 0, 1, 2,...
 *  @param ws a workspace to modify
 */
129
130
131
132
void rebinIntegralWorkspace(Mantid::API::MatrixWorkspace &ws) {
  auto &xs = ws.mutableX(0);
  std::iota(xs.begin(), xs.end(), 0.0);
}
133

134
/// Enumerations to define the rotation plane of the detector.
Antti Soininen's avatar
Antti Soininen committed
135
enum class RotationPlane { horizontal, vertical };
136

137
138
139
/** Calculate the detector position from given parameters.
 *  @param plane rotation plane of the detector
 *  @param distance sample to detector centre distance in meters
Antti Soininen's avatar
Antti Soininen committed
140
 *  @param angle an angle between the Z axis and the detector in degrees
141
142
 *  @return a vector pointing to the new detector centre
 */
Antti Soininen's avatar
Antti Soininen committed
143
144
145
Mantid::Kernel::V3D detectorPosition(const RotationPlane plane,
                                     const double distance,
                                     const double angle) {
146
  const double a = degToRad(angle);
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
  double x, y, z;
  switch (plane) {
  case RotationPlane::horizontal:
    x = distance * std::sin(a);
    y = 0;
    z = distance * std::cos(a);
    break;
  case RotationPlane::vertical:
    x = 0;
    y = distance * std::sin(a);
    z = distance * std::cos(a);
    break;
  }
  return Mantid::Kernel::V3D(x, y, z);
}

163
164
/** Calculates the detector rotation such that it faces the origin.
 *  @param plane rotation plane of the detectorPosition
Antti Soininen's avatar
Antti Soininen committed
165
 *  @param angle an angle between the Z axis and the detector in degrees
166
167
 *  @return the calculated rotation transformation
 */
Antti Soininen's avatar
Antti Soininen committed
168
169
Mantid::Kernel::Quat detectorFaceRotation(const RotationPlane plane,
                                          const double angle) {
170
171
  const Mantid::Kernel::V3D axis = [plane]() {
    double x, y;
Antti Soininen's avatar
Antti Soininen committed
172
    switch (plane) {
173
174
175
176
177
178
179
180
181
182
183
184
185
    case RotationPlane::horizontal:
      x = 0;
      y = 1;
      break;
    case RotationPlane::vertical:
      x = -1;
      y = 0;
      break;
    }
    return Mantid::Kernel::V3D(x, y, 0);
  }();
  return Mantid::Kernel::Quat(angle, axis);
}
186
} // anonymous namespace
Verena Reimund's avatar
Verena Reimund committed
187

188
189
namespace Mantid {
namespace DataHandling {
Yannick Raoul's avatar
Yannick Raoul committed
190

191
192
193
using namespace Kernel;
using namespace API;
using namespace NeXus;
Yannick Raoul's avatar
Yannick Raoul committed
194
195

// Register the algorithm into the AlgorithmFactory
196
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadILLReflectometry)
Yannick Raoul's avatar
Yannick Raoul committed
197

198
199
const double LoadILLReflectometry::PIXEL_CENTER = 127.5;

200
/**
201
 * Return the confidence with this algorithm can load the file
202
203
204
205
 * @param descriptor A descriptor for the file
 * @returns An integer specifying the confidence level. 0 indicates it will not
 * be used
 */
206
207
int LoadILLReflectometry::confidence(
    Kernel::NexusDescriptor &descriptor) const {
208
209

  // fields existent only at the ILL
210
  if ((descriptor.pathExists("/entry0/wavelength") || // ILL D17
211
       descriptor.pathExists("/entry0/theta"))        // ILL FIGARO
LamarMoore's avatar
LamarMoore committed
212
      && descriptor.pathExists("/entry0/experiment_identifier") &&
213
      descriptor.pathExists("/entry0/mode") &&
214
      (descriptor.pathExists("/entry0/instrument/VirtualChopper") || // ILL D17
215
       descriptor.pathExists("/entry0/instrument/Theta")) // ILL FIGARO
LamarMoore's avatar
LamarMoore committed
216
  )
217
    return 80;
Verena Reimund's avatar
Verena Reimund committed
218
  else
219
    return 0;
Verena Reimund's avatar
Verena Reimund committed
220
221
}

222
/// Initialize the algorithm's properties.
223
void LoadILLReflectometry::init() {
224
  declareProperty(std::make_unique<FileProperty>("Filename", std::string(),
Sam Jenkins's avatar
Sam Jenkins committed
225
226
                                                 FileProperty::Load, ".nxs",
                                                 Direction::Input),
Verena Reimund's avatar
Verena Reimund committed
227
                  "Name of the Nexus file to load");
228

229
  declareProperty(std::make_unique<WorkspaceProperty<>>(
230
                      "OutputWorkspace", std::string(), Direction::Output),
Verena Reimund's avatar
Verena Reimund committed
231
                  "Name of the output workspace");
232
233
234
  declareProperty(
      "BeamCentre", EMPTY_DBL(),
      "Beam position in workspace indices (disables peak finding).");
235
  declareProperty(std::make_unique<WorkspaceProperty<ITableWorkspace>>(
Antti Soininen's avatar
Antti Soininen committed
236
237
238
                      "OutputBeamPosition", std::string(), Direction::Output,
                      PropertyMode::Optional),
                  "Name of the fitted beam position output workspace");
239

240
  declareProperty(std::make_unique<WorkspaceProperty<ITableWorkspace>>(
241
                      "DirectBeamPosition", std::string(), Direction::Input,
Antti Soininen's avatar
Antti Soininen committed
242
243
244
                      PropertyMode::Optional),
                  "A workspace defining the beam position; used to calculate "
                  "the Bragg angle");
245

246
247
  declareProperty("BraggAngle", EMPTY_DBL(),
                  "User defined Bragg angle in degrees");
248
  const std::vector<std::string> availableUnits{"Wavelength", "TimeOfFlight"};
249
  declareProperty("XUnit", "Wavelength",
250
                  std::make_shared<StringListValidator>(availableUnits),
251
                  "X unit of the OutputWorkspace");
252
}
253

254
/// Execute the algorithm.
255
256
void LoadILLReflectometry::exec() {
  // open the root node
257
258
  NeXus::NXRoot root(getPropertyValue("Filename"));
  NXEntry firstEntry{root.openFirstEntry()};
Verena Reimund's avatar
Verena Reimund committed
259
260
  // set instrument specific names of Nexus file entries
  initNames(firstEntry);
261
262
  // load Monitor details: n. monitors x monitor contents
  std::vector<std::vector<int>> monitorsData{loadMonitors(firstEntry)};
263
264
  // load Data details (number of tubes, channels, etc)
  loadDataDetails(firstEntry);
265
266
  // initialise workspace
  initWorkspace(monitorsData);
267
  // load the instrument from the IDF if it exists
268
269
  loadInstrument();
  // get properties
270
  loadNexusEntriesIntoProperties();
271
  // load data into the workspace
272
273
274
  loadData(firstEntry, monitorsData, getXValues());
  root.close();
  firstEntry.close();
275
  initPixelWidth();
276
  // Move components.
277
  m_sampleZOffset = sampleHorizontalOffset();
278
  placeSource();
279
  placeDetector();
280
  placeSlits();
281
  // When other components are in-place
282
  convertTofToWavelength();
283
  // Add sample logs loader.two_theta and Facility
284
  addSampleLogs();
285
286
  // Set the output workspace property
  setProperty("OutputWorkspace", m_localWorkspace);
287
288
289
290
291
292
293
294
} // exec

/// Run the Child Algorithm LoadInstrument.
void LoadILLReflectometry::loadInstrument() {
  // execute the Child Algorithm. Catch and log any error, but don't stop.
  g_log.debug("Loading instrument definition...");
  try {
    IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");
295
    const std::string instrumentName =
296
        m_instrument == Supported::D17 ? "D17" : "FIGARO";
297
    loadInst->setPropertyValue("InstrumentName", instrumentName);
298
299
300
301
302
    loadInst->setProperty("RewriteSpectraMap",
                          Mantid::Kernel::OptionalBool(true));
    loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", m_localWorkspace);
    loadInst->executeAsChildAlg();
  } catch (std::runtime_error &e) {
Antti Soininen's avatar
Antti Soininen committed
303
304
305
    g_log.information()
        << "Unable to succesfully run LoadInstrument child algorithm: "
        << e.what() << '\n';
306
  }
307
308
}

Verena Reimund's avatar
Verena Reimund committed
309
/**
LamarMoore's avatar
LamarMoore committed
310
311
312
313
314
 * Init names of sample logs based on instrument specific NeXus file
 * entries
 *
 * @param entry :: the NeXus file entry
 */
Verena Reimund's avatar
Verena Reimund committed
315
void LoadILLReflectometry::initNames(NeXus::NXEntry &entry) {
316
  std::string instrumentNamePath = m_loader.findInstrumentNexusPath(entry);
317
318
  std::string instrumentName =
      entry.getString(instrumentNamePath.append("/name"));
319
  if (instrumentName.empty())
320
321
    throw std::runtime_error(
        "Cannot set the instrument name from the Nexus file!");
322
323
324
325
  boost::to_lower(instrumentName);
  if (instrumentName == "d17") {
    m_instrument = Supported::D17;
  } else if (instrumentName == "figaro") {
326
    m_instrument = Supported::FIGARO;
327
328
  } else {
    std::ostringstream str;
329
    str << "Unsupported instrument: " << instrumentName << '.';
330
331
332
333
    throw std::runtime_error(str.str());
  }
  g_log.debug() << "Instrument name: " << instrumentName << '\n';
  if (m_instrument == Supported::D17) {
Verena Reimund's avatar
Verena Reimund committed
334
335
336
337
    m_detectorAngleName = "dan.value";
    m_offsetFrom = "VirtualChopper";
    m_chopper1Name = "Chopper1";
    m_chopper2Name = "Chopper2";
338
  } else if (m_instrument == Supported::FIGARO) {
339
    m_detectorAngleName = "VirtualAxis.DAN_actual_angle";
340
    m_sampleAngleName = "CollAngle.actual_coll_angle";
Verena Reimund's avatar
Verena Reimund committed
341
    m_offsetFrom = "CollAngle";
342
    // FIGARO: find out which of the four choppers are used
Verena Reimund's avatar
Verena Reimund committed
343
344
345
346
347
348
    NXFloat firstChopper =
        entry.openNXFloat("instrument/ChopperSetting/firstChopper");
    firstChopper.load();
    NXFloat secondChopper =
        entry.openNXFloat("instrument/ChopperSetting/secondChopper");
    secondChopper.load();
349
350
    m_chopper1Name = "CH" + std::to_string(int(firstChopper[0]));
    m_chopper2Name = "CH" + std::to_string(int(secondChopper[0]));
Verena Reimund's avatar
Verena Reimund committed
351
  }
352
353
354
355
356
  // get acquisition mode
  NXInt acqMode = entry.openNXInt("acquisition_mode");
  acqMode.load();
  m_acqMode = acqMode[0];
  m_acqMode ? g_log.debug("TOF mode") : g_log.debug("Monochromatic Mode");
Verena Reimund's avatar
Verena Reimund committed
357
358
}

359
360
361
362
363
364
365
366
367
368
369
370
371
/// Call child algorithm ConvertUnits for conversion from TOF to wavelength
void LoadILLReflectometry::convertTofToWavelength() {
  if (m_acqMode && (getPropertyValue("XUnit") == "Wavelength")) {
    auto convertToWavelength =
        createChildAlgorithm("ConvertUnits", -1, -1, true);
    convertToWavelength->initialize();
    convertToWavelength->setProperty<MatrixWorkspace_sptr>("InputWorkspace",
                                                           m_localWorkspace);
    convertToWavelength->setProperty<MatrixWorkspace_sptr>("OutputWorkspace",
                                                           m_localWorkspace);
    convertToWavelength->setPropertyValue("Target", "Wavelength");
    convertToWavelength->executeAsChildAlg();
  }
372
373
}

374
375
376
377
378
379
/**
 * Creates the workspace and initialises member variables with
 * the corresponding values
 *
 * @param monitorsData :: Monitors data already loaded
 */
Verena Reimund's avatar
Verena Reimund committed
380
void LoadILLReflectometry::initWorkspace(
381
    const std::vector<std::vector<int>> &monitorsData) {
382

383
  g_log.debug() << "Number of monitors: " << monitorsData.size() << '\n';
384
385
  for (size_t i = 0; i < monitorsData.size(); ++i) {
    if (monitorsData[i].size() != m_numberOfChannels)
Antti Soininen's avatar
Antti Soininen committed
386
387
      g_log.debug() << "Data size of monitor ID " << i << " is "
                    << monitorsData[i].size() << '\n';
388
389
390
  }
  // create the workspace
  try {
391
392
393
    m_localWorkspace = DataObjects::create<DataObjects::Workspace2D>(
        m_numberOfHistograms + monitorsData.size(),
        HistogramData::BinEdges(m_numberOfChannels + 1));
394
  } catch (std::out_of_range &) {
Verena Reimund's avatar
Verena Reimund committed
395
    throw std::runtime_error(
396
        "Workspace2D cannot be created, check number of histograms (" +
Antti Soininen's avatar
Antti Soininen committed
397
398
399
        std::to_string(m_numberOfHistograms) + "), monitors (" +
        std::to_string(monitorsData.size()) + "), and channels (" +
        std::to_string(m_numberOfChannels) + '\n');
400
401
  }
  if (m_acqMode)
Antti Soininen's avatar
Antti Soininen committed
402
403
    m_localWorkspace->getAxis(0)->unit() =
        UnitFactory::Instance().create("TOF");
404
  m_localWorkspace->setYUnitLabel("Counts");
405
406
407
408
}

/**
 * Load Data details (number of tubes, channels, etc)
409
 *
410
411
412
 * @param entry First entry of nexus file
 */
void LoadILLReflectometry::loadDataDetails(NeXus::NXEntry &entry) {
413
  // PSD data D17 256 x 1 x 1000
414
  // PSD data FIGARO 1 x 256 x 1000
415

Verena Reimund's avatar
Verena Reimund committed
416
417
418
419
420
421
  if (m_acqMode) {
    NXFloat timeOfFlight = entry.openNXFloat("instrument/PSD/time_of_flight");
    timeOfFlight.load();
    m_channelWidth = static_cast<double>(timeOfFlight[0]);
    m_numberOfChannels = size_t(timeOfFlight[1]);
    m_tofDelay = timeOfFlight[2];
Verena Reimund's avatar
Verena Reimund committed
422
  } else { // monochromatic mode
Verena Reimund's avatar
Verena Reimund committed
423
424
    m_numberOfChannels = 1;
  }
425
426
427
428
429

  NXInt nChannels = entry.openNXInt("instrument/PSD/detsize");
  nChannels.load();
  m_numberOfHistograms = nChannels[0];

Antti Soininen's avatar
Antti Soininen committed
430
431
432
433
  g_log.debug()
      << "Please note that ILL reflectometry instruments have "
         "several tubes, after integration one "
         "tube remains in the Nexus file.\n Number of tubes (banks): 1\n";
434
  g_log.debug() << "Number of pixels per tube (number of detectors and number "
LamarMoore's avatar
LamarMoore committed
435
436
                   "of histograms): "
                << m_numberOfHistograms << '\n';
437
438
439
440
  g_log.debug() << "Number of time channels: " << m_numberOfChannels << '\n';
  g_log.debug() << "Channel width: " << m_channelWidth << " 1e-6 sec\n";
}

441
442
443
444
445
446
447
/**
 * @brief LoadILLReflectometry::doubleFromRun
 * Returns a sample log a single double
 * @param entryName : the name of the log
 * @return the value as double
 * @throws runtime_error : if the log does not exist
 */
448
double LoadILLReflectometry::doubleFromRun(const std::string &entryName) const {
449
450
451
452
453
454
  if (m_localWorkspace->run().hasProperty(entryName)) {
    return m_localWorkspace->run().getPropertyValueAsType<double>(entryName);
  } else {
    throw std::runtime_error("The log with the given name does not exist " +
                             entryName);
  }
455
456
457
}

/**
Verena Reimund's avatar
Verena Reimund committed
458
 * Load single monitor
459
460
 *
 * @param entry :: The Nexus entry
461
 * @param monitor_data :: A std::string containing the Nexus path to the monitor
Verena Reimund's avatar
Verena Reimund committed
462
 *data
Verena Reimund's avatar
Verena Reimund committed
463
 * @return monitor :: A std::vector containing monitor values
464
 */
Verena Reimund's avatar
Verena Reimund committed
465
std::vector<int>
Verena Reimund's avatar
Verena Reimund committed
466
LoadILLReflectometry::loadSingleMonitor(NeXus::NXEntry &entry,
467
                                        const std::string &monitor_data) {
Verena Reimund's avatar
Verena Reimund committed
468
  NXData dataGroup = entry.openNXData(monitor_data);
469
  NXInt data = dataGroup.openIntData();
Verena Reimund's avatar
Verena Reimund committed
470
  // load counts
471
  data.load();
472
  return std::vector<int>(data(), data() + data.size());
Verena Reimund's avatar
Verena Reimund committed
473
474
475
}

/**
476
 * Load monitor data
Verena Reimund's avatar
Verena Reimund committed
477
478
 *
 * @param entry :: The Nexus entry
479
 * @return :: A std::vector of vectors of monitors containing monitor values
Verena Reimund's avatar
Verena Reimund committed
480
481
482
 */
std::vector<std::vector<int>>
LoadILLReflectometry::loadMonitors(NeXus::NXEntry &entry) {
483
  g_log.debug("Read monitor data...");
Verena Reimund's avatar
Verena Reimund committed
484
  // vector of monitors with one entry
485
  const std::vector<std::vector<int>> monitors{
Verena Reimund's avatar
Verena Reimund committed
486
487
      loadSingleMonitor(entry, "monitor1/data"),
      loadSingleMonitor(entry, "monitor2/data")};
488
489
490
491
  return monitors;
}

/**
492
 * Determine x values (unit time-of-flight)
493
 *
494
 * @return :: vector holding the x values
495
 */
496
std::vector<double> LoadILLReflectometry::getXValues() {
Antti Soininen's avatar
Antti Soininen committed
497
  std::vector<double> xVals;             // no initialisation
498
  xVals.reserve(m_numberOfChannels + 1); // reserve memory
Verena Reimund's avatar
Verena Reimund committed
499
  try {
500
    if (m_acqMode) {
501
      if (m_instrument == Supported::FIGARO) {
502
503
504
        if (m_localWorkspace->run().hasProperty(
                "Distance.edelay_delay")) // Valid from 2018.
          m_tofDelay += doubleFromRun("Distance.edelay_delay");
505
        else // Valid before 2018.
506
507
508
          m_tofDelay += doubleFromRun("Theta.edelay_delay");
      }
      g_log.debug() << "TOF delay: " << m_tofDelay << '\n';
509
      std::string chopper{"Chopper"};
Antti Soininen's avatar
Antti Soininen committed
510
      double chop1Speed{0.0}, chop1Phase{0.0}, chop2Speed{0.0}, chop2Phase{0.0};
511
      if (m_instrument == Supported::D17) {
512
        chop1Speed = doubleFromRun("VirtualChopper.chopper1_speed_average");
513
514
        chop1Phase = doubleFromRun("VirtualChopper.chopper1_phase_average");
        chop2Speed = doubleFromRun("VirtualChopper.chopper2_speed_average");
515
        chop2Phase = doubleFromRun("VirtualChopper.chopper2_phase_average");
516
517
518
519
520
        if (chop1Phase > 360.) {
          // This is an ugly workaround for pre-2018 D17 files which have
          // chopper 1 phase and chopper 2 speed swapped.
          std::swap(chop1Phase, chop2Speed);
        }
521
      } else if (m_instrument == Supported::FIGARO) {
Antti Soininen's avatar
Antti Soininen committed
522
        chop1Phase = doubleFromRun(m_chopper1Name + ".phase");
523
        // Chopper 1 phase on FIGARO is set to an arbitrary value (999.9)
524
525
        if (chop1Phase > 360.0)
          chop1Phase = 0.0;
526
      }
527
528
529
530
531
      double POFF;
      try {
        POFF = doubleFromRun(m_offsetFrom + ".poff");
      } catch (std::runtime_error &) {
        try {
532
          POFF = doubleFromRun(m_offsetFrom + ".pickup_offset");
533
534
535
536
537
        } catch (std::runtime_error &) {
          throw std::runtime_error(
              "Unable to find VirtualChopper pickup offset");
        }
      }
538
539
540
541
542
543
      double openOffset;
      if (m_localWorkspace->run().hasProperty(
              m_offsetFrom + ".open_offset")) // Valid from 2018.
        openOffset = doubleFromRun(m_offsetFrom + ".open_offset");
      else // Figaro 2017 / 2018
        openOffset = doubleFromRun(m_offsetFrom + ".openOffset");
544
545
      if (m_instrument == Supported::D17 && chop1Speed != 0.0 &&
          chop2Speed != 0.0 && chop2Phase != 0.0) {
546
547
548
549
        // virtual chopper entries are valid
        chopper = "Virtual chopper";
      } else {
        // use chopper values
550
        chop1Speed = doubleFromRun(m_chopper1Name + ".rotation_speed");
Antti Soininen's avatar
Antti Soininen committed
551
        chop2Speed = doubleFromRun(m_chopper2Name + ".rotation_speed");
552
        chop2Phase = doubleFromRun(m_chopper2Name + ".phase");
553
554
      }
      // logging
555
556
557
558
559
560
      g_log.debug() << "Poff: " << POFF << '\n';
      g_log.debug() << "Open offset: " << openOffset << '\n';
      g_log.debug() << "Chopper 1 phase: " << chop1Phase << '\n';
      g_log.debug() << chopper << " 1 speed: " << chop1Speed << '\n';
      g_log.debug() << chopper << " 2 phase: " << chop2Phase << '\n';
      g_log.debug() << chopper << " 2 speed: " << chop2Speed << '\n';
561

562
      if (chop1Speed <= 0.0) {
Antti Soininen's avatar
Antti Soininen committed
563
564
        g_log.error() << "First chopper velocity " << chop1Speed
                      << ". Check you NeXus file.\n";
565
      }
Antti Soininen's avatar
Antti Soininen committed
566
      const double chopWindow = 45.0;
LamarMoore's avatar
LamarMoore committed
567
568
569
570
      const double t_TOF2 = m_tofDelay - 1.e+6 * 60.0 *
                                             (POFF - chopWindow + chop2Phase -
                                              chop1Phase + openOffset) /
                                             (2.0 * 360 * chop1Speed);
571
      g_log.debug() << "t_TOF2: " << t_TOF2 << '\n';
572
      // compute tof values
Antti Soininen's avatar
Antti Soininen committed
573
      for (int channelIndex = 0;
574
           channelIndex < static_cast<int>(m_numberOfChannels) + 1;
575
           ++channelIndex) {
576
577
        const double t_TOF1 = channelIndex * m_channelWidth;
        xVals.emplace_back(t_TOF1 + t_TOF2);
578
      }
Verena Reimund's avatar
Verena Reimund committed
579
    } else {
580
581
      g_log.debug("Time channel index for axis description \n");
      for (size_t t = 0; t <= m_numberOfChannels; ++t)
582
        xVals.emplace_back(static_cast<double>(t));
Verena Reimund's avatar
Verena Reimund committed
583
584
    }
  } catch (std::runtime_error &e) {
Antti Soininen's avatar
Antti Soininen committed
585
586
    g_log.information() << "Unable to access NeXus file entry: " << e.what()
                        << '\n';
Verena Reimund's avatar
Verena Reimund committed
587
  }
588
  return xVals;
589
590
591
592
593
594
}

/**
 * Load data from nexus file
 *
 * @param entry :: The Nexus file entry
595
 * @param monitorsData :: Monitors data already loaded
596
 * @param xVals :: X values
597
 */
Antti Soininen's avatar
Antti Soininen committed
598
599
600
void LoadILLReflectometry::loadData(
    NeXus::NXEntry &entry, const std::vector<std::vector<int>> &monitorsData,
    const std::vector<double> &xVals) {
601
  g_log.debug("Loading data...");
602
603
604
605
  NXData dataGroup = entry.openNXData("data");
  NXInt data = dataGroup.openIntData();
  // load the counts from the file into memory
  data.load();
606
  const size_t nb_monitors = monitorsData.size();
607
  Progress progress(this, 0, 1, m_numberOfHistograms + nb_monitors);
608

609
  // load data
610
611
  if (!xVals.empty()) {
    HistogramData::BinEdges binEdges(xVals);
612
613
    PARALLEL_FOR_IF(Kernel::threadSafe(*m_localWorkspace))
    for (int j = 0; j < static_cast<int>(m_numberOfHistograms); ++j) {
614
      const int *data_p = &data(0, static_cast<int>(j), 0);
615
      const HistogramData::Counts counts(data_p, data_p + m_numberOfChannels);
616
      m_localWorkspace->setHistogram(j, binEdges, std::move(counts));
617
618
619
620
621
622
623
624
625
626
627
628
      m_localWorkspace->getSpectrum(j).setSpectrumNo(j);
      progress.report();
    }
    for (size_t im = 0; im < nb_monitors; ++im) {
      const int *monitor_p = monitorsData[im].data();
      const HistogramData::Counts monitorCounts(monitor_p,
                                                monitor_p + m_numberOfChannels);
      const size_t spectrum = im + m_numberOfHistograms;
      m_localWorkspace->setHistogram(spectrum, binEdges,
                                     std::move(monitorCounts));
      m_localWorkspace->getSpectrum(spectrum).setSpectrumNo(
          static_cast<specnum_t>(spectrum));
629
630
631
632
      progress.report();
    }
  } else
    g_log.debug("Vector of x values is empty");
633
}
634
635
636

/**
 * Use the LoadHelper utility to load most of the nexus entries into workspace
637
 * sample log properties
638
 */
639
void LoadILLReflectometry::loadNexusEntriesIntoProperties() {
640
  g_log.debug("Building properties...");
641
  // Open NeXus file
642
  const std::string filename{getPropertyValue("Filename")};
643
  NXhandle nxfileID;
644
  NXstatus stat = NXopen(filename.c_str(), NXACC_READ, &nxfileID);
Verena Reimund's avatar
Verena Reimund committed
645
  if (stat == NX_ERROR)
646
    throw Kernel::Exception::FileError("Unable to open File:", filename);
647
  m_loader.addNexusFieldsToWsRun(nxfileID, m_localWorkspace->mutableRun());
648
  NXclose(&nxfileID);
649
650
}

651
/**
LamarMoore's avatar
LamarMoore committed
652
653
654
655
656
 * Gaussian fit to determine peak position if no user position given.
 *
 * @return :: detector position of the peak: Gaussian fit and position
 * of the maximum (serves as start value for the optimization)
 */
657
double LoadILLReflectometry::reflectometryPeak() {
658
659
  if (!isDefault("BeamCentre")) {
    return getProperty("BeamCentre");
660
  }
661
662
  size_t startIndex;
  size_t endIndex;
Antti Soininen's avatar
Antti Soininen committed
663
664
  std::tie(startIndex, endIndex) =
      fitIntegrationWSIndexRange(*m_localWorkspace);
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
  IAlgorithm_sptr integration = createChildAlgorithm("Integration");
  integration->initialize();
  integration->setProperty("InputWorkspace", m_localWorkspace);
  integration->setProperty("OutputWorkspace", "__unused_for_child");
  integration->setProperty("StartWorkspaceIndex", static_cast<int>(startIndex));
  integration->setProperty("EndWorkspaceIndex", static_cast<int>(endIndex));
  integration->execute();
  MatrixWorkspace_sptr integralWS = integration->getProperty("OutputWorkspace");
  IAlgorithm_sptr transpose = createChildAlgorithm("Transpose");
  transpose->initialize();
  transpose->setProperty("InputWorkspace", integralWS);
  transpose->setProperty("OutputWorkspace", "__unused_for_child");
  transpose->execute();
  integralWS = transpose->getProperty("OutputWorkspace");
  rebinIntegralWorkspace(*integralWS);
  // determine initial height: maximum value
  const auto maxValueIt =
      std::max_element(integralWS->y(0).cbegin(), integralWS->y(0).cend());
  const double height = *maxValueIt;
  // determine initial centre: index of the maximum value
  const size_t maxIndex = std::distance(integralWS->y(0).cbegin(), maxValueIt);
686
  const auto centreByMax = static_cast<double>(maxIndex);
687
  g_log.debug() << "Peak maximum position: " << centreByMax << '\n';
688
  // determine sigma
689
  const auto &ys = integralWS->y(0);
690
691
692
  auto lessThanHalfMax = [height](const double x) { return x < 0.5 * height; };
  using IterType = HistogramData::HistogramY::const_iterator;
  std::reverse_iterator<IterType> revMaxValueIt{maxValueIt};
693
694
  auto revMinFwhmIt = std::find_if(revMaxValueIt, ys.crend(), lessThanHalfMax);
  auto maxFwhmIt = std::find_if(maxValueIt, ys.cend(), lessThanHalfMax);
695
  std::reverse_iterator<IterType> revMaxFwhmIt{maxFwhmIt};
696
  if (revMinFwhmIt == ys.crend() || maxFwhmIt == ys.cend()) {
697
698
    g_log.warning() << "Couldn't determine fwhm of beam, using position of max "
                       "value as beam center.\n";
699
700
    return centreByMax;
  }
701
  const auto fwhm =
Antti Soininen's avatar
Antti Soininen committed
702
      static_cast<double>(std::distance(revMaxFwhmIt, revMinFwhmIt) + 1);
703
  g_log.debug() << "Initial fwhm (full width at half maximum): " << fwhm
Antti Soininen's avatar
Antti Soininen committed
704
                << '\n';
705
  // generate Gaussian
706
707
  auto func =
      API::FunctionFactory::Instance().createFunction("CompositeFunction");
708
  auto sum = std::dynamic_pointer_cast<API::CompositeFunction>(func);
709
  func = API::FunctionFactory::Instance().createFunction("Gaussian");
710
  auto gaussian = std::dynamic_pointer_cast<API::IPeakFunction>(func);
711
712
713
714
715
716
717
718
  gaussian->setHeight(height);
  gaussian->setCentre(centreByMax);
  gaussian->setFwhm(fwhm);
  sum->addFunction(gaussian);
  func = API::FunctionFactory::Instance().createFunction("LinearBackground");
  func->setParameter("A0", 0.);
  func->setParameter("A1", 0.);
  sum->addFunction(func);
719
  // call Fit child algorithm
720
721
  API::IAlgorithm_sptr fit = createChildAlgorithm("Fit");
  fit->initialize();
722
  fit->setProperty("Function", std::dynamic_pointer_cast<API::IFunction>(sum));
723
724
725
726
727
728
  fit->setProperty("InputWorkspace", integralWS);
  fit->setProperty("StartX", centreByMax - 3 * fwhm);
  fit->setProperty("EndX", centreByMax + 3 * fwhm);
  fit->execute();
  const std::string fitStatus = fit->getProperty("OutputStatus");
  if (fitStatus != "success") {
729
    g_log.warning("Fit not successful, using position of max value.\n");
730
731
732
733
    return centreByMax;
  }
  const auto centre = gaussian->centre();
  g_log.debug() << "Sigma: " << gaussian->fwhm() << '\n';
734
735
  g_log.debug() << "Estimated peak position: " << centre << '\n';
  return centre;
736
737
}

738
739
740
/// Add sample logs reduction.two_theta and Facility
void LoadILLReflectometry::addSampleLogs() {
  // Add two theta to the sample logs
741
  m_localWorkspace->mutableRun().addProperty("loader.two_theta",
742
743
                                             m_detectorAngle);
  m_localWorkspace->mutableRun()
744
      .getProperty("loader.two_theta")
745
746
747
748
749
      ->setUnits("degree");
  // Add Facility to the sample logs
  m_localWorkspace->mutableRun().addProperty("Facility", std::string("ILL"));
}

750
751
752
/** Compute the detector rotation angle around origin and optionally set the
 *  OutputBeamPosition property.
 *  @return a rotation angle
753
 */
754
double LoadILLReflectometry::detectorRotation() {
755
756
  ITableWorkspace_const_sptr posTable = getProperty("DirectBeamPosition");
  const double peakCentre = reflectometryPeak();
757
758
  g_log.debug() << "Using detector angle (degrees): " << m_detectorAngle
                << '\n';
759
760
  if (!isDefault("OutputBeamPosition")) {
    PeakInfo p;
761
762
    p.detectorAngle = m_detectorAngle;
    p.detectorDistance = m_detectorDistance;
763
764
765
    p.peakCentre = peakCentre;
    setProperty("OutputBeamPosition", createPeakPositionTable(p));
  }
766
  if (!isDefault("BraggAngle")) {
767
    if (posTable) {
768
769
      g_log.notice()
          << "Ignoring DirectBeamPosition, using BraggAngle instead.";
770
    }
771
772
773
774
    const double userAngle = getProperty("BraggAngle");
    const double offset =
        offsetAngle(peakCentre, PIXEL_CENTER, m_detectorDistance);
    m_log.debug() << "Beam offset angle: " << offset << '\n';
775
    return 2. * userAngle - offset;
776
  } else if (!posTable) {
777
    const double deflection = collimationAngle();
778
    if (deflection != 0) {
779
780
      g_log.debug() << "Using incident deflection angle (degrees): "
                    << deflection << '\n';
781
    }
782
    return m_detectorAngle + deflection;
783
  }
784
785
786
787
788
789
790
  const auto dbPeak = parseBeamPositionTable(*posTable);
  const double dbOffset =
      offsetAngle(dbPeak.peakCentre, PIXEL_CENTER, dbPeak.detectorDistance);
  m_log.debug() << "Direct beam offset angle: " << dbOffset << '\n';
  const double detectorAngle =
      m_detectorAngle - dbPeak.detectorAngle - dbOffset;

791
792
  m_log.debug() << "Direct beam calibrated detector angle (degrees): "
                << detectorAngle << '\n';
793
  return detectorAngle;
794
}
Verena Reimund's avatar
Verena Reimund committed
795

796
/// Initialize m_pixelWidth from the IDF and check for NeXus consistency.
797
798
799
800
801
802
void LoadILLReflectometry::initPixelWidth() {
  auto instrument = m_localWorkspace->getInstrument();
  auto detectorPanels = instrument->getAllComponentsWithName("detector");
  if (detectorPanels.size() != 1) {
    throw std::runtime_error("IDF should have a single 'detector' component.");
  }
803
  auto detector =
804
      std::dynamic_pointer_cast<const Geometry::RectangularDetector>(
805
          detectorPanels.front());
806
  double widthInLogs;
807
  if (m_instrument != Supported::FIGARO) {
808
    m_pixelWidth = std::abs(detector->xstep());
809
    widthInLogs = mmToMeter(
810
        m_localWorkspace->run().getPropertyValueAsType<double>("PSD.mppx"));
811
    if (std::abs(widthInLogs - m_pixelWidth) > 1e-10) {
Gagik Vardanyan's avatar
Gagik Vardanyan committed
812
813
814
      m_log.information() << "NeXus pixel width (mppx) " << widthInLogs
                          << " differs from the IDF. Using the IDF value "
                          << m_pixelWidth << '\n';
815
816
817
    }
  } else {
    m_pixelWidth = std::abs(detector->ystep());
818
    widthInLogs = mmToMeter(
819
        m_localWorkspace->run().getPropertyValueAsType<double>("PSD.mppy"));
820
    if (std::abs(widthInLogs - m_pixelWidth) > 1e-10) {
Gagik Vardanyan's avatar
Gagik Vardanyan committed
821
822
823
      m_log.information() << "NeXus pixel width (mppy) " << widthInLogs
                          << " differs from the IDF. Using the IDF value "
                          << m_pixelWidth << '\n';
824
    }
825
826
827
  }
}

828
/// Update detector position according to data file
829
void LoadILLReflectometry::placeDetector() {
Verena Reimund's avatar
Verena Reimund committed
830
  g_log.debug("Move the detector bank \n");
831
  m_detectorDistance = sampleDetectorDistance();
832
  m_detectorAngle = detectorAngle();
833
  g_log.debug() << "Sample-detector distance: " << m_detectorDistance << "m.\n";
834
  const auto detectorRotationAngle = detectorRotation();
835
  const std::string componentName = "detector";
836
  const RotationPlane rotPlane = [this]() {
837
    if (m_instrument != Supported::FIGARO)
Antti Soininen's avatar
Antti Soininen committed
838
839
      return RotationPlane::horizontal;
    else
840
      return RotationPlane::vertical;
841
  }();
Antti Soininen's avatar
Antti Soininen committed
842
  const auto newpos =
843
      detectorPosition(rotPlane, m_detectorDistance, detectorRotationAngle);
844
845
  m_loader.moveComponent(m_localWorkspace, componentName, newpos);
  // apply a local rotation to stay perpendicular to the beam
846
  const auto rotation = detectorFaceRotation(rotPlane, detectorRotationAngle);
847
  m_loader.rotateComponent(m_localWorkspace, componentName, rotation);
848
}