LoadILLTOF2.cpp 14 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
#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"
#include "MantidKernel/UnitFactory.h"

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
const std::vector<std::string> LoadILLTOF2::SUPPORTED_INSTRUMENTS = {"IN4", "IN5", "IN6"};

/**
 * 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") &&
      !descriptor.pathExists(
          "/entry0/dataSD") // This one is for LoadILLIndirect
      &&
      !descriptor.pathExists(
          "/entry0/instrument/VirtualChopper") // This one is for
                                               // LoadILLReflectometry
      ) {
    return 80;
  } else {
    return 0;
  }
}

LoadILLTOF2::LoadILLTOF2() : API::IFileLoader<Kernel::NexusDescriptor>() {
}

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

  declareProperty(make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
                                                   Direction::Output),
                  "The name to use for the output workspace");
}

/**
 * Executes the algorithm
 */
void LoadILLTOF2::exec() {
  // Retrieve filename
72
  const std::string filenameData = getPropertyValue("Filename");
73
74
75
76
77
78
79
80

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

  loadInstrumentDetails(dataFirstEntry);
  loadTimeDetails(dataFirstEntry);

81
  const auto monitors = getMonitorInfo(dataFirstEntry);
82
83
84
85
86
87
88
89
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
117
118
119
120
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
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178

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

  std::vector<std::vector<int>> monitorList;

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

  if (m_instrumentPath == "") {
    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
 */
void LoadILLTOF2::initWorkSpace(NeXus::NXEntry &entry,
                               const std::vector<std::vector<int>> &monitors) {

  // 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());
179
  const size_t numberOfMonitors = monitors.size();
180
181
182
183
184
185
186
187
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

  /**
   * 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
 */
253
void LoadILLTOF2::addAllNexusFieldsAsProperties(const std::string &filename) {
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277

  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();
278
  const double ei = m_loader.calculateEnergy(m_wavelength);
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
  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();
295
296
  double n_pulses = -1;
  double fermiChopperSpeed = -1;
297
298
299
300

  if (m_instrumentName == "IN4") {
    fermiChopperSpeed =
        runDetails.getPropertyAsSingleValue("FC.rotation_speed");
301
    const double bkgChopper1Speed =
302
        runDetails.getPropertyAsSingleValue("BC1.rotation_speed");
303
    const double bkgChopper2Speed =
304
305
306
307
308
309
310
311
312
313
314
        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");
315
    const double suppressorSpeed =
316
317
318
319
320
321
322
        runDetails.getPropertyAsSingleValue("Suppressor.rotation_speed");

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

323
  const double pulseInterval = 60.0 / (2 * fermiChopperSpeed) * n_pulses;
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
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
 */
void LoadILLTOF2::loadDataIntoTheWorkSpace(
    NeXus::NXEntry &entry, const std::vector<std::vector<int>> &monitors) {

  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) {
    X0[i] = m_timeOfFlightDelay +
            m_channelWidth *
                static_cast<double>(i) -
            m_channelWidth /
                2; // to make sure the bin centre is correct
  }

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

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

358
  const auto monitorIDs = instrument->getMonitors();
359
360
361
362
363
364
365
366

  for (const auto &monitor : monitors) {
    m_localWorkspace->setHistogram(spec, m_localWorkspace->binEdges(0),
                                   Counts(monitor.begin(), monitor.end()));
    m_localWorkspace->getSpectrum(spec).setDetectorID(monitorIDs[spec]);
    spec++;
  }

367
368
  const std::vector<detid_t> detectorIDs = instrument->getDetectorIDs(true);
  const size_t numberOfMonitors = monitors.size();
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408

  Progress progress(this, 0, 1, m_numberOfTubes * m_numberOfPixelsPerTube);

  loadSpectra(spec, numberOfMonitors, m_numberOfTubes, detectorIDs, data,
              progress);

  g_log.debug() << "Loading data into the workspace: DONE!\n";

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

    Progress progressRosace(this, 0, 1,
                            numberOfTubes * m_numberOfPixelsPerTube);

    loadSpectra(spec, numberOfMonitors, numberOfTubes, detectorIDs, dataRosace,
                progressRosace);
  }
}

/**
 * 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 numberOfMonitors The number of monitors in the workspace
 * @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
 */
409
410
411
void LoadILLTOF2::loadSpectra(size_t &spec, const size_t numberOfMonitors,
                             const size_t numberOfTubes,
                             const std::vector<detid_t> &detectorIDs, NXInt data,
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
                             Progress progress) {
  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));
      m_localWorkspace->getSpectrum(spec)
          .setDetectorID(detectorIDs[spec - numberOfMonitors]);
      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();
  } catch (...) {
    g_log.information("Cannot load the instrument definition.");
  }
}

} // namespace DataHandling
} // namespace Mantid