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
154
  size_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
359
360
      specToLoad.insert(
          specToLoad.end(),
361
362
          boost::counting_iterator<specnum_t>(m_spec_min),
          boost::counting_iterator<specnum_t>(m_spec_max));
363
364
      specToLoad.insert(specToLoad.end(), m_spec_list.begin(),
                        m_spec_list.end());
365
366
    } else {
      // Load all the spectra
367
      // Start from 1 to N+1 to be consistent with
368
      // the case where spectra are specified
369
      for (int i = 1; i < numDeadTimes / m_numberOfPeriods + 1; i++)
370
        specToLoad.emplace_back(i);
371
    }
372
373

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

380
    } else if (numDeadTimes % m_numberOfSpectra) {
381

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

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

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

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

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

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

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

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

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

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

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

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

459
    int numGroupingEntries = groupingData.dim0();
460

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

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

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

488
489
490
    } 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 "
491
492
                                 "every spectrum in every period",
                                 m_filename);
493

494
    } else {
495

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

516
        TableWorkspace_sptr table =
517
            createDetectorGroupingTable(specToLoad, grouping);
518

519
        if (table->rowCount() != 0)
520
521
          return table;

522
523
524
525
      } else if (numGroupingEntries == m_numberOfSpectra) {
        // Multiple periods - same grouping for each
        specToLoad.clear();
        for (int i = 1; i < m_numberOfSpectra + 1; i++) {
526
          specToLoad.emplace_back(i);
527
528
529
530
531
532
533
534
535
        }
        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;
536
537
538
      } else {
        // More complex case - grouping information for every period

539
        WorkspaceGroup_sptr tableGroup = std::make_shared<WorkspaceGroup>();
540

541
        for (int i = 0; i < m_numberOfPeriods; i++) {
542

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

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

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

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

          return tableGroup;
        }
566
567
      }
    }
568
  }
569
570
571
572
573
574
575
576
577
578
  // 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");
579
    auto dummyGrouping = std::make_shared<Grouping>();
580
581
582
583
584
585
    if (inst->getNumberDetectors() != 0) {
      dummyGrouping = groupLoader.getDummyGrouping();
    } else {
      // Make sure it uses the right number of detectors
      std::ostringstream oss;
      oss << "1-" << m_numberOfSpectra;
586
      dummyGrouping->groups.emplace_back(oss.str());
587
588
589
590
      dummyGrouping->description = "Dummy grouping";
      dummyGrouping->groupNames.emplace_back("all");
    }
    return dummyGrouping->toTable();
591
  }
592
593
594
595
}

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

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

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

  return deadTimeTable;
}

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

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

632
  std::map<int, std::vector<int>> groupingMap;
633

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

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

  return detectorGroupingTable;
}

/** Load in a single spectrum taken from a NeXus file
LamarMoore's avatar
LamarMoore committed
652
653
654
655
656
657
658
659
 *  @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
660
661
662
663
void LoadMuonNexus1::loadData(
    size_t hist, specnum_t &i, specnum_t specNo, MuonNexusReader &nxload,
    const int64_t lengthIn,
    const DataObjects::Workspace2D_sptr &localWorkspace) {
664
665
666
667
  // 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
668

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

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

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

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

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

695
  auto numSpectra = static_cast<int>(localWorkspace->getNumberHistograms());
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
736
737
  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
738
739
void LoadMuonNexus1::runLoadLog(
    const DataObjects::Workspace2D_sptr &localWorkspace) {
740
741
742
743
744
745
746
747
748
749
  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 &) {
750
    g_log.error("Unable to successfully run LoadMuonLog Child Algorithm");
751
  } catch (std::logic_error &) {
752
    g_log.error("Unable to successfully run LoadMuonLog Child Algorithm");
753
754
755
  }

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

  NXRoot root(m_filename);

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

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

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

784
  ISISRunLogs runLogs(run);
785
786
787
788
789
790
791
792
  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
793
794
void LoadMuonNexus1::addPeriodLog(
    const DataObjects::Workspace2D_sptr &localWorkspace, int64_t period) {
795
  auto &run = localWorkspace->mutableRun();
796
  ISISRunLogs runLogs(run);
797
  if (period == 0) {
798
    runLogs.addPeriodLogs(1, run);
799
  } else {
800
801
802
    run.removeLogData("period 1");
    runLogs.addPeriodLog(static_cast<int>(period) + 1, run);
  }
803
804
}

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

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

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

815
    try {
Raquel Alvarez Banos's avatar
Raquel Alvarez Banos committed
816

817
      handle.openPath("run/instrument/beam");
818
819
820
821
822
823
824
825
      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");
      }
826

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

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

834
    } catch (::NeXus::Exception &) {
835
      g_log.warning("Could not read number of good frames");
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
880
881

  } 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
882
883
884
885

  handle.close();
}

886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
/**
 * 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
905
      return 0;
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
  }

  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;
926
    }
927
928
929
930
  } catch (...) {
  }
  return 0;
}
931

932
} // namespace DataHandling
933
} // namespace Mantid