448 lines
17 KiB
C++
448 lines
17 KiB
C++
/***************************************************************************
|
|
|
|
main.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. *
|
|
***************************************************************************/
|
|
|
|
/**
|
|
* @file main.cpp
|
|
* @brief Command-line interface for the h4nexus NeXus HDF4 file reader/writer
|
|
*
|
|
* This file contains the main() function and command-line interface for the
|
|
* h4nexus program. It provides functionality to:
|
|
* - Read and display NeXus HDF4 files
|
|
* - Write NeXus HDF4 files using the PNXdata approach
|
|
* - Calculate dead time corrections for muon detector data
|
|
*
|
|
* **Command-Line Options:**
|
|
* - --fn <file>: Input NeXus HDF4 file (required)
|
|
* - --out <file>: Output NeXus HDF4 file (optional)
|
|
* - --debug, -d: Enable debug output
|
|
* - --dead_time_estimate, -dt: Calculate dead time corrections
|
|
* - --help, -h: Display help message
|
|
* - --version, -v: Display version information
|
|
*
|
|
* @author Andreas Suter
|
|
* @date 2007-2026
|
|
* @copyright GNU General Public License v2
|
|
* @version 1.0
|
|
*
|
|
* @see nxH4::PNeXus
|
|
* @see nxH4::PNeXusDeadTime
|
|
*/
|
|
|
|
#include <cstring>
|
|
#include <iostream>
|
|
#include <string>
|
|
#include <memory>
|
|
#include <ctime>
|
|
|
|
#include <mfhdf.h>
|
|
|
|
#include "PNeXus.h"
|
|
|
|
#ifdef HAVE_CONFIG_H
|
|
#include "config.h"
|
|
#endif
|
|
|
|
#ifdef HAVE_GIT_REV_H
|
|
#include "git-revision.h"
|
|
#endif
|
|
|
|
//-----------------------------------------------------------------------------
|
|
/**
|
|
* @brief Display command-line syntax and help information
|
|
*
|
|
* Prints the usage syntax and available command-line options for the h4nexus
|
|
* program to stdout and exits the program.
|
|
*/
|
|
void h4nexus_syntax() {
|
|
std::cout << std::endl;
|
|
std::cout << "usage: h4nexus [--help | -h] |" << std::endl;
|
|
std::cout << " [--version | -v] |" << std::endl;
|
|
std::cout << " --fn <fln> [--debug | -d]" << std::endl;
|
|
std::cout << " [--dead_time_estimate | -dt]]" << std::endl;
|
|
std::cout << " [--out <fout>]" << std::endl;
|
|
std::cout << std::endl;
|
|
std::cout << "options:" << std::endl;
|
|
std::cout << " --help, -h: this help." << std::endl;
|
|
std::cout << " --version, -v: version of h4nexus." << std::endl;
|
|
std::cout << " --fn <fln>: nexus hdf4 input file name <fn>." << std::endl;
|
|
std::cout << " --dead_time_estimate, -dt: dead time estimate for the read hdf4 nexus file." << std::endl;
|
|
std::cout << " --debug, -d: print additional debug information." << std::endl;
|
|
std::cout << " --out <fout>: write the required datasets of a nexus hdf4 file to file " << std::endl;
|
|
std::cout << " with name <fout>. Only makes sense together with the option --fn <fln>." << std::endl;
|
|
std::cout << std::endl;
|
|
exit(0);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
/**
|
|
* @brief Calculate dead time corrections for muon detector data
|
|
*
|
|
* Estimates dead time corrections for each detector in the NeXus file using
|
|
* the PNeXusDeadTime class and ROOT Minuit2 minimization.
|
|
*
|
|
* @param nxs Pointer to the PNeXus object containing the data
|
|
* @param debug If true, print additional debug information
|
|
*
|
|
* @see nxH4::PNeXusDeadTime
|
|
*/
|
|
void h4nexus_deadTimeEstimate(const nxH4::PNeXus *nxs, bool debug)
|
|
{
|
|
if (debug) {
|
|
std::cout << std::endl;
|
|
std::cout << std::endl << "+++++++++++++++++++";
|
|
std::cout << std::endl << "in deadTimeEstimate";
|
|
std::cout << std::endl << "+++++++++++++++++++";
|
|
std::cout << std::endl;
|
|
}
|
|
|
|
nxH4::PNeXusDeadTime ndt(nxs, debug);
|
|
auto dims = ndt.GetDimensions();
|
|
for (unsigned int i=0; i<dims[1]; i++) {
|
|
ndt.minimize(i);
|
|
}
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
/**
|
|
* @brief Write NeXus HDF4 file using the PNXdata approach
|
|
*
|
|
* Writes all datasets from the PNeXus data map to a new NeXus HDF4 file.
|
|
* This function uses the WriteNexusFile() method to create a complete
|
|
* NeXus file with all groups, datasets, and attributes.
|
|
*
|
|
* @param nxs Pointer to the PNeXus object containing the data to write
|
|
* @param outFileName Output filename for the NeXus HDF4 file
|
|
* @param debug If true, print additional debug information
|
|
*
|
|
* @note Currently only supports IDF version 2 files
|
|
*
|
|
* @see nxH4::PNeXus::WriteNexusFile()
|
|
*/
|
|
void h4nexus_writeTest(const nxH4::PNeXus *nxs, const std::string& outFileName, bool debug)
|
|
{
|
|
if (debug) {
|
|
std::cout << std::endl;
|
|
std::cout << "++++++++++++++++++++" << std::endl;
|
|
std::cout << "Writing NeXus file" << std::endl;
|
|
std::cout << "++++++++++++++++++++" << std::endl;
|
|
std::cout << std::endl;
|
|
}
|
|
|
|
if (nxs->GetIdfVersion() == 1) {
|
|
std::cerr << "Error: IDF v1 write not yet implemented" << std::endl;
|
|
return;
|
|
}
|
|
|
|
// Write using the read object's data
|
|
int result = const_cast<nxH4::PNeXus*>(nxs)->WriteNexusFile(outFileName, nxs->GetIdfVersion());
|
|
|
|
if (result == 0) {
|
|
std::cout << "Successfully wrote: " << outFileName << std::endl;
|
|
} else {
|
|
std::cerr << "Failed to write file: " << outFileName << std::endl;
|
|
}
|
|
|
|
// write data from scratch
|
|
|
|
std::unique_ptr<nxH4::PNeXus> nxs_out = std::make_unique<nxH4::PNeXus>();
|
|
|
|
std::vector<int> ival;
|
|
std::vector<float> fval;
|
|
std::vector<std::string> sval;
|
|
|
|
// ----------
|
|
// raw_data_1
|
|
// ----------
|
|
|
|
// IDF version
|
|
ival.push_back(2);
|
|
nxs_out->AddDataset<int>("/raw_data_1/IDF_version", ival, {1}, nxH4::H4DataType::INT32);
|
|
ival.clear();
|
|
|
|
// add group attribute to '/raw_data_1'
|
|
nxs_out->AddGroupAttribute("/raw_data_1", "NX_class", std::string("NXentry"));
|
|
|
|
// beamline
|
|
sval.push_back("piE3");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/beamline", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// definition
|
|
sval.push_back("muonTD");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/definition", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// run_number
|
|
ival.push_back(1234);
|
|
nxs_out->AddDataset<int>("/raw_data_1/run_number", ival, {1}, nxH4::H4DataType::INT32);
|
|
ival.clear();
|
|
|
|
// title
|
|
sval.push_back("this is the run title.");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/title", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// start time
|
|
sval.push_back("2026-01-01T01:02:03");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/start_time", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// end time
|
|
sval.push_back("2026-01-01T02:03:42");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/end_time", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// experiment_identifier - pgroup for PSI
|
|
sval.push_back("p18324");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/experiment_identifier", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// -------------------
|
|
// detector_1 (NXdata)
|
|
// -------------------
|
|
|
|
// add group attribute to /raw_data_1/instrument
|
|
nxs_out->AddGroupAttribute("/raw_data_1/detector_1", "NX_class", std::string("NXdata"));
|
|
|
|
// counts
|
|
std::vector<int> counts(16*66000, 42); // data 16 histos with length 66000
|
|
nxs_out->AddDataset<int>("/raw_data_1/detector_1/counts", counts, {1, 16, 66000}, nxH4::H4DataType::INT32);
|
|
|
|
// attributes for counts
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "signal", 1);
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "axes", std::string("period_index,spectrum_index,raw_time"));
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "long_name", std::string("positron_counts"));
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "t0_bin", 2741);
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "first_good_bin", 2741);
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "last_good_bin", 66000);
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "units", std::string("counts"));
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "target", std::string("/raw_data_1/instrument/detector_1/counts"));
|
|
|
|
// raw_time
|
|
std::vector<float> raw_time(66000, 0.0);
|
|
for (unsigned int i=0; i<raw_time.size(); i++)
|
|
raw_time[i] = 0.1953125f*1.0e-3*((float)i-2741.0f+0.5f);
|
|
nxs_out->AddDataset<float>("/raw_data_1/detector_1/raw_time", raw_time, {66000}, nxH4::H4DataType::FLOAT32);
|
|
|
|
// attributes raw_time
|
|
nxs_out->AddDatasetAttribute<float>("/raw_data_1/detector_1/raw_time", "units", std::string("microseconds"));
|
|
nxs_out->AddDatasetAttribute<float>("/raw_data_1/detector_1/raw_time", "target", std::string("/raw_data_1/instrument/detector_1/raw_time"));
|
|
|
|
|
|
// ----------
|
|
// instrument
|
|
// ----------
|
|
|
|
// add group attribute to /raw_data_1/instrument
|
|
nxs_out->AddGroupAttribute("/raw_data_1/instrument", "NX_class", std::string("NXinstrument"));
|
|
|
|
// name
|
|
sval.push_back("LEM");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/instrument/name", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// ------
|
|
// source
|
|
// ------
|
|
|
|
// add group attribute to /raw_data_1/instrument/source
|
|
nxs_out->AddGroupAttribute("/raw_data_1/instrument/source", "NX_class", std::string("NXsource"));
|
|
|
|
// name
|
|
sval.push_back("PSI");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/instrument/source/name", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// type
|
|
sval.push_back("continuous muon source");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/instrument/source/types", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// probe
|
|
sval.push_back("postive muons");
|
|
nxs_out->AddDataset<std::string>("/raw_data_1/instrument/source/probe", sval, {1}, nxH4::H4DataType::CHAR8);
|
|
sval.clear();
|
|
|
|
// -----------------------
|
|
// detector_1 (NXdetector)
|
|
// -----------------------
|
|
|
|
// add group attribute to /raw_data_1/instrument/detector_1
|
|
nxs_out->AddGroupAttribute("/raw_data_1/instrument/detector_1", "NX_class", std::string("NXdetector"));
|
|
|
|
// counts
|
|
nxs_out->AddDataset<int>("/raw_data_1/instrument/detector_1/counts", counts, {1, 16, 66000}, nxH4::H4DataType::INT32);
|
|
|
|
// attributes for counts
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "signal", 1);
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "axes", std::string("period_index,spectrum_index,raw_time"));
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "long_name", std::string("positron_counts"));
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "t0_bin", 2741);
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "first_good_bin", 2741);
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "last_good_bin", 66000);
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "units", std::string("counts"));
|
|
nxs_out->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "target", std::string("/raw_data_1/instrument/detector_1/counts"));
|
|
|
|
// raw_time
|
|
nxs_out->AddDataset<float>("/raw_data_1/instrument/detector_1/raw_time", raw_time, {66000}, nxH4::H4DataType::FLOAT32);
|
|
|
|
// attributes raw_time
|
|
nxs_out->AddDatasetAttribute<float>("/raw_data_1/instrument/detector_1/raw_time", "units", std::string("microseconds"));
|
|
nxs_out->AddDatasetAttribute<float>("/raw_data_1/instrument/detector_1/raw_time", "target", std::string("/raw_data_1/instrument/detector_1/raw_time"));
|
|
|
|
// resolution
|
|
fval.push_back(195.3125);
|
|
nxs_out->AddDataset<float>("/raw_data_1/instrument/detector_1/resolution", fval, {1}, nxH4::H4DataType::FLOAT32);
|
|
fval.clear();
|
|
nxs_out->AddDatasetAttribute<float>("/raw_data_1/instrument/detector_1/resolution",
|
|
"units", std::string("picoseconds"));
|
|
|
|
// spectrum_index
|
|
for (unsigned int i=0; i<16; i++)
|
|
ival.push_back(i+1);
|
|
nxs_out->AddDataset<int>("/raw_data_1/instrument/detector_1/spectrum_index", ival, {16}, nxH4::H4DataType::INT32);
|
|
ival.clear();
|
|
|
|
// dead_time
|
|
std::vector<float> deadTime(16, 0.0);
|
|
nxs_out->AddDataset<float>("/raw_data_1/instrument/detector_1/dead_time", deadTime, {16}, nxH4::H4DataType::FLOAT32);
|
|
|
|
// attributes dead_time
|
|
nxs_out->AddDatasetAttribute<float>("/raw_data_1/instrument/detector_1/dead_time", "available", 0);
|
|
nxs_out->AddDatasetAttribute<float>("/raw_data_1/instrument/detector_1/dead_time", "units", std::string("microseconds"));
|
|
nxs_out->AddDatasetAttribute<float>("/raw_data_1/instrument/detector_1/dead_time", "target", std::string("/raw_data_1/instrument/detector_1/dead_time"));
|
|
|
|
// add root attributes
|
|
// file name
|
|
nxs_out->AddRootAttribute("file_name", std::string("_test.nxs"));
|
|
// date-time
|
|
std::time_t time = std::time({});
|
|
char timeString[std::size("yyyy-mm-ddThh:mm:ssZ")];
|
|
std::strftime(std::data(timeString), std::size(timeString),
|
|
"%FT%TZ", std::gmtime(&time));
|
|
nxs_out->AddRootAttribute("file_time", std::string(timeString));
|
|
// NeXus version
|
|
nxs_out->AddRootAttribute("NeXus_Version", std::string("4.3.0"));
|
|
// hdf4 version
|
|
nxs_out->AddRootAttribute("HDF4_Version", std::string(nxs->GetHdf4Version()));
|
|
// creator
|
|
nxs_out->AddRootAttribute("creator", std::string("h4nexus - PSI"));
|
|
|
|
nxs_out->WriteNexusFile("_test.nxs");
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
/**
|
|
* @brief Main entry point for the h4nexus program
|
|
*
|
|
* Parses command-line arguments and performs the requested operations:
|
|
* - Read and display NeXus HDF4 file information
|
|
* - Write NeXus HDF4 files with modified data
|
|
* - Calculate dead time corrections
|
|
*
|
|
* @param argc Number of command-line arguments
|
|
* @param argv Array of command-line argument strings
|
|
*
|
|
* @return 0 on success, non-zero on error
|
|
*
|
|
* @see h4nexus_syntax() for available command-line options
|
|
*/
|
|
int main(int argc, char *argv[])
|
|
{
|
|
std::string fileName{""};
|
|
std::string fileNameOut{""};
|
|
bool printDebug{false};
|
|
bool deadTimeEstimate{false};
|
|
|
|
if (argc == 1)
|
|
h4nexus_syntax();
|
|
|
|
for (int i=1; i<argc; i++) {
|
|
if (!strcmp(argv[i], "-h") || !strcmp(argv[i], "--help")) {
|
|
h4nexus_syntax();
|
|
} else if (!strcmp(argv[i], "-v") || !strcmp(argv[i], "--version")) {
|
|
#ifdef HAVE_CONFIG_H
|
|
#ifdef HAVE_GIT_REV_H
|
|
std::cout << std::endl << "h4nexus version: " << H4NEXUS_VERSION << ", git-branch: " << GIT_BRANCH << ", git-hash: " << GIT_HASH << std::endl << std::endl;
|
|
#else
|
|
std::cout << std::endl << "h4nexus version: " << H4NEXUS_VERSION << std::endl << std::endl;
|
|
#endif
|
|
#else
|
|
#ifdef HAVE_GIT_REV_H
|
|
std::cout << std::endl << "h4nexus git-branch: " << GIT_BRANCH << ", git-hash: " << GIT_HASH << std::endl << std::endl;
|
|
#else
|
|
std::cout << std::endl << "h4nexus version: unknown." << std::endl << std::endl;
|
|
#endif
|
|
#endif
|
|
return 0;
|
|
} else if (!strcmp(argv[i], "--fn")) {
|
|
if (i+1 >= argc) {
|
|
std::cout << std::endl << "**ERROR** found --fn without <fln>." << std::endl;
|
|
h4nexus_syntax();
|
|
}
|
|
i++;
|
|
fileName = argv[i];
|
|
} else if (!strcmp(argv[i], "-dt") || !strcmp(argv[i], "--dead_time_estimate")) {
|
|
deadTimeEstimate = true;
|
|
} else if (!strcmp(argv[i], "-d") || !strcmp(argv[i], "--debug")) {
|
|
printDebug = true;
|
|
} else if (!strcmp(argv[i], "--out")) {
|
|
if (i+1 >= argc) {
|
|
std::cout << std::endl << "**ERROR** found --out without <fout>." << std::endl;
|
|
h4nexus_syntax();
|
|
}
|
|
i++;
|
|
fileNameOut = argv[i];
|
|
} else {
|
|
h4nexus_syntax();
|
|
}
|
|
}
|
|
|
|
if (fileName.empty()) {
|
|
std::cerr << std::endl;
|
|
std::cerr << "**ERROR** <fln> is missing." << std::endl;
|
|
std::cerr << std::endl;
|
|
h4nexus_syntax();
|
|
}
|
|
|
|
std::unique_ptr<nxH4::PNeXus> nxs = std::make_unique<nxH4::PNeXus>(fileName, printDebug);
|
|
|
|
nxs->Dump();
|
|
|
|
if (deadTimeEstimate) {
|
|
h4nexus_deadTimeEstimate(nxs.get(), printDebug);
|
|
}
|
|
|
|
if (!fileNameOut.empty()) {
|
|
h4nexus_writeTest(nxs.get(), fileNameOut, printDebug);
|
|
}
|
|
|
|
return 0;
|
|
}
|