From c3e9c03920f507705aa2b34ecb9c111c8ba49ff0 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 25 Jan 2026 11:15:49 +0100 Subject: [PATCH] added HDF5 ISIS NeXus class handling. --- src/external/nexus/CMakeLists.txt | 8 +- src/external/nexus/PNeXus.cpp | 2721 ++++++++++++++++++++++++++++- src/external/nexus/PNeXus.h | 835 +++++++++ 3 files changed, 3562 insertions(+), 2 deletions(-) diff --git a/src/external/nexus/CMakeLists.txt b/src/external/nexus/CMakeLists.txt index 3e866987..140a82b1 100644 --- a/src/external/nexus/CMakeLists.txt +++ b/src/external/nexus/CMakeLists.txt @@ -28,7 +28,13 @@ target_include_directories( ) #--- add library dependencies ------------------------------------------------- -target_link_libraries(PNeXus ${HDF4_LIBRARIES} ${ROOT_LIBRARIES}) +if (HAVE_HDF4) + set(HDF_LIBS ${HDF4_LIBRARIES} ${HDF5_LIBRARIES}) +else (HAVE_HDF4) + set(HDF_LIBS ${HDF5_LIBRARIES}) +endif (HAVE_HDF4) +message(STATUS "++>> HDF_LIBS: ${HDF_LIBS}") +target_link_libraries(PNeXus ${HDF_LIBS} ${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 32c5c8c3..80ae1950 100644 --- a/src/external/nexus/PNeXus.cpp +++ b/src/external/nexus/PNeXus.cpp @@ -34,6 +34,7 @@ #include #include +#include "Minuit2/MnStrategy.h" #include "Minuit2/MnMinimize.h" #include "Minuit2/FunctionMinimum.h" @@ -77,7 +78,7 @@ nxs::HDFType nxs::checkHDFType(const std::string& filename) { #ifdef HAVE_HDF4 //============================================================================= -// PNeXusDeadTime Constructor +// nxH4::PNeXusDeadTime Constructor //============================================================================= nxH4::PNeXusDeadTime::PNeXusDeadTime(const nxH4::PNeXus *nxs, bool debug) : fDebug(debug) { @@ -2261,3 +2262,2721 @@ int32 nxH4::PNeXus::convertToHdf4Type(H4DataType dataType) } } #endif // HAVE_HDF4 + +// ** HDF5 ******************************************************************** + +//============================================================================= +// nxH5::PNeXusDeadTime Constructor +//============================================================================= +nxH5::PNeXusDeadTime::PNeXusDeadTime(const 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; + + // get t0 + if (counts_data.HasAttribute("t0_bin")) { + try { + auto t0_bin = std::any_cast(counts_data.GetAttribute("t0_bin")); + fT0Bin = t0_bin; + if (fDebug) + std::cout << " fT0Bin: " << fT0Bin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast t0_bin attribute" << std::endl; + } + } + // get fgb + if (counts_data.HasAttribute("first_good_bin")) { + try { + auto first_good_bin = std::any_cast(counts_data.GetAttribute("first_good_bin")); + fFgbBin = first_good_bin; + if (fDebug) + std::cout << " fFgbBin: " << fFgbBin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast first_good_bin attribute" << std::endl; + } + } + // get lgb + if (counts_data.HasAttribute("last_good_bin")) { + try { + auto last_good_bin = std::any_cast(counts_data.GetAttribute("last_good_bin")); + fLgbBin = last_good_bin; + if (fDebug) + std::cout << " fLgbBin: " << fLgbBin << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast last_good_bin attribute" << std::endl; + } + } + } + + // get dead time parameters + if (dataMap.find("/raw_data_1/detector_1/dead_time") != dataMap.end()) { + auto dead_time = std::any_cast>(dataMap["/raw_data_1/detector_1/dead_time"]); + const auto& dt_data = dead_time.GetData(); + fDeadTime = dt_data; + } + fDeadTimeEstimated.resize(fDeadTime.size()); + + // get good_frames + if (dataMap.find("/raw_data_1/good_frames") != dataMap.end()) { + auto intData = std::any_cast>(dataMap["/raw_data_1/good_frames"]); + const auto& ival = intData.GetData(); + fGoodFrames = ival[0]; + if (fDebug) + std::cout << " good_frames: " << fGoodFrames << std::endl; + } + // get resolution with units + if (dataMap.find("/raw_data_1/instrument/detector_1/resolution") != dataMap.end()) { + auto int_data = std::any_cast>(dataMap["/raw_data_1/instrument/detector_1/resolution"]); + const auto& ival = int_data.GetData(); + float scale{1.0}; + if (int_data.HasAttribute("units")) { + try { + auto units = std::any_cast(int_data.GetAttribute("units")); + if (units == "picoseconds") + scale = 1.0e-6; + else if (units == "nanoseconds") + scale = 1.0e-3; + } catch (const std::bad_any_cast& e) { + std::cerr << "Error: Failed to cast units attribute" << std::endl; + } + } + fTimeResolution = ival[0] * scale; + if (fDebug) + std::cout << " time_resolution: " << fTimeResolution << " (us)" << std::endl; + } + } +} + +//============================================================================= +// nxH5::PNeXusDeadTime::operator() +//============================================================================= +/** + * @brief Calculates the negative log-likelihood with dead time correction parameter dtc. + * @param par + * @return + */ +double nxH5::PNeXusDeadTime::operator()(const std::vector &par) const +{ + // par[0]: dtc, par[1]: N0, par[2]: N_bkg + + double nt=0, nd=0, ll=0; + const double tau_mu=2.1969811; + for (unsigned int i=fFgbBin; i::max(); + double tolerance = 0.1; + ROOT::Minuit2::FunctionMinimum min = minimize(maxfcn, tolerance); + + std::unique_ptr fcnMin; ///< Minuit2 function minimum result + + // keep FunctionMinimum object + fcnMin.reset(new ROOT::Minuit2::FunctionMinimum(min)); + + // keep user parameters + if (fcnMin) + params = fcnMin->UserParameters(); + + fDeadTimeEstimated[i] = params.Value(0); + + if (fDebug) + std::cout << "debug> " << i << ": dtc fitted=" << params.Value(0) << ", dtc from file=" << fDeadTime[i] << std::endl; +} + +//============================================================================= +// Case-insensitive string comparison helper +//============================================================================= +/** + * Compares two strings character-by-character using case-insensitive matching. + * The function first checks if the lengths match, then converts each character + * to lowercase using std::tolower before comparison. + * + * @param a First string to compare + * @param b Second string to compare + * @return true if strings are equal (ignoring case), false otherwise + * + * @note Uses static_cast to unsigned char for std::tolower to avoid undefined + * behavior with negative char values + * + * @example + * @code + * caseInsensitiveEquals("NeXus", "nexus") // returns true + * caseInsensitiveEquals("IDF_VERSION", "idf_version") // returns true + * @endcode + */ +bool nxH5::PNeXus::caseInsensitiveEquals(const std::string& a, const std::string& b) +{ + if (a.length() != b.length()) { + return false; + } + + for (size_t i = 0; i < a.length(); ++i) { + if (std::tolower(static_cast(a[i])) != + std::tolower(static_cast(b[i]))) { + return false; + } + } + + return true; +} + +//============================================================================= +// Split HDF5 path into components +//============================================================================= +/** + * Splits an HDF5 path string into individual components using '/' as delimiter. + * The function preserves the information about absolute vs relative paths by + * keeping an empty first component for absolute paths. + * + * @param path HDF5 path string to split + * @return Vector of path components + * + * @note Empty components are filtered out except for the first component + * which indicates an absolute path when empty + * + * @example + * @code + * splitPath("/raw_data_1/IDF_version") // returns ["", "raw_data_1", "IDF_version"] + * splitPath("detector_1/counts") // returns ["detector_1", "counts"] + * splitPath("/") // returns [""] + * @endcode + */ +std::vector nxH5::PNeXus::splitPath(const std::string& path) +{ + std::vector components; + std::string component; + std::istringstream stream(path); + + while (std::getline(stream, component, '/')) { + // Keep empty first component to indicate absolute path + if (!component.empty() || components.empty()) { + components.push_back(component); + } + } + + return components; +} + +//============================================================================= +// Find attribute with case-insensitive name matching +//============================================================================= +/** + * Searches for an attribute by name using case-insensitive matching. + * Lists all attributes on the given HDF5 object and compares each name + * case-insensitively against the requested name. + * + * @param obj HDF5 File object to search for attributes + * @param requestedName Attribute name to find (any case) + * @return The correctly-cased attribute name as it exists in the file + * @throws H5::AttributeIException if no matching attribute is found + * + * @note Useful for handling NeXus files where attribute names may vary in case + * (e.g., "NeXus_version" vs "nexus_version") + * + * @example + * @code + * std::string actualName = findAttributeName(file, "nexus_VERSION"); + * // actualName might be "NeXus_version" if that's how it's stored + * @endcode + */ +std::string nxH5::PNeXus::findAttributeName(H5::H5File& obj, const std::string& requestedName) +{ + int numAttrs = obj.getNumAttrs(); + + for (int i = 0; i < numAttrs; i++) { + H5::Attribute attr = obj.openAttribute(i); + std::string attrName = attr.getName(); + attr.close(); + + if (caseInsensitiveEquals(attrName, requestedName)) { + return attrName; + } + } + + // Not found - throw exception + std::ostringstream msg; + msg << "Could not find attribute matching case-insensitive name: '" << requestedName << "'"; + throw H5::AttributeIException("findAttributeName", msg.str()); +} + +//============================================================================= +// Find group path with case-insensitive matching (File overload) +//============================================================================= +/** + * Finds an HDF5 group path using case-insensitive matching, starting from + * the file root. Traverses the group hierarchy by splitting the path into + * components and matching each component case-insensitively. + * + * @param parent HDF5 File object (file root) + * @param requestedPath Group path to find (any case, absolute or relative) + * @return The correctly-cased group path as it exists in the file + * @throws H5::GroupIException if group not found or path component is not a group + * + * @note Each path component is verified to be a Group using getObjectType() + * + * @example + * @code + * std::string path = findGroupPath(file, "/RAW_DATA_1/detector_1"); + * // path might be "/raw_data_1/detector_1" if that's the actual casing + * @endcode + */ +std::string nxH5::PNeXus::findGroupPath(H5::H5File& parent, const std::string& requestedPath) +{ + std::vector components = splitPath(requestedPath); + + // Handle empty path + if (components.empty()) { + return "/"; + } + + // Check if absolute path + bool isAbsolute = (components[0].empty()); + std::string currentPath = isAbsolute ? "" : ""; + + // Start from index 1 if absolute path (skip empty component) + size_t startIdx = isAbsolute ? 1 : 0; + + // Traverse path components + for (size_t i = startIdx; i < components.size(); ++i) { + const std::string& requestedComponent = components[i]; + + // Build the parent path for querying + std::string parentPath = (currentPath.empty() || currentPath == "") ? "/" : currentPath; + H5::Group currentGroup = (parentPath == "/") ? parent.openGroup("/") : parent.openGroup(parentPath); + + // List object names using iteration + hsize_t numObjs = currentGroup.getNumObjs(); + std::vector objectNames; + for (hsize_t j = 0; j < numObjs; j++) { + objectNames.push_back(currentGroup.getObjnameByIdx(j)); + } + + // Find case-insensitive match + std::string matchedName; + bool found = false; + + for (const auto& objName : objectNames) { + if (caseInsensitiveEquals(objName, requestedComponent)) { + matchedName = objName; + found = true; + break; + } + } + + if (!found) { + std::ostringstream msg; + msg << "Could not find group component matching case-insensitive path: '" + << requestedPath << "' (failed at component: '" << requestedComponent << "')"; + throw H5::GroupIException("findGroupPath", msg.str()); + } + + // Build current path + currentPath += "/" + matchedName; + + // Verify it's a group + H5O_type_t objType = currentGroup.childObjType(matchedName); + if (objType != H5O_TYPE_GROUP) { + std::ostringstream msg; + msg << "Path component '" << matchedName << "' in '" << requestedPath + << "' is not a group"; + throw H5::GroupIException("findGroupPath", msg.str()); + } + } + + return currentPath; +} + +//----------------------------------------------------------------------------- +// Find group path with case-insensitive matching (Group overload) +//----------------------------------------------------------------------------- +std::string nxH5::PNeXus::findGroupPath(H5::Group& parent, const std::string& requestedPath) +{ + std::vector components = splitPath(requestedPath); + + // Handle empty path + if (components.empty()) { + return "/"; + } + + // Check if absolute path + bool isAbsolute = (components[0].empty()); + std::string currentPath = isAbsolute ? "" : ""; + H5::Group currentGroup = parent; + + // Start from index 1 if absolute path (skip empty component) + size_t startIdx = isAbsolute ? 1 : 0; + + // Traverse path components + for (size_t i = startIdx; i < components.size(); ++i) { + const std::string& requestedComponent = components[i]; + + // List object names using iteration + hsize_t numObjs = currentGroup.getNumObjs(); + std::vector objectNames; + for (hsize_t j = 0; j < numObjs; j++) { + objectNames.push_back(currentGroup.getObjnameByIdx(j)); + } + + // Find case-insensitive match + std::string matchedName; + bool found = false; + + for (const auto& objName : objectNames) { + if (caseInsensitiveEquals(objName, requestedComponent)) { + matchedName = objName; + found = true; + break; + } + } + + if (!found) { + std::ostringstream msg; + msg << "Could not find group component matching case-insensitive path: '" + << requestedPath << "' (failed at component: '" << requestedComponent << "')"; + throw H5::GroupIException("findGroupPath", msg.str()); + } + + // Build current path + currentPath += "/" + matchedName; + + // Verify it's a group (unless it's the last component, which is checked by caller) + H5O_type_t objType = currentGroup.childObjType(matchedName); + if (objType != H5O_TYPE_GROUP) { + std::ostringstream msg; + msg << "Path component '" << matchedName << "' in '" << requestedPath + << "' is not a group"; + throw H5::GroupIException("findGroupPath", msg.str()); + } + + // Open the group for next iteration + if (i < components.size() - 1) { + currentGroup = currentGroup.openGroup(matchedName); + } + } + + return currentPath; +} + +//============================================================================= +// Find dataset path with case-insensitive matching (File overload) +//============================================================================= +/** + * Finds an HDF5 dataset path using case-insensitive matching, starting from + * the file root. Traverses the group hierarchy and verifies that the final + * path component is actually a dataset. + * + * @param parent HDF5 File object (file root) + * @param requestedPath Dataset path to find (any case, absolute or relative) + * @return The correctly-cased dataset path as it exists in the file + * @throws H5::DataSetIException if dataset not found, path invalid, or + * final component is not a dataset + * + * @note Intermediate path components must be groups, and the final component + * must be a dataset, otherwise an exception is thrown + * + * @example + * @code + * std::string path = findDatasetPath(file, "/RAW_DATA_1/idf_VERSION"); + * // path might be "/raw_data_1/IDF_version" if that's the actual casing + * std::int data = H5Easy::load(file, path); + * @endcode + */ +std::string nxH5::PNeXus::findDatasetPath(H5::H5File& parent, const std::string& requestedPath) +{ + std::vector components = splitPath(requestedPath); + + // Handle empty path + if (components.empty()) { + throw H5::DataSetIException("findDatasetPath", "Empty dataset path provided"); + } + + // Check if absolute path + bool isAbsolute = (components[0].empty()); + std::string currentPath = isAbsolute ? "" : ""; + + // Start from index 1 if absolute path (skip empty component) + size_t startIdx = isAbsolute ? 1 : 0; + + // Traverse path components + for (size_t i = startIdx; i < components.size(); ++i) { + const std::string& requestedComponent = components[i]; + + // Build the parent path for querying + std::string parentPath = (currentPath.empty() || currentPath == "") ? "/" : currentPath; + H5::Group currentGroup = (parentPath == "/") ? parent.openGroup("/") : parent.openGroup(parentPath); + + // List object names using iteration + hsize_t numObjs = currentGroup.getNumObjs(); + std::vector objectNames; + for (hsize_t j = 0; j < numObjs; j++) { + objectNames.push_back(currentGroup.getObjnameByIdx(j)); + } + + // Find case-insensitive match + std::string matchedName; + bool found = false; + + for (const auto& objName : objectNames) { + if (caseInsensitiveEquals(objName, requestedComponent)) { + matchedName = objName; + found = true; + break; + } + } + + if (!found) { + std::ostringstream msg; + msg << "Could not find dataset matching case-insensitive path: '" + << requestedPath << "' (failed at component: '" << requestedComponent << "')"; + throw H5::DataSetIException("findDatasetPath", msg.str()); + } + + // Build current path + currentPath += "/" + matchedName; + + // Check if this is the last component (should be a dataset) + if (i == components.size() - 1) { + H5O_type_t objType = currentGroup.childObjType(matchedName); + if (objType != H5O_TYPE_DATASET) { + std::ostringstream msg; + msg << "Path '" << requestedPath << "' resolves to '" << currentPath + << "' which is not a dataset"; + throw H5::DataSetIException("findDatasetPath", msg.str()); + } + } else { + // Intermediate component - should be a group + H5O_type_t objType = currentGroup.childObjType(matchedName); + if (objType != H5O_TYPE_GROUP) { + std::ostringstream msg; + msg << "Path component '" << matchedName << "' in '" << requestedPath + << "' is not a group"; + throw H5::DataSetIException("findDatasetPath", msg.str()); + } + } + } + + return currentPath; +} + +//----------------------------------------------------------------------------- +// Find dataset path with case-insensitive matching (Group overload) +//----------------------------------------------------------------------------- +std::string nxH5::PNeXus::findDatasetPath(H5::Group& parent, const std::string& requestedPath) +{ + std::vector components = splitPath(requestedPath); + + // Handle empty path + if (components.empty()) { + throw H5::DataSetIException("findDatasetPath", "Empty dataset path provided"); + } + + // Check if absolute path + bool isAbsolute = (components[0].empty()); + std::string currentPath = isAbsolute ? "" : ""; + H5::Group currentGroup = parent; + + // Start from index 1 if absolute path (skip empty component) + size_t startIdx = isAbsolute ? 1 : 0; + + // Traverse path components + for (size_t i = startIdx; i < components.size(); ++i) { + const std::string& requestedComponent = components[i]; + + // List object names using iteration + hsize_t numObjs = currentGroup.getNumObjs(); + std::vector objectNames; + for (hsize_t j = 0; j < numObjs; j++) { + objectNames.push_back(currentGroup.getObjnameByIdx(j)); + } + + // Find case-insensitive match + std::string matchedName; + bool found = false; + + for (const auto& objName : objectNames) { + if (caseInsensitiveEquals(objName, requestedComponent)) { + matchedName = objName; + found = true; + break; + } + } + + if (!found) { + std::ostringstream msg; + msg << "Could not find dataset matching case-insensitive path: '" + << requestedPath << "' (failed at component: '" << requestedComponent << "')"; + throw H5::DataSetIException("findDatasetPath", msg.str()); + } + + // Build current path + currentPath += "/" + matchedName; + + // Check if this is the last component (should be a dataset) + if (i == components.size() - 1) { + H5O_type_t objType = currentGroup.childObjType(matchedName); + if (objType != H5O_TYPE_DATASET) { + std::ostringstream msg; + msg << "Path '" << requestedPath << "' resolves to '" << currentPath + << "' which is not a dataset"; + throw H5::DataSetIException("findDatasetPath", msg.str()); + } + } else { + // Intermediate component - should be a group + H5O_type_t objType = currentGroup.childObjType(matchedName); + if (objType != H5O_TYPE_GROUP) { + std::ostringstream msg; + msg << "Path component '" << matchedName << "' in '" << requestedPath + << "' is not a group"; + throw H5::DataSetIException("findDatasetPath", msg.str()); + } + // Open the group for next iteration + currentGroup = currentGroup.openGroup(matchedName); + } + } + + return currentPath; +} + +//============================================================================= +// nxH5::PNeXus Constructor +//============================================================================= +nxH5::PNeXus::PNeXus() +{ + // empty +} + +//============================================================================= +// nxH5::PNeXus Constructor +//============================================================================= +/** + * Constructs a PNeXus object and immediately reads the specified NeXus file. + * The constructor initializes the filename member and calls ReadNexusFile() + * to open and parse the HDF5 file. + * + * @param fln Path to the NeXus HDF5 file to read + * @throws H5::FileIException if the file cannot be opened + * @throws H5::AttributeIException if required attributes are missing + * @throws H5::DataSetIException if required datasets are missing + * + * @note The file is read immediately upon construction. If reading fails, + * an exception is thrown and the object is not fully constructed. + * + * @example + * @code + * try { + * PNeXus nexus("data/emu00139040.nxs"); + * std::cout << "Opened: " << nexus.GetFileName() << std::endl; + * } catch (const H5::Exception& e) { + * std::cerr << "Failed to open file: " << e.getDetailMsg() << std::endl; + * } + * @endcode + */ +nxH5::PNeXus::PNeXus(const std::string fln, const bool printDebug) : fFileName(fln), fPrintDebug(printDebug) +{ + ReadNexusFile(); +} + + +//============================================================================= +// Read NeXus File +//============================================================================= +/** + * Reads and parses the NeXus HDF5 file. Opens the file in read-only mode, + * loads the NeXus version attribute, and reads the IDF version dataset. + * Uses case-insensitive path lookup to handle varying path casings. + * + * @return 0 on success, 1 on error + * @throws H5::FileIException if file cannot be opened + * @throws H5::AttributeIException if NeXus_version attribute is missing + * @throws H5::DataSetIException if IDF_version dataset is missing + * + * @note This method is called automatically by the constructor. It performs + * case-insensitive lookups for both attributes and datasets to handle + * NeXus files with varying naming conventions. + * + * @internal + * Current implementation reads: + * - Attribute: NeXus_version (from root "/") + * - Dataset: /raw_data_1/IDF_version + * + * Both paths use case-insensitive lookup, so variations like + * "/RAW_DATA_1/idf_VERSION" will work correctly. + */ +int nxH5::PNeXus::ReadNexusFile() +{ + try { + // Open the HDF5/NeXus file + H5::H5File file(fFileName, H5F_ACC_RDONLY); + + // Turn off the auto-printing when failure occurs so that we can + // handle the errors appropriately + H5::Exception::dontPrint(); + + // Load NeXus version attribute (with case-insensitive lookup) + std::string version; + try { + std::string versionNeXusAttr = findAttributeName(file, "NeXus_version"); + H5::Attribute attr = file.openAttribute(versionNeXusAttr); + H5::DataType dtype = attr.getDataType(); + attr.read(dtype, version); + attr.close(); + fNeXusVersion = version; + if (fPrintDebug) + std::cout << "debug> NeXus version=" << fNeXusVersion << std::endl; + } catch (const H5::AttributeIException& err) { + std::cerr << "Error: Failed to read NeXus_version attribute: " << err.getDetailMsg() << std::endl; + return 1; + } + + // Load IDF version from dataset (with case-insensitive lookup) + int idf_vers{-1}; + try { + std::string idfPath = findDatasetPath(file, "/run/IDF_version"); + H5::DataSet dataset = file.openDataSet(idfPath); + dataset.read(&idf_vers, H5::PredType::NATIVE_INT); + dataset.close(); + fIdfVersion = idf_vers; + } catch (const H5::DataSetIException& err) { + if (fPrintDebug) + std::cout << "Info: not IDF version 1, will check for IDF version 2" << std::endl; + } + + if (fIdfVersion == -1) { + try { + std::string idfPath = findDatasetPath(file, "/raw_data_1/IDF_version"); + H5::DataSet dataset = file.openDataSet(idfPath); + dataset.read(&idf_vers, H5::PredType::NATIVE_INT); + dataset.close(); + fIdfVersion = idf_vers; + } catch (const H5::DataSetIException& err) { + std::cerr << "Error: Failed to read IDF_version dataset: " << err.getDetailMsg() << std::endl; + return 1; + } + + try { + std::string versionHdf5Attr = findAttributeName(file, "HDF5_version"); + H5::Attribute attr = file.openAttribute(versionHdf5Attr); + H5::DataType dtype = attr.getDataType(); + attr.read(dtype, version); + attr.close(); + fHdf5Version = version; + if (fPrintDebug) + std::cout << "debug> HDF5 version=" << fHdf5Version << std::endl; + } catch (const H5::AttributeIException& err) { + std::cerr << "Error: Failed to read HDF5_version attribute: " << err.getDetailMsg() << std::endl; + return 1; + } + } + + if (fPrintDebug) + std::cout << "debug> IDF_version=" << fIdfVersion << std::endl; + + if (fIdfVersion == 1) { + try { + HandleIdfV1(file); + } catch (...) { + + } + } else { // fIdfVersion == 2 + try { + HandleIdfV2(file); + } catch (...) { + + } + } + + return 0; + + } catch (const H5::FileIException& err) { + std::cerr << "Error: Failed to open file '" << fFileName << "': " << err.getDetailMsg() << std::endl; + return 1; + } catch (const H5::Exception& err) { + std::cerr << "Error: HDF5 exception occurred: " << err.getDetailMsg() << std::endl; + return 1; + } catch (const std::exception& err) { + std::cerr << "Error: Unexpected exception: " << err.what() << std::endl; + return 1; + } +} + +//============================================================================= +// Handle NeXus IDF version 1 file +//============================================================================= +void nxH5::PNeXus::HandleIdfV1(H5::H5File &file) +{ + if (fPrintDebug) + std::cout << "debug> in HandleIdfV1 ..." << std::endl; + + std::string attrStr{""}, attrStrName{""}; + // get file_name attribute + try { + attrStrName = findAttributeName(file, "user"); + H5::Attribute attr = file.openAttribute(attrStrName); + H5::DataType dtype = attr.getDataType(); + attr.read(dtype, attrStr); + attr.close(); + fUserV1 = attrStr; + if (fPrintDebug) + std::cout << "debug> user =" << fUserV1 << std::endl; + } catch (const H5::AttributeIException& err) { + std::cerr << "Error: Failed to read user attribute: " << err.getDetailMsg() << std::endl; + } + + // Read the mandatory key datasets and store them in the data map + try { + ReadStringDataset(file, "/run/program_name"); + ReadIntDataset(file, "/run/number"); + ReadStringDataset(file, "/run/title"); + ReadStringDataset(file, "/run/notes"); + ReadStringDataset(file, "/run/analysis"); + ReadStringDataset(file, "/run/lab"); + ReadStringDataset(file, "/run/beamline"); + ReadStringDataset(file, "/run/start_time"); + ReadStringDataset(file, "/run/end_time"); + ReadIntDataset(file, "/run/switching_state"); + ReadStringDataset(file, "/run/user/name"); + ReadStringDataset(file, "/run/user/experiment_number"); + ReadStringDataset(file, "/run/sample/name"); + ReadFloatDataset(file, "/run/sample/temperature"); + ReadFloatDataset(file, "/run/sample/magnetic_field"); + ReadStringDataset(file, "/run/sample/magnetic_field_state"); + ReadFloatDataset(file, "/run/sample/magnetic_field_vector"); + ReadStringDataset(file, "/run/sample/environment"); + ReadStringDataset(file, "/run/instrument/name"); + ReadIntDataset(file, "/run/instrument/detector/number"); + ReadStringDataset(file, "/run/instrument/collimator/type"); + ReadStringDataset(file, "/run/instrument/beam/beamline"); + ReadIntDataset(file, "/run/histogram_data_1/counts"); + ReadIntDataset(file, "/run/histogram_data_1/resolution"); + ReadIntDataset(file, "/run/histogram_data_1/time_zero"); + ReadFloatDataset(file, "/run/histogram_data_1/raw_time"); + ReadFloatDataset(file, "/run/histogram_data_1/corrected_time"); + ReadIntDataset(file, "/run/histogram_data_1/grouping"); + ReadFloatDataset(file, "/run/histogram_data_1/alpha"); + } catch (const H5::Exception& err) { + std::cerr << "Error in HandleIdfV1: " << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Handle NeXus IDF version 2 file +//============================================================================= +void nxH5::PNeXus::HandleIdfV2(H5::H5File &file) +{ + if (fPrintDebug) + std::cout << "debug> in HandleIdfV2 ..." << std::endl; + + std::string attrStr{""}, attrStrName{""}; + // get file_name attribute + try { + attrStrName = findAttributeName(file, "file_name"); + H5::Attribute attr = file.openAttribute(attrStrName); + H5::DataType dtype = attr.getDataType(); + attr.read(dtype, attrStr); + attr.close(); + fFileNameNxs = attrStr; + if (fPrintDebug) + std::cout << "debug> NXS file_name =" << fFileNameNxs << std::endl; + } catch (const H5::AttributeIException& err) { + std::cerr << "Error: Failed to read file_name attribute: " << err.getDetailMsg() << std::endl; + } + + // get file_time attribute + attrStr=""; + try { + attrStrName = findAttributeName(file, "file_time"); + H5::Attribute attr = file.openAttribute(attrStrName); + H5::DataType dtype = attr.getDataType(); + attr.read(dtype, attrStr); + attr.close(); + fFileTimeNxs = attrStr; + if (fPrintDebug) + std::cout << "debug> NXS file_time =" << fFileNameNxs << std::endl; + } catch (const H5::AttributeIException& err) { + std::cerr << "Error: Failed to read file_time attribute: " << err.getDetailMsg() << std::endl; + } + + // Read the mandatory key datasets and store them in the data map + try { + ReadIntDataset(file, "/raw_data_1/IDF_version"); + ReadIntDataset(file, "/raw_data_1/good_frames"); + ReadStringDataset(file, "/raw_data_1/beamline"); + ReadStringDataset(file, "/raw_data_1/definition"); + ReadIntDataset(file, "/raw_data_1/run_number"); + ReadStringDataset(file, "/raw_data_1/title"); + ReadStringDataset(file, "/raw_data_1/start_time"); + ReadStringDataset(file, "/raw_data_1/end_time"); + ReadStringDataset(file, "/raw_data_1/experiment_identifier"); + ReadStringDataset(file, "/raw_data_1/instrument/name"); + ReadStringDataset(file, "/raw_data_1/instrument/source/name"); + ReadStringDataset(file, "/raw_data_1/instrument/source/type"); + ReadStringDataset(file, "/raw_data_1/instrument/source/probe"); + ReadIntDataset(file, "/raw_data_1/instrument/detector_1/resolution"); + ReadIntDataset(file, "/raw_data_1/detector_1/counts"); + ReadFloatDataset(file, "/raw_data_1/detector_1/raw_time"); + ReadIntDataset(file, "/raw_data_1/detector_1/spectrum_index"); + ReadFloatDataset(file, "/raw_data_1/detector_1/dead_time"); + } catch (const H5::Exception& err) { + std::cerr << "Error in HandleIdfV2: " << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Dump hdf5-NeXus file content which was read +//============================================================================= +void nxH5::PNeXus::Dump() +{ + int32_t first_good_bin{0}; + + if (fIdfVersion == 1) { + std::cout << std::endl; + std::cout << std::endl << "hdf5-NeXus file content of file:' " << fFileName << "'"; std::cout << std::endl << "****"; + std::cout << std::endl << "****"; + std::cout << std::endl << "Top Level Attributes:"; + std::cout << std::endl << " NeXus Version: " << fNeXusVersion; + std::cout << std::endl << " user: " << fUserV1; + 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 << std::endl; + } catch (const std::bad_any_cast& e) { + std::cerr << "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/end_time") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/end_time"]); + std::cout << std::endl << " end_time : " << str_data.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("/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 << " user"; + std::cout << std::endl << "----"; + if (fDataMap.find("/run/usr/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/usr/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 << " 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 << 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 << 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 << 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 << 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 << " 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 << std::endl; + } 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 << " 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 << " 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; + } + } + std::cout << std::endl << " collimator"; + std::cout << std::endl << "----"; + if (fDataMap.find("/run/instrument/collimator/type") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(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 << " beam"; + std::cout << std::endl << "----"; + if (fDataMap.find("/run/instrument/beam/beamline") != fDataMap.end()) { + try { + auto str_data = std::any_cast>(fDataMap["/run/instrument/beam/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; + } + } + 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 : " << noOfHistos << 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 << 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 << std::endl << "Error: Failed to cast resolution data" << std::endl; + } + } + + if (fDataMap.find("/run/histogram_data_1/time_zero") != fDataMap.end()) { + try { + auto int_data = std::any_cast>(fDataMap["/run/histogram_data_1/time_zero"]); + std::cout << std::endl << " time_zero : " << int_data.GetData()[0]; + + if (int_data.HasAttribute("units")) { + try { + auto units = std::any_cast(int_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 (int_data.HasAttribute("available")) { + try { + auto available = std::any_cast(int_data.GetAttribute("available")); + std::cout << " available : " << available << std::endl; + } 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 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 << " axis : " << axis << std::endl; + } 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 << " primary : " << primary << std::endl; + } 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 << " units : " << 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 << 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 << " axis : " << axis << std::endl; + } 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 << " units : " << 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 << 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"]); + std::cout << std::endl << " grouping : " << int_data.GetData()[0]; + if (int_data.HasAttribute("available")) { + try { + auto available = std::any_cast(int_data.GetAttribute("available")); + std::cout << " available : " << available << std::endl; + } 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; + } + } + + 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 << " HDF5 Version : " << fHdf5Version; + 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 dims(rank); + dataspace.getSimpleExtentDims(dims.data(), nullptr); + + // Calculate total number of elements + hsize_t numElements = 1; + for (int i = 0; i < rank; i++) { + numElements *= dims[i]; + } + + // Create PNXdata object + PNXdata data(datatype); + data.SetDimensions(dims); + + // Read data + std::vector buffer(numElements); + dataset.read(buffer.data(), H5::PredType::NATIVE_INT); + data.SetData(buffer); + + // Read attributes + ReadDatasetAttributes(dataset, data); + + // Store in map + fDataMap[actualPath] = data; + + dataset.close(); + + if (fPrintDebug) { + std::cout << "debug> Read integer dataset: " << actualPath + << " (dims: "; + for (size_t i = 0; i < dims.size(); i++) { + std::cout << dims[i]; + if (i < dims.size() - 1) std::cout << " x "; + } + std::cout << ")" << std::endl; + } + } catch (const H5::Exception& err) { + std::cerr << "Error reading integer dataset " << path << ": " + << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Read float dataset and store in data map +//============================================================================= +/** + * Reads a float dataset from the HDF5 file and stores it in the data map. + * Handles multi-dimensional datasets by flattening them into a vector. + * + * @param file HDF5 file object + * @param path Path to the dataset + * @throws H5::Exception if reading fails + */ +void nxH5::PNeXus::ReadFloatDataset(H5::H5File& file, const std::string& path) +{ + try { + // Find the actual path (case-insensitive) + std::string actualPath = findDatasetPath(file, path); + + // Open the dataset + H5::DataSet dataset = file.openDataSet(actualPath); + H5::DataType datatype = dataset.getDataType(); + H5::DataSpace dataspace = dataset.getSpace(); + + // Get dimensions + int rank = dataspace.getSimpleExtentNdims(); + std::vector dims(rank); + dataspace.getSimpleExtentDims(dims.data(), nullptr); + + // Calculate total number of elements + hsize_t numElements = 1; + for (int i = 0; i < rank; i++) { + numElements *= dims[i]; + } + + // Create PNXdata object + PNXdata data(datatype); + data.SetDimensions(dims); + + // Read data + std::vector buffer(numElements); + dataset.read(buffer.data(), H5::PredType::NATIVE_FLOAT); + data.SetData(buffer); + + // Read attributes + ReadDatasetAttributes(dataset, data); + + // Store in map + fDataMap[actualPath] = data; + + dataset.close(); + + if (fPrintDebug) { + std::cout << "debug> Read float dataset: " << actualPath + << " (dims: "; + for (size_t i = 0; i < dims.size(); i++) { + std::cout << dims[i]; + if (i < dims.size() - 1) std::cout << " x "; + } + std::cout << ")" << std::endl; + } + } catch (const H5::Exception& err) { + std::cerr << "Error reading float dataset " << path << ": " + << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Read string dataset and store in data map +//============================================================================= +/** + * Reads a string dataset from the HDF5 file and stores it in the data map. + * + * @param file HDF5 file object + * @param path Path to the dataset + * @throws H5::Exception if reading fails + */ +void nxH5::PNeXus::ReadStringDataset(H5::H5File& file, const std::string& path) +{ + try { + // Find the actual path (case-insensitive) + std::string actualPath = findDatasetPath(file, path); + + // Open the dataset + H5::DataSet dataset = file.openDataSet(actualPath); + H5::DataType datatype = dataset.getDataType(); + H5::DataSpace dataspace = dataset.getSpace(); + + // Get dimensions + int rank = dataspace.getSimpleExtentNdims(); + std::vector dims(rank); + dataspace.getSimpleExtentDims(dims.data(), nullptr); + + // Calculate total number of elements + hsize_t numElements = 1; + for (int i = 0; i < rank; i++) { + numElements *= dims[i]; + } + + // Create PNXdata object + PNXdata data(datatype); + data.SetDimensions(dims); + + // Read data + std::vector buffer(numElements); + if (numElements == 1) { + // Single string + std::string value; + dataset.read(value, datatype); + buffer[0] = value; + } else { + // Multiple strings (if needed) + for (hsize_t i = 0; i < numElements; i++) { + dataset.read(buffer[i], datatype); + } + } + data.SetData(buffer); + + // Read attributes + ReadDatasetAttributes(dataset, data); + + // Store in map + fDataMap[actualPath] = data; + + dataset.close(); + + if (fPrintDebug) { + std::cout << "debug> Read string dataset: " << actualPath << std::endl; + } + } catch (const H5::Exception& err) { + std::cerr << "Error reading string dataset " << path << ": " + << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Read dataset attributes (template implementation) +//============================================================================= +/** + * Reads all attributes from an HDF5 dataset and adds them to a PNXdata object. + * Attributes can be integers, floats, or strings. + * + * @tparam T The data type of the PNXdata object + * @param dataset HDF5 dataset object + * @param data PNXdata object to add attributes to + */ +template +void nxH5::PNeXus::ReadDatasetAttributes(H5::DataSet& dataset, PNXdata& data) +{ + int numAttrs = dataset.getNumAttrs(); + + for (int i = 0; i < numAttrs; i++) { + H5::Attribute attr = dataset.openAttribute(i); + std::string attrName = attr.getName(); + H5::DataType attrType = attr.getDataType(); + H5T_class_t typeClass = attrType.getClass(); + + try { + if (typeClass == H5T_INTEGER) { + int32_t value; + attr.read(H5::PredType::NATIVE_INT32, &value); + data.AddAttribute(attrName, value); + } else if (typeClass == H5T_FLOAT) { + float value; + attr.read(H5::PredType::NATIVE_FLOAT, &value); + data.AddAttribute(attrName, value); + } else if (typeClass == H5T_STRING) { + std::string value; + attr.read(attrType, value); + data.AddAttribute(attrName, value); + } + } catch (const H5::Exception& err) { + if (fPrintDebug) { + std::cerr << "Warning: Could not read attribute " << attrName + << ": " << err.getDetailMsg() << std::endl; + } + } + + attr.close(); + } +} + +// Explicit template instantiations +template void nxH5::PNeXus::ReadDatasetAttributes(H5::DataSet&, PNXdata&); +template void nxH5::PNeXus::ReadDatasetAttributes(H5::DataSet&, PNXdata&); +template void nxH5::PNeXus::ReadDatasetAttributes(H5::DataSet&, PNXdata&); + +//============================================================================= +// WRITE METHODS +//============================================================================= + +//============================================================================= +// Create group hierarchy for a given path +//============================================================================= +/** + * Create nested group hierarchy for a given path. Groups are created recursively + * if they don't exist. If a group already exists, it is opened and returned. + * + * @param file HDF5 file object + * @param path Full path (e.g., "/raw_data_1/detector_1") + * @return H5::Group handle to the deepest created/opened group + * @throws H5::GroupIException if group creation fails + */ +H5::Group nxH5::PNeXus::CreateGroupHierarchy(H5::H5File& file, const std::string& path) +{ + std::vector components = splitPath(path); + std::string currentPath = ""; + H5::Group currentGroup; + + for (const auto& component : components) { + if (component.empty()) continue; // Skip empty (root indicator) + + currentPath += "/" + component; + + // Check if group exists + bool exists = false; + try { + H5::Group testGroup = file.openGroup(currentPath); + testGroup.close(); + exists = true; + } catch (const H5::Exception&) { + exists = false; + } + + // Create if doesn't exist + if (!exists) { + if (fPrintDebug) { + std::cout << "debug> Creating group: " << currentPath << std::endl; + } + currentGroup = file.createGroup(currentPath); + + // Write group attributes if any exist for this path + auto attrIt = fGroupAttributes.find(currentPath); + if (attrIt != fGroupAttributes.end()) { + if (fPrintDebug) { + std::cout << "debug> Writing " << attrIt->second.size() + << " attributes to group: " << currentPath << std::endl; + } + WriteGroupAttributes(currentGroup, attrIt->second); + } + + currentGroup.close(); + + if (fPrintDebug) { + std::cout << "debug> Created group: " << currentPath << std::endl; + } + } + } + + // Return an invalid group handle (we just needed to create the hierarchy) + return currentGroup; +} + +//============================================================================= +// Write dataset attributes from PNXdata object +//============================================================================= +/** + * Writes all attributes from a PNXdata object to an HDF5 dataset. + * Supports int32_t, float, and string attribute types. + * + * @tparam T The data type of the PNXdata object + * @param dataset HDF5 dataset object to write attributes to + * @param data PNXdata object containing attributes + */ +template +void nxH5::PNeXus::WriteDatasetAttributes(H5::DataSet& dataset, const PNXdata& data) +{ + const auto& attributes = data.GetAttributes(); + + for (const auto& [attrName, attrValue] : attributes) { + try { + // Try int32_t + if (auto* intVal = std::any_cast(&attrValue)) { + H5::DataSpace attrSpace(H5S_SCALAR); + H5::Attribute attr = dataset.createAttribute( + attrName, H5::PredType::NATIVE_INT32, attrSpace + ); + attr.write(H5::PredType::NATIVE_INT32, intVal); + attr.close(); + attrSpace.close(); + } + // Try float + else if (auto* floatVal = std::any_cast(&attrValue)) { + H5::DataSpace attrSpace(H5S_SCALAR); + H5::Attribute attr = dataset.createAttribute( + attrName, H5::PredType::NATIVE_FLOAT, attrSpace + ); + attr.write(H5::PredType::NATIVE_FLOAT, floatVal); + attr.close(); + attrSpace.close(); + } + // Try string + else if (auto* strVal = std::any_cast(&attrValue)) { + H5::StrType strType(H5::PredType::C_S1, H5T_VARIABLE); + H5::DataSpace attrSpace(H5S_SCALAR); + H5::Attribute attr = dataset.createAttribute( + attrName, strType, attrSpace + ); + attr.write(strType, *strVal); + attr.close(); + attrSpace.close(); + } + else { + if (fPrintDebug) { + std::cerr << "Warning: Unsupported attribute type for " + << attrName << std::endl; + } + } + + } catch (const H5::Exception& err) { + if (fPrintDebug) { + std::cerr << "Warning: Could not write attribute " << attrName + << ": " << err.getDetailMsg() << std::endl; + } + } + } +} + +//============================================================================= +// Write integer dataset with attributes +//============================================================================= +/** + * Writes an integer dataset to the HDF5 file with all its attributes. + * The parent group hierarchy is created automatically if it doesn't exist. + * + * @param file HDF5 file object + * @param path Dataset path (groups created automatically) + * @param data PNXdata object containing data, dimensions, and attributes + * @throws H5::Exception if writing fails + */ +void nxH5::PNeXus::WriteIntDataset(H5::H5File& file, const std::string& path, + const PNXdata& data) +{ + try { + // Extract parent path and dataset name + size_t lastSlash = path.find_last_of('/'); + std::string parentPath = (lastSlash == 0) ? "/" : path.substr(0, lastSlash); + std::string datasetName = path.substr(lastSlash + 1); + + // Create parent group hierarchy + if (!parentPath.empty() && parentPath != "/") { + CreateGroupHierarchy(file, parentPath); + } + + // Get dimensions and data + const auto& dims = data.GetDimensions(); + const auto& buffer = data.GetData(); + + // Validate data consistency + if (buffer.size() != data.GetNumElements()) { + throw std::runtime_error("Data size mismatch with dimensions"); + } + + // Create dataspace + H5::DataSpace dataspace(dims.size(), dims.data()); + + // Create dataset + H5::DataSet dataset = file.createDataSet( + path, + H5::PredType::NATIVE_INT, + dataspace + ); + + // Write data + dataset.write(buffer.data(), H5::PredType::NATIVE_INT); + + // Write attributes + WriteDatasetAttributes(dataset, data); + + // Close resources + dataset.close(); + dataspace.close(); + + if (fPrintDebug) { + std::cout << "debug> Wrote integer dataset: " << path + << " (dims: "; + for (size_t i = 0; i < dims.size(); i++) { + std::cout << dims[i]; + if (i < dims.size() - 1) std::cout << " x "; + } + std::cout << ")" << std::endl; + } + + } catch (const H5::Exception& err) { + std::cerr << "Error writing integer dataset " << path << ": " + << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Write float dataset with attributes +//============================================================================= +/** + * Writes a float dataset to the HDF5 file with all its attributes. + * The parent group hierarchy is created automatically if it doesn't exist. + * + * @param file HDF5 file object + * @param path Dataset path (groups created automatically) + * @param data PNXdata object containing data, dimensions, and attributes + * @throws H5::Exception if writing fails + */ +void nxH5::PNeXus::WriteFloatDataset(H5::H5File& file, const std::string& path, + const PNXdata& data) +{ + try { + // Extract parent path and dataset name + size_t lastSlash = path.find_last_of('/'); + std::string parentPath = (lastSlash == 0) ? "/" : path.substr(0, lastSlash); + std::string datasetName = path.substr(lastSlash + 1); + + // Create parent group hierarchy + if (!parentPath.empty() && parentPath != "/") { + CreateGroupHierarchy(file, parentPath); + } + + // Get dimensions and data + const auto& dims = data.GetDimensions(); + const auto& buffer = data.GetData(); + + // Validate data consistency + if (buffer.size() != data.GetNumElements()) { + throw std::runtime_error("Data size mismatch with dimensions"); + } + + // Create dataspace + H5::DataSpace dataspace(dims.size(), dims.data()); + + // Create dataset + H5::DataSet dataset = file.createDataSet( + path, + H5::PredType::NATIVE_FLOAT, + dataspace + ); + + // Write data + dataset.write(buffer.data(), H5::PredType::NATIVE_FLOAT); + + // Write attributes + WriteDatasetAttributes(dataset, data); + + // Close resources + dataset.close(); + dataspace.close(); + + if (fPrintDebug) { + std::cout << "debug> Wrote float dataset: " << path + << " (dims: "; + for (size_t i = 0; i < dims.size(); i++) { + std::cout << dims[i]; + if (i < dims.size() - 1) std::cout << " x "; + } + std::cout << ")" << std::endl; + } + + } catch (const H5::Exception& err) { + std::cerr << "Error writing float dataset " << path << ": " + << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Write string dataset with attributes +//============================================================================= +/** + * Writes a string dataset to the HDF5 file with all its attributes. + * The parent group hierarchy is created automatically if it doesn't exist. + * Handles both scalar strings and arrays of strings. + * + * @param file HDF5 file object + * @param path Dataset path (groups created automatically) + * @param data PNXdata object containing data, dimensions, and attributes + * @throws H5::Exception if writing fails + */ +void nxH5::PNeXus::WriteStringDataset(H5::H5File& file, const std::string& path, + const PNXdata& data) +{ + try { + // Extract parent path and dataset name + size_t lastSlash = path.find_last_of('/'); + std::string parentPath = (lastSlash == 0) ? "/" : path.substr(0, lastSlash); + std::string datasetName = path.substr(lastSlash + 1); + + // Create parent group hierarchy + if (!parentPath.empty() && parentPath != "/") { + CreateGroupHierarchy(file, parentPath); + } + + // Get data + const auto& buffer = data.GetData(); + const auto& dims = data.GetDimensions(); + + /* //as35 + // Handle single string vs array of strings + if (buffer.size() == 1 && dims.size() == 1 && dims[0] == 1) { + // Single string - use variable-length string type + H5::StrType strType(H5::PredType::C_S1, H5T_VARIABLE); + H5::DataSpace dataspace(H5S_SCALAR); + + H5::DataSet dataset = file.createDataSet(path, strType, dataspace); + dataset.write(buffer[0], strType); + + WriteDatasetAttributes(dataset, data); + dataset.close(); + dataspace.close(); + + } else { + */ //as35 + // Array of strings + H5::StrType strType(H5::PredType::C_S1, H5T_VARIABLE); + H5::DataSpace dataspace(dims.size(), dims.data()); + + H5::DataSet dataset = file.createDataSet(path, strType, dataspace); + + // For arrays, need to create array of C strings + std::vector cStrings; + for (const auto& str : buffer) { + cStrings.push_back(str.c_str()); + } + dataset.write(cStrings.data(), strType); + + WriteDatasetAttributes(dataset, data); + dataset.close(); + dataspace.close(); + /* //as35 + } + */ //as35 + + if (fPrintDebug) { + std::cout << "debug> Wrote string dataset: " << path << std::endl; + } + + } catch (const H5::Exception& err) { + std::cerr << "Error writing string dataset " << path << ": " + << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Write group attributes +//============================================================================= +/** + * Writes attributes to an HDF5 group object. Supports int32_t, float, and + * string attribute types. + * + * @param group HDF5 group object to write attributes to + * @param attributes Map of attribute names to values (stored as std::any) + */ +void nxH5::PNeXus::WriteGroupAttributes(H5::Group& group, const std::map& attributes) +{ + for (const auto& [attrName, attrValue] : attributes) { + try { + // Try int32_t + if (auto* intVal = std::any_cast(&attrValue)) { + H5::DataSpace attrSpace(H5S_SCALAR); + H5::Attribute attr = group.createAttribute( + attrName, H5::PredType::NATIVE_INT32, attrSpace + ); + attr.write(H5::PredType::NATIVE_INT32, intVal); + attr.close(); + attrSpace.close(); + } + // Try float + else if (auto* floatVal = std::any_cast(&attrValue)) { + H5::DataSpace attrSpace(H5S_SCALAR); + H5::Attribute attr = group.createAttribute( + attrName, H5::PredType::NATIVE_FLOAT, attrSpace + ); + attr.write(H5::PredType::NATIVE_FLOAT, floatVal); + attr.close(); + attrSpace.close(); + } + // Try string + else if (auto* strVal = std::any_cast(&attrValue)) { + H5::StrType strType(H5::PredType::C_S1, H5T_VARIABLE); + H5::DataSpace attrSpace(H5S_SCALAR); + H5::Attribute attr = group.createAttribute( + attrName, strType, attrSpace + ); + attr.write(strType, *strVal); + attr.close(); + attrSpace.close(); + } + else { + if (fPrintDebug) { + std::cerr << "Warning: Unsupported attribute type for " + << attrName << std::endl; + } + } + + } catch (const H5::Exception& err) { + if (fPrintDebug) { + std::cerr << "Warning: Failed to write group attribute " + << attrName << ": " << err.getDetailMsg() << std::endl; + } + } + } +} + +//============================================================================= +// Write root-level file attributes +//============================================================================= +/** + * Writes root-level file attributes such as NeXus_version, HDF5_version, + * file_name, and file_time to the HDF5 file. + * + * @param file HDF5 file object + * @throws H5::Exception if writing fails + */ +void nxH5::PNeXus::WriteFileAttributes(H5::H5File& file) +{ + try { + H5::StrType strType(H5::PredType::C_S1, H5T_VARIABLE); + H5::DataSpace attrSpace(H5S_SCALAR); + + // Check if custom root-level attributes exist + auto rootAttrIt = fGroupAttributes.find("/"); + bool hasCustomAttrs = (rootAttrIt != fGroupAttributes.end()); + + // NeXus_version attribute + if (!fNeXusVersion.empty() && + (!hasCustomAttrs || rootAttrIt->second.find("NeXus_version") == rootAttrIt->second.end())) { + H5::Attribute attr = file.createAttribute( + "NeXus_version", strType, attrSpace + ); + attr.write(strType, fNeXusVersion); + attr.close(); + } + + // HDF5_version attribute + if (!fHdf5Version.empty() && + (!hasCustomAttrs || rootAttrIt->second.find("HDF5_version") == rootAttrIt->second.end())) { + H5::Attribute attr = file.createAttribute( + "HDF5_version", strType, attrSpace + ); + attr.write(strType, fHdf5Version); + attr.close(); + } + + // file_name attribute + if (!fFileNameNxs.empty() && + (!hasCustomAttrs || rootAttrIt->second.find("file_name") == rootAttrIt->second.end())) { + H5::Attribute attr = file.createAttribute( + "file_name", strType, attrSpace + ); + attr.write(strType, fFileNameNxs); + attr.close(); + } + + // file_time attribute + if (!fFileTimeNxs.empty() && + (!hasCustomAttrs || rootAttrIt->second.find("file_time") == rootAttrIt->second.end())) { + H5::Attribute attr = file.createAttribute( + "file_time", strType, attrSpace + ); + attr.write(strType, fFileTimeNxs); + attr.close(); + } + + attrSpace.close(); + + // Write any custom root-level attributes from fGroupAttributes["/"] + if (hasCustomAttrs) { + if (fPrintDebug) { + std::cout << "debug> Writing " << rootAttrIt->second.size() + << " custom attributes to root level" << std::endl; + } + + for (const auto& [attrName, attrValue] : rootAttrIt->second) { + try { + // Try int32_t + if (auto* intVal = std::any_cast(&attrValue)) { + H5::DataSpace scalarSpace(H5S_SCALAR); + H5::Attribute attr = file.createAttribute( + attrName, H5::PredType::NATIVE_INT32, scalarSpace + ); + attr.write(H5::PredType::NATIVE_INT32, intVal); + attr.close(); + scalarSpace.close(); + } + // Try float + else if (auto* floatVal = std::any_cast(&attrValue)) { + H5::DataSpace scalarSpace(H5S_SCALAR); + H5::Attribute attr = file.createAttribute( + attrName, H5::PredType::NATIVE_FLOAT, scalarSpace + ); + attr.write(H5::PredType::NATIVE_FLOAT, floatVal); + attr.close(); + scalarSpace.close(); + } + // Try string + else if (auto* strVal = std::any_cast(&attrValue)) { + H5::StrType varStrType(H5::PredType::C_S1, H5T_VARIABLE); + H5::DataSpace scalarSpace(H5S_SCALAR); + H5::Attribute attr = file.createAttribute( + attrName, varStrType, scalarSpace + ); + attr.write(varStrType, *strVal); + attr.close(); + scalarSpace.close(); + } + else { + if (fPrintDebug) { + std::cerr << "Warning: Unsupported attribute type for root-level " + << attrName << std::endl; + } + } + } catch (const H5::Exception& err) { + if (fPrintDebug) { + std::cerr << "Warning: Failed to write root-level attribute " + << attrName << ": " << err.getDetailMsg() << std::endl; + } + } + } + } + + } catch (const H5::Exception& err) { + std::cerr << "Error writing file attributes: " + << err.getDetailMsg() << std::endl; + throw; + } +} + +//============================================================================= +// Write IDF version 2 structure +//============================================================================= +/** + * Writes all datasets from the data map to the HDF5 file following the + * IDF version 2 structure. Iterates through fDataMap and writes each + * dataset based on its type. + * + * @param file HDF5 file object + * @throws H5::Exception if writing fails + */ +void nxH5::PNeXus::WriteIdfV2(H5::H5File& file) +{ + if (fPrintDebug) { + std::cout << "debug> Writing IDF v2 structure..." << std::endl; + } + + // Iterate through all datasets in fDataMap + for (const auto& [path, anyData] : fDataMap) { + try { + // Try PNXdata + try { + auto intData = std::any_cast>(anyData); + WriteIntDataset(file, path, intData); + continue; + } catch (const std::bad_any_cast&) {} + + // Try PNXdata - convert to int + try { + auto int32Data = std::any_cast>(anyData); + + // Convert PNXdata to PNXdata + PNXdata intData(int32Data.GetDataType()); + intData.SetDimensions(int32Data.GetDimensions()); + + std::vector convertedData; + for (auto val : int32Data.GetData()) { + convertedData.push_back(static_cast(val)); + } + intData.SetData(convertedData); + + // Copy attributes + for (const auto& [name, value] : int32Data.GetAttributes()) { + intData.AddAttribute(name, value); + } + + WriteIntDataset(file, path, intData); + continue; + } catch (const std::bad_any_cast&) {} + + // Try PNXdata + try { + auto floatData = std::any_cast>(anyData); + WriteFloatDataset(file, path, floatData); + continue; + } catch (const std::bad_any_cast&) {} + + // Try PNXdata + try { + auto strData = std::any_cast>(anyData); + WriteStringDataset(file, path, strData); + continue; + } catch (const std::bad_any_cast&) {} + + // Unknown type + if (fPrintDebug) { + std::cerr << "Warning: Unknown data type for path " << path + << std::endl; + } + + } catch (const H5::Exception& err) { + std::cerr << "Error writing dataset " << path << ": " + << err.getDetailMsg() << std::endl; + throw; + } + } +} + +//============================================================================= +// Group attribute management methods +//============================================================================= + +/** + * Add or update an attribute for a group. + * + * @param groupPath HDF5 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 + */ +bool nxH5::PNeXus::AddGroupAttribute(const std::string& groupPath, const std::string& attrName, + const std::any& attrValue) +{ + fGroupAttributes[groupPath][attrName] = attrValue; + return true; +} + +/** + * Remove an attribute from a group. + * + * @param groupPath HDF5 path of the group + * @param attrName Attribute name to remove + * @return true if attribute was removed, false if group or attribute not found + */ +bool nxH5::PNeXus::RemoveGroupAttribute(const std::string& groupPath, const std::string& attrName) +{ + auto groupIt = fGroupAttributes.find(groupPath); + if (groupIt == fGroupAttributes.end()) { + return false; + } + + auto attrIt = groupIt->second.find(attrName); + if (attrIt == groupIt->second.end()) { + return false; + } + + groupIt->second.erase(attrIt); + + // Clean up empty group entry + if (groupIt->second.empty()) { + fGroupAttributes.erase(groupIt); + } + + return true; +} + +/** + * Check if a group has a specific attribute. + * + * @param groupPath HDF5 path of the group + * @param attrName Attribute name to check + * @return true if attribute exists, false otherwise + */ +bool nxH5::PNeXus::HasGroupAttribute(const std::string& groupPath, const std::string& attrName) const +{ + auto groupIt = fGroupAttributes.find(groupPath); + if (groupIt == fGroupAttributes.end()) { + return false; + } + + return groupIt->second.find(attrName) != groupIt->second.end(); +} + +/** + * Get an attribute value from a group. + * + * @param groupPath HDF5 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 nxH5::PNeXus::GetGroupAttribute(const std::string& groupPath, const std::string& attrName) const +{ + return fGroupAttributes.at(groupPath).at(attrName); +} + +/** + * Get all attributes for a group. + * + * @param groupPath HDF5 path of the group + * @return Map of attribute names to values + */ +const std::map& nxH5::PNeXus::GetGroupAttributes(const std::string& groupPath) const +{ + static const std::map emptyMap; + auto it = fGroupAttributes.find(groupPath); + if (it == fGroupAttributes.end()) { + return emptyMap; + } + return it->second; +} + +/** + * Clear all attributes from a group. + * + * @param groupPath HDF5 path of the group + * @return true if group was found and attributes cleared, false otherwise + */ +bool nxH5::PNeXus::ClearGroupAttributes(const std::string& groupPath) +{ + auto it = fGroupAttributes.find(groupPath); + if (it == fGroupAttributes.end()) { + return false; + } + + fGroupAttributes.erase(it); + return true; +} + +/** + * Add or update an attribute at the root level "/". + * This is a convenience method that calls AddGroupAttribute("/", attrName, attrValue). + * + * @param attrName Attribute name + * @param attrValue Attribute value (stored as std::any) + * @return true if attribute was added successfully + */ +bool nxH5::PNeXus::AddRootAttribute(const std::string& attrName, const std::any& attrValue) +{ + return AddGroupAttribute("/", attrName, attrValue); +} + +//============================================================================= +// Write NeXus file (main entry point) +//============================================================================= +/** + * Writes the data map contents to a new NeXus HDF5 file. + * This is the main public API for writing NeXus files. + * + * @param filename Path to the output NeXus HDF5 file + * @param idfVersion IDF version to write (default: 2) + * @return 0 on success, 1 on error + * @throws H5::FileIException if file cannot be created + * @throws H5::GroupIException if group creation fails + * @throws H5::DataSetIException if dataset writing fails + */ +int nxH5::PNeXus::WriteNexusFile(const std::string& filename, int idfVersion) +{ + try { + if (fPrintDebug) { + std::cout << std::endl; + std::cout << "debug> Creating NeXus file: " << filename << std::endl; + } + + // Validate that we have data to write + if (fDataMap.empty()) { + std::cerr << "Error: No data to write (data map is empty)" << std::endl; + return 1; + } + + // Create new HDF5 file (truncate if exists) + H5::H5File file(filename, H5F_ACC_TRUNC); + + // Turn off auto-printing for exceptions + H5::Exception::dontPrint(); + + // Write root-level file attributes + WriteFileAttributes(file); + + // Write structure based on IDF version + if (idfVersion == 1) { + std::cerr << "Error: IDF version 1 writing not yet implemented" + << std::endl; + file.close(); + return 1; + } else if (idfVersion == 2) { + WriteIdfV2(file); + } else { + std::cerr << "Error: Unsupported IDF version " << idfVersion + << std::endl; + file.close(); + return 1; + } + + // Close file + file.close(); + + if (fPrintDebug) { + std::cout << "debug> Successfully wrote NeXus file: " << filename + << std::endl; + } + + return 0; + + } catch (const H5::FileIException& err) { + std::cerr << "Error: Failed to create file '" << filename << "': " + << err.getDetailMsg() << std::endl; + return 1; + } catch (const H5::Exception& err) { + std::cerr << "Error: HDF5 exception occurred: " + << err.getDetailMsg() << std::endl; + return 1; + } catch (const std::exception& err) { + std::cerr << "Error: Unexpected exception: " << err.what() << std::endl; + return 1; + } +} + +// Explicit template instantiations for write methods +template void nxH5::PNeXus::WriteDatasetAttributes(H5::DataSet&, const PNXdata&); +template void nxH5::PNeXus::WriteDatasetAttributes(H5::DataSet&, const PNXdata&); +template void nxH5::PNeXus::WriteDatasetAttributes(H5::DataSet&, const PNXdata&); diff --git a/src/external/nexus/PNeXus.h b/src/external/nexus/PNeXus.h index b3212da1..3de538fb 100644 --- a/src/external/nexus/PNeXus.h +++ b/src/external/nexus/PNeXus.h @@ -45,6 +45,8 @@ #include #endif // HAVE_HDF4 +#include + namespace nxs { enum class HDFType { @@ -892,4 +894,837 @@ private: } #endif // HAVE_HDF4 +namespace nxH5 { + +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 + * nxH5::PNeXus nexus("file.nxs"); + * nxH5::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 nxH5::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; } + + std::vector GetDeadTimeEstimated() { return fDeadTimeEstimated; } + +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) from the file + std::vector fDeadTimeEstimated; ///< Dead time values per detector (microseconds) as estimated from the data + 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 +}; + + +/** + * @class PNXdata + * @brief Template class for storing HDF5 dataset content with attributes + * + * The PNXdata class stores data read from an HDF5 dataset along with metadata + * such as dimensions, datatype, and attributes. It is designed to be stored in + * std::map where the key is the HDF5 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 HDF5 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(H5::PredType::NATIVE_INT) {} + + /** + * @brief Constructor with datatype + * @param dataType HDF5 datatype for this dataset + */ + PNXdata(const H5::DataType& dataType) : fDataType(dataType) {} + + /** + * @brief Get the HDF5 DataType + * @return The HDF5 DataType for this dataset + */ + H5::DataType GetDataType() const { return fDataType; } + + /** + * @brief Set the HDF5 DataType + * @param dataType The HDF5 DataType to set + */ + void SetDataType(const H5::DataType& 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: + H5::DataType fDataType; ///< HDF5 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 HDF5 file reader with case-insensitive path lookup + * + * The PNeXus class provides functionality for reading ISIS muon NeXus HDF5 files + * using the HDF5 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 with HDF5 C++ exceptions + * + * @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 HDF5 C++ exceptions on failure + */ +class PNeXus { +public: + PNeXus(); + + /** + * @brief Constructor - creates PNeXus object and reads the NeXus file + * @param fln Filename of the NeXus HDF5 file to read + * @throws H5::FileIException if file cannot be opened + * @throws H5::AttributeIException if required attributes are missing + * @throws H5::DataSetIException if required datasets are missing + */ + PNeXus(const std::string fln, const bool printDebug=false); + + /** + * @brief Get the filename of the NeXus file + * @return The filename as a string + */ + std::string GetFileName() const { return fFileName; } + + /** + * @brief Get the hdf5 version of the NeXus file + * @return The hdf5 version as a string + */ + std::string GetHdf5Version() const { return fHdf5Version; } + + /** + * @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 HDF5 file + * @return 0 on success, 1 on error + * @throws H5::FileIException if file cannot be opened + * @throws H5::AttributeIException if required attributes are missing + * @throws H5::DataSetIException if required datasets are missing + */ + 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 HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 path for the new dataset + * @param data Data vector + * @param dimensions Dimensions vector + * @param dataType HDF5 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 H5::DataType& dataType = H5::PredType::NATIVE_INT) { + 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 HDF5 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 hdf5-NeXus file content which was read + */ + void Dump(); + + /** + * @brief Write the data map contents to a NeXus HDF5 file + * @param filename Path to the output NeXus HDF5 file + * @param idfVersion IDF version to write (default: 2) + * @return 0 on success, 1 on error + * @throws H5::FileIException if file cannot be created + * @throws H5::GroupIException if group creation fails + * @throws H5::DataSetIException if dataset writing fails + */ + int WriteNexusFile(const std::string& filename, int idfVersion = 2); + + /** + * @brief Add or update an attribute for a group + * @param groupPath HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 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 HDF5 filename + int fIdfVersion{-1}; ///< IDF version of the NeXus file + std::string fHdf5Version{""}; ///< HDF5 version of the file + std::string fNeXusVersion{""}; ///< NeXus version of the file + std::string fFileNameNxs{""}; + std::string fFileTimeNxs{""}; + std::string fCreatorNxs{""}; + std::string fUserV1{""}; + std::map fDataMap; ///< Map of HDF5 paths to PNXdata objects + std::map> fGroupAttributes; ///< Map of group paths to their attributes + + /** + * @brief HandleIdfV1 + * @param file + */ + void HandleIdfV1(H5::H5File &file); + + /** + * @brief HandleIdfV2 + * @param file + */ + void HandleIdfV2(H5::H5File &file); + + // ======================================================================== + // Write methods for HDF5 file creation + // ======================================================================== + + /** + * @brief Create nested group hierarchy for a given path + * @param file HDF5 file object + * @param path Full path (e.g., "/raw_data_1/detector_1") + * @return H5::Group handle to the deepest created/opened group + * @throws H5::GroupIException if group creation fails + */ + H5::Group CreateGroupHierarchy(H5::H5File& file, const std::string& path); + + /** + * @brief Write dataset attributes from PNXdata object + * @tparam T The data type of the PNXdata object + * @param dataset HDF5 dataset object to write attributes to + * @param data PNXdata object containing attributes + */ + template + void WriteDatasetAttributes(H5::DataSet& dataset, const PNXdata& data); + + /** + * @brief Write attributes to a group + * @param group HDF5 group object to write attributes to + * @param attributes Map of attribute names to values + */ + void WriteGroupAttributes(H5::Group& group, const std::map& attributes); + + /** + * @brief Write an integer dataset with attributes + * @param file HDF5 file object + * @param path Dataset path (groups created automatically) + * @param data PNXdata object containing data, dimensions, and attributes + * @throws H5::Exception if writing fails + */ + void WriteIntDataset(H5::H5File& file, const std::string& path, + const PNXdata& data); + + /** + * @brief Write a float dataset with attributes + * @param file HDF5 file object + * @param path Dataset path (groups created automatically) + * @param data PNXdata object containing data, dimensions, and attributes + * @throws H5::Exception if writing fails + */ + void WriteFloatDataset(H5::H5File& file, const std::string& path, + const PNXdata& data); + + /** + * @brief Write a string dataset with attributes + * @param file HDF5 file object + * @param path Dataset path (groups created automatically) + * @param data PNXdata object containing data, dimensions, and attributes + * @throws H5::Exception if writing fails + */ + void WriteStringDataset(H5::H5File& file, const std::string& path, + const PNXdata& data); + + /** + * @brief Write root-level file attributes + * @param file HDF5 file object + * @throws H5::Exception if writing fails + */ + void WriteFileAttributes(H5::H5File& file); + + /** + * @brief Write IDF version 2 structure + * @param file HDF5 file object + * @throws H5::Exception if writing fails + */ + void WriteIdfV2(H5::H5File& file); + + // ======================================================================== + // 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 HDF5 path into components + * @param path HDF5 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 obj HDF5 File object to search + * @param requestedName Requested attribute name (any case) + * @return Correctly-cased attribute name as it exists in the file + * @throws H5::AttributeIException if attribute not found + * @example "nexus_version" might resolve to "NeXus_version" + */ + std::string findAttributeName(H5::H5File& obj, const std::string& requestedName); + + /** + * @brief Find group path with case-insensitive matching (File version) + * @param parent HDF5 File object (root) + * @param requestedPath Requested group path (any case) + * @return Correctly-cased group path as it exists in the file + * @throws H5::GroupIException if group not found or path invalid + * @example "/RAW_DATA_1" might resolve to "/raw_data_1" + */ + std::string findGroupPath(H5::H5File& parent, const std::string& requestedPath); + + /** + * @brief Find group path with case-insensitive matching (Group version) + * @param parent HDF5 Group object to start search from + * @param requestedPath Requested group path (any case) + * @return Correctly-cased group path as it exists in the file + * @throws H5::GroupIException if group not found or path invalid + */ + std::string findGroupPath(H5::Group& parent, const std::string& requestedPath); + + /** + * @brief Find dataset path with case-insensitive matching (File version) + * @param parent HDF5 File object (root) + * @param requestedPath Requested dataset path (any case) + * @return Correctly-cased dataset path as it exists in the file + * @throws H5::DataSetIException if dataset not found or path invalid + * @example "/RAW_DATA_1/idf_VERSION" might resolve to "/raw_data_1/IDF_version" + * @note Validates that final path component is actually a dataset + */ + std::string findDatasetPath(H5::H5File& parent, const std::string& requestedPath); + + /** + * @brief Find dataset path with case-insensitive matching (Group version) + * @param parent HDF5 Group object to start search from + * @param requestedPath Requested dataset path (any case) + * @return Correctly-cased dataset path as it exists in the file + * @throws H5::DataSetIException if dataset not found or path invalid + * @note Validates that final path component is actually a dataset + */ + std::string findDatasetPath(H5::Group& parent, const std::string& requestedPath); + + // ======================================================================== + // Dataset reading helper methods + // ======================================================================== + + /** + * @brief Read an integer dataset and store in data map + * @param file HDF5 file object + * @param path Path to the dataset + * @throws H5::Exception if reading fails + */ + void ReadIntDataset(H5::H5File& file, const std::string& path); + + /** + * @brief Read a float dataset and store in data map + * @param file HDF5 file object + * @param path Path to the dataset + * @throws H5::Exception if reading fails + */ + void ReadFloatDataset(H5::H5File& file, const std::string& path); + + /** + * @brief Read a string dataset and store in data map + * @param file HDF5 file object + * @param path Path to the dataset + * @throws H5::Exception if reading fails + */ + void ReadStringDataset(H5::H5File& file, const std::string& path); + + /** + * @brief Read dataset attributes and add to PNXdata object + * @tparam T The data type of the PNXdata object + * @param dataset HDF5 dataset object + * @param data PNXdata object to add attributes to + */ + template + void ReadDatasetAttributes(H5::DataSet& dataset, PNXdata& data); +}; + +} + #endif // _PNEXUS_H_