Skip to content

Commit

Permalink
Refactor code in preparation of having a loading task
Browse files Browse the repository at this point in the history
  • Loading branch information
@peterfpeterson
peterfpeterson committed Jan 9, 2025
1 parent 8e6f386 commit 275a87b
Show file tree
Hide file tree
Showing 2 changed files with 137 additions and 115 deletions.
5 changes: 0 additions & 5 deletions Framework/DataHandling/inc/MantidDataHandling/AlignAndFocusPowderSlim.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,11 +44,6 @@ class MANTID_DATAHANDLING_DLL AlignAndFocusPowderSlim : public API::Algorithm {

void initCalibrationConstants(API::MatrixWorkspace_sptr &wksp);

void loadTOF(std::unique_ptr<std::vector<float>> &data, ::NeXus::File &h5file);
void loadDetid(std::unique_ptr<std::vector<uint32_t>> &data, ::NeXus::File &h5file);
void loadPulseTimes(std::unique_ptr<std::vector<double>> &data, ::NeXus::File &h5file);
void loadEventIndex(std::unique_ptr<std::vector<uint64_t>> &data, ::NeXus::File &h5file);

void loadCalFile(const Mantid::API::Workspace_sptr &inputWS, const std::string &filename);

std::map<detid_t, double> m_calibration; // detid: 1/difc
Expand Down
247 changes: 137 additions & 110 deletions Framework/DataHandling/src/AlignAndFocusPowderSlim.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,134 @@ const std::vector<std::string> AlignAndFocusPowderSlim::seeAlso() const { return

//----------------------------------------------------------------------------------------------
namespace { // anonymous
class NexusLoader {
public:
NexusLoader(const bool is_time_filtered, const size_t pulse_start_index, const size_t pulse_stop_index)
: m_is_time_filtered(is_time_filtered), m_pulse_start_index(pulse_start_index),
m_pulse_stop_index(pulse_stop_index) {}

static void loadPulseTimes(std::unique_ptr<std::vector<double>> &data, ::NeXus::File &h5file) {
// /entry/DASlogs/frequency/time
h5file.openGroup("DASlogs", "NXcollection");
h5file.openGroup("frequency", "NXlog");
h5file.openData("time");

// This is the data size
::NeXus::Info id_info = h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));
data->resize(dim0);

Mantid::NeXus::NeXusIOHelper::readNexusVector<double>(*data, h5file, "time");

// close the sds
h5file.closeData();
h5file.closeGroup();
h5file.closeGroup();
}

void loadTOF(std::unique_ptr<std::vector<float>> &data, ::NeXus::File &h5file,
const std::pair<uint64_t, uint64_t> &eventRange) {
// g_log.information(NxsFieldNames::TIME_OF_FLIGHT);
h5file.openData(NxsFieldNames::TIME_OF_FLIGHT);

// This is the data size
::NeXus::Info id_info = h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));

if (m_is_time_filtered) {
// 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, h5file, NxsFieldNames::TIME_OF_FLIGHT, loadStart, loadSize);
} else {
data->resize(dim0);
Mantid::NeXus::NeXusIOHelper::readNexusVector<float>(*data, h5file, NxsFieldNames::TIME_OF_FLIGHT);
}

// get the units
std::string tof_unit;
h5file.getAttr("units", tof_unit);

// close the sds
h5file.closeData();

// Convert Tof to microseconds
if (tof_unit != MICROSEC)
Kernel::Units::timeConversionVector(*data, tof_unit, MICROSEC);
}

void loadDetid(std::unique_ptr<std::vector<uint32_t>> &data, ::NeXus::File &h5file,
const std::pair<uint64_t, uint64_t> &eventRange) {
// g_log.information(NxsFieldNames::DETID);
h5file.openData(NxsFieldNames::DETID);

// This is the data size
::NeXus::Info id_info = h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));

if (m_is_time_filtered) {
// 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, h5file, NxsFieldNames::DETID, loadStart, loadSize);
} else {
data->resize(dim0);
Mantid::NeXus::NeXusIOHelper::readNexusVector<uint32_t>(*data, h5file, NxsFieldNames::DETID);
}

// close the sds
h5file.closeData();
}

private:
void loadEventIndex(std::unique_ptr<std::vector<uint64_t>> &data, ::NeXus::File &h5file) {
// g_log.information(NxsFieldNames::INDEX_ID);
h5file.openData(NxsFieldNames::INDEX_ID);

// This is the data size
::NeXus::Info id_info = h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));
data->resize(dim0);

Mantid::NeXus::NeXusIOHelper::readNexusVector<uint64_t>(*data, h5file, NxsFieldNames::INDEX_ID);

// close the sds
h5file.closeData();
}

public:
std::pair<uint64_t, uint64_t> getEventIndexRange(::NeXus::File &h5file) {
uint64_t start_event = 0;
uint64_t stop_event = std::numeric_limits<uint64_t>::max();
if (m_is_time_filtered) {
std::unique_ptr<std::vector<uint64_t>> event_index = std::make_unique<std::vector<uint64_t>>();
this->loadEventIndex(event_index, h5file);

start_event = event_index->at(m_pulse_stop_index);
if (m_pulse_stop_index != std::numeric_limits<size_t>::max())
stop_event = event_index->at(m_pulse_stop_index);
}
return {start_event, stop_event};
}

private:
const bool m_is_time_filtered;
const size_t m_pulse_start_index;
const size_t m_pulse_stop_index;
};

class Histogrammer {
public:
Histogrammer(const std::vector<double> *binedges, const double width, const bool linear_bins) : m_binedges(binedges) {
Expand Down Expand Up @@ -296,7 +424,7 @@ void AlignAndFocusPowderSlim::exec() {
is_time_filtered = true;
g_log.information() << "Filtering pulses from " << filter_time_start_sec << " to " << filter_time_stop_sec << "s\n";
std::unique_ptr<std::vector<double>> pulse_times = std::make_unique<std::vector<double>>();
loadPulseTimes(pulse_times, h5file);
NexusLoader::loadPulseTimes(pulse_times, h5file);
g_log.information() << "Pulse times from " << pulse_times->front() << " to " << pulse_times->back()
<< " with length " << pulse_times->size() << '\n';
if (!std::is_sorted(pulse_times->cbegin(), pulse_times->cend())) {
Expand Down Expand Up @@ -356,6 +484,8 @@ void AlignAndFocusPowderSlim::exec() {
}
}

NexusLoader loader(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;
Expand All @@ -365,28 +495,20 @@ void AlignAndFocusPowderSlim::exec() {
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
std::unique_ptr<std::vector<uint64_t>> event_index = std::make_unique<std::vector<uint64_t>>();
g_log.information() << "Loading bank " << entry_name << '\n';
h5file.openGroup(entry_name, "NXevent_data");

if (is_time_filtered) {
const auto startTime = std::chrono::high_resolution_clock::now();
loadEventIndex(event_index, h5file);
addTimer("loadEventIndex" + entry_name, startTime, std::chrono::high_resolution_clock::now());
start_event = event_index->at(pulse_start_index);
if (pulse_stop_index != std::numeric_limits<size_t>::max())
stop_event = event_index->at(pulse_stop_index);
g_log.debug() << "Loading events from " << start_event << " to " << stop_event << '\n';
}
// get filtering range
const auto eventRange = loader.getEventIndexRange(h5file);

{
const auto startTime = std::chrono::high_resolution_clock::now();
loadTOF(event_time_of_flight, h5file);
loader.loadTOF(event_time_of_flight, h5file, eventRange);
addTimer("readTOF" + entry_name, startTime, std::chrono::high_resolution_clock::now());
}
{
const auto startTime = std::chrono::high_resolution_clock::now();
loadDetid(event_detid, h5file);
loader.loadDetid(event_detid, h5file, eventRange);
addTimer("readDetID" + entry_name, startTime, std::chrono::high_resolution_clock::now());
}

Expand All @@ -409,9 +531,10 @@ void AlignAndFocusPowderSlim::exec() {
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), task);
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());
Expand Down Expand Up @@ -442,102 +565,6 @@ void AlignAndFocusPowderSlim::initCalibrationConstants(API::MatrixWorkspace_sptr
}
}

void AlignAndFocusPowderSlim::loadTOF(std::unique_ptr<std::vector<float>> &data, ::NeXus::File &h5file) {
g_log.information(NxsFieldNames::TIME_OF_FLIGHT);
h5file.openData(NxsFieldNames::TIME_OF_FLIGHT);

// This is the data size
::NeXus::Info id_info = h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));

if (is_time_filtered) {
// These are the arguments to getSlab()
loadStart[0] = start_event;
if (stop_event == std::numeric_limits<size_t>::max())
loadSize[0] = dim0 - start_event;
else
loadSize[0] = stop_event - start_event;
data->resize(loadSize[0]);
Mantid::NeXus::NeXusIOHelper::readNexusSlab<float, Mantid::NeXus::NeXusIOHelper::PreventNarrowing>(
*data, h5file, NxsFieldNames::TIME_OF_FLIGHT, loadStart, loadSize);
} else {
data->resize(dim0);
Mantid::NeXus::NeXusIOHelper::readNexusVector<float>(*data, h5file, NxsFieldNames::TIME_OF_FLIGHT);
}

// get the units
std::string tof_unit;
h5file.getAttr("units", tof_unit);

// close the sds
h5file.closeData();

// Convert Tof to microseconds
if (tof_unit != MICROSEC)
Kernel::Units::timeConversionVector(*data, tof_unit, MICROSEC);
}

void AlignAndFocusPowderSlim::loadDetid(std::unique_ptr<std::vector<uint32_t>> &data, ::NeXus::File &h5file) {
g_log.information(NxsFieldNames::DETID);
h5file.openData(NxsFieldNames::DETID);

// This is the data size
::NeXus::Info id_info = h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));

if (is_time_filtered) {
// These are the arguments to getSlab()
loadStart[0] = start_event;
if (stop_event == std::numeric_limits<size_t>::max())
loadSize[0] = dim0 - start_event;
else
loadSize[0] = stop_event - start_event;
data->resize(loadSize[0]);
Mantid::NeXus::NeXusIOHelper::readNexusSlab<uint32_t, Mantid::NeXus::NeXusIOHelper::PreventNarrowing>(
*data, h5file, NxsFieldNames::DETID, loadStart, loadSize);
} else {
data->resize(dim0);
Mantid::NeXus::NeXusIOHelper::readNexusVector<uint32_t>(*data, h5file, NxsFieldNames::DETID);
}

// close the sds
h5file.closeData();
}

void AlignAndFocusPowderSlim::loadPulseTimes(std::unique_ptr<std::vector<double>> &data, ::NeXus::File &h5file) {
// /entry/DASlogs/frequency/time
h5file.openGroup("DASlogs", "NXcollection");
h5file.openGroup("frequency", "NXlog");
h5file.openData("time");

// This is the data size
::NeXus::Info id_info = h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));
data->resize(dim0);

Mantid::NeXus::NeXusIOHelper::readNexusVector<double>(*data, h5file, "time");

// close the sds
h5file.closeData();
h5file.closeGroup();
h5file.closeGroup();
}

void AlignAndFocusPowderSlim::loadEventIndex(std::unique_ptr<std::vector<uint64_t>> &data, ::NeXus::File &h5file) {
g_log.information(NxsFieldNames::INDEX_ID);
h5file.openData(NxsFieldNames::INDEX_ID);

// This is the data size
::NeXus::Info id_info = h5file.getInfo();
const auto dim0 = static_cast<size_t>(LoadBankFromDiskTask::recalculateDataSize(id_info.dims[0]));
data->resize(dim0);

Mantid::NeXus::NeXusIOHelper::readNexusVector<uint64_t>(*data, h5file, NxsFieldNames::INDEX_ID);

// close the sds
h5file.closeData();
}

void AlignAndFocusPowderSlim::loadCalFile(const Mantid::API::Workspace_sptr &inputWS, const std::string &filename) {
auto alg = createChildAlgorithm("LoadDiffCal");
alg->setProperty("InputWorkspace", inputWS);
Expand Down

0 comments on commit 275a87b

Please sign in to comment.