LMU/musrfit
5
1
Fork 0
Code Issues 5 Pull Requests Actions Packages Projects Releases 4 Activity

adopted the PNeXus lib such that the class naming scheme is self-consistent.

Browse Source
This commit is contained in:
suter_a
2026-02-21 07:44:20 +01:00
parent ee3fc2e42a
commit 17fcaa7c69
3 changed files with 8442 additions and 5906 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/classes/PRunDataHandler.cpp
+747 -5
View File
@@ -5067,11 +5067,753 @@ Bool_t PRunDataHandler::WriteNexusFile(TString fln)
if (!fAny2ManyInfo->useStandardOutput)
std::cout << std::endl << ">> PRunDataHandler::WriteNexusFile(): writing a NeXus data file (" << fln.Data() << ") ... " << std::endl;
// create NeXus object
std::unique_ptr<PNeXus> nxs = std::make_unique<PNeXus>();
if (nxs == nullptr) {
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile(): **ERROR** couldn't invoke the NeXus object." << std::endl;
return false;
if (format.Contains("HDF4", TString::kIgnoreCase)) { // HDF4
#ifdef HAVE_HDF4
try {
// create NeXus object
std::unique_ptr<nxH4::PNeXus> nxs = std::make_unique<nxH4::PNeXus>();
if (nxs == nullptr) {
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile(): **ERROR** couldn't invoke the NeXus object." << std::endl;
return false;
}
// set NeXus version
nxs->AddGroupAttribute("/", "NeXus_version", std::string("4.3.0"));
// set HDF4 version
nxs->AddGroupAttribute("/", "HDF_version", nxs->GetHdf4LibVersion());
// set file name
nxs->AddGroupAttribute("/", "file_name", std::string(fln.Data()));
// set creation time
std::string dt = nxs::getIso8601TimestampLocal();
nxs->AddGroupAttribute("/", "file_time", dt);
if (fAny2ManyInfo->idf == 1) { // IDF V1
// set IDF version
nxs->AddDataset<int>("/run/IDF_version", {(int)fAny2ManyInfo->idf}, {1}, nxH4::H4DataType::kINT32);
// set program name
nxs->AddDataset<std::string>("/run/program_name", {"any2many"}, {1}, nxH4::H4DataType::kCHAR8);
str="n/a";
#ifdef HAVE_CONFIG_H
str = PACKAGE_VERSION;
#endif
nxs->AddDatasetAttribute<std::string>("/run/program_name", "version", str);
// set run number
nxs->AddDataset<int>("/run/number", {fData[0].GetRunNumber()}, {1}, nxH4::H4DataType::kINT32);
// set title
nxs->AddDataset<std::string>("/run/title", {fData[0].GetRunTitle()->Data()}, {1}, nxH4::H4DataType::kCHAR8);
// set notes
nxs->AddDataset<std::string>("/run/notes", {std::string("n/a")}, {1}, nxH4::H4DataType::kCHAR8);
// set analysis
nxs->AddDataset<std::string>("/run/analysis", {std::string("muonTD")}, {1}, nxH4::H4DataType::kCHAR8);
// set lab
str = *fData[0].GetLaboratory();
nxs->AddDataset<std::string>("/run/lab", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set beamline
str = *fData[0].GetBeamline();
nxs->AddDataset<std::string>("/run/beamline", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set start time
str = std::string(fData[0].GetStartDate()->Data()) + std::string("T") + std::string(fData[0].GetStartTime()->Data());
nxs->AddDataset<std::string>("/run/start_time", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set stop time
str = std::string(fData[0].GetStopDate()->Data()) + std::string("T") + std::string(fData[0].GetStopTime()->Data());
nxs->AddDataset<std::string>("/run/stop_time", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set switching state
nxs->AddDataset<int>("/run/switching_states", {1}, {1}, nxH4::H4DataType::kINT32);
// set user name
nxs->AddDataset<std::string>("/run/user/name", {std::string("n/a")}, {1}, nxH4::H4DataType::kCHAR8);
// set user experiment_number
nxs->AddDataset<std::string>("/run/user/experiment_number", {std::string("n/a")}, {1}, nxH4::H4DataType::kCHAR8);
// set sample name
nxs->AddDataset<std::string>("/run/sample/name", {fData[0].GetSample()->Data()}, {1}, nxH4::H4DataType::kCHAR8);
// set sample temperature
nxs->AddDataset<float>("/run/sample/temperature", {(float)fData[0].GetTemperature(0)}, {1}, nxH4::H4DataType::kFLOAT32);
nxs->AddDatasetAttribute<float>("/run/sample/temperature", "units", std::string("Kelvin"));
// set magnetic field
nxs->AddDataset<float>("/run/sample/magnetic_field", {(float)fData[0].GetField()}, {1}, nxH4::H4DataType::kFLOAT32);
nxs->AddDatasetAttribute<float>("/run/sample/magnetic_field", "units", std::string("Gauss"));
// set sample environment
nxs->AddDataset<std::string>("/run/sample/environment", {fData[0].GetSetup()->Data()}, {1}, nxH4::H4DataType::kCHAR8);
// set sample shape
nxs->AddDataset<std::string>("/run/sample/shape", {std::string("n/a")}, {1}, nxH4::H4DataType::kCHAR8);
// set magnetic field vector
nxs->AddDataset<float>("/run/sample/magnetic_field_vector", {1.0f, 1.0f, 1.0f}, {3}, nxH4::H4DataType::kFLOAT32);
nxs->AddDatasetAttribute<float>("/run/sample/magnetic_field_vector", "coordinate_system", std::string("cartesian"));
nxs->AddDatasetAttribute<float>("/run/sample/magnetic_field_vector", "units", std::string("Gauss"));
nxs->AddDatasetAttribute<float>("/run/sample/magnetic_field_vector", "available", 0);
// set instrument name
str = *fData[0].GetInstrument();
nxs->AddDataset<std::string>("/run/instrument/name", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set instrument number of detectors
nxs->AddDataset<int>("/run/instrument/detector/number", {(int)fData[0].GetNoOfHistos()}, {1}, nxH4::H4DataType::kINT32);
// set instrument collimator
nxs->AddDataset<std::string>("/run/instrument/collimator/type", {std::string("n/a")}, {1}, nxH4::H4DataType::kCHAR8);
// set instrument beam total number of counts in Mev
// calculate the total number of counts
double total_counts = 0;
PRawRunDataSet *dataSet = nullptr;
for (unsigned int i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=0" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
for (unsigned int j=0; j<dataSet->GetData()->size(); j++)
total_counts += dataSet->GetData()->at(j);
}
float total_counts_mev = (float) total_counts / 1.0e6;
nxs->AddDataset<float>("/run/instrument/beam/total_counts", {total_counts_mev}, {1}, nxH4::H4DataType::kFLOAT32);
nxs->AddDatasetAttribute<float>("/run/instrument/beam/total_counts", "units", std::string("MEv"));
// set time resolution (use FLOAT instead of INT)
float res = (float)(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin*1.0e3);
nxs->AddDataset<float>("/run/histogram_data_1/resolution", {res}, {1}, nxH4::H4DataType::kFLOAT32);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/resolution", "units", std::string("picoseconds"));
// set time zero time to 0. see t0_bin attribute of counts!
nxs->AddDataset<int>("/run/histogram_data_1/time_zero", {0}, {1}, nxH4::H4DataType::kINT32);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/time_zero", "units", std::string("microseconds"));
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/time_zero", "available", 0);
// set raw_time
res = (float)(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin*1.0e-3);
dataSet = fData[0].GetDataSet(0, false); // i.e. the false means, that i is the index and NOT the histo number
unsigned int length = (int)(dataSet->GetData()->size() / fAny2ManyInfo->rebin);
std::vector<float> time;
for (unsigned int i=0; i<length; i++)
time.push_back(((float)i+0.5)*res);
nxs->AddDataset<float>("/run/histogram_data_1/raw_time", time, {length}, nxH4::H4DataType::kFLOAT32);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/raw_time", "axis", 1);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/raw_time", "primary", 0);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/raw_time", "units", std::string("microseconds"));
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/raw_time", "available", 0);
// set corrected_time
nxs->AddDataset<float>("/run/histogram_data_1/corrected_time", time, {length}, nxH4::H4DataType::kFLOAT32);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/corrected_time", "axis", 1);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/corrected_time", "units", std::string("microseconds"));
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/corrected_time", "available", 0);
// set grouping
std::vector<int> grouping(fData[0].GetNoOfHistos(), 0);
nxs->AddDataset<int>("/run/histogram_data_1/grouping", grouping, {(uint32_t)grouping.size()}, nxH4::H4DataType::kINT32);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/grouping", "avaliabe", 0);
// set alpha
int ival=1;
nxs->AddDataset<int>("/run/histogram_data_1/alpha", {ival}, {1}, nxH4::H4DataType::kINT32);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/alpha", "avaliabe", 0);
// set counts
// feed histos
PIntVector data;
UInt_t size = 0;
int noHisto, histoLength;
if (fAny2ManyInfo->rebin == 1) {
noHisto = fData[0].GetNoOfHistos();
for (UInt_t i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
size = dataSet->GetData()->size();
histoLength = size;
for (UInt_t j=0; j<size; j++) {
data.push_back((UInt_t)dataSet->GetData()->at(j));
}
}
} else { // rebin > 1
UInt_t dataRebin = 0;
UInt_t dataCount = 0;
noHisto = fData[0].GetNoOfHistos();
for (UInt_t i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
size = dataSet->GetData()->size();
dataCount = 0;
for (UInt_t j=0; j<size; j++) {
if ((j > 0) && (j % fAny2ManyInfo->rebin == 0)) {
dataCount++;
data.push_back(dataRebin);
dataRebin = 0;
}
dataRebin += static_cast<UInt_t>(dataSet->GetData()->at(j));
}
}
size = dataCount;
}
nxs->AddDataset<int>("/run/histogram_data_1/counts", data, {(uint32_t)noHisto, (uint32_t)size}, nxH4::H4DataType::kINT32);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "units", std::string("counts"));
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "signal", 1);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "number", noHisto);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "length", (int)size);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "t0_bin", 0);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "first_good_bin", 0);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "last_good_bin", 0);
res = (float)(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin*1.0e3)/2.0;
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "offset", res);
} else { // IDF V2
nxs->AddGroupAttribute("/raw_data_1", "NX_class", std::string("NXentry"));
// set IDF version
nxs->AddDataset<int>("/raw_data_1/IDF_version", {(int)fAny2ManyInfo->idf}, {1}, nxH4::H4DataType::kINT32);
// set beamline
str = *fData[0].GetBeamline();
nxs->AddDataset<std::string>("/raw_data_1/beamline", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set definition
nxs->AddDataset<std::string>("/raw_data_1/definition", {std::string("muonTD")}, {1}, nxH4::H4DataType::kCHAR8);
// set run_number
nxs->AddDataset<int>("/raw_data_1/run_number", {fData[0].GetRunNumber()}, {1}, nxH4::H4DataType::kINT32);
// set title
nxs->AddDataset<std::string>("/raw_data_1/title", {fData[0].GetRunTitle()->Data()}, {1}, nxH4::H4DataType::kCHAR8);
// set start time
str = std::string(fData[0].GetStartDate()->Data()) + std::string("T") + std::string(fData[0].GetStartTime()->Data());
nxs->AddDataset<std::string>("/raw_data_1/start_time", {str}, {1}, nxH4::H4DataType::kCHAR8);
nxs->AddDatasetAttribute<std::string>("/raw_data_1/start_time", "units", "ISO8601");
// set end time
str = std::string(fData[0].GetStopDate()->Data()) + std::string("T") + std::string(fData[0].GetStopTime()->Data());
nxs->AddDataset<std::string>("/raw_data_1/end_time", {str}, {1}, nxH4::H4DataType::kCHAR8);
nxs->AddDatasetAttribute<std::string>("/raw_data_1/end_time", "units", "ISO8601");
// set experiment_identifier
str = "n/a";
nxs->AddDataset<std::string>("/raw_data_1/experiment_identifier", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set instrument attribute
nxs->AddGroupAttribute("/raw_data_1/instrument", "NX_class", std::string("NXinstrument"));
// set instrument name
str = *fData[0].GetInstrument();
nxs->AddDataset<std::string>("/raw_data_1/instrument/name", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set instrument/source attribute
nxs->AddGroupAttribute("/raw_data_1/instrument/source", "NX_class", std::string("NXsource"));
// set instrument/source/name
str = fData[0].GetLaboratory()->Data();
nxs->AddDataset<std::string>("/raw_data_1/instrument/source/name", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set instrument/source/type
TString tstr = *fData[0].GetInstrument();
std::string type{"n/a"};
if (tstr.Contains("LEM", TString::kIgnoreCase)) {
type = "low energy muon source";
} else if (tstr.Contains("GPS", TString::kIgnoreCase) || tstr.Contains("GPD", TString::kIgnoreCase) ||
tstr.Contains("LTF", TString::kIgnoreCase) || tstr.Contains("FLAME", TString::kIgnoreCase) ||
tstr.Contains("HAL-9500", TString::kIgnoreCase) || tstr.Contains("DOLLY", TString::kIgnoreCase) ||
tstr.Contains("VMS", TString::kIgnoreCase)) {
type = "quasi-continous muon source";
} else if (tstr.Contains("EMU", TString::kIgnoreCase) || tstr.Contains("MUSR", TString::kIgnoreCase) ||
tstr.Contains("HIFI", TString::kIgnoreCase)) {
type = "pulsed muon source";
}
nxs->AddDataset<std::string>("/raw_data_1/instrument/source/type", {type}, {1}, nxH4::H4DataType::kCHAR8);
// set instrument/source/probe
str = "positive muons";
nxs->AddDataset<std::string>("/raw_data_1/instrument/source/probe", {str}, {1}, nxH4::H4DataType::kCHAR8);
// set instrument/detector info
nxs->AddGroupAttribute("/raw_data_1/instrument/detector_1", "NX_class", std::string("NXdetector"));
// set instrument/detector/spectrum_index
int noHistos = fData[0].GetNoOfHistos();
std::vector<int> spectrum_index(noHistos);
for (unsigned int i=0; i<spectrum_index.size(); i++)
spectrum_index[i] = i+1;
nxs->AddDataset<int>("/raw_data_1/instrument/detector_1/spectrum_index", spectrum_index, {(uint32_t)spectrum_index.size()}, nxH4::H4DataType::kINT32);
// set instrument/detector/raw_time (not useful for quasi-continuous sources)
int ival=0;
nxs->AddDataset<int>("/raw_data_1/instrument/detector_1/raw_time", {ival}, {1}, nxH4::H4DataType::kINT32);
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/raw_time", "available", 0);
// set instrument/detector/counts
// set counts
// feed histos
PRawRunDataSet *dataSet = nullptr;
PIntVector data;
UInt_t size = 0;
int noHisto, histoLength;
if (fAny2ManyInfo->rebin == 1) {
noHisto = fData[0].GetNoOfHistos();
for (UInt_t i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
size = dataSet->GetData()->size();
histoLength = size;
for (UInt_t j=0; j<size; j++) {
data.push_back((UInt_t)dataSet->GetData()->at(j));
}
}
} else { // rebin > 1
UInt_t dataRebin = 0;
UInt_t dataCount = 0;
noHisto = fData[0].GetNoOfHistos();
for (UInt_t i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
size = dataSet->GetData()->size();
dataCount = 0;
for (UInt_t j=0; j<size; j++) {
if ((j > 0) && (j % fAny2ManyInfo->rebin == 0)) {
dataCount++;
data.push_back(dataRebin);
dataRebin = 0;
}
dataRebin += static_cast<UInt_t>(dataSet->GetData()->at(j));
}
}
size = dataCount;
}
nxs->AddDataset<int>("/raw_data_1/instrument/detector_1/counts", data, {(uint32_t)noHisto, (uint32_t)size}, nxH4::H4DataType::kINT32);
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "axes", std::string("spectrum_index,time_bin"));
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "long_name", std::string("positon counts"));
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "signal", 1);
// t0_bin attributes
std::vector<int> t0_bin;
for (unsigned int i=0; i<fData[0].GetNoOfHistos(); i++)
t0_bin.push_back((int)(fData[0].GetT0Bin(i+1)/fAny2ManyInfo->rebin));
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "t0_bin", {t0_bin});
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "units", std::string("counts"));
// set raw_data_1/detector info
nxs->AddGroupAttribute("/raw_data_1/detector_1", "NX_class", std::string("NXdata"));
// set detector/counts
nxs->AddDataset<int>("/raw_data_1/detector_1/counts", data, {(uint32_t)noHisto, (uint32_t)size}, nxH4::H4DataType::kINT32);
nxs->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "axes", std::string("spectrum_index,time_bin"));
nxs->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "long_name", std::string("positon counts"));
nxs->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "signal", 1);
}
int result = nxs->WriteNexusFile(fln.Data(), fAny2ManyInfo->idf);
if (result != 0) {
std::cerr << std::endl << "**ERROR** PRunDataHandler::WriteNexusFile, fln=" << fln << std::endl;
return false;
}
} catch (const std::runtime_error& e) {
std::cerr << std::endl << "HDF4 error: " << e.what() << std::endl;
return false;
}
#endif
} else { // HDF5
try {
// create NeXus object
std::unique_ptr<nxH5::PNeXus> nxs = std::make_unique<nxH5::PNeXus>();
if (nxs == nullptr) {
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile(): **ERROR** couldn't invoke the NeXus object." << std::endl;
return false;
}
// set NeXus version
nxs->AddGroupAttribute("/", "NeXus_version", std::string("4.3.0"));
// set HDF5 version
nxs->AddGroupAttribute("/", "HDF_version", nxs->GetHdf5LibVersion());
// set file name
nxs->AddGroupAttribute("/", "file_name", std::string(fln.Data()));
// set creation time
std::string dt = nxs::getIso8601TimestampLocal();
nxs->AddGroupAttribute("/", "file_time", dt);
if (fAny2ManyInfo->idf == 1) { // IDF V1
// set IDF version
nxs->AddDataset<int>("/run/IDF_version", {(int)fAny2ManyInfo->idf}, {1}, H5::PredType::NATIVE_INT32);
// set program name
nxs->AddDataset<std::string>("/run/program_name", {"any2many"}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
str="n/a";
#ifdef HAVE_CONFIG_H
str = PACKAGE_VERSION;
#endif
nxs->AddDatasetAttribute<std::string>("/run/program_name", "version", str);
// set run number
nxs->AddDataset<int>("/run/number", {fData[0].GetRunNumber()}, {1}, H5::PredType::NATIVE_INT32);
// set title
nxs->AddDataset<std::string>("/run/title", {fData[0].GetRunTitle()->Data()}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set notes
nxs->AddDataset<std::string>("/run/notes", {std::string("n/a")}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set analysis
nxs->AddDataset<std::string>("/run/analysis", {std::string("muonTD")}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set lab
str = *fData[0].GetLaboratory();
nxs->AddDataset<std::string>("/run/lab", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set beamline
str = *fData[0].GetBeamline();
nxs->AddDataset<std::string>("/run/beamline", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set start time
str = std::string(fData[0].GetStartDate()->Data()) + std::string("T") + std::string(fData[0].GetStartTime()->Data());
nxs->AddDataset<std::string>("/run/start_time", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set stop time
str = std::string(fData[0].GetStopDate()->Data()) + std::string("T") + std::string(fData[0].GetStopTime()->Data());
nxs->AddDataset<std::string>("/run/stop_time", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set switching state
nxs->AddDataset<int>("/run/switching_states", {1}, {1}, H5::PredType::NATIVE_INT32);
// set user name
nxs->AddDataset<std::string>("/run/user/name", {std::string("n/a")}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set user experiment_number
nxs->AddDataset<std::string>("/run/user/experiment_number", {std::string("n/a")}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set sample name
nxs->AddDataset<std::string>("/run/sample/name", {fData[0].GetSample()->Data()}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set sample temperature
nxs->AddDataset<float>("/run/sample/temperature", {(float)fData[0].GetTemperature(0)}, {1}, H5::PredType::NATIVE_FLOAT);
nxs->AddDatasetAttribute<float>("/run/sample/temperature", "units", std::string("Kelvin"));
// set magnetic field
nxs->AddDataset<float>("/run/sample/magnetic_field", {(float)fData[0].GetField()}, {1}, H5::PredType::NATIVE_FLOAT);
nxs->AddDatasetAttribute<float>("/run/sample/magnetic_field", "units", std::string("Gauss"));
// set sample environment
nxs->AddDataset<std::string>("/run/sample/environment", {fData[0].GetSetup()->Data()}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set sample shape
nxs->AddDataset<std::string>("/run/sample/shape", {std::string("n/a")}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set magnetic field vector
nxs->AddDataset<float>("/run/sample/magnetic_field_vector", {1.0, 1.0, 1.0}, {3}, H5::PredType::NATIVE_FLOAT);
nxs->AddDatasetAttribute<float>("/run/sample/magnetic_field_vector", "coordinate_system", std::string("cartesian"));
nxs->AddDatasetAttribute<float>("/run/sample/magnetic_field_vector", "units", std::string("Gauss"));
nxs->AddDatasetAttribute<float>("/run/sample/magnetic_field_vector", "available", 0);
// set instrument name
str = *fData[0].GetInstrument();
nxs->AddDataset<std::string>("/run/instrument/name", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set instrument number of detectors
nxs->AddDataset<int>("/run/instrument/detector/number", {(int)fData[0].GetNoOfHistos()}, {1}, H5::PredType::NATIVE_INT32);
// set instrument collimator
nxs->AddDataset<std::string>("/run/instrument/collimator/type", {std::string("n/a")}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set instrument beam total number of counts in Mev
// calculate the total number of counts
double total_counts = 0;
PRawRunDataSet *dataSet = nullptr;
for (unsigned int i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=0" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
for (unsigned int j=0; j<dataSet->GetData()->size(); j++)
total_counts += dataSet->GetData()->at(j);
}
float total_counts_mev = (float) total_counts / 1.0e6;
nxs->AddDataset<float>("/run/instrument/beam/total_counts", {total_counts_mev}, {1}, H5::PredType::NATIVE_FLOAT);
nxs->AddDatasetAttribute<float>("/run/instrument/beam/total_counts", "units", std::string("MEv"));
// set time resolution (use FLOAT instead of INT)
float res = (float)(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin*1.0e3);
nxs->AddDataset<float>("/run/histogram_data_1/resolution", {res}, {1}, H5::PredType::NATIVE_FLOAT);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/resolution", "units", std::string("picoseconds"));
// set time zero time to 0. see t0_bin attribute of counts!
nxs->AddDataset<int>("/run/histogram_data_1/time_zero", {0}, {1}, H5::PredType::NATIVE_INT32);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/time_zero", "units", std::string("microseconds"));
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/time_zero", "available", 0);
// set raw_time
res = (float)(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin*1.0e-3);
dataSet = fData[0].GetDataSet(0, false); // i.e. the false means, that i is the index and NOT the histo number
unsigned int length = (int)(dataSet->GetData()->size() / fAny2ManyInfo->rebin);
std::vector<float> time;
for (unsigned int i=0; i<length; i++)
time.push_back(((float)i+0.5)*res);
nxs->AddDataset<float>("/run/histogram_data_1/raw_time", time, {length}, H5::PredType::NATIVE_FLOAT);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/raw_time", "axis", 1);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/raw_time", "primary", 0);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/raw_time", "units", std::string("microseconds"));
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/raw_time", "available", 0);
// set corrected_time
nxs->AddDataset<float>("/run/histogram_data_1/corrected_time", time, {length}, H5::PredType::NATIVE_FLOAT);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/corrected_time", "axis", 1);
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/corrected_time", "units", std::string("microseconds"));
nxs->AddDatasetAttribute<float>("/run/histogram_data_1/corrected_time", "available", 0);
// set grouping
std::vector<int> grouping(fData[0].GetNoOfHistos(), 0);
nxs->AddDataset<int>("/run/histogram_data_1/grouping", grouping, {grouping.size()}, H5::PredType::NATIVE_INT32);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/grouping", "avaliabe", 0);
// set alpha
int ival=1;
nxs->AddDataset<int>("/run/histogram_data_1/alpha", {ival}, {1}, H5::PredType::NATIVE_INT32);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/alpha", "avaliabe", 0);
// set counts
// feed histos
PIntVector data;
UInt_t size = 0;
int noHisto, histoLength;
if (fAny2ManyInfo->rebin == 1) {
noHisto = fData[0].GetNoOfHistos();
for (UInt_t i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
size = dataSet->GetData()->size();
histoLength = size;
for (UInt_t j=0; j<size; j++) {
data.push_back((UInt_t)dataSet->GetData()->at(j));
}
}
} else { // rebin > 1
UInt_t dataRebin = 0;
UInt_t dataCount = 0;
noHisto = fData[0].GetNoOfHistos();
for (UInt_t i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
size = dataSet->GetData()->size();
dataCount = 0;
for (UInt_t j=0; j<size; j++) {
if ((j > 0) && (j % fAny2ManyInfo->rebin == 0)) {
dataCount++;
data.push_back(dataRebin);
dataRebin = 0;
}
dataRebin += static_cast<UInt_t>(dataSet->GetData()->at(j));
}
}
size = dataCount;
}
nxs->AddDataset<int>("/run/histogram_data_1/counts", data, {(long unsigned int)noHisto, size}, H5::PredType::NATIVE_INT32);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "units", std::string("counts"));
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "signal", 1);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "number", noHisto);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "length", size);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "t0_bin", 0);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "first_good_bin", 0);
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "last_good_bin", 0);
res = (float)(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin*1.0e3)/2.0;
nxs->AddDatasetAttribute<int>("/run/histogram_data_1/counts", "offset", res);
} else { // IDF V2
nxs->AddGroupAttribute("/raw_data_1", "NX_class", std::string("NXentry"));
// set IDF version
nxs->AddDataset<int>("/raw_data_1/IDF_version", {(int)fAny2ManyInfo->idf}, {1}, H5::PredType::NATIVE_INT32);
// set beamline
str = *fData[0].GetBeamline();
nxs->AddDataset<std::string>("/raw_data_1/beamline", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set definition
nxs->AddDataset<std::string>("/raw_data_1/definition", {std::string("muonTD")}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set run_number
nxs->AddDataset<int>("/raw_data_1/run_number", {fData[0].GetRunNumber()}, {1}, H5::PredType::NATIVE_INT32);
// set title
nxs->AddDataset<std::string>("/raw_data_1/title", {fData[0].GetRunTitle()->Data()}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set start time
str = std::string(fData[0].GetStartDate()->Data()) + std::string("T") + std::string(fData[0].GetStartTime()->Data());
nxs->AddDataset<std::string>("/raw_data_1/start_time", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
nxs->AddDatasetAttribute<std::string>("/raw_data_1/start_time", "units", "ISO8601");
// set end time
str = std::string(fData[0].GetStopDate()->Data()) + std::string("T") + std::string(fData[0].GetStopTime()->Data());
nxs->AddDataset<std::string>("/raw_data_1/end_time", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
nxs->AddDatasetAttribute<std::string>("/raw_data_1/end_time", "units", "ISO8601");
// set experiment_identifier
str = "n/a";
nxs->AddDataset<std::string>("/raw_data_1/experiment_identifier", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set instrument attribute
nxs->AddGroupAttribute("/raw_data_1/instrument", "NX_class", std::string("NXinstrument"));
// set instrument name
str = *fData[0].GetInstrument();
nxs->AddDataset<std::string>("/raw_data_1/instrument/name", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set instrument/source attribute
nxs->AddGroupAttribute("/raw_data_1/instrument/source", "NX_class", std::string("NXsource"));
// set instrument/source/name
str = fData[0].GetLaboratory()->Data();
nxs->AddDataset<std::string>("/raw_data_1/instrument/source/name", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set instrument/source/type
TString tstr = *fData[0].GetInstrument();
std::string type{"n/a"};
if (tstr.Contains("LEM", TString::kIgnoreCase)) {
type = "low energy muon source";
} else if (tstr.Contains("GPS", TString::kIgnoreCase) || tstr.Contains("GPD", TString::kIgnoreCase) ||
tstr.Contains("LTF", TString::kIgnoreCase) || tstr.Contains("FLAME", TString::kIgnoreCase) ||
tstr.Contains("HAL-9500", TString::kIgnoreCase) || tstr.Contains("DOLLY", TString::kIgnoreCase) ||
tstr.Contains("VMS", TString::kIgnoreCase)) {
type = "quasi-continous muon source";
} else if (tstr.Contains("EMU", TString::kIgnoreCase) || tstr.Contains("MUSR", TString::kIgnoreCase) ||
tstr.Contains("HIFI", TString::kIgnoreCase)) {
type = "pulsed muon source";
}
nxs->AddDataset<std::string>("/raw_data_1/instrument/source/type", {type}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set instrument/source/probe
str = "positive muons";
nxs->AddDataset<std::string>("/raw_data_1/instrument/source/probe", {str}, {1}, H5::StrType(H5::PredType::C_S1, H5T_VARIABLE));
// set instrument/detector info
nxs->AddGroupAttribute("/raw_data_1/instrument/detector_1", "NX_class", std::string("NXdetector"));
// set instrument/detector/spectrum_index
int noHistos = fData[0].GetNoOfHistos();
std::vector<int> spectrum_index(noHistos);
for (unsigned int i=0; i<spectrum_index.size(); i++)
spectrum_index[i] = i+1;
nxs->AddDataset<int>("/raw_data_1/instrument/detector_1/spectrum_index", spectrum_index, {spectrum_index.size()}, H5::PredType::NATIVE_INT32);
// set instrument/detector/raw_time (not useful for quasi-continuous sources)
int ival=0;
nxs->AddDataset<int>("/raw_data_1/instrument/detector_1/raw_time", {ival}, {1}, H5::PredType::NATIVE_INT32);
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/raw_time", "available", 0);
// set instrument/detector/counts
// set counts
// feed histos
PRawRunDataSet *dataSet = nullptr;
PIntVector data;
UInt_t size = 0;
int noHisto, histoLength;
if (fAny2ManyInfo->rebin == 1) {
noHisto = fData[0].GetNoOfHistos();
for (UInt_t i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
size = dataSet->GetData()->size();
histoLength = size;
for (UInt_t j=0; j<size; j++) {
data.push_back((UInt_t)dataSet->GetData()->at(j));
}
}
} else { // rebin > 1
UInt_t dataRebin = 0;
UInt_t dataCount = 0;
noHisto = fData[0].GetNoOfHistos();
for (UInt_t i=0; i<fData[0].GetNoOfHistos(); i++) {
dataSet = fData[0].GetDataSet(i, false); // i.e. the false means, that i is the index and NOT the histo number
if (dataSet == nullptr) { // something is really wrong
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile: **ERROR** Couldn't get data set (idx=" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
size = dataSet->GetData()->size();
dataCount = 0;
for (UInt_t j=0; j<size; j++) {
if ((j > 0) && (j % fAny2ManyInfo->rebin == 0)) {
dataCount++;
data.push_back(dataRebin);
dataRebin = 0;
}
dataRebin += static_cast<UInt_t>(dataSet->GetData()->at(j));
}
}
size = dataCount;
}
nxs->AddDataset<int>("/raw_data_1/instrument/detector_1/counts", data, {(long unsigned int)noHisto, size}, H5::PredType::NATIVE_INT32);
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "axes", std::string("spectrum_index,time_bin"));
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "long_name", std::string("positon counts"));
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "signal", 1);
// t0_bin attributes
std::vector<int> t0_bin;
for (unsigned int i=0; i<fData[0].GetNoOfHistos(); i++)
t0_bin.push_back((int)(fData[0].GetT0Bin(i+1)/fAny2ManyInfo->rebin));
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "t0_bin", {t0_bin});
nxs->AddDatasetAttribute<int>("/raw_data_1/instrument/detector_1/counts", "units", std::string("counts"));
// set raw_data_1/detector info
nxs->AddGroupAttribute("/raw_data_1/detector_1", "NX_class", std::string("NXdata"));
// set detector/counts
nxs->AddDataset<int>("/raw_data_1/detector_1/counts", data, {(long unsigned int)noHisto, size}, H5::PredType::NATIVE_INT32);
nxs->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "axes", std::string("spectrum_index,time_bin"));
nxs->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "long_name", std::string("positon counts"));
nxs->AddDatasetAttribute<int>("/raw_data_1/detector_1/counts", "signal", 1);
}
int result = nxs->WriteNexusFile(fln.Data(), fAny2ManyInfo->idf);
if (result != 0) {
std::cerr << std::endl << "**ERROR** PRunDataHandler::WriteNexusFile, fln=" << fln << std::endl;
return false;
}
} catch (const H5::Exception& e) {
std::cerr << std::endl << "HDF5 error: " << e.getDetailMsg() << std::endl;
}
}
// set IDF version
src/external/nexus/PNeXus.cpp
Vendored
+5679 -5343
View File
File diff suppressed because it is too large Load Diff
src/external/nexus/PNeXus.h
Vendored
+2016 -558
View File
File diff suppressed because it is too large Load Diff