LoadILLTOF2.cpp 14.2 KB
Newer Older
1
2
3
4
5
6
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//     NScD Oak Ridge National Laboratory, European Spallation Source
//     & Institut Laue - Langevin
// SPDX - License - Identifier: GPL - 3.0 +
7
8
9
10
11
12
13
14
#include "MantidDataHandling/LoadILLTOF2.h"

#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidGeometry/Instrument.h"
Hahn, Steven's avatar
Hahn, Steven committed
15
#include "MantidKernel/OptionalBool.h"
16
17
#include "MantidKernel/UnitFactory.h"

18
19
#include <boost/algorithm/string/predicate.hpp>

20
21
22
23
24
25
26
27
28
29
30
namespace Mantid {
namespace DataHandling {

using namespace Kernel;
using namespace API;
using namespace NeXus;
using namespace HistogramData;

DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadILLTOF2)

/// A vector containing the supported instrument names
Antti Soininen's avatar
Antti Soininen committed
31
32
const std::vector<std::string> LoadILLTOF2::SUPPORTED_INSTRUMENTS = {
    "IN4", "IN5", "IN6"};
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

/**
 * Return the confidence with with this algorithm can load the file
 *
 * @param descriptor A descriptor for the file
 *
 * @return An integer specifying the confidence level. 0 indicates it will not
 * be used
 */
int LoadILLTOF2::confidence(Kernel::NexusDescriptor &descriptor) const {

  // fields existent only at the ILL
  if (descriptor.pathExists("/entry0/wavelength") &&
      descriptor.pathExists("/entry0/experiment_identifier") &&
      descriptor.pathExists("/entry0/mode") &&
Antti Soininen's avatar
Antti Soininen committed
48
49
      !descriptor.pathExists("/entry0/dataSD") // This one is for
                                               // LoadILLIndirect
LamarMoore's avatar
LamarMoore committed
50
51
52
53
      && !descriptor.pathExists(
             "/entry0/instrument/VirtualChopper") // This one is for
                                                  // LoadILLReflectometry
  ) {
54
55
56
57
58
59
    return 80;
  } else {
    return 0;
  }
}

Antti Soininen's avatar
Antti Soininen committed
60
LoadILLTOF2::LoadILLTOF2() : API::IFileLoader<Kernel::NexusDescriptor>() {}
61
62
63
64
65
66
67
68
69

/**
 * Initialises the algorithm
 */
void LoadILLTOF2::init() {
  declareProperty(
      make_unique<FileProperty>("Filename", "", FileProperty::Load, ".nxs"),
      "File path of the Data file to load");

Antti Soininen's avatar
Antti Soininen committed
70
71
  declareProperty(make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
                                                   Direction::Output),
72
73
74
75
76
77
78
79
                  "The name to use for the output workspace");
}

/**
 * Executes the algorithm
 */
void LoadILLTOF2::exec() {
  // Retrieve filename
80
  const std::string filenameData = getPropertyValue("Filename");
81
82
83
84
85
86
87
88

  // open the root node
  NeXus::NXRoot dataRoot(filenameData);
  NXEntry dataFirstEntry = dataRoot.openFirstEntry();

  loadInstrumentDetails(dataFirstEntry);
  loadTimeDetails(dataFirstEntry);

89
  const auto monitors = getMonitorInfo(dataFirstEntry);
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

  initWorkSpace(dataFirstEntry, monitors);

  addAllNexusFieldsAsProperties(filenameData);
  addFacility();

  runLoadInstrument(); // just to get IDF contents

  loadDataIntoTheWorkSpace(dataFirstEntry, monitors);

  addEnergyToRun();
  addPulseInterval();

  // load the instrument from the IDF if it exists
  runLoadInstrument();

  // Set the output workspace property
  setProperty("OutputWorkspace", m_localWorkspace);
}

/**
 * Loads Monitor data into an vector of monitor data
 *
 * @param firstEntry The NeXus entry
 *
 * @return List of monitor data
 */
Antti Soininen's avatar
Antti Soininen committed
117
std::vector<std::vector<int>>
118
119
LoadILLTOF2::getMonitorInfo(NeXus::NXEntry &firstEntry) {

Antti Soininen's avatar
Antti Soininen committed
120
  std::vector<std::vector<int>> monitorList;
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151

  for (std::vector<NXClassInfo>::const_iterator it =
           firstEntry.groups().begin();
       it != firstEntry.groups().end(); ++it) {

    if (it->nxclass == "NXmonitor" ||
        boost::starts_with(it->nxname, "monitor")) {

      g_log.debug() << "Load monitor data from " + it->nxname;

      NXData dataGroup = firstEntry.openNXData(it->nxname + "/data");
      NXInt data = dataGroup.openIntData();
      // load the counts from the file into memory
      data.load();

      std::vector<int> thisMonitor(data(), data() + data.size());
      monitorList.push_back(thisMonitor);
    }
  }
  return monitorList;
}

/**
 * Sets the instrument name along with its path in the nexus file
 *
 * @param firstEntry The NeXus entry
 */
void LoadILLTOF2::loadInstrumentDetails(NeXus::NXEntry &firstEntry) {

  m_instrumentPath = m_loader.findInstrumentNexusPath(firstEntry);

152
  if (m_instrumentPath.empty()) {
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
    throw std::runtime_error(
        "Cannot set the instrument name from the Nexus file!");
  }

  m_instrumentName =
      m_loader.getStringFromNexusPath(firstEntry, m_instrumentPath + "/name");

  if (std::find(SUPPORTED_INSTRUMENTS.begin(), SUPPORTED_INSTRUMENTS.end(),
                m_instrumentName) == SUPPORTED_INSTRUMENTS.end()) {
    std::string message =
        "The instrument " + m_instrumentName + " is not valid for this loader!";
    throw std::runtime_error(message);
  }

  g_log.debug() << "Instrument name set to: " + m_instrumentName << '\n';
}

/**
 * Creates the workspace and initialises member variables with
 * the corresponding values
 *
 * @param entry The NeXus entry
 * @param monitors List of monitor data
 */
Antti Soininen's avatar
Antti Soininen committed
177
178
void LoadILLTOF2::initWorkSpace(NeXus::NXEntry &entry,
                                const std::vector<std::vector<int>> &monitors) {
179
180
181
182
183
184
185
186

  // read in the data
  NXData dataGroup = entry.openNXData("data");
  NXInt data = dataGroup.openIntData();

  m_numberOfTubes = static_cast<size_t>(data.dim0());
  m_numberOfPixelsPerTube = static_cast<size_t>(data.dim1());
  m_numberOfChannels = static_cast<size_t>(data.dim2());
187
  const size_t numberOfMonitors = monitors.size();
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
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260

  /**
   * IN4 : Rosace detector is in a different field.
   */
  size_t numberOfTubesInRosace = 0;
  if (m_instrumentName == "IN4") {
    NXData dataGroupRosace =
        entry.openNXData("instrument/Detector_Rosace/data");
    NXInt dataRosace = dataGroupRosace.openIntData();
    numberOfTubesInRosace += static_cast<size_t>(dataRosace.dim0());
  }

  // dim0 * m_numberOfPixelsPerTube is the total number of detectors
  m_numberOfHistograms =
      (m_numberOfTubes + numberOfTubesInRosace) * m_numberOfPixelsPerTube;

  g_log.debug() << "NumberOfTubes: " << m_numberOfTubes << '\n';
  g_log.debug() << "NumberOfPixelsPerTube: " << m_numberOfPixelsPerTube << '\n';
  g_log.debug() << "NumberOfChannels: " << m_numberOfChannels << '\n';

  // Now create the output workspace
  // total number of spectra + number of monitors,
  // bin boundaries = m_numberOfChannels + 1
  // Z/time dimension
  m_localWorkspace = WorkspaceFactory::Instance().create(
      "Workspace2D", m_numberOfHistograms + numberOfMonitors,
      m_numberOfChannels + 1, m_numberOfChannels);
  m_localWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
  m_localWorkspace->setYUnitLabel("Counts");
}

/**
 * Load the time details from the nexus file.
 *
 * @param entry :: The Nexus entry
 */
void LoadILLTOF2::loadTimeDetails(NeXus::NXEntry &entry) {

  m_wavelength = entry.getFloat("wavelength");

  // Monitor can be monitor (IN5) or monitor1 (IN6)
  std::string monitorName;
  if (entry.containsGroup("monitor"))
    monitorName = "monitor";
  else if (entry.containsGroup("monitor1"))
    monitorName = "monitor1";
  else {
    std::string message("Cannot find monitor[1] in the Nexus file!");
    g_log.error(message);
    throw std::runtime_error(message);
  }

  NXFloat time_of_flight_data =
      entry.openNXFloat(monitorName + "/time_of_flight");
  time_of_flight_data.load();

  // The entry "monitor/time_of_flight", has 3 fields:
  // channel width , number of channels, Time of flight delay
  m_channelWidth = time_of_flight_data[0];
  m_timeOfFlightDelay = time_of_flight_data[2];

  g_log.debug("Nexus Data:");
  g_log.debug() << " ChannelWidth: " << m_channelWidth << '\n';
  g_log.debug() << " TimeOfFlightDealy: " << m_timeOfFlightDelay << '\n';
  g_log.debug() << " Wavelength: " << m_wavelength << '\n';
}

/**
 * Goes through all the fields of the NeXus file and adds them
 * as parameters in the workspace
 *
 * @param filename The NeXus file
 */
261
void LoadILLTOF2::addAllNexusFieldsAsProperties(const std::string &filename) {
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285

  API::Run &runDetails = m_localWorkspace->mutableRun();

  // Open NeXus file
  NXhandle nxfileID;
  NXstatus stat = NXopen(filename.c_str(), NXACC_READ, &nxfileID);

  g_log.debug() << "Starting parsing properties from : " << filename << '\n';
  if (stat == NX_ERROR) {
    g_log.debug() << "convertNexusToProperties: Error loading " << filename;
    throw Kernel::Exception::FileError("Unable to open File:", filename);
  }
  m_loader.addNexusFieldsToWsRun(nxfileID, runDetails);

  g_log.debug() << "End parsing properties from : " << filename << '\n';
}

/**
 * Calculates the incident energy from the wavelength and adds
 * it as sample log 'Ei'
 */
void LoadILLTOF2::addEnergyToRun() {

  API::Run &runDetails = m_localWorkspace->mutableRun();
286
  const double ei = m_loader.calculateEnergy(m_wavelength);
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
  runDetails.addProperty<double>("Ei", ei, true); // overwrite
}

/**
 * Adds facility info to the sample logs
 */
void LoadILLTOF2::addFacility() {
  API::Run &runDetails = m_localWorkspace->mutableRun();
  runDetails.addProperty("Facility", std::string("ILL"));
}

/**
 * Calculates and adds the pulse intervals for the run
 */
void LoadILLTOF2::addPulseInterval() {
  API::Run &runDetails = m_localWorkspace->mutableRun();
303
304
  double n_pulses = -1;
  double fermiChopperSpeed = -1;
305
306
307
308

  if (m_instrumentName == "IN4") {
    fermiChopperSpeed =
        runDetails.getPropertyAsSingleValue("FC.rotation_speed");
309
    const double bkgChopper1Speed =
310
        runDetails.getPropertyAsSingleValue("BC1.rotation_speed");
311
    const double bkgChopper2Speed =
312
313
314
315
316
317
318
319
320
321
322
        runDetails.getPropertyAsSingleValue("BC2.rotation_speed");

    if (std::abs(bkgChopper1Speed - bkgChopper2Speed) > 1) {
      throw std::invalid_argument(
          "Background choppers 1 and 2 have different speeds");
    }

    n_pulses = fermiChopperSpeed / bkgChopper1Speed / 4;
  } else if (m_instrumentName == "IN6") {
    fermiChopperSpeed =
        runDetails.getPropertyAsSingleValue("Fermi.rotation_speed");
323
    const double suppressorSpeed =
324
325
326
327
328
329
330
        runDetails.getPropertyAsSingleValue("Suppressor.rotation_speed");

    n_pulses = fermiChopperSpeed / suppressorSpeed;
  } else {
    return;
  }

331
  const double pulseInterval = 60.0 / (2 * fermiChopperSpeed) * n_pulses;
332
333
334
335
336
337
338
339
340
341
  runDetails.addProperty<double>("pulse_interval", pulseInterval);
}

/**
 * Loads all the spectra into the workspace, including that from the monitor
 *
 * @param entry The Nexus entry
 * @param monitors List of monitor data
 */
void LoadILLTOF2::loadDataIntoTheWorkSpace(
Antti Soininen's avatar
Antti Soininen committed
342
    NeXus::NXEntry &entry, const std::vector<std::vector<int>> &monitors) {
343
344
345
346
347
348
349
350
351
352
353

  g_log.debug() << "Loading data into the workspace...\n";
  // read in the data
  NXData dataGroup = entry.openNXData("data");
  NXInt data = dataGroup.openIntData();
  // load the counts from the file into memory
  data.load();

  // Put tof in an array
  auto &X0 = m_localWorkspace->mutableX(0);
  for (size_t i = 0; i < m_numberOfChannels + 1; ++i) {
Antti Soininen's avatar
Antti Soininen committed
354
355
    X0[i] = m_timeOfFlightDelay + m_channelWidth * static_cast<double>(i) -
            m_channelWidth / 2; // to make sure the bin centre is correct
356
357
358
359
360
361
362
  }

  // The binning for monitors is considered the same as for detectors
  size_t spec = 0;

  auto const &instrument = m_localWorkspace->getInstrument();

363
  const std::vector<detid_t> detectorIDs = instrument->getDetectorIDs(true);
364

365
  Progress progress(this, 0.0, 1.0, m_numberOfTubes * m_numberOfPixelsPerTube);
366

Antti Soininen's avatar
Antti Soininen committed
367
  loadSpectra(spec, m_numberOfTubes, detectorIDs, data, progress);
368

369
  g_log.debug() << "Loading detector data into the workspace: DONE!\n";
370
371
372
373
374
375
376
377
378
379
380
381
382
383

  /**
   * IN4 Rosace detectors are in a different NeXus entry
   */
  if (m_instrumentName == "IN4") {
    g_log.debug() << "Loading data into the workspace: IN4 Rosace!\n";
    // read in the data
    NXData dataGroupRosace =
        entry.openNXData("instrument/Detector_Rosace/data");
    NXInt dataRosace = dataGroupRosace.openIntData();
    auto numberOfTubes = static_cast<size_t>(dataRosace.dim0());
    // load the counts from the file into memory
    dataRosace.load();

384
    Progress progressRosace(this, 0.0, 1.0,
385
386
                            numberOfTubes * m_numberOfPixelsPerTube);

Antti Soininen's avatar
Antti Soininen committed
387
    loadSpectra(spec, numberOfTubes, detectorIDs, dataRosace, progressRosace);
388
  }
389
390
391
392
393
394
395
396
397
398

  const auto monitorIDs = instrument->getMonitors();

  for (size_t i = 0; i < monitors.size(); ++i) {
    const auto &monitor = monitors[i];
    m_localWorkspace->setHistogram(spec, m_localWorkspace->binEdges(0),
                                   Counts(monitor.begin(), monitor.end()));
    m_localWorkspace->getSpectrum(spec).setDetectorID(monitorIDs[i]);
    spec++;
  }
399
400
401
402
403
404
405
406
407
408
409
410
}

/**
 * Loops over all the pixels and loads the correct spectra. Called for each set
 * of detector types in the workspace.
 *
 * @param spec The current spectrum id
 * @param numberOfTubes The number of detector tubes in the workspace
 * @param detectorIDs A list of all of the detector IDs
 * @param data The NeXus data to load into the workspace
 * @param progress The progress monitor
 */
Antti Soininen's avatar
Antti Soininen committed
411
void LoadILLTOF2::loadSpectra(size_t &spec, const size_t numberOfTubes,
Antti Soininen's avatar
Antti Soininen committed
412
                              const std::vector<detid_t> &detectorIDs,
413
                              const NXInt &data, Progress &progress) {
414
415
416
417
418
419
  for (size_t i = 0; i < numberOfTubes; ++i) {
    for (size_t j = 0; j < m_numberOfPixelsPerTube; ++j) {
      int *data_p = &data(static_cast<int>(i), static_cast<int>(j), 0);
      m_localWorkspace->setHistogram(
          spec, m_localWorkspace->binEdges(0),
          Counts(data_p, data_p + m_numberOfChannels));
Antti Soininen's avatar
Antti Soininen committed
420
      m_localWorkspace->getSpectrum(spec).setDetectorID(detectorIDs[spec]);
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
      spec++;
      progress.report();
    }
  }
}

/**
 * Runs LoadInstrument to attach the instrument to the workspace
 */
void LoadILLTOF2::runLoadInstrument() {

  IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");

  // Execute the child algorithm. Catch and log any error, but don't stop.
  try {
    loadInst->setPropertyValue("InstrumentName", m_instrumentName);
    loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", m_localWorkspace);
    loadInst->setProperty("RewriteSpectraMap",
                          Mantid::Kernel::OptionalBool(false));
    loadInst->execute();
Antti Soininen's avatar
Antti Soininen committed
441
  } catch (...) {
442
443
444
445
446
447
    g_log.information("Cannot load the instrument definition.");
  }
}

} // namespace DataHandling
} // namespace Mantid