From 909fa6519da3aadf157fa4c251179026834ba1c6 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 25 Jan 2026 10:51:23 +0100 Subject: [PATCH] added HDF4 ISIS NeXus class handling. --- CMakeLists.txt | 30 +- cmake/FindHDF4.cmake | 97 -- cmake/FindMXML.cmake | 34 - cmake/FindNEXUS.cmake | 45 - src/external/nexus/CMakeLists.txt | 2 +- src/external/nexus/PNeXus.cpp | 2193 +++++++++++++++++++++++++++++ src/external/nexus/PNeXus.h | 848 +++++++++++ 7 files changed, 3071 insertions(+), 178 deletions(-) delete mode 100644 cmake/FindHDF4.cmake delete mode 100644 cmake/FindMXML.cmake delete mode 100644 cmake/FindNEXUS.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index f06d59c4..1c095d2e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -213,7 +213,35 @@ endif (qt_based_tools) if (nexus) find_package(HDF5 COMPONENTS CXX REQUIRED) if (HAVE_HDF4) - find_package(HDF4 REQUIRED) + #--- check for HDF4 ----------------------------------------------------------- + # Find HDF4 manually (pkg-config often doesn't have hdf4) + find_path(HDF4_INCLUDE_DIR + NAMES mfhdf.h + PATHS /usr/include /usr/local/include + PATH_SUFFIXES hdf + ) + + find_library(HDF4_DF_LIBRARY + NAMES df libdf + PATHS /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib + ) + + find_library(HDF4_MFHDF_LIBRARY + NAMES mfhdf libmfhdf + PATHS /usr/lib64 /usr/lib /usr/local/lib64 /usr/local/lib + ) + + if (HDF4_INCLUDE_DIR AND HDF4_DF_LIBRARY AND HDF4_MFHDF_LIBRARY) + set(HDF4_FOUND TRUE) + set(HDF4_INCLUDE_DIRS ${HDF4_INCLUDE_DIR}) + set(HDF4_LIBRARIES ${HDF4_MFHDF_LIBRARY} ${HDF4_DF_LIBRARY}) + message(STATUS "Found HDF4: ${HDF4_INCLUDE_DIR}") + message(STATUS " HDF4 libraries: ${HDF4_LIBRARIES}") + else() + message(FATAL_ERROR "HDF4 library not found. Please install libhdf4-dev or hdf-devel") + endif() + + include_directories(${HDF4_INCLUDE_DIRS}) add_definitions(-DHAVE_HDF4) endif (HAVE_HDF4) add_definitions(-DPNEXUS_ENABLED) diff --git a/cmake/FindHDF4.cmake b/cmake/FindHDF4.cmake deleted file mode 100644 index 9026459a..00000000 --- a/cmake/FindHDF4.cmake +++ /dev/null @@ -1,97 +0,0 @@ -## Process this file with cmake -#============================================================================= -# NeXus - Neutron & X-ray Common Data Format -# -# CMakeLists for building the NeXus library and applications. -# -# Copyright (C) 2011 Stephen Rankin -# -# This library is free software; you can redistribute it and/or modify it -# under the terms of the GNU Lesser General Public License as published by the -# Free Software Foundation; either version 2 of the License, or (at your -# option) any later version. -# -# This library is distributed in the hope that it will be useful, but WITHOUT -# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License -# for more details. -# -# You should have received a copy of the GNU Lesser General Public License -# along with this library; if not, write to the Free Software Foundation, -# Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -# -# For further information, see -# -# -#============================================================================= - - -#------------------------------------------------------------------------------ -# find the runtime binaries of the HDF4 library -#------------------------------------------------------------------------------ -find_library(HDF4_DF_LIBRARY NAMES df hdf - HINTS ENV HDF4_ROOT - PATH_SUFFIXES hdf) - - -if(HDF4_DF_LIBRARY MATCHES HDF4_DF_LIBRARY-NOTFOUND) - message(FATAL_ERROR "Could not find HDF4 DF library!") -else() - get_filename_component(HDF4_LIBRARY_DIRS ${HDF4_DF_LIBRARY} PATH) - message(STATUS "Found HDF4 DF library: ${HDF4_DF_LIBRARY}") - message(STATUS "HDF4 libary path: ${HDF4_LIBRARY_DIRS}") -endif() - -find_library(HDF4_MFHDF_LIBRARY NAMES mfhdf - HINTS ENV HDF4_ROOT - PATH_SUFFIXES hdf) - -if(HDF4_MFHDF_LIBRARY MATCHES HDF4_MFHDF_LIBRARY-NOTFOUND) - message(FATAL_ERROR "Could not find HDF5 MFHDF library!") -else() - message(STATUS "Found HDF4 MFHDF library: ${HDF4_MFHDF_LIBRARY}") -endif() - - -#------------------------------------------------------------------------------ -# find the HDF4 header file -#------------------------------------------------------------------------------ -find_path(HDF4_INCLUDE_DIRS mfhdf.h - HINTS ENV HDF4_ROOT - PATH_SUFFIXES hdf) - -if(HDF4_INCLUDE_DIRS MATCHES HDF4_INCLUDE_DIRS-NOTFOUND) - message(FATAL_ERROR "Could not find HDF4 header files") -else() - message(STATUS "Found HDF4 header files in: ${HDF4_INCLUDE_DIRS}") -endif() - -#------------------------------------------------------------------------------ -# search for additional packages required to link against HDF4 -#------------------------------------------------------------------------------ -find_package(JPEG REQUIRED) - -#------------------------------------------------------------------------------ -# add libraries to the link list for NAPI -#------------------------------------------------------------------------------ -get_filename_component(LIB_EXT ${HDF4_DF_LIBRARY} EXT) -if(LIB_EXT MATCHES .a) - message(STATUS "HDF4 DF library is static") - list(APPEND NAPI_LINK_LIBS "-Wl,-whole-archive" ${HDF4_DF_LIBRARY} "-Wl,-no-whole-archive") -else() - list(APPEND NAPI_LINK_LIBS ${HDF4_DF_LIBRARY}) -endif() - - -get_filename_component(LIB_EXT ${HDF4_MFHDF_LIBRARY} EXT) -if(LIB_EXT MATCHES .a) - message(STATUS "HDF4 MFHDF library is static") - list(APPEND NAPI_LINK_LIBS "-Wl,-whole-archive" ${HDF4_MFHDF_LIBRARY} "-Wl,-no-whole-archive") -else() - list(APPEND NAPI_LINK_LIBS ${HDF4_MFHDF_LIBRARY}) -endif() - -list(APPEND NAPI_LINK_LIBS jpeg) - -include_directories ( SYSTEM ${HDF4_INCLUDE_DIRS} ) -link_directories(${HDF4_LIBRARY_DIRS}) diff --git a/cmake/FindMXML.cmake b/cmake/FindMXML.cmake deleted file mode 100644 index e11a62fe..00000000 --- a/cmake/FindMXML.cmake +++ /dev/null @@ -1,34 +0,0 @@ -# - find MXML -# find the MXML lib and includes -# This module defines -# LIBMXML_INCLUDE_DIR, where to find mxml.h -# LIBMXML_LIBRARY, library to link against -# LIBMXML_FOUND, if false, do not try to use the MXML lib - -find_path(LIBMXML_INCLUDE_DIR mxml.h - HINT "/usr/include" -) -# find position of mxml.h from the end -string(FIND "${LIBMXML_INCLUDE_DIR}" "/mxml.h" pos REVERSE) -# truncate the string -string(SUBSTRING "${LIBMXML_INCLUDE_DIR}" 0 ${pos} substr) -set(LIBMXML_INCLUDE_DIR ${substr}) -unset(substr) - -find_library(LIBMXML_LIBRARY mxml) - -# get version string -# currently do not know from where to get it automatically - -# handle the QUIETLY and REQUIRED arguments and set LIBMXML_FOUND to TRUE if -# all listed variables are TRUE -include(${CMAKE_ROOT}/Modules/FindPackageHandleStandardArgs.cmake) -FIND_PACKAGE_HANDLE_STANDARD_ARGS(MXML - REQUIRED_VARS LIBMXML_LIBRARY LIBMXML_INCLUDE_DIR) - -if (NOT LIBMXML_FOUND) - unset(LIBMXML_LIBRARY) -endif() - -mark_as_advanced(LIBMXML_INCLUDE_DIR LIBMXML_LIBRARY) - diff --git a/cmake/FindNEXUS.cmake b/cmake/FindNEXUS.cmake deleted file mode 100644 index dca59b4a..00000000 --- a/cmake/FindNEXUS.cmake +++ /dev/null @@ -1,45 +0,0 @@ -# - Find NeXus library -# Find the native NEXUS includes and library -# This module defines -# NEXUS_INCLUDE_DIR, where to find NeXus.h, etc. -# NEXUS_LIBRARY, library to link against to use NEXUS -# NEXUS_FOUND, if false, do not try to use NEXUS. - -find_path(NEXUS_INCLUDE_DIR napi.h - HINTS "/usr/local/include" "/opt/nexus/include" "/usr/local/include/nexus" -) -# find position of napi.h from the end -string(FIND "${NEXUS_INCLUDE_DIR}" "/napi.h" pos REVERSE) -# truncate the string -string(SUBSTRING "${NEXUS_INCLUDE_DIR}" 0 ${pos} substr) -set(NEXUS_INCLUDE_DIR ${substr}) -unset(substr) - -find_library(NEXUS_LIBRARY NeXus - HINTS "/usr/lib" "/usr/lib64" "/usr/local/lib" "/usr/local/lib64" "/opt/nexus/lib") - -# get version string -if (NEXUS_INCLUDE_DIR AND EXISTS ${NEXUS_INCLUDE_DIR}/napi.h) - file(STRINGS "${NEXUS_INCLUDE_DIR}/napi.h" NEXUS_version_str - REGEX "^#define[\t ]+NEXUS_VERSION[\t ].*") - - string(REGEX REPLACE "^#define[\t ]+NEXUS_VERSION[\t ]+\"([^\"]*).*" - "\\1" NEXUS_VERSION_STRING "${NEXUS_version_str}") - unset(NEXUS_version_str) -endif() - -# handle the QUIETLY and REQUIRED arguments and set NEXUS_FOUND to TRUE if -# all listed variables are TRUE -include(${CMAKE_ROOT}/Modules/FindPackageHandleStandardArgs.cmake) -FIND_PACKAGE_HANDLE_STANDARD_ARGS(NEXUS - REQUIRED_VARS NEXUS_LIBRARY NEXUS_INCLUDE_DIR - VERSION_VAR NEXUS_VERSION_STRING) - -if (NOT NEXUS_FOUND) - unset(NEXUS_LIBRARY) -endif() - -mark_as_advanced(NEXUS_INCLUDE_DIR NEXUS_LIBRARY) - - - diff --git a/src/external/nexus/CMakeLists.txt b/src/external/nexus/CMakeLists.txt index 2dc2cf1f..3e866987 100644 --- a/src/external/nexus/CMakeLists.txt +++ b/src/external/nexus/CMakeLists.txt @@ -28,7 +28,7 @@ target_include_directories( ) #--- add library dependencies ------------------------------------------------- -target_link_libraries(PNeXus ${NEXUS_LIBRARY}) +target_link_libraries(PNeXus ${HDF4_LIBRARIES} ${ROOT_LIBRARIES}) #--- install PNeXus solib ----------------------------------------------------- install(TARGETS PNeXus DESTINATION lib) diff --git a/src/external/nexus/PNeXus.cpp b/src/external/nexus/PNeXus.cpp index 672170f0..32c5c8c3 100644 --- a/src/external/nexus/PNeXus.cpp +++ b/src/external/nexus/PNeXus.cpp @@ -29,7 +29,13 @@ #include #include +#include +#include #include +#include + +#include "Minuit2/MnMinimize.h" +#include "Minuit2/FunctionMinimum.h" #include "PNeXus.h" @@ -68,3 +74,2190 @@ nxs::HDFType nxs::checkHDFType(const std::string& filename) { return nxs::HDFType::Unknown; } + +#ifdef HAVE_HDF4 +//============================================================================= +// PNeXusDeadTime Constructor +//============================================================================= +nxH4::PNeXusDeadTime::PNeXusDeadTime(const nxH4::PNeXus *nxs, bool debug) : fDebug(debug) +{ + std::map dataMap = nxs->GetDataMap(); + + int idfVersion = nxs->GetIdfVersion(); + if ((idfVersion != 1) && (idfVersion != 2)) { + std::stringstream err; + err << "**ERROR** Found unsupported IDF version '" << idfVersion << "'"; + std::cout << err.str() << std::endl; + fValid = false; + return; + } + + if (idfVersion == 1) { + // not yet implemented + } else { // idfVersion == 2 + // get counts + if (dataMap.find("/raw_data_1/detector_1/counts") != dataMap.end()) { + auto counts_data = std::any_cast>(dataMap["/raw_data_1/detector_1/counts"]); + const auto& counts = counts_data.GetData(); + fCounts = counts; + auto dims = counts_data.GetDimensions(); + if (fDebug) { + std::cout << " counts dimensions: " << dims[0] << " x " + << dims[1] << " x " << dims[2] << std::endl; + } + fDims = dims; + } else { + std::cout << "**ERROR** Couldn't find 'counts' dataset" << std::endl; + fValid = false; + return; + } + + // get time_resolution + if (dataMap.find("/raw_data_1/detector_1/time_resolution") != dataMap.end()) { + auto time_res_data = std::any_cast>(dataMap["/raw_data_1/detector_1/time_resolution"]); + const auto& time_res = time_res_data.GetData(); + if (time_res.size() > 0) { + fTimeResolution = time_res[0]; + if (fDebug) { + std::cout << " time resolution: " << fTimeResolution << " ps" << std::endl; + } + } + } else { + std::cout << "**ERROR** Couldn't find 'time_resolution' dataset" << std::endl; + fValid = false; + return; + } + + // get good_frames + if (dataMap.find("/raw_data_1/good_frames") != dataMap.end()) { + auto good_frames_data = std::any_cast>(dataMap["/raw_data_1/good_frames"]); + const auto& good_frames = good_frames_data.GetData(); + if (good_frames.size() > 0) { + fGoodFrames = good_frames[0]; + if (fDebug) { + std::cout << " good frames: " << fGoodFrames << std::endl; + } + } + } else { + std::cout << "**ERROR** Couldn't find 'good_frames' dataset" << std::endl; + fValid = false; + return; + } + + // get t0_bin, first_good_bin, last_good_bin + if (dataMap.find("/raw_data_1/detector_1/spectrum_index") != dataMap.end()) { + auto spec_idx_data = std::any_cast>(dataMap["/raw_data_1/detector_1/spectrum_index"]); + if (spec_idx_data.HasAttribute("t0_bin")) { + try { + fT0Bin = std::any_cast(spec_idx_data.GetAttribute("t0_bin")); + } catch (...) { + std::cout << "**WARNING** Couldn't read t0_bin attribute" << std::endl; + } + } + + if (spec_idx_data.HasAttribute("first_good_bin")) { + try { + fFgbBin = std::any_cast(spec_idx_data.GetAttribute("first_good_bin")); + } catch (...) { + std::cout << "**WARNING** Couldn't read first_good_bin attribute" << std::endl; + } + } + + if (spec_idx_data.HasAttribute("last_good_bin")) { + try { + fLgbBin = std::any_cast(spec_idx_data.GetAttribute("last_good_bin")); + } catch (...) { + std::cout << "**WARNING** Couldn't read last_good_bin attribute" << std::endl; + } + } + } + + if (fDebug) { + std::cout << " t0_bin: " << fT0Bin << ", first_good_bin: " << fFgbBin + << ", last_good_bin: " << fLgbBin << std::endl; + } + } + + fValid = true; +} + +//============================================================================= +// PNeXusDeadTime::operator() +//============================================================================= +double nxH4::PNeXusDeadTime::operator()(const std::vector &par) const +{ + double chisq = 0.0; + double tau = par[0]; // dead time in microseconds + + // Convert time resolution from ps to microseconds + double dt = fTimeResolution * 1.0e-6; // ps -> us + + // Calculate chi-square for the specified spectrum + int period = 0; // assuming single period + + for (unsigned int bin = fFgbBin; bin <= static_cast(fLgbBin); bin++) { + int idx = period * fDims[1] * fDims[2] + fIdx * fDims[2] + bin; + double n_obs = static_cast(fCounts[idx]); + + if (n_obs > 0) { + // Dead time correction: n_true = n_obs / (1 - n_obs * tau / (good_frames * dt)) + double correction = 1.0 - (n_obs * tau) / (fGoodFrames * dt); + if (correction > 0) { + double n_expected = n_obs / correction; + double diff = n_obs - n_expected; + chisq += (diff * diff) / n_expected; + } + } + } + + return chisq; +} + +//============================================================================= +// PNeXusDeadTime::minimize +//============================================================================= +void nxH4::PNeXusDeadTime::minimize(const int i) +{ + if (i < 0 || i >= static_cast(fDims[1])) { + std::cerr << "**ERROR** Invalid spectrum index: " << i << std::endl; + return; + } + + fIdx = static_cast(i); + + ROOT::Minuit2::MnUserParameters upar; + upar.Add("dead_time", 0.1, 0.01); // initial value, step size + upar.SetLimits("dead_time", 0.0, 10.0); // limits in microseconds + + ROOT::Minuit2::MnMinimize minimize(*this, upar); + ROOT::Minuit2::FunctionMinimum min = minimize(); + + std::cout << "Spectrum " << i << ": "; + if (min.IsValid()) { + double dt = min.UserState().Value("dead_time"); + double dt_err = min.UserState().Error("dead_time"); + std::cout << "Dead time = " << dt << " +/- " << dt_err << " us" << std::endl; + } else { + std::cout << "Minimization failed" << std::endl; + } +} + +//============================================================================= +// PNeXus Default Constructor +//============================================================================= +nxH4::PNeXus::PNeXus() : fSdId(-1), fFileId(-1) +{ +} + +//============================================================================= +// PNeXus Constructor with filename +//============================================================================= +nxH4::PNeXus::PNeXus(const std::string fln, const bool printDebug) + : fPrintDebug(printDebug), fFileName(fln), fSdId(-1), fFileId(-1) +{ + if (ReadNexusFile() != 0) { + throw std::runtime_error("Failed to read NeXus file: " + fln); + } +} + +//============================================================================= +// PNeXus Destructor +//============================================================================= +nxH4::PNeXus::~PNeXus() +{ + if (fSdId != -1) { + SDend(fSdId); + } + if (fFileId != -1) { + Hclose(fFileId); + } +} + +//============================================================================= +// PNeXus::ReadNexusFile +//============================================================================= +int nxH4::PNeXus::ReadNexusFile() +{ + // Open the HDF4 file + fSdId = SDstart(fFileName.c_str(), DFACC_READ); + if (fSdId == FAIL) { + std::cerr << "**ERROR** Failed to open HDF4 file: " << fFileName << std::endl; + return 1; + } + + fFileId = Hopen(fFileName.c_str(), DFACC_READ, 0); + if (fFileId == FAIL) { + std::cerr << "**ERROR** Failed to open HDF4 file (H interface): " << fFileName << std::endl; + SDend(fSdId); + fSdId = -1; + return 1; + } + + // Read file-level attributes + int32 n_datasets, n_file_attrs; + if (SDfileinfo(fSdId, &n_datasets, &n_file_attrs) == FAIL) { + std::cerr << "**ERROR** Failed to get file info" << std::endl; + return 1; + } + + if (fPrintDebug) { + std::cout << "Number of datasets: " << n_datasets << std::endl; + std::cout << "Number of file attributes: " << n_file_attrs << std::endl; + } + + // Read HDF4 version and NeXus version from attributes + char attr_name[H4_MAX_NC_NAME]; + int32 attr_type, attr_count; + + for (int32 i = 0; i < n_file_attrs; i++) { + if (SDattrinfo(fSdId, i, attr_name, &attr_type, &attr_count) == FAIL) { + continue; + } + + if (caseInsensitiveEquals(attr_name, "HDF_version") || + caseInsensitiveEquals(attr_name, "HDF4_version")) { + std::vector buffer(attr_count + 1, '\0'); + if (SDreadattr(fSdId, i, buffer.data()) != FAIL) { + fHdf4Version = std::string(buffer.data()); + } + } else if (caseInsensitiveEquals(attr_name, "NeXus_version")) { + std::vector buffer(attr_count + 1, '\0'); + if (SDreadattr(fSdId, i, buffer.data()) != FAIL) { + fNeXusVersion = std::string(buffer.data()); + } + } else if (caseInsensitiveEquals(attr_name, "file_name")) { + std::vector buffer(attr_count + 1, '\0'); + if (SDreadattr(fSdId, i, buffer.data()) != FAIL) { + fFileNameNxs = std::string(buffer.data()); + } + } else if (caseInsensitiveEquals(attr_name, "file_time")) { + std::vector buffer(attr_count + 1, '\0'); + if (SDreadattr(fSdId, i, buffer.data()) != FAIL) { + fFileTimeNxs = std::string(buffer.data()); + } + } + } + + // Determine IDF version + // Try to find IDF_version dataset - check multiple possible locations + + // First try simple root-level name + int32 idf_idx = SDnametoindex(fSdId, "IDF_version"); + if (idf_idx != FAIL) { + int32 sds_id = SDselect(fSdId, idf_idx); + if (sds_id != FAIL) { + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs; + char name[H4_MAX_NC_NAME]; + if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) != FAIL) { + if (data_type == DFNT_INT32) { + int32 idf_val; + int32 start[H4_MAX_VAR_DIMS] = {0}; + int32 edges[H4_MAX_VAR_DIMS] = {1}; + if (SDreaddata(sds_id, start, nullptr, edges, &idf_val) != FAIL) { + fIdfVersion = idf_val; + } + } + } + SDendaccess(sds_id); + } + } + + // If not found, try /run/IDF_version (IDF version 1 location) + if (fIdfVersion == -1) { + try { + int32 sds_ref = findDatasetRefByPath("/run/IDF_version"); + if (sds_ref != -1) { + int32 sds_idx = SDreftoindex(fSdId, sds_ref); + if (sds_idx != FAIL) { + int32 sds_id = SDselect(fSdId, sds_idx); + if (sds_id != FAIL) { + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs; + char name[H4_MAX_NC_NAME]; + if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) != FAIL) { + if (data_type == DFNT_INT32) { + int32 idf_val; + int32 start[H4_MAX_VAR_DIMS] = {0}; + int32 edges[H4_MAX_VAR_DIMS] = {1}; + if (SDreaddata(sds_id, start, nullptr, edges, &idf_val) != FAIL) { + fIdfVersion = idf_val; + } + } + } + SDendaccess(sds_id); + } + } + } + } catch (...) { + // Path not found, continue to next attempt + } + } + + // If still not found, try /raw_data_1/idf_version (IDF version 2 location) + if (fIdfVersion == -1) { + try { + int32 sds_ref = findDatasetRefByPath("/raw_data_1/idf_version"); + if (sds_ref != -1) { + int32 sds_idx = SDreftoindex(fSdId, sds_ref); + if (sds_idx != FAIL) { + int32 sds_id = SDselect(fSdId, sds_idx); + if (sds_id != FAIL) { + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs; + char name[H4_MAX_NC_NAME]; + if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) != FAIL) { + if (data_type == DFNT_INT32) { + int32 idf_val; + int32 start[H4_MAX_VAR_DIMS] = {0}; + int32 edges[H4_MAX_VAR_DIMS] = {1}; + if (SDreaddata(sds_id, start, nullptr, edges, &idf_val) != FAIL) { + fIdfVersion = idf_val; + } + } + } + SDendaccess(sds_id); + } + } + } + } catch (...) { + // Path not found, will use default + } + } + + if (fIdfVersion == -1) { + // Try alternate approach - check for raw_data_1 structure + // If we find raw_data_1/detector_1, assume IDF v2 + // Otherwise assume IDF v1 + fIdfVersion = 2; // default to version 2 + } + + if (fPrintDebug) { + std::cout << "HDF4 Version: " << fHdf4Version << std::endl; + std::cout << "NeXus Version: " << fNeXusVersion << std::endl; + std::cout << "IDF Version: " << fIdfVersion << std::endl; + } + + // Read datasets based on IDF version + if (fIdfVersion == 1) { + HandleIdfV1(fSdId); + } else if (fIdfVersion == 2) { + HandleIdfV2(fSdId); + } else { + std::cerr << "**ERROR** Unsupported IDF version: " << fIdfVersion << std::endl; + return 1; + } + + return 0; +} + +//============================================================================= +// nxH4::PNeXus::HandleIdfV1 +//============================================================================= +void nxH4::PNeXus::HandleIdfV1(int32 sd_id) +{ + // IDF version 1 handling - not yet implemented + if (fPrintDebug) { + std::cout << "IDF version 1 handling ..." << std::endl; + } + + std::vector int_datasets = { + "/run/IDF_version", + "/run/number", + "/run/switching_state", + "/run/instrument/detector/number", + "/run/histogram_data_1/counts", + "/run/histogram_data_1/resolution", + "/run/histogram_data_1/grouping" + }; + + std::vector float_datasets = { + "/run/sample/temperature", + "/run/sample/magnetic_field", + "/run/sample/magnetic_field_vector", + "/run/instrument/detector/deadtimes", + "/run/histogram_data_1/time_zero", + "/run/histogram_data_1/raw_time", + "/run/histogram_data_1/corrected_time", + "/run/histogram_data_1/alpha" + }; + + std::vector string_datasets = { + "/run/program_name", + "/run/title", + "/run/notes", + "/run/analysis", + "/run/lab", + "/run/beamline", + "/run/start_time", + "/run/stop_time", + "/run/user/name", + "/run/user/experiment_number", + "/run/sample/name", + "/run/sample/magnetic_field_state", + "/run/sample/environment", + "/run/instrument/name", + "/run/instrument/collimator/type" + }; + + // Read integer datasets + for (const auto& path : int_datasets) { + try { + ReadIntDataset(sd_id, path); + } catch (const std::exception& e) { + if (fPrintDebug) { + std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl; + } + } + } + + // Read float datasets + for (const auto& path : float_datasets) { + try { + ReadFloatDataset(sd_id, path); + } catch (const std::exception& e) { + if (fPrintDebug) { + std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl; + } + } + } + + // Read string datasets + for (const auto& path : string_datasets) { + try { + ReadStringDataset(sd_id, path); + } catch (const std::exception& e) { + if (fPrintDebug) { + std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl; + } + } + } +} + +//============================================================================= +// nxH4::PNeXus::HandleIdfV2 +//============================================================================= +void nxH4::PNeXus::HandleIdfV2(int32 sd_id) +{ + // Read common datasets for IDF version 2 + // These are typical datasets found in muon NeXus files + + std::vector int_datasets = { + "/raw_data_1/detector_1/counts", + "/raw_data_1/detector_1/spectrum_index", + "/raw_data_1/good_frames", + "/raw_data_1/detector_1/raw_time", + "/raw_data_1/run_number" + }; + + std::vector float_datasets = { + "/raw_data_1/detector_1/time_resolution", + "/raw_data_1/detector_1/time_zero" + }; + + std::vector string_datasets = { + "/raw_data_1/name", + "/raw_data_1/title", + "/raw_data_1/start_time", + "/raw_data_1/end_time" + }; + + // Read integer datasets + for (const auto& path : int_datasets) { + try { + ReadIntDataset(sd_id, path); + } catch (const std::exception& e) { + if (fPrintDebug) { + std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl; + } + } + } + + // Read float datasets + for (const auto& path : float_datasets) { + try { + ReadFloatDataset(sd_id, path); + } catch (const std::exception& e) { + if (fPrintDebug) { + std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl; + } + } + } + + // Read string datasets + for (const auto& path : string_datasets) { + try { + ReadStringDataset(sd_id, path); + } catch (const std::exception& e) { + if (fPrintDebug) { + std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl; + } + } + } +} + +//============================================================================= +// nxH4::PNeXus::ReadIntDataset +//============================================================================= +void nxH4::PNeXus::ReadIntDataset(int32 sd_id, const std::string& path) +{ + // Extract dataset name from path + std::vector components = splitPath(path); + if (components.empty()) { + throw std::runtime_error("Invalid path: " + path); + } + std::string dataset_name = components.back(); + + // Use VGroup navigation to find the correct dataset reference + int32 sds_ref = -1; + if (components.size() > 1) { + // Try to resolve via VGroup hierarchy + sds_ref = findDatasetRefByPath(path); + } + + int32 sds_idx = -1; + if (sds_ref != -1) { + // Find SDS index from reference + sds_idx = SDreftoindex(sd_id, sds_ref); + } + + // Fallback to name-based search if VGroup navigation failed + if (sds_idx == FAIL || sds_idx == -1) { + sds_idx = findDatasetIndex(sd_id, dataset_name); + } + + int32 sds_id = SDselect(sd_id, sds_idx); + if (sds_id == FAIL) { + throw std::runtime_error("Failed to select dataset: " + dataset_name); + } + + // Get dataset info + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs; + char name[H4_MAX_NC_NAME]; + if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to get dataset info: " + dataset_name); + } + + // Calculate total number of elements + int32 n_elements = 1; + std::vector dimensions; + for (int32 i = 0; i < rank; i++) { + n_elements *= dim_sizes[i]; + dimensions.push_back(static_cast(dim_sizes[i])); + } + + // Read data + std::vector data(n_elements); + int32 start[H4_MAX_VAR_DIMS] = {0}; + if (SDreaddata(sds_id, start, nullptr, dim_sizes, data.data()) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to read dataset: " + dataset_name); + } + + // Convert to int vector + std::vector int_data(data.begin(), data.end()); + + // Create PNXdata object + PNXdata pnx_data(H4DataType::INT32); + pnx_data.SetData(int_data); + pnx_data.SetDimensions(dimensions); + + // Read attributes + ReadDatasetAttributes(sds_id, pnx_data); + + // Store in data map + fDataMap[path] = pnx_data; + + SDendaccess(sds_id); +} + +//============================================================================= +// nxH4::PNeXus::ReadFloatDataset +//============================================================================= +void nxH4::PNeXus::ReadFloatDataset(int32 sd_id, const std::string& path) +{ + // Extract dataset name from path + std::vector components = splitPath(path); + if (components.empty()) { + throw std::runtime_error("Invalid path: " + path); + } + std::string dataset_name = components.back(); + + // Use VGroup navigation to find the correct dataset reference + int32 sds_ref = -1; + if (components.size() > 1) { + // Try to resolve via VGroup hierarchy + sds_ref = findDatasetRefByPath(path); + } + + int32 sds_idx = -1; + if (sds_ref != -1) { + // Find SDS index from reference + sds_idx = SDreftoindex(sd_id, sds_ref); + } + + // Fallback to name-based search if VGroup navigation failed + if (sds_idx == FAIL || sds_idx == -1) { + sds_idx = findDatasetIndex(sd_id, dataset_name); + } + + int32 sds_id = SDselect(sd_id, sds_idx); + if (sds_id == FAIL) { + throw std::runtime_error("Failed to select dataset: " + dataset_name); + } + + // Get dataset info + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs; + char name[H4_MAX_NC_NAME]; + if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to get dataset info: " + dataset_name); + } + + // Calculate total number of elements + int32 n_elements = 1; + std::vector dimensions; + for (int32 i = 0; i < rank; i++) { + n_elements *= dim_sizes[i]; + dimensions.push_back(static_cast(dim_sizes[i])); + } + + // Read data + std::vector data(n_elements); + int32 start[H4_MAX_VAR_DIMS] = {0}; + if (SDreaddata(sds_id, start, nullptr, dim_sizes, data.data()) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to read dataset: " + dataset_name); + } + + // Convert to float vector + std::vector float_data(data.begin(), data.end()); + + // Create PNXdata object + PNXdata pnx_data(H4DataType::FLOAT32); + pnx_data.SetData(float_data); + pnx_data.SetDimensions(dimensions); + + // Read attributes + ReadDatasetAttributes(sds_id, pnx_data); + + // Store in data map + fDataMap[path] = pnx_data; + + SDendaccess(sds_id); +} + +//============================================================================= +// nxH4::PNeXus::ReadStringDataset +//============================================================================= +void nxH4::PNeXus::ReadStringDataset(int32 sd_id, const std::string& path) +{ + // Extract dataset name from path + std::vector components = splitPath(path); + if (components.empty()) { + throw std::runtime_error("Invalid path: " + path); + } + std::string dataset_name = components.back(); + + // Use VGroup navigation to find the correct dataset reference + int32 sds_ref = -1; + if (components.size() > 1) { + // Try to resolve via VGroup hierarchy + sds_ref = findDatasetRefByPath(path); + } + + int32 sds_idx = -1; + if (sds_ref != -1) { + // Find SDS index from reference + sds_idx = SDreftoindex(sd_id, sds_ref); + } + + // Fallback to name-based search if VGroup navigation failed + if (sds_idx == FAIL || sds_idx == -1) { + sds_idx = findDatasetIndex(sd_id, dataset_name); + } + + int32 sds_id = SDselect(sd_id, sds_idx); + if (sds_id == FAIL) { + throw std::runtime_error("Failed to select dataset: " + dataset_name); + } + + // Get dataset info + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs; + char name[H4_MAX_NC_NAME]; + if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to get dataset info: " + dataset_name); + } + + // Calculate total size for string + int32 n_elements = 1; + std::vector dimensions; + for (int32 i = 0; i < rank; i++) { + n_elements *= dim_sizes[i]; + dimensions.push_back(static_cast(dim_sizes[i])); + } + + // Read data as char array + std::vector data(n_elements + 1, '\0'); + int32 start[H4_MAX_VAR_DIMS] = {0}; + if (SDreaddata(sds_id, start, nullptr, dim_sizes, data.data()) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to read dataset: " + dataset_name); + } + + // Create string + std::vector string_data; + string_data.push_back(std::string(data.data())); + + // Create PNXdata object + PNXdata pnx_data(H4DataType::CHAR8); + pnx_data.SetData(string_data); + pnx_data.SetDimensions(dimensions); + + // Read attributes + ReadDatasetAttributes(sds_id, pnx_data); + + // Store in data map + fDataMap[path] = pnx_data; + + SDendaccess(sds_id); +} + +//============================================================================= +// nxH4::PNeXus::ReadDatasetAttributes (template specializations) +//============================================================================= +template +void nxH4::PNeXus::ReadDatasetAttributes(int32 sds_id, PNXdata& data) +{ + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs; + char name[H4_MAX_NC_NAME]; + + if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) == FAIL) { + return; + } + + for (int32 i = 0; i < n_attrs; i++) { + char attr_name[H4_MAX_NC_NAME]; + int32 attr_type, attr_count; + + if (SDattrinfo(sds_id, i, attr_name, &attr_type, &attr_count) == FAIL) { + continue; + } + + // Read attribute based on type - store as primitives like h5nexus + if (attr_type == DFNT_INT32) { + if (attr_count == 1) { + // Single value - store as int directly + int32 value; + if (SDreadattr(sds_id, i, &value) != FAIL) { + data.AddAttribute(attr_name, static_cast(value)); + } + } else { + // Multiple values - store as vector + std::vector attr_data(attr_count); + if (SDreadattr(sds_id, i, attr_data.data()) != FAIL) { + std::vector int_data(attr_data.begin(), attr_data.end()); + data.AddAttribute(attr_name, int_data); + } + } + } else if (attr_type == DFNT_FLOAT32) { + if (attr_count == 1) { + // Single value - store as float directly + float32 value; + if (SDreadattr(sds_id, i, &value) != FAIL) { + data.AddAttribute(attr_name, static_cast(value)); + } + } else { + // Multiple values - store as vector + std::vector attr_data(attr_count); + if (SDreadattr(sds_id, i, attr_data.data()) != FAIL) { + std::vector float_data(attr_data.begin(), attr_data.end()); + data.AddAttribute(attr_name, float_data); + } + } + } else if (attr_type == DFNT_CHAR8 || attr_type == DFNT_UCHAR8) { + // String attributes - always store as string + std::vector attr_data(attr_count + 1, '\0'); + if (SDreadattr(sds_id, i, attr_data.data()) != FAIL) { + data.AddAttribute(attr_name, std::string(attr_data.data())); + } + } + } +} + +//============================================================================= +// nxH4::PNeXus::Dump +//============================================================================= +void nxH4::PNeXus::Dump() +{ + int first_good_bin{0}; + + if (fIdfVersion == 1) { + std::cout << std::endl; + std::cout << "========================================" << std::endl; + std::cout << "NeXus File Dump" << std::endl; + std::cout << "========================================" << std::endl; + std::cout << "Filename: " << fFileName << std::endl; + std::cout << "HDF4 Version: " << fHdf4Version << std::endl; + std::cout << "NeXus Version: " << fNeXusVersion << std::endl; + std::cout << "IDF Version: " << fIdfVersion << std::endl; + std::cout << "----------------------------------------" << std::endl; + + std::cout << std::endl << "++++"; + std::cout << std::endl << "run"; + std::cout << std::endl << "----"; + std::cout << std::endl << " IDF_version: " << fIdfVersion; + if (fDataMap.find("/run/program_name") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/program_name"]); + std::cout << std::endl << " program_name : " << str_data.GetData()[0]; + + // Check for attributes + if (str_data.HasAttribute("version")) { + try { + auto version = std::any_cast(str_data.GetAttribute("version")); + std::cout << " version : " << version; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast version attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast program_name data" << std::endl; + } + } + if (fDataMap.find("/run/number") != fDataMap.end()) { + try { + auto number = std::any_cast>(fDataMap["/run/number"]); + std::cout << std::endl << " number : " << number.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast number data" << std::endl; + } + } + if (fDataMap.find("/run/title") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/title"]); + std::cout << std::endl << " title : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast title data" << std::endl; + } + } + if (fDataMap.find("/run/notes") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/notes"]); + std::cout << std::endl << " notes : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast notes data" << std::endl; + } + } + if (fDataMap.find("/run/analysis") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/analysis"]); + std::cout << std::endl << " analysis : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast analysis data" << std::endl; + } + } + if (fDataMap.find("/run/lab") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/lab"]); + std::cout << std::endl << " lab : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast lab data" << std::endl; + } + } + if (fDataMap.find("/run/beamline") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/beamline"]); + std::cout << std::endl << " beamline : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast beamline data" << std::endl; + } + } + if (fDataMap.find("/run/start_time") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/start_time"]); + std::cout << std::endl << " start_time : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast start_time data" << std::endl; + } + } + if (fDataMap.find("/run/stop_time") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/stop_time"]); + std::cout << std::endl << " stop_time : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast stop_time data" << std::endl; + } + } + if (fDataMap.find("/run/switching_state") != fDataMap.end()) { + try { + auto int_data = std::any_cast>(fDataMap["/run/switching_state"]); + std::cout << std::endl << " switching_state : " << int_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast switching_state data" << std::endl; + } + } + std::cout << std::endl << "----"; + std::cout << std::endl << " user"; + std::cout << std::endl << "----"; + if (fDataMap.find("/run/user/name") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/user/name"]); + std::cout << std::endl << " name : '" << str_data.GetData()[0] << "'"; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast name data" << std::endl; + } + } + if (fDataMap.find("/run/user/experiment_number") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/user/experiment_number"]); + std::cout << std::endl << " experiment_number : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast experiment_number data" << std::endl; + } + } + std::cout << std::endl << "----"; + std::cout << std::endl << " sample"; + std::cout << std::endl << "----"; + if (fDataMap.find("/run/sample/name") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/sample/name"]); + std::cout << std::endl << " name : '" << str_data.GetData()[0] << "'"; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast name data" << std::endl; + } + } + if (fDataMap.find("/run/sample/temperature") != fDataMap.end()) { + try { + auto float_data = std::any_cast>(fDataMap["/run/sample/temperature"]); + std::cout << std::endl << " temperature : " << float_data.GetData()[0]; + + // Check for attributes + if (float_data.HasAttribute("units")) { + try { + auto units = std::any_cast(float_data.GetAttribute("units")); + std::cout << " units : " << units; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast temperature data" << std::endl; + } + } + if (fDataMap.find("/run/sample/magnetic_field") != fDataMap.end()) { + try { + auto float_data = std::any_cast>(fDataMap["/run/sample/magnetic_field"]); + std::cout << std::endl << " magnetic_field : " << float_data.GetData()[0]; + + // Check for attributes + if (float_data.HasAttribute("units")) { + try { + auto units = std::any_cast(float_data.GetAttribute("units")); + std::cout << " units : " << units; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast magnetic_field data" << std::endl; + } + } + if (fDataMap.find("/run/sample/magnetic_field_vector") != fDataMap.end()) { + try { + auto float_data = std::any_cast>(fDataMap["/run/sample/magnetic_field_vector"]); + std::cout << std::endl << " magnetic_field_vector : " << float_data.GetData()[0]; + + // Check for attributes + if (float_data.HasAttribute("coordinate_system")) { + try { + auto coordinate_system = std::any_cast(float_data.GetAttribute("coordinate_system")); + std::cout << std::endl << " coordinate_system : " << coordinate_system << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast coordinate_system attribute" << std::endl; + } + } + if (float_data.HasAttribute("units")) { + try { + auto units = std::any_cast(float_data.GetAttribute("units")); + std::cout << " units : " << units << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl; + } + } + if (float_data.HasAttribute("available")) { + try { + auto available = std::any_cast(float_data.GetAttribute("available")); + std::cout << " available : " << available; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast available attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast magnetic_field_vector data" << std::endl; + } + } + if (fDataMap.find("/run/sample/environment") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/sample/environment"]); + std::cout << std::endl << " environment : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast environment data" << std::endl; + } + } + std::cout << std::endl << "----"; + std::cout << std::endl << " instrument"; + std::cout << std::endl << "----"; + if (fDataMap.find("/run/instrument/name") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/instrument/name"]); + std::cout << std::endl << " name : '" << str_data.GetData()[0] << "'"; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast name data" << std::endl; + } + } + std::cout << std::endl << "----"; + std::cout << std::endl << " detector"; + std::cout << std::endl << "----"; + if (fDataMap.find("/run/instrument/detector/number") != fDataMap.end()) { + try { + auto int_data = std::any_cast>(fDataMap["/run/instrument/detector/number"]); + std::cout << std::endl << " number : " << int_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast number data" << std::endl; + } + } + if (fDataMap.find("/run/instrument/detector/deadtimes") != fDataMap.end()) { + try { + auto float_data = std::any_cast>(fDataMap["/run/instrument/detector/deadtimes"]); + auto data = float_data.GetData(); + unsigned int end{15}; + if (data.size() < end) + end = data.size(); + std::cout << std::endl << " deadtimes: "; + for (unsigned int i=0; i>(fDataMap["/run/instrument/collimator/type"]); + std::cout << std::endl << " type : " << str_data.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast type data" << std::endl; + } + } + + std::cout << std::endl << "----"; + std::cout << std::endl << " histogram_data_1"; + std::cout << std::endl << "----"; + std::cout << std::endl << " counts:"; + std::cout << std::endl; + try { + auto counts_data = std::any_cast>(fDataMap["/run/histogram_data_1/counts"]); + + // Check for attributes + if (counts_data.HasAttribute("units")) { + try { + auto units = std::any_cast(counts_data.GetAttribute("units")); + std::cout << " units : " << units << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl; + } + } + if (counts_data.HasAttribute("signal")) { + try { + auto signal = std::any_cast(counts_data.GetAttribute("signal")); + std::cout << " signal : " << signal << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast signal attribute" << std::endl; + } + } + int noOfHistos{0}; + if (counts_data.HasAttribute("number")) { + try { + noOfHistos = std::any_cast(counts_data.GetAttribute("number")); + std::cout << " number : " << noOfHistos << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast number attribute" << std::endl; + } + } + int histoLength{0}; + if (counts_data.HasAttribute("length")) { + try { + histoLength = std::any_cast(counts_data.GetAttribute("length")); + std::cout << " length : " << histoLength << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast length attribute" << std::endl; + } + } + if (counts_data.HasAttribute("t0_bin")) { + try { + auto t0_bin = std::any_cast(counts_data.GetAttribute("t0_bin")); + std::cout << " t0_bin : " << t0_bin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast t0_bin attribute" << std::endl; + } + } + if (counts_data.HasAttribute("first_good_bin")) { + try { + first_good_bin = std::any_cast(counts_data.GetAttribute("first_good_bin")); + std::cout << " first_good_bin : " << first_good_bin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast first_good_bin attribute" << std::endl; + } + } + if (counts_data.HasAttribute("last_good_bin")) { + try { + auto last_good_bin = std::any_cast(counts_data.GetAttribute("last_good_bin")); + std::cout << " last_good_bin : " << last_good_bin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast last_good_bin attribute" << std::endl; + } + } + if (counts_data.HasAttribute("offset")) { + try { + auto offset = std::any_cast(counts_data.GetAttribute("offset")); + std::cout << " offset : " << offset << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast offset attribute" << std::endl; + } + } + + // dump the first couple of counts of each detector + const auto& data = counts_data.GetData(); + std::cout << std::endl << " first couple of counts of each detector:"; + for (unsigned int i=0; i>(fDataMap["/run/histogram_data_1/resolution"]); + std::cout << std::endl << " resolution : " << int_data.GetData()[0]; + + if (int_data.HasAttribute("units")) { + try { + auto units = std::any_cast(int_data.GetAttribute("units")); + std::cout << " units : " << units; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast resolution data" << std::endl; + } + } + + if (fDataMap.find("/run/histogram_data_1/time_zero") != fDataMap.end()) { + try { + auto float_data = std::any_cast>(fDataMap["/run/histogram_data_1/time_zero"]); + std::cout << std::endl << " time_zero : " << float_data.GetData()[0]; + + if (float_data.HasAttribute("units")) { + try { + auto units = std::any_cast(float_data.GetAttribute("units")); + std::cout << std::endl << " units : " << units; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl; + } + } + if (float_data.HasAttribute("available")) { + try { + auto available = std::any_cast(float_data.GetAttribute("available")); + std::cout << std::endl << " available : " << available; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast available attribute"; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast time_zero data" << std::endl; + } + } + + if (fDataMap.find("/run/histogram_data_1/raw_time") != fDataMap.end()) { + try { + auto float_data = std::any_cast>(fDataMap["/run/histogram_data_1/raw_time"]); + std::cout << std::endl << " raw_time : " << float_data.GetData()[0]; + if (float_data.HasAttribute("axis")) { + try { + auto axis = std::any_cast(float_data.GetAttribute("axis")); + std::cout << std::endl << " axis : " << axis; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast axis attribute" << std::endl; + } + } + if (float_data.HasAttribute("primary")) { + try { + auto primary = std::any_cast(float_data.GetAttribute("primary")); + std::cout << std::endl << " primary : " << primary; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast primary attribute" << std::endl; + } + } + if (float_data.HasAttribute("units")) { + try { + auto units = std::any_cast(float_data.GetAttribute("units")); + std::cout << std::endl << " units : " << units; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast raw_time data" << std::endl; + } + } + if (fDataMap.find("/run/histogram_data_1/corrected_time") != fDataMap.end()) { + try { + auto float_data = std::any_cast>(fDataMap["/run/histogram_data_1/corrected_time"]); + std::cout << std::endl << " corrected_time : " << float_data.GetData()[0]; + if (float_data.HasAttribute("axis")) { + try { + auto axis = std::any_cast(float_data.GetAttribute("axis")); + std::cout << std::endl << " axis : " << axis; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast axis attribute" << std::endl; + } + } + if (float_data.HasAttribute("units")) { + try { + auto units = std::any_cast(float_data.GetAttribute("units")); + std::cout << std::endl << " units : " << units; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast raw_time data" << std::endl; + } + } + if (fDataMap.find("/run/histogram_data_1/grouping") != fDataMap.end()) { + try { + auto int_data = std::any_cast>(fDataMap["/run/histogram_data_1/grouping"]); + auto data = int_data.GetData(); + unsigned int end{15}; + if (data.size() < end) + end = data.size(); + std::cout << std::endl << " grouping : "; + for (unsigned int i=0; i(int_data.GetAttribute("available")); + std::cout << std::endl << " available : " << available; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast available attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast grouping data" << std::endl; + } + } + if (fDataMap.find("/run/histogram_data_1/alpha") != fDataMap.end()) { + try { + auto float_data = std::any_cast>(fDataMap["/run/histogram_data_1/alpha"]); + auto data = float_data.GetData(); + unsigned int end{15}; + if (data.size() < end) + end = data.size(); + std::cout << std::endl << " alpha : "; + for (unsigned int i=0; i(float_data.GetAttribute("available")); + std::cout << std::endl << " available : " << available; + } catch (const std::bad_any_cast& e) { + std::cerr << "**ERROR**: Failed to cast available attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "**ERROR**: Failed to cast alpha data" << std::endl; + } + } + std::cout << std::endl; + std::cout << "========================================" << std::endl; + std::cout << std::endl; + } else { // IDF Version 2 + std::cout << std::endl; + std::cout << std::endl << "hdf5-NeXus file content of file:' " << fFileName << "'"; + std::cout << std::endl << "****"; + std::cout << std::endl << "Top Level Attributes:"; + std::cout << std::endl << " HDF4 Version : " << fHdf4Version; + std::cout << std::endl << " NeXus Version: " << fNeXusVersion; + std::cout << std::endl << " file_name : " << fFileNameNxs; + std::cout << std::endl << " file_time : " << fFileTimeNxs; + std::cout << std::endl << "++++"; + std::cout << std::endl << "raw_data_1"; + std::cout << std::endl << "----"; + std::cout << std::endl << " IDF_version: " << fIdfVersion; + if (fDataMap.find("/raw_data_1/beamline") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/beamline"]); + std::cout << std::endl << " beamline : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast beamline data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/definition") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/definition"]); + std::cout << std::endl << " definition : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast definition data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/run_number") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/run_number"]); + std::cout << std::endl << " run_number : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast run_number data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/title") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/title"]); + std::cout << std::endl << " title : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast title data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/start_time") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/start_time"]); + std::cout << std::endl << " start_time : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast start_time data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/end_time") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/end_time"]); + std::cout << std::endl << " end_time : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast end_time data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/good_frames") != fDataMap.end()) { + try { + auto good_frames = std::any_cast>(fDataMap["/raw_data_1/good_frames"]); + std::cout << std::endl << " good_frames: " << good_frames.GetData()[0]; + } catch(const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast good_frames data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/experiment_identifier") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/experiment_identifier"]); + std::cout << std::endl << " experiment_identifier: " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast experiment_identifier data" << std::endl; + } + } + std::cout << std::endl << "----"; + std::cout << std::endl << " instrument"; + if (fDataMap.find("/raw_data_1/instrument/name") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/instrument/name"]); + std::cout << std::endl << " name : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast instrument/name data" << std::endl; + } + } + std::cout << std::endl << "----"; + std::cout << std::endl << " source"; + if (fDataMap.find("/raw_data_1/instrument/source/name") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/instrument/source/name"]); + std::cout << std::endl << " name : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast instrument/source/name data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/instrument/source/type") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/instrument/source/type"]); + std::cout << std::endl << " type : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast instrument/source/type data" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/instrument/source/probe") != fDataMap.end()) { + try { + auto str = std::any_cast>(fDataMap["/raw_data_1/instrument/source/probe"]); + std::cout << std::endl << " probe : " << str.GetData()[0]; + } catch (const std::bad_any_cast& e) { + std::cerr << std::endl << "Error: Failed to cast instrument/source/probe data" << std::endl; + } + } + std::cout << std::endl << "----"; + std::cout << std::endl << " detector_1"; + std::cout << std::endl << " counts:"; + std::cout << std::endl; + try { + auto counts_data = std::any_cast>(fDataMap["/raw_data_1/detector_1/counts"]); + auto dims = counts_data.GetDimensions(); + std::cout << " counts dimensions: " << dims[0] << " x " + << dims[1] << " x " << dims[2] << std::endl; + std::cout << " total elements: " << counts_data.GetNumElements() << std::endl; + + // Check for attributes + if (counts_data.HasAttribute("signal")) { + try { + auto signal = std::any_cast(counts_data.GetAttribute("signal")); + std::cout << " signal : " << signal << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast signal attribute" << std::endl; + } + } + if (counts_data.HasAttribute("axes")) { + try { + auto axes = std::any_cast(counts_data.GetAttribute("axes")); + std::cout << " axes : " << axes << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast axes attribute" << std::endl; + } + } + if (counts_data.HasAttribute("long_name")) { + try { + auto long_name = std::any_cast(counts_data.GetAttribute("long_name")); + std::cout << " long_name : " << long_name << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast long_name attribute" << std::endl; + } + } + if (counts_data.HasAttribute("t0_bin")) { + try { + auto t0_bin = std::any_cast(counts_data.GetAttribute("t0_bin")); + std::cout << " t0_bin : " << t0_bin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast t0_bin attribute" << std::endl; + } + } + if (counts_data.HasAttribute("first_good_bin")) { + try { + first_good_bin = std::any_cast(counts_data.GetAttribute("first_good_bin")); + std::cout << " first_good_bin : " << first_good_bin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast first_good_bin attribute" << std::endl; + } + } + if (counts_data.HasAttribute("last_good_bin")) { + try { + auto last_good_bin = std::any_cast(counts_data.GetAttribute("last_good_bin")); + std::cout << " last_good_bin : " << last_good_bin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast last_good_bin attribute" << std::endl; + } + } + if (fDataMap.find("/raw_data_1/instrument/detector_1/resolution") != fDataMap.end()) { + try { + auto ivalData = std::any_cast>(fDataMap["/raw_data_1/instrument/detector_1/resolution"]); + std::cout << " resolution : " << ivalData.GetData()[0]; + if (ivalData.HasAttribute("units")) { + try { + auto units = std::any_cast(ivalData.GetAttribute("units")); + std::cout << " " << units << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast units attribute" << std::endl; + } + } + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast resolution attribute" << std::endl; + } + } + + // dump the first couple of counts of each detector + const auto& data = counts_data.GetData(); + std::cout << std::endl << " first couple of counts of each detector:"; + for (unsigned int i=0; i>(fDataMap["/raw_data_1/detector_1/raw_time"]); + const auto& data = raw_time.GetData(); + + // Check for attributes + if (raw_time.HasAttribute("units")) { + try { + auto units = std::any_cast(raw_time.GetAttribute("units")); + std::cout << " units : " << units << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast units attribute" << std::endl; + } + } + + // dump the first couple of raw_times + std::cout << " "; + for (unsigned int i=0; i>(fDataMap["/raw_data_1/detector_1/spectrum_index"]); + const auto& data = spectrum_index.GetData(); + + // dump the first couple of raw_times + std::cout << " "; + for (unsigned int i=0; i>(fDataMap["/raw_data_1/detector_1/dead_time"]); + const auto& data = dead_time.GetData(); + + // dump the first couple of raw_times + std::cout << " "; + for (unsigned int i=0; i(attr_value); + SDsetattr(sd_id, attr_name.c_str(), DFNT_CHAR8, + str_attr.length(), str_attr.c_str()); + } catch (...) { + // Try other types if needed + } + } + } +} + +//============================================================================= +// nxH4::PNeXus::WriteIdfV2 +//============================================================================= +void nxH4::PNeXus::WriteIdfV2(int32 sd_id) +{ + // Write all datasets from data map + for (const auto& [path, data_any] : fDataMap) { + // Try to write as int dataset + try { + auto data = std::any_cast>(data_any); + WriteIntDataset(sd_id, path, data); + continue; + } catch (...) {} + + // Try to write as float dataset + try { + auto data = std::any_cast>(data_any); + WriteFloatDataset(sd_id, path, data); + continue; + } catch (...) {} + + // Try to write as string dataset + try { + auto data = std::any_cast>(data_any); + WriteStringDataset(sd_id, path, data); + continue; + } catch (...) {} + } +} + +//============================================================================= +// nxH4::PNeXus::WriteIntDataset +//============================================================================= +void nxH4::PNeXus::WriteIntDataset(int32 sd_id, const std::string& path, + const PNXdata& data) +{ + // Extract dataset name from path + std::vector components = splitPath(path); + if (components.empty()) { + throw std::runtime_error("Invalid path: " + path); + } + std::string dataset_name = components.back(); + + const auto& dims = data.GetDimensions(); + int32 rank = static_cast(dims.size()); + std::vector dim_sizes(rank); + for (int32 i = 0; i < rank; i++) { + dim_sizes[i] = static_cast(dims[i]); + } + + // Create dataset + int32 sds_id = SDcreate(sd_id, dataset_name.c_str(), DFNT_INT32, rank, dim_sizes.data()); + if (sds_id == FAIL) { + throw std::runtime_error("Failed to create dataset: " + dataset_name); + } + + // Convert data to int32 + const auto& int_data = data.GetData(); + std::vector data32(int_data.begin(), int_data.end()); + + // Write data + int32 start[H4_MAX_VAR_DIMS] = {0}; + if (SDwritedata(sds_id, start, nullptr, dim_sizes.data(), data32.data()) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to write dataset: " + dataset_name); + } + + // Write attributes + WriteDatasetAttributes(sds_id, data); + + SDendaccess(sds_id); +} + +//============================================================================= +// nxH4::PNeXus::WriteFloatDataset +//============================================================================= +void nxH4::PNeXus::WriteFloatDataset(int32 sd_id, const std::string& path, + const PNXdata& data) +{ + // Extract dataset name from path + std::vector components = splitPath(path); + if (components.empty()) { + throw std::runtime_error("Invalid path: " + path); + } + std::string dataset_name = components.back(); + + const auto& dims = data.GetDimensions(); + int32 rank = static_cast(dims.size()); + std::vector dim_sizes(rank); + for (int32 i = 0; i < rank; i++) { + dim_sizes[i] = static_cast(dims[i]); + } + + // Create dataset + int32 sds_id = SDcreate(sd_id, dataset_name.c_str(), DFNT_FLOAT32, rank, dim_sizes.data()); + if (sds_id == FAIL) { + throw std::runtime_error("Failed to create dataset: " + dataset_name); + } + + // Convert data to float32 + const auto& float_data = data.GetData(); + std::vector data32(float_data.begin(), float_data.end()); + + // Write data + int32 start[H4_MAX_VAR_DIMS] = {0}; + if (SDwritedata(sds_id, start, nullptr, dim_sizes.data(), data32.data()) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to write dataset: " + dataset_name); + } + + // Write attributes + WriteDatasetAttributes(sds_id, data); + + SDendaccess(sds_id); +} + +//============================================================================= +// nxH4::PNeXus::WriteStringDataset +//============================================================================= +void nxH4::PNeXus::WriteStringDataset(int32 sd_id, const std::string& path, + const PNXdata& data) +{ + // Extract dataset name from path + std::vector components = splitPath(path); + if (components.empty()) { + throw std::runtime_error("Invalid path: " + path); + } + std::string dataset_name = components.back(); + + const auto& string_data = data.GetData(); + if (string_data.empty()) { + return; + } + + std::string str = string_data[0]; + int32 dims[1] = {static_cast(str.length())}; + + // Create dataset + int32 sds_id = SDcreate(sd_id, dataset_name.c_str(), DFNT_CHAR8, 1, dims); + if (sds_id == FAIL) { + throw std::runtime_error("Failed to create dataset: " + dataset_name); + } + + // Write data - HDF4 requires non-const pointer + std::vector str_buffer(str.begin(), str.end()); + int32 start[1] = {0}; + if (SDwritedata(sds_id, start, nullptr, dims, str_buffer.data()) == FAIL) { + SDendaccess(sds_id); + throw std::runtime_error("Failed to write dataset: " + dataset_name); + } + + // Write attributes + WriteDatasetAttributes(sds_id, data); + + SDendaccess(sds_id); +} + +//============================================================================= +// nxH4::PNeXus::WriteDatasetAttributes (template specializations) +//============================================================================= +template +void nxH4::PNeXus::WriteDatasetAttributes(int32 sds_id, const PNXdata& data) +{ + const auto& attributes = data.GetAttributes(); + + for (const auto& [attr_name, attr_value] : attributes) { + // Try different attribute types - primitives first (like h5nexus) + + // Try scalar int + try { + int32 value = std::any_cast(attr_value); + SDsetattr(sds_id, attr_name.c_str(), DFNT_INT32, 1, &value); + continue; + } catch (...) {} + + // Try scalar float + try { + float32 value = std::any_cast(attr_value); + SDsetattr(sds_id, attr_name.c_str(), DFNT_FLOAT32, 1, &value); + continue; + } catch (...) {} + + // Try string + try { + std::string str = std::any_cast(attr_value); + SDsetattr(sds_id, attr_name.c_str(), DFNT_CHAR8, str.length(), str.c_str()); + continue; + } catch (...) {} + + // Try vector for multi-element attributes + try { + std::vector int_vec = std::any_cast>(attr_value); + std::vector data32(int_vec.begin(), int_vec.end()); + SDsetattr(sds_id, attr_name.c_str(), DFNT_INT32, data32.size(), data32.data()); + continue; + } catch (...) {} + + // Try vector for multi-element attributes + try { + std::vector float_vec = std::any_cast>(attr_value); + std::vector data32(float_vec.begin(), float_vec.end()); + SDsetattr(sds_id, attr_name.c_str(), DFNT_FLOAT32, data32.size(), data32.data()); + continue; + } catch (...) {} + } +} + +//============================================================================= +// Group attribute methods +//============================================================================= +bool nxH4::PNeXus::AddGroupAttribute(const std::string& groupPath, + const std::string& attrName, + const std::any& attrValue) +{ + fGroupAttributes[groupPath][attrName] = attrValue; + return true; +} + +bool nxH4::PNeXus::RemoveGroupAttribute(const std::string& groupPath, + const std::string& attrName) +{ + auto it = fGroupAttributes.find(groupPath); + if (it != fGroupAttributes.end()) { + auto attr_it = it->second.find(attrName); + if (attr_it != it->second.end()) { + it->second.erase(attr_it); + return true; + } + } + return false; +} + +bool nxH4::PNeXus::HasGroupAttribute(const std::string& groupPath, + const std::string& attrName) const +{ + auto it = fGroupAttributes.find(groupPath); + if (it != fGroupAttributes.end()) { + return it->second.find(attrName) != it->second.end(); + } + return false; +} + +std::any nxH4::PNeXus::GetGroupAttribute(const std::string& groupPath, + const std::string& attrName) const +{ + return fGroupAttributes.at(groupPath).at(attrName); +} + +const std::map& nxH4::PNeXus::GetGroupAttributes( + const std::string& groupPath) const +{ + static std::map empty_map; + auto it = fGroupAttributes.find(groupPath); + if (it != fGroupAttributes.end()) { + return it->second; + } + return empty_map; +} + +bool nxH4::PNeXus::ClearGroupAttributes(const std::string& groupPath) +{ + auto it = fGroupAttributes.find(groupPath); + if (it != fGroupAttributes.end()) { + it->second.clear(); + return true; + } + return false; +} + +bool nxH4::PNeXus::AddRootAttribute(const std::string& attrName, const std::any& attrValue) +{ + return AddGroupAttribute("/", attrName, attrValue); +} + +//============================================================================= +// Helper methods +//============================================================================= +bool nxH4::PNeXus::caseInsensitiveEquals(const std::string& a, const std::string& b) +{ + if (a.length() != b.length()) { + return false; + } + return std::equal(a.begin(), a.end(), b.begin(), + [](char a, char b) { + return std::tolower(static_cast(a)) == + std::tolower(static_cast(b)); + }); +} + +std::vector nxH4::PNeXus::splitPath(const std::string& path) +{ + std::vector components; + std::string current; + + for (char c : path) { + if (c == '/') { + if (!current.empty()) { + components.push_back(current); + current.clear(); + } + } else { + current += c; + } + } + + if (!current.empty()) { + components.push_back(current); + } + + return components; +} + +std::string nxH4::PNeXus::findAttributeName(int32 sd_id, const std::string& requestedName) +{ + int32 n_datasets, n_attrs; + if (SDfileinfo(sd_id, &n_datasets, &n_attrs) == FAIL) { + throw std::runtime_error("Failed to get file info"); + } + + char attr_name[H4_MAX_NC_NAME]; + int32 attr_type, attr_count; + + for (int32 i = 0; i < n_attrs; i++) { + if (SDattrinfo(sd_id, i, attr_name, &attr_type, &attr_count) != FAIL) { + if (caseInsensitiveEquals(attr_name, requestedName)) { + return std::string(attr_name); + } + } + } + + throw std::runtime_error("Attribute not found: " + requestedName); +} + +int32 nxH4::PNeXus::findDatasetIndex(int32 sd_id, const std::string& requestedName) +{ + // First try exact match + int32 idx = SDnametoindex(sd_id, requestedName.c_str()); + if (idx != FAIL) { + return idx; + } + + // Try case-insensitive search + int32 n_datasets, n_attrs; + if (SDfileinfo(sd_id, &n_datasets, &n_attrs) == FAIL) { + throw std::runtime_error("Failed to get file info"); + } + + for (int32 i = 0; i < n_datasets; i++) { + int32 sds_id = SDselect(sd_id, i); + if (sds_id != FAIL) { + char name[H4_MAX_NC_NAME]; + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_ds_attrs; + if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_ds_attrs) != FAIL) { + if (caseInsensitiveEquals(name, requestedName)) { + SDendaccess(sds_id); + return i; + } + } + SDendaccess(sds_id); + } + } + + throw std::runtime_error("Dataset not found: " + requestedName); +} + +int32 nxH4::PNeXus::findDatasetRefByPath(const std::string& path) +{ + // Navigate VGroup hierarchy to find the dataset reference + // Path format: /group1/group2/.../dataset_name + + std::vector components = splitPath(path); + if (components.empty()) { + return -1; + } + + // Open file with V interface + if (Vstart(fFileId) == FAIL) { + return -1; + } + + int32 result_ref = -1; + int32 current_vg_ref = -1; + + // Find root VGroup (first component after /) + // For IDF version 1, root is typically "run" + if (components.size() >= 1) { + // Get all lone vgroups + int32 n_vgroups = Vlone(fFileId, nullptr, 0); + if (n_vgroups > 0) { + std::vector vgroup_refs(n_vgroups); + Vlone(fFileId, vgroup_refs.data(), n_vgroups); + + // Find the root VGroup + for (int32 i = 0; i < n_vgroups; i++) { + int32 vg_id = Vattach(fFileId, vgroup_refs[i], "r"); + if (vg_id != FAIL) { + char vg_name[VGNAMELENMAX]; + if (Vgetname(vg_id, vg_name) != FAIL) { + if (caseInsensitiveEquals(vg_name, components[0])) { + current_vg_ref = vgroup_refs[i]; + Vdetach(vg_id); + break; + } + } + Vdetach(vg_id); + } + } + } + } + + // Navigate through intermediate VGroups + for (size_t comp_idx = 1; comp_idx < components.size() - 1; comp_idx++) { + if (current_vg_ref == -1) break; + + int32 vg_id = Vattach(fFileId, current_vg_ref, "r"); + if (vg_id == FAIL) { + current_vg_ref = -1; + break; + } + + int32 n_entries = Vntagrefs(vg_id); + bool found = false; + + for (int32 i = 0; i < n_entries; i++) { + int32 tag, ref; + if (Vgettagref(vg_id, i, &tag, &ref) != FAIL) { + // Check if it's a VGroup + if (tag == DFTAG_VG) { + int32 child_vg_id = Vattach(fFileId, ref, "r"); + if (child_vg_id != FAIL) { + char child_name[VGNAMELENMAX]; + if (Vgetname(child_vg_id, child_name) != FAIL) { + if (caseInsensitiveEquals(child_name, components[comp_idx])) { + current_vg_ref = ref; + found = true; + Vdetach(child_vg_id); + break; + } + } + Vdetach(child_vg_id); + } + } + } + } + + Vdetach(vg_id); + + if (!found) { + current_vg_ref = -1; + break; + } + } + + // Find the dataset in the final VGroup + if (current_vg_ref != -1 && components.size() >= 1) { + int32 vg_id = Vattach(fFileId, current_vg_ref, "r"); + if (vg_id != FAIL) { + int32 n_entries = Vntagrefs(vg_id); + std::string target_name = components.back(); + + for (int32 i = 0; i < n_entries; i++) { + int32 tag, ref; + if (Vgettagref(vg_id, i, &tag, &ref) != FAIL) { + // Check if it's a Numeric Data (SDS) + if (tag == DFTAG_NDG) { + // Try to get the name by converting ref to index + int32 sds_idx = SDreftoindex(fSdId, ref); + if (sds_idx != FAIL) { + int32 sds_id = SDselect(fSdId, sds_idx); + if (sds_id != FAIL) { + char ds_name[H4_MAX_NC_NAME]; + int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs; + if (SDgetinfo(sds_id, ds_name, &rank, dim_sizes, &data_type, &n_attrs) != FAIL) { + if (caseInsensitiveEquals(ds_name, target_name)) { + result_ref = ref; + SDendaccess(sds_id); + break; + } + } + SDendaccess(sds_id); + } + } + } + } + } + + Vdetach(vg_id); + } + } + + Vend(fFileId); + return result_ref; +} + +nxH4::H4DataType nxH4::PNeXus::convertHdf4Type(int32 hdf4_type) +{ + switch (hdf4_type) { + case DFNT_INT32: return H4DataType::INT32; + case DFNT_FLOAT32: return H4DataType::FLOAT32; + case DFNT_FLOAT64: return H4DataType::FLOAT64; + case DFNT_CHAR8: return H4DataType::CHAR8; + case DFNT_UINT32: return H4DataType::UINT32; + case DFNT_INT16: return H4DataType::INT16; + case DFNT_UINT16: return H4DataType::UINT16; + case DFNT_INT8: return H4DataType::INT8; + case DFNT_UINT8: return H4DataType::UINT8; + default: return H4DataType::INT32; + } +} + +int32 nxH4::PNeXus::convertToHdf4Type(H4DataType dataType) +{ + switch (dataType) { + case H4DataType::INT32: return DFNT_INT32; + case H4DataType::FLOAT32: return DFNT_FLOAT32; + case H4DataType::FLOAT64: return DFNT_FLOAT64; + case H4DataType::CHAR8: return DFNT_CHAR8; + case H4DataType::UINT32: return DFNT_UINT32; + case H4DataType::INT16: return DFNT_INT16; + case H4DataType::UINT16: return DFNT_UINT16; + case H4DataType::INT8: return DFNT_INT8; + case H4DataType::UINT8: return DFNT_UINT8; + default: return DFNT_INT32; + } +} +#endif // HAVE_HDF4 diff --git a/src/external/nexus/PNeXus.h b/src/external/nexus/PNeXus.h index c855900a..b3212da1 100644 --- a/src/external/nexus/PNeXus.h +++ b/src/external/nexus/PNeXus.h @@ -31,6 +31,19 @@ #define _PNEXUS_H_ #include +#include +#include +#include +#include +#include +#include + +#include "Minuit2/FCNBase.h" + +#ifdef HAVE_HDF4 +#include +#include +#endif // HAVE_HDF4 namespace nxs { @@ -44,4 +57,839 @@ HDFType checkHDFType(const std::string& filename); } // end namespace nxs +#ifdef HAVE_HDF4 +namespace nxH4 { + +class PNeXus; + +/** + * @class PNeXusDeadTime + * @brief Dead time correction calculator for muon detector data + * + * The PNeXusDeadTime class calculates dead time corrections for muon detector + * count data using ROOT Minuit2 minimization. It inherits from ROOT::Minuit2::FCNBase + * to implement a chi-square minimization function. + * + * Dead time is the period after detecting an event during which the detector + * cannot register another event. This class estimates the dead time parameter + * for each detector spectrum by minimizing the deviation from expected count rates. + * + * **Typical Usage:** + * @code + * nxH4::PNeXus nexus("file.nxs"); + * nxH4::PNeXusDeadTime ndt(&nexus, false); // false = no debug output + * + * if (ndt.IsValid()) { + * auto dims = ndt.GetDimensions(); + * // Minimize for each spectrum + * for (unsigned int i = 0; i < dims[1]; i++) { + * ndt.minimize(i); + * } + * } + * @endcode + * + * @note Requires a valid PNeXus object with count data + * @note Uses ROOT Minuit2 for minimization + * + * @see ROOT::Minuit2::FCNBase + * @see nxH4::PNeXus + */ +class PNeXusDeadTime : public ROOT::Minuit2::FCNBase +{ +public: + /** + * @brief Constructor - initializes dead time calculator from NeXus data + * * Reads count data, time resolution, and bin information from the NeXus file + * to set up the dead time calculation. + * * @param nxs Pointer to PNeXus object containing the detector data + * @param debug If true, print additional debug information during calculations + * * @note After construction, check IsValid() to ensure data was loaded successfully + */ + PNeXusDeadTime(const PNeXus *nxs, bool debug=false); + + /** + * @brief Check if the dead time calculator was initialized successfully + * * @return true if required datasets were found and loaded, false otherwise + */ + bool IsValid() { return fValid; } + + /** + * @brief Get the error definition for the minimizer (FCNBase requirement) + * * Returns the UP parameter which defines how the minimizer estimates errors. + * For chi-square minimization, UP = 0.5 corresponds to 1-sigma errors. + * * @return Error definition parameter (0.5 for chi-square) + */ + double Up() const { return fUp; } + + /** + * @brief Function call operator - calculates chi-square for given parameters + * * This is the objective function that ROOT::Minuit2 minimizes. It calculates + * the chi-square deviation between observed counts and expected counts + * corrected for dead time. + * * @param par Vector of parameters (dead time values in microseconds) + * @return Chi-square value for the given parameters + */ + double operator()(const std::vector &par) const; + + /** + * @brief Minimize dead time for a specific detector spectrum + * * Performs Minuit2 minimization to find the optimal dead time parameter + * for the specified spectrum index. Results are printed to stdout. + * * @param i Spectrum index to minimize (0-based index into second dimension) + * * @note Prints minimization results including dead time estimate and errors + * @note The index i must be less than dims[1] (number of spectra) + */ + void minimize(const int i); + + /** + * @brief Get the dimensions of the count dataset + * * Returns the dimensions [periods, spectra, bins] of the count data + * being used for dead time calculation. + * * @return Vector of dimension sizes (typically 3D: [periods, spectra, bins]) + */ + const std::vector& GetDimensions() const { return fDims; } + +private: + bool fDebug{false}; ///< Debug flag - if true, print additional diagnostic information + bool fValid; ///< Validity flag - true if required data was loaded successfully + int fGoodFrames{-1}; ///< Number of good time frames for analysis + float fTimeResolution{0.0}; ///< Time resolution in picoseconds + std::vectorfDims = {0, 0, 0}; ///< Dimensions of count data [periods, spectra, bins] + int fT0Bin{-1}; ///< T0 bin number (time zero reference) + int fFgbBin{-1}; ///< First good bin for analysis + int fLgbBin{-1}; ///< Last good bin for analysis + std::vector fDeadTime; ///< Dead time values per detector (microseconds) + std::vector fCounts; ///< Count data for minimization + double fUp{0.5}; ///< UP parameter for Minuit2 (0.5 for chi-square) + unsigned int fIdx{0}; ///< Current spectrum index being minimized +}; + + +/** + * @brief HDF4 data type enumeration + * + * This enum maps to HDF4 numeric types (DFNT_*) for type-safe dataset handling + */ +enum class H4DataType { + INT32, ///< 32-bit signed integer (DFNT_INT32) + FLOAT32, ///< 32-bit floating point (DFNT_FLOAT32) + FLOAT64, ///< 64-bit floating point (DFNT_FLOAT64) + CHAR8, ///< 8-bit character (DFNT_CHAR8) + UINT32, ///< 32-bit unsigned integer (DFNT_UINT32) + INT16, ///< 16-bit signed integer (DFNT_INT16) + UINT16, ///< 16-bit unsigned integer (DFNT_UINT16) + INT8, ///< 8-bit signed integer (DFNT_INT8) + UINT8 ///< 8-bit unsigned integer (DFNT_UINT8) +}; + +/** + * @class PNXdata + * @brief Template class for storing HDF4 dataset content with attributes + * + * The PNXdata class stores data read from an HDF4 dataset along with metadata + * such as dimensions, datatype, and attributes. It is designed to be stored in + * std::map where the key is the HDF4 path. + * + * @tparam T The element type of the dataset (int, float, std::string, etc.) + * + * Key features: + * - Stores multi-dimensional data as flattened std::vector + * - Preserves dimension information (shape) + * - Stores HDF4 DataType for type safety + * - Supports attributes as nested PNXdata objects + * - Compatible with std::any for type-erased storage + * + * @example + * @code + * // Create PNXdata for a 3D integer dataset (1 x 96 x 2048) + * PNXdata counts; + * counts.SetDimensions({1, 96, 2048}); + * counts.SetData(dataVector); + * + * // Store in map + * std::map dataMap; + * dataMap["/raw_data_1/detector_1/counts"] = counts; + * @endcode + */ +template class PNXdata { +public: + /** + * @brief Default constructor + */ + PNXdata() : fDataType(H4DataType::INT32) {} + + /** + * @brief Constructor with datatype + * @param dataType HDF4 datatype for this dataset + */ + PNXdata(const H4DataType& dataType) : fDataType(dataType) {} + + /** + * @brief Get the HDF4 DataType + * @return The HDF4 DataType for this dataset + */ + H4DataType GetDataType() const { return fDataType; } + + /** + * @brief Set the HDF4 DataType + * @param dataType The HDF4 DataType to set + */ + void SetDataType(const H4DataType& dataType) { fDataType = dataType; } + + /** + * @brief Get the data as a vector + * @return Reference to the data vector + */ + const std::vector& GetData() const { return fData; } + + /** + * @brief Get mutable reference to the data vector + * @return Reference to the data vector + */ + std::vector& GetData() { return fData; } + + /** + * @brief Set the data vector + * @param data Vector of data to store + */ + void SetData(const std::vector& data) { fData = data; } + + /** + * @brief Get the dimensions of the dataset + * @return Vector of dimension sizes + */ + const std::vector& GetDimensions() const { return fDimensions; } + + /** + * @brief Set the dimensions of the dataset + * @param dims Vector of dimension sizes + */ + void SetDimensions(const std::vector& dims) { fDimensions = dims; } + + /** + * @brief Get the number of dimensions + * @return Number of dimensions (rank) + */ + size_t GetRank() const { return fDimensions.size(); } + + /** + * @brief Get total number of elements + * @return Product of all dimensions + */ + size_t GetNumElements() const { + if (fDimensions.empty()) return 0; + size_t total = 1; + for (auto dim : fDimensions) { + total *= dim; + } + return total; + } + + /** + * @brief Add an attribute to this dataset + * @param name Attribute name + * @param value Attribute value (stored as std::any to support different types) + */ + void AddAttribute(const std::string& name, const std::any& value) { + fAttributes[name] = value; + } + + /** + * @brief Get an attribute by name + * @param name Attribute name + * @return The attribute value as std::any + * @throws std::out_of_range if attribute doesn't exist + */ + std::any GetAttribute(const std::string& name) const { + return fAttributes.at(name); + } + + /** + * @brief Check if an attribute exists + * @param name Attribute name + * @return true if attribute exists, false otherwise + */ + bool HasAttribute(const std::string& name) const { + return fAttributes.find(name) != fAttributes.end(); + } + + /** + * @brief Get all attributes + * @return Map of attribute names to values + */ + const std::map& GetAttributes() const { + return fAttributes; + } + + /** + * @brief Get mutable reference to all attributes + * @return Mutable map of attribute names to values + */ + std::map& GetAttributes() { + return fAttributes; + } + +private: + H4DataType fDataType; ///< HDF4 datatype of the dataset + std::vector fData; ///< Data storage (flattened multi-dimensional array) + std::vector fDimensions; ///< Dimensions of the dataset + std::map fAttributes; ///< Attributes associated with this dataset +}; + +/** + * @class PNeXus + * @brief NeXus HDF4 file reader with case-insensitive path lookup + * + * The PNeXus class provides functionality for reading ISIS muon NeXus HDF4 files + * using the HDF4 C API. It implements case-insensitive path lookup for + * datasets, groups, and attributes to handle files with varying path casings. + * + * Key features: + * - Automatic file reading upon construction + * - Case-insensitive path resolution for datasets and groups + * - Case-insensitive attribute name matching + * - Comprehensive exception handling + * + * @example + * @code + * PNeXus nexus("data/experiment.nxs"); + * // Case-insensitive paths work: "/RAW_DATA_1/IDF_version" -> "/raw_data_1/IDF_version" + * @endcode + * + * @note All path lookups throw appropriate exceptions on failure + */ +class PNeXus { +public: + PNeXus(); + + /** + * @brief Constructor - creates PNeXus object and reads the NeXus file + * @param fln Filename of the NeXus HDF4 file to read + * @param printDebug Enable debug output + * @throws std::runtime_error if file cannot be opened or read + */ + PNeXus(const std::string fln, const bool printDebug=false); + + /** + * @brief Destructor - closes HDF4 file if open + */ + ~PNeXus(); + + /** + * @brief Get the filename of the NeXus file + * @return The filename as a string + */ + std::string GetFileName() const { return fFileName; } + + /** + * @brief Get the hdf4 version of the NeXus file + * @return The hdf4 version as a string + */ + std::string GetHdf4Version() const { return fHdf4Version; } + + /** + * @brief Get the NeXus version of the file + * @return The NeXus version as a string + */ + std::string GetNeXusVersion() const { return fNeXusVersion; } + + /** + * @brief Get the Idf version of the file + * @return The Idf version tag + */ + int GetIdfVersion() const { return fIdfVersion; } + + /** + * @brief Read and parse the NeXus HDF4 file + * @return 0 on success, 1 on error + * @throws std::runtime_error if file cannot be opened or read + */ + int ReadNexusFile(); + + /** + * @brief Get the data map containing all datasets + * @return Reference to the data map (path -> PNXdata stored in std::any) + */ + const std::map& GetDataMap() const { return fDataMap; } + + /** + * @brief Get mutable reference to the data map + * @return Reference to the data map + */ + std::map& GetDataMap() { return fDataMap; } + + /** + * @brief Check if a dataset path exists in the data map + * @param path HDF4 path to check + * @return true if path exists, false otherwise + */ + bool HasDataset(const std::string& path) const { return fDataMap.find(path) != fDataMap.end(); } + + /** + * @brief Add or update a dataset in the data map + * @tparam T Data type of the PNXdata object + * @param path HDF4 path for the dataset + * @param data PNXdata object to store + */ + template + void SetDataset(const std::string& path, const PNXdata& data) { + fDataMap[path] = data; + } + + /** + * @brief Get a dataset from the data map + * @tparam T Data type of the PNXdata object + * @param path HDF4 path of the dataset + * @return PNXdata object + * @throws std::out_of_range if path doesn't exist + * @throws std::bad_any_cast if type T doesn't match stored type + */ + template + PNXdata GetDataset(const std::string& path) const { + return std::any_cast>(fDataMap.at(path)); + } + + /** + * @brief Get a mutable reference to a dataset in the data map + * @tparam T Data type of the PNXdata object + * @param path HDF4 path of the dataset + * @return Reference to PNXdata object + * @throws std::out_of_range if path doesn't exist + * @throws std::bad_any_cast if type T doesn't match stored type + */ + template + PNXdata& GetDatasetRef(const std::string& path) { + return std::any_cast&>(fDataMap.at(path)); + } + + /** + * @brief Remove a dataset from the data map + * @param path HDF4 path of the dataset to remove + * @return true if dataset was removed, false if path didn't exist + */ + bool RemoveDataset(const std::string& path) { + auto it = fDataMap.find(path); + if (it != fDataMap.end()) { + fDataMap.erase(it); + return true; + } + return false; + } + + /** + * @brief Clear all datasets from the data map + */ + void ClearDataMap() { fDataMap.clear(); } + + /** + * @brief Get the number of datasets in the data map + * @return Number of datasets stored + */ + size_t GetNumDatasets() const { return fDataMap.size(); } + + /** + * @brief Update the data of an existing dataset + * @tparam T Data type of the PNXdata object + * @param path HDF4 path of the dataset + * @param newData New data vector to set + * @return true if update succeeded, false if path doesn't exist + */ + template + bool UpdateDatasetData(const std::string& path, const std::vector& newData) { + if (!HasDataset(path)) return false; + try { + auto& dataset = std::any_cast&>(fDataMap.at(path)); + dataset.SetData(newData); + return true; + } catch (const std::bad_any_cast&) { + return false; + } + } + + /** + * @brief Update the dimensions of an existing dataset + * @tparam T Data type of the PNXdata object + * @param path HDF4 path of the dataset + * @param newDimensions New dimensions vector to set + * @return true if update succeeded, false if path doesn't exist + */ + template + bool UpdateDatasetDimensions(const std::string& path, const std::vector& newDimensions) { + if (!HasDataset(path)) return false; + try { + auto& dataset = std::any_cast&>(fDataMap.at(path)); + dataset.SetDimensions(newDimensions); + return true; + } catch (const std::bad_any_cast&) { + return false; + } + } + + /** + * @brief Add or update an attribute for a dataset + * @tparam T Data type of the PNXdata object + * @param path HDF4 path of the dataset + * @param attrName Attribute name + * @param attrValue Attribute value (stored as std::any) + * @return true if attribute was added, false if dataset doesn't exist + */ + template + bool AddDatasetAttribute(const std::string& path, const std::string& attrName, + const std::any& attrValue) { + if (!HasDataset(path)) return false; + try { + auto& dataset = std::any_cast&>(fDataMap.at(path)); + dataset.AddAttribute(attrName, attrValue); + return true; + } catch (const std::bad_any_cast&) { + return false; + } + } + + /** + * @brief Remove an attribute from a dataset + * @tparam T Data type of the PNXdata object + * @param path HDF4 path of the dataset + * @param attrName Attribute name to remove + * @return true if attribute was removed, false if dataset doesn't exist or attribute not found + */ + template + bool RemoveDatasetAttribute(const std::string& path, const std::string& attrName) { + if (!HasDataset(path)) return false; + try { + auto& dataset = std::any_cast&>(fDataMap.at(path)); + auto& attrs = dataset.GetAttributes(); + auto it = attrs.find(attrName); + if (it != attrs.end()) { + attrs.erase(it); + return true; + } + return false; + } catch (const std::bad_any_cast&) { + return false; + } + } + + /** + * @brief Create and add a new dataset with data, dimensions, and optional attributes + * @tparam T Data type of the PNXdata object + * @param path HDF4 path for the new dataset + * @param data Data vector + * @param dimensions Dimensions vector + * @param dataType HDF4 DataType (optional, uses default if not provided) + * @return true if dataset was created, false if path already exists + */ + template + bool AddDataset(const std::string& path, const std::vector& data, + const std::vector& dimensions, + const H4DataType& dataType = H4DataType::INT32) { + if (HasDataset(path)) return false; + PNXdata dataset(dataType); + dataset.SetData(data); + dataset.SetDimensions(dimensions); + fDataMap[path] = dataset; + return true; + } + + /** + * @brief Modify an existing dataset's data and dimensions + * @tparam T Data type of the PNXdata object + * @param path HDF4 path of the dataset + * @param newData New data vector + * @param newDimensions New dimensions vector + * @return true if modification succeeded, false if path doesn't exist + */ + template + bool ModifyDataset(const std::string& path, const std::vector& newData, + const std::vector& newDimensions) { + if (!HasDataset(path)) return false; + try { + auto& dataset = std::any_cast&>(fDataMap.at(path)); + dataset.SetData(newData); + dataset.SetDimensions(newDimensions); + return true; + } catch (const std::bad_any_cast&) { + return false; + } + } + + /** + * @brief Dump the read hdf4-NeXus file content which was read + */ + void Dump(); + + /** + * @brief Write the data map contents to a NeXus HDF4 file + * @param filename Path to the output NeXus HDF4 file + * @param idfVersion IDF version to write (default: 2) + * @return 0 on success, 1 on error + * @throws std::runtime_error if file cannot be created or written + */ + int WriteNexusFile(const std::string& filename, int idfVersion = 2); + + /** + * @brief Add or update an attribute for a group + * @param groupPath HDF4 path of the group (e.g., "/raw_data_1") + * @param attrName Attribute name + * @param attrValue Attribute value (stored as std::any) + * @return true if attribute was added successfully + * @example + * @code + * nexus.AddGroupAttribute("/raw_data_1", "NX_class", std::string("NXentry")); + * @endcode + */ + bool AddGroupAttribute(const std::string& groupPath, const std::string& attrName, + const std::any& attrValue); + + /** + * @brief Remove an attribute from a group + * @param groupPath HDF4 path of the group + * @param attrName Attribute name to remove + * @return true if attribute was removed, false if group or attribute not found + */ + bool RemoveGroupAttribute(const std::string& groupPath, const std::string& attrName); + + /** + * @brief Check if a group has a specific attribute + * @param groupPath HDF4 path of the group + * @param attrName Attribute name to check + * @return true if attribute exists, false otherwise + */ + bool HasGroupAttribute(const std::string& groupPath, const std::string& attrName) const; + + /** + * @brief Get an attribute value from a group + * @param groupPath HDF4 path of the group + * @param attrName Attribute name + * @return The attribute value as std::any + * @throws std::out_of_range if group or attribute doesn't exist + */ + std::any GetGroupAttribute(const std::string& groupPath, const std::string& attrName) const; + + /** + * @brief Get all attributes for a group + * @param groupPath HDF4 path of the group + * @return Map of attribute names to values + */ + const std::map& GetGroupAttributes(const std::string& groupPath) const; + + /** + * @brief Clear all attributes from a group + * @param groupPath HDF4 path of the group + * @return true if group was found and attributes cleared, false otherwise + */ + bool ClearGroupAttributes(const std::string& groupPath); + + /** + * @brief Add or update an attribute at the root level "/" + * @param attrName Attribute name + * @param attrValue Attribute value (stored as std::any) + * @return true if attribute was added successfully + * @note This is a convenience method that calls AddGroupAttribute("/", attrName, attrValue) + * @example + * @code + * nexus.AddRootAttribute("experiment_id", std::string("EXP-2024-001")); + * nexus.AddRootAttribute("temperature", 293.15f); + * nexus.AddRootAttribute("scan_number", static_cast(42)); + * @endcode + */ + bool AddRootAttribute(const std::string& attrName, const std::any& attrValue); + +private: + bool fPrintDebug{false}; ///< if true print additional debug information + std::string fFileName{""}; ///< NeXus HDF4 filename + int fIdfVersion{-1}; ///< IDF version of the NeXus file + std::string fHdf4Version{""}; ///< HDF4 version of the file + std::string fNeXusVersion{""}; ///< NeXus version of the file + std::string fFileNameNxs{""}; + std::string fFileTimeNxs{""}; + std::string fCreatorNxs{""}; + std::string fUserV1{""}; + int32 fSdId{-1}; ///< HDF4 SD interface identifier + int32 fFileId{-1}; ///< HDF4 file identifier + std::map fDataMap; ///< Map of HDF4 paths to PNXdata objects + std::map> fGroupAttributes; ///< Map of group paths to their attributes + + /** + * @brief HandleIdfV1 + * @param sd_id HDF4 SD interface ID + */ + void HandleIdfV1(int32 sd_id); + + /** + * @brief HandleIdfV2 + * @param sd_id HDF4 SD interface ID + */ + void HandleIdfV2(int32 sd_id); + + // ======================================================================== + // Write methods for HDF4 file creation + // ======================================================================== + + /** + * @brief Write dataset attributes from PNXdata object + * @tparam T The data type of the PNXdata object + * @param sds_id HDF4 dataset ID to write attributes to + * @param data PNXdata object containing attributes + */ + template + void WriteDatasetAttributes(int32 sds_id, const PNXdata& data); + + /** + * @brief Write attributes to a group + * @param sd_id HDF4 SD interface ID + * @param groupPath Group path + * @param attributes Map of attribute names to values + */ + void WriteGroupAttributes(int32 sd_id, const std::string& groupPath, + const std::map& attributes); + + /** + * @brief Write an integer dataset with attributes + * @param sd_id HDF4 SD interface ID + * @param path Dataset path + * @param data PNXdata object containing data, dimensions, and attributes + * @throws std::runtime_error if writing fails + */ + void WriteIntDataset(int32 sd_id, const std::string& path, + const PNXdata& data); + + /** + * @brief Write a float dataset with attributes + * @param sd_id HDF4 SD interface ID + * @param path Dataset path + * @param data PNXdata object containing data, dimensions, and attributes + * @throws std::runtime_error if writing fails + */ + void WriteFloatDataset(int32 sd_id, const std::string& path, + const PNXdata& data); + + /** + * @brief Write a string dataset with attributes + * @param sd_id HDF4 SD interface ID + * @param path Dataset path + * @param data PNXdata object containing data, dimensions, and attributes + * @throws std::runtime_error if writing fails + */ + void WriteStringDataset(int32 sd_id, const std::string& path, + const PNXdata& data); + + /** + * @brief Write root-level file attributes + * @param sd_id HDF4 SD interface ID + * @throws std::runtime_error if writing fails + */ + void WriteFileAttributes(int32 sd_id); + + /** + * @brief Write IDF version 2 structure + * @param sd_id HDF4 SD interface ID + * @throws std::runtime_error if writing fails + */ + void WriteIdfV2(int32 sd_id); + + // ======================================================================== + // Case-insensitive lookup helper methods + // ======================================================================== + + /** + * @brief Compare two strings case-insensitively + * @param a First string to compare + * @param b Second string to compare + * @return true if strings are equal (ignoring case), false otherwise + * @note Uses std::tolower for character-by-character comparison + */ + static bool caseInsensitiveEquals(const std::string& a, const std::string& b); + + /** + * @brief Split an HDF4 path into components + * @param path HDF4 path string (e.g., "/raw_data_1/IDF_version") + * @return Vector of path components + * @note Empty first component indicates absolute path + * @example "/raw_data_1/IDF_version" -> ["", "raw_data_1", "IDF_version"] + */ + static std::vector splitPath(const std::string& path); + + /** + * @brief Find attribute name with case-insensitive matching + * @param sd_id HDF4 SD interface ID + * @param requestedName Requested attribute name (any case) + * @return Correctly-cased attribute name as it exists in the file + * @throws std::runtime_error if attribute not found + * @example "nexus_version" might resolve to "NeXus_version" + */ + std::string findAttributeName(int32 sd_id, const std::string& requestedName); + + /** + * @brief Find dataset index with case-insensitive matching + * @param sd_id HDF4 SD interface ID + * @param requestedName Requested dataset name (any case) + * @return Dataset index + * @throws std::runtime_error if dataset not found + */ + int32 findDatasetIndex(int32 sd_id, const std::string& requestedName); + + /** + * @brief Find dataset reference by navigating VGroup hierarchy + * @param path Full hierarchical path (e.g., "/run/sample/name") + * @return Dataset reference number, or -1 if not found + * @note Uses VGroup navigation to resolve paths in HDF4 files with duplicate dataset names + */ + int32 findDatasetRefByPath(const std::string& path); + + // ======================================================================== + // Dataset reading helper methods + // ======================================================================== + + /** + * @brief Read an integer dataset and store in data map + * @param sd_id HDF4 SD interface ID + * @param path Path to the dataset + * @throws std::runtime_error if reading fails + */ + void ReadIntDataset(int32 sd_id, const std::string& path); + + /** + * @brief Read a float dataset and store in data map + * @param sd_id HDF4 SD interface ID + * @param path Path to the dataset + * @throws std::runtime_error if reading fails + */ + void ReadFloatDataset(int32 sd_id, const std::string& path); + + /** + * @brief Read a string dataset and store in data map + * @param sd_id HDF4 SD interface ID + * @param path Path to the dataset + * @throws std::runtime_error if reading fails + */ + void ReadStringDataset(int32 sd_id, const std::string& path); + + /** + * @brief Read dataset attributes and add to PNXdata object + * @tparam T The data type of the PNXdata object + * @param sds_id HDF4 dataset ID + * @param data PNXdata object to add attributes to + */ + template + void ReadDatasetAttributes(int32 sds_id, PNXdata& data); + + /** + * @brief Convert HDF4 data type to H4DataType enum + * @param hdf4_type HDF4 numeric type constant (DFNT_*) + * @return H4DataType enum value + */ + static H4DataType convertHdf4Type(int32 hdf4_type); + + /** + * @brief Convert H4DataType enum to HDF4 data type + * @param dataType H4DataType enum value + * @return HDF4 numeric type constant (DFNT_*) + */ + static int32 convertToHdf4Type(H4DataType dataType); +}; + +} +#endif // HAVE_HDF4 + #endif // _PNEXUS_H_