/***************************************************************************

  PNeXus.cpp

  Author: Andreas Suter
  e-mail: andreas.suter@psi.ch

***************************************************************************/

/***************************************************************************
 *   Copyright (C) 2007-2026 by Andreas Suter                              *
 *   andreas.suter@psi.ch                                                  *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include <iostream>
#include <fstream>
#include <sstream>
#include <stdexcept>
#include <cstring>
#include <cmath>

#include "Minuit2/MnMinimize.h"
#include "Minuit2/FunctionMinimum.h"

#include "PNeXus.h"

nxs::HDFType nxs::checkHDFType(const std::string& filename) {

  std::ifstream file(filename, std::ios::binary);
  if (!file.is_open()) {
    std::cerr << "Error: Cannot open file '" << filename << "'" << std::endl;
    return nxs::HDFType::Unknown;
  }

  // Read first 8 bytes (enough for both signatures)
  unsigned char header[8];
  file.read(reinterpret_cast<char*>(header), 8);

  if (file.gcount() < 4) {
    std::cerr << "**Error**: File too small to be HDF4 or HDF5" << std::endl;
    return nxs::HDFType::Unknown;
  }

  // HDF5 signature: 0x89 'H' 'D' 'F' 0x0d 0x0a 0x1a 0x0a
  const unsigned char hdf5_sig[8] = {0x89, 0x48, 0x44, 0x46, 0x0d, 0x0a, 0x1a, 0x0a};

  // HDF4 signature: 0x0e 0x03 0x13 0x01
  const unsigned char hdf4_sig[4] = {0x0e, 0x03, 0x13, 0x01};

  // Check HDF5 first (needs 8 bytes)
  if (file.gcount() >= 8 && std::memcmp(header, hdf5_sig, 8) == 0) {
    return nxs::HDFType::HDF5;
  }

  // Check HDF4 (needs 4 bytes)
  if (std::memcmp(header, hdf4_sig, 4) == 0) {
    return nxs::HDFType::HDF4;
  }

  return nxs::HDFType::Unknown;
}

#ifdef HAVE_HDF4
//=============================================================================
// PNeXusDeadTime Constructor
//=============================================================================
nxH4::PNeXusDeadTime::PNeXusDeadTime(const nxH4::PNeXus *nxs, bool debug) : fDebug(debug)
{
  std::map<std::string, std::any> 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<PNXdata<int>>(dataMap["/raw_data_1/detector_1/counts"]);
      const auto& counts = counts_data.GetData();
      fCounts = counts;
      auto dims = counts_data.GetDimensions();
      if (fDebug) {
        std::cout << " counts dimensions: " << dims[0] << " x "
                  << dims[1] << " x " << dims[2] << std::endl;
      }
      fDims = dims;
    } else {
      std::cout << "**ERROR** Couldn't find 'counts' dataset" << std::endl;
      fValid = false;
      return;
    }

    // get time_resolution
    if (dataMap.find("/raw_data_1/detector_1/time_resolution") != dataMap.end()) {
      auto time_res_data = std::any_cast<PNXdata<float>>(dataMap["/raw_data_1/detector_1/time_resolution"]);
      const auto& time_res = time_res_data.GetData();
      if (time_res.size() > 0) {
        fTimeResolution = time_res[0];
        if (fDebug) {
          std::cout << " time resolution: " << fTimeResolution << " ps" << std::endl;
        }
      }
    } else {
      std::cout << "**ERROR** Couldn't find 'time_resolution' dataset" << std::endl;
      fValid = false;
      return;
    }

    // get good_frames
    if (dataMap.find("/raw_data_1/good_frames") != dataMap.end()) {
      auto good_frames_data = std::any_cast<PNXdata<int>>(dataMap["/raw_data_1/good_frames"]);
      const auto& good_frames = good_frames_data.GetData();
      if (good_frames.size() > 0) {
        fGoodFrames = good_frames[0];
        if (fDebug) {
          std::cout << " good frames: " << fGoodFrames << std::endl;
        }
      }
    } else {
      std::cout << "**ERROR** Couldn't find 'good_frames' dataset" << std::endl;
      fValid = false;
      return;
    }

    // get t0_bin, first_good_bin, last_good_bin
    if (dataMap.find("/raw_data_1/detector_1/spectrum_index") != dataMap.end()) {
      auto spec_idx_data = std::any_cast<PNXdata<int>>(dataMap["/raw_data_1/detector_1/spectrum_index"]);
      if (spec_idx_data.HasAttribute("t0_bin")) {
        try {
          fT0Bin = std::any_cast<int>(spec_idx_data.GetAttribute("t0_bin"));
        } catch (...) {
          std::cout << "**WARNING** Couldn't read t0_bin attribute" << std::endl;
        }
      }

      if (spec_idx_data.HasAttribute("first_good_bin")) {
        try {
          fFgbBin = std::any_cast<int>(spec_idx_data.GetAttribute("first_good_bin"));
        } catch (...) {
          std::cout << "**WARNING** Couldn't read first_good_bin attribute" << std::endl;
        }
      }

      if (spec_idx_data.HasAttribute("last_good_bin")) {
        try {
          fLgbBin = std::any_cast<int>(spec_idx_data.GetAttribute("last_good_bin"));
        } catch (...) {
          std::cout << "**WARNING** Couldn't read last_good_bin attribute" << std::endl;
        }
      }
    }

    if (fDebug) {
      std::cout << " t0_bin: " << fT0Bin << ", first_good_bin: " << fFgbBin
                << ", last_good_bin: " << fLgbBin << std::endl;
    }
  }

  fValid = true;
}

//=============================================================================
// PNeXusDeadTime::operator()
//=============================================================================
double nxH4::PNeXusDeadTime::operator()(const std::vector<double> &par) const
{
  double chisq = 0.0;
  double tau = par[0]; // dead time in microseconds

  // Convert time resolution from ps to microseconds
  double dt = fTimeResolution * 1.0e-6; // ps -> us

  // Calculate chi-square for the specified spectrum
  int period = 0; // assuming single period

  for (unsigned int bin = fFgbBin; bin <= static_cast<unsigned int>(fLgbBin); bin++) {
    int idx = period * fDims[1] * fDims[2] + fIdx * fDims[2] + bin;
    double n_obs = static_cast<double>(fCounts[idx]);

    if (n_obs > 0) {
      // Dead time correction: n_true = n_obs / (1 - n_obs * tau / (good_frames * dt))
      double correction = 1.0 - (n_obs * tau) / (fGoodFrames * dt);
      if (correction > 0) {
        double n_expected = n_obs / correction;
        double diff = n_obs - n_expected;
        chisq += (diff * diff) / n_expected;
      }
    }
  }

  return chisq;
}

//=============================================================================
// PNeXusDeadTime::minimize
//=============================================================================
void nxH4::PNeXusDeadTime::minimize(const int i)
{
  if (i < 0 || i >= static_cast<int>(fDims[1])) {
    std::cerr << "**ERROR** Invalid spectrum index: " << i << std::endl;
    return;
  }

  fIdx = static_cast<unsigned int>(i);

  ROOT::Minuit2::MnUserParameters upar;
  upar.Add("dead_time", 0.1, 0.01); // initial value, step size
  upar.SetLimits("dead_time", 0.0, 10.0); // limits in microseconds

  ROOT::Minuit2::MnMinimize minimize(*this, upar);
  ROOT::Minuit2::FunctionMinimum min = minimize();

  std::cout << "Spectrum " << i << ": ";
  if (min.IsValid()) {
    double dt = min.UserState().Value("dead_time");
    double dt_err = min.UserState().Error("dead_time");
    std::cout << "Dead time = " << dt << " +/- " << dt_err << " us" << std::endl;
  } else {
    std::cout << "Minimization failed" << std::endl;
  }
}

//=============================================================================
// PNeXus Default Constructor
//=============================================================================
nxH4::PNeXus::PNeXus() : fSdId(-1), fFileId(-1)
{
}

//=============================================================================
// PNeXus Constructor with filename
//=============================================================================
nxH4::PNeXus::PNeXus(const std::string fln, const bool printDebug)
  : fPrintDebug(printDebug), fFileName(fln), fSdId(-1), fFileId(-1)
{
  if (ReadNexusFile() != 0) {
    throw std::runtime_error("Failed to read NeXus file: " + fln);
  }
}

//=============================================================================
// PNeXus Destructor
//=============================================================================
nxH4::PNeXus::~PNeXus()
{
  if (fSdId != -1) {
    SDend(fSdId);
  }
  if (fFileId != -1) {
    Hclose(fFileId);
  }
}

//=============================================================================
// PNeXus::ReadNexusFile
//=============================================================================
int nxH4::PNeXus::ReadNexusFile()
{
  // Open the HDF4 file
  fSdId = SDstart(fFileName.c_str(), DFACC_READ);
  if (fSdId == FAIL) {
    std::cerr << "**ERROR** Failed to open HDF4 file: " << fFileName << std::endl;
    return 1;
  }

  fFileId = Hopen(fFileName.c_str(), DFACC_READ, 0);
  if (fFileId == FAIL) {
    std::cerr << "**ERROR** Failed to open HDF4 file (H interface): " << fFileName << std::endl;
    SDend(fSdId);
    fSdId = -1;
    return 1;
  }

  // Read file-level attributes
  int32 n_datasets, n_file_attrs;
  if (SDfileinfo(fSdId, &n_datasets, &n_file_attrs) == FAIL) {
    std::cerr << "**ERROR** Failed to get file info" << std::endl;
    return 1;
  }

  if (fPrintDebug) {
    std::cout << "Number of datasets: " << n_datasets << std::endl;
    std::cout << "Number of file attributes: " << n_file_attrs << std::endl;
  }

  // Read HDF4 version and NeXus version from attributes
  char attr_name[H4_MAX_NC_NAME];
  int32 attr_type, attr_count;

  for (int32 i = 0; i < n_file_attrs; i++) {
    if (SDattrinfo(fSdId, i, attr_name, &attr_type, &attr_count) == FAIL) {
      continue;
    }

    if (caseInsensitiveEquals(attr_name, "HDF_version") ||
        caseInsensitiveEquals(attr_name, "HDF4_version")) {
      std::vector<char> buffer(attr_count + 1, '\0');
      if (SDreadattr(fSdId, i, buffer.data()) != FAIL) {
        fHdf4Version = std::string(buffer.data());
      }
    } else if (caseInsensitiveEquals(attr_name, "NeXus_version")) {
      std::vector<char> buffer(attr_count + 1, '\0');
      if (SDreadattr(fSdId, i, buffer.data()) != FAIL) {
        fNeXusVersion = std::string(buffer.data());
      }
    } else if (caseInsensitiveEquals(attr_name, "file_name")) {
      std::vector<char> buffer(attr_count + 1, '\0');
      if (SDreadattr(fSdId, i, buffer.data()) != FAIL) {
        fFileNameNxs = std::string(buffer.data());
      }
    } else if (caseInsensitiveEquals(attr_name, "file_time")) {
      std::vector<char> buffer(attr_count + 1, '\0');
      if (SDreadattr(fSdId, i, buffer.data()) != FAIL) {
        fFileTimeNxs = std::string(buffer.data());
      }
    }
  }

  // Determine IDF version
  // Try to find IDF_version dataset - check multiple possible locations

  // First try simple root-level name
  int32 idf_idx = SDnametoindex(fSdId, "IDF_version");
  if (idf_idx != FAIL) {
    int32 sds_id = SDselect(fSdId, idf_idx);
    if (sds_id != FAIL) {
      int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs;
      char name[H4_MAX_NC_NAME];
      if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) != FAIL) {
        if (data_type == DFNT_INT32) {
          int32 idf_val;
          int32 start[H4_MAX_VAR_DIMS] = {0};
          int32 edges[H4_MAX_VAR_DIMS] = {1};
          if (SDreaddata(sds_id, start, nullptr, edges, &idf_val) != FAIL) {
            fIdfVersion = idf_val;
          }
        }
      }
      SDendaccess(sds_id);
    }
  }

  // If not found, try /run/IDF_version (IDF version 1 location)
  if (fIdfVersion == -1) {
    try {
      int32 sds_ref = findDatasetRefByPath("/run/IDF_version");
      if (sds_ref != -1) {
        int32 sds_idx = SDreftoindex(fSdId, sds_ref);
        if (sds_idx != FAIL) {
          int32 sds_id = SDselect(fSdId, sds_idx);
          if (sds_id != FAIL) {
            int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs;
            char name[H4_MAX_NC_NAME];
            if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) != FAIL) {
              if (data_type == DFNT_INT32) {
                int32 idf_val;
                int32 start[H4_MAX_VAR_DIMS] = {0};
                int32 edges[H4_MAX_VAR_DIMS] = {1};
                if (SDreaddata(sds_id, start, nullptr, edges, &idf_val) != FAIL) {
                  fIdfVersion = idf_val;
                }
              }
            }
            SDendaccess(sds_id);
          }
        }
      }
    } catch (...) {
      // Path not found, continue to next attempt
    }
  }

  // If still not found, try /raw_data_1/idf_version (IDF version 2 location)
  if (fIdfVersion == -1) {
    try {
      int32 sds_ref = findDatasetRefByPath("/raw_data_1/idf_version");
      if (sds_ref != -1) {
        int32 sds_idx = SDreftoindex(fSdId, sds_ref);
        if (sds_idx != FAIL) {
          int32 sds_id = SDselect(fSdId, sds_idx);
          if (sds_id != FAIL) {
            int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs;
            char name[H4_MAX_NC_NAME];
            if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) != FAIL) {
              if (data_type == DFNT_INT32) {
                int32 idf_val;
                int32 start[H4_MAX_VAR_DIMS] = {0};
                int32 edges[H4_MAX_VAR_DIMS] = {1};
                if (SDreaddata(sds_id, start, nullptr, edges, &idf_val) != FAIL) {
                  fIdfVersion = idf_val;
                }
              }
            }
            SDendaccess(sds_id);
          }
        }
      }
    } catch (...) {
      // Path not found, will use default
    }
  }

  if (fIdfVersion == -1) {
    // Try alternate approach - check for raw_data_1 structure
    // If we find raw_data_1/detector_1, assume IDF v2
    // Otherwise assume IDF v1
    fIdfVersion = 2; // default to version 2
  }

  if (fPrintDebug) {
    std::cout << "HDF4 Version: " << fHdf4Version << std::endl;
    std::cout << "NeXus Version: " << fNeXusVersion << std::endl;
    std::cout << "IDF Version: " << fIdfVersion << std::endl;
  }

  // Read datasets based on IDF version
  if (fIdfVersion == 1) {
    HandleIdfV1(fSdId);
  } else if (fIdfVersion == 2) {
    HandleIdfV2(fSdId);
  } else {
    std::cerr << "**ERROR** Unsupported IDF version: " << fIdfVersion << std::endl;
    return 1;
  }

  return 0;
}

//=============================================================================
// nxH4::PNeXus::HandleIdfV1
//=============================================================================
void nxH4::PNeXus::HandleIdfV1(int32 sd_id)
{
  // IDF version 1 handling - not yet implemented
  if (fPrintDebug) {
    std::cout << "IDF version 1 handling ..." << std::endl;
  }

  std::vector<std::string> int_datasets = {
    "/run/IDF_version",
    "/run/number",
    "/run/switching_state",
    "/run/instrument/detector/number",
    "/run/histogram_data_1/counts",
    "/run/histogram_data_1/resolution",
    "/run/histogram_data_1/grouping"
  };

  std::vector<std::string> float_datasets = {
    "/run/sample/temperature",
    "/run/sample/magnetic_field",
    "/run/sample/magnetic_field_vector",
    "/run/instrument/detector/deadtimes",
    "/run/histogram_data_1/time_zero",
    "/run/histogram_data_1/raw_time",
    "/run/histogram_data_1/corrected_time",
    "/run/histogram_data_1/alpha"
  };

  std::vector<std::string> string_datasets = {
    "/run/program_name",
    "/run/title",
    "/run/notes",
    "/run/analysis",
    "/run/lab",
    "/run/beamline",
    "/run/start_time",
    "/run/stop_time",
    "/run/user/name",
    "/run/user/experiment_number",
    "/run/sample/name",
    "/run/sample/magnetic_field_state",
    "/run/sample/environment",
    "/run/instrument/name",
    "/run/instrument/collimator/type"
  };

  // Read integer datasets
  for (const auto& path : int_datasets) {
    try {
      ReadIntDataset(sd_id, path);
    } catch (const std::exception& e) {
      if (fPrintDebug) {
        std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl;
      }
    }
  }

  // Read float datasets
  for (const auto& path : float_datasets) {
    try {
      ReadFloatDataset(sd_id, path);
    } catch (const std::exception& e) {
      if (fPrintDebug) {
        std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl;
      }
    }
  }

  // Read string datasets
  for (const auto& path : string_datasets) {
    try {
      ReadStringDataset(sd_id, path);
    } catch (const std::exception& e) {
      if (fPrintDebug) {
        std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl;
      }
    }
  }
}

//=============================================================================
// nxH4::PNeXus::HandleIdfV2
//=============================================================================
void nxH4::PNeXus::HandleIdfV2(int32 sd_id)
{
  // Read common datasets for IDF version 2
  // These are typical datasets found in muon NeXus files

  std::vector<std::string> int_datasets = {
    "/raw_data_1/detector_1/counts",
    "/raw_data_1/detector_1/spectrum_index",
    "/raw_data_1/good_frames",
    "/raw_data_1/detector_1/raw_time",
    "/raw_data_1/run_number"
  };

  std::vector<std::string> float_datasets = {
    "/raw_data_1/detector_1/time_resolution",
    "/raw_data_1/detector_1/time_zero"
  };

  std::vector<std::string> string_datasets = {
    "/raw_data_1/name",
    "/raw_data_1/title",
    "/raw_data_1/start_time",
    "/raw_data_1/end_time"
  };

  // Read integer datasets
  for (const auto& path : int_datasets) {
    try {
      ReadIntDataset(sd_id, path);
    } catch (const std::exception& e) {
      if (fPrintDebug) {
        std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl;
      }
    }
  }

  // Read float datasets
  for (const auto& path : float_datasets) {
    try {
      ReadFloatDataset(sd_id, path);
    } catch (const std::exception& e) {
      if (fPrintDebug) {
        std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl;
      }
    }
  }

  // Read string datasets
  for (const auto& path : string_datasets) {
    try {
      ReadStringDataset(sd_id, path);
    } catch (const std::exception& e) {
      if (fPrintDebug) {
        std::cout << "Note: Could not read " << path << ": " << e.what() << std::endl;
      }
    }
  }
}

//=============================================================================
// nxH4::PNeXus::ReadIntDataset
//=============================================================================
void nxH4::PNeXus::ReadIntDataset(int32 sd_id, const std::string& path)
{
  // Extract dataset name from path
  std::vector<std::string> components = splitPath(path);
  if (components.empty()) {
    throw std::runtime_error("Invalid path: " + path);
  }
  std::string dataset_name = components.back();

  // Use VGroup navigation to find the correct dataset reference
  int32 sds_ref = -1;
  if (components.size() > 1) {
    // Try to resolve via VGroup hierarchy
    sds_ref = findDatasetRefByPath(path);
  }

  int32 sds_idx = -1;
  if (sds_ref != -1) {
    // Find SDS index from reference
    sds_idx = SDreftoindex(sd_id, sds_ref);
  }

  // Fallback to name-based search if VGroup navigation failed
  if (sds_idx == FAIL || sds_idx == -1) {
    sds_idx = findDatasetIndex(sd_id, dataset_name);
  }

  int32 sds_id = SDselect(sd_id, sds_idx);
  if (sds_id == FAIL) {
    throw std::runtime_error("Failed to select dataset: " + dataset_name);
  }

  // Get dataset info
  int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs;
  char name[H4_MAX_NC_NAME];
  if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to get dataset info: " + dataset_name);
  }

  // Calculate total number of elements
  int32 n_elements = 1;
  std::vector<uint32_t> dimensions;
  for (int32 i = 0; i < rank; i++) {
    n_elements *= dim_sizes[i];
    dimensions.push_back(static_cast<uint32_t>(dim_sizes[i]));
  }

  // Read data
  std::vector<int32> data(n_elements);
  int32 start[H4_MAX_VAR_DIMS] = {0};
  if (SDreaddata(sds_id, start, nullptr, dim_sizes, data.data()) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to read dataset: " + dataset_name);
  }

  // Convert to int vector
  std::vector<int> int_data(data.begin(), data.end());

  // Create PNXdata object
  PNXdata<int> pnx_data(H4DataType::INT32);
  pnx_data.SetData(int_data);
  pnx_data.SetDimensions(dimensions);

  // Read attributes
  ReadDatasetAttributes(sds_id, pnx_data);

  // Store in data map
  fDataMap[path] = pnx_data;

  SDendaccess(sds_id);
}

//=============================================================================
// nxH4::PNeXus::ReadFloatDataset
//=============================================================================
void nxH4::PNeXus::ReadFloatDataset(int32 sd_id, const std::string& path)
{
  // Extract dataset name from path
  std::vector<std::string> components = splitPath(path);
  if (components.empty()) {
    throw std::runtime_error("Invalid path: " + path);
  }
  std::string dataset_name = components.back();

  // Use VGroup navigation to find the correct dataset reference
  int32 sds_ref = -1;
  if (components.size() > 1) {
    // Try to resolve via VGroup hierarchy
    sds_ref = findDatasetRefByPath(path);
  }

  int32 sds_idx = -1;
  if (sds_ref != -1) {
    // Find SDS index from reference
    sds_idx = SDreftoindex(sd_id, sds_ref);
  }

  // Fallback to name-based search if VGroup navigation failed
  if (sds_idx == FAIL || sds_idx == -1) {
    sds_idx = findDatasetIndex(sd_id, dataset_name);
  }

  int32 sds_id = SDselect(sd_id, sds_idx);
  if (sds_id == FAIL) {
    throw std::runtime_error("Failed to select dataset: " + dataset_name);
  }

  // Get dataset info
  int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs;
  char name[H4_MAX_NC_NAME];
  if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to get dataset info: " + dataset_name);
  }

  // Calculate total number of elements
  int32 n_elements = 1;
  std::vector<uint32_t> dimensions;
  for (int32 i = 0; i < rank; i++) {
    n_elements *= dim_sizes[i];
    dimensions.push_back(static_cast<uint32_t>(dim_sizes[i]));
  }

  // Read data
  std::vector<float32> data(n_elements);
  int32 start[H4_MAX_VAR_DIMS] = {0};
  if (SDreaddata(sds_id, start, nullptr, dim_sizes, data.data()) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to read dataset: " + dataset_name);
  }

  // Convert to float vector
  std::vector<float> float_data(data.begin(), data.end());

  // Create PNXdata object
  PNXdata<float> pnx_data(H4DataType::FLOAT32);
  pnx_data.SetData(float_data);
  pnx_data.SetDimensions(dimensions);

  // Read attributes
  ReadDatasetAttributes(sds_id, pnx_data);

  // Store in data map
  fDataMap[path] = pnx_data;

  SDendaccess(sds_id);
}

//=============================================================================
// nxH4::PNeXus::ReadStringDataset
//=============================================================================
void nxH4::PNeXus::ReadStringDataset(int32 sd_id, const std::string& path)
{
  // Extract dataset name from path
  std::vector<std::string> components = splitPath(path);
  if (components.empty()) {
    throw std::runtime_error("Invalid path: " + path);
  }
  std::string dataset_name = components.back();

  // Use VGroup navigation to find the correct dataset reference
  int32 sds_ref = -1;
  if (components.size() > 1) {
    // Try to resolve via VGroup hierarchy
    sds_ref = findDatasetRefByPath(path);
  }

  int32 sds_idx = -1;
  if (sds_ref != -1) {
    // Find SDS index from reference
    sds_idx = SDreftoindex(sd_id, sds_ref);
  }

  // Fallback to name-based search if VGroup navigation failed
  if (sds_idx == FAIL || sds_idx == -1) {
    sds_idx = findDatasetIndex(sd_id, dataset_name);
  }

  int32 sds_id = SDselect(sd_id, sds_idx);
  if (sds_id == FAIL) {
    throw std::runtime_error("Failed to select dataset: " + dataset_name);
  }

  // Get dataset info
  int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs;
  char name[H4_MAX_NC_NAME];
  if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to get dataset info: " + dataset_name);
  }

  // Calculate total size for string
  int32 n_elements = 1;
  std::vector<uint32_t> dimensions;
  for (int32 i = 0; i < rank; i++) {
    n_elements *= dim_sizes[i];
    dimensions.push_back(static_cast<uint32_t>(dim_sizes[i]));
  }

  // Read data as char array
  std::vector<char> data(n_elements + 1, '\0');
  int32 start[H4_MAX_VAR_DIMS] = {0};
  if (SDreaddata(sds_id, start, nullptr, dim_sizes, data.data()) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to read dataset: " + dataset_name);
  }

  // Create string
  std::vector<std::string> string_data;
  string_data.push_back(std::string(data.data()));

  // Create PNXdata object
  PNXdata<std::string> pnx_data(H4DataType::CHAR8);
  pnx_data.SetData(string_data);
  pnx_data.SetDimensions(dimensions);

  // Read attributes
  ReadDatasetAttributes(sds_id, pnx_data);

  // Store in data map
  fDataMap[path] = pnx_data;

  SDendaccess(sds_id);
}

//=============================================================================
// nxH4::PNeXus::ReadDatasetAttributes (template specializations)
//=============================================================================
template <typename T>
void nxH4::PNeXus::ReadDatasetAttributes(int32 sds_id, PNXdata<T>& data)
{
  int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs;
  char name[H4_MAX_NC_NAME];

  if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_attrs) == FAIL) {
    return;
  }

  for (int32 i = 0; i < n_attrs; i++) {
    char attr_name[H4_MAX_NC_NAME];
    int32 attr_type, attr_count;

    if (SDattrinfo(sds_id, i, attr_name, &attr_type, &attr_count) == FAIL) {
      continue;
    }

    // Read attribute based on type - store as primitives like h5nexus
    if (attr_type == DFNT_INT32) {
      if (attr_count == 1) {
        // Single value - store as int directly
        int32 value;
        if (SDreadattr(sds_id, i, &value) != FAIL) {
          data.AddAttribute(attr_name, static_cast<int>(value));
        }
      } else {
        // Multiple values - store as vector
        std::vector<int32> attr_data(attr_count);
        if (SDreadattr(sds_id, i, attr_data.data()) != FAIL) {
          std::vector<int> int_data(attr_data.begin(), attr_data.end());
          data.AddAttribute(attr_name, int_data);
        }
      }
    } else if (attr_type == DFNT_FLOAT32) {
      if (attr_count == 1) {
        // Single value - store as float directly
        float32 value;
        if (SDreadattr(sds_id, i, &value) != FAIL) {
          data.AddAttribute(attr_name, static_cast<float>(value));
        }
      } else {
        // Multiple values - store as vector
        std::vector<float32> attr_data(attr_count);
        if (SDreadattr(sds_id, i, attr_data.data()) != FAIL) {
          std::vector<float> float_data(attr_data.begin(), attr_data.end());
          data.AddAttribute(attr_name, float_data);
        }
      }
    } else if (attr_type == DFNT_CHAR8 || attr_type == DFNT_UCHAR8) {
      // String attributes - always store as string
      std::vector<char> attr_data(attr_count + 1, '\0');
      if (SDreadattr(sds_id, i, attr_data.data()) != FAIL) {
        data.AddAttribute(attr_name, std::string(attr_data.data()));
      }
    }
  }
}

//=============================================================================
// nxH4::PNeXus::Dump
//=============================================================================
void nxH4::PNeXus::Dump()
{
  int first_good_bin{0};

  if (fIdfVersion == 1) {
    std::cout << std::endl;
    std::cout << "========================================" << std::endl;
    std::cout << "NeXus File Dump" << std::endl;
    std::cout << "========================================" << std::endl;
    std::cout << "Filename: " << fFileName << std::endl;
    std::cout << "HDF4 Version: " << fHdf4Version << std::endl;
    std::cout << "NeXus Version: " << fNeXusVersion << std::endl;
    std::cout << "IDF Version: " << fIdfVersion << std::endl;
    std::cout << "----------------------------------------" << std::endl;

    std::cout << std::endl << "++++";
    std::cout << std::endl << "run";
    std::cout << std::endl << "----";
    std::cout << std::endl << " IDF_version: " << fIdfVersion;
    if (fDataMap.find("/run/program_name") != fDataMap.end()) {
      try {
        auto str_data = std::any_cast<PNXdata<std::string>>(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<std::string>(str_data.GetAttribute("version"));
            std::cout << "   version : " << version;
          } catch (const std::bad_any_cast& e) {
            std::cerr << std::endl << "**ERROR**: Failed to cast version attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast program_name data" << std::endl;
      }
    }
    if (fDataMap.find("/run/number") != fDataMap.end()) {
      try {
        auto number = std::any_cast<PNXdata<int>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(fDataMap["/run/start_time"]);
        std::cout << std::endl << " start_time : " << str_data.GetData()[0];
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast start_time data" << std::endl;
      }
    }
    if (fDataMap.find("/run/stop_time") != fDataMap.end()) {
      try {
        auto str_data = std::any_cast<PNXdata<std::string>>(fDataMap["/run/stop_time"]);
        std::cout << std::endl << " stop_time  : " << str_data.GetData()[0];
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast stop_time data" << std::endl;
      }
    }
    if (fDataMap.find("/run/switching_state") != fDataMap.end()) {
      try {
        auto int_data = std::any_cast<PNXdata<std::string>>(fDataMap["/run/switching_state"]);
        std::cout << std::endl << " switching_state : " << int_data.GetData()[0];
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast switching_state data" << std::endl;
      }
    }
    std::cout << std::endl << "----";
    std::cout << std::endl << " user";
    std::cout << std::endl << "----";
    if (fDataMap.find("/run/user/name") != fDataMap.end()) {
      try {
        auto str_data = std::any_cast<PNXdata<std::string>>(fDataMap["/run/user/name"]);
        std::cout << std::endl << "  name : '" << str_data.GetData()[0] << "'";
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast name data" << std::endl;
      }
    }
    if (fDataMap.find("/run/user/experiment_number") != fDataMap.end()) {
      try {
        auto str_data = std::any_cast<PNXdata<std::string>>(fDataMap["/run/user/experiment_number"]);
        std::cout << std::endl << "  experiment_number : " << str_data.GetData()[0];
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast experiment_number data" << std::endl;
      }
    }
    std::cout << std::endl << "----";
    std::cout << std::endl << " sample";
    std::cout << std::endl << "----";
    if (fDataMap.find("/run/sample/name") != fDataMap.end()) {
      try {
        auto str_data = std::any_cast<PNXdata<std::string>>(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<PNXdata<float>>(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<std::string>(float_data.GetAttribute("units"));
            std::cout << "   units : " << units;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast temperature data" << std::endl;
      }
    }
    if (fDataMap.find("/run/sample/magnetic_field") != fDataMap.end()) {
      try {
        auto float_data = std::any_cast<PNXdata<float>>(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<std::string>(float_data.GetAttribute("units"));
            std::cout << "   units : " << units;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast magnetic_field data" << std::endl;
      }
    }
    if (fDataMap.find("/run/sample/magnetic_field_vector") != fDataMap.end()) {
      try {
        auto float_data = std::any_cast<PNXdata<float>>(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<std::string>(float_data.GetAttribute("coordinate_system"));
            std::cout << std::endl << "   coordinate_system : " << coordinate_system << std::endl;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast coordinate_system attribute" << std::endl;
          }
        }
        if (float_data.HasAttribute("units")) {
          try {
            auto units = std::any_cast<std::string>(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<int>(float_data.GetAttribute("available"));
            std::cout << "   available : " << available;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast available attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast magnetic_field_vector data" << std::endl;
      }
    }
    if (fDataMap.find("/run/sample/environment") != fDataMap.end()) {
      try {
        auto str_data = std::any_cast<PNXdata<std::string>>(fDataMap["/run/sample/environment"]);
        std::cout << std::endl << "  environment : " << str_data.GetData()[0];
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast environment data" << std::endl;
      }
    }
    std::cout << std::endl << "----";
    std::cout << std::endl << " instrument";
    std::cout << std::endl << "----";
    if (fDataMap.find("/run/instrument/name") != fDataMap.end()) {
      try {
        auto str_data = std::any_cast<PNXdata<std::string>>(fDataMap["/run/instrument/name"]);
        std::cout << std::endl << "  name : '" << str_data.GetData()[0] << "'";
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast name data" << std::endl;
      }
    }
    std::cout << std::endl << "----";
    std::cout << std::endl << "  detector";
    std::cout << std::endl << "----";
    if (fDataMap.find("/run/instrument/detector/number") != fDataMap.end()) {
      try {
        auto int_data = std::any_cast<PNXdata<int>>(fDataMap["/run/instrument/detector/number"]);
        std::cout << std::endl << "   number : " << int_data.GetData()[0];
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast number data" << std::endl;
      }
    }
    if (fDataMap.find("/run/instrument/detector/deadtimes") != fDataMap.end()) {
      try {
        auto float_data = std::any_cast<PNXdata<float>>(fDataMap["/run/instrument/detector/deadtimes"]);
        auto data = float_data.GetData();
        unsigned int end{15};
        if (data.size() < end)
          end = data.size();
        std::cout << std::endl << "   deadtimes: ";
        for (unsigned int i=0; i<end; i++)
          std::cout << data[i] << ", ";
        std::cout << "...";
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast deadtimes data" << std::endl;
      }
    }
    std::cout << 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<PNXdata<std::string>>(fDataMap["/run/instrument/collimator/type"]);
        std::cout << std::endl << "   type : " << str_data.GetData()[0];
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast type data" << std::endl;
      }
    }

    std::cout << std::endl << "----";
    std::cout << std::endl << " histogram_data_1";
    std::cout << std::endl << "----";
    std::cout << std::endl << "  counts:";
    std::cout << std::endl;
    try {
      auto counts_data = std::any_cast<PNXdata<int>>(fDataMap["/run/histogram_data_1/counts"]);

      // Check for attributes
      if (counts_data.HasAttribute("units")) {
        try {
          auto units = std::any_cast<std::string>(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<int>(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<int>(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<int>(counts_data.GetAttribute("length"));
          std::cout << "   length : " << histoLength << std::endl;
        } catch (const std::bad_any_cast& e) {
          std::cerr << "**ERROR**: Failed to cast length attribute" << std::endl;
        }
      }
      if (counts_data.HasAttribute("t0_bin")) {
        try {
          auto t0_bin = std::any_cast<int>(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<int>(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<int>(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<float>(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<noOfHistos; i++) {
        if (i<9)
          std::cout << std::endl << "   " << i+1 << ": ";
        else
          std::cout << std::endl << "  " << i+1 << ": ";
        for (unsigned int j=0; j<first_good_bin+5; j++)
          std::cout << data[i*histoLength+j] << ", ";
        std::cout << "...";
      }
    } catch (const std::bad_any_cast& e) {
      std::cerr << "**ERROR**: Failed to cast counts data" << std::endl;
    }

    if (fDataMap.find("/run/histogram_data_1/resolution") != fDataMap.end()) {
      try {
        auto int_data = std::any_cast<PNXdata<int>>(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<std::string>(int_data.GetAttribute("units"));
            std::cout << "   units : " << units;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast resolution data" << std::endl;
      }
    }

    if (fDataMap.find("/run/histogram_data_1/time_zero") != fDataMap.end()) {
      try {
        auto float_data = std::any_cast<PNXdata<float>>(fDataMap["/run/histogram_data_1/time_zero"]);
        std::cout << std::endl << "   time_zero : " << float_data.GetData()[0];

        if (float_data.HasAttribute("units")) {
          try {
            auto units = std::any_cast<std::string>(float_data.GetAttribute("units"));
            std::cout << std::endl << "    units : " << units;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl;
          }
        }
        if (float_data.HasAttribute("available")) {
          try {
            auto available = std::any_cast<int>(float_data.GetAttribute("available"));
            std::cout << std::endl << "    available : " << available;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast available attribute";
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast time_zero data" << std::endl;
      }
    }

    if (fDataMap.find("/run/histogram_data_1/raw_time") != fDataMap.end()) {
      try {
        auto float_data = std::any_cast<PNXdata<float>>(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<int>(float_data.GetAttribute("axis"));
            std::cout << std::endl << "    axis : " << axis;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast axis attribute" << std::endl;
          }
        }
        if (float_data.HasAttribute("primary")) {
          try {
            auto primary = std::any_cast<int>(float_data.GetAttribute("primary"));
            std::cout << std::endl << "    primary : " << primary;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast primary attribute" << std::endl;
          }
        }
        if (float_data.HasAttribute("units")) {
          try {
            auto units = std::any_cast<std::string>(float_data.GetAttribute("units"));
            std::cout << std::endl << "    units : " << units;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast raw_time data" << std::endl;
      }
    }
    if (fDataMap.find("/run/histogram_data_1/corrected_time") != fDataMap.end()) {
      try {
        auto float_data = std::any_cast<PNXdata<float>>(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<int>(float_data.GetAttribute("axis"));
            std::cout << std::endl << "   axis : " << axis;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast axis attribute" << std::endl;
          }
        }
        if (float_data.HasAttribute("units")) {
          try {
            auto units = std::any_cast<std::string>(float_data.GetAttribute("units"));
            std::cout << std::endl << "    units : " << units;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast units attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast raw_time data" << std::endl;
      }
    }
    if (fDataMap.find("/run/histogram_data_1/grouping") != fDataMap.end()) {
      try {
        auto int_data = std::any_cast<PNXdata<int>>(fDataMap["/run/histogram_data_1/grouping"]);
        auto data = int_data.GetData();
        unsigned int end{15};
        if (data.size() < end)
          end = data.size();
        std::cout << std::endl << "   grouping : ";
        for (unsigned int i=0; i<end; i++)
          std::cout << data[i] << ", ";
        std::cout << "...";
        if (int_data.HasAttribute("available")) {
          try {
            auto available = std::any_cast<int>(int_data.GetAttribute("available"));
            std::cout << std::endl << "    available : " << available;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast available attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast grouping data" << std::endl;
      }
    }
    if (fDataMap.find("/run/histogram_data_1/alpha") != fDataMap.end()) {
      try {
        auto float_data = std::any_cast<PNXdata<float>>(fDataMap["/run/histogram_data_1/alpha"]);
        auto data = float_data.GetData();
        unsigned int end{15};
        if (data.size() < end)
          end = data.size();
        std::cout << std::endl << "   alpha : ";
        for (unsigned int i=0; i<end; i++)
          std::cout << data[i] << ", ";
        std::cout << "...";
        if (float_data.HasAttribute("available")) {
          try {
            auto available = std::any_cast<int>(float_data.GetAttribute("available"));
            std::cout << std::endl << "    available : " << available;
          } catch (const std::bad_any_cast& e) {
            std::cerr << "**ERROR**: Failed to cast available attribute" << std::endl;
          }
        }
      } catch (const std::bad_any_cast& e) {
        std::cerr << std::endl << "**ERROR**: Failed to cast alpha data" << std::endl;
      }
    }
    std::cout << std::endl;
    std::cout << "========================================" << std::endl;
    std::cout << std::endl;
  } else { // IDF Version 2
    std::cout << std::endl;
    std::cout << std::endl << "hdf5-NeXus file content of file:' " << fFileName << "'";
    std::cout << std::endl << "****";
    std::cout << std::endl << "Top Level Attributes:";
    std::cout << std::endl << " HDF4 Version : " << fHdf4Version;
    std::cout << std::endl << " NeXus Version: " << fNeXusVersion;
    std::cout << std::endl << " file_name    : " << fFileNameNxs;
    std::cout << std::endl << " file_time    : " << fFileTimeNxs;
    std::cout << std::endl << "++++";
    std::cout << std::endl << "raw_data_1";
    std::cout << std::endl << "----";
    std::cout << std::endl << " IDF_version: " << fIdfVersion;
    if (fDataMap.find("/raw_data_1/beamline") != fDataMap.end()) {
      try {
        auto str = std::any_cast<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<int>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<int>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<std::string>>(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<PNXdata<int>>(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<int>(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<std::string>(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<std::string>(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<int32_t>(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<int32_t>(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<int32_t>(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<PNXdata<int32_t>>(fDataMap["/raw_data_1/instrument/detector_1/resolution"]);
          std::cout << "   resolution     : " << ivalData.GetData()[0];
          if (ivalData.HasAttribute("units")) {
            try {
              auto units = std::any_cast<std::string>(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<dims[1]; i++) {
        if (i<9)
          std::cout << std::endl << "   " << i+1 << ": ";
        else
          std::cout << std::endl << "  " << i+1 << ": ";
        for (unsigned int j=0; j<first_good_bin+5; j++)
          std::cout << data[i*dims[2]+j] << ", ";
        std::cout << "...";
      }
    } catch (const std::bad_any_cast& e) {
      std::cerr << "Error: Failed to cast counts data" << std::endl;
    }

    std::cout << std::endl;
    std::cout << std::endl << "  raw_time:";
    std::cout << std::endl;
    try {
      auto raw_time = std::any_cast<PNXdata<float>>(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<std::string>(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<first_good_bin+5; i++)
        std::cout << data[i] << ", ";
      std::cout << "...";
    } catch (const std::bad_any_cast& e) {
      std::cerr << "Error: Failed to cast raw_time data" << std::endl;
    }

    std::cout << std::endl;
    std::cout << std::endl << "  spectrum_index:";
    std::cout << std::endl;
    try {
      auto spectrum_index = std::any_cast<PNXdata<int>>(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<first_good_bin+5; i++)
        std::cout << data[i] << ", ";
      std::cout << "...";
    } catch (const std::bad_any_cast& e) {
      std::cerr << "Error: Failed to cast spectrum_index data" << std::endl;
    }

    std::cout << std::endl;
    std::cout << std::endl << "  dead_time:";
    std::cout << std::endl;
    try {
      auto dead_time = std::any_cast<PNXdata<float>>(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<first_good_bin+5; i++)
        std::cout << data[i] << ", ";
      std::cout << "...";
    } catch (const std::bad_any_cast& e) {
      std::cerr << "Error: Failed to cast spectrum_index data" << std::endl;
    }

    std::cout << std::endl;
  }
}

//=============================================================================
// nxH4::PNeXus::WriteNexusFile
//=============================================================================
int nxH4::PNeXus::WriteNexusFile(const std::string& filename, int idfVersion)
{
  // Create new HDF4 file
  int32 sd_id = SDstart(filename.c_str(), DFACC_CREATE);
  if (sd_id == FAIL) {
    std::cerr << "**ERROR** Failed to create HDF4 file: " << filename << std::endl;
    return 1;
  }

  // Write file attributes
  WriteFileAttributes(sd_id);

  // Write datasets based on IDF version
  if (idfVersion == 2) {
    WriteIdfV2(sd_id);
  } else {
    std::cerr << "**ERROR** Unsupported IDF version for writing: " << idfVersion << std::endl;
    SDend(sd_id);
    return 1;
  }

  SDend(sd_id);
  return 0;
}

//=============================================================================
// nxH4::PNeXus::WriteFileAttributes
//=============================================================================
void nxH4::PNeXus::WriteFileAttributes(int32 sd_id)
{
  // Write HDF4 version
  if (!fHdf4Version.empty()) {
    SDsetattr(sd_id, "HDF4_version", DFNT_CHAR8,
              fHdf4Version.length(), fHdf4Version.c_str());
  }

  // Write NeXus version
  if (!fNeXusVersion.empty()) {
    SDsetattr(sd_id, "NeXus_version", DFNT_CHAR8,
              fNeXusVersion.length(), fNeXusVersion.c_str());
  }

  // Write file name
  if (!fFileNameNxs.empty()) {
    SDsetattr(sd_id, "file_name", DFNT_CHAR8,
              fFileNameNxs.length(), fFileNameNxs.c_str());
  }

  // Write file time
  if (!fFileTimeNxs.empty()) {
    SDsetattr(sd_id, "file_time", DFNT_CHAR8,
              fFileTimeNxs.length(), fFileTimeNxs.c_str());
  }

  // Write root level attributes from fGroupAttributes
  if (fGroupAttributes.find("/") != fGroupAttributes.end()) {
    const auto& attrs = fGroupAttributes.at("/");
    for (const auto& [attr_name, attr_value] : attrs) {
      // Handle different attribute types
      try {
        auto str_attr = std::any_cast<std::string>(attr_value);
        SDsetattr(sd_id, attr_name.c_str(), DFNT_CHAR8,
                  str_attr.length(), str_attr.c_str());
      } catch (...) {
        // Try other types if needed
      }
    }
  }
}

//=============================================================================
// nxH4::PNeXus::WriteIdfV2
//=============================================================================
void nxH4::PNeXus::WriteIdfV2(int32 sd_id)
{
  // Write all datasets from data map
  for (const auto& [path, data_any] : fDataMap) {
    // Try to write as int dataset
    try {
      auto data = std::any_cast<PNXdata<int>>(data_any);
      WriteIntDataset(sd_id, path, data);
      continue;
    } catch (...) {}

    // Try to write as float dataset
    try {
      auto data = std::any_cast<PNXdata<float>>(data_any);
      WriteFloatDataset(sd_id, path, data);
      continue;
    } catch (...) {}

    // Try to write as string dataset
    try {
      auto data = std::any_cast<PNXdata<std::string>>(data_any);
      WriteStringDataset(sd_id, path, data);
      continue;
    } catch (...) {}
  }
}

//=============================================================================
// nxH4::PNeXus::WriteIntDataset
//=============================================================================
void nxH4::PNeXus::WriteIntDataset(int32 sd_id, const std::string& path,
                             const PNXdata<int>& data)
{
  // Extract dataset name from path
  std::vector<std::string> components = splitPath(path);
  if (components.empty()) {
    throw std::runtime_error("Invalid path: " + path);
  }
  std::string dataset_name = components.back();

  const auto& dims = data.GetDimensions();
  int32 rank = static_cast<int32>(dims.size());
  std::vector<int32> dim_sizes(rank);
  for (int32 i = 0; i < rank; i++) {
    dim_sizes[i] = static_cast<int32>(dims[i]);
  }

  // Create dataset
  int32 sds_id = SDcreate(sd_id, dataset_name.c_str(), DFNT_INT32, rank, dim_sizes.data());
  if (sds_id == FAIL) {
    throw std::runtime_error("Failed to create dataset: " + dataset_name);
  }

  // Convert data to int32
  const auto& int_data = data.GetData();
  std::vector<int32> data32(int_data.begin(), int_data.end());

  // Write data
  int32 start[H4_MAX_VAR_DIMS] = {0};
  if (SDwritedata(sds_id, start, nullptr, dim_sizes.data(), data32.data()) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to write dataset: " + dataset_name);
  }

  // Write attributes
  WriteDatasetAttributes(sds_id, data);

  SDendaccess(sds_id);
}

//=============================================================================
// nxH4::PNeXus::WriteFloatDataset
//=============================================================================
void nxH4::PNeXus::WriteFloatDataset(int32 sd_id, const std::string& path,
                               const PNXdata<float>& data)
{
  // Extract dataset name from path
  std::vector<std::string> components = splitPath(path);
  if (components.empty()) {
    throw std::runtime_error("Invalid path: " + path);
  }
  std::string dataset_name = components.back();

  const auto& dims = data.GetDimensions();
  int32 rank = static_cast<int32>(dims.size());
  std::vector<int32> dim_sizes(rank);
  for (int32 i = 0; i < rank; i++) {
    dim_sizes[i] = static_cast<int32>(dims[i]);
  }

  // Create dataset
  int32 sds_id = SDcreate(sd_id, dataset_name.c_str(), DFNT_FLOAT32, rank, dim_sizes.data());
  if (sds_id == FAIL) {
    throw std::runtime_error("Failed to create dataset: " + dataset_name);
  }

  // Convert data to float32
  const auto& float_data = data.GetData();
  std::vector<float32> data32(float_data.begin(), float_data.end());

  // Write data
  int32 start[H4_MAX_VAR_DIMS] = {0};
  if (SDwritedata(sds_id, start, nullptr, dim_sizes.data(), data32.data()) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to write dataset: " + dataset_name);
  }

  // Write attributes
  WriteDatasetAttributes(sds_id, data);

  SDendaccess(sds_id);
}

//=============================================================================
// nxH4::PNeXus::WriteStringDataset
//=============================================================================
void nxH4::PNeXus::WriteStringDataset(int32 sd_id, const std::string& path,
                                const PNXdata<std::string>& data)
{
  // Extract dataset name from path
  std::vector<std::string> components = splitPath(path);
  if (components.empty()) {
    throw std::runtime_error("Invalid path: " + path);
  }
  std::string dataset_name = components.back();

  const auto& string_data = data.GetData();
  if (string_data.empty()) {
    return;
  }

  std::string str = string_data[0];
  int32 dims[1] = {static_cast<int32>(str.length())};

  // Create dataset
  int32 sds_id = SDcreate(sd_id, dataset_name.c_str(), DFNT_CHAR8, 1, dims);
  if (sds_id == FAIL) {
    throw std::runtime_error("Failed to create dataset: " + dataset_name);
  }

  // Write data - HDF4 requires non-const pointer
  std::vector<char> str_buffer(str.begin(), str.end());
  int32 start[1] = {0};
  if (SDwritedata(sds_id, start, nullptr, dims, str_buffer.data()) == FAIL) {
    SDendaccess(sds_id);
    throw std::runtime_error("Failed to write dataset: " + dataset_name);
  }

  // Write attributes
  WriteDatasetAttributes(sds_id, data);

  SDendaccess(sds_id);
}

//=============================================================================
// nxH4::PNeXus::WriteDatasetAttributes (template specializations)
//=============================================================================
template <typename T>
void nxH4::PNeXus::WriteDatasetAttributes(int32 sds_id, const PNXdata<T>& data)
{
  const auto& attributes = data.GetAttributes();

  for (const auto& [attr_name, attr_value] : attributes) {
    // Try different attribute types - primitives first (like h5nexus)

    // Try scalar int
    try {
      int32 value = std::any_cast<int>(attr_value);
      SDsetattr(sds_id, attr_name.c_str(), DFNT_INT32, 1, &value);
      continue;
    } catch (...) {}

    // Try scalar float
    try {
      float32 value = std::any_cast<float>(attr_value);
      SDsetattr(sds_id, attr_name.c_str(), DFNT_FLOAT32, 1, &value);
      continue;
    } catch (...) {}

    // Try string
    try {
      std::string str = std::any_cast<std::string>(attr_value);
      SDsetattr(sds_id, attr_name.c_str(), DFNT_CHAR8, str.length(), str.c_str());
      continue;
    } catch (...) {}

    // Try vector<int> for multi-element attributes
    try {
      std::vector<int> int_vec = std::any_cast<std::vector<int>>(attr_value);
      std::vector<int32> data32(int_vec.begin(), int_vec.end());
      SDsetattr(sds_id, attr_name.c_str(), DFNT_INT32, data32.size(), data32.data());
      continue;
    } catch (...) {}

    // Try vector<float> for multi-element attributes
    try {
      std::vector<float> float_vec = std::any_cast<std::vector<float>>(attr_value);
      std::vector<float32> data32(float_vec.begin(), float_vec.end());
      SDsetattr(sds_id, attr_name.c_str(), DFNT_FLOAT32, data32.size(), data32.data());
      continue;
    } catch (...) {}
  }
}

//=============================================================================
// Group attribute methods
//=============================================================================
bool nxH4::PNeXus::AddGroupAttribute(const std::string& groupPath,
                           const std::string& attrName,
                           const std::any& attrValue)
{
  fGroupAttributes[groupPath][attrName] = attrValue;
  return true;
}

bool nxH4::PNeXus::RemoveGroupAttribute(const std::string& groupPath,
                                  const std::string& attrName)
{
  auto it = fGroupAttributes.find(groupPath);
  if (it != fGroupAttributes.end()) {
    auto attr_it = it->second.find(attrName);
    if (attr_it != it->second.end()) {
      it->second.erase(attr_it);
      return true;
    }
  }
  return false;
}

bool nxH4::PNeXus::HasGroupAttribute(const std::string& groupPath,
                               const std::string& attrName) const
{
  auto it = fGroupAttributes.find(groupPath);
  if (it != fGroupAttributes.end()) {
    return it->second.find(attrName) != it->second.end();
  }
  return false;
}

std::any nxH4::PNeXus::GetGroupAttribute(const std::string& groupPath,
                           const std::string& attrName) const
{
  return fGroupAttributes.at(groupPath).at(attrName);
}

const std::map<std::string, std::any>& nxH4::PNeXus::GetGroupAttributes(
    const std::string& groupPath) const
{
  static std::map<std::string, std::any> empty_map;
  auto it = fGroupAttributes.find(groupPath);
  if (it != fGroupAttributes.end()) {
    return it->second;
  }
  return empty_map;
}

bool nxH4::PNeXus::ClearGroupAttributes(const std::string& groupPath)
{
  auto it = fGroupAttributes.find(groupPath);
  if (it != fGroupAttributes.end()) {
    it->second.clear();
    return true;
    }
  return false;
}

bool nxH4::PNeXus::AddRootAttribute(const std::string& attrName, const std::any& attrValue)
{
  return AddGroupAttribute("/", attrName, attrValue);
}

//=============================================================================
// Helper methods
//=============================================================================
bool nxH4::PNeXus::caseInsensitiveEquals(const std::string& a, const std::string& b)
{
  if (a.length() != b.length()) {
    return false;
  }
  return std::equal(a.begin(), a.end(), b.begin(),
                    [](char a, char b) {
                      return std::tolower(static_cast<unsigned char>(a)) ==
                             std::tolower(static_cast<unsigned char>(b));
                    });
}

std::vector<std::string> nxH4::PNeXus::splitPath(const std::string& path)
{
  std::vector<std::string> components;
  std::string current;

  for (char c : path) {
    if (c == '/') {
      if (!current.empty()) {
        components.push_back(current);
        current.clear();
      }
    } else {
      current += c;
    }
  }

  if (!current.empty()) {
    components.push_back(current);
  }

  return components;
}

std::string nxH4::PNeXus::findAttributeName(int32 sd_id, const std::string& requestedName)
{
  int32 n_datasets, n_attrs;
  if (SDfileinfo(sd_id, &n_datasets, &n_attrs) == FAIL) {
    throw std::runtime_error("Failed to get file info");
  }

  char attr_name[H4_MAX_NC_NAME];
  int32 attr_type, attr_count;

  for (int32 i = 0; i < n_attrs; i++) {
    if (SDattrinfo(sd_id, i, attr_name, &attr_type, &attr_count) != FAIL) {
      if (caseInsensitiveEquals(attr_name, requestedName)) {
        return std::string(attr_name);
      }
    }
  }

  throw std::runtime_error("Attribute not found: " + requestedName);
}

int32 nxH4::PNeXus::findDatasetIndex(int32 sd_id, const std::string& requestedName)
{
  // First try exact match
  int32 idx = SDnametoindex(sd_id, requestedName.c_str());
  if (idx != FAIL) {
    return idx;
  }

  // Try case-insensitive search
  int32 n_datasets, n_attrs;
  if (SDfileinfo(sd_id, &n_datasets, &n_attrs) == FAIL) {
    throw std::runtime_error("Failed to get file info");
  }

  for (int32 i = 0; i < n_datasets; i++) {
    int32 sds_id = SDselect(sd_id, i);
    if (sds_id != FAIL) {
      char name[H4_MAX_NC_NAME];
      int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_ds_attrs;
      if (SDgetinfo(sds_id, name, &rank, dim_sizes, &data_type, &n_ds_attrs) != FAIL) {
        if (caseInsensitiveEquals(name, requestedName)) {
          SDendaccess(sds_id);
          return i;
        }
      }
      SDendaccess(sds_id);
    }
  }

  throw std::runtime_error("Dataset not found: " + requestedName);
}

int32 nxH4::PNeXus::findDatasetRefByPath(const std::string& path)
{
  // Navigate VGroup hierarchy to find the dataset reference
  // Path format: /group1/group2/.../dataset_name

  std::vector<std::string> components = splitPath(path);
  if (components.empty()) {
    return -1;
  }

  // Open file with V interface
  if (Vstart(fFileId) == FAIL) {
    return -1;
  }

  int32 result_ref = -1;
  int32 current_vg_ref = -1;

  // Find root VGroup (first component after /)
  // For IDF version 1, root is typically "run"
  if (components.size() >= 1) {
    // Get all lone vgroups
    int32 n_vgroups = Vlone(fFileId, nullptr, 0);
    if (n_vgroups > 0) {
      std::vector<int32> vgroup_refs(n_vgroups);
      Vlone(fFileId, vgroup_refs.data(), n_vgroups);

      // Find the root VGroup
      for (int32 i = 0; i < n_vgroups; i++) {
        int32 vg_id = Vattach(fFileId, vgroup_refs[i], "r");
        if (vg_id != FAIL) {
          char vg_name[VGNAMELENMAX];
          if (Vgetname(vg_id, vg_name) != FAIL) {
            if (caseInsensitiveEquals(vg_name, components[0])) {
              current_vg_ref = vgroup_refs[i];
              Vdetach(vg_id);
              break;
            }
          }
          Vdetach(vg_id);
        }
      }
    }
  }

  // Navigate through intermediate VGroups
  for (size_t comp_idx = 1; comp_idx < components.size() - 1; comp_idx++) {
    if (current_vg_ref == -1) break;

    int32 vg_id = Vattach(fFileId, current_vg_ref, "r");
    if (vg_id == FAIL) {
      current_vg_ref = -1;
      break;
    }

    int32 n_entries = Vntagrefs(vg_id);
    bool found = false;

    for (int32 i = 0; i < n_entries; i++) {
      int32 tag, ref;
      if (Vgettagref(vg_id, i, &tag, &ref) != FAIL) {
        // Check if it's a VGroup
        if (tag == DFTAG_VG) {
          int32 child_vg_id = Vattach(fFileId, ref, "r");
          if (child_vg_id != FAIL) {
            char child_name[VGNAMELENMAX];
            if (Vgetname(child_vg_id, child_name) != FAIL) {
              if (caseInsensitiveEquals(child_name, components[comp_idx])) {
                current_vg_ref = ref;
                found = true;
                Vdetach(child_vg_id);
                break;
              }
            }
            Vdetach(child_vg_id);
          }
        }
      }
    }

    Vdetach(vg_id);

    if (!found) {
      current_vg_ref = -1;
      break;
    }
  }

  // Find the dataset in the final VGroup
  if (current_vg_ref != -1 && components.size() >= 1) {
    int32 vg_id = Vattach(fFileId, current_vg_ref, "r");
    if (vg_id != FAIL) {
      int32 n_entries = Vntagrefs(vg_id);
      std::string target_name = components.back();

      for (int32 i = 0; i < n_entries; i++) {
        int32 tag, ref;
        if (Vgettagref(vg_id, i, &tag, &ref) != FAIL) {
          // Check if it's a Numeric Data (SDS)
          if (tag == DFTAG_NDG) {
            // Try to get the name by converting ref to index
            int32 sds_idx = SDreftoindex(fSdId, ref);
            if (sds_idx != FAIL) {
              int32 sds_id = SDselect(fSdId, sds_idx);
              if (sds_id != FAIL) {
                char ds_name[H4_MAX_NC_NAME];
                int32 rank, dim_sizes[H4_MAX_VAR_DIMS], data_type, n_attrs;
                if (SDgetinfo(sds_id, ds_name, &rank, dim_sizes, &data_type, &n_attrs) != FAIL) {
                  if (caseInsensitiveEquals(ds_name, target_name)) {
                    result_ref = ref;
                    SDendaccess(sds_id);
                    break;
                  }
                }
                SDendaccess(sds_id);
              }
            }
          }
        }
      }

      Vdetach(vg_id);
    }
  }

  Vend(fFileId);
  return result_ref;
}

nxH4::H4DataType nxH4::PNeXus::convertHdf4Type(int32 hdf4_type)
{
  switch (hdf4_type) {
  case DFNT_INT32:   return H4DataType::INT32;
  case DFNT_FLOAT32: return H4DataType::FLOAT32;
  case DFNT_FLOAT64: return H4DataType::FLOAT64;
  case DFNT_CHAR8:   return H4DataType::CHAR8;
  case DFNT_UINT32:  return H4DataType::UINT32;
  case DFNT_INT16:   return H4DataType::INT16;
  case DFNT_UINT16:  return H4DataType::UINT16;
  case DFNT_INT8:    return H4DataType::INT8;
  case DFNT_UINT8:   return H4DataType::UINT8;
  default:           return H4DataType::INT32;
  }
}

int32 nxH4::PNeXus::convertToHdf4Type(H4DataType dataType)
{
  switch (dataType) {
  case H4DataType::INT32:   return DFNT_INT32;
  case H4DataType::FLOAT32: return DFNT_FLOAT32;
  case H4DataType::FLOAT64: return DFNT_FLOAT64;
  case H4DataType::CHAR8:   return DFNT_CHAR8;
  case H4DataType::UINT32:  return DFNT_UINT32;
  case H4DataType::INT16:   return DFNT_INT16;
  case H4DataType::UINT16:  return DFNT_UINT16;
  case H4DataType::INT8:    return DFNT_INT8;
  case H4DataType::UINT8:   return DFNT_UINT8;
  default:                  return DFNT_INT32;
  }
}
#endif // HAVE_HDF4