diff --git a/.github/workflows/build/osx/action.yml b/.github/workflows/build/osx/action.yml index 47ba37a..19d1b03 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 @@ -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 -DLOCAL_STATIC_HDF5:bool=true -DHDF5_DIR:path=${HDF5_DIR} ninja - name: Create Zip diff --git a/CMakeLists.txt b/CMakeLists.txt index 645517c..22c0bd3 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( @@ -140,8 +141,7 @@ if(BUILD_HDF) endif(BUILD_SZIP) else(BUILD_HDF) if(LOCAL_STATIC_HDF5) - list(APPEND LINK_LIBS "${HDF5_DIR}/lib/libhdf5_cpp.a" - "${HDF5_DIR}/lib/libhdf5.a" + list(APPEND LINK_LIBS "${HDF5_DIR}/lib/libhdf5.a" "${HDF5_DIR}/lib/libszip-static.a") if(WIN32) list(APPEND LINK_LIBS "${HDF5_DIR}/lib/libzlibstatic.a") diff --git a/flake.lock b/flake.lock index 8554147..1d5cd11 100644 --- a/flake.lock +++ b/flake.lock @@ -149,16 +149,16 @@ }, "nixpkgs_4": { "locked": { - "lastModified": 1688392541, - "narHash": "sha256-lHrKvEkCPTUO+7tPfjIcb7Trk6k31rz18vkyqmkeJfY=", + "lastModified": 1704290814, + "narHash": "sha256-LWvKHp7kGxk/GEtlrGYV68qIvPHkU9iToomNFGagixU=", "owner": "NixOS", "repo": "nixpkgs", - "rev": "ea4c80b39be4c09702b0cb3b42eab59e2ba4f24b", + "rev": "70bdadeb94ffc8806c0570eb5c2695ad29f0e421", "type": "github" }, "original": { "owner": "NixOS", - "ref": "nixos-22.11", + "ref": "nixos-23.05", "repo": "nixpkgs", "type": "github" } diff --git a/flake.nix b/flake.nix index abefff5..92707e1 100644 --- a/flake.nix +++ b/flake.nix @@ -1,8 +1,8 @@ { inputs = { - nixpkgs.url = "github:NixOS/nixpkgs/nixos-22.11"; - future.url = "github:NixOS/nixpkgs/nixos-unstable"; + nixpkgs.url = "github:NixOS/nixpkgs/nixos-23.05"; outdated.url = "github:NixOS/nixpkgs/nixos-21.05"; + future.url = "github:NixOS/nixpkgs/nixos-unstable"; }; outputs = { self, nixpkgs, future, outdated, flake-utils, bundlers }: @@ -34,13 +34,11 @@ name = "np-shell"; buildInputs = base_libs pkgs ++ (with pkgs; [ - (pkgs.clang-tools.override { - llvmPackages = pkgs.llvmPackages_13; - }) ccache ccls cmake-format cmake-language-server + clang-tools distcc gdb valgrind diff --git a/np.cpp b/np.cpp index f609454..a463b98 100644 --- a/np.cpp +++ b/np.cpp @@ -26,46 +26,142 @@ int main(int argc, char **argv) double windowDelta_{0.0}; // Number of slices to partition window in to int windowSlices_{1}; - // Target spectrum for event get (optional) - std::optional spectrumId_; + // Integer value for processing mode (if required) + int processingInt_; // Define and parse CLI arguments - CLI::App app("NeXuS Processor (np), Copyright (C) 2024 Jared Swift and Tristan Youngs."); + CLI::App app("NeXuS Processor (np), Copyright (C) 2024-2025 Jared Swift and Tristan Youngs.\n\nNotes:\n- Detector and " + "monitor spectrum indices start at 1.\n"); // -- 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)") + "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")->required(); - 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_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 rather than 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"); + 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 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 id) + { + if (processingMode_ != Processors::ProcessingMode::None) + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + processingMode_ = Processors::ProcessingMode::DumpEvents; + processingInt_ = id; + }, + "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; + processingInt_ = id; + }, + "Print all events for specified detector index") + ->group("Processing"); + app.add_option_function( + "--count-detector", + [&](int id) + { + if (processingMode_ != Processors::ProcessingMode::None) + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + processingMode_ = Processors::ProcessingMode::CountDetector; + processingInt_ = id; + }, + "Count events in specified detector histogram") + ->group("Processing"); + app.add_option_function( + "--count-monitor", + [&](int id) + { + if (processingMode_ != Processors::ProcessingMode::None) + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + processingMode_ = Processors::ProcessingMode::CountMonitor; + processingInt_ = id; + }, + "Count events in specified monitor histogram") + ->group("Processing"); + app.add_option_function( + "--dump-detector", + [&](int id) + { + if (processingMode_ != Processors::ProcessingMode::None) + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + processingMode_ = Processors::ProcessingMode::DumpDetector; + processingInt_ = id; + }, + "Dump specified detector histogram") + ->group("Processing"); + app.add_option_function( + "--dump-detector-from-events", + [&](int id) + { + if (processingMode_ != Processors::ProcessingMode::None) + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + processingMode_ = Processors::ProcessingMode::DumpDetectorFromEvents; + processingInt_ = id; + }, + "Dump specified detector histogram after constructing it from event data") + ->group("Processing"); + app.add_option_function( + "--dump-monitor", + [&](int id) + { + if (processingMode_ != Processors::ProcessingMode::None) + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + processingMode_ = Processors::ProcessingMode::DumpMonitor; + processingInt_ = id; + }, + "Dump specified monitor histogram") + ->group("Processing"); app.add_flag_callback( "--summed", [&]() { - if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::Summed; - 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"); @@ -73,17 +169,28 @@ int main(int argc, char **argv) "--individual", [&]() { - if (processingMode_ == Processors::ProcessingMode::None) - processingMode_ = Processors::ProcessingMode::Individual; - 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"); - app.add_option("-l,--slices", windowSlices_, "Number of slices to split window definition in to (default = 1, no slicing)") + app.add_option_function( + "--resize-detectors", + [&](int id) + { + if (processingMode_ != Processors::ProcessingMode::None) + { + fmt::print("Error: Multiple processing modes given.\n"); + throw(CLI::RuntimeError()); + } + processingMode_ = Processors::ProcessingMode::ResizeDetectors; + processingInt_ = id; + }, + "Dump all events for specified detector index") ->group("Processing"); // -- Post Processing app.add_flag_callback( @@ -94,56 +201,86 @@ int main(int argc, char **argv) "--scale-detectors", [&]() { Processors::postProcessingMode_ = Processors::PostProcessingMode::ScaleDetectors; }, "Scale detector counts in final output to match the number of frames used for monitor counts") ->group("Post-Processing"); + // -- Other Options + app.add_flag_callback( + "-v,--verbose", [&]() { NeXuSFile::verbose = true; }, "Enable verbose output") + ->group("Other Options"); 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"); + fmt::print("No processing mode specified. Basic data from files will be shown.\n"); + for (const auto &file : inputFiles_) + { + NeXuSFile nxs(file, true); + nxs.prepareSpectraSpace(true); + } + break; + case (Processors::ProcessingMode::DumpEvents): + case (Processors::ProcessingMode::PrintEvents): + Processors::dumpEventTimesEpoch(inputFiles_, processingInt_, + processingMode_ == Processors::ProcessingMode::PrintEvents); + break; + case (Processors::ProcessingMode::CountDetector): + Processors::countDetector(inputFiles_, processingInt_); + break; + case (Processors::ProcessingMode::CountMonitor): + Processors::countMonitor(inputFiles_, processingInt_); + break; + case (Processors::ProcessingMode::DumpDetector): + Processors::dumpDetector(inputFiles_, processingInt_); break; - case (Processors::ProcessingMode::Individual): - Processors::processIndividual(inputFiles_, outputDirectory_, window, windowSlices_, windowDelta_); + case (Processors::ProcessingMode::DumpDetectorFromEvents): + Processors::dumpDetectorFromEvents(inputFiles_, processingInt_); + break; + case (Processors::ProcessingMode::DumpMonitor): + Processors::dumpMonitor(inputFiles_, processingInt_); + break; + 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()); + 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::partitionEventsIndividual(inputFiles_, outputDirectory_, + {windowName_, windowStartTime_ + windowOffset_, windowWidth_}, + windowSlices_, windowDelta_); + else + Processors::partitionEventsSummed(inputFiles_, outputDirectory_, + {windowName_, windowStartTime_ + windowOffset_, windowWidth_}, windowSlices_, + windowDelta_); break; - case (Processors::ProcessingMode::Summed): - Processors::processSummed(inputFiles_, outputDirectory_, window, windowSlices_, windowDelta_); + case (Processors::ProcessingMode::ResizeDetectors): + Processors::resizeDetectors(inputFiles_, processingInt_); break; default: throw(std::runtime_error("Unhandled processing mode.\n")); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2b200fa..bfe565e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,16 +1,27 @@ -add_library(nexusProcess +add_library( + nexusProcess + countDetector.cpp + countMonitor.cpp + dumpDetector.cpp + dumpDetectorFromEvents.cpp + dumpMonitor.cpp + frameChunker.cpp getEvents.cpp + histogram.cpp nexusFile.cpp + partitionEventsIndividual.cpp + partitionEventsSummed.cpp processCommon.cpp - processIndividual.cpp - processSummed.cpp + resizeDetectors.cpp window.cpp + frameChunker.h + histogram.h nexusFile.h processors.h - window.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) diff --git a/src/countDetector.cpp b/src/countDetector.cpp new file mode 100644 index 0000000..1a09404 --- /dev/null +++ b/src/countDetector.cpp @@ -0,0 +1,29 @@ +#include "nexusFile.h" +#include "processors.h" +#include +#include + +namespace Processors +{ +void countDetector(const std::vector &inputNeXusFiles, int detectorIndex) +{ + /* + * Count events in the specified detector histogram + */ + + fmt::print("Counting events for detector index {}...\n", detectorIndex); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and load in detector counts + NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); + nxs.loadDetectorCounts(); + + const auto spectrumId = nxs.spectrumForDetector(detectorIndex); + fmt::print("{}\n", nxs.detectorHistograms().at(spectrumId).sum()); + } +} + +} // namespace Processors diff --git a/src/countMonitor.cpp b/src/countMonitor.cpp new file mode 100644 index 0000000..86d2647 --- /dev/null +++ b/src/countMonitor.cpp @@ -0,0 +1,30 @@ +#include "nexusFile.h" +#include "processors.h" +#include +#include + +namespace Processors +{ +void countMonitor(const std::vector &inputNeXusFiles, int monitorIndex) +{ + /* + * Count events in the specified monitor histogram + */ + + fmt::print("Counting events for monitor index {}...\n", monitorIndex); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and load in monitor counts + NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); + nxs.loadMonitorCounts(); + + const auto &counts = nxs.monitorCounts().at(monitorIndex); + auto sum = std::accumulate(counts.begin(), counts.end(), 0); + fmt::print("{}\n", sum); + } +} + +} // namespace Processors diff --git a/src/dumpDetector.cpp b/src/dumpDetector.cpp new file mode 100644 index 0000000..c4f11bb --- /dev/null +++ b/src/dumpDetector.cpp @@ -0,0 +1,41 @@ +#include "nexusFile.h" +#include "processors.h" +#include +#include +#include + +namespace Processors +{ +void dumpDetector(const std::vector &inputNeXusFiles, int detectorIndex) +{ + /* + * Dump histograms for the specified detector index + */ + + fmt::print("Dumping histogram for detector index {}...\n", detectorIndex); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and load in detector counts + NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); + nxs.loadDetectorCounts(); + + const auto spectrumId = nxs.spectrumForDetector(detectorIndex); + + // Open the output file + auto filename = std::string(std::filesystem::path(nxsFileName).filename().c_str()); + std::ofstream output(fmt::format("{}.det.{}", filename, detectorIndex).c_str()); + output << fmt::format("# TCB/us Counts [detector index {}, spectrum index = {}]\n", detectorIndex, spectrumId); + auto bin = 0; + const auto &histo = nxs.detectorHistograms().at(spectrumId); + for (auto tof : nxs.tofBoundaries()) + { + output << fmt::format("{} {}\n", tof, histo.value(bin++)); + } + output.close(); + } +} + +} // namespace Processors diff --git a/src/dumpDetectorFromEvents.cpp b/src/dumpDetectorFromEvents.cpp new file mode 100644 index 0000000..50144a3 --- /dev/null +++ b/src/dumpDetectorFromEvents.cpp @@ -0,0 +1,61 @@ +#include "frameChunker.h" +#include "nexusFile.h" +#include "processors.h" +#include +#include +#include + +namespace Processors +{ +void dumpDetectorFromEvents(const std::vector &inputNeXusFiles, int detectorIndex) +{ + /* + * Dump histograms for the specified detector index after constructing it from event data + */ + + fmt::print("Dumping histogram for detector index {}...\n", detectorIndex); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and prepare spectra space for histograms + NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); + + // Get spectrum index and histogram + const auto spectrumId = nxs.spectrumForDetector(detectorIndex); + auto &histo = nxs.detectorHistograms().at(spectrumId); + + // Read in chunked frames + FrameChunker frameChunker(nxs); + while (frameChunker.getNextFrameData()) + { + auto &frameData = frameChunker.frameData(); + for (const auto &frame : frameData) + { + const auto &eventIndices = frame.detectorIndices; + const auto &eventTimes = frame.times; + + // Loop over events in this frame + for (auto i = 0; i < eventIndices.size(); ++i) + { + if (eventIndices[i] == spectrumId) + histo.bin(eventTimes[i]); + } + } + } + + // Open the output file + auto filename = std::string(std::filesystem::path(nxsFileName).filename().c_str()); + std::ofstream output(fmt::format("{}.detev.{}", filename, detectorIndex).c_str()); + output << fmt::format("# TCB/us Counts [detector index {}, spectrum index = {}]\n", detectorIndex, spectrumId); + auto bin = 0; + for (auto tof : nxs.tofBoundaries()) + { + output << fmt::format("{} {}\n", tof, histo.value(bin++)); + } + output.close(); + } +} + +} // namespace Processors diff --git a/src/dumpMonitor.cpp b/src/dumpMonitor.cpp new file mode 100644 index 0000000..2584a5e --- /dev/null +++ b/src/dumpMonitor.cpp @@ -0,0 +1,40 @@ +#include "nexusFile.h" +#include "processors.h" +#include +#include +#include + +namespace Processors +{ +void dumpMonitor(const std::vector &inputNeXusFiles, int monitorIndex) +{ + /* + * Dump histograms for the specified monitor index + */ + + fmt::print("Dumping histogram for monitor index {}...\n", monitorIndex); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and load in monitor counts + NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); + nxs.loadMonitorCounts(); + + // Open the output file + auto filename = std::string(std::filesystem::path(nxsFileName).filename().c_str()); + std::ofstream output(fmt::format("{}.mon.{}", filename, monitorIndex).c_str()); + output << "# TCB/us Counts\n"; + + const auto &counts = nxs.monitorCounts().at(monitorIndex); + auto bin = 0; + for (auto tof : nxs.tofBoundaries()) + { + output << fmt::format("{} {}\n", tof, counts[bin++]); + } + output.close(); + } +} + +} // namespace Processors diff --git a/src/frameChunker.cpp b/src/frameChunker.cpp new file mode 100644 index 0000000..2c6887c --- /dev/null +++ b/src/frameChunker.cpp @@ -0,0 +1,107 @@ +#include "frameChunker.h" +#include +#include + +FrameChunker::FrameChunker(NeXuSFile &source) : neXuSFile_(source) +{ + // Initialise frame chunk storage + frameData_.reserve(frameChunkSize_); + nextFrameIndex_ = 0; + + // Open input NeXuS file in read only mode. + fileHandle_ = H5::H5File(neXuSFile_.filename(), H5F_ACC_RDONLY); + + // Get event counts per frame + auto &&[eventsPerFrameID, eventsPerFrameDimension] = + NeXuSFile::get1DDataset(fileHandle_, "raw_data_1/framelog/events_log", "value"); + eventsPerFrame_.resize(eventsPerFrameDimension); + H5Dread(eventsPerFrameID.getId(), H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, eventsPerFrame_.data()); + + // Get first event indices per frame + auto &&[frameFirstIndices, frameFirstIndicesDimension] = + NeXuSFile::get1DDataset(fileHandle_, "raw_data_1/detector_1_events", "event_index"); + frameFirstIndices_.resize(frameFirstIndicesDimension); + H5Dread(frameFirstIndices.getId(), H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, frameFirstIndices_.data()); + + // Get total number of events + auto &&[totalCountsID, totalCountsDimension] = + NeXuSFile::get1DDataset(fileHandle_, "raw_data_1/detector_1_events", "total_counts"); + std::array totalCountsBuffer; + H5Dread(totalCountsID.getId(), H5T_NATIVE_LONG, H5S_ALL, H5S_ALL, H5P_DEFAULT, totalCountsBuffer.data()); + totalEvents_ = totalCountsBuffer[0]; + + // Read in good frames + auto &&[goodFramesID, goodFramesDimension] = NeXuSFile::get1DDataset(fileHandle_, "raw_data_1", "good_frames"); + std::array goodFramesTemp; + H5Dread(goodFramesID.getId(), H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, goodFramesTemp.data()); + totalFrames_ = goodFramesTemp[0]; + + // Read in frame offsets. + auto &&[frameZeroID, frameZerosDimension] = + NeXuSFile::get1DDataset(fileHandle_, "raw_data_1/detector_1_events", "event_time_zero"); + frameZero_.resize(frameZerosDimension); + H5Dread(frameZeroID.getId(), H5T_NATIVE_LDOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, frameZero_.data()); + + fmt::print("There are {} total events.\n", totalEvents_); +} + +/* + * Current Frame Chunk + */ + +// Move to next frame data +bool FrameChunker::getNextFrameData() +{ + // If our next frame index is out of range, return false as we are done + if (nextFrameIndex_ >= totalFrames_) + return false; + + // Clear current data + frameData_.clear(); + + // Get dataset handles + auto &&[indicesDataset, indicesDimension] = + NeXuSFile::get1DDataset(fileHandle_, "raw_data_1/detector_1_events", "event_id"); + auto &&[timesDataset, timesDimension] = + NeXuSFile::get1DDataset(fileHandle_, "raw_data_1/detector_1_events", "event_time_offset"); + + // Read in events for the current frame + auto nextFrameLimit = std::min(nextFrameIndex_ + frameChunkSize_, totalFrames_); + for (auto i = nextFrameIndex_; i < nextFrameLimit; ++i) + { + // Push a new frame data + auto &frame = frameData_.emplace_back(); + frame.index = i; + frame.timeZero = frameZero_[i]; + + // Set vector sizes + frame.detectorIndices.resize(eventsPerFrame_[i]); + frame.times.resize(eventsPerFrame_[i]); + + // Set up dataspaces to read in only events for this frame. Note that the order of parameters given to selectHyperslab() + // is different in the actual API - online manual lists [start, stride, count, block] but it is actually [count, start, + // stride, block]. + std::array start = {(hsize_t)frameFirstIndices_[i]}, stride = {1}, block = {(hsize_t)eventsPerFrame_[i]}, + count = {1}; + H5::DataSpace indicesSpace = indicesDataset.getSpace(); + indicesSpace.selectHyperslab(H5S_SELECT_SET, count.data(), start.data(), stride.data(), block.data()); + H5::DataSpace timesSpace = timesDataset.getSpace(); + timesSpace.selectHyperslab(H5S_SELECT_SET, count.data(), start.data(), stride.data(), block.data()); + + // Also need a memory dataspace of the same size + auto memSpace = H5Screate_simple(1, block.data(), block.data()); + + // Read the partial data + H5Dread(indicesDataset.getId(), H5T_STD_I32LE, memSpace, indicesSpace.getId(), H5P_DEFAULT, + frame.detectorIndices.data()); + H5Dread(timesDataset.getId(), H5T_NATIVE_DOUBLE, memSpace, timesSpace.getId(), H5P_DEFAULT, frame.times.data()); + } + + // Set frame index for next call + nextFrameIndex_ = nextFrameLimit; + + return true; +} + +// Return current frame data +std::vector &FrameChunker::frameData() { return frameData_; } \ No newline at end of file diff --git a/src/frameChunker.h b/src/frameChunker.h new file mode 100644 index 0000000..a54269c --- /dev/null +++ b/src/frameChunker.h @@ -0,0 +1,59 @@ +#pragma once + +#include "nexusFile.h" +#include + +struct FrameData +{ + // Index of this frame + int index; + // Time zero of the frame, in seconds + double timeZero; + // Event detector indices (detector_1_events/event_id) - (M+1 -> N+M+1) + std::vector detectorIndices; + // Event times (detector_1_events/event_time_offset) in us relative to frame start + std::vector times; +}; + +class FrameChunker +{ + public: + FrameChunker(NeXuSFile &source); + ~FrameChunker() = default; + + /* + * File and Basic Info + */ + private: + // Source NeXuS file + NeXuSFile &neXuSFile_; + // HDF5 file handle + H5::H5File fileHandle_; + // Total number of events + long long totalEvents_; + // First event index of each frame (detector_1_events/event_index) + std::vector frameFirstIndices_; + // Number of events per frame (framelog/events_log/value) + std::vector eventsPerFrame_; + // Frame start times (detector_1_events/event_time_zero) - seconds, relative to start time since epoch + std::vector frameZero_; + // Total number of frames + int totalFrames_{0}; + // Frame chunk size + const int frameChunkSize_{1024}; + + /* + * Current Frame Chunk + */ + private: + // Next frame chunk to distil into frame data + int nextFrameIndex_{0}; + // Current chunk of frames + std::vector frameData_; + + public: + // Move to next frame data + bool getNextFrameData(); + // Return current frame data + std::vector &frameData(); +}; diff --git a/src/getEvents.cpp b/src/getEvents.cpp index 0d99163..505104b 100644 --- a/src/getEvents.cpp +++ b/src/getEvents.cpp @@ -1,71 +1,84 @@ +#include "frameChunker.h" #include "nexusFile.h" #include "processors.h" #include "window.h" +#include #include +#include +#include #include namespace Processors { -// Get events from specified spectrum, returning seconds since epoch for each -std::map> getEvents(const std::vector &inputNeXusFiles, int spectrumId, 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 */ - printf("Get events...\n"); - fmt::print("Target detector spectrum is {}\n", spectrumId); + char timeBuffer[20]; - std::map> eventMap; - std::optional lastSecondsSinceEpoch; + fmt::print("Retrieving all events from detector index {}...\n", detectorIndex); - // Loop over input Nexus files + // 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(); - nxs.loadTimes(); - fmt::print("... file '{}' has {} events...\n", nxsFileName, nxs.eventTimes().size()); + nxs.prepareSpectraSpace(); - 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(); + std::optional lastSecondsSinceEpoch; - // 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]; + const auto spectrumId = nxs.spectrumForDetector(detectorIndex); + fmt::print("NeXuS file spectrum ID for detector index {} is {}.\n", detectorIndex, spectrumId); + + std::ofstream fileOutput; + if (!toStdOut) + fileOutput.open(fmt::format("{}.events.{}", nxsFileName, detectorIndex).c_str()); + std::ostream &output = toStdOut ? std::cout : fileOutput; + + // Prepare the frame chunker + FrameChunker frameChunker(nxs); - for (auto k = eventStart; k < eventEnd; ++k) + // Write header + output << fmt::format("# {:20s} {:20s} {:20s} {:20s} {}\n", "frame_offset(us)", "start_time_offset(s)", + "epoch_offset(s)", "local time", "delta(s)"); + + // Read in chunked frames + while (frameChunker.getNextFrameData()) + { + auto &frameData = frameChunker.frameData(); + for (const auto &frame : frameData) { - if (eventIndices[k] == spectrumId) + const auto &eventIndices = frame.detectorIndices; + const auto &eventTimes = frame.times; + + // Loop over events in this frame + for (auto i = 0; i < eventIndices.size(); ++i) { - auto eMicroSeconds = eventTimes[k]; - auto eSeconds = eMicroSeconds * 0.000001; - auto eSecondsSinceEpoch = eSeconds + frameZero + nxs.startSinceEpoch(); + if (eventIndices[i] != spectrumId) + continue; + + auto eSeconds = eventTimes[i] * 0.000001; + auto eSecondsSinceEpoch = eSeconds + frame.timeZero + nxs.startSinceEpoch(); + auto convertedSeconds = time_t(eSecondsSinceEpoch); + strftime(timeBuffer, 20, "%d/%m/%y %H:%M:%S", std::localtime(&convertedSeconds)); if (lastSecondsSinceEpoch) - fmt::print("{:20.6f} {:20.10f} {:20.5f} {}\n", eMicroSeconds, eSeconds + frameZero, - eSecondsSinceEpoch, eSecondsSinceEpoch - *lastSecondsSinceEpoch); + output << fmt::format("{:20.6f} {:20.10f} {:20.5f} {:20s} {}\n", eventTimes[i], + eSeconds + frame.timeZero, eSecondsSinceEpoch, timeBuffer, + eSecondsSinceEpoch - *lastSecondsSinceEpoch); else - fmt::print("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch); - eventMap[spectrumId].push_back(eSecondsSinceEpoch); + output << fmt::format("{:20.6f} {:20.10f} {:20.5f} {:20s}\n", eventTimes[i], + eSeconds + frame.timeZero, eSecondsSinceEpoch, timeBuffer); + lastSecondsSinceEpoch = eSecondsSinceEpoch; - if (firstOnly) - return eventMap; } } - - // Update start event index - eventStart = eventEnd; } - } - return eventMap; + if (!toStdOut) + fileOutput.close(); + } } } // namespace Processors diff --git a/src/histogram.cpp b/src/histogram.cpp new file mode 100644 index 0000000..0c9ce68 --- /dev/null +++ b/src/histogram.cpp @@ -0,0 +1,42 @@ +#include "histogram.h" +#include +#include + +// Initialise +void IntegerHistogram::initialise(const std::vector &boundaries) +{ + boundaries_ = boundaries; + bins_.resize(boundaries_.size() - 1); + std::fill(bins_.begin(), bins_.end(), 0); +} + +// Increment relevant bin +void IntegerHistogram::bin(double x) +{ + // Out of initial range? + if (x < boundaries_.front()) + return; + + // Find relevant bin + for (auto i = 1; i < boundaries_.size(); ++i) + if (x < boundaries_[i]) + { + ++bins_[i - 1]; + return; + } +} + +// Add to specific bin +void IntegerHistogram::add(int index, long int value) { bins_[index] += value; } + +// Return specified bin value +long int IntegerHistogram::value(int index) const { return bins_[index]; } + +// Return the sum of all bins +long int IntegerHistogram::sum() const { return std::accumulate(bins_.begin(), bins_.end(), (long int)0); } + +// Scale bin values +void IntegerHistogram::scale(double factor) +{ + std::transform(bins_.begin(), bins_.end(), bins_.begin(), [factor](const long int value) { return value * factor; }); +} diff --git a/src/histogram.h b/src/histogram.h new file mode 100644 index 0000000..68542d2 --- /dev/null +++ b/src/histogram.h @@ -0,0 +1,30 @@ +#pragma once + +#include + +class IntegerHistogram +{ + public: + IntegerHistogram() = default; + ~IntegerHistogram() = default; + + private: + // Bins + std::vector bins_; + // Boundaries + std::vector boundaries_; + + public: + // Initialise + void initialise(const std::vector &boundaries); + // Increment relevant bin + void bin(double x); + // Add to specific bin + void add(int index, long int value); + // Return specified bin value + long int value(int index) const; + // Return the sum of all bins + long int sum() const; + // Scale bin values + void scale(double factor); +}; diff --git a/src/modex.cpp b/src/modex.cpp deleted file mode 100644 index a767bee..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 diff --git a/src/nexusFile.cpp b/src/nexusFile.cpp index fe170d8..f765b22 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -5,14 +5,17 @@ #include #include +// Whether verbose output is enabled +bool NeXuSFile::verbose = false; + // Basic paths required when copying / creating a NeXuS file 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,59 +24,219 @@ 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, bool loadEvents) : filename_(filename) +NeXuSFile::NeXuSFile(std::string filename, bool printInfo) : filename_(filename) { - if (loadEvents) + if (!filename_.empty()) { - fmt::print("Loading event data from file '{}'...\n", filename_); - loadFrameCounts(); - loadEventData(); - loadTimes(); - fmt::print("... file '{}' has {} goodframes and {} events...\n", filename_, nGoodFrames_, eventTimes_.size()); + fmt::print("Opening NeXuS file '{}'...\n", filename_); + + loadBasicData(printInfo); } } +void NeXuSFile::operator=(NeXuSFile &source) { copy(source); } + +NeXuSFile::NeXuSFile(const 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_; + nGoodFrames_ = source.nGoodFrames_; + startSinceEpoch_ = source.startSinceEpoch_; + endSinceEpoch_ = source.endSinceEpoch_; + tofBoundaries_ = source.tofBoundaries_; + monitorCounts_ = source.monitorCounts_; + detectorHistograms_ = source.detectorHistograms_; +} + +// Clear all data and arrays +void NeXuSFile::clear() +{ + filename_ = {}; + detectorSpectrumIndices_.clear(); + nMonitorSpectra_ = 0; + nMonitorFrames_ = 0; + nGoodFrames_ = 0; + startSinceEpoch_ = 0; + endSinceEpoch_ = 0; + tofBoundaries_.clear(); + monitorCounts_.clear(); + detectorHistograms_.clear(); +} + /* * I/O */ // Return handle and (simple) dimension for named leaf dataset -std::pair NeXuSFile::find1DDataset(H5::H5File file, H5std_string groupName, H5std_string datasetName) +std::pair NeXuSFile::get1DDataset(H5::H5File file, H5std_string groupName, H5std_string datasetName) { if (!file.nameExists(groupName)) + { + fmt::print("ERROR: No group named '{}' exists.\n", groupName); return {}; + } H5::Group group = file.openGroup(groupName); if (!group.nameExists(datasetName)) + { + fmt::print("ERROR: No dataset named '{}' exists in the group '{}' exists.\n", datasetName, groupName); return {}; + } H5::DataSet dataset = group.openDataSet(datasetName); H5::DataSpace space = dataset.getSpace(); - hsize_t spaceNDims = space.getSimpleExtentNdims(); + hsize_t spaceNDimensions = space.getSimpleExtentNdims(); + + std::vector spaceDimensions(spaceNDimensions); + space.getSimpleExtentDims(spaceDimensions.data()); + hsize_t nPoints = space.getSimpleExtentNpoints(); - hsize_t *spaceDims = new hsize_t[spaceNDims]; - space.getSimpleExtentDims(spaceDims); + // Create a string of the dimensions + std::string dims; + for (auto d : spaceDimensions) + dims += dims.empty() ? fmt::format("{}", d) : fmt::format(",{}", d); - return {dataset, spaceDims[0]}; + if (verbose) + fmt::print("Got dataset '{}' in group '{}', rank/dims/N = {}/{}/{}\n", datasetName, groupName, spaceNDimensions, dims, + nPoints); + return {dataset, spaceDimensions[0]}; +} + +// Resize 1D dataset +void NeXuSFile::resize1DDataset(H5::DataSet dataset, std::vector dimensions) +{ + H5::DataSpace space = dataset.getSpace(); + hsize_t spaceNDimensions = space.getSimpleExtentNdims(); + std::vector spaceDimensions(spaceNDimensions); + space.getSimpleExtentDims(spaceDimensions.data()); + + // Create strings of the dimensions + std::string oldDims, newDims; + for (auto d : spaceDimensions) + oldDims += oldDims.empty() ? fmt::format("{}", d) : fmt::format(",{}", d); + for (auto d : dimensions) + newDims += newDims.empty() ? fmt::format("{}", d) : fmt::format(",{}", d); + + hsize_t existingNPoints = space.getSimpleExtentNpoints(); + if (verbose) + fmt::print("Checking dataset size - existing rank/dims = {}/{}, new = {}/{}...\n", spaceNDimensions, oldDims, + dimensions.size(), newDims); + if (spaceNDimensions == existingNPoints && spaceDimensions == dimensions) + { + if (verbose) + fmt::print(" ... Extent has not changed - no change will be made.\n"); + } + + else + { + if (verbose) + fmt::print(" ... Extent has changed - resizing dataset...\n"); + // space.setExtentSimple(dimensions.size(), dimensions.data(), dimensions.data()); + H5Dset_extent(dataset.getId(), dimensions.data()); + } } // Return filename std::string NeXuSFile::filename() const { return filename_; } -// Template basic paths from the referenceFile, and make ready for histogram binning -void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) +// Load basic information from the NeXuS file +void NeXuSFile::loadBasicData(bool printInfo) { - filename_ = outputFile; + 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(referenceFile, H5F_ACC_RDONLY); + // Open input NeXuS file in read only mode. + H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); - // Create new Nexus file for output. - H5::H5File output = H5::H5File(filename_, H5F_ACC_TRUNC); + // Get the run title + auto &&[titleID, titleDimension] = NeXuSFile::get1DDataset(input, "raw_data_1", "title"); + H5Dread(titleID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, 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; + + auto &&[startTimeID, startTimeDimension] = NeXuSFile::get1DDataset(input, "raw_data_1", "start_time"); + H5Dread(startTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, charBuffer); + + 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; + stime.tm_mday = d; + stime.tm_hour = h; + stime.tm_min = m; + stime.tm_sec = s; + startSinceEpoch_ = (int)mktime(&stime); + if (printInfo) + fmt::print("... run started at {} ({} s since epoch).\n", charBuffer, startSinceEpoch_); + + // Read in end time in Unix time. + auto &&[endTimeID, endTimeDimension] = NeXuSFile::get1DDataset(input, "raw_data_1", "end_time"); + H5Dread(endTimeID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, charBuffer); + + 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; + etime.tm_mday = d; + etime.tm_hour = h; + etime.tm_min = m; + etime.tm_sec = s; + endSinceEpoch_ = (int)mktime(&etime); + if (printInfo) + { + 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 + auto &&[goodFramesID, goodFramesDimension] = NeXuSFile::get1DDataset(input, "raw_data_1", "good_frames"); + std::array goodFramesTemp; + H5Dread(goodFramesID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, goodFramesTemp.data()); + nGoodFrames_ = goodFramesTemp[0]; + if (printInfo) + fmt::print("... there were {} good frames.\n", nGoodFrames_); + + input.close(); +} + +// Template a new NeXusFile from that specified +void NeXuSFile::templateTo(std::string sourceFilename, std::string newFilename) +{ + // 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(newFilename, H5F_ACC_TRUNC); + + if (verbose) + 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) @@ -93,141 +256,92 @@ 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()); + input.close(); + 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 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()); + // Read in detector spectra information + auto &&[detSpecIndices, spectraDimension] = NeXuSFile::get1DDataset(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; - // Set up detector histograms - for (auto spec : spectra_) + if (printInfo) { - detectorHistograms_[spec] = gsl_histogram_alloc(tofBins_.size() - 1); - gsl_histogram_set_ranges(detectorHistograms_[spec], tofBins_.data(), tofBins_.size()); + fmt::print("... total number of detector spectra is {}.\n", detectorSpectrumIndices_.size()); + fmt::print("... inferred number of monitor spectra is {}.\n", nMonitorSpectra_); } - // 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) - { - auto &&[monitorSpectrum, monitorSpectrumDimension] = - NeXuSFile::find1DDataset(input, "/raw_data_1/monitor_" + std::to_string(i), "data"); - if (monitorSpectrum.getId() <= 0) - break; - - monitorCounts_[i].resize(tofBinsDimension); - H5Dread(monitorSpectrum.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, monitorCounts_[i].data()); + // Read in TOF boundary information. + auto &&[tofBoundariesID, tofBoundariesDimension] = + NeXuSFile::get1DDataset(input, "raw_data_1/detector_1", "time_of_flight"); + tofBoundaries_.resize(tofBoundariesDimension); + H5Dread(tofBoundariesID.getId(), H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBoundaries_.data()); - ++i; - } - - // Read in 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]; + // Set up detector histograms and straight counts vectors + for (auto spec : detectorSpectrumIndices_) + appendEmptyDetector(spec); input.close(); - output.close(); } -// Load frame counts -void NeXuSFile::loadFrameCounts() +// Load in monitor histograms +void NeXuSFile::loadMonitorCounts() { - printf("Load frame counts...\n"); + printf("Load monitor 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 - 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"); + // Read in monitor data - start from index 1 and end when we fail to find the named dataset with this suffix + for (auto i = 1; i <= nMonitorSpectra_; ++i) + { + auto &&[monitorSpectrum, monitorSpectrumDimension] = + NeXuSFile::get1DDataset(input, "/raw_data_1/monitor_" + std::to_string(i), "data"); - // Open our Nexus file in read only mode. - H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); + monitorCounts_[i].resize(tofBoundaries_.size()); + H5Dread(monitorSpectrum.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, monitorCounts_[i].data()); + } - // Read in event indices. - auto &&[eventIndicesID, eventIndicesDimension] = - NeXuSFile::find1DDataset(input, "raw_data_1/detector_1_events", "event_id"); - eventIndices_.resize(eventIndicesDimension); - H5Dread(eventIndicesID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, eventIndices_.data()); - - // Read in events. - auto &&[eventTimesID, eventTimesDimension] = - NeXuSFile::find1DDataset(input, "raw_data_1/detector_1_events", "event_time_offset"); - eventTimes_.resize(eventTimesDimension); - H5Dread(eventTimesID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, eventTimes_.data()); - - // Read in event counts per frame - auto &&[eventsPerFrameID, eventsPerFrameDimension] = - NeXuSFile::find1DDataset(input, "raw_data_1/framelog/events_log", "value"); - eventsPerFrame_.resize(eventsPerFrameDimension); - H5Dread(eventsPerFrameID.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, eventsPerFrame_.data()); - - // Read in frame offsets. - auto &&[frameOffsetsID, frameOffsetsDimension] = - NeXuSFile::find1DDataset(input, "raw_data_1/detector_1_events", "event_time_zero"); - frameOffsets_.resize(frameOffsetsDimension); - H5Dread(frameOffsetsID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, frameOffsets_.data()); + // 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::get1DDataset(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(); } -// Load start/end times -void NeXuSFile::loadTimes() +// Load detector counts from the file +void NeXuSFile::loadDetectorCounts() { - printf("Load times....\n"); + 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); - // 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); + const auto nSpec = detectorSpectrumIndices_.size(); + const auto nTOFBins = tofBoundaries_.size() - 1; - 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); + // Note to Future Me: While you might think that the buffer into which we will read must match the datatype of that stored + // (a 32-bit int) using a std::vector completely breaks the read and gives nothing but zeroes. + std::vector countsBuffer; + countsBuffer.resize(nSpec * nTOFBins); // Need contiguous memory. This is a pain. + auto &&[counts, detectorCountsDimension] = NeXuSFile::get1DDataset(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 &histo = detectorHistograms_[detectorSpectrumIndices_[i]]; + for (auto j = 0; j < nTOFBins; ++j) + histo.add(j, countsBuffer[i * nTOFBins + j]); + } input.close(); } @@ -235,35 +349,51 @@ void NeXuSFile::loadTimes() // 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); + const auto nSpec = detectorSpectrumIndices_.size(); + const auto nTOFBins = tofBoundaries_.size() - 1; + // Write good frames - std::array framesBuffer{0}; - framesBuffer[0] = nDetectorFrames_; - auto &&[goodFrames, goodFramesDimension] = NeXuSFile::find1DDataset(output, "raw_data_1", "good_frames"); + std::array framesBuffer; + framesBuffer[0] = nGoodFrames_; + auto &&[goodFrames, goodFramesDimension] = NeXuSFile::get1DDataset(output, "raw_data_1", "good_frames"); goodFrames.write(framesBuffer.data(), H5::PredType::STD_I32LE); // Write monitors for (auto &&[index, counts] : monitorCounts_) { auto &&[monitorCounts, monitorCountsDimension] = - NeXuSFile::find1DDataset(output, "raw_data_1/monitor_" + std::to_string(index), "data"); + NeXuSFile::get1DDataset(output, "raw_data_1/monitor_" + std::to_string(index), "data"); monitorCounts.write(counts.data(), H5::PredType::STD_I32LE); } - // Write detector counts - const auto nSpec = spectra_.size(); - const auto nTofBins = tofBins_.size() - 1; + // Write spectrum indices - we have to pass in any adjusted DataSpace otherwise we end up with the same sized array as + // before. + auto &&[detSpecIndices, spectraDimension] = NeXuSFile::get1DDataset(output, "raw_data_1/detector_1", "spectrum_index"); + resize1DDataset(detSpecIndices, {detectorSpectrumIndices_.size()}); + detSpecIndices.write(detectorSpectrumIndices_.data(), H5::PredType::STD_I32LE, detSpecIndices.getSpace(), + detSpecIndices.getSpace()); - auto *countsBuffer = new int[nSpec * nTofBins]; // HDF5 expects contiguous memory. This is a pain. + // Write detector counts + // Note to Future Me: Similar to the note above, using a long int here for the countsBuffer corrupts the resulting + // data in the file. Presumably this has something to so with the DataSpace type / extent in memory vs how it is + // to be stored in the HDF5 file, but honestly who cares at this point. + std::vector countsBuffer; + countsBuffer.resize(nSpec * nTOFBins); // Need 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); - auto &&[counts, detectorCountsDimension] = NeXuSFile::find1DDataset(output, "raw_data_1/detector_1", "counts"); - counts.write(countsBuffer, H5::PredType::STD_I32LE); + { + auto specID = detectorSpectrumIndices_[i]; + for (auto j = 0; j < nTOFBins; ++j) + countsBuffer[i * nTOFBins + j] = detectorHistograms_[specID].value(j); + } - delete[](countsBuffer); + // Get and resize the counts dataset - we have to pass in any adjusted DataSpace otherwise we end up with the same + // sized array as before. + auto &&[counts, detectorCountsDimension] = NeXuSFile::get1DDataset(output, "raw_data_1/detector_1", "counts"); + resize1DDataset(counts, {1, nSpec, nTOFBins}); + counts.write(countsBuffer.data(), H5::PredType::STD_I32LE, counts.getSpace(), counts.getSpace()); output.close(); @@ -275,33 +405,55 @@ bool NeXuSFile::saveModifiedData() */ int NeXuSFile::nGoodFrames() const { return nGoodFrames_; } +void NeXuSFile::zeroGoodFrames() { nGoodFrames_ = 0; } +void NeXuSFile::incrementGoodFrames(int delta) { nGoodFrames_ += delta; } int NeXuSFile::nMonitorFrames() const { return nMonitorFrames_; } -int NeXuSFile::nDetectorFrames() const { return nDetectorFrames_; } -void NeXuSFile::incrementDetectorFrameCount(int delta) { nDetectorFrames_ += delta; } int NeXuSFile::startSinceEpoch() const { return startSinceEpoch_; } int NeXuSFile::endSinceEpoch() const { return endSinceEpoch_; } -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::map> &NeXuSFile::monitorCounts() const { return monitorCounts_; } -std::map &NeXuSFile::detectorHistograms() { return detectorHistograms_; } -const std::map> &NeXuSFile::partitions() const { return partitions_; } +const std::vector &NeXuSFile::tofBoundaries() const { return tofBoundaries_; } +const int NeXuSFile::spectrumForDetector(int detectorId) const { return detectorSpectrumIndices_.at(detectorId - 1); } +const int NeXuSFile::nDetectors() const { return detectorSpectrumIndices_.size(); } +const std::map> &NeXuSFile::monitorCounts() const { return monitorCounts_; } +std::map &NeXuSFile::detectorHistograms() { return detectorHistograms_; } /* * Manipulation */ +// Remove the last detector spectrum, returning the index that was removed +int NeXuSFile::removeLastDetector() +{ + // Get target spectrum index + auto specID = detectorSpectrumIndices_.back(); + detectorHistograms_.erase(specID); + detectorSpectrumIndices_.pop_back(); + + return specID; +} + +// Append an empty detector +int NeXuSFile::appendEmptyDetector(int specID) +{ + // Generate and push target spectrum index if we need to + if (specID == -1) + { + specID = detectorSpectrumIndices_.empty() ? 0 : (detectorSpectrumIndices_.back() + 1); + detectorSpectrumIndices_.push_back(specID); + } + + // Initialise the histogram + detectorHistograms_[specID].initialise(tofBoundaries_); + + return specID; +} + // Scale monitors by specified factor void NeXuSFile::scaleMonitors(double factor) { for (auto &&[index, counts] : monitorCounts_) { for (auto &bin : counts) - { bin *= factor; - } } nMonitorFrames_ *= factor; @@ -311,15 +463,15 @@ void NeXuSFile::scaleMonitors(double factor) // Scale detectors by specified factor void NeXuSFile::scaleDetectors(double factor) { - auto oldSum = 0, newSum = 0; - for (auto i : spectra_) + unsigned long long oldSum = 0, newSum = 0; + for (auto i : detectorSpectrumIndices_) { - oldSum += gsl_histogram_sum(detectorHistograms_[i]); - gsl_histogram_scale(detectorHistograms_[i], factor); - newSum += gsl_histogram_sum(detectorHistograms_[i]); + oldSum += detectorHistograms_[i].sum(); + detectorHistograms_[i].scale(factor); + newSum += detectorHistograms_[i].sum(); } fmt::print(" ... Old counts was {}, now scaled to {} (ratio = {}).\n", oldSum, newSum, double(oldSum) / double(newSum)); - nDetectorFrames_ *= factor; - fmt::print(" ... New number of effective contributing detector frames is {}.\n", nDetectorFrames_); + nGoodFrames_ *= factor; + fmt::print(" ... New number of effective contributing detector frames is {}.\n", nGoodFrames_); } diff --git a/src/nexusFile.h b/src/nexusFile.h index be4922b..1452339 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -1,16 +1,24 @@ #pragma once +#include "histogram.h" #include -#include #include +#include #include #include class NeXuSFile { public: - NeXuSFile(std::string filename = "", bool loadEvents = false); - ~NeXuSFile() = default; + NeXuSFile(std::string filename = "", bool printInfo = false); + void operator=(NeXuSFile &source); + NeXuSFile(const 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 @@ -19,21 +27,29 @@ class NeXuSFile // Filename std::string filename_; - private: - // Return handle and (simple) dimension for named leaf dataset - static std::pair find1DDataset(H5::H5File file, H5std_string terminal, H5std_string datasetName); + public: + // Return handle and (simple) dimension for named leaf dataset, resizing if requested + static std::pair get1DDataset(H5::H5File file, H5std_string terminal, H5std_string datasetName); + // Resize 1D dataset + static void resize1DDataset(H5::DataSet dataset, std::vector dimensions); + + public: + // Whether verbose output is enabled + static bool verbose; public: // Return filename std::string filename() const; - // Template basic paths from the referenceFile, and make ready for histogram binning - void templateFile(std::string referenceFile, std::string outputFile); - // Load frame counts - void loadFrameCounts(); - // Load event data - void loadEventData(); - // Load start/end times - void loadTimes(); + // Load basic information from the NeXuS file + 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 detector counts from the file + void loadDetectorCounts(); // Save key modified data back to the file bool saveModifiedData(); @@ -41,41 +57,55 @@ class NeXuSFile * Data */ private: - std::vector spectra_; + /* + * Detector spectrum indices (detector_1/spectrum_index) + * These are detector, not monitor, indices, with indices running from M+1 -> N+M+1 where M is the + * number of monitors and N is the number of detectors. If there were no detectors the range would + * therefore be 1 -> N, for NIMROD's 9 monitors the first detector is 10, etc. The ordering of the + * indices (which is normally continuous?) reflects the order of data in the detector_1/counts + * array, while the indices themselves are used as keys for per-detector lookup in the + * detectorHistograms_ map. + */ + std::vector detectorSpectrumIndices_; + // Number of monitor spectra, determined from the first index in detectorSpectrumIndices_ + int nMonitorSpectra_{0}; + // Number of frames contributing to monitor counts int nMonitorFrames_{0}; - int nDetectorFrames_{0}; + // Number of frames contributing to detector counts int nGoodFrames_{0}; + // Start time in seconds since epoch (raw_data_1/start_time) int startSinceEpoch_{0}; + // End time in seconds since epoch (raw_data_1/end_time) int endSinceEpoch_{0}; - std::vector eventIndices_; - std::vector eventTimes_; - std::vector eventsPerFrame_; - std::vector frameOffsets_; - std::vector tofBins_; - std::map> monitorCounts_; - std::map detectorHistograms_; - std::map> partitions_; + // Time-of-flight bin boundaries in detector histograms (detector_1/time_of_flight) - us + std::vector tofBoundaries_; + // Monitor counts (TOF bins), mapped by monitor index (monitor_M/data) + std::map> monitorCounts_; + // Detector counts (TOF bins), mapped by detector index (M+1 -> N+M+1) (detector_1/counts) + std::map detectorHistograms_; public: [[nodiscard]] int nGoodFrames() const; + void zeroGoodFrames(); + void incrementGoodFrames(int delta = 1); [[nodiscard]] int nMonitorFrames() const; - [[nodiscard]] int nDetectorFrames() const; - void incrementDetectorFrameCount(int delta = 1); [[nodiscard]] int startSinceEpoch() const; [[nodiscard]] int endSinceEpoch() const; - [[nodiscard]] const std::vector &eventIndices() const; - [[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::map> &monitorCounts() const; - std::map &detectorHistograms(); - [[nodiscard]] const std::map> &partitions() const; + [[nodiscard]] const std::vector &tofBoundaries() const; + [[nodiscard]] const int spectrumForDetector(int detectorId) const; + [[nodiscard]] const int nDetectors() const; + [[nodiscard]] const std::map> &monitorCounts() const; + [[nodiscard]] const std::map> &detectorCounts() const; + std::map &detectorHistograms(); /* * Manipulation */ public: + // Remove the last detector spectrum, returning the index that was removed + int removeLastDetector(); + // Append an empty detector + int appendEmptyDetector(int specID = -1); // Scale monitors by specified factor void scaleMonitors(double factor); // Scale detectors by specified factor diff --git a/src/partitionEventsIndividual.cpp b/src/partitionEventsIndividual.cpp new file mode 100644 index 0000000..7e77751 --- /dev/null +++ b/src/partitionEventsIndividual.cpp @@ -0,0 +1,113 @@ +#include "frameChunker.h" +#include "nexusFile.h" +#include "processors.h" +#include "window.h" +#include +#include + +namespace Processors +{ +// 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("Partitioning events into individual windows/slices...\n"); + + // Copy the master window definition since we will be modifying it + auto window = windowDefinition; + std::vector> slices; + auto sliceIt = slices.end(); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and get its event data + NeXuSFile nxs(nxsFileName, true); + nxs.prepareSpectraSpace(); + + // Read in chunked frames + FrameChunker frameChunker(nxs); + auto tenPercentsDone = -1; + while (frameChunker.getNextFrameData()) + { + auto &frameData = frameChunker.frameData(); + auto currentTenPercents = int(100 * (frameData.front().index / double(nxs.nGoodFrames()))) / 10; + if (currentTenPercents > tenPercentsDone) + { + tenPercentsDone = currentTenPercents; + fmt::print("Processing frame / event data... {}%\n", tenPercentsDone * 10); + } + + for (const auto &frame : frameData) + { + const auto &eventIndices = frame.detectorIndices; + const auto &eventTimes = frame.times; + auto frameZero = frame.timeZero + nxs.startSinceEpoch(); + + // If the current slice end time is less than the frame zero, iterate the slices. + while (sliceIt != slices.end() && sliceIt->first.endTime() < frameZero) + { + sliceIt++; + + // If we have run out of slices propagate the set forward. + if (sliceIt == slices.end()) + { + // Save current data + postProcess(slices); + saveSlices(slices); + + // Empty current slices + slices.clear(); + sliceIt = slices.end(); + } + } + + // Do we need to generate new window / slices? + if (slices.empty()) + { + // Set new window start time + while (window.endTime() < frameZero) + { + window.shiftStartTime(windowDelta); + printf("Propagated window forwards... new start time is %16.2f\n", sliceIt->first.startTime()); + } + + // Create the new slices + slices = prepareSlices(window, nSlices, inputNeXusFiles[0], outputFilePath); + sliceIt = slices.begin(); + } + + // If this frame zero is greater than or equal to the start time of the current window slice we can process + // events + if (frameZero >= sliceIt->first.startTime()) + { + // Sanity check! + if (frameZero > sliceIt->first.endTime()) + throw(std::runtime_error("Somebody's done something wrong here....\n")); + + // Grab the destination histograms for this slice and bin events + auto &destinationHistograms = sliceIt->second.detectorHistograms(); + for (auto i = 0; i < eventIndices.size(); ++i) + if (eventIndices[i] > 0) + destinationHistograms[eventIndices[i]].bin(eventTimes[i]); + + // Increment the frame counter for this slice + sliceIt->second.incrementGoodFrames(); + } + } + } + + fmt::print("Processing frame / event data... Done!\n"); + } + + // Perform post-processing on any slices we still have and save + postProcess(slices); + saveSlices(slices); +} + +} // namespace Processors diff --git a/src/partitionEventsSummed.cpp b/src/partitionEventsSummed.cpp new file mode 100644 index 0000000..9a2b632 --- /dev/null +++ b/src/partitionEventsSummed.cpp @@ -0,0 +1,97 @@ +#include "frameChunker.h" +#include "nexusFile.h" +#include "processors.h" +#include "window.h" +#include +#include + +namespace Processors +{ +// 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. + */ + + 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); + + // Initialise the slice iterator and window slice / NeXuSFile references + auto sliceIt = slices.begin(); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and get its event data + NeXuSFile nxs(nxsFileName, true); + + // Read in chunked frames + FrameChunker frameChunker(nxs); + auto tenPercentsDone = -1; + while (frameChunker.getNextFrameData()) + { + auto &frameData = frameChunker.frameData(); + auto currentTenPercents = int(100 * (frameData.front().index / double(nxs.nGoodFrames()))) / 10; + if (currentTenPercents > tenPercentsDone) + { + tenPercentsDone = currentTenPercents; + fmt::print("Processing frame / event data... {}%\n", tenPercentsDone * 10); + } + + for (const auto &frame : frameData) + { + const auto &eventIndices = frame.detectorIndices; + const auto &eventTimes = frame.times; + auto frameZero = frame.timeZero + nxs.startSinceEpoch(); + + // If the current slice end time is less than the frame zero, iterate the slices. + while (sliceIt->first.endTime() < frameZero) + { + sliceIt++; + + // If we have run out of slices propagate the set forward. + if (sliceIt == slices.end()) + { + sliceIt = slices.begin(); + for (auto &&[slice, _unused] : slices) + slice.shiftStartTime(windowDelta); + printf("Propagated window forwards... new start time is %16.2f\n", sliceIt->first.startTime()); + } + } + + // If this frame zero is greater than or equal to the start time of the current window slice we can process + // events + if (frameZero >= sliceIt->first.startTime()) + { + // Sanity check! + if (frameZero > sliceIt->first.endTime()) + throw(std::runtime_error("Somebody's done something wrong here....\n")); + + // Grab the destination histograms for this slice and bin events + auto &destinationHistograms = sliceIt->second.detectorHistograms(); + for (auto i = 0; i < eventIndices.size(); ++i) + if (eventIndices[i] > 0) + destinationHistograms[eventIndices[i]].bin(eventTimes[i]); + + // Increment the frame counter for this slice + sliceIt->second.incrementGoodFrames(); + } + } + } + + fmt::print("Processing frame / event data... Done!\n"); + } + + // Perform post-processing + postProcess(slices); + + // Save slices + saveSlices(slices); +} + +} // namespace Processors diff --git a/src/processCommon.cpp b/src/processCommon.cpp index 51f6aa7..2bc5327 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,15 @@ 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(); + nexus.loadMonitorCounts(); + nexus.zeroGoodFrames(); sliceStartTime += sliceDuration; } @@ -56,18 +65,18 @@ void postProcess(std::vector> &slices) for (auto &&[slice, outputNeXuSFile] : slices) { fmt::print("Output '{}' ({} -> {}) has {} detector frames and {} monitor frames.\n", std::string(slice.id()).c_str(), - slice.startTime(), slice.endTime(), outputNeXuSFile.nDetectorFrames(), outputNeXuSFile.nMonitorFrames()); + slice.startTime(), slice.endTime(), outputNeXuSFile.nGoodFrames(), outputNeXuSFile.nMonitorFrames()); switch (postProcessingMode_) { case (Processors::PostProcessingMode::None): break; case (Processors::PostProcessingMode::ScaleMonitors): - factor = (double)outputNeXuSFile.nDetectorFrames() / (double)outputNeXuSFile.nMonitorFrames(); + factor = (double)outputNeXuSFile.nGoodFrames() / (double)outputNeXuSFile.nMonitorFrames(); fmt::print(" --> Scaling monitors by processed detector-to-monitor frame ratio ({}).\n", factor); outputNeXuSFile.scaleMonitors(factor); break; case (Processors::PostProcessingMode::ScaleDetectors): - factor = (double)outputNeXuSFile.nMonitorFrames() / (double)outputNeXuSFile.nDetectorFrames(); + factor = (double)outputNeXuSFile.nMonitorFrames() / (double)outputNeXuSFile.nGoodFrames(); fmt::print(" --> Scaling detectors by processed monitor-to-detector frame ratio ({}).\n", factor); outputNeXuSFile.scaleDetectors(factor); break; diff --git a/src/processIndividual.cpp b/src/processIndividual.cpp deleted file mode 100644 index bb26620..0000000 --- a/src/processIndividual.cpp +++ /dev/null @@ -1,107 +0,0 @@ -#include "nexusFile.h" -#include "processors.h" -#include "window.h" -#include -#include - -namespace Processors -{ -// Perform summed processing -void processIndividual(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"); - - // Copy the master window definition since we will be modifying it - auto window = windowDefinition; - std::vector> slices; - auto sliceIt = slices.end(); - - // Loop over input Nexus files - for (auto &nxsFileName : inputNeXusFiles) - { - // Open the NeXuS file and get its event data - NeXuSFile nxs(nxsFileName, true); - - 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 - auto eventStart = 0, eventEnd = 0; - 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] + nxs.startSinceEpoch(); - - // If the current slice end time is less than the frame zero, iterate the slices. - while (sliceIt != slices.end() && sliceIt->first.endTime() < frameZero) - { - sliceIt++; - - // If we have run out of slices propagate the set forward. - if (sliceIt == slices.end()) - { - // Save current data - postProcess(slices); - saveSlices(slices); - - // Empty current slices - slices.clear(); - sliceIt = slices.end(); - } - } - - // Do we need to generate new window / slices? - if (slices.empty()) - { - // Set new window start time - while (window.endTime() < frameZero) - { - window.shiftStartTime(windowDelta); - printf("Propagated window forwards... new start time is %16.2f\n", sliceIt->first.startTime()); - } - - // Create the new slices - slices = prepareSlices(window, nSlices, inputNeXusFiles[0], outputFilePath); - sliceIt = slices.begin(); - } - - // If this frame zero is greater than or equal to the start time of the current window slice we can process events - if (frameZero >= sliceIt->first.startTime()) - { - // Sanity check! - if (frameZero > sliceIt->first.endTime()) - throw(std::runtime_error("Somebody's done something wrong here....\n")); - - // Grab the destination datafile for this slice and bin events - auto &destinationHistograms = sliceIt->second.detectorHistograms(); - for (int k = eventStart; k < eventEnd; ++k) - { - auto id = eventIndices[k]; - if (id > 0) - gsl_histogram_accumulate(destinationHistograms[id], eventTimes[k], 1.0); - } - - // Increment the frame counter for this slice - sliceIt->second.incrementDetectorFrameCount(); - } - - // Update start event index - eventStart = eventEnd; - } - } - - // Perform post-processing on any slices we still have and save - postProcess(slices); - saveSlices(slices); -} - -} // namespace Processors diff --git a/src/processSummed.cpp b/src/processSummed.cpp deleted file mode 100644 index 6255fdd..0000000 --- a/src/processSummed.cpp +++ /dev/null @@ -1,91 +0,0 @@ -#include "nexusFile.h" -#include "processors.h" -#include "window.h" -#include - -namespace Processors -{ -// Perform summed processing -void processSummed(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"); - - // Generate a new set of window "slices" and associated output NeXuS files to sum data into - auto slices = prepareSlices(windowDefinition, nSlices, inputNeXusFiles[0], outputFilePath); - - // Initialise the slice iterator and window slice / NeXuSFile references - auto sliceIt = slices.begin(); - - // Loop over input Nexus files - for (auto &nxsFileName : inputNeXusFiles) - { - // Open the NeXuS file and get its event data - NeXuSFile nxs(nxsFileName, true); - - 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 - auto eventStart = 0, eventEnd = 0; - 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] + nxs.startSinceEpoch(); - - // If the current slice end time is less than the frame zero, iterate the slices. - while (sliceIt->first.endTime() < frameZero) - { - sliceIt++; - - // If we have run out of slices propagate the set forward. - if (sliceIt == slices.end()) - { - sliceIt = slices.begin(); - for (auto &&[slice, _unused] : slices) - slice.shiftStartTime(windowDelta); - printf("Propagated window forwards... new start time is %16.2f\n", sliceIt->first.startTime()); - } - } - - // If this frame zero is greater than or equal to the start time of the current window slice we can process events - if (frameZero >= sliceIt->first.startTime()) - { - // Sanity check! - if (frameZero > sliceIt->first.endTime()) - throw(std::runtime_error("Somebody's done something wrong here....\n")); - - // Grab the destination datafile for this slice and bin events - auto &destinationHistograms = sliceIt->second.detectorHistograms(); - for (int k = eventStart; k < eventEnd; ++k) - { - auto id = eventIndices[k]; - if (id > 0) - gsl_histogram_accumulate(destinationHistograms[id], eventTimes[k], 1.0); - } - - // Increment the frame counter for this slice - sliceIt->second.incrementDetectorFrameCount(); - } - - // Update start event index - eventStart = eventEnd; - } - } - - // Perform post-processing - postProcess(slices); - - // Save slices - saveSlices(slices); -} - -} // namespace Processors diff --git a/src/processors.h b/src/processors.h index 6ed0831..3a8489c 100644 --- a/src/processors.h +++ b/src/processors.h @@ -13,8 +13,16 @@ namespace Processors enum class ProcessingMode { None, - Individual, - Summed + CountDetector, + CountMonitor, + DumpDetector, + DumpDetectorFromEvents, + DumpEvents, + DumpMonitor, + PartitionEventsIndividual, + PartitionEventsSummed, + PrintEvents, + ResizeDetectors }; // Processing Direction @@ -51,13 +59,24 @@ void saveSlices(std::vector> &slices); * Processors */ -// Get Events -std::map> getEvents(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); +// Count detector histogram +void countDetector(const std::vector &inputNeXusFiles, int detectorIndex); +// Count monitor histogram +void countMonitor(const std::vector &inputNeXusFiles, int monitorIndex); +// 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); +// Dump detector histogram after constructing it from event data +void dumpDetectorFromEvents(const std::vector &inputNeXusFiles, int detectorIndex); +// Dump monitor histogram +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); +// Partition events into summed windows / slices +void partitionEventsSummed(const std::vector &inputNeXusFiles, std::string_view outputFilePath, + const Window &windowDefinition, int nSlices, double windowDelta); +// Resize detector array +void resizeDetectors(const std::vector &inputNeXusFiles, int newNDetectors); }; // namespace Processors diff --git a/src/resizeDetectors.cpp b/src/resizeDetectors.cpp new file mode 100644 index 0000000..a345064 --- /dev/null +++ b/src/resizeDetectors.cpp @@ -0,0 +1,52 @@ +#include "nexusFile.h" +#include "processors.h" +#include +#include + +namespace Processors +{ +void resizeDetectors(const std::vector &inputNeXusFiles, int totalDetectorCount) +{ + /* + * Resize detector space, deleting extra or adding empty as necessary + */ + + fmt::print("Resizing number of detector spectra to {}...\n", totalDetectorCount); + + // Loop over input NeXuS files + for (auto &nxsFileName : inputNeXusFiles) + { + // Open the NeXuS file and load in detector counts + NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); + nxs.loadDetectorCounts(); + + if (nxs.nDetectors() == totalDetectorCount) + { + fmt::print("Skipping - file already contains the correct number of detectors.\n"); + continue; + } + + // Adjust detector space size + while (nxs.nDetectors() != totalDetectorCount) + { + if (nxs.nDetectors() > totalDetectorCount) + { + // Remove the last (highest index) detector spectrum + auto specID = nxs.removeLastDetector(); + fmt::print(" ... Removed detector {} (spectrum {})\n", specID, specID - 1); + } + else + { + // Add a new detector on the end + auto specID = nxs.appendEmptyDetector(); + fmt::print(" ... Appended empty detector {} (spectrum {})\n", specID, specID - 1); + } + } + + // Write the new data + nxs.saveModifiedData(); + } +} + +} // namespace Processors