fix dead time correction for Red/Green (period) NeXus IDF V2 data
Build and Deploy Documentation / build-and-deploy (push) Successful in 34s

Since the Red/Green mode handling was added for NeXus HDF5 IDF V2, the
histogram numbers of the routed histograms carry the period offset, i.e.
histoNo = period*10000 + histo. PRunBase::DeadTimeCorrection() was still
addressing the dead time parameter vector with 'histoNo-1', hence for any
period > 0 it read far beyond the end of the vector. The garbage picked up
there scaled the counts, so every invocation of musrfit gave a different
(sometimes nan) result for the affected run blocks.

The dead time parameters are stored contiguously, period-by-period, and
therefore need to be addressed as period*noOfHistosPerPeriod + histo - 1.
Introduce PERIOD_HISTO_OFFSET for the period encoding, keep the number of
histos per period in PRawRunData, and add PRawRunData::GetDeadTimeParam(histoNo)
which does the mapping and guards against out-of-range access.

Cross-checked against the HDF4 IDF V1 twin of the same run, where the two
periods are flattened into histos 1-192: both readers now yield bit-identical
chisq, with and without dead time correction.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
2026-07-26 11:10:53 +02:00
co-authored by Claude Opus 5
parent 68a8fae1e5
commit 0577d5d211
5 changed files with 60 additions and 3 deletions
+8
View File
@@ -15,6 +15,14 @@ https://bitbucket.org/muonspin/musrfit
Release of V1.12.0, 2026/07/22
==============================
Bug fix:
Fix the dead time correction for Red/Green (period) data of NeXus HDF5 IDF V2 files.
The histogram numbers of the routed histograms carry the period offset (10000),
whereas the dead time parameters are stored contiguously, period-by-period. The
dead time parameter vector was addressed with 'histoNo-1' which, for period > 0,
read beyond the end of the vector. As a consequence any fit of such data gave a
different (random) result at each invocation.
Refactoring:
Replace ROOT's TString::Tokenize() + TObjArray/TObjString token handling with
the dependency-free PStringUtils::Split() across all classes. This reduces the
+34
View File
@@ -700,6 +700,40 @@ const Double_t PRawRunData::GetRingAnode(const UInt_t idx)
return fRingAnode[idx];
}
//--------------------------------------------------------------------------
// GetDeadTimeParam (public)
//--------------------------------------------------------------------------
/**
* <p>Returns the dead time parameter belonging to the histogram number histoNo.
*
* <p>The histogram numbers of Red/Green (period) data carry the period offset,
* i.e. histoNo = period*PERIOD_HISTO_OFFSET + histo, whereas the dead time
* parameters are stored contiguously, period-by-period, and hence need to be
* addressed as period*fNoOfHistosPerPeriod + histo - 1.
*
* <b>return:</b>
* - dead time parameter, if it can be addressed
* - 0.0, otherwise, i.e. no dead time correction will be applied
*
* \param histoNo histogram number as given in the msr-file
*/
const Double_t PRawRunData::GetDeadTimeParam(const UInt_t histoNo)
{
UInt_t idx = histoNo - 1;
if (fNoOfHistosPerPeriod > 0) // Red/Green (period) data
idx = (histoNo / PERIOD_HISTO_OFFSET) * fNoOfHistosPerPeriod + (histoNo % PERIOD_HISTO_OFFSET) - 1;
if (idx >= fDeadTimeParam.size()) {
std::cerr << std::endl << ">> PRawRunData::GetDeadTimeParam: **WARNING** histoNo=" << histoNo << " maps onto idx=" << idx;
std::cerr << " which is out of range [0," << fDeadTimeParam.size() << "[. Will not dead time correct this histogram.";
std::cerr << std::endl;
return 0.0;
}
return fDeadTimeParam[idx];
}
//--------------------------------------------------------------------------
// GetDataSet (public)
//--------------------------------------------------------------------------
+2 -2
View File
@@ -210,9 +210,9 @@ void PRunBase::DeadTimeCorrection(std::vector<PDoubleVector> &histos, PUIntVecto
// Dead time correction: n_true = n_obs / (1 - n_obs * t_dt / (good_frames * dt))
Double_t n_true;
Int_t gf = runData->GetNumberOfGoodFrames();
std::vector<float> t_dt = runData->GetDeadTimeParam();
Double_t t_dt = runData->GetDeadTimeParam(histoNo[i]);
for (UInt_t j=0; j<histos[i].size(); j++) {
n_true = histos[i][j] / (1.0 - histos[i][j] * t_dt[histoNo[i]-1] / (gf * fTimeResolution));
n_true = histos[i][j] / (1.0 - histos[i][j] * t_dt / (gf * fTimeResolution));
histos[i][j] = n_true;
}
}
+13
View File
@@ -264,6 +264,15 @@
/// Offset added to function indices for parameter parsing
#define MSR_PARAM_FUN_OFFSET 20000
//-------------------------------------------------------------
/**
* <p>Offset used to encode the period (Red/Green mode) into the histogram
* number, i.e. histoNo = period*PERIOD_HISTO_OFFSET + histo, with histo
* starting at 1 and period starting at 0. Currently used for NeXus IDF V2
* (ISIS) data which are organized as counts[period][histo][bin].
*/
#define PERIOD_HISTO_OFFSET 10000
//-------------------------------------------------------------
/**
* <p>Fourier transform unit tags.
@@ -880,7 +889,9 @@ class PRawRunData {
virtual const Bool_t DeadTimeCorrectionReady();
virtual const Int_t GetNumberOfGoodFrames() { return fNumberOfGoodFrames; }
virtual const std::vector<float> GetDeadTimeParam() { return fDeadTimeParam; }
virtual const Double_t GetDeadTimeParam(const UInt_t histoNo);
virtual const UInt_t GetNoOfHistos() { return fData.Size(); }
virtual const UInt_t GetNoOfHistosPerPeriod() { return fNoOfHistosPerPeriod; }
virtual PRawRunDataSet* GetDataSet(const UInt_t idx, Bool_t wantHistoNo = true);
virtual const PDoubleVector* GetDataBin(const UInt_t histoNo) { return fData.GetData(histoNo); }
virtual const PNonMusrRawRunData* GetDataNonMusr() { return &fDataNonMusr; }
@@ -923,6 +934,7 @@ class PRawRunData {
virtual void SetRedGreenOffset(PIntVector &ivec) { fRedGreenOffset = ivec; }
virtual void SetNumberOfGoodFrames(Int_t ival) { fNumberOfGoodFrames = ival; }
virtual void SetDeadTimeParam(std::vector<float> dvec) { fDeadTimeParam = dvec; }
virtual void SetNoOfHistosPerPeriod(UInt_t ival) { fNoOfHistosPerPeriod = ival; }
virtual void SetDataSet(PRawRunDataSet &dataSet, UInt_t idx=-1) { fData.Set(dataSet, idx); }
PNonMusrRawRunData fDataNonMusr; ///< keeps all ascii- or db-file info in case of nonMusr fit
@@ -964,6 +976,7 @@ class PRawRunData {
PIntVector fRedGreenOffset; ///< keeps the Red/Green offsets
std::vector<float> fDeadTimeParam; ///< dead time parameter vector needed for pulsed sources
Int_t fNumberOfGoodFrames{0}; ///< needed to correct dead times at pulsed sources
UInt_t fNoOfHistosPerPeriod{0}; ///< number of histos per period (Red/Green mode); 0 means the run has no periods
PRawRunDataVector fData; ///< keeps the histos together with the histo related properties such as T0, first good bin, etc.
};
+3 -1
View File
@@ -876,6 +876,8 @@ Bool_t PRunDataHandler::ReadNexusFileIdf2(T& nxs_file)
noOfPeriods = dims[0];
noOfHistos = dims[1];
histoLength = dims[2];
// needed to properly address period related information like the dead time parameters
runData.SetNoOfHistosPerPeriod(noOfHistos);
// get all necessary attributes
if (count_ds.HasAttribute("t0_bin"))
t0_bin = std::any_cast<int>(count_ds.GetAttribute("t0_bin"));
@@ -893,7 +895,7 @@ Bool_t PRunDataHandler::ReadNexusFileIdf2(T& nxs_file)
for (int i=0; i<noOfPeriods; i++) {
for (int j=0; j<noOfHistos; j++) {
dataSet.Clear();
dataSet.SetHistoNo(i*10000+j+1); // i.e. histo numbers start with 1
dataSet.SetHistoNo(i*PERIOD_HISTO_OFFSET+j+1); // i.e. histo numbers start with 1
dataSet.SetTimeZeroBin(t0_bin);
dataSet.SetFirstGoodBin(fgb);
dataSet.SetLastGoodBin(lgb);