LoadILLTOF2.cpp 14.1 KB
Newer Older
1
2
// Mantid Repository : https://github.com/mantidproject/mantid
//
3
// Copyright © 2019 ISIS Rutherford Appleton Laboratory UKRI,
4
5
6
//     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
namespace {
/// An array containing the supported instrument names
const std::array<std::string, 4> SUPPORTED_INSTRUMENTS = {
    {"IN4", "IN5", "IN6", "PANTHER"}};
} // namespace

26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
namespace Mantid {
namespace DataHandling {

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

DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadILLTOF2)

/**
 * 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
50
51
      !descriptor.pathExists("/entry0/dataSD") // This one is for
                                               // LoadILLIndirect
LamarMoore's avatar
LamarMoore committed
52
53
54
55
      && !descriptor.pathExists(
             "/entry0/instrument/VirtualChopper") // This one is for
                                                  // LoadILLReflectometry
  ) {
56
57
58
59
60
61
    return 80;
  } else {
    return 0;
  }
}

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

/**
 * 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
72
73
  declareProperty(make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
                                                   Direction::Output),
74
75
76
77
78
79
80
81
                  "The name to use for the output workspace");
}

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

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

  loadInstrumentDetails(dataFirstEntry);
  loadTimeDetails(dataFirstEntry);

91
  const auto monitors = getMonitorInfo(dataFirstEntry);
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
117
118

  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
119
std::vector<std::vector<int>>
120
121
LoadILLTOF2::getMonitorInfo(NeXus::NXEntry &firstEntry) {

Antti Soininen's avatar
Antti Soininen committed
122
  std::vector<std::vector<int>> monitorList;
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
152
153

  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);

154
  if (m_instrumentPath.empty()) {
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
    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
179
180
void LoadILLTOF2::initWorkSpace(NeXus::NXEntry &entry,
                                const std::vector<std::vector<int>> &monitors) {
181
182
183
184
185
186
187
188

  // 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());
189
  const size_t numberOfMonitors = monitors.size();
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
261
262

  /**
   * 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
 */
263
void LoadILLTOF2::addAllNexusFieldsAsProperties(const std::string &filename) {
264
265
266
267
268
269
270
271
272
273
274
275
276

  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);
277
  stat = NXclose(&nxfileID);
278
279
280
281
282
283
284
285
286
287
  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();
288
  const double ei = m_loader.calculateEnergy(m_wavelength);
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
  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();
305
306
  double n_pulses = -1;
  double fermiChopperSpeed = -1;
307
308
309
310

  if (m_instrumentName == "IN4") {
    fermiChopperSpeed =
        runDetails.getPropertyAsSingleValue("FC.rotation_speed");
311
    const double bkgChopper1Speed =
312
        runDetails.getPropertyAsSingleValue("BC1.rotation_speed");
313
    const double bkgChopper2Speed =
314
315
316
317
318
319
320
321
322
323
324
        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");
325
    const double suppressorSpeed =
326
327
328
329
330
331
332
        runDetails.getPropertyAsSingleValue("Suppressor.rotation_speed");

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

333
  const double pulseInterval = 60.0 / (2 * fermiChopperSpeed) * n_pulses;
334
335
336
337
338
339
340
341
342
343
  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
344
    NeXus::NXEntry &entry, const std::vector<std::vector<int>> &monitors) {
345
346
347
348
349
350
351
352
353
354
355

  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
356
357
    X0[i] = m_timeOfFlightDelay + m_channelWidth * static_cast<double>(i) -
            m_channelWidth / 2; // to make sure the bin centre is correct
358
359
360
361
362
363
364
  }

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

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

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

367
  Progress progress(this, 0.0, 1.0, m_numberOfTubes * m_numberOfPixelsPerTube);
368

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

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

  /**
   * 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();

386
    Progress progressRosace(this, 0.0, 1.0,
387
388
                            numberOfTubes * m_numberOfPixelsPerTube);

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

  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++;
  }
401
402
403
404
405
406
407
408
409
410
411
412
}

/**
 * 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
413
void LoadILLTOF2::loadSpectra(size_t &spec, const size_t numberOfTubes,
Antti Soininen's avatar
Antti Soininen committed
414
                              const std::vector<detid_t> &detectorIDs,
415
                              const NXInt &data, Progress &progress) {
416
417
418
419
420
421
  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
422
      m_localWorkspace->getSpectrum(spec).setDetectorID(detectorIDs[spec]);
423
424
425
426
427
428
429
430
431
432
433
434
435
      spec++;
      progress.report();
    }
  }
}

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

  IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");

436
437
438
439
440
  loadInst->setPropertyValue("InstrumentName", m_instrumentName);
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", m_localWorkspace);
  loadInst->setProperty("RewriteSpectraMap",
                        Mantid::Kernel::OptionalBool(false));
  loadInst->execute();
441
442
443
444
}

} // namespace DataHandling
} // namespace Mantid