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

adopted to new hdf4/5 approach of NeXus.

Browse Source
This commit is contained in:
suter_a
2026-03-02 13:46:47 +01:00
parent c0ad630740
commit 8cb7da2e2a
13 changed files with 672 additions and 673 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/classes/CMakeLists.txt
+6 -2
View File
@@ -6,6 +6,7 @@ set(MUSRFIT_INC ${CMAKE_SOURCE_DIR}/src/include)
# Hence, target_include_directories cannot be used here because, targets are
# setup only afterwards.
include_directories(${MUSRFIT_INC})
include_directories(${FFTW3_INC})
root_generate_dictionary(
PFourierCanvasDict
@@ -81,7 +82,7 @@ set(prefix "${CMAKE_INSTALL_PREFIX}")
set(exec_prefix "\$\{prefix\}")
set(libdir "\$\{exec_prefix\}/lib")
set(includedir "\$\{prefix\}/include")
set(MUSR_VERSION "1.7.0")
set(MUSR_VERSION "1.8.0")
set(MUSR_LIBRARY_NAME "PMusr")
configure_file("PMusr.pc.in" "PMusr.pc" @ONLY)
set(USERFCN_LIBRARY_NAME "PUserFcnBase")
@@ -132,13 +133,15 @@ target_include_directories(
$<BUILD_INTERFACE:${DKS_INCLUDE_DIR}>
$<BUILD_INTERFACE:${Boost_INCLUDE_DIRS}>
$<BUILD_INTERFACE:${FFTW3_INCLUDE}>
$<BUILD_INTERFACE:${NEXUS_INCLUDE_DIR}>
$<BUILD_INTERFACE:${HDF4_INCLUDE_DIRS}>
$<BUILD_INTERFACE:${HDF5_INCLUDE_DIRS}>
$<BUILD_INTERFACE:${MUSRFIT_INC}>
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src/external/MusrRoot>
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src/external/TLemRunHeader>
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src/external/MuSR_software/Class_MuSR_PSI>
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src/external/mud/src>
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src/external/nexus>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}>
)
add_library(PUserFcnBase SHARED
@@ -152,6 +155,7 @@ add_library(PRgeHandler SHARED
)
#--- set target properties, e.g. version --------------------------------------
target_compile_options(PMusr BEFORE PRIVATE "-DHAVE_CONFIG_H")
set_target_properties(PMusr
PROPERTIES
VERSION ${MUSR_VERSION}
src/classes/PMsrHandler.cpp
+42 -4
View File
@@ -814,6 +814,10 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
fout.width(16);
fout << std::left << "packing";
fout << fGlobal.GetPacking() << std::endl;
} else if (sstr.BeginsWith("deadtime-cor")) {
fout.width(16);
fout << std::left << "deadtime-cor";
fout << fGlobal.GetDeadTimeCorrection() << std::endl;
} else {
fout << str.Data() << std::endl;
}
@@ -1157,6 +1161,10 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
fout.width(16);
fout << std::left << "packing";
fout << fRuns[runNo].GetPacking() << std::endl;
} else if (sstr.BeginsWith("deadtime-cor")) {
fout.width(16);
fout << std::left << "deadtime-cor";
fout << fRuns[runNo].GetDeadTimeCorrection() << std::endl;
} else {
fout << str.Data() << std::endl;
}
@@ -3334,6 +3342,20 @@ Bool_t PMsrHandler::HandleGlobalEntry(PMsrLines &lines)
error = true;
}
}
} else if (iter->fLine.BeginsWith("deadtime-cor", TString::kIgnoreCase)) { // deadtime correction
if (tokens->GetEntries() < 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (!str.CompareTo("no", TString::kIgnoreCase) ||
!str.CompareTo("file", TString::kIgnoreCase) ||
!str.CompareTo("estimate", TString::kIgnoreCase)) {
global.SetDeadTimeCorrection(str);
} else {
error = true;
}
}
}
// clean up
@@ -3931,6 +3953,26 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
}
}
// deadtime-correction -----------------------------------
if (iter->fLine.BeginsWith("deadtime-cor", TString::kIgnoreCase)) { // deadtime correction
runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following
if (tokens->GetEntries() < 2) {
error = true;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (!str.CompareTo("no", TString::kIgnoreCase) ||
!str.CompareTo("file", TString::kIgnoreCase) ||
!str.CompareTo("estimate", TString::kIgnoreCase)) {
param.SetDeadTimeCorrection(str);
} else {
error = true;
}
}
}
// xy-data -----------------------------------------------
if (line.BeginsWith("xy-data", TString::kIgnoreCase)) {
@@ -6413,10 +6455,6 @@ Bool_t PMsrHandler::CheckMaps()
for (UInt_t i=0; i<map->size(); i++)
if (map->at(i) != 0)
fNoOfMaps++;
/*as
if (fNoOfMaps == 0)
fNoOfMaps = -1;
*/
}
// clean up
src/classes/PMusr.cpp
+16
View File
@@ -721,6 +721,22 @@ PRawRunDataSet* PRawRunData::GetDataSet(const UInt_t idx, Bool_t wantHistoNo)
return fData.GetSet(idx);
}
//--------------------------------------------------------------------------
// DeadTimeCorrectionReady (public)
//--------------------------------------------------------------------------
/**
* <p>Checks if deadtime correction information is sufficient to apply it.
* This means that fNumberOfGoodFrames must be present and the deadtime
* parameter vector.
*
* @return true if ready to apply deadtime correctio, false otherwise
*/
const Bool_t PRawRunData::DeadTimeCorrectionReady()
{
if ((fNumberOfGoodFrames > 0) && (fDeadTimeParam.size() > 0))
return true;
return false;
}
//--------------------------------------------------------------------------
// SetRingAnode (public)
src/classes/PRunDataHandler.cpp
+96 -640
View File
@@ -408,9 +408,9 @@ Bool_t PRunDataHandler::WriteData(TString fileName)
break;
case A2M_NEXUS:
if (fAny2ManyInfo->outFileName.Length() == 0)
success = WriteNexusFile(fileName);
success = WriteNexusFile(fAny2ManyInfo->outFormat, fileName);
else
success = WriteNexusFile(fAny2ManyInfo->outFileName);
success = WriteNexusFile(fAny2ManyInfo->outFormat, fAny2ManyInfo->outFileName);
break;
default:
break;
@@ -684,9 +684,9 @@ Bool_t PRunDataHandler::ReadWriteFilesList()
break;
case A2M_NEXUS:
if (fAny2ManyInfo->outFileName.Length() == 0)
success = WriteNexusFile();
success = WriteNexusFile(fAny2ManyInfo->outFormat);
else
success = WriteNexusFile(fAny2ManyInfo->outFileName);
success = WriteNexusFile(fAny2ManyInfo->outFormat, fAny2ManyInfo->outFileName);
break;
case A2M_MUD:
if (fAny2ManyInfo->outFileName.Length() == 0)
@@ -787,7 +787,7 @@ Bool_t PRunDataHandler::ReadWriteFilesList()
success = WritePsiBinFile(fln);
break;
case A2M_NEXUS:
success = WriteNexusFile(fln);
success = WriteNexusFile(fAny2ManyInfo->outFormat, fln);
break;
case A2M_MUD:
success = WriteMudFile(fln);
@@ -2193,348 +2193,98 @@ Bool_t PRunDataHandler::ReadRootFile()
Bool_t PRunDataHandler::ReadNexusFile()
{
#ifdef PNEXUS_ENABLED
std::cout << std::endl << ">> PRunDataHandler::ReadNexusFile(): Will read nexus file " << fRunPathName.Data() << " ...";
std::cout << std::endl << ">> PRunDataHandler::ReadNexusFile(): Will read nexus file " << fRunPathName.Data() << " ...";
PDoubleVector histoData;
PRawRunData runData;
PRawRunDataSet dataSet;
TString str;
std::string sstr;
Double_t dval;
bool ok;
nxs::HDFType type = nxs::checkHDFType(fRunPathName.Data());
std::unique_ptr<PNeXus> nxs_file = std::make_unique<PNeXus>(fRunPathName.Data());
if (!nxs_file->IsValid()) {
std::cerr << std::endl << ">> PRunDataHandler::ReadNexusFile(): Not a valid NeXus file.";
std::cerr << std::endl << ">> Error Message: " << nxs_file->GetErrorMsg().data() << std::endl;
return false;
}
if (nxs_file->GetIdfVersion() == 1) {
if (!nxs_file->IsValid()) {
std::cout << std::endl << "**ERROR** invalid NeXus IDF 2 version file found." << std::endl;
return false;
// check for type errors, missing enabled HDF4
switch (type) {
case nxs::HDFType::HDF4:
std::cout << std::endl << ">> PRunDataHandler::ReadNexusFile(): HDF4 file." << std::endl;
#ifndef HAVE_HDF4
std::cerr << std::endl << ">> PRunDataHandler::ReadNexusFile(): **ERROR**, HDF4 is not enabled." << std::endl;
return false;
#endif
break;
case nxs::HDFType::HDF5:
std::cout << std::endl << ">> PRunDataHandler::ReadNexusFile(): HDF5 file." << std::endl;
break;
case nxs::HDFType::Unknown:
std::cerr << std::endl << ">> PRunDataHandler::ReadNexusFile(): Not a valid NeXus file." << std::endl;
return false;
}
// get header information
// get/set laboratory
str = TString(nxs_file->GetEntryIdf1()->GetLaboratory());
runData.SetLaboratory(str);
// get/set beamline
str = TString(nxs_file->GetEntryIdf1()->GetBeamline());
runData.SetBeamline(str);
// get/set instrument
str = TString(nxs_file->GetEntryIdf1()->GetInstrument()->GetName());
runData.SetInstrument(str);
// get/set run title
str = TString(nxs_file->GetEntryIdf1()->GetTitle());
runData.SetRunTitle(str);
// get/set run number
runData.SetRunNumber(nxs_file->GetEntryIdf1()->GetRunNumber());
// get/set temperature
dval = nxs_file->GetEntryIdf1()->GetSample()->GetPhysPropValue("temperature", ok);
if (ok)
runData.SetTemperature(0, dval, 0.0);
// get/set field
dval = nxs_file->GetEntryIdf1()->GetSample()->GetPhysPropValue("magnetic_field", ok);
nxs_file->GetEntryIdf1()->GetSample()->GetPhysPropUnit("magnetic_field", sstr, ok);
str = sstr;
// since field has to be given in Gauss, check the units
Double_t factor=1.0;
if (!str.CompareTo("gauss", TString::kIgnoreCase))
factor=1.0;
else if (!str.CompareTo("tesla", TString::kIgnoreCase))
factor=1.0e4;
else
factor=1.0;
runData.SetField(factor*dval);
// get/set implantation energy
runData.SetEnergy(PMUSR_UNDEFINED);
// get/set moderator HV
runData.SetTransport(PMUSR_UNDEFINED);
// get/set RA HV's (LEM specific)
for (UInt_t i=0; i<4; i++)
runData.SetRingAnode(i, PMUSR_UNDEFINED);
// get/set setup
runData.SetSetup(nxs_file->GetEntryIdf1()->GetNotes());
// get/set sample
runData.SetSample(nxs_file->GetEntryIdf1()->GetSample()->GetName());
// get/set orientation
runData.SetOrientation("??");
// get/set time resolution (ns)
runData.SetTimeResolution(nxs_file->GetEntryIdf1()->GetData()->GetTimeResolution("ns"));
// get/set start/stop time
sstr = nxs_file->GetEntryIdf1()->GetStartTime();
str = sstr;
TString date, time;
SplitTimeDate(str, time, date, ok);
if (ok) {
runData.SetStartTime(time);
runData.SetStartDate(date);
}
sstr = nxs_file->GetEntryIdf1()->GetStopTime();
str = sstr;
SplitTimeDate(str, time, date, ok);
if (ok) {
runData.SetStopTime(time);
runData.SetStopDate(date);
}
// get/set t0, firstGoodBin, lastGoodBin
std::vector<unsigned int> *t0 = nxs_file->GetEntryIdf1()->GetData()->GetT0s();
std::vector<unsigned int> *fgb = nxs_file->GetEntryIdf1()->GetData()->GetFirstGoodBins();
std::vector<unsigned int> *lgb = nxs_file->GetEntryIdf1()->GetData()->GetLastGoodBins();
// get/set data
std::vector<unsigned int> *pdata;
unsigned int max=0, binMax=0;
PDoubleVector data;
for (UInt_t i=0; i<nxs_file->GetEntryIdf1()->GetData()->GetNoOfHistos(); i++) {
pdata = nxs_file->GetEntryIdf1()->GetData()->GetHisto(i);
for (UInt_t j=0; j<pdata->size(); j++) {
data.push_back(pdata->at(j));
if (pdata->at(j) > max) {
max = pdata->at(j);
binMax = j;
}
}
// fill data set
dataSet.Clear();
dataSet.SetHistoNo(i+1); // i.e. histo numbers start with 1
// set time zero bin
if (i<t0->size())
dataSet.SetTimeZeroBin(t0->at(i));
else
dataSet.SetTimeZeroBin(t0->at(0));
// set time zero bin estimate
dataSet.SetTimeZeroBinEstimated(binMax);
// set first good bin
if (i<fgb->size())
dataSet.SetFirstGoodBin(fgb->at(i));
else
dataSet.SetFirstGoodBin(fgb->at(0));
// set last good bin
if (i<lgb->size())
dataSet.SetFirstGoodBin(lgb->at(i));
else
dataSet.SetFirstGoodBin(lgb->at(0));
dataSet.SetData(data);
runData.SetDataSet(dataSet);
data.clear();
}
// keep run name from the msr-file
runData.SetRunName(fRunName);
// keep the information
fData.push_back(runData);
} else if (nxs_file->GetIdfVersion() == 2) {
if (!nxs_file->IsValid()) {
std::cout << std::endl << "**ERROR** invalid NeXus IDF 2 version file found." << std::endl;
return false;
}
// get header information
// get/set laboratory
str = TString(nxs_file->GetEntryIdf2()->GetInstrument()->GetSource()->GetName());
runData.SetLaboratory(str);
// get/set beamline
str = TString(nxs_file->GetEntryIdf2()->GetInstrument()->GetName());
runData.SetBeamline(str);
// get/set instrument
str = TString(nxs_file->GetEntryIdf2()->GetInstrument()->GetName());
runData.SetInstrument(str);
// get/set muon source
str = TString(nxs_file->GetEntryIdf2()->GetInstrument()->GetSource()->GetType());
runData.SetMuonSource(str);
// get/set muon species
str = TString(nxs_file->GetEntryIdf2()->GetInstrument()->GetSource()->GetProbe());
runData.SetMuonSpecies(str);
// get/set run title
str = TString(nxs_file->GetEntryIdf2()->GetTitle());
runData.SetRunTitle(str);
// get/set run number
runData.SetRunNumber(nxs_file->GetEntryIdf2()->GetRunNumber());
// get/set temperature
dval = nxs_file->GetEntryIdf2()->GetSample()->GetPhysPropValue("temperature", ok);
if (ok)
runData.SetTemperature(0, dval, 0.0);
// get/set field
dval = nxs_file->GetEntryIdf2()->GetSample()->GetPhysPropValue("magnetic_field", ok);
nxs_file->GetEntryIdf2()->GetSample()->GetPhysPropUnit("magnetic_field", sstr, ok);
str = sstr;
// since field has to be given in Gauss, check the units
Double_t factor=1.0;
if (!str.CompareTo("gauss", TString::kIgnoreCase))
factor=1.0;
else if (!str.CompareTo("tesla", TString::kIgnoreCase))
factor=1.0e4;
else
factor=1.0;
runData.SetField(factor*dval);
// get/set implantation energy
runData.SetEnergy(PMUSR_UNDEFINED);
// get/set moderator HV
runData.SetTransport(PMUSR_UNDEFINED);
// get/set RA HV's (LEM specific)
for (UInt_t i=0; i<4; i++)
runData.SetRingAnode(i, PMUSR_UNDEFINED);
// get/set setup take NXsample/temperature_1_env and NXsample/magnetic_field_1_env
sstr = nxs_file->GetEntryIdf2()->GetSample()->GetEnvironmentTemp() + std::string("/");
sstr += nxs_file->GetEntryIdf2()->GetSample()->GetEnvironmentField();
str = sstr;
runData.SetSetup(str);
// get/set sample
runData.SetSample(nxs_file->GetEntryIdf2()->GetSample()->GetName());
// get/set orientation
runData.SetOrientation("??");
// get/set time resolution (ns)
runData.SetTimeResolution(nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetTimeResolution("ns"));
// get/set start/stop time
sstr = nxs_file->GetEntryIdf2()->GetStartTime();
str = sstr;
TString date, time;
SplitTimeDate(str, time, date, ok);
if (ok) {
runData.SetStartTime(time);
runData.SetStartDate(date);
}
sstr = nxs_file->GetEntryIdf2()->GetStopTime();
str = sstr;
SplitTimeDate(str, time, date, ok);
if (ok) {
runData.SetStopTime(time);
runData.SetStopDate(date);
}
// get/set data, t0, fgb, lgb
PDoubleVector data;
PDoubleVector histoData;
PRawRunData runData;
PRawRunDataSet dataSet;
UInt_t histoNo = 0;
Int_t ival;
int *histos = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetHistos();
if (nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfPeriods() > 0) { // counts[][][]
for (int i=0; i<nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfPeriods(); i++) {
for (int j=0; j<nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfSpectra(); j++) {
for (int k=0; k<nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfBins(); k++) {
data.push_back(*(histos+i*nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfSpectra()+j*nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfBins()+k));
}
dataSet.Clear();
dataSet.SetHistoNo(++histoNo); // i.e. histo numbers start with 1
// get t0
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetT0(i,j);
if (ival == -1) // i.e. single value only
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetT0();
dataSet.SetTimeZeroBin(ival);
// get first good bin
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetFirstGoodBin(i,j);
if (ival == -1) // i.e. single value only
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetFirstGoodBin();
dataSet.SetFirstGoodBin(ival);
// get last good bin
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetLastGoodBin(i,j);
if (ival == -1) // i.e. single value only
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetLastGoodBin();
dataSet.SetLastGoodBin(ival);
TString str;
std::string sstr;
Int_t ival, idf{-1};
Double_t dval, factor;
bool ok;
dataSet.SetData(data);
runData.SetDataSet(dataSet);
data.clear();
if (type == nxs::HDFType::HDF4) {
#ifdef HAVE_HDF4
std::unique_ptr<nxH4::PNeXus> nxs_file = std::make_unique<nxH4::PNeXus>(fRunPathName.Data());
if (nxs_file == nullptr) {
std::cerr << std::endl << "**ERROR** allocation of nxH4::PNeXus object failed." << std::endl;
return true;
}
// check for IDF_version
if (nxs_file->HasDataset("/run/IDF_version")) {
idf = nxs_file->GetDataset<int>("/run/IDF_version").GetData()[0];
std::cout << ">> PRunDataHandler::ReadNexusFile(): IDF V" << idf << std::endl;
}
if (idf == -1) { // IDF_version not found
if (nxs_file->HasDataset("/raw_data_1/IDF_version")) {
idf = nxs_file->GetDataset<int>("/raw_data_1/IDF_version").GetData()[0];
std::cout << ">> PRunDataHandler::ReadNexusFile(): IDF V" << idf << std::endl;
}
}
} else {
if (nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfSpectra() > 0) { // counts[][]
for (int i=0; i<nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfSpectra(); i++) {
for (int j=0; j<nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfBins(); j++) {
data.push_back(*(histos+i*nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfBins()+j));
}
dataSet.Clear();
dataSet.SetHistoNo(++histoNo); // i.e. histo numbers start with 1
// get t0
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetT0(-1,i);
if (ival == -1) // i.e. single value only
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetT0();
dataSet.SetTimeZeroBin(ival);
// get first good bin
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetFirstGoodBin(-1,i);
if (ival == -1) // i.e. single value only
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetFirstGoodBin();
dataSet.SetFirstGoodBin(ival);
// get last good bin
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetLastGoodBin(-1,i);
if (ival == -1) // i.e. single value only
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetLastGoodBin();
dataSet.SetLastGoodBin(ival);
if ((idf != 1) && (idf != 2)) {
std::cerr << std::endl << ">> PRunDataHandler::ReadNexusFile(): a NeXus file with an invalid IDF V" << idf << std::endl;
return false;
}
dataSet.SetData(data);
runData.SetDataSet(dataSet);
data.clear();
}
} else { // counts[]
for (int i=0; i<nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfBins(); i++) {
data.push_back(*(histos+i));
}
dataSet.Clear();
dataSet.SetHistoNo(++histoNo); // i.e. histo numbers start with 1
// get t0
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetT0();
dataSet.SetTimeZeroBin(ival);
// get first good bin
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetFirstGoodBin();
dataSet.SetFirstGoodBin(ival);
// get last good bin
ival = nxs_file->GetEntryIdf2()->GetInstrument()->GetDetector()->GetLastGoodBin();
dataSet.SetLastGoodBin(ival);
if (idf == 1) { // HDF4 IDF V1
if (!ReadNexusFileIdf1(nxs_file))
return false;
} else { // HDF4 IDF V2
// not yet implemented
}
#endif
} else { // HDF5
std::unique_ptr<nxH5::PNeXus> nxs_file = std::make_unique<nxH5::PNeXus>(fRunPathName.Data());
if (nxs_file == nullptr) {
std::cerr << std::endl << "**ERROR** allocation of nxH5::PNeXus object failed." << std::endl;
return true;
}
dataSet.SetData(data);
runData.SetDataSet(dataSet);
data.clear();
// check for IDF_version
if (nxs_file->HasDataset("/run/IDF_version")) {
idf = nxs_file->GetDataset<int>("/run/IDF_version").GetData()[0];
std::cout << ">> PRunDataHandler::ReadNexusFile(): IDF V" << idf << std::endl;
}
if (idf == -1) { // IDF_version not found
if (nxs_file->HasDataset("/raw_data_1/IDF_version")) {
idf = nxs_file->GetDataset<int>("/raw_data_1/IDF_version").GetData()[0];
std::cout << ">> PRunDataHandler::ReadNexusFile(): IDF V" << idf << std::endl;
}
}
if ((idf != 1) && (idf != 2)) {
std::cerr << std::endl << ">> PRunDataHandler::ReadNexusFile(): a NeXus file with an invalid IDF V" << idf << std::endl;
return false;
}
if (idf == 1) { // HDF5 IDF V1
if (!ReadNexusFileIdf1(nxs_file))
return false;
} else { // HDF5 IDF V2
if (!ReadNexusFileIdf2(nxs_file))
return false;
}
}
// keep run name from the msr-file
runData.SetRunName(fRunName);
// keep the information
fData.push_back(runData);
} else {
std::cout << std::endl << ">> PRunDataHandler::ReadNexusFile(): IDF version " << nxs_file->GetIdfVersion() << ", not implemented." << std::endl;
}
#else
std::cout << std::endl << ">> PRunDataHandler::ReadNexusFile(): Sorry, not enabled at configuration level, i.e. --enable-NeXus when executing configure" << std::endl << std::endl;
#endif
@@ -4956,7 +4706,7 @@ Bool_t PRunDataHandler::WriteRootFile(TString fln)
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=0" << i << ")";
std::cerr << std::endl << ">> PRunDataHandler::WriteRootFile: **ERROR** Couldn't get data set (idx=0" << i << ")";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
@@ -5050,15 +4800,14 @@ Bool_t PRunDataHandler::WriteRootFile(TString fln)
/**
* <p> Write the nexus-file format.
*
* <b>return:</b>
* - true on successful writting,
* - otherwise false.
*
* \param fln file name. If empty, the routine will try to construct one
* \return true on successful writting, otherwise false.
*/
Bool_t PRunDataHandler::WriteNexusFile(TString fln)
Bool_t PRunDataHandler::WriteNexusFile(TString format, TString fln)
{
#ifdef PNEXUS_ENABLED
std::string str{""};
Bool_t ok = false;
fln = GenerateOutputFileName(fln, ".nxs", ok);
if (!ok)
@@ -5097,9 +4846,9 @@ Bool_t PRunDataHandler::WriteNexusFile(TString fln)
// set program name
nxs->AddDataset<std::string>("/run/program_name", {"any2many"}, {1}, nxH4::H4DataType::kCHAR8);
str="n/a";
#ifdef HAVE_CONFIG_H
#ifdef HAVE_CONFIG_H
str = PACKAGE_VERSION;
#endif
#endif
nxs->AddDatasetAttribute<std::string>("/run/program_name", "version", str);
// set run number
@@ -5471,9 +5220,9 @@ Bool_t PRunDataHandler::WriteNexusFile(TString fln)
// 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
#ifdef HAVE_CONFIG_H
str = PACKAGE_VERSION;
#endif
#endif
nxs->AddDatasetAttribute<std::string>("/run/program_name", "version", str);
// set run number
@@ -5816,302 +5565,9 @@ Bool_t PRunDataHandler::WriteNexusFile(TString fln)
}
}
// set IDF version
nxs->SetIdfVersion(fAny2ManyInfo->idf);
if (fAny2ManyInfo->idf == 1) {
// fill necessary data structures
nxs->SetFileName(fln.Data());
// set file creating time
time_t now;
struct tm *tm;
time(&now);
tm = localtime(&now);
std::string str("");
char cstr[128];
strftime(cstr, sizeof(cstr), "%FT%T", tm);
str = std::string(cstr);
nxs->SetFileTime(str);
nxs->GetEntryIdf1()->SetProgramName("any2many");
nxs->GetEntryIdf1()->SetProgramVersion("$Id$");
nxs->GetEntryIdf1()->SetRunNumber(fData[0].GetRunNumber());
nxs->GetEntryIdf1()->SetTitle(fData[0].GetRunTitle()->Data());
nxs->GetEntryIdf1()->SetNotes("n/a");
nxs->GetEntryIdf1()->SetAnalysis("muonTD");
if (*fData[0].GetLaboratory() != "n/a")
nxs->GetEntryIdf1()->SetLaboratory(fData[0].GetLaboratory()->Data());
if (*fData[0].GetBeamline() != "n/a")
nxs->GetEntryIdf1()->SetBeamline(fData[0].GetBeamline()->Data());
str = std::string(fData[0].GetStartDate()->Data()) + std::string("T") + std::string(fData[0].GetStartTime()->Data());
nxs->GetEntryIdf1()->SetStartTime(str);
str = std::string(fData[0].GetStopDate()->Data()) + std::string("T") + std::string(fData[0].GetStopTime()->Data());
nxs->GetEntryIdf1()->SetStopTime(str);
nxs->GetEntryIdf1()->SetSwitchingState(1);
nxs->GetEntryIdf1()->GetUser()->SetName("n/a");
nxs->GetEntryIdf1()->GetUser()->SetExperimentNumber("n/a");
nxs->GetEntryIdf1()->GetSample()->SetName(fData[0].GetSample()->Data());
nxs->GetEntryIdf1()->GetSample()->SetPhysProp("temperature", fData[0].GetTemperature(0), "Kelvin");
nxs->GetEntryIdf1()->GetSample()->SetPhysProp("magnetic_field", fData[0].GetField(), "Gauss");
nxs->GetEntryIdf1()->GetSample()->SetEnvironment(fData[0].GetSetup()->Data());
nxs->GetEntryIdf1()->GetSample()->SetShape("n/a");
nxs->GetEntryIdf1()->GetSample()->SetMagneticFieldVectorAvailable(0);
if (*fData[0].GetInstrument() != "n/a")
nxs->GetEntryIdf1()->GetInstrument()->SetName(fData[0].GetInstrument()->Data());
nxs->GetEntryIdf1()->GetInstrument()->GetDetector()->SetNumber(fData[0].GetNoOfHistos());
nxs->GetEntryIdf1()->GetInstrument()->GetCollimator()->SetType("n/a");
// 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);
}
double total_counts_mev = (double) total_counts / 1.0e6;
nxs->GetEntryIdf1()->GetInstrument()->GetBeam()->SetTotalCounts(total_counts_mev);
nxs->GetEntryIdf1()->GetInstrument()->GetBeam()->SetUnits("Mev");
nxs->GetEntryIdf1()->GetData()->SetTimeResolution(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin, "ns");
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
nxs->GetEntryIdf1()->GetData()->SetT0(static_cast<Int_t>(dataSet->GetTimeZeroBin()/fAny2ManyInfo->rebin), i);
nxs->GetEntryIdf1()->GetData()->SetFirstGoodBin(static_cast<Int_t>(dataSet->GetFirstGoodBin()/fAny2ManyInfo->rebin), i);
nxs->GetEntryIdf1()->GetData()->SetLastGoodBin(static_cast<Int_t>(dataSet->GetLastGoodBin()/fAny2ManyInfo->rebin), i);
}
// feed histos
PUIntVector data;
UInt_t size = 0;
if (fAny2ManyInfo->rebin == 1) {
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();
for (UInt_t j=0; j<size; j++) {
data.push_back((UInt_t)dataSet->GetData()->at(j));
}
nxs->GetEntryIdf1()->GetData()->SetHisto(data, i);
data.clear();
}
} else { // rebin > 1
UInt_t dataRebin = 0;
UInt_t dataCount = 0;
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));
}
nxs->GetEntryIdf1()->GetData()->SetHisto(data, i);
data.clear();
}
}
} else if (fAny2ManyInfo->idf == 2) {
// fill necessary data structures
nxs->SetFileName(fln.Data());
// set file creating time
time_t now;
struct tm *tm;
time(&now);
tm = localtime(&now);
std::string str("");
char cstr[128];
strftime(cstr, sizeof(cstr), "%FT%T", tm);
str = std::string(cstr);
nxs->SetFileTime(str);
// NXroot info
nxs->SetCreator("PSI: any2many");
// NXentry info
nxs->GetEntryIdf2()->SetDefinition("muonTD");
nxs->GetEntryIdf2()->SetProgramName("any2many");
nxs->GetEntryIdf2()->SetProgramVersion("$Id$");
nxs->GetEntryIdf2()->SetRunNumber(fData[0].GetRunNumber());
nxs->GetEntryIdf2()->SetTitle(fData[0].GetRunTitle()->Data());
str = std::string(fData[0].GetStartDate()->Data()) + std::string("T") + std::string(fData[0].GetStartTime()->Data());
nxs->GetEntryIdf2()->SetStartTime(str);
str = std::string(fData[0].GetStopDate()->Data()) + std::string("T") + std::string(fData[0].GetStopTime()->Data());
nxs->GetEntryIdf2()->SetStopTime(str);
nxs->GetEntryIdf2()->SetExperimentIdentifier("n/a");
// NXuser info
nxs->GetEntryIdf2()->GetUser()->SetName("n/a");
// NXsample info
nxs->GetEntryIdf2()->GetSample()->SetName(fData[0].GetSample()->Data());
nxs->GetEntryIdf2()->GetSample()->SetDescription("n/a");
nxs->GetEntryIdf2()->GetSample()->SetPhysProp("temperature_1", fData[0].GetTemperature(0), "Kelvin");
nxs->GetEntryIdf2()->GetSample()->SetPhysProp("magnetic_field_1", fData[0].GetField(), "Gauss");
nxs->GetEntryIdf2()->GetSample()->SetEnvironmentTemp(fData[0].GetSetup()->Data());
nxs->GetEntryIdf2()->GetSample()->SetEnvironmentField("n/a");
// here would be the information for NXinstrument. Currently there are not much information to feed this
nxs->GetEntryIdf2()->GetInstrument()->SetName(fData[0].GetInstrument()->Data());
// NXinstrument/NXsource
nxs->GetEntryIdf2()->GetInstrument()->GetSource()->SetName(fData[0].GetLaboratory()->Data());
nxs->GetEntryIdf2()->GetInstrument()->GetSource()->SetType(fData[0].GetMuonSource()->Data());
nxs->GetEntryIdf2()->GetInstrument()->GetSource()->SetProbe(fData[0].GetMuonSpecies()->Data());
// NXinstrument/NXbeamline
nxs->GetEntryIdf2()->GetInstrument()->GetBeamline()->SetName(fData[0].GetBeamline()->Data());
// NXinstrument/NXdetector
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetDescription(fData[0].GetInstrument()->Data()); // assume that this should be the instrument name
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetNoOfPeriods(0); // currently red/green is not distinguished
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetNoOfSpectra(fData[0].GetNoOfHistos());
PRawRunDataSet *dataSet = fData[0].GetDataSet(0, 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)";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetNoOfBins((unsigned int)(dataSet->GetData()->size() / fAny2ManyInfo->rebin));
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetTimeResolution(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin, "ns");
int *histo = nullptr;
int idx = 0;
if (fAny2ManyInfo->rebin == 1) {
histo = new int[fData[0].GetNoOfHistos()*dataSet->GetData()->size()];
idx = 0;
for (int i=0; i<nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfSpectra(); 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;
}
for (unsigned int j=0; j<dataSet->GetData()->size(); j++) {
*(histo+idx++) = (int) dataSet->GetData()->at(j);
}
}
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetHistos(histo);
// clean up
if (histo) {
delete [] histo;
histo = nullptr;
}
} else { // rebin > 1
histo = new int[fData[0].GetNoOfHistos()*(int)(dataSet->GetData()->size()/fAny2ManyInfo->rebin)];
int counts = 0;
idx = 0;
for (int i=0; i<nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfSpectra(); 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;
}
for (unsigned int j=0; j<dataSet->GetData()->size(); j++) {
if ((j>0) && (j % fAny2ManyInfo->rebin == 0)) {
*(histo+idx++) = counts;
counts = 0;
}
counts += (int) dataSet->GetData()->at(j);
}
}
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetHistos(histo);
// clean up
if (histo) {
delete [] histo;
histo = nullptr;
}
}
// handle spectrum index
for (int i=0; i<nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->GetNoOfSpectra(); i++)
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetSpectrumIndex(i+1);
// handle histogram resolution
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetTimeResolution(fData[0].GetTimeResolution()*fAny2ManyInfo->rebin, "ns");
// handle raw time
std::vector<double> raw_time;
UInt_t size = (unsigned int)(dataSet->GetData()->size() / fAny2ManyInfo->rebin);
for (unsigned int i=0; i<size; i++) {
raw_time.push_back((double)i * fData[0].GetTimeResolution() * fAny2ManyInfo->rebin * 1.0e-3); // since time resolution is given in ns, the factor 1.0e-3 is needed to convert to us
}
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetRawTime(raw_time);
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetRawTimeUnit("micro.second");
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetRawTimeName("time");
raw_time.clear();
// handle t0
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetT0Tag(2); // i.e. t0[#histo] format
int *t0 = new int[fData[0].GetNoOfHistos()];
int *fgb = new int[fData[0].GetNoOfHistos()];
int *lgb = new int[fData[0].GetNoOfHistos()];
if ((t0==0) || (fgb==0) || (lgb==0)) {
std::cerr << std::endl << ">> PRunDataHandler::WriteNeXusFile: **ERROR** Couldn't allocate memory for t0, fgb, lgb" << std::endl;
return false;
}
for (unsigned int i=0; i<fData[0].GetNoOfHistos(); i++) {
PRawRunDataSet *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)";
std::cerr << std::endl << ">> something is really wrong!" << std::endl;
return false;
}
t0[i] = (int)(dataSet->GetTimeZeroBin() / fAny2ManyInfo->rebin);
fgb[i] = (int)(dataSet->GetFirstGoodBin() / fAny2ManyInfo->rebin);
lgb[i] = (int)(dataSet->GetLastGoodBin() / fAny2ManyInfo->rebin);
}
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetT0(t0);
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetFirstGoodBin(fgb);
nxs->GetEntryIdf2()->GetInstrument()->GetDetector()->SetLastGoodBin(lgb);
// clean up
if (t0) delete [] t0;
if (fgb) delete [] fgb;
if (lgb) delete [] lgb;
} else {
return false;
}
// filter out the proper file type, i.e. HDF4, HDF5, or XML
char fileType[32];
memset(fileType, '\0', 32);
if (!fAny2ManyInfo->outFormat.CompareTo("nexus1-hdf4", TString::kIgnoreCase) || !fAny2ManyInfo->outFormat.CompareTo("nexus2-hdf4", TString::kIgnoreCase))
strncpy(fileType, "hdf4", sizeof(fileType));
else if (!fAny2ManyInfo->outFormat.CompareTo("nexus1-hdf5", TString::kIgnoreCase) || !fAny2ManyInfo->outFormat.CompareTo("nexus2-hdf5", TString::kIgnoreCase))
strncpy(fileType, "hdf5", sizeof(fileType));
else if (!fAny2ManyInfo->outFormat.CompareTo("nexus1-xml", TString::kIgnoreCase) || !fAny2ManyInfo->outFormat.CompareTo("nexus2-xml", TString::kIgnoreCase))
strncpy(fileType, "xml", sizeof(fileType));
else {
std::cerr << std::endl << ">> PRunDataHandler::WriteNexusFile(): **ERROR** undefined output NeXus format " << fAny2ManyInfo->outFormat.Data() << " found.";
std::cerr << std::endl << ">> Allowed are: hdf4, hdf5, xml" << std::endl;
return false;
}
// write file
nxs->WriteFile(fln, fileType, fAny2ManyInfo->idf);
return true;
#else
std::cout << std::endl << ">> PRunDataHandler::WriteNexusFile(): Sorry, not enabled at configuration level, i.e. --enable-NeXus when executing configure" << std::endl << std::endl;
std::cout << std::endl << ">> PRunDataHandler::WriteNexusFile(): Sorry, not enabled at configuration level, i.e. -Dnexus=1 when executing configure" << std::endl << std::endl;
#endif
return true;