/***************************************************************************

  PRunSingleHistoRRF.cpp

  Author: Andreas Suter
  e-mail: andreas.suter@psi.ch

***************************************************************************/

/***************************************************************************
 *   Copyright (C) 2007-2025 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.             *
 ***************************************************************************/

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef HAVE_GOMP
#include <omp.h>
#endif

#include <cmath>
#include <iostream>
#include <fstream>
#include <memory>

#include <TString.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TH1F.h>

#include "PMusr.h"
#include "PFourier.h"
#include "PRunSingleHistoRRF.h"

//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
 * \brief Default constructor for RRF single histogram fitting class.
 *
 * Initializes all member variables to safe default values:
 * - fNoOfFitBins = 0 (no bins to fit)
 * - fBackground = 0 (will be estimated or set from MSR file)
 * - fBkgErr = 1.0 (default error estimate)
 * - fRRFPacking = -1 (invalid until set from GLOBAL block)
 * - fTheoAsData = false (high-resolution theory grid)
 * - fGoodBins[0,1] = -1 (calculated from data range)
 * - fN0EstimateEndTime = 1.0 μs (default N₀ estimation window)
 *
 * \warning This constructor creates an invalid object that cannot be used
 *          until properly initialized with MSR file data. Use the full
 *          constructor for normal operation.
 *
 * \see PRunSingleHistoRRF(PMsrHandler*, PRunDataHandler*, UInt_t, EPMusrHandleTag, Bool_t)
 */
PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase()
{
  fNoOfFitBins  = 0;
  fBackground = 0.0;
  fBkgErr = 1.0;
  fRRFPacking = -1;
  fTheoAsData = false;

  // the 2 following variables are need in case fit range is given in bins, and since
  // the fit range can be changed in the command block, these variables need to be accessible
  fGoodBins[0] = -1;
  fGoodBins[1] = -1;
  
  fN0EstimateEndTime = 1.0; // end time in (us) over which N0 is estimated.
}

//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
 * \brief Main constructor for RRF single histogram fitting and viewing.
 *
 * Constructs a fully initialized RRF single histogram run object by:
 * -# Validating GLOBAL block presence (mandatory for RRF analysis)
 * -# Validating RRF frequency specification (rrf_freq in GLOBAL block)
 * -# Validating RRF packing specification (rrf_packing in GLOBAL block)
 * -# Calling PrepareData() to load histogram and perform RRF transformation
 *
 * <b>GLOBAL Block Requirements:</b>
 * The RRF fit type requires the following entries in the GLOBAL block:
 * - \c rrf_freq: Rotation frequency with unit (e.g., "13.554 T", "183.7 MHz")
 * - \c rrf_packing: Number of bins to average (integer)
 * - \c rrf_phase: (optional) Initial phase in degrees
 *
 * <b>Error Handling:</b>
 * If any validation fails, the constructor:
 * - Outputs detailed error message to stderr
 * - Sets fValid = false
 * - Returns immediately (PrepareData() is not called)
 *
 * \param msrInfo Pointer to MSR file handler (NOT owned, must outlive this object)
 * \param rawData Pointer to raw run data handler (NOT owned, must outlive this object)
 * \param runNo Zero-based index of the RUN block in the MSR file
 * \param tag Operation mode: kFit (fitting) or kView (viewing/plotting)
 * \param theoAsData If true, theory calculated only at data points (for viewing);
 *                   if false, theory uses finer time grid (8× data resolution)
 *
 * \warning GLOBAL block with RRF parameters is MANDATORY for this fit type.
 *          Always check IsValid() after construction.
 *
 * \note After construction, check IsValid() to ensure initialization succeeded.
 *       Common failure modes:
 *       - Missing GLOBAL block
 *       - Missing rrf_freq specification
 *       - Missing rrf_packing specification
 *       - Data file not found or histogram missing
 *
 * \see PrepareData(), PrepareFitData(), PrepareViewData()
 */
PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) :
  PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData)
{
  fNoOfFitBins  = 0;

  PMsrGlobalBlock *global = msrInfo->GetMsrGlobal();

  if (!global->IsPresent()) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no GLOBAL-block present!";
    std::cerr << std::endl << ">> For Single Histo RRF the GLOBAL-block is mandatory! Please fix this first.";
    std::cerr << std::endl;
    fValid = false;
    return;
  }

  if (!global->GetRRFUnit().CompareTo("??")) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Frequency found!";
    std::cerr << std::endl;
    fValid = false;
    return;
  }

  fRRFPacking = global->GetRRFPacking();
  if (fRRFPacking == -1) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Packing found!";
    std::cerr << std::endl;
    fValid = false;
    return;
  }

  // the 2 following variables are need in case fit range is given in bins, and since
  // the fit range can be changed in the command block, these variables need to be accessible
  fGoodBins[0] = -1;
  fGoodBins[1] = -1;

  fN0EstimateEndTime = 1.0; // end time in (us) over which N0 is estimated.
 
  if (!PrepareData()) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: Couldn't prepare data for fitting!";
    std::cerr << std::endl << ">> This is very bad :-(, will quit ...";
    std::cerr << std::endl;
    fValid = false;
  }
}

//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
 * \brief Destructor for RRF single histogram fitting class.
 *
 * Cleans up dynamically allocated memory:
 * - Clears the forward histogram data vector (fForward)
 * - Other vectors (fM, fMerr, fW, fAerr) are local to PrepareFitData
 *   and cleared automatically
 *
 * Base class destructor (PRunBase) handles cleanup of:
 * - Theory objects
 * - Function value arrays
 * - Other shared resources
 */
PRunSingleHistoRRF::~PRunSingleHistoRRF()
{
  fForward.clear();
}

//--------------------------------------------------------------------------
// CalcChiSquare (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates χ² between RRF-transformed data and theory (least-squares fit metric).
 *
 * Computes the standard chi-square goodness-of-fit statistic for RRF asymmetry:
 * \f[
 * \chi^2 = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{\sigma_i^2}
 * \f]
 *
 * Unlike standard single histogram fitting, no explicit N₀ or exponential decay
 * factors appear since the RRF transformation already produces dimensionless
 * asymmetry with properly propagated errors.
 *
 * <b>Algorithm:</b>
 * -# Evaluate all user-defined functions from FUNCTIONS block
 * -# Pre-evaluate theory at t=1.0 to initialize any stateful functions
 *    (e.g., LF relaxation, user functions with internal state)
 * -# Loop over fit range bins [fStartTimeBin, fEndTimeBin)
 * -# For each bin: calculate time, evaluate theory, accumulate χ²
 *
 * <b>OpenMP Parallelization:</b>
 * When compiled with OpenMP (HAVE_GOMP defined):
 * - Dynamic scheduling with chunk size = max(10, N_bins / N_processors)
 * - Private variables per thread: i, time, diff
 * - Reduction performed on chisq accumulator
 * - Thread-safe due to pre-evaluation of theory at t=1.0
 *
 * \param par Parameter vector from MINUIT minimizer, containing current
 *            estimates of all fit parameters
 *
 * \return χ² value (sum over all bins in fit range). Minimize this value
 *         during fitting to find optimal parameters.
 *
 * \note The theory function is evaluated in the RRF frame. The THEORY block
 *       should describe the low-frequency RRF signal, not the laboratory frame
 *       precession.
 *
 * \see CalcChiSquareExpected(), CalcMaxLikelihood()
 */
Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
{
  Double_t chisq     = 0.0;
  Double_t diff      = 0.0;

  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    UInt_t funcNo = fMsrInfo->GetFuncNo(i);
    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
  }

  // calculate chi square
  Double_t time(1.0);
  Int_t i;

  // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
  // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
  // for a given set of parameters---which should be done outside of the parallelized loop.
  // For all other functions it means a tiny and acceptable overhead.
  time = fTheory->Func(time, par, fFuncValues);

  #ifdef HAVE_GOMP
  Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
  if (chunk < 10)
    chunk = 10;
  #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
  #endif
  for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
    time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
    diff = fData.GetValue()->at(i) - fTheory->Func(time, par, fFuncValues);
    chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
  }

  return chisq;
}

//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates expected χ² using theory variance instead of data variance.
 *
 * Computes the expected chi-square where the error estimate in the denominator
 * comes from the theory prediction rather than the data:
 * \f[
 * \chi^2_{\rm exp} = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{A_{\rm RRF}^{\rm theory}(t_i)}
 * \f]
 *
 * This metric is useful for:
 * - Diagnostic purposes to assess fit quality
 * - Detecting systematic deviations from the model
 * - Comparing with standard χ² to identify error estimation issues
 *
 * <b>Algorithm:</b>
 * Same as CalcChiSquare() but uses theory value as variance estimate instead
 * of measured error bars. OpenMP parallelization is applied when available.
 *
 * \param par Parameter vector from MINUIT minimizer
 *
 * \return Expected χ² value. For a good fit, this should be approximately
 *         equal to the number of degrees of freedom (N_bins - N_params).
 *
 * \warning Theory values must be positive for valid variance estimate.
 *          Negative theory values can lead to incorrect χ² calculation.
 *
 * \see CalcChiSquare()
 */
Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
  Double_t chisq = 0.0;
  Double_t diff  = 0.0;
  Double_t theo  = 0.0;

  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    UInt_t funcNo = fMsrInfo->GetFuncNo(i);
    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
  }

  // calculate chi square
  Double_t time(1.0);
  Int_t i;

  // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
  // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
  // for a given set of parameters---which should be done outside of the parallelized loop.
  // For all other functions it means a tiny and acceptable overhead.
  time = fTheory->Func(time, par, fFuncValues);

  #ifdef HAVE_GOMP
  Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
  if (chunk < 10)
    chunk = 10;
  #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
  #endif
  for (i=fStartTimeBin; i < fEndTimeBin; ++i) {
    time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
    theo = fTheory->Func(time, par, fFuncValues);
    diff = fData.GetValue()->at(i) - theo;
    chisq += diff*diff / theo;
  }

  return chisq;
}

//--------------------------------------------------------------------------
// CalcMaxLikelihood (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates maximum likelihood for RRF data (NOT YET IMPLEMENTED).
 *
 * Maximum likelihood estimation for RRF single histogram data is more complex
 * than for raw histograms due to the non-linear transformation from
 * Poisson-distributed counts to RRF asymmetry.
 *
 * <b>Theoretical Background:</b>
 * For raw histogram data, the likelihood is:
 * \f[
 * \mathcal{L} = \prod_i \frac{\mu_i^{n_i} e^{-\mu_i}}{n_i!}
 * \f]
 * where \f$\mu_i\f$ is the expected count and \f$n_i\f$ is the observed count.
 *
 * For RRF-transformed data, the error propagation through the transformation
 * must be properly accounted for in the likelihood function.
 *
 * <b>Current Implementation:</b>
 * Returns 0.0 (not implemented). Use χ² minimization (CalcChiSquare) instead.
 *
 * \param par Parameter vector from MINUIT minimizer (unused)
 *
 * \return 0.0 (not implemented)
 *
 * \todo Implement proper maximum likelihood for RRF data by:
 *       -# Deriving the likelihood function for transformed asymmetry
 *       -# Accounting for error propagation through RRF transformation
 *       -# Including correlations introduced by packing
 *
 * \see CalcChiSquare() for currently supported fit metric
 */
Double_t PRunSingleHistoRRF::CalcMaxLikelihood(const std::vector<Double_t>& par)
{
  // not yet implemented

  return 0.0;
}

//--------------------------------------------------------------------------
// CalcTheory (public)
//--------------------------------------------------------------------------
/**
 * \brief Evaluates theory function at all data points for viewing/plotting.
 *
 * Calculates the theoretical RRF asymmetry using the current MSR parameter
 * values and stores results in fData for display. This method is called
 * after fitting to generate the theory curve overlay.
 *
 * <b>Algorithm:</b>
 * -# Extract parameter values from MSR parameter list
 * -# Evaluate all user-defined functions from FUNCTIONS block
 * -# Loop over data points (size matches RRF-packed data)
 * -# Calculate time: t = dataTimeStart + i × dataTimeStep
 * -# Evaluate theory: P(t) = Func(t, par, funcValues)
 * -# Store results via fData.AppendTheoryValue()
 *
 * <b>Theory Function:</b>
 * The theory is evaluated directly in the RRF frame. The THEORY block should
 * specify the low-frequency RRF signal (after transformation), not the
 * laboratory-frame high-frequency precession.
 *
 * Example THEORY block for RRF analysis:
 * \code
 * THEORY
 *   asymmetry    1        (amplitude)
 *   simpleGss    2        (Gaussian relaxation)
 *   TFieldCos   map1 fun1 (frequency shift, phase)
 * \endcode
 *
 * where the frequency in TFieldCos is the difference frequency (ω - ω_RRF).
 *
 * \note Unlike CalcChiSquare(), this method does not return a value.
 *       Results are stored internally in fData.fTheory vector.
 *
 * \see CalcChiSquare(), PrepareViewData()
 */
void PRunSingleHistoRRF::CalcTheory()
{
  // feed the parameter vector
  std::vector<Double_t> par;
  PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
  for (UInt_t i=0; i<paramList->size(); i++)
    par.push_back((*paramList)[i].fValue);

  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
  }

  // calculate theory
  UInt_t size = fData.GetValue()->size();
  Double_t start = fData.GetDataTimeStart();
  Double_t resolution = fData.GetDataTimeStep();
  Double_t time;
  for (UInt_t i=0; i<size; i++) {
    time = start + static_cast<Double_t>(i)*resolution;
    fData.AppendTheoryValue(fTheory->Func(time, par, fFuncValues));
  }

  // clean up
  par.clear();
}

//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
 * \brief Returns the number of bins included in the current fit range.
 *
 * Triggers CalcNoOfFitBins() to ensure the bin count is current based on
 * the latest fit range settings, then returns the cached value.
 *
 * The number of fit bins is needed for:
 * - Calculating degrees of freedom: ν = N_bins - N_params
 * - Reduced χ²: χ²_red = χ² / ν
 * - Statistical diagnostics and fit quality assessment
 *
 * \return Number of RRF-packed bins within [fFitStartTime, fFitEndTime].
 *         This accounts for RRF packing: fewer bins than raw data.
 *
 * \note The fit range may be modified during fitting by COMMANDS block
 *       instructions. Always call this method to get the current count.
 *
 * \see CalcNoOfFitBins(), SetFitRangeBin()
 */
UInt_t PRunSingleHistoRRF::GetNoOfFitBins()
{
  CalcNoOfFitBins();

  return fNoOfFitBins;
}

//--------------------------------------------------------------------------
// SetFitRangeBin (public)
//--------------------------------------------------------------------------
/**
 * \brief Sets fit range using bin-offset syntax from COMMANDS block.
 *
 * Allows dynamic modification of the fit range during fitting, as specified
 * in the COMMANDS block. This is used for systematic studies where the fit
 * range needs to be varied.
 *
 * <b>Syntax:</b>
 * \code
 * FIT_RANGE fgb[+n0] lgb[-n1]                          # single range for all runs
 * FIT_RANGE fgb+n00 lgb-n01 fgb+n10 lgb-n11 ...        # per-run ranges
 * \endcode
 *
 * where:
 * - \c fgb = first good bin (symbolic, replaced by actual value)
 * - \c lgb = last good bin (symbolic, replaced by actual value)
 * - \c +n0 = positive offset added to fgb
 * - \c -n1 = positive offset subtracted from lgb
 *
 * <b>Conversion to Time:</b>
 * \f[
 * t_{\rm start} = ({\rm fgb} + n_0 - t_0) \times \Delta t
 * \f]
 * \f[
 * t_{\rm end} = ({\rm lgb} - n_1 - t_0) \times \Delta t
 * \f]
 *
 * where Δt is the raw time resolution (before RRF packing).
 *
 * <b>Multiple Run Handling:</b>
 * - Single pair: Applied to all runs
 * - Multiple pairs: Must match number of RUN blocks; each run uses its own range
 * - Run selection: Uses (2 × runNo + 1) to index into token array
 *
 * <b>Example:</b>
 * \code
 * COMMANDS
 *   FIT_RANGE fgb+100 lgb-500    # start 100 bins after fgb, end 500 bins before lgb
 * \endcode
 *
 * \param fitRange String containing FIT_RANGE specification from COMMANDS block
 *
 * \warning Invalid syntax produces error message but does not throw exception.
 *          The previous fit range values are retained on error.
 *
 * \see GetProperFitRange(), CalcNoOfFitBins()
 */
void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange)
{
  TObjArray *tok = nullptr;
  TObjString *ostr = nullptr;
  TString str;
  Ssiz_t idx = -1;
  Int_t offset = 0;

  tok = fitRange.Tokenize(" \t");

  if (tok->GetEntries() == 3) { // structure FIT_RANGE fgb+n0 lgb-n1
    // handle fgb+n0 entry
    ostr = dynamic_cast<TObjString*>(tok->At(1));
    str = ostr->GetString();
    // check if there is an offset present
    idx = str.First("+");
    if (idx != -1) { // offset present
      str.Remove(0, idx+1);
      if (str.IsFloat()) // if str is a valid number, convert is to an integer
        offset = str.Atoi();
    }
    fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution;

    // handle lgb-n1 entry
    ostr = dynamic_cast<TObjString*>(tok->At(2));
    str = ostr->GetString();
    // check if there is an offset present
    idx = str.First("-");
    if (idx != -1) { // offset present
      str.Remove(0, idx+1);
      if (str.IsFloat()) // if str is a valid number, convert is to an integer
        offset = str.Atoi();
    }
    fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution;
  } else if ((tok->GetEntries() > 3) && (tok->GetEntries() % 2 == 1)) { // structure FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]]
    Int_t pos = 2*(fRunNo+1)-1;

    if (pos + 1 >= tok->GetEntries()) {
      std::cerr << std::endl << ">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
      std::cerr << std::endl << ">> will ignore it. Sorry ..." << std::endl;
    } else {
      // handle fgb+n0 entry
      ostr = dynamic_cast<TObjString*>(tok->At(pos));
      str = ostr->GetString();
      // check if there is an offset present
      idx = str.First("+");
      if (idx != -1) { // offset present
        str.Remove(0, idx+1);
        if (str.IsFloat()) // if str is a valid number, convert is to an integer
          offset = str.Atoi();
      }
      fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution;

      // handle lgb-n1 entry
      ostr = dynamic_cast<TObjString*>(tok->At(pos+1));
      str = ostr->GetString();
      // check if there is an offset present
      idx = str.First("-");
      if (idx != -1) { // offset present
        str.Remove(0, idx+1);
        if (str.IsFloat()) // if str is a valid number, convert is to an integer
          offset = str.Atoi();
      }
      fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution;
    }
  } else { // error
    std::cerr << std::endl << ">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
    std::cerr << std::endl << ">> will ignore it. Sorry ..." << std::endl;
  }

  // clean up
  if (tok) {
    delete tok;
  }
}

//--------------------------------------------------------------------------
// CalcNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates start/end bin indices from fit time range.
 *
 * Converts the fit range from time (μs) to RRF-packed bin indices.
 * This method is called whenever the fit range may have changed.
 *
 * <b>Conversion Formulas:</b>
 * \f[
 * i_{\rm start} = \lceil \frac{t_{\rm start} - t_{\rm data,0}}{\Delta t_{\rm data}} \rceil
 * \f]
 * \f[
 * i_{\rm end} = \lfloor \frac{t_{\rm end} - t_{\rm data,0}}{\Delta t_{\rm data}} \rfloor + 1
 * \f]
 *
 * where:
 * - \f$t_{\rm data,0}\f$ = fData.GetDataTimeStart() (first RRF-packed bin center)
 * - \f$\Delta t_{\rm data}\f$ = fData.GetDataTimeStep() (RRF-packed bin width)
 *
 * <b>Bounds Checking:</b>
 * - fStartTimeBin clamped to [0, data size)
 * - fEndTimeBin clamped to [0, data size]
 * - fNoOfFitBins = 0 if fEndTimeBin <= fStartTimeBin
 *
 * <b>Side Effects:</b>
 * Updates member variables:
 * - fStartTimeBin: First bin index in fit range (inclusive)
 * - fEndTimeBin: Last bin index in fit range (exclusive)
 * - fNoOfFitBins: Number of bins = fEndTimeBin - fStartTimeBin
 *
 * \note Time step includes RRF packing: dataTimeStep = rawTimeRes × fRRFPacking
 *
 * \see GetNoOfFitBins(), SetFitRangeBin()
 */
void PRunSingleHistoRRF::CalcNoOfFitBins()
{
  // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly
  fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
  if (fStartTimeBin < 0)
    fStartTimeBin = 0;
  fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
  if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
    fEndTimeBin = fData.GetValue()->size();

  if (fEndTimeBin > fStartTimeBin)
    fNoOfFitBins = fEndTimeBin - fStartTimeBin;
  else
    fNoOfFitBins = 0;
}

//--------------------------------------------------------------------------
// PrepareData (protected)
//--------------------------------------------------------------------------
/**
 * \brief Main data preparation orchestrator for RRF single histogram analysis.
 *
 * Coordinates the loading and preprocessing of histogram data before RRF
 * transformation. This method handles all operations common to both fitting
 * and viewing modes.
 *
 * <b>Processing Steps:</b>
 * -# Validate object state (return false if already marked invalid)
 * -# Get GLOBAL block reference for settings
 * -# Load raw run data from data handler using run name from MSR file
 * -# Extract metadata from data file:
 *    - Magnetic field (fMetaData.fField)
 *    - Beam energy (fMetaData.fEnergy)
 *    - Sample temperature(s) (fMetaData.fTemp)
 * -# Collect histogram numbers from RUN block forward specification
 * -# Validate all specified histograms exist in data file
 * -# Determine time resolution: ns → μs conversion
 * -# Determine t0 values via GetProperT0()
 * -# Initialize forward histogram from first group
 * -# Handle addrun (co-add multiple runs with t0 alignment)
 * -# Handle grouping (combine multiple detectors with t0 alignment)
 * -# Determine data range via GetProperDataRange()
 * -# Determine fit range via GetProperFitRange()
 * -# Dispatch to PrepareFitData() or PrepareViewData() based on tag
 *
 * <b>Run Addition (addrun):</b>
 * When multiple runs are specified in the RUN block, histograms are co-added
 * with t0 alignment:
 * \code
 * fForward[i] += addRunData[i + addT0 - mainT0]
 * \endcode
 *
 * <b>Detector Grouping:</b>
 * When multiple forward histograms are specified, they are summed with
 * t0 alignment to form a single combined histogram.
 *
 * \return true if data preparation succeeds, false on any error:
 *         - Run data not found
 *         - Histogram not present in data file
 *         - Invalid t0 determination
 *         - Data range validation failure
 *         - Fit/view data preparation failure
 *
 * \see PrepareFitData(), PrepareViewData(), GetProperT0(), GetProperDataRange()
 */
Bool_t PRunSingleHistoRRF::PrepareData()
{
  Bool_t success = true;

  if (!fValid)
    return false;

  // keep the Global block info
  PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal();

  // get the proper run
  PRawRunData* runData = fRawData->GetRunData(*fRunInfo->GetRunName());
  if (!runData) { // couldn't get run
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get run " << fRunInfo->GetRunName()->Data() << "!";
    std::cerr << std::endl;
    return false;
  }

  // keep the field from the meta-data from the data-file
  fMetaData.fField = runData->GetField();

  // keep the energy from the meta-data from the data-file
  fMetaData.fEnergy = runData->GetEnergy();

  // keep the temperature(s) from the meta-data from the data-file
  for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
    fMetaData.fTemp.push_back(runData->GetTemperature(i));

  // collect histogram numbers
  PUIntVector histoNo; // histoNo = msr-file forward + redGreen_offset - 1
  for (UInt_t i=0; i<fRunInfo->GetForwardHistoNoSize(); i++) {
    histoNo.push_back(fRunInfo->GetForwardHistoNo(i));

    if (!runData->IsPresent(histoNo[i])) {
      std::cerr << std::endl << ">> PRunSingleHistoRRF::PrepareData(): **PANIC ERROR**:";
      std::cerr << std::endl << ">> histoNo found = " << histoNo[i] << ", which is NOT present in the data file!?!?";
      std::cerr << std::endl << ">> Will quit :-(";
      std::cerr << std::endl;
      histoNo.clear();
      return false;
    }
  }

  // keep the time resolution in (us)
  fTimeResolution = runData->GetTimeResolution()/1.0e3;
  std::cout.precision(10);
  std::cout << std::endl << ">> PRunSingleHisto::PrepareData(): time resolution=" << std::fixed << runData->GetTimeResolution() << "(ns)" << std::endl;

  // get all the proper t0's and addt0's for the current RUN block
  if (!GetProperT0(runData, globalBlock, histoNo)) {
    return false;
  }

  // keep the histo of each group at this point (addruns handled below)
  std::vector<PDoubleVector> forward;
  forward.resize(histoNo.size());   // resize to number of groups
  for (UInt_t i=0; i<histoNo.size(); i++) {
    forward[i].resize(runData->GetDataBin(histoNo[i])->size());
    forward[i] = *runData->GetDataBin(histoNo[i]);
  }

  // check if there are runs to be added to the current one
  if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present
    PRawRunData *addRunData;
    for (UInt_t i=1; i<fRunInfo->GetRunNameSize(); i++) {

      // get run to be added to the main one
      addRunData = fRawData->GetRunData(*fRunInfo->GetRunName(i));
      if (addRunData == nullptr) { // couldn't get run
        std::cerr << std::endl << ">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!";
        std::cerr << std::endl;
        return false;
      }

      // add forward run
      UInt_t addRunSize;
      for (UInt_t k=0; k<histoNo.size(); k++) { // fill each group
        addRunSize = addRunData->GetDataBin(histoNo[k])->size();
        for (UInt_t j=0; j<addRunData->GetDataBin(histoNo[k])->size(); j++) { // loop over the bin indices
          // make sure that the index stays in the proper range
          if ((static_cast<Int_t>(j)+static_cast<Int_t>(fAddT0s[i-1][k])-static_cast<Int_t>(fT0s[k]) >= 0) &&
              (j+static_cast<Int_t>(fAddT0s[i-1][k])-static_cast<Int_t>(fT0s[k]) < addRunSize)) {
            forward[k][j] += addRunData->GetDataBin(histoNo[k])->at(j+static_cast<Int_t>(fAddT0s[i-1][k])-static_cast<Int_t>(fT0s[k]));
          }
        }
      }
    }
  }

  // set forward histo data of the first group
  fForward.resize(forward[0].size());
  for (UInt_t i=0; i<fForward.size(); i++) {
    fForward[i] = forward[0][i];
  }

  // group histograms, add all the remaining forward histograms of the group
  for (UInt_t i=1; i<histoNo.size(); i++) { // loop over the groupings
    for (UInt_t j=0; j<runData->GetDataBin(histoNo[i])->size(); j++) { // loop over the bin indices
      // make sure that the index stays within proper range
      if ((static_cast<Int_t>(j)+static_cast<Int_t>(fT0s[i])-static_cast<Int_t>(fT0s[0]) >= 0) &&
          (j+static_cast<Int_t>(fT0s[i])-static_cast<Int_t>(fT0s[0]) < runData->GetDataBin(histoNo[i])->size())) {
        fForward[j] += forward[i][j+static_cast<Int_t>(fT0s[i])-static_cast<Int_t>(fT0s[0])];
      }
    }
  }

  // get the data range (fgb/lgb) for the current RUN block
  if (!GetProperDataRange()) {
    return false;
  }

  // get the fit range for the current RUN block
  GetProperFitRange(globalBlock);

  // do the more fit/view specific stuff
  if (fHandleTag == kFit)
    success = PrepareFitData(runData, histoNo[0]);
  else if (fHandleTag == kView)
    success = PrepareViewData(runData, histoNo[0]);
  else
    success = false;

  // cleanup
  histoNo.clear();

  return success;
}

//--------------------------------------------------------------------------
// PrepareFitData (protected)
//--------------------------------------------------------------------------
/**
 * \brief Performs full RRF transformation on histogram data for fitting.
 *
 * Takes the pre-processed histogram (grouping and addrun already applied) and
 * transforms it into RRF asymmetry suitable for fitting. This is the core
 * method implementing the RRF analysis technique.
 *
 * <b>Processing Steps:</b>
 *
 * <b>1. Frequency Analysis:</b>
 * - Calls GetMainFrequency() on raw data to find dominant precession frequency
 * - Used to determine optimal N₀ estimation window (integer oscillation cycles)
 * - Prints "optimal packing" suggestion for user reference
 *
 * <b>2. Background Handling:</b>
 * - If fixed background specified: use directly from RUN block
 * - If background range given: estimate via EstimateBkg()
 * - If neither: estimate range as [0.1×t0, 0.6×t0] with warning
 * - Subtract background: N(t) → N(t) - B
 *
 * <b>3. Lifetime Correction:</b>
 * \f[
 * M(t) = [N(t) - B] \cdot e^{+t/\tau_\mu}
 * \f]
 * - Removes exponential decay from signal
 * - Error: \f$\sigma_M = e^{+t/\tau_\mu} \sqrt{N(t) + \sigma_B^2}\f$
 * - Weights: \f$w = 1/\sigma_M^2\f$
 *
 * <b>4. N₀ Estimation:</b>
 * - Call EstimateN0() with frequency information
 * - Uses average of M(t) over integer number of oscillation cycles
 * - Returns N₀ and its error σ_N₀
 *
 * <b>5. Asymmetry Extraction:</b>
 * \f[
 * A(t) = \frac{M(t)}{N_0} - 1
 * \f]
 * Error propagation:
 * \f[
 * \sigma_A(t) = \frac{e^{+t/\tau_\mu}}{N_0} \sqrt{N(t) + \left(\frac{N(t)-B}{N_0}\right)^2 \sigma_{N_0}^2}
 * \f]
 *
 * <b>6. RRF Rotation:</b>
 * \f[
 * A_{\rm RRF}(t) = 2 \cdot A(t) \cdot \cos(\omega_{\rm RRF} t + \phi_{\rm RRF})
 * \f]
 * - ω_RRF from GLOBAL block (rrf_freq converted to angular frequency)
 * - φ_RRF from GLOBAL block (rrf_phase in radians)
 * - Factor 2 compensates for discarded counter-rotating component
 *
 * <b>7. RRF Packing:</b>
 * - Average over fRRFPacking consecutive bins
 * - Data: \f$A_{\rm packed} = \frac{1}{n}\sum_{i=1}^{n} A_{\rm RRF}(t_i)\f$
 * - Error: \f$\sigma_{\rm packed} = \frac{\sqrt{2}}{n}\sqrt{\sum_{i=1}^{n} \sigma_A^2(t_i)}\f$
 * - √2 factor accounts for doubling in RRF rotation
 *
 * <b>8. Time Grid Setup:</b>
 * - Data time start: accounts for packing offset (center of first packed bin)
 * - Data time step: rawTimeResolution × fRRFPacking
 *
 * \param runData Raw run data handler (for background estimation)
 * \param histoNo Forward histogram number (0-based, for background error messages)
 *
 * \return true on success, false on error:
 *         - Background estimation failure
 *         - Invalid bin ranges
 *
 * \see PrepareViewData(), GetMainFrequency(), EstimateN0(), EstimateBkg()
 */
Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t histoNo)
{
  // keep the raw data for the RRF asymmetry error estimate for later
  PDoubleVector rawNt;
  for (UInt_t i=0; i<fForward.size(); i++) {
    rawNt.push_back(fForward[i]);  // N(t) without any corrections
  }
  Double_t freqMax = GetMainFrequency(rawNt);
  std::cout << "info> freqMax=" << freqMax << " (MHz)" << std::endl;

  // "optimal packing"
  Double_t optNoPoints = 8;
  if (freqMax < 271.0) // < 271 MHz, i.e ~ 2T
    optNoPoints = 5;
  std::cout << "info> optimal packing: " << static_cast<Int_t>(1.0 / (fTimeResolution*(freqMax - fMsrInfo->GetMsrGlobal()->GetRRFFreq("MHz"))) / optNoPoints);

  // initially fForward is the "raw data set" (i.e. grouped histo and raw runs already added) to be fitted. This means fForward = N(t) at this point.

  // 1) check how the background shall be handled
  // subtract background from histogramms ------------------------------------------
  if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given
    if (fRunInfo->GetBkgRange(0) >= 0) {
      if (!EstimateBkg(histoNo))
        return false;
    } else { // no background given to do the job, try estimate
      fRunInfo->SetBkgRange(static_cast<Int_t>(fT0s[0]*0.1), 0);
      fRunInfo->SetBkgRange(static_cast<Int_t>(fT0s[0]*0.6), 1);
      std::cerr << std::endl << ">> PRunSingleHistoRRF::PrepareFitData(): **WARNING** Neither fix background nor background bins are given!";
      std::cerr << std::endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1);
      std::cerr << std::endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ...";
      std::cerr << std::endl;
      if (!EstimateBkg(histoNo))
        return false;
    }
    // subtract background from fForward
    for (UInt_t i=0; i<fForward.size(); i++)
      fForward[i] -= fBackground;
  } else { // fixed background given
    for (UInt_t i=0; i<fForward.size(); i++) {
      fForward[i] -= fRunInfo->GetBkgFix(0);
    }
    fBackground = fRunInfo->GetBkgFix(0);
  }
  // here fForward = N(t) - Nbkg

  Int_t t0 = static_cast<Int_t>(fT0s[0]);

  // 2) N(t) - Nbkg -> exp(+t/tau) [N(t)-Nbkg]
  Double_t startTime = fTimeResolution * (static_cast<Double_t>(fGoodBins[0]) - static_cast<Double_t>(t0));

  Double_t time_tau=0.0;
  Double_t exp_t_tau=0.0;
  for (Int_t i=fGoodBins[0]; i<fGoodBins[1]; i++) {
    time_tau = (startTime + fTimeResolution * (i - fGoodBins[0])) / PMUON_LIFETIME;
    exp_t_tau = exp(time_tau);
    fForward[i] *= exp_t_tau;
    fM.push_back(fForward[i]);   // i.e. M(t) = [N(t)-Nbkg] exp(+t/tau); needed to estimate N0 later on
    fMerr.push_back(exp_t_tau*sqrt(rawNt[i]+fBkgErr*fBkgErr));
  }

  // calculate weights
  for (UInt_t i=0; i<fMerr.size(); i++) {
    if (fMerr[i] > 0.0)
      fW.push_back(1.0/(fMerr[i]*fMerr[i]));
    else
      fW.push_back(1.0);
  }
  // now fForward = exp(+t/tau) [N(t)-Nbkg] = M(t)

  // 3) estimate N0
  Double_t errN0 = 0.0;
  Double_t n0 = EstimateN0(errN0, freqMax);

  // 4a) A(t) = exp(+t/tau) [N(t)-Nbkg] / N0 - 1.0
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    fForward[i] = fForward[i] / n0 - 1.0;
  }

  // 4b) error estimate of A(t): errA(t) = exp(+t/tau)/N0 sqrt( N(t) + ([N(t)-N_bkg]/N0)^2 errN0^2 )
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    time_tau = (startTime + fTimeResolution * (i - fGoodBins[0])) / PMUON_LIFETIME;
    exp_t_tau = exp(time_tau);
    fAerr.push_back(exp_t_tau/n0*sqrt(rawNt[i]+pow(((rawNt[i]-fBackground)/n0)*errN0,2.0)));
  }

  // 5) rotate A(t): A(t) -> 2* A(t) * cos(wRRF t + phiRRF), the factor 2.0 is needed since the high frequency part is suppressed.
  PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal();
  Double_t wRRF = globalBlock->GetRRFFreq("Mc");
  Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0;
  Double_t time = 0.0;
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    time = startTime + fTimeResolution * (static_cast<Double_t>(i) - static_cast<Double_t>(fGoodBins[0]));
    fForward[i] *= 2.0*cos(wRRF * time + phaseRRF);
  }

  // 6) RRF packing
  Double_t dval=0.0;
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    if (fRRFPacking == 1) {
      fData.AppendValue(fForward[i]);
    } else { // RRF packing > 1
      if (((i-fGoodBins[0]) % fRRFPacking == 0) && (i != fGoodBins[0])) { // fill data
        dval /= fRRFPacking;
        fData.AppendValue(dval);
        // reset dval
        dval = 0.0;
      }
      dval += fForward[i];
    }
  }

  // 7) estimate packed RRF errors (see log-book p.204)
  //    the error estimate of the unpacked RRF asymmetry is: errA_RRF(t) \simeq exp(t/tau)/N0 sqrt( [N(t) + ((N(t)-N_bkg)/N0)^2 errN0^2] )
  dval = 0.0;
  // the packed RRF asymmetry error
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    if (((i-fGoodBins[0]) % fRRFPacking == 0) && (i != fGoodBins[0])) { // fill data
      fData.AppendErrorValue(sqrt(2.0*dval)/fRRFPacking); // the factor 2.0 is needed since the high frequency part is suppressed.
      dval = 0.0;
    }
    dval += fAerr[i-fGoodBins[0]]*fAerr[i-fGoodBins[0]];
  }

  // set start time and time step
  fData.SetDataTimeStart(fTimeResolution*(static_cast<Double_t>(fGoodBins[0])-static_cast<Double_t>(t0)+static_cast<Double_t>(fRRFPacking-1)/2.0));
  fData.SetDataTimeStep(fTimeResolution*fRRFPacking);

  CalcNoOfFitBins();

  return true;
}

//--------------------------------------------------------------------------
// PrepareViewData (protected)
//--------------------------------------------------------------------------
/**
 * \brief Prepares RRF-transformed data for viewing/plotting.
 *
 * Takes the pre-processed histogram and prepares both data and theory curves
 * for display. This method is used when the operation mode is kView rather
 * than kFit.
 *
 * <b>Processing Steps:</b>
 *
 * <b>1. Data Preparation:</b>
 * - Calls PrepareFitData() to perform full RRF transformation
 * - Data is ready for display after this step
 *
 * <b>2. View Packing Check:</b>
 * - Checks if additional view packing is specified in PLOT block
 * - If viewPacking > fRRFPacking: additional packing would be applied (TODO)
 * - If viewPacking < fRRFPacking: warning issued and view packing ignored
 *   (cannot unpack already-packed RRF data)
 *
 * <b>3. Theory Grid Setup:</b>
 * Determines theory evaluation resolution:
 * - If fTheoAsData = true: theory calculated only at data points
 * - If fTheoAsData = false: theory calculated on 8× finer grid for smooth curves
 *
 * Time grid:
 * - Theory time start = data time start
 * - Theory time step = data time step / factor (where factor = 1 or 8)
 *
 * <b>4. Theory Evaluation:</b>
 * - Extract parameter values from MSR parameter list
 * - Evaluate all user-defined functions from FUNCTIONS block
 * - Loop over theory grid points
 * - Evaluate theory: P(t) = Func(t, par, funcValues)
 * - Apply sanity check: |P(t)| > 10 → set to 0 (dirty hack for edge cases)
 * - Store results via fData.AppendTheoryValue()
 *
 * <b>Theory Function:</b>
 * The theory is evaluated directly in the RRF frame. The THEORY block should
 * specify the low-frequency signal after RRF transformation, not the
 * laboratory-frame high-frequency precession.
 *
 * \param runData Raw run data handler (passed to PrepareFitData)
 * \param histoNo Forward histogram number (passed to PrepareFitData)
 *
 * \return true on success, false on error (typically from PrepareFitData)
 *
 * \note The 8× theory resolution provides smoother curves for visualization
 *       and better Fourier transforms without affecting fit performance.
 *
 * \see PrepareFitData(), CalcTheory()
 */
Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t histoNo)
{
  // --------------
  // prepare data
  // --------------

  // prepare RRF single histo
  PrepareFitData(runData, histoNo);

  // check for view packing
  Int_t viewPacking = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking;
  if (viewPacking > 0) {
    if (viewPacking < fRRFPacking) {
      std::cerr << ">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** Found View Packing (" << viewPacking << ") < RRF Packing (" << fRRFPacking << ").";
      std::cerr << ">> Will ignore View Packing." << std::endl;
    } else {
      // STILL MISSING
    }
  }

  // --------------
  // prepare theory
  // --------------

  // feed the parameter vector
  std::vector<Double_t> par;
  PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
  for (UInt_t i=0; i<paramList->size(); i++)
    par.push_back((*paramList)[i].fValue);

  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
  }

  // check if a finer binning for the theory is needed
  UInt_t size = fForward.size();
  Int_t factor = 8; // 8 times more points for the theory (if fTheoAsData == false)
  fData.SetTheoryTimeStart(fData.GetDataTimeStart());
  if (fTheoAsData) { // calculate theory only at the data points
    fData.SetTheoryTimeStep(fData.GetDataTimeStep());
  } else {
    // finer binning for the theory (8 times as many points = factor)
    size *= factor;
    fData.SetTheoryTimeStep(fData.GetDataTimeStep()/(Double_t)factor);
  }

  // calculate theory
  Double_t time = 0.0;
  Double_t theoryValue = 0.0;
  for (UInt_t i=0; i<size; i++) {
    time = fData.GetTheoryTimeStart() + static_cast<Double_t>(i)*fData.GetTheoryTimeStep();
    theoryValue = fTheory->Func(time, par, fFuncValues);
    if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!!
      theoryValue = 0.0;
    }
    fData.AppendTheoryValue(theoryValue);
  }

  return true;
}

//--------------------------------------------------------------------------
// GetProperT0 (private)
//--------------------------------------------------------------------------
/**
 * \brief Determines and validates t0 values for all histograms in the run.
 *
 * Searches for time-zero (t0) values through a priority cascade and performs
 * validation. t0 marks the muon arrival time and is critical for proper
 * timing alignment.
 *
 * <b>t0 Search Priority:</b>
 * -# RUN block: t0 values specified in the MSR file RUN block
 * -# GLOBAL block: fallback t0 values from GLOBAL block
 * -# Data file: t0 values stored in the raw data file header
 * -# Automatic estimation: estimated from histogram shape (with warning)
 *
 * <b>Algorithm:</b>
 * -# Initialize fT0s vector with size = number of forward histograms (grouping)
 * -# Set all t0 values to -1.0 (unset marker)
 * -# Fill from RUN block if present
 * -# Fill remaining -1.0 values from GLOBAL block
 * -# Fill remaining -1.0 values from data file
 * -# Fill remaining -1.0 values from automatic estimation (with warning)
 * -# Validate all t0 values are within histogram bounds
 *
 * <b>Addrun Handling:</b>
 * When multiple runs are co-added (addrun), each additional run needs its own
 * t0 value for proper time alignment:
 * -# Initialize fAddT0s[runIdx] for each addrun
 * -# Apply same priority cascade for each addrun
 * -# Validate addrun t0 values
 *
 * <b>Validation:</b>
 * - t0 must be ≥ 0
 * - t0 must be < histogram length
 * - Error message and return false on validation failure
 *
 * \param runData Pointer to raw run data for histogram access and data file t0
 * \param globalBlock Pointer to GLOBAL block for fallback t0 values
 * \param histoNo Vector of histogram indices (0-based) for forward grouping
 *
 * \return true if valid t0 found for all histograms, false on:
 *         - t0 out of histogram bounds
 *         - addrun data not accessible
 *         - addrun t0 validation failure
 *
 * \warning Automatic t0 estimation is unreliable for LEM and other specialized
 *          setups. Always specify t0 explicitly for best results.
 *
 * \see GetProperDataRange(), PrepareData()
 */
Bool_t PRunSingleHistoRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo)
{
  // feed all T0's
  // first init T0's, T0's are stored as (forward T0, backward T0, etc.)
  fT0s.clear();
  fT0s.resize(histoNo.size());
  for (UInt_t i=0; i<fT0s.size(); i++) {
    fT0s[i] = -1.0;
  }

  // fill in the T0's from the msr-file (if present)
  for (UInt_t i=0; i<fRunInfo->GetT0BinSize(); i++) {
    fT0s[i] = fRunInfo->GetT0Bin(i);
  }

  // fill in the T0's from the GLOBAL block section (if present)
  for (UInt_t i=0; i<globalBlock->GetT0BinSize(); i++) {
    if (fT0s[i] == -1.0) { // i.e. not given in the RUN block section
      fT0s[i] = globalBlock->GetT0Bin(i);
    }
  }

  // fill in the T0's from the data file, if not already present in the msr-file
  for (UInt_t i=0; i<histoNo.size(); i++) {
    if (fT0s[i] == -1.0) { // i.e. not present in the msr-file, try the data file
      if (runData->GetT0Bin(histoNo[i]) > 0.0) {
        fT0s[i] = runData->GetT0Bin(histoNo[i]);
        fRunInfo->SetT0Bin(fT0s[i], i); // keep value for the msr-file
      }
    }
  }

  // fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file
  for (UInt_t i=0; i<histoNo.size(); i++) {
    if (fT0s[i] == -1.0) { // i.e. not present in the msr-file and data file, use the estimated T0
      fT0s[i] = runData->GetT0BinEstimated(histoNo[i]);
      fRunInfo->SetT0Bin(fT0s[i], i); // keep value for the msr-file

      std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
      std::cerr << std::endl << ">> run: " << fRunInfo->GetRunName()->Data();
      std::cerr << std::endl << ">> will try the estimated one: forward t0 = " << runData->GetT0BinEstimated(histoNo[i]);
      std::cerr << std::endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
      std::cerr << std::endl;
    }
  }

  // check if t0 is within proper bounds
  for (UInt_t i=0; i<fRunInfo->GetForwardHistoNoSize(); i++) {
    if ((fT0s[i] < 0.0) || (fT0s[i] > static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size()))) {
      std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** t0 data bin (" << fT0s[i] << ") doesn't make any sense!";
      std::cerr << std::endl;
      return false;
    }
  }

  // check if there are runs to be added to the current one. If yes keep the needed t0's
  if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present
    PRawRunData *addRunData;
    fAddT0s.resize(fRunInfo->GetRunNameSize()-1); // resize to the number of addruns
    for (UInt_t i=1; i<fRunInfo->GetRunNameSize(); i++) {

      // get run to be added to the main one
      addRunData = fRawData->GetRunData(*fRunInfo->GetRunName(i));
      if (addRunData == nullptr) { // couldn't get run
        std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!";
        std::cerr << std::endl;
        return false;
      }

      // feed all T0's
      // first init T0's, T0's are stored as (forward T0, backward T0, etc.)
      fAddT0s[i-1].resize(histoNo.size());
      for (UInt_t j=0; j<fAddT0s[i-1].size(); j++) {
        fAddT0s[i-1][j] = -1.0;
      }

      // fill in the T0's from the msr-file (if present)
      for (UInt_t j=0; j<fRunInfo->GetT0BinSize(); j++) {
        fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i-1,j); // addRunIdx starts at 0
      }

      // fill in the T0's from the data file, if not already present in the msr-file
      for (UInt_t j=0; j<histoNo.size(); j++) {
        if (fAddT0s[i-1][j] == -1.0) // i.e. not present in the msr-file, try the data file
          if (addRunData->GetT0Bin(histoNo[j]) > 0.0) {
            fAddT0s[i-1][j] = addRunData->GetT0Bin(histoNo[j]);
            fRunInfo->SetAddT0Bin(fAddT0s[i-1][j], i-1, j); // keep value for the msr-file
          }
      }

      // fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file
      for (UInt_t j=0; j<histoNo.size(); j++) {
        if (fAddT0s[i-1][j] == -1.0) { // i.e. not present in the msr-file and data file, use the estimated T0
          fAddT0s[i-1][j] = addRunData->GetT0BinEstimated(histoNo[j]);
          fRunInfo->SetAddT0Bin(fAddT0s[i-1][j], i-1, j); // keep value for the msr-file

          std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
          std::cerr << std::endl << ">> run: " << fRunInfo->GetRunName(i)->Data();
          std::cerr << std::endl << ">> will try the estimated one: forward t0 = " << addRunData->GetT0BinEstimated(histoNo[j]);
          std::cerr << std::endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
          std::cerr << std::endl;
        }
      }

      // check if t0 is within proper bounds
      for (UInt_t j=0; j<fRunInfo->GetForwardHistoNoSize(); j++) {
        if ((fAddT0s[i-1][j] < 0) || (fAddT0s[i-1][j] > static_cast<Int_t>(addRunData->GetDataBin(histoNo[j])->size()))) {
          std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** addt0 data bin (" << fAddT0s[i-1][j] << ") doesn't make any sense!";
          std::cerr << std::endl;
          return false;
        }
      }
    }
  }

  return true;
}

//--------------------------------------------------------------------------
// GetProperDataRange (private)
//--------------------------------------------------------------------------
/**
 * \brief Determines valid data range (first/last good bins) for analysis.
 *
 * Establishes the boundaries of usable histogram data through a priority
 * cascade. The data range defines which bins contain valid signal (after t0,
 * before noise/overflow).
 *
 * <b>Data Range Search Priority:</b>
 * -# RUN block: data range specified in MSR file RUN block
 * -# GLOBAL block: fallback data range from GLOBAL block
 * -# Automatic estimation: estimate from t0 (with warning)
 *    - Start: t0 + 10 ns offset
 *    - End: histogram length
 *
 * <b>Validation:</b>
 * -# Check start < end (swap if reversed)
 * -# Check start ≥ 0 and start < histogram size
 * -# Check end ≥ 0 (adjust if > histogram size with warning)
 *
 * <b>Result:</b>
 * Sets fGoodBins[0] (first good bin) and fGoodBins[1] (last good bin).
 * These are used for:
 * - RRF transformation range
 * - Fit range specification in bins (fgb+n0 lgb-n1)
 * - Data validity checking
 *
 * \return true if valid data range determined, false on:
 *         - Start bin out of bounds
 *         - End bin negative
 *         - Other range validation failures
 *
 * \note The data range is typically wider than the fit range. Data range
 *       defines where valid data exists; fit range defines what is fitted.
 *
 * \see GetProperFitRange(), GetProperT0()
 */
Bool_t PRunSingleHistoRRF::GetProperDataRange()
{
  // get start/end data
  Int_t start;
  Int_t end;
  start = fRunInfo->GetDataRange(0);
  end   = fRunInfo->GetDataRange(1);

  // check if data range has been given in the RUN block, if not try to get it from the GLOBAL block
  if (start < 0) {
    start = fMsrInfo->GetMsrGlobal()->GetDataRange(0);
  }
  if (end < 0) {
    end = fMsrInfo->GetMsrGlobal()->GetDataRange(1);
  }

  // check if data range has been provided, and if not try to estimate them
  if (start < 0) {
    Int_t offset = static_cast<Int_t>(10.0e-3/fTimeResolution);
    start = static_cast<Int_t>(fT0s[0])+offset;
    fRunInfo->SetDataRange(start, 0);
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start << ".";
    std::cerr << std::endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE.";
    std::cerr << std::endl;
  }
  if (end < 0) {
    end = fForward.size();
    fRunInfo->SetDataRange(end, 1);
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range end = " << end << ".";
    std::cerr << std::endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE.";
    std::cerr << std::endl;
  }

  // check if start and end make any sense
  // 1st check if start and end are in proper order
  if (end < start) { // need to swap them
    Int_t keep = end;
    end = start;
    start = keep;
  }
  // 2nd check if start is within proper bounds
  if ((start < 0) || (start > static_cast<Int_t>(fForward.size()))) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** start data bin (" << start << ") doesn't make any sense!";
    std::cerr << std::endl;
    return false;
  }
  // 3rd check if end is within proper bounds
  if (end < 0) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** end data bin (" << end << ") doesn't make any sense!";
    std::cerr << std::endl;
    return false;
  }
  if (end > static_cast<Int_t>(fForward.size())) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** end data bin (" << end << ") > histo length (" << fForward.size() << ").";
    std::cerr << std::endl << ">>    Will set end = (histo length - 1). Consider to change it in the msr-file." << std::endl;
    std::cerr << std::endl;
    end = static_cast<Int_t>(fForward.size()-1);
  }

  // keep good bins for potential later use
  fGoodBins[0] = start;
  fGoodBins[1] = end;

  // make sure that fGoodBins are in proper range for fForward
  if (fGoodBins[0] < 0)
    fGoodBins[0]=0;
  if (fGoodBins[1] > fForward.size()) {
    std::cerr << std::endl << ">> PRunSingleHisto::GetProperDataRange **WARNING** needed to shift forward lgb,";
    std::cerr << std::endl << ">> from " << fGoodBins[1] << " to " << fForward.size()-1 << std::endl;
    fGoodBins[1]=fForward.size()-1;
  }

  return true;
}

//--------------------------------------------------------------------------
// GetProperFitRange (private)
//--------------------------------------------------------------------------
/**
 * \brief Determines fit time range from MSR file settings.
 *
 * Extracts the fit range (time window for parameter extraction) through a
 * priority cascade. Supports both time-based and bin-based specifications.
 *
 * <b>Fit Range Formats:</b>
 * - Time-based: \c fit \c 0.1 \c 8.0 (in microseconds relative to t0)
 * - Bin-based: \c fit \c fgb+100 \c lgb-500 (bins relative to good bins)
 *
 * <b>Search Priority:</b>
 * -# RUN block time-based fit range
 * -# RUN block bin-based fit range (converted to time)
 * -# GLOBAL block time-based fit range
 * -# GLOBAL block bin-based fit range (converted to time)
 * -# Default to data range (fgb to lgb) with warning
 *
 * <b>Bin-to-Time Conversion:</b>
 * \f[
 * t_{\rm start} = ({\rm fgb} + n_0 - t_0) \times \Delta t
 * \f]
 * \f[
 * t_{\rm end} = ({\rm lgb} - n_1 - t_0) \times \Delta t
 * \f]
 *
 * where Δt = fTimeResolution (raw, before RRF packing).
 *
 * <b>Result:</b>
 * Sets fFitStartTime and fFitEndTime in microseconds relative to t0.
 * These define the range used in χ² calculation.
 *
 * \param globalBlock Pointer to GLOBAL block for fallback fit range settings
 *
 * \note The converted time values are written back to the MSR data structures
 *       for log file output consistency.
 *
 * \warning If no fit range is specified, the full data range is used with
 *          a warning message. This may not be appropriate for all analyses.
 *
 * \see GetProperDataRange(), SetFitRangeBin(), CalcNoOfFitBins()
 */
void PRunSingleHistoRRF::GetProperFitRange(PMsrGlobalBlock *globalBlock)
{
  // set fit start/end time; first check RUN Block
  fFitStartTime = fRunInfo->GetFitRange(0);
  fFitEndTime   = fRunInfo->GetFitRange(1);
  // if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now
  if (fRunInfo->IsFitRangeInBin()) {
    fFitStartTime = (fGoodBins[0] + fRunInfo->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt
    fFitEndTime = (fGoodBins[1] - fRunInfo->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution;   // (lgb-n1-t0)*dt
    // write these times back into the data structure. This way it is available when writting the log-file
    fRunInfo->SetFitRange(fFitStartTime, 0);
    fRunInfo->SetFitRange(fFitEndTime, 1);
  }
  if (fFitStartTime == PMUSR_UNDEFINED) { // fit start/end NOT found in the RUN block, check GLOBAL block
    fFitStartTime = globalBlock->GetFitRange(0);
    fFitEndTime   = globalBlock->GetFitRange(1);
    // if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now
    if (globalBlock->IsFitRangeInBin()) {
      fFitStartTime = (fGoodBins[0] + globalBlock->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt
      fFitEndTime = (fGoodBins[1] - globalBlock->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution;   // (lgb-n1-t0)*dt
      // write these times back into the data structure. This way it is available when writting the log-file
      globalBlock->SetFitRange(fFitStartTime, 0);
      globalBlock->SetFitRange(fFitEndTime, 1);
    }
  }
  if ((fFitStartTime == PMUSR_UNDEFINED) || (fFitEndTime == PMUSR_UNDEFINED)) {
    fFitStartTime = (fGoodBins[0] - fT0s[0]) * fTimeResolution; // (fgb-t0)*dt
    fFitEndTime = (fGoodBins[1] - fT0s[0]) * fTimeResolution;   // (lgb-t0)*dt
    std::cerr << ">> PRunSingleHistoRRF::GetProperFitRange(): **WARNING** Couldn't get fit start/end time!" << std::endl;
    std::cerr << ">>    Will set it to fgb/lgb which given in time is: " << fFitStartTime << "..." << fFitEndTime << " (usec)" << std::endl;
  }
}

//--------------------------------------------------------------------------
// GetMainFrequency (private)
//--------------------------------------------------------------------------
/**
 * \brief Finds the dominant precession frequency in raw histogram data.
 *
 * Performs Fourier analysis on the raw histogram to identify the main muon
 * spin precession frequency. This is essential for:
 * - Determining optimal N₀ estimation window (integer oscillation cycles)
 * - Calculating suggested "optimal packing" for user information
 * - Validating that the RRF frequency is appropriate
 *
 * <b>Algorithm:</b>
 * -# Create ROOT TH1F histogram from data in good bin range [fGoodBins[0], fGoodBins[1]]
 * -# Set histogram binning to match time structure
 * -# Apply Fourier transform using PFourier class
 * -# Use strong apodization (windowing) to reduce spectral leakage
 * -# Search power spectrum for maximum above 10 MHz (ignores DC component)
 * -# Return frequency at maximum power
 *
 * <b>Frequency Search Constraints:</b>
 * - Ignores frequencies below 10 MHz (DC and low-frequency noise)
 * - Searches for local maximum (rising then falling power)
 * - Returns 0 if no clear maximum found
 *
 * <b>Output Information:</b>
 * The method prints diagnostic information:
 * - Detected maximum frequency (MHz)
 * - Suggested optimal packing for ~5-8 points per cycle
 *
 * \param data Raw histogram data vector (counts per bin)
 *
 * \return Maximum frequency in MHz, or 0 if no maximum found above 10 MHz.
 *
 * \note The frequency is used to determine the N₀ estimation window:
 *       window = ceil(fN0EstimateEndTime × freqMax / 2π) × (2π / freqMax)
 *       This ensures an integer number of complete oscillation cycles.
 *
 * \see EstimateN0(), PrepareFitData()
 */
Double_t PRunSingleHistoRRF::GetMainFrequency(PDoubleVector &data)
{
  Double_t freqMax = 0.0;

  // create histo
  Double_t startTime = (fGoodBins[0]-fT0s[0]) * fTimeResolution;
  Int_t noOfBins = fGoodBins[1]-fGoodBins[0]+1;
  std::unique_ptr<TH1F> histo = std::make_unique<TH1F>("data", "data", noOfBins, startTime-fTimeResolution/2.0, startTime+fTimeResolution/2.0+noOfBins*fTimeResolution);
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    histo->SetBinContent(i-fGoodBins[0]+1, data[i]);
  }

  // Fourier transform
  std::unique_ptr<PFourier> ft = std::make_unique<PFourier>(histo.get(), FOURIER_UNIT_FREQ);
  ft->Transform(F_APODIZATION_STRONG);

  // find frequency maximum
  TH1F *power = ft->GetPowerFourier();
  Double_t maxFreqVal = 0.0;
  for (Int_t i=1; i<power->GetNbinsX(); i++) {
    // ignore dc part at 0 frequency
    if (i<power->GetNbinsX()-1) {
      if (power->GetBinContent(i)>power->GetBinContent(i+1))
        continue;
    }
    // ignore everything below 10 MHz
    if (power->GetBinCenter(i) < 10.0)
      continue;
    // check for maximum
    if (power->GetBinContent(i) > maxFreqVal) {
      maxFreqVal = power->GetBinContent(i);
      freqMax = power->GetBinCenter(i);
    }
  }

  // clean up
  if (power)
    delete power;

  return freqMax;
}

//--------------------------------------------------------------------------
// EstimateN0 (private)
//--------------------------------------------------------------------------
/**
 * \brief Estimates initial normalization N₀ from lifetime-corrected histogram data.
 *
 * Calculates N₀ by averaging the lifetime-corrected signal M(t) over an initial
 * time window. The window is chosen to span an integer number of complete
 * precession cycles to minimize bias from oscillations.
 *
 * <b>N₀ Estimation Formula:</b>
 * \f[
 * N_0 = \frac{1}{n} \sum_{i=0}^{n-1} M(t_i)
 * \f]
 *
 * where M(t) = [N(t) - B] × exp(+t/τ_μ) is the lifetime-corrected histogram
 * stored in fM vector.
 *
 * <b>Window Determination:</b>
 * The end bin is calculated to include an integer number of oscillation cycles:
 * \f[
 * n_{\rm end} = {\rm round}\left( \left\lceil \frac{T \cdot f_{\rm max}}{2\pi} \right\rceil \cdot \frac{2\pi}{f_{\rm max} \cdot \Delta t} \right)
 * \f]
 *
 * where:
 * - T = fN0EstimateEndTime (default 1.0 μs)
 * - f_max = main precession frequency from GetMainFrequency()
 * - Δt = fTimeResolution
 *
 * <b>Error Estimation:</b>
 * \f[
 * \sigma_{N_0} = \frac{\sqrt{\sum_{i=0}^{n-1} w_i^2 \sigma_{M_i}^2}}{\sum_{i=0}^{n-1} w_i}
 * \f]
 *
 * where w_i = 1/σ_M_i² are the weights stored in fW vector.
 *
 * <b>Output Information:</b>
 * Prints diagnostic message: "N₀ = value (error)"
 *
 * \param errN0 [out] Reference to store the estimated error on N₀
 * \param freqMax Main precession frequency (MHz) from GetMainFrequency()
 *
 * \return Estimated N₀ value (counts × exp(+t/τ))
 *
 * \note The simple average (rather than weighted average) is used for N₀ itself,
 *       while the error uses the full weight information.
 *
 * \see GetMainFrequency(), PrepareFitData()
 */
Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax)
{
  // endBin is estimated such that the number of full cycles (according to the maximum frequency of the data)
  // is approximately the time fN0EstimateEndTime.
  Int_t endBin = static_cast<Int_t>(round(ceil(fN0EstimateEndTime*freqMax/TMath::TwoPi()) * (TMath::TwoPi()/freqMax) / fTimeResolution));

  Double_t n0 = 0.0;
  Double_t wN = 0.0;
  for (Int_t i=0; i<endBin; i++) {
//    n0 += fW[i]*fM[i];
    n0 += fM[i];
    wN += fW[i];
  }
//  n0 /= wN;
  n0 /= endBin;

  errN0 = 0.0;
  for (Int_t i=0; i<endBin; i++) {
    errN0 += fW[i]*fW[i]*fMerr[i]*fMerr[i];
  }
  errN0 = sqrt(errN0)/wN;

  std::cout << "info> PRunSingleHistoRRF::EstimateN0(): N0=" << n0 << "(" << errN0 << ")" << std::endl;

  return n0;
}

//--------------------------------------------------------------------------
// EstimatBkg (private)
//--------------------------------------------------------------------------
/**
 * \brief Estimates background level and error from pre-t0 histogram bins.
 *
 * Calculates the background (random coincidences, accidentals) from a
 * specified range of histogram bins, typically before t0. The background
 * range is adjusted to span an integer number of accelerator RF cycles
 * for proper averaging.
 *
 * <b>Accelerator RF Periods:</b>
 * - PSI: ACCEL_PERIOD_PSI ns
 * - RAL: ACCEL_PERIOD_RAL ns
 * - TRIUMF: ACCEL_PERIOD_TRIUMF ns
 * - Unknown: no RF adjustment
 *
 * <b>Algorithm:</b>
 * -# Get background range [start, end] from RUN block
 * -# Swap start/end if in wrong order (with message)
 * -# Determine accelerator RF period from institute name
 * -# Adjust end to span integer number of RF cycles
 * -# Validate range is within histogram bounds
 * -# Calculate mean background:
 *    \f$ B = \frac{1}{n}\sum_{i={\rm start}}^{{\rm end}} N(i) \f$
 * -# Calculate background error (standard deviation):
 *    \f$ \sigma_B = \sqrt{\frac{1}{n-1}\sum_{i={\rm start}}^{{\rm end}} (N(i) - B)^2} \f$
 * -# Store results in fBackground and fBkgErr
 * -# Update RUN block with estimated background value
 *
 * <b>Output Information:</b>
 * Prints diagnostic messages:
 * - Adjusted background range (if RF adjustment applied)
 * - Background value and error: "fBackground = value (error)"
 *
 * \param histoNo Forward histogram number (for error message context)
 *
 * \return true on success, false on error:
 *         - Background start out of histogram bounds
 *         - Background end out of histogram bounds
 *
 * \note The RF cycle adjustment ensures proper averaging over the
 *       accelerator bunch structure, reducing systematic bias.
 *
 * \see PrepareFitData()
 */
Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo)
{
  Double_t beamPeriod = 0.0;

  // check if data are from PSI, RAL, or TRIUMF
  if (fRunInfo->GetInstitute()->Contains("psi"))
    beamPeriod = ACCEL_PERIOD_PSI;
  else if (fRunInfo->GetInstitute()->Contains("ral"))
    beamPeriod = ACCEL_PERIOD_RAL;
  else if (fRunInfo->GetInstitute()->Contains("triumf"))
    beamPeriod = ACCEL_PERIOD_TRIUMF;
  else
    beamPeriod = 0.0;

  // check if start and end are in proper order
  UInt_t start = fRunInfo->GetBkgRange(0);
  UInt_t end   = fRunInfo->GetBkgRange(1);
  if (end < start) {
    std::cout << std::endl << "PRunSingleHistoRRF::EstimatBkg(): end = " << end << " > start = " << start << "! Will swap them!";
    UInt_t keep = end;
    end = start;
    start = keep;
  }

  // calculate proper background range
  if (beamPeriod != 0.0) {
    Double_t timeBkg = static_cast<Double_t>(end-start)*fTimeResolution; // length of the background intervall in time
    UInt_t fullCycles = static_cast<UInt_t>(timeBkg/beamPeriod); // how many proton beam cylces can be placed within the proposed background intervall
    // correct the end of the background intervall such that the background is as close as possible to a multiple of the proton cylce
    end = start + static_cast<UInt_t>((fullCycles*beamPeriod)/fTimeResolution);
    std::cout << std::endl << "PRunSingleHistoRRF::EstimatBkg(): Background " << start << ", " << end;
    if (end == start)
      end = fRunInfo->GetBkgRange(1);
  }

  // check if start is within histogram bounds
  if (start >= fForward.size()) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!";
    std::cerr << std::endl << ">> histo lengths    = " << fForward.size();
    std::cerr << std::endl << ">> background start = " << start;
    std::cerr << std::endl;
    return false;
  }

  // check if end is within histogram bounds
  if (end >= fForward.size()) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!";
    std::cerr << std::endl << ">> histo lengths  = " << fForward.size();
    std::cerr << std::endl << ">> background end = " << end;
    std::cerr << std::endl;
    return false;
  }

  // calculate background
  Double_t bkg    = 0.0;

  // forward
  for (UInt_t i=start; i<end; i++)
    bkg += fForward[i];
  bkg /= static_cast<Double_t>(end - start + 1);

  fBackground = bkg;  // keep background (per bin)

  bkg = 0.0;
  for (UInt_t i=start; i<end; i++)
    bkg += pow(fForward[i]-fBackground, 2.0);
  fBkgErr = sqrt(bkg/(static_cast<Double_t>(end - start)));

  std::cout << std::endl << "info> fBackground=" << fBackground << "(" << fBkgErr << ")" << std::endl;

  fRunInfo->SetBkgEstimated(fBackground, 0);

  return true;
}