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..9808db5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -44,6 +44,7 @@ if(LOCAL_STATIC_GSL) else(LOCAL_STATIC_GSL) find_package(GSL REQUIRED) list(APPEND LINK_LIBS "${GSL_LIBRARIES}") + include_directories("${GSL_INCLUDE_DIRS}") endif(LOCAL_STATIC_GSL) option( diff --git a/np.cpp b/np.cpp index f609454..550f6e6 100644 --- a/np.cpp +++ b/np.cpp @@ -27,45 +27,98 @@ int main(int argc, char **argv) // Number of slices to partition window in to int windowSlices_{1}; // Target spectrum for event get (optional) - std::optional spectrumId_; + int targetIndex_; // Define and parse CLI arguments CLI::App app("NeXuS Processor (np), Copyright (C) 2024 Jared Swift and Tristan Youngs."); // -- Window Definition app.add_option("-n,--name", windowName_, "Name of the window, used as a prefix to all output files") - ->group("Window Definition") - ->required(); + ->group("Window Definition"); app.add_option("-s,--start", windowStartTime_, - "Start time of the window (relative to first input file start time unless --absolute-start is given)") + "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; + targetIndex_ = 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; + targetIndex_ = id; + }, + "Print all events for specified detector index") + ->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; + targetIndex_ = id; + }, + "Dump specified detector histogram") + ->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; + targetIndex_ = 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,18 +126,15 @@ 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)") - ->group("Processing"); // -- Post Processing app.add_flag_callback( "--scale-monitors", [&]() { Processors::postProcessingMode_ = Processors::PostProcessingMode::ScaleMonitors; }, @@ -97,53 +147,67 @@ int main(int argc, char **argv) CLI11_PARSE(app, argc, argv); - // Sanity check - if ((windowWidth_ + windowOffset_) > windowDelta_) - { - fmt::print("Error: Window width (including any optional offset) is greater than window delta.\n"); - return 1; - } - if (windowSlices_ < 1) - { - fmt::print("Error: Invalid number of window slices provided ({}).\n", windowSlices_); - return 1; - } - - // Perform pre-processing if requested - if (spectrumId_) - { - Processors::getEvents(inputFiles_, *spectrumId_); - } - - // Construct the master window definition - if (relativeStartTime_) - { - // Need to query first NeXuS file to get its start time - if (inputFiles_.empty()) - { - fmt::print("Error: Need at least one input NeXuS file."); - return 1; - } - NeXuSFile firstFile(inputFiles_.front()); - firstFile.loadTimes(); - fmt::print("Window start time converted from relative to absolute time: {} => {} (= {} + {})\n", windowStartTime_, - windowStartTime_ + firstFile.startSinceEpoch(), firstFile.startSinceEpoch(), windowStartTime_); - windowStartTime_ += firstFile.startSinceEpoch(); - } - Window window(windowName_, windowStartTime_ + windowOffset_, windowWidth_); - fmt::print("Window start time (including any offset) is {}.\n", window.startTime()); - // Perform processing switch (processingMode_) { case (Processors::ProcessingMode::None): - fmt::print("No processing mode specified. We are done.\n"); + 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_, targetIndex_, + processingMode_ == Processors::ProcessingMode::PrintEvents); break; - case (Processors::ProcessingMode::Individual): - Processors::processIndividual(inputFiles_, outputDirectory_, window, windowSlices_, windowDelta_); + case (Processors::ProcessingMode::DumpDetector): + Processors::dumpDetector(inputFiles_, targetIndex_); break; - case (Processors::ProcessingMode::Summed): - Processors::processSummed(inputFiles_, outputDirectory_, window, windowSlices_, windowDelta_); + case (Processors::ProcessingMode::DumpMonitor): + Processors::dumpMonitor(inputFiles_, targetIndex_); + 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; default: throw(std::runtime_error("Unhandled processing mode.\n")); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2b200fa..80c607b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,16 +1,19 @@ -add_library(nexusProcess +add_library( + nexusProcess + dumpDetector.cpp + dumpMonitor.cpp getEvents.cpp nexusFile.cpp + partitionEventsIndividual.cpp + partitionEventsSummed.cpp processCommon.cpp - processIndividual.cpp - processSummed.cpp window.cpp 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/dumpDetector.cpp b/src/dumpDetector.cpp new file mode 100644 index 0000000..b5c90a2 --- /dev/null +++ b/src/dumpDetector.cpp @@ -0,0 +1,39 @@ +#include "nexusFile.h" +#include "processors.h" +#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 + std::ofstream output(fmt::format("{}.det.{}", nxsFileName, detectorIndex).c_str()); + output << fmt::format("# TCB/us Counts [detector index {}, spectrum index = {}]\n", detectorIndex, spectrumId); + auto bin = 0; + const auto &counts = nxs.detectorCounts().at(spectrumId); + for (auto tof : nxs.tofBoundaries()) + { + output << fmt::format("{} {}\n", tof, counts[bin++]); + } + output.close(); + } +} + +} // namespace Processors diff --git a/src/dumpMonitor.cpp b/src/dumpMonitor.cpp new file mode 100644 index 0000000..972c67d --- /dev/null +++ b/src/dumpMonitor.cpp @@ -0,0 +1,37 @@ +#include "nexusFile.h" +#include "processors.h" +#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 + std::ofstream output(fmt::format("{}.mon.{}", nxsFileName, monitorIndex).c_str()); + output << "# TCB/us Counts\n"; + auto bin = 0; + const auto &counts = nxs.monitorCounts().at(monitorIndex); + for (auto tof : nxs.tofBoundaries()) + { + output << fmt::format("{} {}\n", tof, counts[bin++]); + } + output.close(); + } +} + +} // namespace Processors diff --git a/src/getEvents.cpp b/src/getEvents.cpp index 0d99163..4b518ae 100644 --- a/src/getEvents.cpp +++ b/src/getEvents.cpp @@ -2,39 +2,47 @@ #include "processors.h" #include "window.h" #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); + fmt::print("Retrieving all events from detector index {}...\n", detectorIndex); - std::map> eventMap; - std::optional lastSecondsSinceEpoch; - - // Loop over input Nexus files + // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) { - // Open the Nexus file ready for use + // Open the NeXuS file ready for use NeXuSFile nxs(nxsFileName); + nxs.prepareSpectraSpace(); nxs.loadEventData(); - nxs.loadTimes(); - fmt::print("... file '{}' has {} events...\n", nxsFileName, nxs.eventTimes().size()); + std::optional lastSecondsSinceEpoch; 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(); + 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; + output << fmt::format("# {:20s} {:20s} {:20s} {}\n", "frame_offset(us)", "start_time_offset(s)", "epoch_offset(s)", + "delta(s)"); - // Loop over frames in the Nexus file + // 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 @@ -49,23 +57,23 @@ std::map> getEvents(const std::vector &inp auto eSeconds = eMicroSeconds * 0.000001; auto eSecondsSinceEpoch = eSeconds + frameZero + nxs.startSinceEpoch(); if (lastSecondsSinceEpoch) - fmt::print("{:20.6f} {:20.10f} {:20.5f} {}\n", eMicroSeconds, eSeconds + frameZero, - eSecondsSinceEpoch, eSecondsSinceEpoch - *lastSecondsSinceEpoch); + output << fmt::format("{:20.6f} {:20.10f} {:20.5f} {}\n", eMicroSeconds, eSeconds + frameZero, + eSecondsSinceEpoch, eSecondsSinceEpoch - *lastSecondsSinceEpoch); else - fmt::print("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, eSecondsSinceEpoch); - eventMap[spectrumId].push_back(eSecondsSinceEpoch); + output << fmt::format("{:20.6f} {:20.10f} {:20.5f}\n", eMicroSeconds, eSeconds + frameZero, + eSecondsSinceEpoch); + lastSecondsSinceEpoch = eSecondsSinceEpoch; - if (firstOnly) - return eventMap; } } // Update start event index eventStart = eventEnd; } - } - return eventMap; + if (!toStdOut) + fileOutput.close(); + } } } // namespace Processors 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..e5c7221 100644 --- a/src/nexusFile.cpp +++ b/src/nexusFile.cpp @@ -9,10 +9,10 @@ std::vector neXuSBasicPaths_ = {"/raw_data_1/title", "/raw_data_1/user_1/name", "/raw_data_1/start_time", + "/raw_data_1/end_time", "/raw_data_1/good_frames", "/raw_data_1/raw_frames", "/raw_data_1/monitor_1/data", - "/raw_data_1/monitor_1/time_of_flight", "/raw_data_1/monitor_2/data", "/raw_data_1/monitor_3/data", "/raw_data_1/monitor_4/data", @@ -21,18 +21,87 @@ std::vector neXuSBasicPaths_ = {"/raw_data_1/title", "/raw_data_1/monitor_7/data", "/raw_data_1/monitor_8/data", "/raw_data_1/monitor_9/data", - "/raw_data_1/detector_1/counts"}; + "/raw_data_1/detector_1/counts", + "/raw_data_1/detector_1/spectrum_index", + "/raw_data_1/detector_1/time_of_flight"}; -NeXuSFile::NeXuSFile(std::string filename, 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_; + nDetectorFrames_ = source.nDetectorFrames_; + nGoodFrames_ = source.nGoodFrames_; + startSinceEpoch_ = source.startSinceEpoch_; + endSinceEpoch_ = source.endSinceEpoch_; + eventIndices_ = source.eventIndices_; + eventTimes_ = source.eventTimes_; + eventsPerFrame_ = source.eventsPerFrame_; + frameOffsets_ = source.frameOffsets_; + tofBoundaries_ = source.tofBoundaries_; + monitorCounts_ = source.monitorCounts_; + detectorCounts_ = source.detectorCounts_; + if (deepCopyHistograms) + { + detectorHistograms_.clear(); + for (auto &[specId, histogram] : detectorHistograms_) + detectorHistograms_[specId] = gsl_histogram_clone(histogram); } + else + detectorHistograms_ = source.detectorHistograms_; +} + +// Clear all data and arrays +void NeXuSFile::clear() +{ + filename_ = {}; + detectorSpectrumIndices_.clear(); + nMonitorSpectra_ = 0; + nMonitorFrames_ = 0; + nDetectorFrames_ = 0; + nGoodFrames_ = 0; + startSinceEpoch_ = 0; + endSinceEpoch_ = 0; + eventIndices_.clear(); + eventTimes_.clear(); + eventsPerFrame_.clear(); + frameOffsets_.clear(); + tofBoundaries_.clear(); + monitorCounts_.clear(); + detectorCounts_.clear(); + // Free GSL histograms + for (auto &[specId, histogram] : detectorHistograms_) + gsl_histogram_free(histogram); + detectorHistograms_.clear(); } /* @@ -62,18 +131,81 @@ std::pair NeXuSFile::find1DDataset(H5::H5File file, H5std // Return filename std::string NeXuSFile::filename() const { return filename_; } -// Template basic paths from the referenceFile, and make ready for histogram binning -void NeXuSFile::templateFile(std::string referenceFile, std::string outputFile) +// Load basic information from the NeXuS file +void NeXuSFile::loadBasicData(bool printInfo) +{ + hid_t memType = H5Tcopy(H5T_C_S1); + H5Tset_size(memType, UCHAR_MAX); + char charBuffer[UCHAR_MAX]; + + // Open input NeXuS file in read only mode. + H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); + + // Get the run title + auto &&[titleID, titleDimension] = NeXuSFile::find1DDataset(input, "raw_data_1", "title"); + H5Dread(titleID.getId(), memType, H5S_ALL, H5S_ALL, H5P_DEFAULT, charBuffer); + 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::find1DDataset(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::find1DDataset(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::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]; + 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) { - filename_ = outputFile; + // Open this NeXuS file in read only mode. + H5::H5File input = H5::H5File(sourceFilename, H5F_ACC_RDONLY); - // Open input Nexus file in read only mode. - H5::H5File input = H5::H5File(referenceFile, H5F_ACC_RDONLY); + // Create new NeXuS file for output. + H5::H5File output = H5::H5File(newFilename, H5F_ACC_TRUNC); - // Create new Nexus file for output. - H5::H5File output = H5::H5File(filename_, H5F_ACC_TRUNC); + printf("Templating file '%s' to '%s'...\n", sourceFilename.c_str(), newFilename.c_str()); - printf("Templating file '%s' to '%s'...\n", referenceFile.c_str(), filename_.c_str()); hid_t ocpl_id, lcpl_id; ocpl_id = H5Pcreate(H5P_OBJECT_COPY); if (ocpl_id < 0) @@ -93,61 +225,73 @@ 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::find1DDataset(input, "raw_data_1/detector_1", "spectrum_index"); + detectorSpectrumIndices_.resize(spectraDimension); + H5Dread(detSpecIndices.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, detectorSpectrumIndices_.data()); + nMonitorSpectra_ = detectorSpectrumIndices_.front() - 1; - // 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; + // Read in TOF boundary information. + auto &&[tofBoundariesID, tofBoundariesDimension] = + NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "time_of_flight"); + tofBoundaries_.resize(tofBoundariesDimension); + H5Dread(tofBoundariesID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, tofBoundaries_.data()); - monitorCounts_[i].resize(tofBinsDimension); - H5Dread(monitorSpectrum.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, monitorCounts_[i].data()); + // Set up detector histograms and straight counts vectors + for (auto spec : detectorSpectrumIndices_) + { + // For event re-binning + detectorHistograms_[spec] = gsl_histogram_alloc(tofBoundaries_.size() - 1); + gsl_histogram_set_ranges(detectorHistograms_[spec], tofBoundaries_.data(), tofBoundaries_.size()); - ++i; + // For histogram manipulation + detectorCounts_[spec].resize(tofBoundaries_.size() - 1); + std::fill(detectorCounts_[spec].begin(), detectorCounts_[spec].end(), 0); } - // 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]; - 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 + // Read in monitor data - start from index 1 and end when we fail to find the named dataset with this suffix + for (auto i = 0; i > nMonitorSpectra_; ++i) + { + auto &&[monitorSpectrum, monitorSpectrumDimension] = + NeXuSFile::find1DDataset(input, "/raw_data_1/monitor_" + std::to_string(i + 1), "data"); + + monitorCounts_[i].resize(tofBoundaries_.size()); + H5Dread(monitorSpectrum.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, monitorCounts_[i].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::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]; + nMonitorFrames_ = goodFramesTemp[0]; input.close(); } @@ -155,9 +299,9 @@ void NeXuSFile::loadFrameCounts() // Load event data void NeXuSFile::loadEventData() { - printf("Load event data...\n"); + printf("Loading event data...\n"); - // Open our Nexus file in read only mode. + // Open our NeXuS file in read only mode. H5::H5File input = H5::H5File(filename_, H5F_ACC_RDONLY); // Read in event indices. @@ -184,58 +328,38 @@ void NeXuSFile::loadEventData() frameOffsets_.resize(frameOffsetsDimension); H5Dread(frameOffsetsID.getId(), H5T_IEEE_F64LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, frameOffsets_.data()); + fmt::print("... there are {} events...\n", eventTimes_.size()); + input.close(); } -// 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); - - 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); + 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 etime = {0}; - etime.tm_year = y - 1900; - etime.tm_mon = M - 1; - etime.tm_mday = d; - etime.tm_hour = h; - etime.tm_min = m; - etime.tm_sec = s; - - endSinceEpoch_ = (int)mktime(&etime); - - input.close(); + std::vector countsBuffer; + countsBuffer.resize(nSpec * nTOFBins); // HDF5 expects contiguous memory. This is a pain. + auto &&[counts, detectorCountsDimension] = NeXuSFile::find1DDataset(input, "raw_data_1/detector_1", "counts"); + H5Dread(counts.getId(), H5T_STD_I32LE, H5S_ALL, H5S_ALL, H5P_DEFAULT, countsBuffer.data()); + for (auto i = 0; i < nSpec; ++i) + { + auto &counts = detectorCounts_[detectorSpectrumIndices_[i]]; + for (auto j = 0; j < nTOFBins; ++j) + counts[j] = countsBuffer[i * nTOFBins + j]; + } } // Save key modified data back to the file bool NeXuSFile::saveModifiedData() { - // Open Nexus file in read/write mode. + // Open NeXuS file in read/write mode. H5::H5File output = H5::H5File(filename_, H5F_ACC_RDWR); // Write good frames @@ -253,13 +377,13 @@ bool NeXuSFile::saveModifiedData() } // Write detector counts - const auto nSpec = spectra_.size(); - const auto nTofBins = tofBins_.size() - 1; + const auto nSpec = detectorSpectrumIndices_.size(); + const auto ntofBoundaries = tofBoundaries_.size() - 1; - auto *countsBuffer = new int[nSpec * nTofBins]; // HDF5 expects contiguous memory. This is a pain. + auto *countsBuffer = new int[nSpec * ntofBoundaries]; // HDF5 expects contiguous memory. This is a pain. for (auto i = 0; i < nSpec; ++i) - for (auto j = 0; j < nTofBins; ++j) - countsBuffer[i * nTofBins + j] = gsl_histogram_get(detectorHistograms_[spectra_[i]], j); + for (auto j = 0; j < ntofBoundaries; ++j) + countsBuffer[i * ntofBoundaries + j] = gsl_histogram_get(detectorHistograms_[detectorSpectrumIndices_[i]], j); auto &&[counts, detectorCountsDimension] = NeXuSFile::find1DDataset(output, "raw_data_1/detector_1", "counts"); counts.write(countsBuffer, H5::PredType::STD_I32LE); @@ -284,10 +408,11 @@ const std::vector &NeXuSFile::eventIndices() const { return eventIndices_; const std::vector &NeXuSFile::eventTimes() const { return eventTimes_; } const std::vector &NeXuSFile::eventsPerFrame() const { return eventsPerFrame_; } const std::vector &NeXuSFile::frameOffsets() const { return frameOffsets_; } -const std::vector &NeXuSFile::tofBins() const { return tofBins_; } +const std::vector &NeXuSFile::tofBoundaries() const { return tofBoundaries_; } +const int NeXuSFile::spectrumForDetector(int detectorId) const { return detectorSpectrumIndices_.at(detectorId - 1); } const std::map> &NeXuSFile::monitorCounts() const { return monitorCounts_; } +const std::map> &NeXuSFile::detectorCounts() const { return detectorCounts_; } std::map &NeXuSFile::detectorHistograms() { return detectorHistograms_; } -const std::map> &NeXuSFile::partitions() const { return partitions_; } /* * Manipulation @@ -312,7 +437,7 @@ void NeXuSFile::scaleMonitors(double factor) void NeXuSFile::scaleDetectors(double factor) { auto oldSum = 0, newSum = 0; - for (auto i : spectra_) + for (auto i : detectorSpectrumIndices_) { oldSum += gsl_histogram_sum(detectorHistograms_[i]); gsl_histogram_scale(detectorHistograms_[i], factor); diff --git a/src/nexusFile.h b/src/nexusFile.h index be4922b..0052bb3 100644 --- a/src/nexusFile.h +++ b/src/nexusFile.h @@ -9,8 +9,15 @@ 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 @@ -26,14 +33,18 @@ class NeXuSFile 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 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 event data void loadEventData(); - // Load start/end times - void loadTimes(); + // Load detector counts from the file + void loadDetectorCounts(); // Save key modified data back to the file bool saveModifiedData(); @@ -41,7 +52,8 @@ class NeXuSFile * Data */ private: - std::vector spectra_; + std::vector detectorSpectrumIndices_; + int nMonitorSpectra_{0}; int nMonitorFrames_{0}; int nDetectorFrames_{0}; int nGoodFrames_{0}; @@ -51,10 +63,10 @@ class NeXuSFile std::vector eventTimes_; std::vector eventsPerFrame_; std::vector frameOffsets_; - std::vector tofBins_; + std::vector tofBoundaries_; std::map> monitorCounts_; + std::map> detectorCounts_; std::map detectorHistograms_; - std::map> partitions_; public: [[nodiscard]] int nGoodFrames() const; @@ -67,10 +79,11 @@ class NeXuSFile [[nodiscard]] const std::vector &eventTimes() const; [[nodiscard]] const std::vector &eventsPerFrame() const; [[nodiscard]] const std::vector &frameOffsets() const; - [[nodiscard]] const std::vector &tofBins() const; + [[nodiscard]] const std::vector &tofBoundaries() const; + [[nodiscard]] const int spectrumForDetector(int detectorId) const; [[nodiscard]] const std::map> &monitorCounts() const; + [[nodiscard]] const std::map> &detectorCounts() const; std::map &detectorHistograms(); - [[nodiscard]] const std::map> &partitions() const; /* * Manipulation diff --git a/src/processIndividual.cpp b/src/partitionEventsIndividual.cpp similarity index 88% rename from src/processIndividual.cpp rename to src/partitionEventsIndividual.cpp index bb26620..aca4e74 100644 --- a/src/processIndividual.cpp +++ b/src/partitionEventsIndividual.cpp @@ -6,34 +6,36 @@ namespace Processors { -// Perform summed processing -void processIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, - const Window &windowDefinition, int nSlices, double windowDelta) +// Partition events into individual windows / slices +void partitionEventsIndividual(const std::vector &inputNeXusFiles, std::string_view outputFilePath, + const Window &windowDefinition, int nSlices, double windowDelta) { /* * From our main windowDefinition we will continually propagate it forwards in time (by the window delta) splitting it into * nSlices and until we go over the end time of the current file. */ - fmt::print("Processing in INDIVIDUAL mode...\n"); + fmt::print("Partitioning events into individual windows/slices...\n"); // Copy the master window definition since we will be modifying it auto window = windowDefinition; std::vector> slices; auto sliceIt = slices.end(); - // Loop over input Nexus files + // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) { // Open the NeXuS file and get its event data NeXuSFile nxs(nxsFileName, true); + nxs.prepareSpectraSpace(); + nxs.loadEventData(); 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 + // Loop over frames in the NeXuS file auto eventStart = 0, eventEnd = 0; for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) { diff --git a/src/processSummed.cpp b/src/partitionEventsSummed.cpp similarity index 87% rename from src/processSummed.cpp rename to src/partitionEventsSummed.cpp index 6255fdd..92e6f6d 100644 --- a/src/processSummed.cpp +++ b/src/partitionEventsSummed.cpp @@ -1,20 +1,21 @@ #include "nexusFile.h" #include "processors.h" #include "window.h" +#include #include namespace Processors { -// Perform summed processing -void processSummed(const std::vector &inputNeXusFiles, std::string_view outputFilePath, - const Window &windowDefinition, int nSlices, double windowDelta) +// Partition events into summed windows / slices +void partitionEventsSummed(const std::vector &inputNeXusFiles, std::string_view outputFilePath, + const Window &windowDefinition, int nSlices, double windowDelta) { /* * From our main windowDefinition we will continually propagate it forwards in time (by the window delta) splitting it into * nSlices and until we go over the end time of the current file. */ - printf("Processing in SUMMED mode...\n"); + fmt::print("Partitioning events in summed windows/slices...\n"); // Generate a new set of window "slices" and associated output NeXuS files to sum data into auto slices = prepareSlices(windowDefinition, nSlices, inputNeXusFiles[0], outputFilePath); @@ -22,18 +23,19 @@ void processSummed(const std::vector &inputNeXusFiles, std::string_ // Initialise the slice iterator and window slice / NeXuSFile references auto sliceIt = slices.begin(); - // Loop over input Nexus files + // Loop over input NeXuS files for (auto &nxsFileName : inputNeXusFiles) { // Open the NeXuS file and get its event data NeXuSFile nxs(nxsFileName, true); + nxs.loadEventData(); 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 + // Loop over frames in the NeXuS file auto eventStart = 0, eventEnd = 0; for (auto frameIndex = 0; frameIndex < nxs.eventsPerFrame().size(); ++frameIndex) { diff --git a/src/processCommon.cpp b/src/processCommon.cpp index 51f6aa7..1aa7214 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,14 @@ 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(); sliceStartTime += sliceDuration; } diff --git a/src/processors.h b/src/processors.h index 6ed0831..0765491 100644 --- a/src/processors.h +++ b/src/processors.h @@ -13,8 +13,12 @@ namespace Processors enum class ProcessingMode { None, - Individual, - Summed + DumpDetector, + DumpEvents, + DumpMonitor, + PartitionEventsIndividual, + PartitionEventsSummed, + PrintEvents }; // Processing Direction @@ -51,13 +55,16 @@ 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); +// 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 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); }; // namespace Processors