diff --git a/ChangeLog b/ChangeLog index bf509f5f..122587a2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -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 diff --git a/src/classes/PMusr.cpp b/src/classes/PMusr.cpp index 3a928880..19a21174 100644 --- a/src/classes/PMusr.cpp +++ b/src/classes/PMusr.cpp @@ -700,6 +700,40 @@ const Double_t PRawRunData::GetRingAnode(const UInt_t idx) return fRingAnode[idx]; } +//-------------------------------------------------------------------------- +// GetDeadTimeParam (public) +//-------------------------------------------------------------------------- +/** + *

Returns the dead time parameter belonging to the histogram number histoNo. + * + *

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. + * + * return: + * - 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) //-------------------------------------------------------------------------- diff --git a/src/classes/PRunBase.cpp b/src/classes/PRunBase.cpp index 4bb89453..b9c15cf8 100644 --- a/src/classes/PRunBase.cpp +++ b/src/classes/PRunBase.cpp @@ -210,9 +210,9 @@ void PRunBase::DeadTimeCorrection(std::vector &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 t_dt = runData->GetDeadTimeParam(); + Double_t t_dt = runData->GetDeadTimeParam(histoNo[i]); for (UInt_t j=0; jOffset 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 + //------------------------------------------------------------- /** *

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 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 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 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. }; diff --git a/src/include/PRunDataHandler.h b/src/include/PRunDataHandler.h index fa2fa2c8..32331c32 100644 --- a/src/include/PRunDataHandler.h +++ b/src/include/PRunDataHandler.h @@ -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(count_ds.GetAttribute("t0_bin")); @@ -893,7 +895,7 @@ Bool_t PRunDataHandler::ReadNexusFileIdf2(T& nxs_file) for (int i=0; i