Skip to content

Commit

Permalink
Add event filtering and loading logs
Browse files Browse the repository at this point in the history
  • Loading branch information
@peterfpeterson
peterfpeterson committed Jan 13, 2025
1 parent df055a9 commit d1c9962
Showing 1 changed file with 30 additions and 109 deletions.
139 changes: 30 additions & 109 deletions Framework/DataHandling/src/AlignAndFocusPowderSlim.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,8 @@
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/Sample.h"
#include "MantidDataHandling/LoadEventNexus.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/MaskWorkspace.h"
Expand All @@ -18,6 +20,7 @@
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidKernel/ArrayLengthValidator.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/Timer.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidNexus/H5Util.h"
Expand All @@ -38,6 +41,7 @@ using Mantid::DataObjects::Workspace2D;
using Mantid::Kernel::ArrayLengthValidator;
using Mantid::Kernel::ArrayProperty;
using Mantid::Kernel::Direction;
using Mantid::Kernel::TimeSeriesProperty;

namespace { // anonymous namespace
namespace PropertyNames {
Expand Down Expand Up @@ -103,30 +107,15 @@ class NexusLoader {
const std::pair<uint64_t, uint64_t> &eventRange) {
// g_log.information(NxsFieldNames::TIME_OF_FLIGHT);
auto tof_SDS = event_group.openDataSet(NxsFieldNames::TIME_OF_FLIGHT);
// TODO probably should resize data array
/*
// This is the data size
::NeXus::Info id_info = m_h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));
*/

// H5Util resizes the vector
if (m_is_time_filtered) {
throw std::runtime_error("filtering not implemented");
// TODO sort this out H5Util doesn't have slab read yet
/*
// These are the arguments to getSlab()
std::vector<int64_t> loadStart(1, eventRange.first);
int64_t slabsize = (eventRange.second == std::numeric_limits<size_t>::max())
? dim0 - eventRange.first
: eventRange.second - eventRange.first;
std::vector<int64_t> loadSize(1, slabsize);
data->resize(loadSize[0]);
Mantid::NeXus::NeXusIOHelper::readNexusSlab<float, Mantid::NeXus::NeXusIOHelper::PreventNarrowing>(
*data, m_h5file, NxsFieldNames::TIME_OF_FLIGHT, loadStart, loadSize);
*/
const size_t offset = eventRange.first;
const size_t slabsize = (eventRange.second == std::numeric_limits<size_t>::max())
? std::numeric_limits<size_t>::max()
: eventRange.second - eventRange.first;
NeXus::H5Util::readArray1DCoerce(tof_SDS, *data, slabsize, offset);
} else {
// TODO probably should resize data array
NeXus::H5Util::readArray1DCoerce(tof_SDS, *data);
}

Expand All @@ -143,30 +132,15 @@ class NexusLoader {
const std::pair<uint64_t, uint64_t> &eventRange) {
// g_log.information(NxsFieldNames::DETID);
auto detID_SDS = event_group.openDataSet(NxsFieldNames::DETID);
// TODO probably should resize data array
/*
// This is the data size
::NeXus::Info id_info = m_h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));
*/

// H5Util resizes the vector
if (m_is_time_filtered) {
throw std::runtime_error("filtering not implemented");
// TODO sort this out H5Util doesn't have slab read yet
/*
// These are the arguments to getSlab()
std::vector<int64_t> loadStart(1, eventRange.first);
int64_t slabsize = (eventRange.second == std::numeric_limits<size_t>::max())
? dim0 - eventRange.first
: eventRange.second - eventRange.first;
std::vector<int64_t> loadSize(1, slabsize);
data->resize(loadSize[0]);
Mantid::NeXus::NeXusIOHelper::readNexusSlab<uint32_t, Mantid::NeXus::NeXusIOHelper::PreventNarrowing>(
*data, m_h5file, NxsFieldNames::DETID, loadStart, loadSize);
*/
const size_t offset = eventRange.first;
const size_t slabsize = (eventRange.second == std::numeric_limits<size_t>::max())
? std::numeric_limits<size_t>::max()
: eventRange.second - eventRange.first;
NeXus::H5Util::readArray1DCoerce(detID_SDS, *data, slabsize, offset);
} else {
// TODO probably should resize data array
NeXus::H5Util::readArray1DCoerce(detID_SDS, *data);
}
}
Expand Down Expand Up @@ -452,7 +426,6 @@ void AlignAndFocusPowderSlim::exec() {

// Load the instrument
// prog->doReport("Loading instrument"); TODO add progress bar stuff
// LoadEventNexus::loadInstrument<MatrixWorkspace_sptr>(filename, wksp, "entry", this, &descriptor);
LoadEventNexus::loadInstrument<MatrixWorkspace_sptr>(filename, wksp, ENTRY_TOP_LEVEL, this, &descriptor);

const std::string cal_filename = getPropertyValue(PropertyNames::CAL_FILE);
Expand All @@ -463,25 +436,12 @@ void AlignAndFocusPowderSlim::exec() {
}

/*
// load run metadata
// prog->doReport("Loading metadata"); TODO add progress bar stuff
try {
LoadEventNexus::loadEntryMetadata(filename, WS, "entry", descriptor);
} catch (std::exception &e) {
g_log.warning() << "Error while loading meta data: " << e.what() << '\n';
}
// create IndexInfo
// prog->doReport("Creating IndexInfo"); TODO add progress bar stuff
const std::vector<int32_t> range;
LoadEventNexusIndexSetup indexSetup(WS, EMPTY_INT(), EMPTY_INT(), range);
auto indexInfo = indexSetup.makeIndexInfo();
const size_t numHist = indexInfo.size();
// make output workspace with correct number of histograms
MatrixWorkspace_sptr outWS = WorkspaceFactory::Instance().create(WS, numHist, 2, 1);
// set spectrum index information
outWS->setIndexInfo(indexInfo);
*/

// load the events
Expand Down Expand Up @@ -561,63 +521,24 @@ void AlignAndFocusPowderSlim::exec() {
m_calibration, m_masked, binWidth, linearBins);
constexpr size_t GRAINSIZE_BANK{2};
tbb::parallel_for(tbb::blocked_range<size_t>(0, bankEntryNames.size(), GRAINSIZE_BANK), task);
}

/*
NexusLoader loader(h5file, is_time_filtered, pulse_start_index, pulse_stop_index);
size_t specnum = 0;
for (const std::string &entry_name : bankEntryNames) {
std::cout << "ENTRY: " << entry_name << std::endl;
const auto startTimeBank = std::chrono::high_resolution_clock::now();
// TODO should re-use vectors to save malloc/free calls
std::unique_ptr<std::vector<uint32_t>> event_detid = std::make_unique<std::vector<uint32_t>>();
std::unique_ptr<std::vector<float>> event_time_of_flight = std::make_unique<std::vector<float>>();
// TODO std::unique_ptr<std::vector<float>> event_weight; some other time
g_log.information() << "Loading bank " << entry_name << '\n';
h5file.openGroup(entry_name, "NXevent_data");
// get filtering range
const auto eventRange = loader.getEventIndexRange(h5file);
loader.loadTOF(event_time_of_flight, eventRange);
loader.loadDetid(event_detid, eventRange);
if (event_time_of_flight->empty() || event_detid->empty()) {
g_log.warning() << "No data for bank " << entry_name << '\n';
h5file.closeGroup();
continue;
}
const auto startTimeSetup = std::chrono::high_resolution_clock::now();
const auto [minval, maxval] = parallel_minmax(event_detid.get());
BankCalibration calibration(static_cast<detid_t>(minval), static_cast<detid_t>(maxval), m_calibration);
auto &spectrum = wksp->getSpectrum(specnum);
Histogrammer histogrammer(&spectrum.readX(), binWidth, linearBins);
const auto numEvent = event_time_of_flight->size();
// std::atomic allows for multi-threaded accumulation and who cares about floats when you are just
// counting things
std::vector<std::atomic_uint32_t> y_temp(spectrum.dataY().size());
addTimer("setup" + entry_name, startTimeSetup, std::chrono::high_resolution_clock::now());
const auto startTimeProcess = std::chrono::high_resolution_clock::now();
constexpr size_t GRAINSIZE_EVENT{2000};
ProcessEventsTask task(&histogrammer, event_detid.get(), event_time_of_flight.get(), &calibration, &y_temp,
&m_masked);
tbb::parallel_for(tbb::blocked_range<size_t>(0, numEvent, GRAINSIZE_EVENT), task);
auto &y_values = spectrum.dataY();
std::copy(y_temp.cbegin(), y_temp.cend(), y_values.begin());
addTimer("proc" + entry_name, startTimeProcess, std::chrono::high_resolution_clock::now());
addTimer(entry_name, startTimeBank, std::chrono::high_resolution_clock::now());
// close the file so child algorithms can do their thing
h5file.close();

h5file.closeGroup();
specnum++;
}
*/
// load run metadata
// prog->doReport("Loading metadata"); TODO add progress bar stuff
try {
LoadEventNexus::loadEntryMetadata(filename, wksp, ENTRY_TOP_LEVEL, descriptor);
} catch (std::exception &e) {
g_log.warning() << "Error while loading meta data: " << e.what() << '\n';
}

// TODO load logs
// load logs
auto periodLog = std::make_unique<const TimeSeriesProperty<int>>("period_log"); // not used
int nPeriods{1};
LoadEventNexus::runLoadNexusLogs<MatrixWorkspace_sptr>(filename, wksp, *this, false, nPeriods, periodLog);

wksp->setYUnit("Counts");
wksp->getAxis(0)->setUnit("DSpacing");
setProperty(PropertyNames::OUTPUT_WKSP, std::move(wksp));
Expand Down

0 comments on commit d1c9962

Please sign in to comment.