LoadILLTOF2.cpp 15.5 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

/**
 * Initialises the algorithm
 */
void LoadILLTOF2::init() {
Sam Jenkins's avatar
Sam Jenkins committed
68
69
70
  declareProperty(std::make_unique<FileProperty>("Filename", "",
                                                 FileProperty::Load, ".nxs"),
                  "File path of the Data file to load");
71

72
  declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
Sam Jenkins's avatar
Sam Jenkins committed
73
                                                        Direction::Output),
74
                  "The name to use for the output workspace");
75
  declareProperty("ConvertToTOF", false, Direction::Input);
76
77
78
79
80
81
82
}

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

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

  loadInstrumentDetails(dataFirstEntry);
  loadTimeDetails(dataFirstEntry);

93
  const auto monitors = getMonitorInfo(dataFirstEntry);
94
95
96
97
98
99
100
101

  initWorkSpace(dataFirstEntry, monitors);

  addAllNexusFieldsAsProperties(filenameData);
  addFacility();

  runLoadInstrument(); // just to get IDF contents

102
  loadDataIntoTheWorkSpace(dataFirstEntry, monitors, convertToTOF);
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

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

Antti Soininen's avatar
Antti Soininen committed
124
  std::vector<std::vector<int>> monitorList;
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

  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());
141
      monitorList.emplace_back(thisMonitor);
142
143
144
145
146
147
148
149
150
151
152
153
154
155
    }
  }
  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);

156
  if (m_instrumentPath.empty()) {
157
158
159
160
161
162
163
164
165
166
167
168
169
170
    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);
  }

171
172
173
174
175
176
177
178
179
180
  // Monitor can be monitor (IN5, PANTHER) or monitor1 (IN6)
  if (firstEntry.containsGroup("monitor"))
    m_monitorName = "monitor";
  else if (firstEntry.containsGroup("monitor1"))
    m_monitorName = "monitor1";
  else {
    std::string message("Cannot find monitor[1] in the Nexus file!");
    throw std::runtime_error(message);
  }

181
182
183
184
185
186
187
188
189
190
  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
191
192
void LoadILLTOF2::initWorkSpace(NeXus::NXEntry &entry,
                                const std::vector<std::vector<int>> &monitors) {
193
194
195
196
197
198
199
200

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

  /**
   * 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);
229
230
231

  NXClass monitor = entry.openNXGroup(m_monitorName);
  if (monitor.containsDataSet("time_of_flight")) {
232
233
234
235
    m_localWorkspace->getAxis(0)->unit() =
        UnitFactory::Instance().create("TOF");
    m_localWorkspace->setYUnitLabel("Counts");
  } else {
236
    g_log.debug("PANTHER diffraction mode");
237
238
239
240
    m_localWorkspace->getAxis(0)->unit() =
        UnitFactory::Instance().create("Wavelength");
    m_localWorkspace->setYUnitLabel("Counts");
  }
241
242
243
244
245
246
247
248
249
250
251
}

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

  m_wavelength = entry.getFloat("wavelength");

252
  NeXus::NXClass monitorEntry = entry.openNXGroup(m_monitorName);
253

254
  if (monitorEntry.containsDataSet("time_of_flight")) {
255
256

    NXFloat time_of_flight_data =
257
        entry.openNXFloat(m_monitorName + "/time_of_flight");
258
    time_of_flight_data.load();
259

260
261
262
263
    // 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];
264

265
266
    g_log.debug("Nexus Data:");
    g_log.debug() << " ChannelWidth: " << m_channelWidth << '\n';
267
    g_log.debug() << " TimeOfFlightDelay: " << m_timeOfFlightDelay << '\n';
268
    g_log.debug() << " Wavelength: " << m_wavelength << '\n';
269
270
  } // the other case is the diffraction mode for PANTHER, where nothing is
    // needed here
271
272
273
274
275
276
277
278
}

/**
 * Goes through all the fields of the NeXus file and adds them
 * as parameters in the workspace
 *
 * @param filename The NeXus file
 */
279
void LoadILLTOF2::addAllNexusFieldsAsProperties(const std::string &filename) {
280
281
282
283
284
285
286
287
288
289
290
291
292

  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);
293
  NXclose(&nxfileID);
294
295
296
297
298
299
300
301
302
303
  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();
304
  const double ei = m_loader.calculateEnergy(m_wavelength);
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
  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();
321
322
  double n_pulses = -1;
  double fermiChopperSpeed = -1;
323
324
325
326

  if (m_instrumentName == "IN4") {
    fermiChopperSpeed =
        runDetails.getPropertyAsSingleValue("FC.rotation_speed");
327
    const double bkgChopper1Speed =
328
        runDetails.getPropertyAsSingleValue("BC1.rotation_speed");
329
    const double bkgChopper2Speed =
330
331
332
333
334
335
336
337
338
339
340
        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");
341
    const double suppressorSpeed =
342
343
344
345
346
347
348
        runDetails.getPropertyAsSingleValue("Suppressor.rotation_speed");

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

349
  const double pulseInterval = 60.0 / (2 * fermiChopperSpeed) * n_pulses;
350
351
352
353
354
355
356
357
  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
358
359
 * @param convertToTOF Should the bin edges be converted to time of flight or
 * keep the channel indexes
360
361
 */
void LoadILLTOF2::loadDataIntoTheWorkSpace(
362
363
    NeXus::NXEntry &entry, const std::vector<std::vector<int>> &monitors,
    bool convertToTOF) {
364
365
366
367
368
369
370
371

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

372
373
  NXClass monitor = entry.openNXGroup(m_monitorName);

374
375
  // Put tof in an array
  auto &X0 = m_localWorkspace->mutableX(0);
376
  if (monitor.containsDataSet("time_of_flight")) {
377
    if (convertToTOF) {
378
379
380
381
382
383
384
385
      for (size_t i = 0; i < m_numberOfChannels + 1; ++i) {
        X0[i] = m_timeOfFlightDelay + m_channelWidth * static_cast<double>(i) -
                m_channelWidth / 2; // to make sure the bin centre is correct
      }
    } else {
      for (size_t i = 0; i < m_numberOfChannels + 1; ++i) {
        X0[i] = static_cast<double>(i); // just take the channel index
      }
386
387
388
389
390
    }
  } else {
    // Diffraction PANTHER
    X0[0] = m_wavelength * 0.9;
    X0[1] = m_wavelength * 1.1;
391
392
393
394
395
396
  }
  // The binning for monitors is considered the same as for detectors
  size_t spec = 0;

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

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

399
  Progress progress(this, 0.0, 1.0, m_numberOfTubes * m_numberOfPixelsPerTube);
400

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

403
  g_log.debug() << "Loading detector data into the workspace: DONE!\n";
404
405
406
407
408
409
410
411
412
413
414
415
416
417

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

418
    Progress progressRosace(this, 0.0, 1.0,
419
420
                            numberOfTubes * m_numberOfPixelsPerTube);

Antti Soininen's avatar
Antti Soininen committed
421
    loadSpectra(spec, numberOfTubes, detectorIDs, dataRosace, progressRosace);
422
  }
423
424
425
426
427
428
429
430
431
432

  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++;
  }
433
434
435
436
437
438
439
440
441
442
443
444
}

/**
 * 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
445
void LoadILLTOF2::loadSpectra(size_t &spec, const size_t numberOfTubes,
Antti Soininen's avatar
Antti Soininen committed
446
                              const std::vector<detid_t> &detectorIDs,
447
                              const NXInt &data, Progress &progress) {
448
449
  PARALLEL_FOR_IF(Kernel::threadSafe(*m_localWorkspace))
  for (int i = 0; i < static_cast<int>(numberOfTubes); ++i) {
450
451
    for (size_t j = 0; j < m_numberOfPixelsPerTube; ++j) {
      int *data_p = &data(static_cast<int>(i), static_cast<int>(j), 0);
452
      const size_t currentSpectrum = spec + i * m_numberOfPixelsPerTube + j;
453
      m_localWorkspace->setHistogram(
454
          currentSpectrum, m_localWorkspace->binEdges(0),
455
          Counts(data_p, data_p + m_numberOfChannels));
456
457
      m_localWorkspace->getSpectrum(currentSpectrum)
          .setDetectorID(detectorIDs[currentSpectrum]);
458
459
460
      progress.report();
    }
  }
461
  spec += numberOfTubes * m_numberOfPixelsPerTube;
462
463
464
465
466
467
468
469
470
}

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

  IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");

471
472
473
474
475
  loadInst->setPropertyValue("InstrumentName", m_instrumentName);
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", m_localWorkspace);
  loadInst->setProperty("RewriteSpectraMap",
                        Mantid::Kernel::OptionalBool(false));
  loadInst->execute();
476
477
478
479
}

} // namespace DataHandling
} // namespace Mantid