LoadMuonNexus1.cpp 31.4 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
8
#include "MantidDataHandling/LoadMuonNexus1.h"

9
#include "MantidAPI/Axis.h"
10
#include "MantidAPI/FileProperty.h"
11
#include "MantidAPI/GroupingLoader.h"
12
#include "MantidAPI/Progress.h"
13
#include "MantidAPI/RegisterFileLoader.h"
14
#include "MantidAPI/SpectrumDetectorMapping.h"
15
#include "MantidAPI/TableRow.h"
16
#include "MantidAPI/WorkspaceFactory.h"
17
#include "MantidAPI/WorkspaceGroup.h"
18
#include "MantidDataHandling/ISISRunLogs.h"
19
20
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidDataObjects/Workspace2D.h"
21
22
23
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
24
#include "MantidKernel/ConfigService.h"
25
#include "MantidKernel/ListValidator.h"
26
#include "MantidKernel/TimeSeriesProperty.h"
27
#include "MantidKernel/Unit.h"
28
#include "MantidKernel/UnitFactory.h"
29
#include "MantidKernel/UnitLabelTypes.h"
30
31
#include "MantidNexus/MuonNexusReader.h"
#include "MantidNexus/NexusClasses.h"
32
// clang-format off
Raquel Alvarez Banos's avatar
Raquel Alvarez Banos committed
33
#include <nexus/NeXusFile.hpp>
34
#include <nexus/NeXusException.hpp>
35
// clang-format on
36

37
#include <boost/iterator/counting_iterator.hpp>
38
#include <memory>
Raquel Alvarez Banos's avatar
Raquel Alvarez Banos committed
39
#include <boost/scoped_array.hpp>
40

41
42
43
44
45
#include <Poco/Path.h>

#include <cmath>
#include <limits>

46
47
48
49
50
namespace Mantid {
namespace DataHandling {
using namespace DataObjects;

// Register the algorithm into the algorithm factory
51
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadMuonNexus1)
52
53
54
55

using namespace Kernel;
using namespace API;
using namespace Mantid::NeXus;
LamarMoore's avatar
LamarMoore committed
56
57
using HistogramData::BinEdges;
using HistogramData::Counts;
58
59
60
61
62

/// Empty default constructor
LoadMuonNexus1::LoadMuonNexus1() : LoadMuonNexus() {}

/** Executes the algorithm. Reading in the file and creating and populating
LamarMoore's avatar
LamarMoore committed
63
64
65
66
67
68
 *  the output workspace
 *
 *  @throw Exception::FileError If the Nexus file cannot be found/opened
 *  @throw std::invalid_argument If the optional properties are set to invalid
 *values
 */
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
void LoadMuonNexus1::exec() {
  // Retrieve the filename from the properties
  m_filename = getPropertyValue("Filename");
  // Retrieve the entry number
  m_entrynumber = getProperty("EntryNumber");

  NXRoot root(m_filename);
  NXEntry entry = root.openEntry("run/histogram_data_1");
  try {
    NXInfo info = entry.getDataSetInfo("time_zero");
    if (info.stat != NX_ERROR) {
      double dum = root.getFloat("run/histogram_data_1/time_zero");
      setProperty("TimeZero", dum);
    }
  } catch (...) {
  }

  try {
    NXInfo infoResolution = entry.getDataSetInfo("resolution");
    NXInt counts = root.openNXInt("run/histogram_data_1/counts");
    std::string firstGoodBin = counts.attributes("first_good_bin");
    if (!firstGoodBin.empty() && infoResolution.stat != NX_ERROR) {
      double resolution;

      switch (infoResolution.type) {
      case NX_FLOAT32:
        resolution = static_cast<double>(entry.getFloat("resolution"));
        break;
      case NX_INT32:
        resolution = static_cast<double>(entry.getInt("resolution"));
        break;
      default:
        throw std::runtime_error("Unsupported data type for resolution");
      }

104
      auto bin = static_cast<double>(boost::lexical_cast<int>(firstGoodBin));
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
      double bin_size = resolution / 1000000.0;

      setProperty("FirstGoodData", bin * bin_size);
    }
  } catch (std::exception &e) {
    g_log.warning() << "Error while loading the FirstGoodData value: "
                    << e.what() << "\n";
  }

  NXEntry nxRun = root.openEntry("run");
  std::string title;
  std::string notes;
  try {
    title = nxRun.getString("title");
    notes = nxRun.getString("notes");
  } catch (...) {
  }
  std::string run_num;
  try {
124
    run_num = std::to_string(nxRun.getInt("number"));
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
  } catch (...) {
  }

  MuonNexusReader nxload;
  nxload.readFromFile(m_filename);

  // Read in the instrument name from the Nexus file
  m_instrument_name = nxload.getInstrumentName();
  // Read in the number of spectra in the Nexus file
  m_numberOfSpectra = nxload.t_nsp1;
  if (m_entrynumber != 0) {
    m_numberOfPeriods = 1;
    if (m_entrynumber > nxload.t_nper) {
      throw std::invalid_argument("Invalid Entry Number:Enter a valid number");
    }
  } else {
    // Read the number of periods in this file
    m_numberOfPeriods = nxload.t_nper;
  }
144

145
  bool autoGroup = getProperty("AutoGroup");
146

147
148
  // Grouping info should be returned if user has set the property
  bool returnGrouping = !getPropertyValue("DetectorGroupingTable").empty();
149

150
151
  // Call private method to validate the optional parameters, if set
  checkOptionalProperties();
152
153
  // Calculate the size of a workspace, given its number of periods & spectra to
  // read
Stephen's avatar
Stephen committed
154
  int64_t total_specs;
155
156
  if (m_interval || m_list) {
    // Remove from list possible duplicate specs
157
158
    for (auto it = m_spec_list.begin(); it != m_spec_list.end();) {
      if ((*it >= m_spec_min) && (*it <= m_spec_max)) {
159
160
161
162
163
164
165
166
167
168
169
170
171
172
        it = m_spec_list.erase(it);
      } else {
        ++it;
      }
    }
    total_specs = m_spec_list.size();
    if (m_interval) {
      total_specs += (m_spec_max - m_spec_min + 1);
      m_spec_max += 1;
    }
  } else {
    total_specs = m_numberOfSpectra;
    // for nexus return all spectra
    m_spec_min = 1;
173
    m_spec_max = m_numberOfSpectra + 1; // Add +1 to iterate
174
175
  }

176
177
178
179
180
  // Read the number of time channels (i.e. bins) from the Nexus file
  const int channelsPerSpectrum = nxload.t_ntc1;
  // Read in the time bin boundaries
  const int lengthIn = channelsPerSpectrum + 1;

181
182
183
  // Try to load dead time info
  loadDeadTimes(root);

184
185
  // Create the 2D workspace for the output
  DataObjects::Workspace2D_sptr localWorkspace =
186
      std::dynamic_pointer_cast<DataObjects::Workspace2D>(
187
188
189
190
191
192
          WorkspaceFactory::Instance().create("Workspace2D", total_specs,
                                              lengthIn, lengthIn - 1));
  localWorkspace->setTitle(title);
  localWorkspace->setComment(notes);
  localWorkspace->mutableRun().addLogData(
      new PropertyWithValue<std::string>("run_number", run_num));
193
194

  // Add 'FirstGoodData' to list of logs if possible
195
  if (existsProperty("FirstGoodData") && existsProperty("TimeZero")) {
196
197
    double fgd = getProperty("FirstGoodData");
    double tz = getProperty("TimeZero");
198
199
    localWorkspace->mutableRun().addLogData(
        new PropertyWithValue<double>("FirstGoodData", fgd - tz));
200
  }
201
202
  // Set the unit on the workspace to muon time, for now in the form of a Label
  // Unit
203
204
  std::shared_ptr<Kernel::Units::Label> lblUnit =
      std::dynamic_pointer_cast<Kernel::Units::Label>(
205
206
207
208
209
210
211
212
          UnitFactory::Instance().create("Label"));
  lblUnit->setLabel("Time", Units::Symbol::Microsecond);
  localWorkspace->getAxis(0)->unit() = lblUnit;
  // Set y axis unit
  localWorkspace->setYUnit("Counts");

  WorkspaceGroup_sptr wsGrpSptr = WorkspaceGroup_sptr(new WorkspaceGroup);

213
  API::Progress progress(this, 0.0, 1.0, m_numberOfPeriods * total_specs);
214
215
216
217
218
219
220
221
  // Loop over the number of periods in the Nexus file, putting each period in a
  // separate workspace
  for (int64_t period = 0; period < m_numberOfPeriods; ++period) {
    if (m_entrynumber != 0) {
      period = m_entrynumber - 1;
      if (period != 0) {
        loadRunDetails(localWorkspace);
        runLoadInstrument(localWorkspace);
222
      }
223
    }
224

225
226
227
228
229
230
231
232
    if (period == 0) {
      // Only run the Child Algorithms once
      loadRunDetails(localWorkspace);
      runLoadInstrument(localWorkspace);
      runLoadLog(localWorkspace);
      localWorkspace->populateInstrumentParameters();
    } else // We are working on a higher period of a multiperiod raw file
    {
233
      localWorkspace = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
234
          WorkspaceFactory::Instance().create(localWorkspace));
235
236
237
      localWorkspace->setTitle(title);
      localWorkspace->setComment(notes);
    }
238
239
    addPeriodLog(localWorkspace, period);
    addGoodFrames(localWorkspace, period, nxload.t_nper);
240

241
    size_t counter = 0;
242
    for (auto i = m_spec_min; i < m_spec_max; ++i) {
243
      // Shift the histogram to read if we're not in the first period
244
      auto histToRead = static_cast<specnum_t>(i - 1 + period * nxload.t_nsp1);
245
      auto specNo = i;
246
      loadData(counter, histToRead, specNo, nxload, lengthIn - 1,
247
248
249
250
251
252
               localWorkspace); // added -1 for NeXus
      counter++;
      progress.report();
    }
    // Read in the spectra in the optional list parameter, if set
    if (m_list) {
253
      for (auto specid : m_spec_list) {
254
        auto histToRead =
255
            static_cast<specnum_t>(specid - 1 + period * nxload.t_nsp1);
256
        auto specNo = static_cast<specnum_t>(specid);
257
        loadData(counter, histToRead, specNo, nxload, lengthIn - 1,
258
                 localWorkspace);
259
260
        counter++;
        progress.report();
261
262
      }
    }
263
264
    // Just a sanity check
    assert(counter == size_t(total_specs));
265

266
    Workspace_sptr outWs;
267

268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
    Workspace_sptr loadedGrouping;

    // Try to load detector grouping info, if needed for auto-grouping or user
    // requested it
    if (autoGroup || returnGrouping) {
      loadedGrouping =
          loadDetectorGrouping(root, localWorkspace->getInstrument());

      if (loadedGrouping && returnGrouping) {
        // Return loaded grouping, if requested
        setProperty("DetectorGroupingTable", loadedGrouping);
      }

      if (!loadedGrouping && autoGroup) {
        // If autoGroup requested and no grouping in the file - show a warning
        g_log.warning(
            "Unable to load grouping from the file. Grouping not applied.");
      }
    }

288
289
    if (autoGroup && loadedGrouping) {
      TableWorkspace_sptr groupingTable;
290

291
      if (auto table =
292
              std::dynamic_pointer_cast<TableWorkspace>(loadedGrouping)) {
293
        groupingTable = table;
294
      } else if (auto group = std::dynamic_pointer_cast<WorkspaceGroup>(
295
296
                     loadedGrouping)) {
        groupingTable =
297
            std::dynamic_pointer_cast<TableWorkspace>(group->getItem(period));
298
      }
299
      std::vector<int> specIDs, detecIDs;
300
      for (size_t i = 0; i < localWorkspace->getNumberHistograms(); i++) {
301
302
        specIDs.emplace_back(localWorkspace->getSpectrum(i).getSpectrumNo());
        detecIDs.emplace_back(localWorkspace->getSpectrum(i).getSpectrumNo());
303
      }
304
      API::SpectrumDetectorMapping mapping(specIDs, detecIDs);
305
      localWorkspace->updateSpectraUsing(mapping);
306

307
308
309
310
      Algorithm_sptr groupDet = createChildAlgorithm("MuonGroupDetectors");
      groupDet->setProperty("InputWorkspace", localWorkspace);
      groupDet->setProperty("DetectorGroupingTable", groupingTable);
      groupDet->execute();
311

312
      MatrixWorkspace_sptr groupedWs = groupDet->getProperty("OutputWorkspace");
313

314
315
316
      outWs = groupedWs;
    } else {
      outWs = localWorkspace;
317
318
    }

319
320
321
322
323
324
325
326
327
328
329
330
    if (m_numberOfPeriods == 1)
      setProperty("OutputWorkspace", outWs);
    else
      // In case of multiple periods, just add workspace to the group, and we
      // will return the
      // group later
      wsGrpSptr->addWorkspace(outWs);

  } // loop over periods

  if (m_numberOfPeriods > 1) {
    setProperty("OutputWorkspace",
331
                std::dynamic_pointer_cast<Workspace>(wsGrpSptr));
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
  }
}

/**
 * Loads dead time table for the detector.
 * @param root :: Root entry of the Nexus to read dead times from
 */
void LoadMuonNexus1::loadDeadTimes(NXRoot &root) {
  // If dead times workspace name is empty - caller doesn't need dead times
  if (getPropertyValue("DeadTimeTable").empty())
    return;

  NXEntry detector = root.openEntry("run/instrument/detector");

  NXInfo infoDeadTimes = detector.getDataSetInfo("deadtimes");
  if (infoDeadTimes.stat != NX_ERROR) {
    NXFloat deadTimesData = detector.openNXFloat("deadtimes");
    deadTimesData.load();

    int numDeadTimes = deadTimesData.dim0();

353
    std::vector<int> specToLoad;
354
    std::vector<double> deadTimes;
355
356

    // Set the spectrum list that should be loaded
357
    if (m_interval || m_list) {
358
      // Load only selected spectra
Stephen's avatar
Stephen committed
359
360
361
      specToLoad.insert(specToLoad.end(),
                        boost::counting_iterator<specnum_t>(m_spec_min),
                        boost::counting_iterator<specnum_t>(m_spec_max));
362
363
      specToLoad.insert(specToLoad.end(), m_spec_list.begin(),
                        m_spec_list.end());
364
365
    } else {
      // Load all the spectra
366
      // Start from 1 to N+1 to be consistent with
367
      // the case where spectra are specified
368
      for (int i = 1; i < numDeadTimes / m_numberOfPeriods + 1; i++)
369
        specToLoad.emplace_back(i);
370
    }
371
372

    if (numDeadTimes < m_numberOfSpectra) {
373
      // Check number of dead time entries match the number of
374
      // spectra in the nexus file
375
      throw Exception::FileError(
376
377
          "Number of dead times specified is less than number of spectra",
          m_filename);
378

379
    } else if (numDeadTimes % m_numberOfSpectra) {
380

381
382
      // At least, number of dead times should cover the number of spectra
      throw Exception::FileError(
383
384
          "Number of dead times doesn't cover every spectrum in every period",
          m_filename);
385
    } else {
386

387
      if (m_numberOfPeriods == 1) {
388
        // Simplest case - one dead time for one detector
389
390

        // Populate deadTimes
391
        deadTimes.reserve(specToLoad.size());
392
393
        for (auto &spectra : specToLoad) {
          deadTimes.emplace_back(deadTimesData[spectra - 1]);
394
395
        }
        // Load into table
396
397
        TableWorkspace_sptr table = createDeadTimeTable(specToLoad, deadTimes);
        setProperty("DeadTimeTable", table);
398

399
400
401
402
403
      } else if (numDeadTimes == m_numberOfSpectra) {
        // Multiple periods, but the same dead times for each

        specToLoad.clear();
        for (int i = 1; i < numDeadTimes + 1; i++) {
404
          specToLoad.emplace_back(i);
405
        }
406
        deadTimes.reserve(specToLoad.size());
407
408
409
410
411
412
        for (const auto &spectrum : specToLoad) {
          deadTimes.emplace_back(deadTimesData[spectrum - 1]);
        }
        // Load into table
        TableWorkspace_sptr table = createDeadTimeTable(specToLoad, deadTimes);
        setProperty("DeadTimeTable", table);
413
414
      } else {
        // More complex case - different dead times for different periods
415
        WorkspaceGroup_sptr tableGroup = std::make_shared<WorkspaceGroup>();
416

417
        for (int64_t i = 0; i < m_numberOfPeriods; i++) {
418
419

          // Populate deadTimes
420
          deadTimes.reserve(specToLoad.size());
421
          for (auto spec : specToLoad) {
422
            auto index = static_cast<int>(spec - 1 + i * m_numberOfSpectra);
423
            deadTimes.emplace_back(deadTimesData[index]);
424
425
426
          }

          // Load into table
427
428
          TableWorkspace_sptr table =
              createDeadTimeTable(specToLoad, deadTimes);
429

430
431
432
433
434
          tableGroup->addWorkspace(table);
        }

        setProperty("DeadTimeTable", tableGroup);
      }
435
    }
436
437
438
439
440
  }
  // It is expected that file might not contain any dead times, so not finding
  // them is not an
  // error
}
441

442
443
/**
 * Loads detector grouping.
444
 * If no entry in NeXus file for grouping, load it from the IDF.
445
 * @param root :: Root entry of the Nexus file to read from
446
447
 * @param inst :: Pointer to instrument (to use if IDF needed)
 * @returns :: Grouping table - or tables, if per period
448
 */
David Fairbrother's avatar
David Fairbrother committed
449
450
Workspace_sptr LoadMuonNexus1::loadDetectorGrouping(
    NXRoot &root, const Geometry::Instrument_const_sptr &inst) {
451
  NXEntry dataEntry = root.openEntry("run/histogram_data_1");
452

453
454
455
456
  NXInfo infoGrouping = dataEntry.getDataSetInfo("grouping");
  if (infoGrouping.stat != NX_ERROR) {
    NXInt groupingData = dataEntry.openNXInt("grouping");
    groupingData.load();
457

458
    int numGroupingEntries = groupingData.dim0();
459

460
    std::vector<int> specToLoad;
461
    std::vector<int> grouping;
462
463

    // Set the spectrum list that should be loaded
464
    if (m_interval || m_list) {
465
      // Load only selected spectra
Stephen's avatar
Stephen committed
466
467
468
      specToLoad.insert(specToLoad.end(),
                        boost::counting_iterator<specnum_t>(m_spec_min),
                        boost::counting_iterator<specnum_t>(m_spec_max));
469
470
      specToLoad.insert(specToLoad.end(), m_spec_list.begin(),
                        m_spec_list.end());
471
472
    } else {
      // Load all the spectra
473
      // Start from 1 to N+1 to be consistent with
474
      // the case where spectra are specified
475
      for (int i = 1; i < m_numberOfSpectra + 1; i++)
476
        specToLoad.emplace_back(i);
477
    }
478

479
    if (numGroupingEntries < m_numberOfSpectra) {
480
      // Check number of dead time entries match the number of
481
      // spectra in the nexus file
482
483
484
      throw Exception::FileError(
          "Number of grouping entries is less than number of spectra",
          m_filename);
485

486
487
488
    } else if (numGroupingEntries % m_numberOfSpectra) {
      // At least the number of entries should cover all the spectra
      throw Exception::FileError("Number of grouping entries doesn't cover "
489
490
                                 "every spectrum in every period",
                                 m_filename);
491

492
    } else {
493

494
      if (m_numberOfPeriods == 1) {
495
        // Simplest case - one grouping entry per spectrum
496
        grouping.reserve(grouping.size() + specToLoad.size());
497
498
499
        if (!m_entrynumber) {
          // m_entrynumber = 0 && m_numberOfPeriods = 1 means that user did not
          // select
500
          // any periods but there is only one in the dataset
501
502
          for (auto spec : specToLoad) {
            grouping.emplace_back(groupingData[spec - 1]);
503
504
505
          }
        } else {
          // User selected an entry number
506
          for (auto &spec : specToLoad) {
LamarMoore's avatar
LamarMoore committed
507
508
509
            int index =
                spec - 1 +
                static_cast<int>((m_entrynumber - 1) * m_numberOfSpectra);
510
            grouping.emplace_back(groupingData[index]);
511
          }
512
513
        }

514
        TableWorkspace_sptr table =
515
            createDetectorGroupingTable(specToLoad, grouping);
516

517
        if (table->rowCount() != 0)
518
519
          return table;

520
521
522
523
      } else if (numGroupingEntries == m_numberOfSpectra) {
        // Multiple periods - same grouping for each
        specToLoad.clear();
        for (int i = 1; i < m_numberOfSpectra + 1; i++) {
524
          specToLoad.emplace_back(i);
525
526
527
528
529
530
531
532
533
        }
        for (const auto &spectrum : specToLoad) {
          grouping.emplace_back(groupingData[spectrum - 1]);
        }
        // Load into table
        TableWorkspace_sptr table =
            createDetectorGroupingTable(specToLoad, grouping);
        if (table->rowCount() != 0)
          return table;
534
535
536
      } else {
        // More complex case - grouping information for every period

537
        WorkspaceGroup_sptr tableGroup = std::make_shared<WorkspaceGroup>();
538

539
        for (int i = 0; i < m_numberOfPeriods; i++) {
540

541
          // Get the grouping
542
          grouping.clear();
543
544
545
          for (auto &spec : specToLoad) {
            int index = spec - 1 + i * static_cast<int>(m_numberOfSpectra);
            grouping.emplace_back(groupingData[index]);
546
547
548
549
          }

          // Set table for period i
          TableWorkspace_sptr table =
550
              createDetectorGroupingTable(specToLoad, grouping);
551
552
553
554

          // Add table to group
          if (table->rowCount() != 0)
            tableGroup->addWorkspace(table);
555
556
557
558
559
        }

        if (tableGroup->size() != 0) {
          if (tableGroup->size() != static_cast<size_t>(m_numberOfPeriods))
            throw Exception::FileError("Zero grouping for some of the periods",
560
                                       m_filename);
561
562
563

          return tableGroup;
        }
564
565
      }
    }
566
  }
567
568
569
570
571
572
573
574
575
576
  // If we reach this point, no/zero grouping found.
  // Try to load from IDF instead
  const std::string mainFieldDirection = getProperty("MainFieldDirection");
  API::GroupingLoader groupLoader(inst, mainFieldDirection);
  try {
    const auto idfGrouping = groupLoader.getGroupingFromIDF();
    g_log.warning("Loading grouping from IDF");
    return idfGrouping->toTable();
  } catch (const std::runtime_error &) {
    g_log.warning("Loading dummy grouping");
577
    auto dummyGrouping = std::make_shared<Grouping>();
578
579
580
581
582
583
    if (inst->getNumberDetectors() != 0) {
      dummyGrouping = groupLoader.getDummyGrouping();
    } else {
      // Make sure it uses the right number of detectors
      std::ostringstream oss;
      oss << "1-" << m_numberOfSpectra;
584
      dummyGrouping->groups.emplace_back(oss.str());
585
586
587
588
      dummyGrouping->description = "Dummy grouping";
      dummyGrouping->groupNames.emplace_back("all");
    }
    return dummyGrouping->toTable();
589
  }
590
591
592
593
}

/**
 * Creates Dead Time Table using all the data between begin and end.
594
595
 * @param specToLoad :: vector containing the spectrum numbers to load
 * @param deadTimes :: vector containing the corresponding dead times
596
597
598
 * @return Dead Time Table create using the data
 */
TableWorkspace_sptr
599
600
LoadMuonNexus1::createDeadTimeTable(std::vector<int> specToLoad,
                                    std::vector<double> deadTimes) {
601
602
  TableWorkspace_sptr deadTimeTable = std::dynamic_pointer_cast<TableWorkspace>(
      WorkspaceFactory::Instance().createTable("TableWorkspace"));
603
604
605
606

  deadTimeTable->addColumn("int", "spectrum");
  deadTimeTable->addColumn("double", "dead-time");

607
  for (size_t i = 0; i < specToLoad.size(); i++) {
608
    TableRow row = deadTimeTable->appendRow();
609
    row << specToLoad[i] << deadTimes[i];
610
611
612
613
614
615
616
617
  }

  return deadTimeTable;
}

/**
 * Creates Detector Grouping Table using all the data between begin and end.
 *
618
619
 * @param specToLoad :: Vector containing the spectrum list to load
 * @param grouping :: Vector containing corresponding grouping
620
621
 * @return Detector Grouping Table create using the data
 */
622
623
624
TableWorkspace_sptr
LoadMuonNexus1::createDetectorGroupingTable(std::vector<int> specToLoad,
                                            std::vector<int> grouping) {
625
  auto detectorGroupingTable = std::dynamic_pointer_cast<TableWorkspace>(
626
627
628
629
      WorkspaceFactory::Instance().createTable("TableWorkspace"));

  detectorGroupingTable->addColumn("vector_int", "Detectors");

630
  std::map<int, std::vector<int>> groupingMap;
631

632
  for (size_t i = 0; i < specToLoad.size(); i++) {
633
634
    // Add detector ID to the list of group detectors. Detector ID is always
    // spectra index + 1
635
    groupingMap[grouping[i]].emplace_back(specToLoad[i]);
636
637
  }

Hahn, Steven's avatar
Hahn, Steven committed
638
639
  for (auto &group : groupingMap) {
    if (group.first != 0) // Skip 0 group
640
    {
641
      TableRow newRow = detectorGroupingTable->appendRow();
Hahn, Steven's avatar
Hahn, Steven committed
642
      newRow << group.second;
643
644
645
646
647
648
649
    }
  }

  return detectorGroupingTable;
}

/** Load in a single spectrum taken from a NeXus file
LamarMoore's avatar
LamarMoore committed
650
651
652
653
654
655
656
657
 *  @param hist ::     The workspace index
 *  @param i ::        The spectrum number
 *  @param specNo ::   The spectrum number
 *  @param nxload ::   A reference to the MuonNeXusReader object
 *  @param lengthIn :: The number of elements in a spectrum
 *  @param localWorkspace :: A pointer to the workspace in which the data will
 * be stored
 */
David Fairbrother's avatar
David Fairbrother committed
658
659
660
661
void LoadMuonNexus1::loadData(
    size_t hist, specnum_t &i, specnum_t specNo, MuonNexusReader &nxload,
    const int64_t lengthIn,
    const DataObjects::Workspace2D_sptr &localWorkspace) {
662
663
664
665
  // Read in a spectrum
  // Put it into a vector, discarding the 1st entry, which is rubbish
  // But note that the last (overflow) bin is kept
  // For Nexus, not sure if above is the case, hence give all data for now
666

667
  // Create and fill another vector for the X axis
668
  auto timeChannels = new float[lengthIn + 1]();
669
  nxload.getTimeChannels(timeChannels, static_cast<int>(lengthIn + 1));
670
671
  // Put the read in array into a vector (inside a shared pointer)

LamarMoore's avatar
LamarMoore committed
672
673
674
675
676
  localWorkspace->setHistogram(
      hist, BinEdges(timeChannels, timeChannels + lengthIn + 1),
      Counts(nxload.counts + i * lengthIn,
             nxload.counts + i * lengthIn + lengthIn));

677
  localWorkspace->getSpectrum(hist).setSpectrumNo(specNo);
678
679
  // Muon v1 files: always a one-to-one mapping between spectra and detectors
  localWorkspace->getSpectrum(hist).setDetectorID(static_cast<detid_t>(specNo));
680
681
  // Clean up
  delete[] timeChannels;
682
683
684
}

/**  Log the run details from the file
LamarMoore's avatar
LamarMoore committed
685
686
 * @param localWorkspace :: The workspace details to use
 */
687
void LoadMuonNexus1::loadRunDetails(
David Fairbrother's avatar
David Fairbrother committed
688
    const DataObjects::Workspace2D_sptr &localWorkspace) {
689
690
691
692
  API::Run &runDetails = localWorkspace->mutableRun();

  runDetails.addProperty("run_title", localWorkspace->getTitle(), true);

693
  auto numSpectra = static_cast<int>(localWorkspace->getNumberHistograms());
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
  runDetails.addProperty("nspectra", numSpectra);

  NXRoot root(m_filename);
  try {
    std::string start_time = root.getString("run/start_time");
    runDetails.addProperty("run_start", start_time);
  } catch (std::runtime_error &) {
    g_log.warning("run/start_time is not available, run_start log not added.");
  }

  try {
    std::string stop_time = root.getString("run/stop_time");
    runDetails.addProperty("run_end", stop_time);
  } catch (std::runtime_error &) {
    g_log.warning("run/stop_time is not available, run_end log not added.");
  }

  try {
    std::string dur = root.getString("run/duration");
    runDetails.addProperty("dur", dur);
    runDetails.addProperty("durunits", 1); // 1 means second here
    runDetails.addProperty("dur_secs", dur);
  } catch (std::runtime_error &) {
    g_log.warning("run/duration is not available, dur log not added.");
  }

  // Get sample parameters
  NXEntry runSample = root.openEntry("run/sample");

  if (runSample.containsDataSet("temperature")) {
    float temperature = runSample.getFloat("temperature");
    runDetails.addProperty("sample_temp", static_cast<double>(temperature));
  }

  if (runSample.containsDataSet("magnetic_field")) {
    float magn_field = runSample.getFloat("magnetic_field");
    runDetails.addProperty("sample_magn_field",
                           static_cast<double>(magn_field));
  }
}

/// Run the LoadLog Child Algorithm
David Fairbrother's avatar
David Fairbrother committed
736
737
void LoadMuonNexus1::runLoadLog(
    const DataObjects::Workspace2D_sptr &localWorkspace) {
738
739
740
741
742
743
744
745
746
747
  IAlgorithm_sptr loadLog = createChildAlgorithm("LoadMuonLog");
  // Pass through the same input filename
  loadLog->setPropertyValue("Filename", m_filename);
  // Set the workspace property to be the same one filled above
  loadLog->setProperty<MatrixWorkspace_sptr>("Workspace", localWorkspace);

  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
  try {
    loadLog->execute();
  } catch (std::runtime_error &) {
748
    g_log.error("Unable to successfully run LoadMuonLog Child Algorithm");
749
  } catch (std::logic_error &) {
750
    g_log.error("Unable to successfully run LoadMuonLog Child Algorithm");
751
752
753
  }

  if (!loadLog->isExecuted())
754
    g_log.error("Unable to successfully run LoadMuonLog Child Algorithm");
755
756
757

  NXRoot root(m_filename);

758
759
  // Get main field direction
  std::string mainFieldDirection = "Longitudinal"; // default
760
761
762
763
764
765
  try {
    NXChar orientation = root.openNXChar("run/instrument/detector/orientation");
    // some files have no data there
    orientation.load();

    if (orientation[0] == 't') {
766
      auto p =
767
          std::make_unique<Kernel::TimeSeriesProperty<double>>("fromNexus");
768
769
      std::string start_time = root.getString("run/start_time");
      p->addValue(start_time, -90.0);
770
771
      localWorkspace->mutableRun().addLogData(std::move(p));
      mainFieldDirection = "Transverse";
772
773
    }
  } catch (...) {
774
    // no data - assume main field was longitudinal
775
  }
776

777
  // set output property and add to workspace logs
778
  auto &run = localWorkspace->mutableRun();
779
780
781
  setProperty("MainFieldDirection", mainFieldDirection);
  run.addProperty("main_field_direction", mainFieldDirection);

782
  ISISRunLogs runLogs(run);
783
784
785
786
787
788
789
790
  runLogs.addStatusLog(run);
}

/**
 * Add the 'period i' log to a workspace.
 * @param localWorkspace A workspace to add the log to.
 * @param period A period for this workspace.
 */
David Fairbrother's avatar
David Fairbrother committed
791
792
void LoadMuonNexus1::addPeriodLog(
    const DataObjects::Workspace2D_sptr &localWorkspace, int64_t period) {
793
  auto &run = localWorkspace->mutableRun();
794
  ISISRunLogs runLogs(run);
795
  if (period == 0) {
796
    runLogs.addPeriodLogs(1, run);
797
  } else {
798
799
800
    run.removeLogData("period 1");
    runLogs.addPeriodLog(static_cast<int>(period) + 1, run);
  }
801
802
}

David Fairbrother's avatar
David Fairbrother committed
803
804
805
void LoadMuonNexus1::addGoodFrames(
    const DataObjects::Workspace2D_sptr &localWorkspace, int64_t period,
    int nperiods) {
Raquel Alvarez Banos's avatar
Raquel Alvarez Banos committed
806
807
808
809

  // Get handle to nexus file
  ::NeXus::File handle(m_filename, NXACC_READ);

810
  // For single-period datasets, read /run/instrument/beam/frames_good
811
  if (nperiods == 1) {
Raquel Alvarez Banos's avatar
Raquel Alvarez Banos committed
812

813
    try {
Raquel Alvarez Banos's avatar
Raquel Alvarez Banos committed
814

815
      handle.openPath("run/instrument/beam");
816
817
818
819
820
821
822
823
      try {
        handle.openData("frames_good");
      } catch (::NeXus::Exception &) {
        // If it's not there, read "frames" instead and assume they are good
        g_log.warning("Could not read /run/instrument/beam/frames_good");
        handle.openData("frames");
        g_log.warning("Using run/instrument/beam/frames instead");
      }
824

825
826
827
828
829
      // read frames_period_daq
      boost::scoped_array<int> dataVals(new int[1]);
      handle.getData(dataVals.get());

      auto &run = localWorkspace->mutableRun();
830
831
      run.addProperty("goodfrm", dataVals[0]);

832
    } catch (::NeXus::Exception &) {
833
      g_log.warning("Could not read number of good frames");
834
    }
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879

  } else {
    // For multi-period datasets, read entries in
    // /run/instrument/beam/frames_period_daq
    try {

      handle.openPath("run/instrument/beam/");
      handle.openData("frames_period_daq");

      ::NeXus::Info info = handle.getInfo();
      // Check that frames_period_daq contains values for
      // every period
      if (period >= info.dims[0]) {
        std::ostringstream error;
        error << "goodfrm not found for period " << period;
        throw std::runtime_error(error.str());
      }
      if (nperiods != info.dims[0]) {
        std::ostringstream error;
        error << "Inconsistent number of period entries found (";
        error << info.dims[0];
        error << "!=" << nperiods << ")";
        throw std::runtime_error(error.str());
      }

      // read frames_period_daq
      boost::scoped_array<int> dataVals(new int[info.dims[0]]);
      handle.getData(dataVals.get());

      auto &run = localWorkspace->mutableRun();
      if (period == 0) {
        // If this is the first period
        // localWorkspace will not contain goodfrm
        run.addProperty("goodfrm", dataVals[0]);

      } else {
        // If period > 0, we need to remove
        // previous goodfrm log value
        run.removeLogData("goodfrm");
        run.addProperty("goodfrm", dataVals[period]);
      }
    } catch (::NeXus::Exception &) {
      g_log.warning("Could not read /run/instrument/beam/frames_period_daq");
    }
  } // else
Raquel Alvarez Banos's avatar
Raquel Alvarez Banos committed
880
881
882
883

  handle.close();
}

884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
/**
 * Return the confidence with with this algorithm can load the file
 * @param descriptor A descriptor for the file
 * @returns An integer specifying the confidence level. 0 indicates it will not
 * be used
 */
int LoadMuonNexus1::confidence(Kernel::NexusDescriptor &descriptor) const {
  const auto &firstEntryNameType = descriptor.firstEntryNameType();
  const std::string root = "/" + firstEntryNameType.first;
  if (!descriptor.pathExists(root + "/analysis"))
    return 0;

  bool upperIDF(true);
  if (descriptor.pathExists(root + "/IDF_version"))
    upperIDF = true;
  else {
    if (descriptor.pathExists(root + "/idf_version"))
      upperIDF = false;
    else
903
      return 0;
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
  }

  try {
    std::string versionField = "idf_version";
    if (upperIDF)
      versionField = "IDF_version";

    auto &file = descriptor.data();
    file.openPath(root + "/" + versionField);
    int32_t version = 0;
    file.getData(&version);
    if (version != 1)
      return 0;

    file.openPath(root + "/analysis");
    std::string def = file.getStrData();
    if (def == "muonTD" || def == "pulsedTD") {
      // If all this succeeded then we'll assume this is an ISIS Muon NeXus file
      // version 2
      return 81;
924
    }
925
926
927
928
  } catch (...) {
  }
  return 0;
}
929

930
} // namespace DataHandling
931
} // namespace Mantid