LoadNGEM.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2019 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidDataHandling/LoadNGEM.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/MultipleFileProperty.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/Run.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/OptionalBool.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/UnitFactory.h"
 
#include <boost/algorithm/string.hpp>
#include <fstream>
 
namespace Mantid::DataHandling {
 
DECLARE_FILELOADER_ALGORITHM(LoadNGEM)
 
// Constants and helper functions.
namespace {
constexpr int NUM_OF_SPECTRA = 16384;
 
uint64_t swapUint64(uint64_t word) {
  word = ((word << 8) & 0xFF00FF00FF00FF00ULL) | ((word >> 8) & 0x00FF00FF00FF00FFULL);
  word = ((word << 16) & 0xFFFF0000FFFF0000ULL) | ((word >> 16) & 0x0000FFFF0000FFFFULL);
  return (word << 32) | (word >> 32);
}
 
void correctForBigEndian(const EventUnion &bigEndian, EventUnion &smallEndian) {
  smallEndian.splitWord.words[0] = swapUint64(bigEndian.splitWord.words[1]);
  smallEndian.splitWord.words[1] = swapUint64(bigEndian.splitWord.words[0]);
}
 
void addFrameToOutputWorkspace(int &rawFrames, int &goodFrames, const int &eventCountInFrame, const int &minEventsReq,
                               const int &maxEventsReq, MantidVec &frameEventCounts,
                               std::vector<DataObjects::EventList> &events,
                               std::vector<DataObjects::EventList> &eventsInFrame) {
  ++rawFrames;
  if (eventCountInFrame >= minEventsReq && eventCountInFrame <= maxEventsReq) {
    // Add number of event counts to workspace.
    frameEventCounts.emplace_back(eventCountInFrame);
    ++goodFrames;
 
    PARALLEL_FOR_NO_WSP_CHECK()
    // Add events that match parameters to workspace
    for (auto i = 0; i < NUM_OF_SPECTRA; ++i) {
      if (eventsInFrame[i].getNumberEvents() > 0) {
        events[i] += eventsInFrame[i];
        eventsInFrame[i].clear();
      }
    }
  }
}
 
void createEventWorkspace(const double &maxToF, const double &binWidth,
                          const std::vector<DataObjects::EventList> &events,
                          DataObjects::EventWorkspace_sptr &dataWorkspace) {
  // Round up number of bins needed
  std::vector<double> xAxis(int(std::ceil(maxToF / binWidth)));
  std::generate(xAxis.begin(), xAxis.end(), [i = 0, &binWidth]() mutable { return binWidth * i++; });
 
  dataWorkspace = DataObjects::create<DataObjects::EventWorkspace>(
      NUM_OF_SPECTRA, HistogramData::Histogram(HistogramData::BinEdges(xAxis)));
  PARALLEL_FOR_NO_WSP_CHECK()
  for (int i = 0; i < NUM_OF_SPECTRA; ++i) {
    dataWorkspace->getSpectrum(i) = events[i];
    dataWorkspace->getSpectrum(i).setSpectrumNo(i + 1);
    dataWorkspace->getSpectrum(i).setDetectorID(i + 1);
  }
  dataWorkspace->setAllX(HistogramData::BinEdges{xAxis});
  dataWorkspace->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("TOF");
  dataWorkspace->setYUnit("Counts");
}
 
template <typename ValueType>
void addToSampleLog(const std::string &name, const ValueType &value, DataObjects::EventWorkspace &ws) {
  ws.mutableRun().addProperty(name, value, false);
}
 
} // namespace
 
int LoadNGEM::confidence(Kernel::FileDescriptor &descriptor) const {
  if (descriptor.extension() == ".edb") {
    return 95;
  } else {
    return 0;
  }
}
 
void LoadNGEM::init() {
  // Filename property.
  const std::vector<std::string> extentions{".edb"};
  declareProperty(std::make_unique<API::MultipleFileProperty>("Filename", extentions),
                  "The name of the nGEM file to load. Selecting multiple files will "
                  "combine them into one workspace.");
  // Output workspace
  declareProperty(
      std::make_unique<API::WorkspaceProperty<API::Workspace>>("OutputWorkspace", "", Kernel::Direction::Output),
      "The output workspace");
 
  auto mustBePositive = std::make_shared<Kernel::BoundedValidator<int>>();
  mustBePositive->setLower(0);
 
  auto mustBePositiveDbl = std::make_shared<Kernel::BoundedValidator<double>>();
  mustBePositiveDbl->setLower(0.0);
 
  // Bin Width
  declareProperty("BinWidth", 10.0, mustBePositiveDbl, "The width of the time bins in the output.");
 
  declareProperty("MinEventsPerFrame", 0, mustBePositive,
                  "The minimum number of events required in a frame before a "
                  "it is considered 'good'.");
  declareProperty("MaxEventsPerFrame", EMPTY_INT(), mustBePositive,
                  "The maximum number of events allowed in a frame to be "
                  "considered 'good'.");
  declareProperty(
      std::make_unique<Kernel::PropertyWithValue<bool>>("GenerateEventsPerFrame", false, Kernel::Direction::Input),
      "Generate a workspace to show the number of events captured by each "
      "frame. (optional, default False).");
}
 
void LoadNGEM::exec() {
  progress(0);
 
  std::vector<std::vector<std::string>> filePaths = getProperty("Filename");
 
  const int minEventsReq(getProperty("MinEventsPerFrame"));
  const int maxEventsReq(getProperty("MaxEventsPerFrame"));
 
  double maxToF{-std::numeric_limits<double>::max()}, minToF{std::numeric_limits<double>::max()};
  const double binWidth(getProperty("BinWidth"));
 
  int rawFrames = 0;
  int goodFrames = 0;
  std::vector<double> frameEventCounts;
  int eventCountInFrame = 0;
 
  std::vector<DataObjects::EventList> events, eventsInFrame;
  events.resize(NUM_OF_SPECTRA);
  eventsInFrame.resize(NUM_OF_SPECTRA);
  progress(0.04);
 
  size_t totalFilePaths(filePaths.size());
  int counter(1);
  for (const auto &filePath : filePaths) {
    loadSingleFile(filePath, eventCountInFrame, maxToF, minToF, rawFrames, goodFrames, minEventsReq, maxEventsReq,
                   frameEventCounts, events, eventsInFrame, totalFilePaths, counter);
  }
  // Add the final frame of events (as they are not followed by a T0 event)
  addFrameToOutputWorkspace(rawFrames, goodFrames, eventCountInFrame, minEventsReq, maxEventsReq, frameEventCounts,
                            events, eventsInFrame);
  progress(0.90);
 
  DataObjects::EventWorkspace_sptr dataWorkspace;
  createEventWorkspace(maxToF, binWidth, events, dataWorkspace);
 
  addToSampleLog("raw_frames", rawFrames, *dataWorkspace);
  addToSampleLog("good_frames", goodFrames, *dataWorkspace);
  addToSampleLog("max_ToF", maxToF, *dataWorkspace);
  addToSampleLog("min_ToF", minToF, *dataWorkspace);
 
  loadInstrument(dataWorkspace);
 
  setProperty("OutputWorkspace", dataWorkspace);
  if (this->getProperty("GenerateEventsPerFrame")) {
    createCountWorkspace(frameEventCounts);
  }
  progress(1.00);
}
 
void LoadNGEM::loadSingleFile(const std::vector<std::string> &filePath, int &eventCountInFrame, double &maxToF,
                              double &minToF, int &rawFrames, int &goodFrames, const int &minEventsReq,
                              const int &maxEventsReq, MantidVec &frameEventCounts,
                              std::vector<DataObjects::EventList> &events,
                              std::vector<DataObjects::EventList> &eventsInFrame, const size_t &totalFilePaths,
                              int &fileCount) {
  // Create file reader
  if (filePath.size() > 1) {
    throw std::runtime_error("Invalid filename parameter.");
  }
  std::ifstream file(filePath[0].c_str(), std::ifstream::binary);
  if (!file.is_open()) {
    throw std::runtime_error("File could not be found.");
  }
 
  const size_t totalNumEvents = verifyFileSize(file) / 16;
  constexpr size_t SKIP_WORD_SIZE = 4;
  size_t numProcessedEvents = 0;
  size_t numWordsSkipped = 0;
 
  while (true) {
    // Load an event into the variable.
    // Occasionally we may get a file where the first event has been chopped,
    // so we seek to the start of a valid event.
    // Chopping only seems to occur on a 4 byte word, hence seekg() of 4
    EventUnion event, eventBigEndian;
    do {
      file.read(reinterpret_cast<char *>(&eventBigEndian), sizeof(eventBigEndian));
      // Correct for the big endian format of nGEM datafile.
      correctForBigEndian(eventBigEndian, event);
    } while (!event.generic.check() && !file.seekg(SKIP_WORD_SIZE, std::ios_base::cur).eof() && ++numWordsSkipped);
    if (file.eof()) {
      break; // we have either not read an event, or only read part of one
    }
    if (event.coincidence.check()) { // Check for coincidence event.
      ++eventCountInFrame;
      uint64_t pixel = event.coincidence.getPixel();
      // Convert to microseconds (us)
      const double tof = event.coincidence.timeOfFlight / 1000.0;
 
      if (tof > maxToF) {
        maxToF = tof;
      } else if (tof < minToF) {
        minToF = tof;
      }
      eventsInFrame[pixel].addEventQuickly(Types::Event::TofEvent(tof));
 
    } else if (event.tZero.check()) { // Check for T0 event.
      addFrameToOutputWorkspace(rawFrames, goodFrames, eventCountInFrame, minEventsReq, maxEventsReq, frameEventCounts,
                                events, eventsInFrame);
 
      if (reportProgressAndCheckCancel(numProcessedEvents, eventCountInFrame, totalNumEvents, totalFilePaths,
                                       fileCount)) {
        return;
      }
    } else if (event.generic.check()) { // match all other events and notify.
      g_log.warning() << "Unexpected event type ID=" << event.generic.id << " loaded.\n";
    } else { // if we were to get to here, must be a corrupt event
      g_log.warning() << "Corrupt event detected.\n";
    }
  }
  if (numWordsSkipped > 0) {
    g_log.warning() << SKIP_WORD_SIZE * numWordsSkipped
                    << " bytes of file data were skipped when locating valid events.\n";
  }
  g_log.information() << "Finished loading a file.\n";
  ++fileCount;
}
 
size_t LoadNGEM::verifyFileSize(std::ifstream &file) {
  // Check that the file fits into 16 byte sections.
  file.seekg(0, file.end);
  size_t size = file.tellg();
  if (size % 16 != 0) {
    g_log.warning() << "Invalid file size. File is size is " << size
                    << " bytes which is not a multiple of 16. There may be some bytes "
                       "missing from the data. \n";
  }
  file.seekg(0);
  return size;
}
 
bool LoadNGEM::reportProgressAndCheckCancel(size_t &numProcessedEvents, int &eventCountInFrame,
                                            const size_t &totalNumEvents, const size_t &totalFilePaths,
                                            const int &fileCount) {
  numProcessedEvents += eventCountInFrame;
  std::string message(std::to_string(fileCount) + "/" + std::to_string(totalFilePaths));
  progress(double(numProcessedEvents) / double(totalNumEvents) / 1.11111, message);
  eventCountInFrame = 0;
  // Check for cancel flag.
  return this->getCancel();
}
 
void LoadNGEM::createCountWorkspace(const std::vector<double> &frameEventCounts) {
  std::vector<double> xAxisCounts(frameEventCounts.size() + 1);
  std::generate(xAxisCounts.begin(), xAxisCounts.end(), [n = 0.0]() mutable { return ++n; });
 
  DataObjects::Workspace2D_sptr countsWorkspace =
      DataObjects::create<DataObjects::Workspace2D>(1, HistogramData::Histogram(HistogramData::BinEdges(xAxisCounts)));
 
  countsWorkspace->mutableY(0) = frameEventCounts;
  std::string countsWorkspaceName(this->getProperty("OutputWorkspace"));
  countsWorkspaceName.append("_event_counts");
  countsWorkspace->setYUnit("Counts");
  std::shared_ptr<Kernel::Units::Label> XLabel =
      std::dynamic_pointer_cast<Kernel::Units::Label>(Kernel::UnitFactory::Instance().create("Label"));
  XLabel->setLabel("Frame");
  countsWorkspace->getAxis(0)->unit() = XLabel;
 
  this->declareProperty(std::make_unique<API::WorkspaceProperty<API::Workspace>>("CountsWorkspace", countsWorkspaceName,
                                                                                 Kernel::Direction::Output),
                        "Counts of events per frame.");
  progress(1.00);
  this->setProperty("CountsWorkspace", countsWorkspace);
}
 
void LoadNGEM::loadInstrument(DataObjects::EventWorkspace_sptr &dataWorkspace) {
  auto loadInstrument = this->createChildAlgorithm("LoadInstrument");
  loadInstrument->setPropertyValue("InstrumentName", "NGEM");
  loadInstrument->setProperty<API::MatrixWorkspace_sptr>("Workspace", dataWorkspace);
  loadInstrument->setProperty("RewriteSpectraMap", Kernel::OptionalBool(false));
  loadInstrument->execute();
}
 
std::map<std::string, std::string> LoadNGEM::validateInputs() {
  std::map<std::string, std::string> results;
 
  int MinEventsPerFrame = getProperty("MinEventsPerFrame");
  int MaxEventsPerFrame = getProperty("MaxEventsPerFrame");
 
  if (MaxEventsPerFrame < MinEventsPerFrame) {
    results["MaxEventsPerFrame"] = "MaxEventsPerFrame is less than MinEvents per frame";
  }
  return results;
}
 
} // namespace Mantid::DataHandling