From d242a5491212dd3dc25668c8e29ef9e49b842e32 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 12 Mar 2024 12:59:54 +0000 Subject: [PATCH 01/35] Clarify event partitioning modes and only check / create window definition if selected. --- np.cpp | 79 +++++++++++++++++++++++++----------------------- src/processors.h | 4 +-- 2 files changed, 44 insertions(+), 39 deletions(-) diff --git a/np.cpp b/np.cpp index f609454..bfb9c85 100644 --- a/np.cpp +++ b/np.cpp @@ -60,7 +60,7 @@ int main(int argc, char **argv) [&]() { if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::Summed; + processingMode_ = Processors::ProcessingMode::PartitionEventsSummed; else { fmt::print("Error: Multiple processing modes given.\n"); @@ -74,7 +74,7 @@ int main(int argc, char **argv) [&]() { if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::Individual; + processingMode_ = Processors::ProcessingMode::PartitionEventsIndividual; else { fmt::print("Error: Multiple processing modes given.\n"); @@ -97,53 +97,58 @@ int main(int argc, char **argv) CLI11_PARSE(app, argc, argv); - // Sanity check - if ((windowWidth_ + windowOffset_) > windowDelta_) - { - fmt::print("Error: Window width (including any optional offset) is greater than window delta.\n"); - return 1; - } - if (windowSlices_ < 1) - { - fmt::print("Error: Invalid number of window slices provided ({}).\n", windowSlices_); - return 1; - } - // Perform pre-processing if requested if (spectrumId_) { Processors::getEvents(inputFiles_, *spectrumId_); } - // Construct the master window definition - if (relativeStartTime_) - { - // Need to query first NeXuS file to get its start time - if (inputFiles_.empty()) - { - fmt::print("Error: Need at least one input NeXuS file."); - return 1; - } - NeXuSFile firstFile(inputFiles_.front()); - firstFile.loadTimes(); - fmt::print("Window start time converted from relative to absolute time: {} => {} (= {} + {})\n", windowStartTime_, - windowStartTime_ + firstFile.startSinceEpoch(), firstFile.startSinceEpoch(), windowStartTime_); - windowStartTime_ += firstFile.startSinceEpoch(); - } - Window window(windowName_, windowStartTime_ + windowOffset_, windowWidth_); - fmt::print("Window start time (including any offset) is {}.\n", window.startTime()); - // Perform processing switch (processingMode_) { case (Processors::ProcessingMode::None): fmt::print("No processing mode specified. We are done.\n"); break; - case (Processors::ProcessingMode::Individual): - Processors::processIndividual(inputFiles_, outputDirectory_, window, windowSlices_, windowDelta_); - break; - case (Processors::ProcessingMode::Summed): - Processors::processSummed(inputFiles_, outputDirectory_, window, windowSlices_, windowDelta_); + case (Processors::ProcessingMode::PartitionEventsIndividual): + case (Processors::ProcessingMode::PartitionEventsSummed): + // Sanity check + if ((windowWidth_ + windowOffset_) > windowDelta_) + { + fmt::print("Error: Window width (including any optional offset) is greater than window delta.\n"); + return 1; + } + if (windowSlices_ < 1) + { + fmt::print("Error: Invalid number of window slices provided ({}).\n", windowSlices_); + return 1; + } + + // Construct the master window definition + if (relativeStartTime_) + { + // Need to query first NeXuS file to get its start time + if (inputFiles_.empty()) + { + fmt::print("Error: Need at least one input NeXuS file."); + return 1; + } + NeXuSFile firstFile(inputFiles_.front()); + firstFile.loadTimes(); + fmt::print("Window start time converted from relative to absolute time: {} => {} (= {} + {})\n", + windowStartTime_, windowStartTime_ + firstFile.startSinceEpoch(), firstFile.startSinceEpoch(), + windowStartTime_); + windowStartTime_ += firstFile.startSinceEpoch(); + } + fmt::print("Window start time (including any offset) is {}.\n", windowStartTime_ + windowOffset_); + + if (processingMode_ == Processors::ProcessingMode::PartitionEventsIndividual) + Processors::processIndividual(inputFiles_, outputDirectory_, + {windowName_, windowStartTime_ + windowOffset_, windowWidth_}, windowSlices_, + windowDelta_); + else + Processors::processSummed(inputFiles_, outputDirectory_, + {windowName_, windowStartTime_ + windowOffset_, windowWidth_}, windowSlices_, + windowDelta_); break; default: throw(std::runtime_error("Unhandled processing mode.\n")); diff --git a/src/processors.h b/src/processors.h index 6ed0831..a8cc952 100644 --- a/src/processors.h +++ b/src/processors.h @@ -13,8 +13,8 @@ namespace Processors enum class ProcessingMode { None, - Individual, - Summed + PartitionEventsIndividual, + PartitionEventsSummed }; // Processing Direction From bd188b4a06f5bb04023ed2f8982fe7b5710c4cb0 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 12 Mar 2024 13:00:53 +0000 Subject: [PATCH 02/35] Window arguments are now optional. --- np.cpp | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/np.cpp b/np.cpp index bfb9c85..2faff21 100644 --- a/np.cpp +++ b/np.cpp @@ -33,19 +33,16 @@ int main(int argc, char **argv) CLI::App app("NeXuS Processor (np), Copyright (C) 2024 Jared Swift and Tristan Youngs."); // -- Window Definition app.add_option("-n,--name", windowName_, "Name of the window, used as a prefix to all output files") - ->group("Window Definition") - ->required(); + ->group("Window Definition"); app.add_option("-s,--start", windowStartTime_, "Start time of the window (relative to first input file start time unless --absolute-start is given)") ->group("Window Definition"); - app.add_option("-w,--width", windowWidth_, "Window width in seconds)")->group("Window Definition")->required(); + app.add_option("-w,--width", windowWidth_, "Window width in seconds)")->group("Window Definition"); app.add_flag( "--relative-start", relativeStartTime_, "Flag that the given window start time is relative to the first run start time, not absolute (seconds since epoch)") ->group("Window Definition"); - app.add_option("-d,--delta", windowDelta_, "Time between window occurrences, in seconds") - ->group("Window Definition") - ->required(); + app.add_option("-d,--delta", windowDelta_, "Time between window occurrences, in seconds")->group("Window Definition"); app.add_option("--offset", windowOffset_, "Time after start time, in seconds, that the window begins.") ->group("Window Definition"); // -- Input Files From 2cb203e83687b146b51f0af5afbc33ce43fade0d Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 12 Mar 2024 13:15:23 +0000 Subject: [PATCH 03/35] Keep all window options together. --- np.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/np.cpp b/np.cpp index 2faff21..938333f 100644 --- a/np.cpp +++ b/np.cpp @@ -45,6 +45,8 @@ int main(int argc, char **argv) app.add_option("-d,--delta", windowDelta_, "Time between window occurrences, in seconds")->group("Window Definition"); app.add_option("--offset", windowOffset_, "Time after start time, in seconds, that the window begins.") ->group("Window Definition"); + app.add_option("-l,--slices", windowSlices_, "Number of slices to split window definition in to (default = 1, no slicing)") + ->group("Window Definition"); // -- Input Files app.add_option("-f,--files", inputFiles_, "List of NeXuS files to process")->group("Input Files")->required(); // -- Output Files @@ -80,8 +82,6 @@ int main(int argc, char **argv) }, "Output NeXuS files for each window / slice") ->group("Processing"); - app.add_option("-l,--slices", windowSlices_, "Number of slices to split window definition in to (default = 1, no slicing)") - ->group("Processing"); // -- Post Processing app.add_flag_callback( "--scale-monitors", [&]() { Processors::postProcessingMode_ = Processors::PostProcessingMode::ScaleMonitors; }, From 70c5651ce38e75035f34d3c35108ef12db96c49e Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 12 Mar 2024 13:15:32 +0000 Subject: [PATCH 04/35] Create DumpEvents option. --- np.cpp | 28 +++++++++++++++++++--------- src/processors.h | 1 + 2 files changed, 20 insertions(+), 9 deletions(-) diff --git a/np.cpp b/np.cpp index 938333f..73bed81 100644 --- a/np.cpp +++ b/np.cpp @@ -27,7 +27,7 @@ int main(int argc, char **argv) // Number of slices to partition window in to int windowSlices_{1}; // Target spectrum for event get (optional) - std::optional spectrumId_; + int spectrumId_; // Define and parse CLI arguments CLI::App app("NeXuS Processor (np), Copyright (C) 2024 Jared Swift and Tristan Youngs."); @@ -51,9 +51,22 @@ int main(int argc, char **argv) app.add_option("-f,--files", inputFiles_, "List of NeXuS files to process")->group("Input Files")->required(); // -- Output Files app.add_option("--output-dir", outputDirectory_, "Output directory for generated NeXuS files.")->group("Output Files"); - // -- Pre Processing - app.add_option("-g,--get", spectrumId_, "Get all events from specified spectrum index")->group("Pre-Processing"); // -- Processing Modes + app.add_option_function( + "--dump-events", + [&](int specId) + { + if (processingMode_ == Processors::ProcessingMode::None) + processingMode_ = Processors::ProcessingMode::DumpEvents; + else + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + spectrumId_ = specId; + }, + "Dump all events from specified spectrum index") + ->group("Processing"); app.add_flag_callback( "--summed", [&]() @@ -94,18 +107,15 @@ int main(int argc, char **argv) CLI11_PARSE(app, argc, argv); - // Perform pre-processing if requested - if (spectrumId_) - { - Processors::getEvents(inputFiles_, *spectrumId_); - } - // Perform processing switch (processingMode_) { case (Processors::ProcessingMode::None): fmt::print("No processing mode specified. We are done.\n"); break; + case (Processors::ProcessingMode::DumpEvents): + Processors::getEvents(inputFiles_, spectrumId_); + break; case (Processors::ProcessingMode::PartitionEventsIndividual): case (Processors::ProcessingMode::PartitionEventsSummed): // Sanity check diff --git a/src/processors.h b/src/processors.h index a8cc952..9191d9d 100644 --- a/src/processors.h +++ b/src/processors.h @@ -13,6 +13,7 @@ namespace Processors enum class ProcessingMode { None, + DumpEvents, PartitionEventsIndividual, PartitionEventsSummed }; From eb432098ef08968772bd1d5db0490364259fed7b Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 12 Mar 2024 13:17:47 +0000 Subject: [PATCH 05/35] Rename function, adjust output text. --- np.cpp | 2 +- src/getEvents.cpp | 7 +++---- src/processors.h | 4 ++-- 3 files changed, 6 insertions(+), 7 deletions(-) diff --git a/np.cpp b/np.cpp index 73bed81..c5e2276 100644 --- a/np.cpp +++ b/np.cpp @@ -114,7 +114,7 @@ int main(int argc, char **argv) fmt::print("No processing mode specified. We are done.\n"); break; case (Processors::ProcessingMode::DumpEvents): - Processors::getEvents(inputFiles_, spectrumId_); + Processors::dumpEvents(inputFiles_, spectrumId_); break; case (Processors::ProcessingMode::PartitionEventsIndividual): case (Processors::ProcessingMode::PartitionEventsSummed): diff --git a/src/getEvents.cpp b/src/getEvents.cpp index 0d99163..b1d993d 100644 --- a/src/getEvents.cpp +++ b/src/getEvents.cpp @@ -6,15 +6,14 @@ namespace Processors { -// Get events from specified spectrum, returning seconds since epoch for each -std::map> getEvents(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) +// Dump events from specified spectrum, returning seconds since epoch for each +std::map> dumpEvents(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) { /* * Dump all events for the specified detector spectrum */ - printf("Get events...\n"); - fmt::print("Target detector spectrum is {}\n", spectrumId); + fmt::print("Dumping all events from detector spectrum {}...\n", spectrumId); std::map> eventMap; std::optional lastSecondsSinceEpoch; diff --git a/src/processors.h b/src/processors.h index 9191d9d..d33bb31 100644 --- a/src/processors.h +++ b/src/processors.h @@ -53,8 +53,8 @@ void saveSlices(std::vector> &slices); */ // Get Events -std::map> getEvents(const std::vector &inputNeXusFiles, int detectorId, - bool firstOnly = false); +std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorId, + bool firstOnly = false); // Perform individual processing void processIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, const Window &windowDefinition, int nSlices, double windowDelta); From 8443bdf8360032d41a0e06a95047434476a331a4 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 12 Mar 2024 13:21:03 +0000 Subject: [PATCH 06/35] Rename event partitioning functions. --- np.cpp | 12 ++++++------ src/CMakeLists.txt | 6 +++--- src/{getEvents.cpp => dumpEvents.cpp} | 0 ...sIndividual.cpp => partitionEventsIndividual.cpp} | 8 ++++---- src/{processSummed.cpp => partitionEventsSummed.cpp} | 9 +++++---- src/processors.h | 12 ++++++------ 6 files changed, 24 insertions(+), 23 deletions(-) rename src/{getEvents.cpp => dumpEvents.cpp} (100%) rename src/{processIndividual.cpp => partitionEventsIndividual.cpp} (91%) rename src/{processSummed.cpp => partitionEventsSummed.cpp} (90%) diff --git a/np.cpp b/np.cpp index c5e2276..a253df6 100644 --- a/np.cpp +++ b/np.cpp @@ -149,13 +149,13 @@ int main(int argc, char **argv) fmt::print("Window start time (including any offset) is {}.\n", windowStartTime_ + windowOffset_); if (processingMode_ == Processors::ProcessingMode::PartitionEventsIndividual) - Processors::processIndividual(inputFiles_, outputDirectory_, - {windowName_, windowStartTime_ + windowOffset_, windowWidth_}, windowSlices_, - windowDelta_); + Processors::partitionEventsIndividual(inputFiles_, outputDirectory_, + {windowName_, windowStartTime_ + windowOffset_, windowWidth_}, + windowSlices_, windowDelta_); else - Processors::processSummed(inputFiles_, outputDirectory_, - {windowName_, windowStartTime_ + windowOffset_, windowWidth_}, windowSlices_, - windowDelta_); + Processors::partitionEventsSummed(inputFiles_, outputDirectory_, + {windowName_, windowStartTime_ + windowOffset_, windowWidth_}, windowSlices_, + windowDelta_); break; default: throw(std::runtime_error("Unhandled processing mode.\n")); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2b200fa..2baa8e8 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,9 +1,9 @@ add_library(nexusProcess - getEvents.cpp + dumpEvents.cpp nexusFile.cpp + partitionEventsIndividual.cpp + partitionEventsSummed.cpp processCommon.cpp - processIndividual.cpp - processSummed.cpp window.cpp nexusFile.h processors.h diff --git a/src/getEvents.cpp b/src/dumpEvents.cpp similarity index 100% rename from src/getEvents.cpp rename to src/dumpEvents.cpp diff --git a/src/processIndividual.cpp b/src/partitionEventsIndividual.cpp similarity index 91% rename from src/processIndividual.cpp rename to src/partitionEventsIndividual.cpp index bb26620..ae52ecb 100644 --- a/src/processIndividual.cpp +++ b/src/partitionEventsIndividual.cpp @@ -6,16 +6,16 @@ namespace Processors { -// Perform summed processing -void processIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, - const Window &windowDefinition, int nSlices, double windowDelta) +// Partition events into individual windows / slices +void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, + const Window &windowDefinition, int nSlices, double windowDelta) { /* * From our main windowDefinition we will continually propagate it forwards in time (by the window delta) splitting it into * nSlices and until we go over the end time of the current file. */ - fmt::print("Processing in INDIVIDUAL mode...\n"); + fmt::print("Partitioning events into individual windows/slices...\n"); // Copy the master window definition since we will be modifying it auto window = windowDefinition; diff --git a/src/processSummed.cpp b/src/partitionEventsSummed.cpp similarity index 90% rename from src/processSummed.cpp rename to src/partitionEventsSummed.cpp index 6255fdd..2d90eee 100644 --- a/src/processSummed.cpp +++ b/src/partitionEventsSummed.cpp @@ -1,20 +1,21 @@ #include "nexusFile.h" #include "processors.h" #include "window.h" +#include #include namespace Processors { -// Perform summed processing -void processSummed(const std::vector &inputNeXusFiles, std::string_view outputFilePath, - const Window &windowDefinition, int nSlices, double windowDelta) +// Partition events into summed windows / slices +void partitionEventsSummed(const std::vector &inputNeXusFiles, std::string_view outputFilePath, + const Window &windowDefinition, int nSlices, double windowDelta) { /* * From our main windowDefinition we will continually propagate it forwards in time (by the window delta) splitting it into * nSlices and until we go over the end time of the current file. */ - printf("Processing in SUMMED mode...\n"); + fmt::print("Partitioning events in summed windows/slices...\n"); // Generate a new set of window "slices" and associated output NeXuS files to sum data into auto slices = prepareSlices(windowDefinition, nSlices, inputNeXusFiles[0], outputFilePath); diff --git a/src/processors.h b/src/processors.h index d33bb31..2f0acd9 100644 --- a/src/processors.h +++ b/src/processors.h @@ -55,10 +55,10 @@ void saveSlices(std::vector> &slices); // Get Events std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorId, bool firstOnly = false); -// Perform individual processing -void processIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, - const Window &windowDefinition, int nSlices, double windowDelta); -// Perform summed processing -void processSummed(const std::vector &inputNeXusFiles, std::string_view outputFilePath, - const Window &windowDefinition, int nSlices, double windowDelta); +// Partition events into individual windows / slices +void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, + const Window &windowDefinition, int nSlices, double windowDelta); +// Partition events into summed windows / slices +void partitionEventsSummed(const std::vector &inputNeXusFiles, std::string_view outputFilePath, + const Window &windowDefinition, int nSlices, double windowDelta); }; // namespace Processors From e55e095d83c780f1a7509d58a3bac3499d854f59 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 15 Oct 2024 10:07:50 +0100 Subject: [PATCH 07/35] Adding --dump-histogram option. --- np.cpp | 18 ++++++++++++++++ src/CMakeLists.txt | 1 + src/dumpEvents.cpp | 2 +- src/modex.cpp | 2 +- src/nexusFile.cpp | 35 ++++++++++++++++++++++++++----- src/nexusFile.h | 3 +++ src/partitionEventsIndividual.cpp | 2 +- src/partitionEventsSummed.cpp | 2 +- src/processors.h | 3 +++ 9 files changed, 59 insertions(+), 9 deletions(-) diff --git a/np.cpp b/np.cpp index a253df6..7b98010 100644 --- a/np.cpp +++ b/np.cpp @@ -67,6 +67,21 @@ int main(int argc, char **argv) }, "Dump all events from specified spectrum index") ->group("Processing"); + app.add_option_function( + "--dump-histogram", + [&](int specId) + { + if (processingMode_ == Processors::ProcessingMode::None) + processingMode_ = Processors::ProcessingMode::DumpHistogram; + else + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + spectrumId_ = specId; + }, + "Dump specified histogram index") + ->group("Processing"); app.add_flag_callback( "--summed", [&]() @@ -116,6 +131,9 @@ int main(int argc, char **argv) case (Processors::ProcessingMode::DumpEvents): Processors::dumpEvents(inputFiles_, spectrumId_); break; + case (Processors::ProcessingMode::DumpHistogram): + Processors::dumpHistogram(inputFiles_, spectrumId_); + break; case (Processors::ProcessingMode::PartitionEventsIndividual): case (Processors::ProcessingMode::PartitionEventsSummed): // Sanity check diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2baa8e8..546ba1c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,6 @@ add_library(nexusProcess dumpEvents.cpp + dumpHistogram.cpp nexusFile.cpp partitionEventsIndividual.cpp partitionEventsSummed.cpp diff --git a/src/dumpEvents.cpp b/src/dumpEvents.cpp index b1d993d..a0022d7 100644 --- a/src/dumpEvents.cpp +++ b/src/dumpEvents.cpp @@ -18,7 +18,7 @@ std::map> dumpEvents(const std::vector &in std::map> eventMap; std::optional lastSecondsSinceEpoch; - // Loop over input Nexus files + // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) { // Open the Nexus file ready for use diff --git a/src/modex.cpp b/src/modex.cpp index a767bee..9a07fc6 100644 --- a/src/modex.cpp +++ b/src/modex.cpp @@ -63,7 +63,7 @@ bool ModEx::process() return false; } - // Loop over input Nexus files + // Loop over input NeXuS files for (auto &nxsFileName : cfg.runs) { // Open the Nexus file ready for use diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index fe170d8..a7a85dc 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -67,7 +67,7 @@ void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) { filename_ = outputFile; - // Open input Nexus file in read only mode. + // Open input NeXuS file in read only mode. H5::H5File input = H5::H5File(referenceFile, H5F_ACC_RDONLY); // Create new Nexus file for output. @@ -103,11 +103,16 @@ void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) tofBins_.resize(tofBinsDimension); H5Dread(tofBinsID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBins_.data()); - // Set up detector histograms + // Set up detector histograms and straight counts vectors for (auto spec : spectra_) { + // For event re-binning detectorHistograms_[spec] = gsl_histogram_alloc(tofBins_.size() - 1); gsl_histogram_set_ranges(detectorHistograms_[spec], tofBins_.data(), tofBins_.size()); + + // For histogram manipulation + detectorCounts_[spec].resize(tofBins_.size() - 1); + std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); } // Read in monitor data - start from index 1 and end when we fail to find the named dataset with this suffix @@ -140,7 +145,7 @@ void NeXuSFile::loadFrameCounts() { printf("Load frame counts...\n"); - // Open our Nexus file in read only mode. + // Open our NeXuS file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); // Read in good frames @@ -157,7 +162,7 @@ void NeXuSFile::loadEventData() { printf("Load event data...\n"); - // Open our Nexus file in read only mode. + // Open our NeXuS file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); // Read in event indices. @@ -192,7 +197,7 @@ void NeXuSFile::loadTimes() { printf("Load times....\n"); - // Open our Nexus file in read only mode. + // Open our NeXuS file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); // Read in start time in Unix time. @@ -232,6 +237,26 @@ void NeXuSFile::loadTimes() input.close(); } +// Load detector counts from the file +void NeXuSFile::loadDetectorCounts() +{ + printf("Load detector events....\n"); + + // Open our Nexus file in read only mode. + H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); + + const auto nSpec = spectra_.size(); + const auto nTofBins = tofBins_.size() - 1; + + std::vector countsBuffer; + countsBuffer.resize(nSpec * nTofBins); // HDF5 expects contiguous memory. This is a pain. + auto &&[counts, detectorCountsDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "counts"); + H5Dread(counts.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, countsBuffer.data()); + for (auto i = 0; i < nSpec; ++i) + for (auto j = 0; j < nTofBins; ++j) + gsl_histogram_se(detectorHistograms_[spectra_[i]], j); countsBuffer[i * nTofBins + j] +} + // Save key modified data back to the file bool NeXuSFile::saveModifiedData() { diff --git a/src/nexusFile.h b/src/nexusFile.h index be4922b..d5121a7 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -34,6 +34,8 @@ class NeXuSFile void loadEventData(); // Load start/end times void loadTimes(); + // Load detector counts from the file + void loadDetectorCounts(); // Save key modified data back to the file bool saveModifiedData(); @@ -53,6 +55,7 @@ class NeXuSFile std::vector frameOffsets_; std::vector tofBins_; std::map> monitorCounts_; + std::map> detectorCounts_; std::map detectorHistograms_; std::map> partitions_; diff --git a/src/partitionEventsIndividual.cpp b/src/partitionEventsIndividual.cpp index ae52ecb..172ee16 100644 --- a/src/partitionEventsIndividual.cpp +++ b/src/partitionEventsIndividual.cpp @@ -22,7 +22,7 @@ void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::vector> slices; auto sliceIt = slices.end(); - // Loop over input Nexus files + // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) { // Open the NeXuS file and get its event data diff --git a/src/partitionEventsSummed.cpp b/src/partitionEventsSummed.cpp index 2d90eee..7f64a87 100644 --- a/src/partitionEventsSummed.cpp +++ b/src/partitionEventsSummed.cpp @@ -23,7 +23,7 @@ void partitionEventsSummed(const std::vector &inputNeXusFiles, std: // Initialise the slice iterator and window slice / NeXuSFile references auto sliceIt = slices.begin(); - // Loop over input Nexus files + // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) { // Open the NeXuS file and get its event data diff --git a/src/processors.h b/src/processors.h index 2f0acd9..8b46fb3 100644 --- a/src/processors.h +++ b/src/processors.h @@ -13,6 +13,7 @@ namespace Processors enum class ProcessingMode { None, + DumpHistogram, DumpEvents, PartitionEventsIndividual, PartitionEventsSummed @@ -55,6 +56,8 @@ void saveSlices(std::vector> &slices); // Get Events std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorId, bool firstOnly = false); +// Dump histogram from specified spectrum +void dumpHistogram(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly = false); // Partition events into individual windows / slices void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, const Window &windowDefinition, int nSlices, double windowDelta); From 24481d011035ec7779d6296c9f70eeca50ca5e6c Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 15 Oct 2024 11:09:31 +0100 Subject: [PATCH 08/35] Add missing source file. --- src/dumpHistogram.cpp | 68 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 src/dumpHistogram.cpp diff --git a/src/dumpHistogram.cpp b/src/dumpHistogram.cpp new file mode 100644 index 0000000..95b174c --- /dev/null +++ b/src/dumpHistogram.cpp @@ -0,0 +1,68 @@ +#include "nexusFile.h" +#include "processors.h" +#include "window.h" +#include +#include + +namespace Processors +{ +// Dump histogram from specified spectrum +void dumpHistogram(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) +{ + /* + * Dump all events for the specified detector spectrum + */ + + fmt::print("Dumping histogram from detector spectrum {}...\n", spectrumId); + + std::map> eventMap; + std::optional lastSecondsSinceEpoch; + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the Nexus file ready for use + NeXuSFile nxs(nxsFileName); + nxs.loadDetectorCounts(); + + auto eventStart = 0, eventEnd = 0; + const auto &eventsPerFrame = nxs.eventsPerFrame(); + const auto &eventIndices = nxs.eventIndices(); + const auto &eventTimes = nxs.eventTimes(); + const auto &frameOffsets = nxs.frameOffsets(); + + // Loop over frames in the Nexus file + for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) + { + // Set new end event index and get zero for frame + eventEnd += eventsPerFrame[frameIndex]; + auto frameZero = frameOffsets[frameIndex]; + + for (auto k = eventStart; k < eventEnd; ++k) + { + if (eventIndices[k] == spectrumId) + { + auto eMicroSeconds = eventTimes[k]; + auto eSeconds = eMicroSeconds * 0.000001; + auto eSecondsSinceEpoch = eSeconds + frameZero + nxs.startSinceEpoch(); + if (lastSecondsSinceEpoch) + fmt::print("{:20.6f} {:20.10f} {:20.5f} {}\n", eMicroSeconds, eSeconds + frameZero, + eSecondsSinceEpoch, eSecondsSinceEpoch - *lastSecondsSinceEpoch); + else + fmt::print("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch); + eventMap[spectrumId].push_back(eSecondsSinceEpoch); + lastSecondsSinceEpoch = eSecondsSinceEpoch; + if (firstOnly) + return eventMap; + } + } + + // Update start event index + eventStart = eventEnd; + } + } + + return eventMap; +} + +} // namespace Processors From e40ab283b9338c825ea3321e85b81542e7f07a21 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 15 Oct 2024 11:46:22 +0100 Subject: [PATCH 09/35] Dump histogram working. Separate array creation and monitor loading. --- src/CMakeLists.txt | 23 +++---------- src/dumpHistogram.cpp | 52 ++++++----------------------- src/nexusFile.cpp | 76 +++++++++++++++++++++++++++++-------------- src/nexusFile.h | 7 ++-- 4 files changed, 72 insertions(+), 86 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 546ba1c..4d27f0a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,20 +1,7 @@ -add_library(nexusProcess - dumpEvents.cpp - dumpHistogram.cpp - nexusFile.cpp - partitionEventsIndividual.cpp - partitionEventsSummed.cpp - processCommon.cpp - window.cpp - nexusFile.h - processors.h - window.h -) +add_library(nexusProcess dumpEvents.cpp dumpHistogram.cpp nexusFile.cpp partitionEventsIndividual.cpp partitionEventsSummed + .cpp processCommon.cpp window.cpp nexusFile.h processors.h window.h) -target_include_directories(nexusProcess PRIVATE ${PROJECT_SOURCE_DIR}/src ${CONAN_INCLUDE_DIRS}) + target_include_directories(nexusProcess PRIVATE ${PROJECT_SOURCE_DIR} / src ${CONAN_INCLUDE_DIRS}) -if(CONAN) - target_link_libraries(nexusProcess PUBLIC CONAN_PKG::fmt) -else(CONAN) - target_link_libraries(nexusProcess PUBLIC fmt::fmt) -endif(CONAN) + if (CONAN) target_link_libraries(nexusProcess PUBLIC CONAN_PKG::fmt) else(CONAN) + target_link_libraries(nexusProcess PUBLIC fmt::fmt) endif(CONAN) diff --git a/src/dumpHistogram.cpp b/src/dumpHistogram.cpp index 95b174c..2426498 100644 --- a/src/dumpHistogram.cpp +++ b/src/dumpHistogram.cpp @@ -1,8 +1,7 @@ #include "nexusFile.h" #include "processors.h" -#include "window.h" #include -#include +#include namespace Processors { @@ -13,56 +12,25 @@ void dumpHistogram(const std::vector &inputNeXusFiles, int spectrum * Dump all events for the specified detector spectrum */ - fmt::print("Dumping histogram from detector spectrum {}...\n", spectrumId); - - std::map> eventMap; - std::optional lastSecondsSinceEpoch; + fmt::print("Dumping histogram for detector spectrum {}...\n", spectrumId); // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) { - // Open the Nexus file ready for use + // Open the NeXuS file and load in detector counts NeXuSFile nxs(nxsFileName); nxs.loadDetectorCounts(); - auto eventStart = 0, eventEnd = 0; - const auto &eventsPerFrame = nxs.eventsPerFrame(); - const auto &eventIndices = nxs.eventIndices(); - const auto &eventTimes = nxs.eventTimes(); - const auto &frameOffsets = nxs.frameOffsets(); - - // Loop over frames in the Nexus file - for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) + // Open the output file + std::ofstream output(fmt::format("{}.{}", nxsFileName, spectrumId).c_str()); + auto bin = 0; + const auto &counts = nxs.detectorCounts().at(spectrumId); + for (auto tof : nxs.tofBins()) { - // Set new end event index and get zero for frame - eventEnd += eventsPerFrame[frameIndex]; - auto frameZero = frameOffsets[frameIndex]; - - for (auto k = eventStart; k < eventEnd; ++k) - { - if (eventIndices[k] == spectrumId) - { - auto eMicroSeconds = eventTimes[k]; - auto eSeconds = eMicroSeconds * 0.000001; - auto eSecondsSinceEpoch = eSeconds + frameZero + nxs.startSinceEpoch(); - if (lastSecondsSinceEpoch) - fmt::print("{:20.6f} {:20.10f} {:20.5f} {}\n", eMicroSeconds, eSeconds + frameZero, - eSecondsSinceEpoch, eSecondsSinceEpoch - *lastSecondsSinceEpoch); - else - fmt::print("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch); - eventMap[spectrumId].push_back(eSecondsSinceEpoch); - lastSecondsSinceEpoch = eSecondsSinceEpoch; - if (firstOnly) - return eventMap; - } - } - - // Update start event index - eventStart = eventEnd; + output << fmt::format("{} {}\n", tof, counts[bin++]); } + output.close(); } - - return eventMap; } } // namespace Processors diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index a7a85dc..3b2b9bb 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -25,6 +25,33 @@ std::vector neXuSBasicPaths_ = {"/raw_data_1/title", NeXuSFile::NeXuSFile(std::string filename, bool loadEvents) : filename_(filename) { + // Open input NeXuS file in read only mode. + H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); + + // Read in detector spectra information + auto &&[spectraID, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); + spectra_.resize(spectraDimension); + H5Dread(spectraID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, spectra_.data()); + + // Read in TOF bin information. + auto &&[tofBinsID, tofBinsDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/monitor_1", "time_of_flight"); + tofBins_.resize(tofBinsDimension); + H5Dread(tofBinsID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBins_.data()); + + // Set up detector histograms and straight counts vectors + for (auto spec : spectra_) + { + // For event re-binning + detectorHistograms_[spec] = gsl_histogram_alloc(tofBins_.size() - 1); + gsl_histogram_set_ranges(detectorHistograms_[spec], tofBins_.data(), tofBins_.size()); + + // For histogram manipulation + detectorCounts_[spec].resize(tofBins_.size() - 1); + std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); + } + + input.close(); + if (loadEvents) { fmt::print("Loading event data from file '{}'...\n", filename_); @@ -35,6 +62,13 @@ NeXuSFile::NeXuSFile(std::string filename, bool loadEvents) : filename_(filename } } +NeXuSFile::~NeXuSFile() +{ + // Free GSL histograms + for (auto &[specId, histogram] : detectorHistograms_) + gsl_histogram_free(histogram); +} + /* * I/O */ @@ -93,27 +127,17 @@ void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) H5Pclose(ocpl_id); H5Pclose(lcpl_id); - // Read in detector spectra information - auto &&[spectraID, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); - spectra_.resize(spectraDimension); - H5Dread(spectraID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, spectra_.data()); - - // Read in TOF bin information. - auto &&[tofBinsID, tofBinsDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/monitor_1", "time_of_flight"); - tofBins_.resize(tofBinsDimension); - H5Dread(tofBinsID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBins_.data()); + input.close(); + output.close(); +} - // Set up detector histograms and straight counts vectors - for (auto spec : spectra_) - { - // For event re-binning - detectorHistograms_[spec] = gsl_histogram_alloc(tofBins_.size() - 1); - gsl_histogram_set_ranges(detectorHistograms_[spec], tofBins_.data(), tofBins_.size()); +// Load in monitor histograms +void NeXuSFile::loadMonitorCounts() +{ + printf("Load monitor counts...\n"); - // For histogram manipulation - detectorCounts_[spec].resize(tofBins_.size() - 1); - std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); - } + // Open our NeXuS file in read only mode. + H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); // Read in monitor data - start from index 1 and end when we fail to find the named dataset with this suffix auto i = 1; @@ -124,20 +148,20 @@ void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) if (monitorSpectrum.getId() <= 0) break; - monitorCounts_[i].resize(tofBinsDimension); + monitorCounts_[i].resize(tofBins_.size()); H5Dread(monitorSpectrum.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, monitorCounts_[i].data()); ++i; } - // Read in good frames - this will reflect our current monitor frame count since we copied those histograms in full + // Read in number of good frames - this will reflect our current monitor frame count since we copied those histograms in + // full auto &&[goodFramesID, goodFramesDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "good_frames"); auto goodFramesTemp = new int[(long int)goodFramesDimension]; H5Dread(goodFramesID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, goodFramesTemp); nMonitorFrames_ = goodFramesTemp[0]; input.close(); - output.close(); } // Load frame counts @@ -240,7 +264,7 @@ void NeXuSFile::loadTimes() // Load detector counts from the file void NeXuSFile::loadDetectorCounts() { - printf("Load detector events....\n"); + printf("Load detector counts....\n"); // Open our Nexus file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); @@ -253,8 +277,11 @@ void NeXuSFile::loadDetectorCounts() auto &&[counts, detectorCountsDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "counts"); H5Dread(counts.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, countsBuffer.data()); for (auto i = 0; i < nSpec; ++i) + { + auto &counts = detectorCounts_[spectra_[i]]; for (auto j = 0; j < nTofBins; ++j) - gsl_histogram_se(detectorHistograms_[spectra_[i]], j); countsBuffer[i * nTofBins + j] + counts[j] = countsBuffer[i * nTofBins + j]; + } } // Save key modified data back to the file @@ -311,6 +338,7 @@ const std::vector &NeXuSFile::eventsPerFrame() const { return eventsPerFram const std::vector &NeXuSFile::frameOffsets() const { return frameOffsets_; } const std::vector &NeXuSFile::tofBins() const { return tofBins_; } const std::map> &NeXuSFile::monitorCounts() const { return monitorCounts_; } +const std::map> &NeXuSFile::detectorCounts() const { return detectorCounts_; } std::map &NeXuSFile::detectorHistograms() { return detectorHistograms_; } const std::map> &NeXuSFile::partitions() const { return partitions_; } diff --git a/src/nexusFile.h b/src/nexusFile.h index d5121a7..b0a8e5f 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -10,7 +10,7 @@ class NeXuSFile { public: NeXuSFile(std::string filename = "", bool loadEvents = false); - ~NeXuSFile() = default; + ~NeXuSFile(); /* * I/O @@ -26,8 +26,10 @@ class NeXuSFile public: // Return filename std::string filename() const; - // Template basic paths from the referenceFile, and make ready for histogram binning + // Template basic paths from the referenceFile void templateFile(std::string referenceFile, std::string outputFile); + // Load in monitor histograms + void loadMonitorCounts(); // Load frame counts void loadFrameCounts(); // Load event data @@ -72,6 +74,7 @@ class NeXuSFile [[nodiscard]] const std::vector &frameOffsets() const; [[nodiscard]] const std::vector &tofBins() const; [[nodiscard]] const std::map> &monitorCounts() const; + [[nodiscard]] const std::map> &detectorCounts() const; std::map &detectorHistograms(); [[nodiscard]] const std::map> &partitions() const; From 629c01e269b0bc07de118485369a90bdeaed364d Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 15 Oct 2024 11:51:29 +0100 Subject: [PATCH 10/35] Rename to dump-detector. --- np.cpp | 8 +++---- src/CMakeLists.txt | 24 ++++++++++++++++----- src/{dumpHistogram.cpp => dumpDetector.cpp} | 2 +- src/nexusFile.cpp | 2 +- src/processors.h | 4 ++-- 5 files changed, 27 insertions(+), 13 deletions(-) rename src/{dumpHistogram.cpp => dumpDetector.cpp} (90%) diff --git a/np.cpp b/np.cpp index 7b98010..ab56e32 100644 --- a/np.cpp +++ b/np.cpp @@ -68,11 +68,11 @@ int main(int argc, char **argv) "Dump all events from specified spectrum index") ->group("Processing"); app.add_option_function( - "--dump-histogram", + "--dump-detector", [&](int specId) { if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::DumpHistogram; + processingMode_ = Processors::ProcessingMode::DumpDetector; else { fmt::print("Error: Multiple processing modes given.\n"); @@ -131,8 +131,8 @@ int main(int argc, char **argv) case (Processors::ProcessingMode::DumpEvents): Processors::dumpEvents(inputFiles_, spectrumId_); break; - case (Processors::ProcessingMode::DumpHistogram): - Processors::dumpHistogram(inputFiles_, spectrumId_); + case (Processors::ProcessingMode::DumpDetector): + Processors::DumpDetector(inputFiles_, spectrumId_); break; case (Processors::ProcessingMode::PartitionEventsIndividual): case (Processors::ProcessingMode::PartitionEventsSummed): diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4d27f0a..23e3091 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,7 +1,21 @@ -add_library(nexusProcess dumpEvents.cpp dumpHistogram.cpp nexusFile.cpp partitionEventsIndividual.cpp partitionEventsSummed - .cpp processCommon.cpp window.cpp nexusFile.h processors.h window.h) +add_library( + nexusProcess + dumpDetector.cpp + dumpEvents.cpp + nexusFile.cpp + partitionEventsIndividual.cpp + partitionEventsSummed.cpp + processCommon.cpp + window.cpp + nexusFile.h + processors.h + window.h) - target_include_directories(nexusProcess PRIVATE ${PROJECT_SOURCE_DIR} / src ${CONAN_INCLUDE_DIRS}) +target_include_directories(nexusProcess PRIVATE ${PROJECT_SOURCE_DIR} / src + ${CONAN_INCLUDE_DIRS}) - if (CONAN) target_link_libraries(nexusProcess PUBLIC CONAN_PKG::fmt) else(CONAN) - target_link_libraries(nexusProcess PUBLIC fmt::fmt) endif(CONAN) +if(CONAN) + target_link_libraries(nexusProcess PUBLIC CONAN_PKG::fmt) +else(CONAN) + target_link_libraries(nexusProcess PUBLIC fmt::fmt) +endif(CONAN) diff --git a/src/dumpHistogram.cpp b/src/dumpDetector.cpp similarity index 90% rename from src/dumpHistogram.cpp rename to src/dumpDetector.cpp index 2426498..fd60190 100644 --- a/src/dumpHistogram.cpp +++ b/src/dumpDetector.cpp @@ -6,7 +6,7 @@ namespace Processors { // Dump histogram from specified spectrum -void dumpHistogram(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) +void DumpDetector(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) { /* * Dump all events for the specified detector spectrum diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 3b2b9bb..7c7717c 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -96,7 +96,7 @@ std::pair NeXuSFile::find1DDataset(H5::H5File file, H5std // Return filename std::string NeXuSFile::filename() const { return filename_; } -// Template basic paths from the referenceFile, and make ready for histogram binning +// Template basic paths from the referenceFile void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) { filename_ = outputFile; diff --git a/src/processors.h b/src/processors.h index 8b46fb3..783a69c 100644 --- a/src/processors.h +++ b/src/processors.h @@ -13,7 +13,7 @@ namespace Processors enum class ProcessingMode { None, - DumpHistogram, + DumpDetector, DumpEvents, PartitionEventsIndividual, PartitionEventsSummed @@ -57,7 +57,7 @@ void saveSlices(std::vector> &slices); std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorId, bool firstOnly = false); // Dump histogram from specified spectrum -void dumpHistogram(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly = false); +void DumpDetector(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly = false); // Partition events into individual windows / slices void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, const Window &windowDefinition, int nSlices, double windowDelta); From 2333d7d5dff5f9f758a55e3f441d281e1a479c42 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 15 Oct 2024 15:16:29 +0100 Subject: [PATCH 11/35] Add --dump-monitor. --- np.cpp | 20 +++++++++++++++++++- src/CMakeLists.txt | 1 + src/dumpDetector.cpp | 6 +++--- src/dumpMonitor.cpp | 36 ++++++++++++++++++++++++++++++++++++ src/nexusFile.cpp | 5 ++++- src/processors.h | 5 ++++- 6 files changed, 67 insertions(+), 6 deletions(-) create mode 100644 src/dumpMonitor.cpp diff --git a/np.cpp b/np.cpp index ab56e32..9dd0032 100644 --- a/np.cpp +++ b/np.cpp @@ -80,7 +80,22 @@ int main(int argc, char **argv) } spectrumId_ = specId; }, - "Dump specified histogram index") + "Dump specified detector histogram") + ->group("Processing"); + app.add_option_function( + "--dump-monitor", + [&](int specId) + { + if (processingMode_ == Processors::ProcessingMode::None) + processingMode_ = Processors::ProcessingMode::DumpMonitor; + else + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + spectrumId_ = specId; + }, + "Dump specified monitor histogram") ->group("Processing"); app.add_flag_callback( "--summed", @@ -134,6 +149,9 @@ int main(int argc, char **argv) case (Processors::ProcessingMode::DumpDetector): Processors::DumpDetector(inputFiles_, spectrumId_); break; + case (Processors::ProcessingMode::DumpMonitor): + Processors::DumpMonitor(inputFiles_, spectrumId_); + break; case (Processors::ProcessingMode::PartitionEventsIndividual): case (Processors::ProcessingMode::PartitionEventsSummed): // Sanity check diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 23e3091..449c104 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,7 @@ add_library( nexusProcess dumpDetector.cpp dumpEvents.cpp + dumpMonitor.cpp nexusFile.cpp partitionEventsIndividual.cpp partitionEventsSummed.cpp diff --git a/src/dumpDetector.cpp b/src/dumpDetector.cpp index fd60190..f1d849c 100644 --- a/src/dumpDetector.cpp +++ b/src/dumpDetector.cpp @@ -5,11 +5,10 @@ namespace Processors { -// Dump histogram from specified spectrum void DumpDetector(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) { /* - * Dump all events for the specified detector spectrum + * Dump histograms for the specified detector spectrum */ fmt::print("Dumping histogram for detector spectrum {}...\n", spectrumId); @@ -22,7 +21,8 @@ void DumpDetector(const std::vector &inputNeXusFiles, int spectrumI nxs.loadDetectorCounts(); // Open the output file - std::ofstream output(fmt::format("{}.{}", nxsFileName, spectrumId).c_str()); + std::ofstream output(fmt::format("{}.det.{}", nxsFileName, spectrumId).c_str()); + output << "# TCB/usec Counts\n"; auto bin = 0; const auto &counts = nxs.detectorCounts().at(spectrumId); for (auto tof : nxs.tofBins()) diff --git a/src/dumpMonitor.cpp b/src/dumpMonitor.cpp new file mode 100644 index 0000000..8637207 --- /dev/null +++ b/src/dumpMonitor.cpp @@ -0,0 +1,36 @@ +#include "nexusFile.h" +#include "processors.h" +#include +#include + +namespace Processors +{ +void DumpMonitor(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) +{ + /* + * Dump histograms for the specified monitor spectrum + */ + + fmt::print("Dumping histogram for monitor spectrum {}...\n", spectrumId); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and load in monitor counts + NeXuSFile nxs(nxsFileName); + nxs.loadMonitorCounts(); + + // Open the output file + std::ofstream output(fmt::format("{}.mon.{}", nxsFileName, spectrumId).c_str()); + output << "# TCB/usec Counts\n"; + auto bin = 0; + const auto &counts = nxs.monitorCounts().at(spectrumId); + for (auto tof : nxs.tofBins()) + { + output << fmt::format("{} {}\n", tof, counts[bin++]); + } + output.close(); + } +} + +} // namespace Processors diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 7c7717c..141c7d3 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -25,6 +25,8 @@ std::vector neXuSBasicPaths_ = {"/raw_data_1/title", NeXuSFile::NeXuSFile(std::string filename, bool loadEvents) : filename_(filename) { + fmt::print("Opening NeXuS file '{}'...\n", filename_); + // Open input NeXuS file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); @@ -32,6 +34,7 @@ NeXuSFile::NeXuSFile(std::string filename, bool loadEvents) : filename_(filename auto &&[spectraID, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); spectra_.resize(spectraDimension); H5Dread(spectraID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, spectra_.data()); + fmt::print("... total number of detector spectra is {}.\n", spectra_.size()); // Read in TOF bin information. auto &&[tofBinsID, tofBinsDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/monitor_1", "time_of_flight"); @@ -54,7 +57,7 @@ NeXuSFile::NeXuSFile(std::string filename, bool loadEvents) : filename_(filename if (loadEvents) { - fmt::print("Loading event data from file '{}'...\n", filename_); + fmt::print("... Loading event data...\n"); loadFrameCounts(); loadEventData(); loadTimes(); diff --git a/src/processors.h b/src/processors.h index 783a69c..661642c 100644 --- a/src/processors.h +++ b/src/processors.h @@ -15,6 +15,7 @@ enum class ProcessingMode None, DumpDetector, DumpEvents, + DumpMonitor, PartitionEventsIndividual, PartitionEventsSummed }; @@ -56,8 +57,10 @@ void saveSlices(std::vector> &slices); // Get Events std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorId, bool firstOnly = false); -// Dump histogram from specified spectrum +// Dump detector histogram void DumpDetector(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly = false); +// Dump monitor histogram +void DumpMonitor(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly = false); // Partition events into individual windows / slices void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, const Window &windowDefinition, int nSlices, double windowDelta); From d4c3dd66a1c3c2ce18ee7b74234bee50d6bd3a87 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 15 Oct 2024 15:50:22 +0100 Subject: [PATCH 12/35] Rename some variables for clarity, and index monitors and detectors from 1 from the user perspective. --- np.cpp | 14 ++++----- src/dumpDetector.cpp | 14 +++++---- src/dumpEvents.cpp | 6 ++-- src/dumpMonitor.cpp | 12 ++++---- src/nexusFile.cpp | 67 ++++++++++++++++++++++---------------------- src/nexusFile.h | 8 ++++-- src/processors.h | 6 ++-- 7 files changed, 66 insertions(+), 61 deletions(-) diff --git a/np.cpp b/np.cpp index 9dd0032..2f01145 100644 --- a/np.cpp +++ b/np.cpp @@ -27,7 +27,7 @@ int main(int argc, char **argv) // Number of slices to partition window in to int windowSlices_{1}; // Target spectrum for event get (optional) - int spectrumId_; + int targetIndex_; // Define and parse CLI arguments CLI::App app("NeXuS Processor (np), Copyright (C) 2024 Jared Swift and Tristan Youngs."); @@ -63,7 +63,7 @@ int main(int argc, char **argv) fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } - spectrumId_ = specId; + targetIndex_ = specId; }, "Dump all events from specified spectrum index") ->group("Processing"); @@ -78,7 +78,7 @@ int main(int argc, char **argv) fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } - spectrumId_ = specId; + targetIndex_ = specId; }, "Dump specified detector histogram") ->group("Processing"); @@ -93,7 +93,7 @@ int main(int argc, char **argv) fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } - spectrumId_ = specId; + targetIndex_ = specId; }, "Dump specified monitor histogram") ->group("Processing"); @@ -144,13 +144,13 @@ int main(int argc, char **argv) fmt::print("No processing mode specified. We are done.\n"); break; case (Processors::ProcessingMode::DumpEvents): - Processors::dumpEvents(inputFiles_, spectrumId_); + Processors::dumpEvents(inputFiles_, targetIndex_); break; case (Processors::ProcessingMode::DumpDetector): - Processors::DumpDetector(inputFiles_, spectrumId_); + Processors::DumpDetector(inputFiles_, targetIndex_); break; case (Processors::ProcessingMode::DumpMonitor): - Processors::DumpMonitor(inputFiles_, spectrumId_); + Processors::DumpMonitor(inputFiles_, targetIndex_); break; case (Processors::ProcessingMode::PartitionEventsIndividual): case (Processors::ProcessingMode::PartitionEventsSummed): diff --git a/src/dumpDetector.cpp b/src/dumpDetector.cpp index f1d849c..40f99c4 100644 --- a/src/dumpDetector.cpp +++ b/src/dumpDetector.cpp @@ -5,13 +5,13 @@ namespace Processors { -void DumpDetector(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) +void DumpDetector(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly) { /* - * Dump histograms for the specified detector spectrum + * Dump histograms for the specified detector index */ - fmt::print("Dumping histogram for detector spectrum {}...\n", spectrumId); + fmt::print("Dumping histogram for detector index {}...\n", detectorIndex); // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) @@ -20,12 +20,14 @@ void DumpDetector(const std::vector &inputNeXusFiles, int spectrumI NeXuSFile nxs(nxsFileName); nxs.loadDetectorCounts(); + const auto spectrumId = nxs.spectrumForDetector(detectorIndex); + // Open the output file - std::ofstream output(fmt::format("{}.det.{}", nxsFileName, spectrumId).c_str()); - output << "# TCB/usec Counts\n"; + std::ofstream output(fmt::format("{}.det.{}", nxsFileName, detectorIndex).c_str()); + output << fmt::format("# TCB/usec Counts [detector index {}, spectrum index = {}]\n", detectorIndex, spectrumId); auto bin = 0; const auto &counts = nxs.detectorCounts().at(spectrumId); - for (auto tof : nxs.tofBins()) + for (auto tof : nxs.tofBoundaries()) { output << fmt::format("{} {}\n", tof, counts[bin++]); } diff --git a/src/dumpEvents.cpp b/src/dumpEvents.cpp index a0022d7..f1513b4 100644 --- a/src/dumpEvents.cpp +++ b/src/dumpEvents.cpp @@ -7,13 +7,14 @@ namespace Processors { // Dump events from specified spectrum, returning seconds since epoch for each -std::map> dumpEvents(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) +std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorIndex, + bool firstOnly) { /* * Dump all events for the specified detector spectrum */ - fmt::print("Dumping all events from detector spectrum {}...\n", spectrumId); + fmt::print("Dumping all events from detector index {}...\n", detectorIndex); std::map> eventMap; std::optional lastSecondsSinceEpoch; @@ -32,6 +33,7 @@ std::map> dumpEvents(const std::vector &in const auto &eventIndices = nxs.eventIndices(); const auto &eventTimes = nxs.eventTimes(); const auto &frameOffsets = nxs.frameOffsets(); + const auto spectrumId = nxs.spectrumForDetector(detectorIndex); // Loop over frames in the Nexus file for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) diff --git a/src/dumpMonitor.cpp b/src/dumpMonitor.cpp index 8637207..699e7dd 100644 --- a/src/dumpMonitor.cpp +++ b/src/dumpMonitor.cpp @@ -5,13 +5,13 @@ namespace Processors { -void DumpMonitor(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly) +void DumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex, bool firstOnly) { /* - * Dump histograms for the specified monitor spectrum + * Dump histograms for the specified monitor index */ - fmt::print("Dumping histogram for monitor spectrum {}...\n", spectrumId); + fmt::print("Dumping histogram for monitor index {}...\n", monitorIndex); // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) @@ -21,11 +21,11 @@ void DumpMonitor(const std::vector &inputNeXusFiles, int spectrumId nxs.loadMonitorCounts(); // Open the output file - std::ofstream output(fmt::format("{}.mon.{}", nxsFileName, spectrumId).c_str()); + std::ofstream output(fmt::format("{}.mon.{}", nxsFileName, monitorIndex).c_str()); output << "# TCB/usec Counts\n"; auto bin = 0; - const auto &counts = nxs.monitorCounts().at(spectrumId); - for (auto tof : nxs.tofBins()) + const auto &counts = nxs.monitorCounts().at(monitorIndex); + for (auto tof : nxs.tofBoundaries()) { output << fmt::format("{} {}\n", tof, counts[bin++]); } diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 141c7d3..0486184 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -31,25 +31,28 @@ NeXuSFile::NeXuSFile(std::string filename, bool loadEvents) : filename_(filename H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); // Read in detector spectra information - auto &&[spectraID, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); - spectra_.resize(spectraDimension); - H5Dread(spectraID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, spectra_.data()); - fmt::print("... total number of detector spectra is {}.\n", spectra_.size()); - - // Read in TOF bin information. - auto &&[tofBinsID, tofBinsDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/monitor_1", "time_of_flight"); - tofBins_.resize(tofBinsDimension); - H5Dread(tofBinsID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBins_.data()); + auto &&[detSpecIndices, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); + detectorSpectrumIndices_.resize(spectraDimension); + H5Dread(detSpecIndices.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, detectorSpectrumIndices_.data()); + fmt::print("... total number of detector spectra is {}.\n", detectorSpectrumIndices_.size()); + nMonitorSpectra_ = detectorSpectrumIndices_.front() - 1; + fmt::print("... inferred number of monitor spectra is {}.\n", nMonitorSpectra_); + + // Read in TOF boundary information. + auto &&[tofBoundariesID, tofBoundariesDimension] = + NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "time_of_flight"); + tofBoundaries_.resize(tofBoundariesDimension); + H5Dread(tofBoundariesID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBoundaries_.data()); // Set up detector histograms and straight counts vectors - for (auto spec : spectra_) + for (auto spec : detectorSpectrumIndices_) { // For event re-binning - detectorHistograms_[spec] = gsl_histogram_alloc(tofBins_.size() - 1); - gsl_histogram_set_ranges(detectorHistograms_[spec], tofBins_.data(), tofBins_.size()); + detectorHistograms_[spec] = gsl_histogram_alloc(tofBoundaries_.size() - 1); + gsl_histogram_set_ranges(detectorHistograms_[spec], tofBoundaries_.data(), tofBoundaries_.size()); // For histogram manipulation - detectorCounts_[spec].resize(tofBins_.size() - 1); + detectorCounts_[spec].resize(tofBoundaries_.size() - 1); std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); } @@ -143,18 +146,13 @@ void NeXuSFile::loadMonitorCounts() H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); // Read in monitor data - start from index 1 and end when we fail to find the named dataset with this suffix - auto i = 1; - while (true) + for (auto i = 0; i > nMonitorSpectra_; ++i) { auto &&[monitorSpectrum, monitorSpectrumDimension] = - NeXuSFile::find1DDataset(input, "/raw_data_1/monitor_" + std::to_string(i), "data"); - if (monitorSpectrum.getId() <= 0) - break; + NeXuSFile::find1DDataset(input, "/raw_data_1/monitor_" + std::to_string(i + 1), "data"); - monitorCounts_[i].resize(tofBins_.size()); + monitorCounts_[i].resize(tofBoundaries_.size()); H5Dread(monitorSpectrum.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, monitorCounts_[i].data()); - - ++i; } // Read in number of good frames - this will reflect our current monitor frame count since we copied those histograms in @@ -272,18 +270,18 @@ void NeXuSFile::loadDetectorCounts() // Open our Nexus file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); - const auto nSpec = spectra_.size(); - const auto nTofBins = tofBins_.size() - 1; + const auto nSpec = detectorSpectrumIndices_.size(); + const auto nTOFBins = tofBoundaries_.size() - 1; std::vector countsBuffer; - countsBuffer.resize(nSpec * nTofBins); // HDF5 expects contiguous memory. This is a pain. + countsBuffer.resize(nSpec * nTOFBins); // HDF5 expects contiguous memory. This is a pain. auto &&[counts, detectorCountsDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "counts"); H5Dread(counts.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, countsBuffer.data()); for (auto i = 0; i < nSpec; ++i) { - auto &counts = detectorCounts_[spectra_[i]]; - for (auto j = 0; j < nTofBins; ++j) - counts[j] = countsBuffer[i * nTofBins + j]; + auto &counts = detectorCounts_[detectorSpectrumIndices_[i]]; + for (auto j = 0; j < nTOFBins; ++j) + counts[j] = countsBuffer[i * nTOFBins + j]; } } @@ -308,13 +306,13 @@ bool NeXuSFile::saveModifiedData() } // Write detector counts - const auto nSpec = spectra_.size(); - const auto nTofBins = tofBins_.size() - 1; + const auto nSpec = detectorSpectrumIndices_.size(); + const auto ntofBoundaries = tofBoundaries_.size() - 1; - auto *countsBuffer = new int[nSpec * nTofBins]; // HDF5 expects contiguous memory. This is a pain. + auto *countsBuffer = new int[nSpec * ntofBoundaries]; // HDF5 expects contiguous memory. This is a pain. for (auto i = 0; i < nSpec; ++i) - for (auto j = 0; j < nTofBins; ++j) - countsBuffer[i * nTofBins + j] = gsl_histogram_get(detectorHistograms_[spectra_[i]], j); + for (auto j = 0; j < ntofBoundaries; ++j) + countsBuffer[i * ntofBoundaries + j] = gsl_histogram_get(detectorHistograms_[detectorSpectrumIndices_[i]], j); auto &&[counts, detectorCountsDimension] = NeXuSFile::find1DDataset(output, "raw_data_1/detector_1", "counts"); counts.write(countsBuffer, H5::PredType::STD_I32LE); @@ -339,7 +337,8 @@ const std::vector &NeXuSFile::eventIndices() const { return eventIndices_; const std::vector &NeXuSFile::eventTimes() const { return eventTimes_; } const std::vector &NeXuSFile::eventsPerFrame() const { return eventsPerFrame_; } const std::vector &NeXuSFile::frameOffsets() const { return frameOffsets_; } -const std::vector &NeXuSFile::tofBins() const { return tofBins_; } +const std::vector &NeXuSFile::tofBoundaries() const { return tofBoundaries_; } +const int NeXuSFile::spectrumForDetector(int detectorId) const { return detectorSpectrumIndices_[detectorId - 1]; } const std::map> &NeXuSFile::monitorCounts() const { return monitorCounts_; } const std::map> &NeXuSFile::detectorCounts() const { return detectorCounts_; } std::map &NeXuSFile::detectorHistograms() { return detectorHistograms_; } @@ -368,7 +367,7 @@ void NeXuSFile::scaleMonitors(double factor) void NeXuSFile::scaleDetectors(double factor) { auto oldSum = 0, newSum = 0; - for (auto i : spectra_) + for (auto i : detectorSpectrumIndices_) { oldSum += gsl_histogram_sum(detectorHistograms_[i]); gsl_histogram_scale(detectorHistograms_[i], factor); diff --git a/src/nexusFile.h b/src/nexusFile.h index b0a8e5f..01b5f90 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -45,7 +45,8 @@ class NeXuSFile * Data */ private: - std::vector spectra_; + std::vector detectorSpectrumIndices_; + int nMonitorSpectra_{0}; int nMonitorFrames_{0}; int nDetectorFrames_{0}; int nGoodFrames_{0}; @@ -55,7 +56,7 @@ class NeXuSFile std::vector eventTimes_; std::vector eventsPerFrame_; std::vector frameOffsets_; - std::vector tofBins_; + std::vector tofBoundaries_; std::map> monitorCounts_; std::map> detectorCounts_; std::map detectorHistograms_; @@ -72,7 +73,8 @@ class NeXuSFile [[nodiscard]] const std::vector &eventTimes() const; [[nodiscard]] const std::vector &eventsPerFrame() const; [[nodiscard]] const std::vector &frameOffsets() const; - [[nodiscard]] const std::vector &tofBins() const; + [[nodiscard]] const std::vector &tofBoundaries() const; + [[nodiscard]] const int spectrumForDetector(int detectorId) const; [[nodiscard]] const std::map> &monitorCounts() const; [[nodiscard]] const std::map> &detectorCounts() const; std::map &detectorHistograms(); diff --git a/src/processors.h b/src/processors.h index 661642c..e19e5bb 100644 --- a/src/processors.h +++ b/src/processors.h @@ -55,12 +55,12 @@ void saveSlices(std::vector> &slices); */ // Get Events -std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorId, +std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly = false); // Dump detector histogram -void DumpDetector(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly = false); +void DumpDetector(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly = false); // Dump monitor histogram -void DumpMonitor(const std::vector &inputNeXusFiles, int spectrumId, bool firstOnly = false); +void DumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex, bool firstOnly = false); // Partition events into individual windows / slices void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, const Window &windowDefinition, int nSlices, double windowDelta); From 216d1701d50a44f949f08fc7a10cb83259ca52a9 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Tue, 15 Oct 2024 16:37:26 +0100 Subject: [PATCH 13/35] Another variable rename. --- np.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/np.cpp b/np.cpp index 2f01145..153d83e 100644 --- a/np.cpp +++ b/np.cpp @@ -54,7 +54,7 @@ int main(int argc, char **argv) // -- Processing Modes app.add_option_function( "--dump-events", - [&](int specId) + [&](int id) { if (processingMode_ == Processors::ProcessingMode::None) processingMode_ = Processors::ProcessingMode::DumpEvents; @@ -63,13 +63,13 @@ int main(int argc, char **argv) fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } - targetIndex_ = specId; + targetIndex_ = id; }, "Dump all events from specified spectrum index") ->group("Processing"); app.add_option_function( "--dump-detector", - [&](int specId) + [&](int id) { if (processingMode_ == Processors::ProcessingMode::None) processingMode_ = Processors::ProcessingMode::DumpDetector; @@ -78,13 +78,13 @@ int main(int argc, char **argv) fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } - targetIndex_ = specId; + targetIndex_ = id; }, "Dump specified detector histogram") ->group("Processing"); app.add_option_function( "--dump-monitor", - [&](int specId) + [&](int id) { if (processingMode_ == Processors::ProcessingMode::None) processingMode_ = Processors::ProcessingMode::DumpMonitor; @@ -93,7 +93,7 @@ int main(int argc, char **argv) fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } - targetIndex_ = specId; + targetIndex_ = id; }, "Dump specified monitor histogram") ->group("Processing"); From 1352c83244acf3df8bbb631f8888775fcd130915 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:19:18 +0100 Subject: [PATCH 14/35] Reorder logic. --- np.cpp | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/np.cpp b/np.cpp index 153d83e..a46c3b2 100644 --- a/np.cpp +++ b/np.cpp @@ -56,13 +56,12 @@ int main(int argc, char **argv) "--dump-events", [&](int id) { - if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::DumpEvents; - else + if (processingMode_ != Processors::ProcessingMode::None) { fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } + processingMode_ = Processors::ProcessingMode::DumpEvents; targetIndex_ = id; }, "Dump all events from specified spectrum index") @@ -71,13 +70,12 @@ int main(int argc, char **argv) "--dump-detector", [&](int id) { - if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::DumpDetector; - else + if (processingMode_ != Processors::ProcessingMode::None) { fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } + processingMode_ = Processors::ProcessingMode::DumpDetector; targetIndex_ = id; }, "Dump specified detector histogram") @@ -86,13 +84,12 @@ int main(int argc, char **argv) "--dump-monitor", [&](int id) { - if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::DumpMonitor; - else + if (processingMode_ != Processors::ProcessingMode::None) { fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } + processingMode_ = Processors::ProcessingMode::DumpMonitor; targetIndex_ = id; }, "Dump specified monitor histogram") @@ -101,13 +98,12 @@ int main(int argc, char **argv) "--summed", [&]() { - if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::PartitionEventsSummed; - else + if (processingMode_ != Processors::ProcessingMode::None) { fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } + processingMode_ = Processors::ProcessingMode::PartitionEventsSummed; }, "Sum windows / slices over all files and output NeXuS files per-slice") ->group("Processing"); @@ -115,13 +111,12 @@ int main(int argc, char **argv) "--individual", [&]() { - if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::PartitionEventsIndividual; - else + if (processingMode_ != Processors::ProcessingMode::None) { fmt::print("Error: Multiple processing modes given.\n"); throw(CLI::RuntimeError()); } + processingMode_ = Processors::ProcessingMode::PartitionEventsIndividual; }, "Output NeXuS files for each window / slice") ->group("Processing"); From b5b5a4d2338c908214c1d1d2ef23cc29b25c17c9 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:19:29 +0100 Subject: [PATCH 15/35] Dump events to file rather than stdout. --- src/dumpEvents.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/dumpEvents.cpp b/src/dumpEvents.cpp index f1513b4..fb8731e 100644 --- a/src/dumpEvents.cpp +++ b/src/dumpEvents.cpp @@ -3,6 +3,7 @@ #include "window.h" #include #include +#include namespace Processors { @@ -35,6 +36,9 @@ std::map> dumpEvents(const std::vector &in const auto &frameOffsets = nxs.frameOffsets(); const auto spectrumId = nxs.spectrumForDetector(detectorIndex); + std::ofstream output(fmt::format("{}.events.{}", nxsFileName, detectorIndex).c_str()); + output << fmt::format("# event(us) event(relative) epoch(s) delta(s)"); + // Loop over frames in the Nexus file for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) { @@ -50,10 +54,10 @@ std::map> dumpEvents(const std::vector &in auto eSeconds = eMicroSeconds * 0.000001; auto eSecondsSinceEpoch = eSeconds + frameZero + nxs.startSinceEpoch(); if (lastSecondsSinceEpoch) - fmt::print("{:20.6f} {:20.10f} {:20.5f} {}\n", eMicroSeconds, eSeconds + frameZero, + output << fmt::format("{:20.6f} {:20.10f} {:20.5f} {}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch, eSecondsSinceEpoch - *lastSecondsSinceEpoch); else - fmt::print("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch); + output << fmt::format("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch); eventMap[spectrumId].push_back(eSecondsSinceEpoch); lastSecondsSinceEpoch = eSecondsSinceEpoch; if (firstOnly) @@ -64,6 +68,8 @@ std::map> dumpEvents(const std::vector &in // Update start event index eventStart = eventEnd; } + + output.close(); } return eventMap; From c9fb32d61c256c4c755e66d304995cf9c28727b3 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:29:20 +0100 Subject: [PATCH 16/35] Always load in times. Add loadBasicData() function. --- np.cpp | 1 - src/dumpEvents.cpp | 8 +-- src/nexusFile.cpp | 148 ++++++++++++++++++++++----------------------- src/nexusFile.h | 4 +- 4 files changed, 78 insertions(+), 83 deletions(-) diff --git a/np.cpp b/np.cpp index a46c3b2..fed45a6 100644 --- a/np.cpp +++ b/np.cpp @@ -171,7 +171,6 @@ int main(int argc, char **argv) return 1; } NeXuSFile firstFile(inputFiles_.front()); - firstFile.loadTimes(); fmt::print("Window start time converted from relative to absolute time: {} => {} (= {} + {})\n", windowStartTime_, windowStartTime_ + firstFile.startSinceEpoch(), firstFile.startSinceEpoch(), windowStartTime_); diff --git a/src/dumpEvents.cpp b/src/dumpEvents.cpp index fb8731e..5fba0d1 100644 --- a/src/dumpEvents.cpp +++ b/src/dumpEvents.cpp @@ -2,8 +2,8 @@ #include "processors.h" #include "window.h" #include -#include #include +#include namespace Processors { @@ -26,7 +26,6 @@ std::map> dumpEvents(const std::vector &in // Open the Nexus file ready for use NeXuSFile nxs(nxsFileName); nxs.loadEventData(); - nxs.loadTimes(); fmt::print("... file '{}' has {} events...\n", nxsFileName, nxs.eventTimes().size()); auto eventStart = 0, eventEnd = 0; @@ -55,9 +54,10 @@ std::map> dumpEvents(const std::vector &in auto eSecondsSinceEpoch = eSeconds + frameZero + nxs.startSinceEpoch(); if (lastSecondsSinceEpoch) output << fmt::format("{:20.6f} {:20.10f} {:20.5f} {}\n", eMicroSeconds, eSeconds + frameZero, - eSecondsSinceEpoch, eSecondsSinceEpoch - *lastSecondsSinceEpoch); + eSecondsSinceEpoch, eSecondsSinceEpoch - *lastSecondsSinceEpoch); else - output << fmt::format("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch); + output << fmt::format("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, + eSecondsSinceEpoch); eventMap[spectrumId].push_back(eSecondsSinceEpoch); lastSecondsSinceEpoch = eSecondsSinceEpoch; if (firstOnly) diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 0486184..a048616 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -27,43 +27,13 @@ NeXuSFile::NeXuSFile(std::string filename, bool loadEvents) : filename_(filename { fmt::print("Opening NeXuS file '{}'...\n", filename_); - // Open input NeXuS file in read only mode. - H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); - - // Read in detector spectra information - auto &&[detSpecIndices, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); - detectorSpectrumIndices_.resize(spectraDimension); - H5Dread(detSpecIndices.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, detectorSpectrumIndices_.data()); - fmt::print("... total number of detector spectra is {}.\n", detectorSpectrumIndices_.size()); - nMonitorSpectra_ = detectorSpectrumIndices_.front() - 1; - fmt::print("... inferred number of monitor spectra is {}.\n", nMonitorSpectra_); - - // Read in TOF boundary information. - auto &&[tofBoundariesID, tofBoundariesDimension] = - NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "time_of_flight"); - tofBoundaries_.resize(tofBoundariesDimension); - H5Dread(tofBoundariesID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBoundaries_.data()); - - // Set up detector histograms and straight counts vectors - for (auto spec : detectorSpectrumIndices_) - { - // For event re-binning - detectorHistograms_[spec] = gsl_histogram_alloc(tofBoundaries_.size() - 1); - gsl_histogram_set_ranges(detectorHistograms_[spec], tofBoundaries_.data(), tofBoundaries_.size()); - - // For histogram manipulation - detectorCounts_[spec].resize(tofBoundaries_.size() - 1); - std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); - } - - input.close(); + loadBasicData(); if (loadEvents) { fmt::print("... Loading event data...\n"); loadFrameCounts(); loadEventData(); - loadTimes(); fmt::print("... file '{}' has {} goodframes and {} events...\n", filename_, nGoodFrames_, eventTimes_.size()); } } @@ -102,6 +72,77 @@ std::pair NeXuSFile::find1DDataset(H5::H5File file, H5std // Return filename std::string NeXuSFile::filename() const { return filename_; } +// Load basic information from the NeXuS file +void NeXuSFile::loadBasicData() +{ + // Open input NeXuS file in read only mode. + H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); + + // Read in start time in Unix time. + hid_t memType = H5Tcopy(H5T_C_S1); + H5Tset_size(memType, UCHAR_MAX); + char timeBuffer[UCHAR_MAX]; + int y = 0, M = 0, d = 0, h = 0, m = 0, s = 0; + + auto &&[startTimeID, startTimeDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "start_time"); + H5Dread(startTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, timeBuffer); + + sscanf(timeBuffer, "%d-%d-%dT%d:%d:%d", &y, &M, &d, &h, &m, &s); + std::tm stime = {0}; + stime.tm_year = y - 1900; + stime.tm_mon = M - 1; + stime.tm_mday = d; + stime.tm_hour = h; + stime.tm_min = m; + stime.tm_sec = s; + startSinceEpoch_ = (int)mktime(&stime); + fmt::print("... run started at {} ({} s since epoch).\n", timeBuffer, startSinceEpoch_); + + // Read in end time in Unix time. + auto &&[endTimeID, endTimeDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "end_time"); + H5Dread(endTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, timeBuffer); + + sscanf(timeBuffer, "%d-%d-%dT%d:%d:%d", &y, &M, &d, &h, &m, &s); + std::tm etime = {0}; + etime.tm_year = y - 1900; + etime.tm_mon = M - 1; + etime.tm_mday = d; + etime.tm_hour = h; + etime.tm_min = m; + etime.tm_sec = s; + endSinceEpoch_ = (int)mktime(&etime); + fmt::print("... run ended at {} ({} s since epoch).\n", timeBuffer, endSinceEpoch_); + fmt::print("... run duration was {} s.\n", endSinceEpoch_ - startSinceEpoch_); + + // Read in detector spectra information + auto &&[detSpecIndices, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); + detectorSpectrumIndices_.resize(spectraDimension); + H5Dread(detSpecIndices.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, detectorSpectrumIndices_.data()); + fmt::print("... total number of detector spectra is {}.\n", detectorSpectrumIndices_.size()); + nMonitorSpectra_ = detectorSpectrumIndices_.front() - 1; + fmt::print("... inferred number of monitor spectra is {}.\n", nMonitorSpectra_); + + // Read in TOF boundary information. + auto &&[tofBoundariesID, tofBoundariesDimension] = + NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "time_of_flight"); + tofBoundaries_.resize(tofBoundariesDimension); + H5Dread(tofBoundariesID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBoundaries_.data()); + + // Set up detector histograms and straight counts vectors + for (auto spec : detectorSpectrumIndices_) + { + // For event re-binning + detectorHistograms_[spec] = gsl_histogram_alloc(tofBoundaries_.size() - 1); + gsl_histogram_set_ranges(detectorHistograms_[spec], tofBoundaries_.data(), tofBoundaries_.size()); + + // For histogram manipulation + detectorCounts_[spec].resize(tofBoundaries_.size() - 1); + std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); + } + + input.close(); +} + // Template basic paths from the referenceFile void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) { @@ -217,51 +258,6 @@ void NeXuSFile::loadEventData() input.close(); } -// Load start/end times -void NeXuSFile::loadTimes() -{ - printf("Load times....\n"); - - // Open our NeXuS file in read only mode. - H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); - - // Read in start time in Unix time. - hid_t memType = H5Tcopy(H5T_C_S1); - H5Tset_size(memType, UCHAR_MAX); - char timeBuffer[UCHAR_MAX]; - int y = 0, M = 0, d = 0, h = 0, m = 0, s = 0; - - auto &&[startTimeID, startTimeDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "start_time"); - H5Dread(startTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, timeBuffer); - - sscanf(timeBuffer, "%d-%d-%dT%d:%d:%d", &y, &M, &d, &h, &m, &s); - std::tm stime = {0}; - stime.tm_year = y - 1900; - stime.tm_mon = M - 1; - stime.tm_mday = d; - stime.tm_hour = h; - stime.tm_min = m; - stime.tm_sec = s; - startSinceEpoch_ = (int)mktime(&stime); - - // Read in end time in Unix time. - auto &&[endTimeID, endTimeDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "end_time"); - H5Dread(endTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, timeBuffer); - - sscanf(timeBuffer, "%d-%d-%dT%d:%d:%d", &y, &M, &d, &h, &m, &s); - std::tm etime = {0}; - etime.tm_year = y - 1900; - etime.tm_mon = M - 1; - etime.tm_mday = d; - etime.tm_hour = h; - etime.tm_min = m; - etime.tm_sec = s; - - endSinceEpoch_ = (int)mktime(&etime); - - input.close(); -} - // Load detector counts from the file void NeXuSFile::loadDetectorCounts() { diff --git a/src/nexusFile.h b/src/nexusFile.h index 01b5f90..b33e2f9 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -26,6 +26,8 @@ class NeXuSFile public: // Return filename std::string filename() const; + // Load basic information from the NeXuS file + void loadBasicData(); // Template basic paths from the referenceFile void templateFile(std::string referenceFile, std::string outputFile); // Load in monitor histograms @@ -34,8 +36,6 @@ class NeXuSFile void loadFrameCounts(); // Load event data void loadEventData(); - // Load start/end times - void loadTimes(); // Load detector counts from the file void loadDetectorCounts(); // Save key modified data back to the file From 68ae7a5295543be918e8acaef1179b1d85e16708 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:35:47 +0100 Subject: [PATCH 17/35] Simplify ctor, always load goodframes in basic data. --- src/nexusFile.cpp | 43 +++++++++++-------------------- src/nexusFile.h | 4 +-- src/partitionEventsIndividual.cpp | 3 ++- src/partitionEventsSummed.cpp | 3 ++- 4 files changed, 20 insertions(+), 33 deletions(-) diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index a048616..718a7ee 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -23,18 +23,13 @@ std::vector neXuSBasicPaths_ = {"/raw_data_1/title", "/raw_data_1/monitor_9/data", "/raw_data_1/detector_1/counts"}; -NeXuSFile::NeXuSFile(std::string filename, bool loadEvents) : filename_(filename) +NeXuSFile::NeXuSFile(std::string filename) : filename_(filename) { - fmt::print("Opening NeXuS file '{}'...\n", filename_); - - loadBasicData(); - - if (loadEvents) + if (!filename_.empty()) { - fmt::print("... Loading event data...\n"); - loadFrameCounts(); - loadEventData(); - fmt::print("... file '{}' has {} goodframes and {} events...\n", filename_, nGoodFrames_, eventTimes_.size()); + fmt::print("Opening NeXuS file '{}'...\n", filename_); + + loadBasicData(); } } @@ -140,6 +135,13 @@ void NeXuSFile::loadBasicData() std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); } + // Read in good frames + auto &&[goodFramesID, goodFramesDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "good_frames"); + auto goodFramesTemp = new int[(long int)goodFramesDimension]; + H5Dread(goodFramesID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, goodFramesTemp); + nGoodFrames_ = goodFramesTemp[0]; + fmt::print("... there were {} good frames.\n", nGoodFrames_); + input.close(); } @@ -206,27 +208,10 @@ void NeXuSFile::loadMonitorCounts() input.close(); } -// Load frame counts -void NeXuSFile::loadFrameCounts() -{ - printf("Load frame counts...\n"); - - // Open our NeXuS file in read only mode. - H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); - - // Read in good frames - auto &&[goodFramesID, goodFramesDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "good_frames"); - auto goodFramesTemp = new int[(long int)goodFramesDimension]; - H5Dread(goodFramesID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, goodFramesTemp); - nGoodFrames_ = goodFramesTemp[0]; - - input.close(); -} - // Load event data void NeXuSFile::loadEventData() { - printf("Load event data...\n"); + printf("Loading event data...\n"); // Open our NeXuS file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); @@ -255,6 +240,8 @@ void NeXuSFile::loadEventData() frameOffsets_.resize(frameOffsetsDimension); H5Dread(frameOffsetsID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, frameOffsets_.data()); + fmt::print("... there are {} events...\n", eventTimes_.size()); + input.close(); } diff --git a/src/nexusFile.h b/src/nexusFile.h index b33e2f9..a4ebbf7 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -9,7 +9,7 @@ class NeXuSFile { public: - NeXuSFile(std::string filename = "", bool loadEvents = false); + NeXuSFile(std::string filename = ""); ~NeXuSFile(); /* @@ -32,8 +32,6 @@ class NeXuSFile void templateFile(std::string referenceFile, std::string outputFile); // Load in monitor histograms void loadMonitorCounts(); - // Load frame counts - void loadFrameCounts(); // Load event data void loadEventData(); // Load detector counts from the file diff --git a/src/partitionEventsIndividual.cpp b/src/partitionEventsIndividual.cpp index 172ee16..3bd62b2 100644 --- a/src/partitionEventsIndividual.cpp +++ b/src/partitionEventsIndividual.cpp @@ -26,7 +26,8 @@ void partitionEventsIndividual(const std::vector &inputNeXusFiles, for (auto &nxsFileName : inputNeXusFiles) { // Open the NeXuS file and get its event data - NeXuSFile nxs(nxsFileName, true); + NeXuSFile nxs(nxsFileName); + nxs.loadEventData(); const auto &eventsPerFrame = nxs.eventsPerFrame(); const auto &eventIndices = nxs.eventIndices(); diff --git a/src/partitionEventsSummed.cpp b/src/partitionEventsSummed.cpp index 7f64a87..88487de 100644 --- a/src/partitionEventsSummed.cpp +++ b/src/partitionEventsSummed.cpp @@ -27,7 +27,8 @@ void partitionEventsSummed(const std::vector &inputNeXusFiles, std: for (auto &nxsFileName : inputNeXusFiles) { // Open the NeXuS file and get its event data - NeXuSFile nxs(nxsFileName, true); + NeXuSFile nxs(nxsFileName); + nxs.loadEventData(); const auto &eventsPerFrame = nxs.eventsPerFrame(); const auto &eventIndices = nxs.eventIndices(); From f21d366546b58c6ef7d09a8b366032c10b610b30 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:40:16 +0100 Subject: [PATCH 18/35] Use at() to access array, fix output messages. --- src/dumpEvents.cpp | 2 +- src/nexusFile.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dumpEvents.cpp b/src/dumpEvents.cpp index 5fba0d1..5223e5b 100644 --- a/src/dumpEvents.cpp +++ b/src/dumpEvents.cpp @@ -26,7 +26,6 @@ std::map> dumpEvents(const std::vector &in // Open the Nexus file ready for use NeXuSFile nxs(nxsFileName); nxs.loadEventData(); - fmt::print("... file '{}' has {} events...\n", nxsFileName, nxs.eventTimes().size()); auto eventStart = 0, eventEnd = 0; const auto &eventsPerFrame = nxs.eventsPerFrame(); @@ -34,6 +33,7 @@ std::map> dumpEvents(const std::vector &in const auto &eventTimes = nxs.eventTimes(); const auto &frameOffsets = nxs.frameOffsets(); const auto spectrumId = nxs.spectrumForDetector(detectorIndex); + fmt::print("NeXuS file spectrum ID for detector index {} is {}.\n", detectorIndex, spectrumId); std::ofstream output(fmt::format("{}.events.{}", nxsFileName, detectorIndex).c_str()); output << fmt::format("# event(us) event(relative) epoch(s) delta(s)"); diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 718a7ee..2fa5731 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -321,7 +321,7 @@ const std::vector &NeXuSFile::eventTimes() const { return eventTimes_; } const std::vector &NeXuSFile::eventsPerFrame() const { return eventsPerFrame_; } const std::vector &NeXuSFile::frameOffsets() const { return frameOffsets_; } const std::vector &NeXuSFile::tofBoundaries() const { return tofBoundaries_; } -const int NeXuSFile::spectrumForDetector(int detectorId) const { return detectorSpectrumIndices_[detectorId - 1]; } +const int NeXuSFile::spectrumForDetector(int detectorId) const { return detectorSpectrumIndices_.at(detectorId - 1); } const std::map> &NeXuSFile::monitorCounts() const { return monitorCounts_; } const std::map> &NeXuSFile::detectorCounts() const { return detectorCounts_; } std::map &NeXuSFile::detectorHistograms() { return detectorHistograms_; } From 51ee8c09f19067d72f0524ee4617e7db38bcc6bd Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:41:49 +0100 Subject: [PATCH 19/35] Remove old files. --- src/modex.cpp | 475 -------------------------------------------------- src/modex.h | 41 ----- 2 files changed, 516 deletions(-) delete mode 100644 src/modex.cpp delete mode 100644 src/modex.h diff --git a/src/modex.cpp b/src/modex.cpp deleted file mode 100644 index 9a07fc6..0000000 --- a/src/modex.cpp +++ /dev/null @@ -1,475 +0,0 @@ -#include -#include -#include -#include -#include - -#include "modex.h" -#include "nexus.h" - -ModEx::ModEx(Config cfg_) : cfg(cfg_) -{ - diagnosticFile = std::ofstream(diagnosticPath, std::ofstream::out | std::ofstream::trunc); -} - -ModEx::~ModEx() -{ - if (diagnosticFile) - diagnosticFile.close(); -} - -bool ModEx::process() -{ - if (cfg.extrapolationMode == NONE) - { - epochPulses(cfg.rawPulses); - binPulsesToRuns(cfg.rawPulses); - totalPulses = cfg.rawPulses.size(); - for (auto &pulse : cfg.rawPulses) - { - printf(" Processing %f -> %f\n", pulse.start, pulse.end); - processPulse(pulse); - } - } - else if (cfg.extrapolationMode == FORWARDS_SUMMED) - { - // Our period ("sequence of pulses") contains exactly one pulse that we're - // interested in (enforced in config.cpp) Extrapolate defined pulse forwards - // in time to generate all pulses we need to bin from over the entire run - Period superPeriod; - if (!createSuperPeriod(superPeriod, cfg.summedNSlices)) - return false; - - // Determine slice multiplier - auto sliceDuration = superPeriod.pulses.front().end - superPeriod.pulses.front().start; - auto sliceMultiplier = 3600 / (int)sliceDuration; - printf("Slice multiplier for events is %i (based on a slice duration of %e)\n", sliceMultiplier, sliceDuration); - - // Template our output file(s) from the first run - std::map outputFiles; - for (auto n = 0; n < cfg.summedNSlices; ++n) - { - // std::string outpath = cfg.outputDir + "/sum-" + std::to_string((int) - // cfg.periodBegin); - std::stringstream ss; - ss << cfg.outputDir << "/sum-" << std::to_string((int)cfg.periodBegin); - if (!cfg.periodDefinition.pulseDefinitions.front().label.empty()) - ss << "-" << cfg.periodDefinition.pulseDefinitions.front().label; - if (cfg.summedNSlices > 1) - ss << "-" << std::setw(3) << std::setfill('0') << (n + 1); - ss << ".nxs"; - auto data = outputFiles.emplace(n, NeXuSFile(cfg.runs[0], ss.str())); - if (!data.first->second.createEmpty(cfg.nxsDefinitionPaths)) - return false; - } - - // Loop over input NeXuS files - for (auto &nxsFileName : cfg.runs) - { - // Open the Nexus file ready for use - NeXuSFile nxs(nxsFileName); - nxs.load(true); - printf("Nexus file has %i goodframes...\n", nxs.totalGoodFrames()); - - // Zero frame counters in all pulses - for (auto &p : superPeriod.pulses) - p.frameCounter = 0; - - // Get first and last pulses which this file might contribute to - auto beginPulseIt = - std::find_if(superPeriod.pulses.begin(), superPeriod.pulses.end(), - [&nxs](const auto &p) { return p.end > nxs.startSinceEpoch() && p.start < nxs.endSinceEpoch(); }); - if (beginPulseIt == superPeriod.pulses.end()) - { - printf("!!! No pulses fall into the time range of this file - moving " - "on to the next...\n"); - continue; - } - auto endPulseIt = std::find_if(superPeriod.pulses.begin(), superPeriod.pulses.end(), - [&nxs](const auto &p) { return p.start > nxs.endSinceEpoch(); }); - - // Loop over frames in the Nexus file - auto pulseIt = beginPulseIt; - auto eventStart = 0, eventEnd = 0; - const auto &eventsPerFrame = nxs.eventsPerFrame(); - const auto &eventIndices = nxs.eventIndices(); - const auto &events = nxs.events(); - const auto &frameOffsets = nxs.frameOffsets(); - - for (auto i = 0; i < nxs.eventsPerFrame().size(); ++i) - { - // Set new end event index and get zero for frame - eventEnd += eventsPerFrame[i]; - auto frameZero = frameOffsets[i] + nxs.startSinceEpoch(); - - // If the frame zero is less than the start time of the current pulse, - // move on - if (frameZero < pulseIt->start) - continue; - - // If the frame zero is less than the end time of the current pulse, bin - // the events. Otherwise, increase the pulse iterator - if (frameZero < pulseIt->end) - { - // Grab the destination datafile for this pulse and bin events - auto &destinationNexus = outputFiles[pulseIt->sliceIndex]; - auto &destinationHistograms = destinationNexus.detectorHistograms(); - for (int k = eventStart; k < eventEnd; ++k) - { - auto id = eventIndices[k]; - auto event = events[k]; - if (id > 0) - gsl_histogram_accumulate(destinationHistograms[id], event, sliceMultiplier); - } - - // Increment the goodframes counter for this pulse - pulseIt->frameCounter += sliceMultiplier; - } - else - ++pulseIt; - - // If we have no more pulses, we can stop processing frames - if (pulseIt == endPulseIt) - break; - - // Update start event index - eventStart = eventEnd; - } - - // For each pulse we just added to, increase the goodFrames, and monitors - // by the correct fractional amount - for (auto it = beginPulseIt; it < endPulseIt; ++it) - { - if (it->frameCounter == 0) - continue; - auto &destinationNexus = outputFiles[it->sliceIndex]; - destinationNexus.nProcessedGoodFrames() += it->frameCounter; - nxs.addMonitors((double)it->frameCounter / it->frameCounter, destinationNexus); - } - - // If we have no more pulses, we can stop processing nexus files - if (pulseIt == endPulseIt) - break; - } - - // Save output files - for (auto &data : outputFiles) - { - auto sliceIndex = data.first + 1; - auto &nexus = data.second; - - printf("Output nexus for slice %i has %i good frames in total.\n", sliceIndex, nexus.nProcessedGoodFrames()); - - if (!nexus.output(cfg.nxsDefinitionPaths)) - return false; - diagnosticFile << nexus.getOutpath() << " " << nexus.nProcessedGoodFrames() << std::endl; - std::cout << "Finished processing: " << nexus.getOutpath() << std::endl; - } - } - else - { - std::vector periods; - - extrapolatePeriods(periods); - binPeriodsToRuns(periods); - totalPulses = periods.size() * cfg.periodDefinition.pulseDefinitions.size(); - for (auto &period : periods) - { - if (period.isComplete()) - { - for (auto &pulse : period.pulses) - { - printf(" Processing %f -> %f\n", pulse.start, pulse.end); - processPulse(pulse); - } - } - else - { - currentPulse += cfg.periodDefinition.pulseDefinitions.size(); - progress = (double)currentPulse / (double)totalPulses; - progress *= 100; - diagnosticFile << "Period " << period.start << " to " << period.end << " was ignored (incomplete period)." - << std::endl; - } - } - } - return true; -} - -bool ModEx::extrapolatePeriods(std::vector &periods) -{ - - std::cout << "Extrapolating periods (in seconds since epoch)\n"; - NeXuSFile firstRunNXS(cfg.runs[0]); - firstRunNXS.load(); - NeXuSFile lastRunNXS(cfg.runs[cfg.runs.size() - 1]); - lastRunNXS.load(); - - const int expStart = firstRunNXS.startSinceEpoch(); - const int expEnd = lastRunNXS.endSinceEpoch(); - double startPeriod = double(expStart) + cfg.periodBegin; - double periodBegin = 0; - - // First period - std::vector firstPeriodPulses; - for (auto &p : cfg.periodDefinition.pulseDefinitions) - { - firstPeriodPulses.push_back(Pulse(p, startPeriod + p.periodOffset, startPeriod + p.periodOffset + p.duration)); - } - - periods.push_back( - Period(cfg.periodDefinition, startPeriod, startPeriod + cfg.periodDefinition.duration, firstPeriodPulses)); - - if (cfg.extrapolationMode == BACKWARDS || cfg.extrapolationMode == BI_DIRECTIONAL) - { - periodBegin = startPeriod - cfg.periodDefinition.duration; - while (periodBegin > expStart) - { - std::vector pulses; - for (auto &p : cfg.periodDefinition.pulseDefinitions) - { - pulses.push_back(Pulse(p, periodBegin + p.periodOffset, periodBegin + p.periodOffset + p.duration)); - } - periods.push_back(Period(cfg.periodDefinition, periodBegin, periodBegin + cfg.periodDefinition.duration, pulses)); - periodBegin -= cfg.periodDefinition.duration; - } - } - if (cfg.extrapolationMode == FORWARDS || cfg.extrapolationMode == BI_DIRECTIONAL) - { - periodBegin = startPeriod + cfg.periodDefinition.duration; - while (periodBegin < expEnd) - { - std::vector pulses; - for (auto &p : cfg.periodDefinition.pulseDefinitions) - { - pulses.push_back(Pulse(p, periodBegin + p.periodOffset, periodBegin + p.periodOffset + p.duration)); - } - periods.push_back(Period(cfg.periodDefinition, periodBegin, periodBegin + cfg.periodDefinition.duration, pulses)); - periodBegin += cfg.periodDefinition.duration; - } - } - - std::sort(periods.begin(), periods.end(), [](const Period a, const Period b) { return a.start < b.end; }); - - printf("Extrapolated %li periods:\n", periods.size()); - for (const auto &p : periods) - printf(" %f -> %f\n", p.start, p.end); - return true; -} - -bool ModEx::createSuperPeriod(Period &period, int nSlices) -{ - std::cout << "Extrapolating periods (in seconds since epoch)\n"; - NeXuSFile firstRunNXS(cfg.runs[0]); - if (!firstRunNXS.load()) - return false; - NeXuSFile lastRunNXS(cfg.runs[cfg.runs.size() - 1]); - if (!lastRunNXS.load()) - return false; - - // Get limiting times of experiment (seconds since epoch values) - const int expStart = firstRunNXS.startSinceEpoch(); - const int expEnd = lastRunNXS.endSinceEpoch(); - - // Set absolute start time of first period, equal to the experiment start plus - // first defined pulse time - double periodStart = double(expStart) + cfg.periodBegin; - - // We have enforced that only a single master PulseDefinition exists - grab it - auto &masterPulseDef = cfg.periodDefinition.pulseDefinitions.front(); - - // Create copies of the master pulse extrapolating forwards in time by the - // period duration - std::vector pulses; - auto nWholePulses = 0; - while (periodStart < expEnd) - { - // Check that the actual start of the pulse is before the experimental end - // time - if ((periodStart + masterPulseDef.periodOffset) >= expEnd) - break; - - // Take the pulse and split into the number of slices - auto sliceDelta = (double)masterPulseDef.duration / (double)nSlices; - auto sliceStart = periodStart + masterPulseDef.periodOffset; - for (auto n = 0; n < nSlices; ++n) - { - pulses.push_back(Pulse(masterPulseDef, sliceStart, sliceStart + sliceDelta, n)); - sliceStart += sliceDelta; - } - - ++nWholePulses; - periodStart += cfg.periodDefinition.duration; - } - - // Create the superperiod - auto superPeriodStart = double(expStart) + cfg.periodBegin; - auto superPeriodEnd = superPeriodStart + cfg.periodDefinition.duration * pulses.size(); - period = Period(cfg.periodDefinition, superPeriodStart, superPeriodEnd, pulses); - - printf("Extrapolated %li pulses into superperiod (from %i whole pulses split " - "into %i slices each):\n", - period.pulses.size(), nWholePulses, nSlices); - for (const auto &p : period.pulses) - printf(" %f -> %f\n", p.start, p.end); - return true; -} - -bool ModEx::binPeriodsToRuns(std::vector &periods) -{ - - for (auto &period : periods) - { - binPulsesToRuns(period.pulses); - } - - return true; -} - -bool ModEx::binPulsesToRuns(std::vector &pulses) -{ - - std::vector>> runBoundaries; - - for (const auto &run : cfg.runs) - { - NeXuSFile nxs(run); - nxs.load(); - runBoundaries.push_back(std::make_pair(run, std::make_pair(nxs.startSinceEpoch(), nxs.endSinceEpoch()))); - } - - for (auto &pulse : pulses) - { - for (int i = 0; i < runBoundaries.size(); ++i) - { - if ((pulse.start >= runBoundaries[i].second.first) && (pulse.start < runBoundaries[i].second.second)) - { - pulse.startRun = runBoundaries[i].first; - if (i < runBoundaries.size() - 1) - { - pulse.endRun = runBoundaries[i + 1].first; - } - else - { - pulse.endRun = pulse.startRun; - } - } - if ((pulse.end >= runBoundaries[i].second.first) && (pulse.end < runBoundaries[i].second.second)) - { - pulse.endRun = runBoundaries[i].first; - if (pulse.startRun.empty()) - { - if (i > 0) - { - pulse.startRun = runBoundaries[i - 1].first; - } - else - { - pulse.startRun = pulse.endRun; - } - } - break; - } - } - } - return true; -} - -bool ModEx::processPulse(Pulse &pulse) -{ - if (pulse.startRun == pulse.endRun) - { - std::string outpath; - if (!pulse.definition.label.empty()) - outpath = cfg.outputDir + "/" + std::to_string((int)pulse.start) + "-" + pulse.definition.label + ".nxs"; - else - outpath = cfg.outputDir + "/" + std::to_string((int)pulse.start) + ".nxs"; - NeXuSFile nxs = NeXuSFile(pulse.startRun, outpath); - if (!nxs.load(true)) - return false; - nxs.nProcessedGoodFrames() = nxs.createHistogram(pulse, nxs.startSinceEpoch()); - if (!nxs.reduceMonitors((double)nxs.nProcessedGoodFrames() / (double)nxs.totalGoodFrames())) - return false; - if (!nxs.output(cfg.nxsDefinitionPaths)) - return false; - diagnosticFile << outpath << " " << nxs.nProcessedGoodFrames() << std::endl; - progress = (double)currentPulse / (double)totalPulses; - progress *= 100; - ++currentPulse; - std::cout << "Finished processing: " << outpath << std::endl; - std::cout << "Progress: " << progress << "%" << std::endl; - return true; - } - else - { - std::string outpath; - if (!pulse.definition.label.empty()) - outpath = cfg.outputDir + "/" + std::to_string((int)pulse.start) + "-" + pulse.definition.label + ".nxs"; - else - outpath = cfg.outputDir + "/" + std::to_string((int)pulse.start) + ".nxs"; - std::cout << outpath << std::endl; - std::cout << pulse.startRun << std::endl; - std::cout << pulse.endRun << std::endl; - NeXuSFile startNxs = NeXuSFile(pulse.startRun); - NeXuSFile endNxs = NeXuSFile(pulse.endRun, outpath); - - if (!startNxs.load(true)) - return false; - if (!endNxs.load(true)) - return false; - - std::cout << pulse.start << " " << pulse.end << std::endl; - Pulse firstPulse(pulse.start, startNxs.endSinceEpoch()); - Pulse secondPulse(endNxs.startSinceEpoch(), pulse.end); - std::cout << "Pulse duration: " << secondPulse.end - firstPulse.start << std::endl; - std::map> monitors = startNxs.monitorCounts(); - std::cout << "Sum monitors" << std::endl; - for (auto &pair : monitors) - { - for (int i = 0; i < pair.second.size(); ++i) - { - pair.second[i] += endNxs.monitorCounts().at(pair.first)[i]; - } - } - - if (!startNxs.createHistogram(firstPulse, startNxs.startSinceEpoch())) - return false; - if (!endNxs.createHistogram(secondPulse, startNxs.detectorHistograms(), endNxs.startSinceEpoch())) - return false; - int availableGoodFrames = startNxs.totalGoodFrames() + endNxs.totalGoodFrames(); - int totalProcessedFrames = startNxs.nProcessedGoodFrames() + endNxs.nProcessedGoodFrames(); - if (!endNxs.reduceMonitors((double)totalProcessedFrames / (double)availableGoodFrames)) - return false; - if (!endNxs.output(cfg.nxsDefinitionPaths)) - return false; - diagnosticFile << outpath << " " << startNxs.nProcessedGoodFrames() + endNxs.nProcessedGoodFrames() << std::endl; - progress = (double)currentPulse / (double)totalPulses; - progress *= 100; - ++currentPulse; - std::cout << "Finished processing: " << outpath << std::endl; - std::cout << "Progress: " << progress << "%" << std::endl; - return true; - } -} - -bool ModEx::epochPulses(std::vector &pulses) -{ - - // Assume runs are ordered. - std::string firstRun = cfg.runs[0]; - // Load the first run. - NeXuSFile firstRunNXS(firstRun); - firstRunNXS.load(); - - // Apply offset. - - const int expStart = firstRunNXS.startSinceEpoch(); - - for (int i = 0; i < pulses.size(); ++i) - { - pulses[i].start += expStart; - pulses[i].end += expStart; - } - - return true; -} diff --git a/src/modex.h b/src/modex.h deleted file mode 100644 index b3dc250..0000000 --- a/src/modex.h +++ /dev/null @@ -1,41 +0,0 @@ -#ifndef MODEX_H -#define MODEX_H - -#include "config.h" -#include "nexus.h" -#include "period.h" -#include "pulse.h" -#include -#include - -class ModEx -{ - - private: - int currentPulse = 1; - - public: - Config cfg; - std::string out; - std::string dataDir; - std::string diagnosticPath = "modex.diagnostics"; - std::ofstream diagnosticFile; - int expStart; - double progress; - int totalPulses = 0; - - ModEx(Config cfg_); - ModEx() = default; - ~ModEx(); - - bool process(); - bool processPulse(Pulse &pulse); - bool epochPulses(std::vector &pulses); - bool extrapolatePeriods(std::vector &periods); - bool createSuperPeriod(Period &period, int nSlices); - bool processPeriod(Period &period); - bool binPulsesToRuns(std::vector &pulses); - bool binPeriodsToRuns(std::vector &periods); -}; - -#endif // MODEX_H From 9c98d427a5b6852069d51bdaa6ac2349402a099e Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:42:09 +0100 Subject: [PATCH 20/35] Capitalisation. --- src/dumpEvents.cpp | 4 ++-- src/nexusFile.cpp | 6 +++--- src/partitionEventsIndividual.cpp | 2 +- src/partitionEventsSummed.cpp | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/dumpEvents.cpp b/src/dumpEvents.cpp index 5223e5b..e9ae888 100644 --- a/src/dumpEvents.cpp +++ b/src/dumpEvents.cpp @@ -23,7 +23,7 @@ std::map> dumpEvents(const std::vector &in // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) { - // Open the Nexus file ready for use + // Open the NeXuS file ready for use NeXuSFile nxs(nxsFileName); nxs.loadEventData(); @@ -38,7 +38,7 @@ std::map> dumpEvents(const std::vector &in std::ofstream output(fmt::format("{}.events.{}", nxsFileName, detectorIndex).c_str()); output << fmt::format("# event(us) event(relative) epoch(s) delta(s)"); - // Loop over frames in the Nexus file + // Loop over frames in the NeXuS file for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) { // Set new end event index and get zero for frame diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 2fa5731..618096a 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -153,7 +153,7 @@ void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) // Open input NeXuS file in read only mode. H5::H5File input = H5::H5File(referenceFile, H5F_ACC_RDONLY); - // Create new Nexus file for output. + // Create new NeXuS file for output. H5::H5File output = H5::H5File(filename_, H5F_ACC_TRUNC); printf("Templating file '%s' to '%s'...\n", referenceFile.c_str(), filename_.c_str()); @@ -250,7 +250,7 @@ void NeXuSFile::loadDetectorCounts() { printf("Load detector counts....\n"); - // Open our Nexus file in read only mode. + // Open our NeXuS file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); const auto nSpec = detectorSpectrumIndices_.size(); @@ -271,7 +271,7 @@ void NeXuSFile::loadDetectorCounts() // Save key modified data back to the file bool NeXuSFile::saveModifiedData() { - // Open Nexus file in read/write mode. + // Open NeXuS file in read/write mode. H5::H5File output = H5::H5File(filename_, H5F_ACC_RDWR); // Write good frames diff --git a/src/partitionEventsIndividual.cpp b/src/partitionEventsIndividual.cpp index 3bd62b2..1a6260c 100644 --- a/src/partitionEventsIndividual.cpp +++ b/src/partitionEventsIndividual.cpp @@ -34,7 +34,7 @@ void partitionEventsIndividual(const std::vector &inputNeXusFiles, const auto &eventTimes = nxs.eventTimes(); const auto &frameOffsets = nxs.frameOffsets(); - // Loop over frames in the Nexus file + // Loop over frames in the NeXuS file auto eventStart = 0, eventEnd = 0; for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) { diff --git a/src/partitionEventsSummed.cpp b/src/partitionEventsSummed.cpp index 88487de..eb91a20 100644 --- a/src/partitionEventsSummed.cpp +++ b/src/partitionEventsSummed.cpp @@ -35,7 +35,7 @@ void partitionEventsSummed(const std::vector &inputNeXusFiles, std: const auto &eventTimes = nxs.eventTimes(); const auto &frameOffsets = nxs.frameOffsets(); - // Loop over frames in the Nexus file + // Loop over frames in the NeXuS file auto eventStart = 0, eventEnd = 0; for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) { From 91d61ca0d18c71fdd3743375f38e58b69dcc06e3 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:44:39 +0100 Subject: [PATCH 21/35] Process function capitalisation. --- np.cpp | 4 ++-- src/dumpDetector.cpp | 2 +- src/dumpMonitor.cpp | 2 +- src/processors.h | 4 ++-- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/np.cpp b/np.cpp index fed45a6..c8ffd26 100644 --- a/np.cpp +++ b/np.cpp @@ -142,10 +142,10 @@ int main(int argc, char **argv) Processors::dumpEvents(inputFiles_, targetIndex_); break; case (Processors::ProcessingMode::DumpDetector): - Processors::DumpDetector(inputFiles_, targetIndex_); + Processors::dumpDetector(inputFiles_, targetIndex_); break; case (Processors::ProcessingMode::DumpMonitor): - Processors::DumpMonitor(inputFiles_, targetIndex_); + Processors::dumpMonitor(inputFiles_, targetIndex_); break; case (Processors::ProcessingMode::PartitionEventsIndividual): case (Processors::ProcessingMode::PartitionEventsSummed): diff --git a/src/dumpDetector.cpp b/src/dumpDetector.cpp index 40f99c4..5da08a5 100644 --- a/src/dumpDetector.cpp +++ b/src/dumpDetector.cpp @@ -5,7 +5,7 @@ namespace Processors { -void DumpDetector(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly) +void dumpDetector(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly) { /* * Dump histograms for the specified detector index diff --git a/src/dumpMonitor.cpp b/src/dumpMonitor.cpp index 699e7dd..4d266ea 100644 --- a/src/dumpMonitor.cpp +++ b/src/dumpMonitor.cpp @@ -5,7 +5,7 @@ namespace Processors { -void DumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex, bool firstOnly) +void dumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex, bool firstOnly) { /* * Dump histograms for the specified monitor index diff --git a/src/processors.h b/src/processors.h index e19e5bb..16bb192 100644 --- a/src/processors.h +++ b/src/processors.h @@ -58,9 +58,9 @@ void saveSlices(std::vector> &slices); std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly = false); // Dump detector histogram -void DumpDetector(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly = false); +void dumpDetector(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly = false); // Dump monitor histogram -void DumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex, bool firstOnly = false); +void dumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex, bool firstOnly = false); // Partition events into individual windows / slices void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, const Window &windowDefinition, int nSlices, double windowDelta); From 0e15ef7c4d55390b96631bb41435a661fbcb7841 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:52:49 +0100 Subject: [PATCH 22/35] Display basic info on files if no processing mode is selected. --- np.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/np.cpp b/np.cpp index c8ffd26..ebde104 100644 --- a/np.cpp +++ b/np.cpp @@ -136,7 +136,9 @@ int main(int argc, char **argv) switch (processingMode_) { case (Processors::ProcessingMode::None): - fmt::print("No processing mode specified. We are done.\n"); + fmt::print("No processing mode specified. Basic data from files will be shown.\n"); + for (const auto &file : inputFiles_) + NeXuSFile nxs(file); break; case (Processors::ProcessingMode::DumpEvents): Processors::dumpEvents(inputFiles_, targetIndex_); From 278eaf431ffe67bbad60a5e09fdc5a38ddb7d49d Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:53:03 +0100 Subject: [PATCH 23/35] Rename var, load and display run title. --- src/nexusFile.cpp | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 618096a..1ab43c2 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -70,19 +70,25 @@ std::string NeXuSFile::filename() const { return filename_; } // Load basic information from the NeXuS file void NeXuSFile::loadBasicData() { + hid_t memType = H5Tcopy(H5T_C_S1); + H5Tset_size(memType, UCHAR_MAX); + char charBuffer[UCHAR_MAX]; + // Open input NeXuS file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); + // Get the run title + auto &&[titleID, titleDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "title"); + H5Dread(titleID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, charBuffer); + fmt::print("... run title was '{}'\n", charBuffer); + // Read in start time in Unix time. - hid_t memType = H5Tcopy(H5T_C_S1); - H5Tset_size(memType, UCHAR_MAX); - char timeBuffer[UCHAR_MAX]; int y = 0, M = 0, d = 0, h = 0, m = 0, s = 0; auto &&[startTimeID, startTimeDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "start_time"); - H5Dread(startTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, timeBuffer); + H5Dread(startTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, charBuffer); - sscanf(timeBuffer, "%d-%d-%dT%d:%d:%d", &y, &M, &d, &h, &m, &s); + sscanf(charBuffer, "%d-%d-%dT%d:%d:%d", &y, &M, &d, &h, &m, &s); std::tm stime = {0}; stime.tm_year = y - 1900; stime.tm_mon = M - 1; @@ -91,13 +97,13 @@ void NeXuSFile::loadBasicData() stime.tm_min = m; stime.tm_sec = s; startSinceEpoch_ = (int)mktime(&stime); - fmt::print("... run started at {} ({} s since epoch).\n", timeBuffer, startSinceEpoch_); + fmt::print("... run started at {} ({} s since epoch).\n", charBuffer, startSinceEpoch_); // Read in end time in Unix time. auto &&[endTimeID, endTimeDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "end_time"); - H5Dread(endTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, timeBuffer); + H5Dread(endTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, charBuffer); - sscanf(timeBuffer, "%d-%d-%dT%d:%d:%d", &y, &M, &d, &h, &m, &s); + sscanf(charBuffer, "%d-%d-%dT%d:%d:%d", &y, &M, &d, &h, &m, &s); std::tm etime = {0}; etime.tm_year = y - 1900; etime.tm_mon = M - 1; @@ -106,7 +112,7 @@ void NeXuSFile::loadBasicData() etime.tm_min = m; etime.tm_sec = s; endSinceEpoch_ = (int)mktime(&etime); - fmt::print("... run ended at {} ({} s since epoch).\n", timeBuffer, endSinceEpoch_); + fmt::print("... run ended at {} ({} s since epoch).\n", charBuffer, endSinceEpoch_); fmt::print("... run duration was {} s.\n", endSinceEpoch_ - startSinceEpoch_); // Read in detector spectra information From 053e66336fb692cd8c84b02a22e5dbea497f1a9a Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Wed, 16 Oct 2024 15:54:30 +0100 Subject: [PATCH 24/35] Clarify output. --- src/nexusFile.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 1ab43c2..788f724 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -113,7 +113,7 @@ void NeXuSFile::loadBasicData() etime.tm_sec = s; endSinceEpoch_ = (int)mktime(&etime); fmt::print("... run ended at {} ({} s since epoch).\n", charBuffer, endSinceEpoch_); - fmt::print("... run duration was {} s.\n", endSinceEpoch_ - startSinceEpoch_); + fmt::print("... literal run duration was {} s.\n", endSinceEpoch_ - startSinceEpoch_); // Read in detector spectra information auto &&[detSpecIndices, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); From d2b91fe9419fe04e0129877edf67f0bc2acaff39 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Thu, 17 Oct 2024 09:10:24 +0100 Subject: [PATCH 25/35] Clarify header. --- src/dumpEvents.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dumpEvents.cpp b/src/dumpEvents.cpp index e9ae888..c2cc3a4 100644 --- a/src/dumpEvents.cpp +++ b/src/dumpEvents.cpp @@ -36,7 +36,7 @@ std::map> dumpEvents(const std::vector &in fmt::print("NeXuS file spectrum ID for detector index {} is {}.\n", detectorIndex, spectrumId); std::ofstream output(fmt::format("{}.events.{}", nxsFileName, detectorIndex).c_str()); - output << fmt::format("# event(us) event(relative) epoch(s) delta(s)"); + output << fmt::format("# frame_offset(us) start_time_offset(s) epoch_offset(s) delta(s)"); // Loop over frames in the NeXuS file for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) From 6e5be1495f8e0c8bf0192838045e4e7543c5b027 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Thu, 17 Oct 2024 09:45:37 +0100 Subject: [PATCH 26/35] Add print events. --- np.cpp | 20 +++++++++++++++-- src/CMakeLists.txt | 2 +- src/dumpDetector.cpp | 4 ++-- src/dumpMonitor.cpp | 4 ++-- src/{dumpEvents.cpp => getEvents.cpp} | 32 +++++++++++++-------------- src/processors.h | 12 +++++----- 6 files changed, 45 insertions(+), 29 deletions(-) rename src/{dumpEvents.cpp => getEvents.cpp} (71%) diff --git a/np.cpp b/np.cpp index ebde104..2c90a17 100644 --- a/np.cpp +++ b/np.cpp @@ -64,7 +64,21 @@ int main(int argc, char **argv) processingMode_ = Processors::ProcessingMode::DumpEvents; targetIndex_ = id; }, - "Dump all events from specified spectrum index") + "Dump all events for specified detector index") + ->group("Processing"); + app.add_option_function( + "--print-events", + [&](int id) + { + if (processingMode_ != Processors::ProcessingMode::None) + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + processingMode_ = Processors::ProcessingMode::PrintEvents; + targetIndex_ = id; + }, + "Print all events for specified detector index") ->group("Processing"); app.add_option_function( "--dump-detector", @@ -141,7 +155,9 @@ int main(int argc, char **argv) NeXuSFile nxs(file); break; case (Processors::ProcessingMode::DumpEvents): - Processors::dumpEvents(inputFiles_, targetIndex_); + case (Processors::ProcessingMode::PrintEvents): + Processors::dumpEventTimesEpoch(inputFiles_, targetIndex_, + processingMode_ == Processors::ProcessingMode::PrintEvents); break; case (Processors::ProcessingMode::DumpDetector): Processors::dumpDetector(inputFiles_, targetIndex_); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 449c104..80c607b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,8 +1,8 @@ add_library( nexusProcess dumpDetector.cpp - dumpEvents.cpp dumpMonitor.cpp + getEvents.cpp nexusFile.cpp partitionEventsIndividual.cpp partitionEventsSummed.cpp diff --git a/src/dumpDetector.cpp b/src/dumpDetector.cpp index 5da08a5..55c7a10 100644 --- a/src/dumpDetector.cpp +++ b/src/dumpDetector.cpp @@ -5,7 +5,7 @@ namespace Processors { -void dumpDetector(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly) +void dumpDetector(const std::vector &inputNeXusFiles, int detectorIndex) { /* * Dump histograms for the specified detector index @@ -24,7 +24,7 @@ void dumpDetector(const std::vector &inputNeXusFiles, int detectorI // Open the output file std::ofstream output(fmt::format("{}.det.{}", nxsFileName, detectorIndex).c_str()); - output << fmt::format("# TCB/usec Counts [detector index {}, spectrum index = {}]\n", detectorIndex, spectrumId); + output << fmt::format("# TCB/us Counts [detector index {}, spectrum index = {}]\n", detectorIndex, spectrumId); auto bin = 0; const auto &counts = nxs.detectorCounts().at(spectrumId); for (auto tof : nxs.tofBoundaries()) diff --git a/src/dumpMonitor.cpp b/src/dumpMonitor.cpp index 4d266ea..929d751 100644 --- a/src/dumpMonitor.cpp +++ b/src/dumpMonitor.cpp @@ -5,7 +5,7 @@ namespace Processors { -void dumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex, bool firstOnly) +void dumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex) { /* * Dump histograms for the specified monitor index @@ -22,7 +22,7 @@ void dumpMonitor(const std::vector &inputNeXusFiles, int monitorInd // Open the output file std::ofstream output(fmt::format("{}.mon.{}", nxsFileName, monitorIndex).c_str()); - output << "# TCB/usec Counts\n"; + output << "# TCB/us Counts\n"; auto bin = 0; const auto &counts = nxs.monitorCounts().at(monitorIndex); for (auto tof : nxs.tofBoundaries()) diff --git a/src/dumpEvents.cpp b/src/getEvents.cpp similarity index 71% rename from src/dumpEvents.cpp rename to src/getEvents.cpp index c2cc3a4..8412d21 100644 --- a/src/dumpEvents.cpp +++ b/src/getEvents.cpp @@ -3,22 +3,18 @@ #include "window.h" #include #include +#include #include namespace Processors { -// Dump events from specified spectrum, returning seconds since epoch for each -std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorIndex, - bool firstOnly) +void dumpEventTimesEpoch(const std::vector &inputNeXusFiles, int detectorIndex, bool toStdOut) { /* - * Dump all events for the specified detector spectrum + * Get all events for the specified detector spectrum, returning seconds since epoch for each */ - fmt::print("Dumping all events from detector index {}...\n", detectorIndex); - - std::map> eventMap; - std::optional lastSecondsSinceEpoch; + fmt::print("Retrieving all events from detector index {}...\n", detectorIndex); // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) @@ -27,6 +23,7 @@ std::map> dumpEvents(const std::vector &in NeXuSFile nxs(nxsFileName); nxs.loadEventData(); + std::optional lastSecondsSinceEpoch; auto eventStart = 0, eventEnd = 0; const auto &eventsPerFrame = nxs.eventsPerFrame(); const auto &eventIndices = nxs.eventIndices(); @@ -35,8 +32,14 @@ std::map> dumpEvents(const std::vector &in const auto spectrumId = nxs.spectrumForDetector(detectorIndex); fmt::print("NeXuS file spectrum ID for detector index {} is {}.\n", detectorIndex, spectrumId); - std::ofstream output(fmt::format("{}.events.{}", nxsFileName, detectorIndex).c_str()); - output << fmt::format("# frame_offset(us) start_time_offset(s) epoch_offset(s) delta(s)"); + std::ofstream fileOutput; + + if (!toStdOut) + fileOutput.open(fmt::format("{}.events.{}", nxsFileName, detectorIndex).c_str()); + + std::ostream &output = toStdOut ? std::cout : fileOutput; + output << fmt::format("# {:20s} {:20s} {:20s} {}\n", "frame_offset(us)", "start_time_offset(s)", "epoch_offset(s)", + "delta(s)"); // Loop over frames in the NeXuS file for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) @@ -58,10 +61,8 @@ std::map> dumpEvents(const std::vector &in else output << fmt::format("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch); - eventMap[spectrumId].push_back(eSecondsSinceEpoch); + lastSecondsSinceEpoch = eSecondsSinceEpoch; - if (firstOnly) - return eventMap; } } @@ -69,10 +70,9 @@ std::map> dumpEvents(const std::vector &in eventStart = eventEnd; } - output.close(); + if (!toStdOut) + fileOutput.close(); } - - return eventMap; } } // namespace Processors diff --git a/src/processors.h b/src/processors.h index 16bb192..0765491 100644 --- a/src/processors.h +++ b/src/processors.h @@ -17,7 +17,8 @@ enum class ProcessingMode DumpEvents, DumpMonitor, PartitionEventsIndividual, - PartitionEventsSummed + PartitionEventsSummed, + PrintEvents }; // Processing Direction @@ -54,13 +55,12 @@ void saveSlices(std::vector> &slices); * Processors */ -// Get Events -std::map> dumpEvents(const std::vector &inputNeXusFiles, int detectorIndex, - bool firstOnly = false); +// Dump all events for the specified detector spectrum, returning seconds since epoch for each +void dumpEventTimesEpoch(const std::vector &inputNeXusFiles, int detectorIndex, bool toStdOut = false); // Dump detector histogram -void dumpDetector(const std::vector &inputNeXusFiles, int detectorIndex, bool firstOnly = false); +void dumpDetector(const std::vector &inputNeXusFiles, int detectorIndex); // Dump monitor histogram -void dumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex, bool firstOnly = false); +void dumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex); // Partition events into individual windows / slices void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, const Window &windowDefinition, int nSlices, double windowDelta); From 566c6468eed2a2804b75be97d900b1f6de5253c4 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 10:39:00 +0100 Subject: [PATCH 27/35] Clarify start time input arguments. --- np.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/np.cpp b/np.cpp index 2c90a17..354c0a3 100644 --- a/np.cpp +++ b/np.cpp @@ -35,12 +35,13 @@ int main(int argc, char **argv) app.add_option("-n,--name", windowName_, "Name of the window, used as a prefix to all output files") ->group("Window Definition"); app.add_option("-s,--start", windowStartTime_, - "Start time of the window (relative to first input file start time unless --absolute-start is given)") + "Absolute start time of the window (i.e. seconds since epoch). Specify --relative-time if using a time " + "relative to the run start.") ->group("Window Definition"); app.add_option("-w,--width", windowWidth_, "Window width in seconds)")->group("Window Definition"); - app.add_flag( - "--relative-start", relativeStartTime_, - "Flag that the given window start time is relative to the first run start time, not absolute (seconds since epoch)") + app.add_flag("--relative-start", relativeStartTime_, + "Flag that the given window start time is relative to the first run start time rather than absolute (seconds " + "since epoch)") ->group("Window Definition"); app.add_option("-d,--delta", windowDelta_, "Time between window occurrences, in seconds")->group("Window Definition"); app.add_option("--offset", windowOffset_, "Time after start time, in seconds, that the window begins.") From 29da49e1309d86ecb4fcddc8f6ea0132bf8fdd4d Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 11:39:39 +0100 Subject: [PATCH 28/35] Proper handling of data copying and file templating. --- np.cpp | 5 +- src/dumpDetector.cpp | 1 + src/dumpMonitor.cpp | 1 + src/getEvents.cpp | 1 + src/nexusFile.cpp | 169 ++++++++++++++++++++++-------- src/nexusFile.h | 19 ++-- src/partitionEventsIndividual.cpp | 1 + src/processCommon.cpp | 11 +- 8 files changed, 155 insertions(+), 53 deletions(-) diff --git a/np.cpp b/np.cpp index 354c0a3..550f6e6 100644 --- a/np.cpp +++ b/np.cpp @@ -153,7 +153,10 @@ int main(int argc, char **argv) case (Processors::ProcessingMode::None): fmt::print("No processing mode specified. Basic data from files will be shown.\n"); for (const auto &file : inputFiles_) - NeXuSFile nxs(file); + { + NeXuSFile nxs(file, true); + nxs.prepareSpectraSpace(true); + } break; case (Processors::ProcessingMode::DumpEvents): case (Processors::ProcessingMode::PrintEvents): diff --git a/src/dumpDetector.cpp b/src/dumpDetector.cpp index 55c7a10..b5c90a2 100644 --- a/src/dumpDetector.cpp +++ b/src/dumpDetector.cpp @@ -18,6 +18,7 @@ void dumpDetector(const std::vector &inputNeXusFiles, int detectorI { // Open the NeXuS file and load in detector counts NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); nxs.loadDetectorCounts(); const auto spectrumId = nxs.spectrumForDetector(detectorIndex); diff --git a/src/dumpMonitor.cpp b/src/dumpMonitor.cpp index 929d751..972c67d 100644 --- a/src/dumpMonitor.cpp +++ b/src/dumpMonitor.cpp @@ -18,6 +18,7 @@ void dumpMonitor(const std::vector &inputNeXusFiles, int monitorInd { // Open the NeXuS file and load in monitor counts NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); nxs.loadMonitorCounts(); // Open the output file diff --git a/src/getEvents.cpp b/src/getEvents.cpp index 8412d21..4b518ae 100644 --- a/src/getEvents.cpp +++ b/src/getEvents.cpp @@ -21,6 +21,7 @@ void dumpEventTimesEpoch(const std::vector &inputNeXusFiles, int de { // Open the NeXuS file ready for use NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); nxs.loadEventData(); std::optional lastSecondsSinceEpoch; diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index 788f724..f5d5051 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -9,10 +9,10 @@ std::vector neXuSBasicPaths_ = {"/raw_data_1/title", "/raw_data_1/user_1/name", "/raw_data_1/start_time", + "/raw_data_1/end_time", "/raw_data_1/good_frames", "/raw_data_1/raw_frames", "/raw_data_1/monitor_1/data", - "/raw_data_1/monitor_1/time_of_flight", "/raw_data_1/monitor_2/data", "/raw_data_1/monitor_3/data", "/raw_data_1/monitor_4/data", @@ -21,23 +21,87 @@ std::vector neXuSBasicPaths_ = {"/raw_data_1/title", "/raw_data_1/monitor_7/data", "/raw_data_1/monitor_8/data", "/raw_data_1/monitor_9/data", - "/raw_data_1/detector_1/counts"}; + "/raw_data_1/detector_1/counts", + "/raw_data_1/detector_1/spectrum_index", + "/raw_data_1/detector_1/time_of_flight"}; -NeXuSFile::NeXuSFile(std::string filename) : filename_(filename) +NeXuSFile::NeXuSFile(std::string filename, bool printInfo) : filename_(filename) { if (!filename_.empty()) { fmt::print("Opening NeXuS file '{}'...\n", filename_); - loadBasicData(); + loadBasicData(printInfo); } } -NeXuSFile::~NeXuSFile() +void NeXuSFile::operator=(NeXuSFile &source) { copy(source); } + +NeXuSFile::NeXuSFile(NeXuSFile &source) { copy(source); } + +NeXuSFile::NeXuSFile(NeXuSFile &&source) +{ + // Copy data, but don't deep copy histograms (just copy pointers) + copy(source, false); + + // Clear the detectorHistograms_ map here so we don't try to delete the histograms (which we just copied the pointers for) + source.detectorHistograms_.clear(); + + source.clear(); +} + +NeXuSFile::~NeXuSFile() { clear(); } + +// Copy data from specified source +void NeXuSFile::copy(const NeXuSFile &source, bool deepCopyHistograms) +{ + filename_ = source.filename_; + detectorSpectrumIndices_ = source.detectorSpectrumIndices_; + nMonitorSpectra_ = source.nMonitorSpectra_; + nMonitorFrames_ = source.nMonitorFrames_; + nDetectorFrames_ = source.nDetectorFrames_; + nGoodFrames_ = source.nGoodFrames_; + startSinceEpoch_ = source.startSinceEpoch_; + endSinceEpoch_ = source.endSinceEpoch_; + eventIndices_ = source.eventIndices_; + eventTimes_ = source.eventTimes_; + eventsPerFrame_ = source.eventsPerFrame_; + frameOffsets_ = source.frameOffsets_; + tofBoundaries_ = source.tofBoundaries_; + monitorCounts_ = source.monitorCounts_; + detectorCounts_ = source.detectorCounts_; + if (deepCopyHistograms) + { + detectorHistograms_.clear(); + for (auto &[specId, histogram] : detectorHistograms_) + detectorHistograms_[specId] = gsl_histogram_clone(histogram); + } + else + detectorHistograms_ = source.detectorHistograms_; +} + +// Clear all data and arrays +void NeXuSFile::clear() { + filename_ = {}; + detectorSpectrumIndices_.clear(); + nMonitorSpectra_ = 0; + nMonitorFrames_ = 0; + nDetectorFrames_ = 0; + nGoodFrames_ = 0; + startSinceEpoch_ = 0; + endSinceEpoch_ = 0; + eventIndices_.clear(); + eventTimes_.clear(); + eventsPerFrame_.clear(); + frameOffsets_.clear(); + tofBoundaries_.clear(); + monitorCounts_.clear(); + detectorCounts_.clear(); // Free GSL histograms for (auto &[specId, histogram] : detectorHistograms_) gsl_histogram_free(histogram); + detectorHistograms_.clear(); } /* @@ -68,7 +132,7 @@ std::pair NeXuSFile::find1DDataset(H5::H5File file, H5std std::string NeXuSFile::filename() const { return filename_; } // Load basic information from the NeXuS file -void NeXuSFile::loadBasicData() +void NeXuSFile::loadBasicData(bool printInfo) { hid_t memType = H5Tcopy(H5T_C_S1); H5Tset_size(memType, UCHAR_MAX); @@ -80,7 +144,8 @@ void NeXuSFile::loadBasicData() // Get the run title auto &&[titleID, titleDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "title"); H5Dread(titleID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, charBuffer); - fmt::print("... run title was '{}'\n", charBuffer); + if (printInfo) + fmt::print("... run title was '{}'\n", charBuffer); // Read in start time in Unix time. int y = 0, M = 0, d = 0, h = 0, m = 0, s = 0; @@ -97,7 +162,8 @@ void NeXuSFile::loadBasicData() stime.tm_min = m; stime.tm_sec = s; startSinceEpoch_ = (int)mktime(&stime); - fmt::print("... run started at {} ({} s since epoch).\n", charBuffer, startSinceEpoch_); + if (printInfo) + fmt::print("... run started at {} ({} s since epoch).\n", charBuffer, startSinceEpoch_); // Read in end time in Unix time. auto &&[endTimeID, endTimeDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "end_time"); @@ -112,33 +178,10 @@ void NeXuSFile::loadBasicData() etime.tm_min = m; etime.tm_sec = s; endSinceEpoch_ = (int)mktime(&etime); - fmt::print("... run ended at {} ({} s since epoch).\n", charBuffer, endSinceEpoch_); - fmt::print("... literal run duration was {} s.\n", endSinceEpoch_ - startSinceEpoch_); - - // Read in detector spectra information - auto &&[detSpecIndices, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); - detectorSpectrumIndices_.resize(spectraDimension); - H5Dread(detSpecIndices.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, detectorSpectrumIndices_.data()); - fmt::print("... total number of detector spectra is {}.\n", detectorSpectrumIndices_.size()); - nMonitorSpectra_ = detectorSpectrumIndices_.front() - 1; - fmt::print("... inferred number of monitor spectra is {}.\n", nMonitorSpectra_); - - // Read in TOF boundary information. - auto &&[tofBoundariesID, tofBoundariesDimension] = - NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "time_of_flight"); - tofBoundaries_.resize(tofBoundariesDimension); - H5Dread(tofBoundariesID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBoundaries_.data()); - - // Set up detector histograms and straight counts vectors - for (auto spec : detectorSpectrumIndices_) + if (printInfo) { - // For event re-binning - detectorHistograms_[spec] = gsl_histogram_alloc(tofBoundaries_.size() - 1); - gsl_histogram_set_ranges(detectorHistograms_[spec], tofBoundaries_.data(), tofBoundaries_.size()); - - // For histogram manipulation - detectorCounts_[spec].resize(tofBoundaries_.size() - 1); - std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); + fmt::print("... run ended at {} ({} s since epoch).\n", charBuffer, endSinceEpoch_); + fmt::print("... literal run duration was {} s.\n", endSinceEpoch_ - startSinceEpoch_); } // Read in good frames @@ -146,23 +189,23 @@ void NeXuSFile::loadBasicData() auto goodFramesTemp = new int[(long int)goodFramesDimension]; H5Dread(goodFramesID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, goodFramesTemp); nGoodFrames_ = goodFramesTemp[0]; - fmt::print("... there were {} good frames.\n", nGoodFrames_); + if (printInfo) + fmt::print("... there were {} good frames.\n", nGoodFrames_); input.close(); } -// Template basic paths from the referenceFile -void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) +// Template a new NeXusFile from that specified +void NeXuSFile::templateTo(std::string sourceFilename, std::string newFilename) { - filename_ = outputFile; - - // Open input NeXuS file in read only mode. - H5::H5File input = H5::H5File(referenceFile, H5F_ACC_RDONLY); + // Open this NeXuS file in read only mode. + H5::H5File input = H5::H5File(sourceFilename, H5F_ACC_RDONLY); // Create new NeXuS file for output. - H5::H5File output = H5::H5File(filename_, H5F_ACC_TRUNC); + H5::H5File output = H5::H5File(newFilename, H5F_ACC_TRUNC); + + printf("Templating file '%s' to '%s'...\n", sourceFilename.c_str(), newFilename.c_str()); - printf("Templating file '%s' to '%s'...\n", referenceFile.c_str(), filename_.c_str()); hid_t ocpl_id, lcpl_id; ocpl_id = H5Pcreate(H5P_OBJECT_COPY); if (ocpl_id < 0) @@ -186,6 +229,45 @@ void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) output.close(); } +// Prepare spectra storage, including loading TOF boundaries etc. +void NeXuSFile::prepareSpectraSpace(bool printInfo) +{ + // Open input NeXuS file in read only mode. + H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); + + // Read in detector spectra information + auto &&[detSpecIndices, spectraDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); + detectorSpectrumIndices_.resize(spectraDimension); + H5Dread(detSpecIndices.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, detectorSpectrumIndices_.data()); + nMonitorSpectra_ = detectorSpectrumIndices_.front() - 1; + + if (printInfo) + { + fmt::print("... total number of detector spectra is {}.\n", detectorSpectrumIndices_.size()); + fmt::print("... inferred number of monitor spectra is {}.\n", nMonitorSpectra_); + } + + // Read in TOF boundary information. + auto &&[tofBoundariesID, tofBoundariesDimension] = + NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "time_of_flight"); + tofBoundaries_.resize(tofBoundariesDimension); + H5Dread(tofBoundariesID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBoundaries_.data()); + + // Set up detector histograms and straight counts vectors + for (auto spec : detectorSpectrumIndices_) + { + // For event re-binning + detectorHistograms_[spec] = gsl_histogram_alloc(tofBoundaries_.size() - 1); + gsl_histogram_set_ranges(detectorHistograms_[spec], tofBoundaries_.data(), tofBoundaries_.size()); + + // For histogram manipulation + detectorCounts_[spec].resize(tofBoundaries_.size() - 1); + std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); + } + + input.close(); +} + // Load in monitor histograms void NeXuSFile::loadMonitorCounts() { @@ -331,7 +413,6 @@ const int NeXuSFile::spectrumForDetector(int detectorId) const { return detector const std::map> &NeXuSFile::monitorCounts() const { return monitorCounts_; } const std::map> &NeXuSFile::detectorCounts() const { return detectorCounts_; } std::map &NeXuSFile::detectorHistograms() { return detectorHistograms_; } -const std::map> &NeXuSFile::partitions() const { return partitions_; } /* * Manipulation diff --git a/src/nexusFile.h b/src/nexusFile.h index a4ebbf7..005c19a 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -9,8 +9,15 @@ class NeXuSFile { public: - NeXuSFile(std::string filename = ""); + NeXuSFile(std::string filename = "", bool printInfo = false); + void operator=(NeXuSFile &source); + NeXuSFile(NeXuSFile &source); + NeXuSFile(NeXuSFile &&source); ~NeXuSFile(); + // Clear all data and arrays + void clear(); + // Copy data from specified source + void copy(const NeXuSFile &source, bool deepCopyHistograms = true); /* * I/O @@ -27,9 +34,11 @@ class NeXuSFile // Return filename std::string filename() const; // Load basic information from the NeXuS file - void loadBasicData(); - // Template basic paths from the referenceFile - void templateFile(std::string referenceFile, std::string outputFile); + void loadBasicData(bool printInfo = false); + // Template a new NeXusFile from that specified + static void templateTo(std::string sourceFilename, std::string newFilename); + // Prepare spectra storage, including loading TOF boundaries etc. + void prepareSpectraSpace(bool printInfo = false); // Load in monitor histograms void loadMonitorCounts(); // Load event data @@ -58,7 +67,6 @@ class NeXuSFile std::map> monitorCounts_; std::map> detectorCounts_; std::map detectorHistograms_; - std::map> partitions_; public: [[nodiscard]] int nGoodFrames() const; @@ -76,7 +84,6 @@ class NeXuSFile [[nodiscard]] const std::map> &monitorCounts() const; [[nodiscard]] const std::map> &detectorCounts() const; std::map &detectorHistograms(); - [[nodiscard]] const std::map> &partitions() const; /* * Manipulation diff --git a/src/partitionEventsIndividual.cpp b/src/partitionEventsIndividual.cpp index 1a6260c..5a337ac 100644 --- a/src/partitionEventsIndividual.cpp +++ b/src/partitionEventsIndividual.cpp @@ -27,6 +27,7 @@ void partitionEventsIndividual(const std::vector &inputNeXusFiles, { // Open the NeXuS file and get its event data NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); nxs.loadEventData(); const auto &eventsPerFrame = nxs.eventsPerFrame(); diff --git a/src/processCommon.cpp b/src/processCommon.cpp index 51f6aa7..1fb7554 100644 --- a/src/processCommon.cpp +++ b/src/processCommon.cpp @@ -28,6 +28,9 @@ std::vector> prepareSlices(const Window &window, in fmt::print("Individual slice duration is {} seconds.\n", sliceDuration); } + // Open the source NeXuS file + NeXuSFile sourceNxs(templatingSourceFilename); + for (auto i = 0; i < nSlices; ++i) { std::stringstream outputFileName; @@ -39,9 +42,13 @@ std::vector> prepareSlices(const Window &window, in std::stringstream sliceName; sliceName << window.id() << i + 1; - auto &[newWin, nexus] = slices.emplace_back(Window(sliceName.str(), sliceStartTime, sliceDuration), NeXuSFile()); + // Template the NeXuS file + NeXuSFile::templateTo(templatingSourceFilename, outputFileName.str()); - nexus.templateFile(templatingSourceFilename, outputFileName.str()); + auto &[newWin, nexus] = + slices.emplace_back(Window(sliceName.str(), sliceStartTime, sliceDuration), NeXuSFile(outputFileName.str())); + nexus.loadBasicData(); + nexus.prepareSpectraSpace(true); sliceStartTime += sliceDuration; } From 8efa2169044bcfc488dfa4fce038d225e75a283a Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 13:18:37 +0100 Subject: [PATCH 29/35] Adjust output verbosity. --- src/partitionEventsIndividual.cpp | 2 +- src/partitionEventsSummed.cpp | 2 +- src/processCommon.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/partitionEventsIndividual.cpp b/src/partitionEventsIndividual.cpp index 5a337ac..aca4e74 100644 --- a/src/partitionEventsIndividual.cpp +++ b/src/partitionEventsIndividual.cpp @@ -26,7 +26,7 @@ void partitionEventsIndividual(const std::vector &inputNeXusFiles, for (auto &nxsFileName : inputNeXusFiles) { // Open the NeXuS file and get its event data - NeXuSFile nxs(nxsFileName); + NeXuSFile nxs(nxsFileName, true); nxs.prepareSpectraSpace(); nxs.loadEventData(); diff --git a/src/partitionEventsSummed.cpp b/src/partitionEventsSummed.cpp index eb91a20..92e6f6d 100644 --- a/src/partitionEventsSummed.cpp +++ b/src/partitionEventsSummed.cpp @@ -27,7 +27,7 @@ void partitionEventsSummed(const std::vector &inputNeXusFiles, std: for (auto &nxsFileName : inputNeXusFiles) { // Open the NeXuS file and get its event data - NeXuSFile nxs(nxsFileName); + NeXuSFile nxs(nxsFileName, true); nxs.loadEventData(); const auto &eventsPerFrame = nxs.eventsPerFrame(); diff --git a/src/processCommon.cpp b/src/processCommon.cpp index 1fb7554..704ca2e 100644 --- a/src/processCommon.cpp +++ b/src/processCommon.cpp @@ -48,7 +48,7 @@ std::vector> prepareSlices(const Window &window, in auto &[newWin, nexus] = slices.emplace_back(Window(sliceName.str(), sliceStartTime, sliceDuration), NeXuSFile(outputFileName.str())); nexus.loadBasicData(); - nexus.prepareSpectraSpace(true); + nexus.prepareSpectraSpace(); sliceStartTime += sliceDuration; } From 78ecddb18ad666a530cf132155491b536098f13c Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 13:53:50 +0100 Subject: [PATCH 30/35] Don't forget to load monitor counts when preparing slices. --- src/processCommon.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/processCommon.cpp b/src/processCommon.cpp index 704ca2e..1aa7214 100644 --- a/src/processCommon.cpp +++ b/src/processCommon.cpp @@ -49,6 +49,7 @@ std::vector> prepareSlices(const Window &window, in slices.emplace_back(Window(sliceName.str(), sliceStartTime, sliceDuration), NeXuSFile(outputFileName.str())); nexus.loadBasicData(); nexus.prepareSpectraSpace(); + nexus.loadMonitorCounts(); sliceStartTime += sliceDuration; } From 5e7bdd511ae5f35219b996628b9da2c442b98ef7 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 16:05:48 +0100 Subject: [PATCH 31/35] Willingly break system packages on the Mac. --- .github/workflows/build/osx/action.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build/osx/action.yml b/.github/workflows/build/osx/action.yml index 47ba37a..444d30d 100644 --- a/.github/workflows/build/osx/action.yml +++ b/.github/workflows/build/osx/action.yml @@ -15,7 +15,7 @@ runs: - name: Install Python Dependencies shell: bash run: | - pip3 install --user conan==1.* + pip3 install --user conan==1.* --break-system-packages - name: Download HDF5 Artifacts shell: bash From 303c4db9eec70a932d4657f807d100a87d2c3fe4 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 16:12:39 +0100 Subject: [PATCH 32/35] Don't hardcode gsl path. --- .github/workflows/build/osx/action.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.github/workflows/build/osx/action.yml b/.github/workflows/build/osx/action.yml index 444d30d..70ff175 100644 --- a/.github/workflows/build/osx/action.yml +++ b/.github/workflows/build/osx/action.yml @@ -29,11 +29,10 @@ runs: set -ex export PATH="$(python3 -m site --user-base)/bin:$PATH" HDF5_DIR="$(pwd)/${{ env.hdf5tag }}" - GSL_DIR="/usr/local/Cellar/gsl/2.7.1" mkdir build && cd build conan install ../ - cmake ../ -G Ninja -DCMAKE_Fortran_COMPILER:string="gfortran-11" -DLOCAL_STATIC_HDF5:bool=true -DHDF5_DIR:path=${HDF5_DIR} -DLOCAL_STATIC_GSL:bool=true -DGSL_DIR:path=${GSL_DIR} + cmake ../ -G Ninja -DCMAKE_Fortran_COMPILER:string="gfortran-11" -DLOCAL_STATIC_HDF5:bool=true -DHDF5_DIR:path=${HDF5_DIR} ninja - name: Create Zip From a0658820ce86b5c65572a3b1bd1e3235e9dc215b Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 16:17:42 +0100 Subject: [PATCH 33/35] No need for a Fortran compiler. --- .github/workflows/build/osx/action.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/build/osx/action.yml b/.github/workflows/build/osx/action.yml index 70ff175..19d1b03 100644 --- a/.github/workflows/build/osx/action.yml +++ b/.github/workflows/build/osx/action.yml @@ -32,7 +32,7 @@ runs: mkdir build && cd build conan install ../ - cmake ../ -G Ninja -DCMAKE_Fortran_COMPILER:string="gfortran-11" -DLOCAL_STATIC_HDF5:bool=true -DHDF5_DIR:path=${HDF5_DIR} + cmake ../ -G Ninja -DLOCAL_STATIC_HDF5:bool=true -DHDF5_DIR:path=${HDF5_DIR} ninja - name: Create Zip From 7dad3847be1209a5ed40cdc85932c994e2c7aa69 Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 16:19:59 +0100 Subject: [PATCH 34/35] Add GSL include dirs to search path. --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 645517c..9808db5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -44,6 +44,7 @@ if(LOCAL_STATIC_GSL) else(LOCAL_STATIC_GSL) find_package(GSL REQUIRED) list(APPEND LINK_LIBS "${GSL_LIBRARIES}") + include_directories("${GSL_INCLUDE_DIRS}") endif(LOCAL_STATIC_GSL) option( From d4a648d03180895ef1f55bc202c30532b5c3c64f Mon Sep 17 00:00:00 2001 From: Tristan Youngs Date: Mon, 21 Oct 2024 16:30:34 +0100 Subject: [PATCH 35/35] Const? --- src/nexusFile.cpp | 2 +- src/nexusFile.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index f5d5051..e5c7221 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -37,7 +37,7 @@ NeXuSFile::NeXuSFile(std::string filename, bool printInfo) : filename_(filename) void NeXuSFile::operator=(NeXuSFile &source) { copy(source); } -NeXuSFile::NeXuSFile(NeXuSFile &source) { copy(source); } +NeXuSFile::NeXuSFile(const NeXuSFile &source) { copy(source); } NeXuSFile::NeXuSFile(NeXuSFile &&source) { diff --git a/src/nexusFile.h b/src/nexusFile.h index 005c19a..0052bb3 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -11,7 +11,7 @@ class NeXuSFile public: NeXuSFile(std::string filename = "", bool printInfo = false); void operator=(NeXuSFile &source); - NeXuSFile(NeXuSFile &source); + NeXuSFile(const NeXuSFile &source); NeXuSFile(NeXuSFile &&source); ~NeXuSFile(); // Clear all data and arrays