musrfit/src/classes/PRunSingleHisto.cpp

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

  PRunSingleHisto.cpp

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

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

/***************************************************************************
 *   Copyright (C) 2007-2026 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 <TString.h>
#include <TObjArray.h>
#include <TObjString.h>

#include "PMusr.h"
#include "PRunSingleHisto.h"

//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
 * \brief Default constructor for single histogram fitting class.
 *
 * Initializes all member variables to safe default values:
 * - fScaleN0AndBkg = true (normalize N₀ and background to 1/ns)
 * - fPacking = -1 (invalid until set from MSR file)
 * - fBackground = 0 (will be estimated or set from MSR file)
 * - fStartTimeBin / fEndTimeBin = -1 (calculated from fit range)
 * - fGoodBins[0,1] = -1 (calculated from data range)
 *
 * \warning This constructor creates an invalid object until initialized
 *          with MSR file data. Use the full constructor for normal operation.
 */
PRunSingleHisto::PRunSingleHisto() : PRunBase()
{
  fScaleN0AndBkg = true;
  fNoOfFitBins  = 0;
  fBackground = 0;
  fPacking = -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;

  fStartTimeBin = -1;
  fEndTimeBin   = -1;
}

//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
 * \brief Main constructor for single histogram fitting and viewing.
 *
 * Constructs a fully initialized single histogram run object by:
 * -# Extracting packing value from RUN block (or falling back to GLOBAL block)
 * -# Determining if N₀ and background should be scaled to 1/ns
 * -# Calling PrepareData() to load and process histogram data
 * -# Setting up fit ranges and background estimation
 *
 * \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 is calculated only at data points (for viewing);
 *                   if false, theory uses finer time grid (8× data resolution)
 *
 * \warning Packing MUST be specified either in the RUN block or GLOBAL block.
 *          If packing is not found, the constructor sets fValid=false and returns.
 *
 * \note After construction, check IsValid() to ensure initialization succeeded.
 *
 * \see PrepareData(), IsScaleN0AndBkg()
 */
PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) :
  PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData)
{
  fScaleN0AndBkg = IsScaleN0AndBkg();
  fNoOfFitBins  = 0;
  fBackground = 0;

  fPacking = fRunInfo->GetPacking();
  if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block
    fPacking = fMsrInfo->GetMsrGlobal()->GetPacking();
  }
  if (fPacking == -1) { // this should NOT happen, somethin is severely wrong
    std::cerr << std::endl << ">> PRunSingleHisto::PRunSingleHisto: **SEVERE ERROR**: Couldn't find any packing information!";
    std::cerr << std::endl << ">> This is very bad :-(, will quit ...";
    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;

  fStartTimeBin = -1;
  fEndTimeBin   = -1;

  if (!PrepareData()) {
    std::cerr << std::endl << ">> PRunSingleHisto::PRunSingleHisto: **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 single histogram fitting class.
 *
 * Cleans up dynamically allocated memory:
 * - Clears the forward histogram data vector
 * - Base class destructor handles theory objects and other shared resources
 */
PRunSingleHisto::~PRunSingleHisto()
{
  fForward.clear();
}

//--------------------------------------------------------------------------
// CalcChiSquare (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates χ² between data and theory (least-squares fit metric).
 *
 * Computes the standard chi-square goodness-of-fit statistic:
 * \f[
 * \chi^2 = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[N_i - N_{\rm theo}(t_i)]^2}{\sigma_i^2}
 * \f]
 *
 * where the theory function is:
 * \f[
 * N_{\rm theo}(t) = N_0 e^{-t/\tau_\mu} [1 + P(t)] + B
 * \f]
 *
 * <b>Algorithm:</b>
 * -# Extract N₀ from parameter vector or evaluate as a function
 * -# Extract muon lifetime τ (defaults to PMUON_LIFETIME if not fitted)
 * -# Extract background B (from fit parameter, fixed value, or estimated range)
 * -# Evaluate all user-defined functions in FUNCTIONS block
 * -# Pre-calculate theory at t=1.0 to initialize LF/user functions (thread-safe)
 * -# Loop over fit range bins [fStartTimeBin, fEndTimeBin) using OpenMP parallelization
 * -# Accumulate χ² with reduction across threads
 * -# Apply correction factor if fScaleN0AndBkg is true
 *
 * <b>N₀ Parameter vs. Function Handling:</b>
 * - If norm parameter number < MSR_PARAM_FUN_OFFSET: N₀ is a fit parameter
 * - If norm parameter number ≥ MSR_PARAM_FUN_OFFSET: N₀ is a user-defined function
 *
 * <b>OpenMP Parallelization:</b>
 * - Dynamic scheduling with chunk size = (N_bins / N_processors), minimum 10
 * - Private variables per thread: i, time, diff
 * - Reduction performed on chisq sum
 *
 * <b>Scaling Correction:</b>
 * If fScaleN0AndBkg is true, χ² is multiplied by:
 * \f[
 * \text{correction} = \text{packing} \times (t_{\rm res} \times 1000)
 * \f]
 * This accounts for the fact that data scales like pack×t_res, but errors
 * scale like √(pack×t_res), ensuring correct χ² when normalizing to 1/ns.
 *
 * \param par Parameter vector from MINUIT2 optimizer (1-based indexing in MSR file,
 *            but 0-based in this vector)
 *
 * \return Chi-square value for the current parameter set
 *
 * \see CalcChiSquareExpected(), CalcMaxLikelihood(), PTheory::Func()
 */
Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
{
  Double_t chisq     = 0.0;
  Double_t diff      = 0.0;

  Double_t N0 = 0.0;

  // check if norm is a parameter or a function
  if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
    N0 = par[fRunInfo->GetNormParamNo()-1];
  } else { // norm is a function
    // get function number
    UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
    // evaluate function
    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
  }

  // get tau
  Double_t tau;
  if (fRunInfo->GetLifetimeParamNo() != -1)
    tau = par[fRunInfo->GetLifetimeParamNo()-1];
  else
    tau = PMUON_LIFETIME;

  // get background
  Double_t bkg;
  if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted
    if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
      bkg = fBackground;
    } else { // fixed bkg given
      bkg = fRunInfo->GetBkgFix(0);
    }
  } else { // bkg fitted
    bkg = par[fRunInfo->GetBkgFitParamNo()-1];
  }

  // 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) -
          (N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg);
    chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
  }

  // the correction factor is need since the data scales like pack*t_res,
  // whereas the error scales like sqrt(pack*t_res)
  if (fScaleN0AndBkg)
    chisq *= fPacking * (fTimeResolution * 1.0e3);

  return chisq;
}

//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates expected χ² using theory as variance (alternative fit metric).
 *
 * Computes chi-square using the expected variance (theory value) instead of
 * observed variance. This is sometimes called the "Neyman χ²" or "expected χ²":
 * \f[
 * \chi^2_{\rm exp} = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[N_i - N_{\rm theo}(t_i)]^2}{N_{\rm theo}(t_i)}
 * \f]
 *
 * <b>Difference from Standard χ²:</b>
 * - Standard χ²: variance = σ²ᵢ (from observed data)
 * - Expected χ²: variance = N_theo(tᵢ) (from theory prediction)
 *
 * This metric can be useful when:
 * - Theory predictions are more reliable than data errors
 * - Data contains zero or very low counts (standard χ² undefined)
 * - Testing model consistency against expected distribution
 *
 * <b>Algorithm:</b>
 * -# Extract N₀ from parameter vector or evaluate as a function
 * -# Extract muon lifetime τ (defaults to PMUON_LIFETIME if not fitted)
 * -# Extract background B (from fit parameter, fixed value, or estimated range)
 * -# Evaluate all user-defined functions in FUNCTIONS block
 * -# Pre-calculate theory at t=1.0 to initialize LF/user functions (thread-safe)
 * -# Loop over fit range bins [fStartTimeBin, fEndTimeBin) using OpenMP parallelization
 * -# Accumulate χ²_exp with reduction across threads
 * -# Apply correction factor if fScaleN0AndBkg is true
 *
 * <b>OpenMP Parallelization:</b>
 * - Dynamic scheduling with chunk size = (N_bins / N_processors), minimum 10
 * - Private variables per thread: i, time, theo, diff
 * - Reduction performed on chisq sum
 *
 * \param par Parameter vector from MINUIT2 optimizer
 *
 * \return Expected chi-square value for the current parameter set
 *
 * \see CalcChiSquare(), CalcMaxLikelihood()
 */
Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
  Double_t chisq = 0.0;
  Double_t diff  = 0.0;
  Double_t theo  = 0.0;

  Double_t N0 = 0.0;

  // check if norm is a parameter or a function
  if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
    N0 = par[fRunInfo->GetNormParamNo()-1];
  } else { // norm is a function
    // get function number
    UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
    // evaluate function
    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
  }

  // get tau
  Double_t tau;
  if (fRunInfo->GetLifetimeParamNo() != -1)
    tau = par[fRunInfo->GetLifetimeParamNo()-1];
  else
    tau = PMUON_LIFETIME;

  // get background
  Double_t bkg;
  if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted
    if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
      bkg = fBackground;
    } else { // fixed bkg given
      bkg = fRunInfo->GetBkgFix(0);
    }
  } else { // bkg fitted
    bkg = par[fRunInfo->GetBkgFitParamNo()-1];
  }

  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    Int_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,theo,diff) schedule(dynamic,chunk) reduction(+:chisq)
  #endif
  for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
    time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
    theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg;
    diff = fData.GetValue()->at(i) - theo;
    chisq += diff*diff / theo;
  }

  // the correction factor is need since the data scales like pack*t_res,
  // whereas the error scales like sqrt(pack*t_res)
  if (fScaleN0AndBkg)
    chisq *= fPacking * (fTimeResolution * 1.0e3);

  return chisq;
}

//--------------------------------------------------------------------------
// CalcMaxLikelihood (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates -2 log(maximum likelihood) for Poisson-distributed histogram data.
 *
 * Computes the negative log-likelihood assuming Poisson statistics for each bin.
 * This is the preferred fit metric for low-count data where Gaussian approximations
 * break down. The likelihood function is:
 * \f[
 * -2\ln\mathcal{L} = 2 \sum_{i} \left[ N_{\rm theo}(t_i) - N_i + N_i \ln\frac{N_i}{N_{\rm theo}(t_i)} \right]
 * \f]
 *
 * This is derived from the Poisson probability:
 * \f[
 * P(N_i | N_{\rm theo}) = \frac{N_{\rm theo}^{N_i} e^{-N_{\rm theo}}}{N_i!}
 * \f]
 *
 * The factor of 2 makes -2ln(L) asymptotically distributed as χ² for large N,
 * allowing use of standard error estimation from MINUIT.
 *
 * <b>Algorithm:</b>
 * -# Extract N₀ from parameter vector or evaluate as a function
 * -# Extract muon lifetime τ (defaults to PMUON_LIFETIME if not fitted)
 * -# Extract background B (from fit parameter, fixed value, or estimated range)
 * -# Evaluate all user-defined functions in FUNCTIONS block
 * -# Pre-calculate theory at t=1.0 to initialize LF/user functions (thread-safe)
 * -# Calculate normalizer = packing × t_res × 1000 (if fScaleN0AndBkg is true)
 * -# Loop over fit range bins [fStartTimeBin, fEndTimeBin) using OpenMP parallelization
 * -# For each bin:
 *    - Calculate theory N_theo(t)
 *    - If N_theo ≤ 0: skip bin with warning (negative theory is unphysical)
 *    - If N_data > 10⁻⁹: add (theo - data) + data×ln(data/theo)
 *    - If N_data ≈ 0: add (theo - data) only (limit as data→0)
 * -# Accumulate -2ln(L) with reduction across threads
 * -# Apply normalizer scaling
 *
 * <b>Edge Cases:</b>
 * - <b>Zero data (Nᵢ = 0):</b> Uses limit: -2ln(L) → 2×N_theo
 * - <b>Negative theory:</b> Skips bin and prints warning (should not occur with valid parameters)
 * - <b>Data threshold:</b> Uses 10⁻⁹ to distinguish zero from non-zero data
 *
 * <b>OpenMP Parallelization:</b>
 * - Dynamic scheduling with chunk size = (N_bins / N_processors), minimum 10
 * - Private variables per thread: i, time, theo, data
 * - Reduction performed on mllh sum (reduction(+:mllh))
 *
 * <b>When to Use Maximum Likelihood vs. χ²:</b>
 * - <b>Use likelihood:</b> Low count rates (< 100 counts/bin), asymmetric errors
 * - <b>Use χ²:</b> High count rates (> 100 counts/bin), Gaussian regime
 *
 * \param par Parameter vector from MINUIT2 optimizer
 *
 * \return -2 × log(maximum likelihood) for the current parameter set
 *
 * \see CalcChiSquare(), CalcMaxLikelihoodExpected()
 * \see PDG Review of Particle Physics: Statistics section (http://pdg.lbl.gov)
 */
Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
{
  Double_t mllh = 0.0; // maximum log likelihood assuming poisson distribution for the single bin

  Double_t N0;

  // check if norm is a parameter or a function
  if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
    N0 = par[fRunInfo->GetNormParamNo()-1];
  } else { // norm is a function
    // get function number
    Int_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
    // evaluate function
    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
  }

  // get tau
  Double_t tau;
  if (fRunInfo->GetLifetimeParamNo() != -1)
    tau = par[fRunInfo->GetLifetimeParamNo()-1];
  else
    tau = PMUON_LIFETIME;

  // get background
  Double_t bkg;
  if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted
    if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
      bkg = fBackground;
    } else { // fixed bkg given
      bkg = fRunInfo->GetBkgFix(0);
    }
  } else { // bkg fitted
    bkg = par[fRunInfo->GetBkgFitParamNo()-1];
  }

  // 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 maximum log likelihood
  Double_t theo;
  Double_t data;
  Double_t time(1.0);
  Int_t i;

  // norm is needed since there is no simple scaling like in chisq case to get the correct Max.Log.Likelihood value when normlizing N(t) to 1/ns
  Double_t normalizer = 1.0;

  if (fScaleN0AndBkg)
    normalizer = fPacking * (fTimeResolution * 1.0e3);

  // 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,theo,data) schedule(dynamic,chunk) reduction(+:mllh)
  #endif
  for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
    time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
    // calculate theory for the given parameter set
    theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg;

    data = fData.GetValue()->at(i);

    if (theo <= 0.0) {
      std::cerr << ">> PRunSingleHisto::CalcMaxLikelihood: **WARNING** NEGATIVE theory!!" << std::endl;
      continue;
    }

    if (data > 1.0e-9) {
      mllh += (theo-data) + data*log(data/theo);
    } else {
      mllh += (theo-data);
    }
  }

  return normalizer*2.0*mllh;
}

//--------------------------------------------------------------------------
// CalcMaxLikelihoodExpected (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates expected -2 log(maximum likelihood) using G-test formulation.
 *
 * Computes an alternative form of the Poisson likelihood using only the data×ln(data/theo)
 * term. This is related to the G-test (likelihood ratio test) and represents the
 * "expected" contribution to the likelihood:
 * \f[
 * -2\ln\mathcal{L}_{\rm exp} = 2 \sum_{i} N_i \ln\frac{N_i}{N_{\rm theo}(t_i)}
 * \f]
 *
 * <b>Difference from CalcMaxLikelihood():</b>
 * - <b>Full likelihood:</b> includes (theo - data) + data×ln(data/theo)
 * - <b>Expected likelihood:</b> includes only data×ln(data/theo)
 *
 * The omitted (theo - data) term represents the "prior" expectation and is
 * constant for a given theory. This formulation is sometimes used in:
 * - G-test for goodness-of-fit (likelihood ratio test)
 * - Comparing relative likelihoods between models
 *
 * <b>Algorithm:</b>
 * -# Extract N₀, τ, and background B (same as CalcMaxLikelihood)
 * -# Evaluate all user-defined functions in FUNCTIONS block
 * -# Pre-calculate theory at t=1.0 to initialize LF/user functions (thread-safe)
 * -# Calculate normalizer = packing × t_res × 1000 (if fScaleN0AndBkg is true)
 * -# Loop over fit range bins using OpenMP parallelization
 * -# For each bin with N_data > 10⁻⁹:
 *    - Calculate theory N_theo(t)
 *    - Add data × ln(data/theo) to likelihood sum
 * -# Skip bins with N_data ≈ 0 (zero contribution to expected likelihood)
 * -# Apply normalizer × 2.0 scaling
 *
 * \warning The comment "is this correct?? needs to be checked. See G-test"
 *          in the code indicates this implementation may need verification.
 *
 * <b>OpenMP Parallelization:</b>
 * - Dynamic scheduling with chunk size = (N_bins / N_processors), minimum 10
 * - Private variables per thread: i, time, theo, data
 * - Reduction performed on mllh sum
 *
 * \param par Parameter vector from MINUIT2 optimizer
 *
 * \return -2 × log(expected likelihood) for the current parameter set
 *
 * \see CalcMaxLikelihood(), G-test (likelihood ratio test)
 */
Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector<Double_t>& par)
{
  Double_t mllh = 0.0; // maximum log likelihood assuming poisson distribution for the single bin

  Double_t N0;

  // check if norm is a parameter or a function
  if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
    N0 = par[fRunInfo->GetNormParamNo()-1];
  } else { // norm is a function
    // get function number
    Int_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
    // evaluate function
    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
  }

  // get tau
  Double_t tau;
  if (fRunInfo->GetLifetimeParamNo() != -1)
    tau = par[fRunInfo->GetLifetimeParamNo()-1];
  else
    tau = PMUON_LIFETIME;

  // get background
  Double_t bkg;
  if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted
    if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
      bkg = fBackground;
    } else { // fixed bkg given
      bkg = fRunInfo->GetBkgFix(0);
    }
  } else { // bkg fitted
    bkg = par[fRunInfo->GetBkgFitParamNo()-1];
  }

  // 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 maximum log likelihood
  Double_t theo;
  Double_t data;
  Double_t time(1.0);
  Int_t i;

  // norm is needed since there is no simple scaling like in chisq case to get the correct Max.Log.Likelihood value when normlizing N(t) to 1/ns
  Double_t normalizer = 1.0;

  if (fScaleN0AndBkg)
    normalizer = fPacking * (fTimeResolution * 1.0e3);

  // 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,theo,data) schedule(dynamic,chunk) reduction(+:mllh)
  #endif
  for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
    time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
    // calculate theory for the given parameter set
    theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg;

    data = fData.GetValue()->at(i);

    if (theo <= 0.0) {
      std::cerr << ">> PRunSingleHisto::CalcMaxLikelihood: **WARNING** NEGATIVE theory!!" << std::endl;
      continue;
    }

    if (data > 1.0e-9) { // is this correct?? needs to be checked. See G-test
      mllh += data*log(data/theo);
    }
  }

  return normalizer*2.0*mllh;
}

//--------------------------------------------------------------------------
// CalcTheory (public)
//--------------------------------------------------------------------------
/**
 * \brief Calculates theory curve N(t) for the current parameter values.
 *
 * Evaluates the single histogram theory function:
 * \f[
 * N(t) = N_0 e^{-t/\tau_\mu} [1 + P(t)] + B
 * \f]
 *
 * for all time bins in the data set, storing results in fData.fTheory.
 * This is used for:
 * - Displaying fitted theory curves in plots
 * - Calculating residuals (data - theory)
 * - Exporting theory predictions
 *
 * <b>Algorithm:</b>
 * -# Extract current parameter values from MSR parameter list
 * -# Determine N₀ (from parameter or function evaluation)
 * -# Determine muon lifetime τ (from parameter or default PMUON_LIFETIME)
 * -# Determine background B (from fit parameter, fixed value, or estimate)
 * -# Evaluate all user-defined functions in FUNCTIONS block
 * -# Loop over all data bins (not just fit range):
 *    - Calculate time t for bin i
 *    - Evaluate P(t) = fTheory->Func(t, par, fFuncValues)
 *    - Calculate N(t) and append to theory vector
 * -# Clean up temporary parameter vector
 *
 * <b>Time Grid:</b>
 * - Start time: fData.GetDataTimeStart()
 * - Time step: fData.GetDataTimeStep()
 * - Number of points: fData.GetValue()->size()
 *
 * \note Theory is calculated for the entire data range, not just the fit range,
 *       to enable full visualization of the model.
 *
 * \see PRunDataHandler::AppendTheoryValue(), PTheory::Func()
 */
void PRunSingleHisto::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 asymmetry
  Double_t N0;
  // check if norm is a parameter or a function
  if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
    N0 = par[fRunInfo->GetNormParamNo()-1];
  } else { // norm is a function
    // get function number
    Int_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
    // evaluate function
    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
  }

  // get tau
  Double_t tau;
  if (fRunInfo->GetLifetimeParamNo() != -1)
    tau = par[fRunInfo->GetLifetimeParamNo()-1];
  else
    tau = PMUON_LIFETIME;

  // get background
  Double_t bkg;
  if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted
    if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
      bkg = fBackground;
    } else { // fixed bkg given
      bkg = fRunInfo->GetBkgFix(0);
    }
  } else { // bkg fitted
    bkg = par[fRunInfo->GetBkgFitParamNo()-1];
  }

  // 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(N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg);
  }

  // clean up
  par.clear();
}

//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
 * \brief Returns the number of bins in the current fit range.
 *
 * Calculates (if not already done) and returns the number of data bins
 * that will be included in the χ² or likelihood calculation. This is
 * determined by the fit range [fFitStartTime, fFitEndTime] and the
 * data time grid.
 *
 * The calculation is performed by CalcNoOfFitBins(), which sets:
 * - fStartTimeBin: first bin index in fit range
 * - fEndTimeBin: one past last bin index in fit range
 * - fNoOfFitBins = fEndTimeBin - fStartTimeBin
 *
 * \return Number of bins in the fit range (degrees of freedom = N_bins - N_params)
 *
 * \see CalcNoOfFitBins(), SetFitRangeBin()
 */
UInt_t PRunSingleHisto::GetNoOfFitBins()
{
  CalcNoOfFitBins();

  return fNoOfFitBins;
}

//--------------------------------------------------------------------------
// SetFitRangeBin (public)
//--------------------------------------------------------------------------
/**
 * \brief Dynamically changes the fit range from COMMAND block instructions.
 *
 * Parses and applies a FIT_RANGE command to modify the fit range on the fly,
 * typically used during interactive fitting sessions or systematic scans.
 *
 * <b>Syntax (in COMMAND block):</b>
 * \code
 * FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]]
 * \endcode
 *
 * where:
 * - <b>fgb:</b> First good bin (start of fit range)
 * - <b>lgb:</b> Last good bin (end of fit range)
 * - <b>+nXY / -nXY:</b> Optional offsets to shift the range (+ extends, - contracts)
 * - <b>Multiple pairs:</b> If N+1 pairs given, they apply to each of N RUN blocks
 *
 * <b>Two modes:</b>
 * -# <b>Single pair:</b> `FIT_RANGE fgb lgb` applies to all runs
 * -# <b>Per-run pairs:</b> `FIT_RANGE fgb₀ lgb₀ fgb₁ lgb₁ ...` applies
 *    pair i to RUN block i
 *
 * <b>Algorithm:</b>
 * -# Tokenize the fitRange string by spaces/tabs
 * -# If 3 tokens (FIT_RANGE + 2 values): apply to this run
 * -# If >3 tokens and odd number: extract pair for this run's index (fRunNo)
 * -# Parse offsets from + or - characters in fgb/lgb strings
 * -# Calculate new fFitStartTime and fFitEndTime:
 *    - fFitStartTime = (fGoodBins[0] + offset - t0) × t_res
 *    - fFitEndTime = (fGoodBins[1] - offset - t0) × t_res
 *
 * <b>Example:</b>
 * \code
 * FIT_RANGE 100+10 500-20  # Fit from bin 110 to bin 480 (applying offsets)
 * \endcode
 *
 * \param fitRange String from COMMAND block containing FIT_RANGE specification
 *
 * \note Errors in parsing (wrong number of tokens) are reported to std::cerr
 *       and the command is ignored.
 *
 * \see CalcNoOfFitBins(), GetProperFitRange()
 */
void PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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 the number of bins in the fit range and caches bin indices.
 *
 * Converts the fit time range [fFitStartTime, fFitEndTime] to bin indices
 * [fStartTimeBin, fEndTimeBin) and computes the total number of fit bins.
 *
 * <b>Algorithm:</b>
 * -# Calculate start bin: \f$ \lceil \frac{t_{\rm start} - t_{\rm data,0}}{\Delta t} \rceil \f$
 * -# Clamp fStartTimeBin to [0, N_data)
 * -# Calculate end bin: \f$ \lfloor \frac{t_{\rm end} - t_{\rm data,0}}{\Delta t} \rfloor + 1 \f$
 * -# Clamp fEndTimeBin to [0, N_data]
 * -# Compute fNoOfFitBins = fEndTimeBin - fStartTimeBin (or 0 if invalid)
 *
 * where:
 * - t_data,0 = fData.GetDataTimeStart() (time of first data bin)
 * - Δt = fData.GetDataTimeStep() (time bin width after packing)
 *
 * <b>Edge Cases:</b>
 * - If fStartTimeBin < 0: clamped to 0
 * - If fEndTimeBin > N_data: clamped to N_data
 * - If fEndTimeBin ≤ fStartTimeBin: fNoOfFitBins = 0 (invalid range)
 *
 * \note This method is called automatically by GetNoOfFitBins() and by
 *       PrepareData() after setting up the data arrays.
 *
 * \see GetNoOfFitBins(), SetFitRangeBin()
 */
void PRunSingleHisto::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 preprocessing pipeline for single histogram runs.
 *
 * Orchestrates the complete data loading and preprocessing workflow:
 * -# <b>Load raw data:</b> Fetch run from PRunDataHandler using run name
 * -# <b>Extract metadata:</b> Magnetic field, beam energy, temperature(s)
 * -# <b>Validate histograms:</b> Check that forward histogram numbers exist in data file
 * -# <b>Get time resolution:</b> Extract bin width (typically 0.1-10 ns)
 * -# <b>Determine t0:</b> Call GetProperT0() for muon arrival times
 * -# <b>Load histogram data:</b> Copy forward histogram bins from raw data
 * -# <b>Add runs (ADDRUN):</b> If multiple runs specified, add them with t0 alignment
 * -# <b>Group histograms:</b> Sum multiple detectors within a group (with t0 alignment)
 * -# <b>Get data range (fgb/lgb):</b> Call GetProperDataRange() for good bin limits
 * -# <b>Get fit range:</b> Call GetProperFitRange() for fit time window
 * -# <b>Check lifetime correction:</b> Determine if exponential decay should be removed (for viewing)
 * -# <b>Dispatch to preparation:</b>
 *    - kFit → PrepareFitData(): packing, background subtraction
 *    - kView (no lifetime corr.) → PrepareRawViewData(): packing, theory calculation
 *    - kView (with lifetime corr.) → PrepareViewData(): lifetime removal, theory
 *
 * <b>ADDRUN t0 Alignment:</b>
 * When adding runs, histograms are aligned by their t0 values:
 * \code
 * forward[k][j] += addRunData[k]->at(j + addT0[k] - mainT0[k])
 * \endcode
 * This ensures muon arrival times coincide across added runs.
 *
 * <b>Grouping t0 Alignment:</b>
 * When grouping histograms, they are aligned to the first histogram's t0:
 * \code
 * fForward[j] += forward[i][j + t0[i] - t0[0]]
 * \endcode
 *
 * \return true if all preprocessing steps succeeded, false otherwise
 *
 * \note If any step fails (missing data file, invalid histogram numbers, t0 errors),
 *       this method returns false and error messages are printed to std::cerr.
 *
 * \see GetProperT0(), GetProperDataRange(), GetProperFitRange(),
 *      PrepareFitData(), PrepareRawViewData(), PrepareViewData()
 */
Bool_t PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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 a dead time correction has to be done
  // this will be done automatically in the function itself, which also
  // checks in the global and run section
  DeadTimeCorrection(forward, histoNo);

  // check if there are runs to be added to the current one
  if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present
    PRawRunData *addRunData;
    std::vector<PDoubleVector> addForward;
    for (UInt_t i=1; i<fRunInfo->GetRunNameSize(); i++) { // loop over all ADDRUN's

      // 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 << ">> PRunSingleHisto::PrepareData(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!";
        std::cerr << std::endl;
        return false;
      }

      addForward.clear();
      addForward.resize(histoNo.size());   // resize to number of groups
      for (UInt_t j=0; j<histoNo.size(); j++) {
        addForward[j].resize(addRunData->GetDataBin(histoNo[j])->size());
        addForward[j] = *addRunData->GetDataBin(histoNo[j]);
      }
      DeadTimeCorrection(addForward, histoNo);

      // add forward run
      UInt_t addRunSize;
      for (UInt_t k=0; k<histoNo.size(); k++) { // fill each group
        addRunSize = addForward[k].size();
        for (UInt_t j=0; j<addRunSize; 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] += addForward[k][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)+fT0s[i]-fT0s[0] >= 0) && (j+fT0s[i]-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);

  // get the lifetimecorrection flag
  Bool_t lifetimecorrection = false;
  PMsrPlotList *plot = fMsrInfo->GetMsrPlotList();
  lifetimecorrection = plot->at(0).fLifeTimeCorrection;

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

  // cleanup
  histoNo.clear();

  return success;
}

//--------------------------------------------------------------------------
// PrepareFitData (protected)
//--------------------------------------------------------------------------
/**
 * \brief Prepares histogram data for fitting (kFit mode).
 *
 * Performs final data transformations after PrepareData() has loaded and grouped
 * the raw histogram data:
 * -# <b>Estimate N₀ (optional):</b> If MSR file requests it, call EstimateN0()
 * -# <b>Handle background:</b>
 *    - If background is fitted: leave data unchanged
 *    - If fixed background given: subtract it from all bins
 *    - If background range given: call EstimateBkg() and subtract estimate
 *    - If nothing specified: auto-estimate from bins [0.1×t0, 0.6×t0] with warning
 * -# <b>Packing (rebinning):</b> Combine consecutive bins to improve statistics:
 *    - If packing = 1: copy bins directly
 *    - If packing > 1: sum every 'packing' bins into one
 * -# <b>Normalization:</b> If fScaleN0AndBkg is true, divide by (packing × t_res × 1000)
 *    to normalize counts to 1/ns
 * -# <b>Error calculation:</b>
 *    - If N > 0: σ = √N (Poisson statistics)
 *    - If N = 0: σ = 1/normalizer (avoid division by zero in χ²)
 * -# <b>Set time grid:</b>
 *    - Data start time: (fgb - 0.5 + pack/2 - t0) × t_res
 *    - Data time step: pack × t_res
 * -# <b>Calculate fit bins:</b> Call CalcNoOfFitBins() to set fStartTimeBin, fEndTimeBin
 *
 * <b>Packing Algorithm:</b>
 * \code
 * for (i = fgb; i < lgb; i++) {
 *   value += forward[i];
 *   if ((i-fgb) % packing == 0 && i != fgb) {
 *     data.push_back(value / normalizer);
 *     error.push_back(sqrt(value) / normalizer);
 *     value = 0;
 *   }
 * }
 * \endcode
 *
 * <b>Background Handling Priority:</b>
 * -# Check if background is fitted (bkgFitParamNo ≠ -1) → leave data as-is
 * -# Check if fixed background given (bkgFix ≠ PMUSR_UNDEFINED) → subtract fixed value
 * -# Check if background range given (bkgRange[0] ≥ 0) → estimate and subtract
 * -# Fallback: auto-estimate from [0.1×t0, 0.6×t0] → print warning
 *
 * \param runData Pointer to raw run data handler (for metadata access)
 * \param histoNo Forward histogram number (for background estimation)
 *
 * \return true if preparation succeeded, false if EstimateBkg() failed
 *
 * \note This method populates fData (PRunData object) with packed data ready for fitting.
 *
 * \see PrepareData(), EstimateBkg(), EstimateN0(), CalcNoOfFitBins()
 */
Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoNo)
{
  if (fMsrInfo->EstimateN0()) {
    EstimateN0();
  }

  // transform raw histo data. This is done the following way (for details see the manual):
  // for the single histo fit, just the rebinned raw data are copied

  // check how the background shall be handled
  if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg shall **NOT** be fitted
    // 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 << ">> PRunSingleHisto::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;
      }
    } else { // fixed background given
      for (UInt_t i=0; i<fForward.size(); i++) {
        fForward[i] -= fRunInfo->GetBkgFix(0);
      }
    }
  }

  // everything looks fine, hence fill data set
  Int_t t0 = static_cast<Int_t>(fT0s[0]);
  Double_t value = 0.0;
  Double_t normalizer = 1.0;
  // in order that after rebinning the fit does not need to be redone (important for plots)
  // the value is normalize to per 1 nsec if scaling is whished
  if (fScaleN0AndBkg)
    normalizer = fPacking * (fTimeResolution * 1.0e3); // fTimeResolution us->ns
  // data start at data_start-t0
  // time shifted so that packing is included correctly, i.e. t0 == t0 after packing
  fData.SetDataTimeStart(fTimeResolution*((static_cast<Double_t>(fGoodBins[0])-0.5) + static_cast<Double_t>(fPacking)/2.0 - static_cast<Double_t>(t0)));
  fData.SetDataTimeStep(fTimeResolution*fPacking);
  for (Int_t i=fGoodBins[0]; i<fGoodBins[1]; i++) {
    if (fPacking == 1) {
      value = fForward[i];
      value /= normalizer;
      fData.AppendValue(value);
      if (value == 0.0)
        fData.AppendErrorValue(1.0/normalizer);
      else
        fData.AppendErrorValue(TMath::Sqrt(value));
    } else { // packed data, i.e. fPacking > 1
      if (((i-fGoodBins[0]) % fPacking == 0) && (i != fGoodBins[0])) { // fill data
        value /= normalizer;
        fData.AppendValue(value);
        if (value == 0.0)
          fData.AppendErrorValue(1.0/normalizer);
        else
          fData.AppendErrorValue(TMath::Sqrt(value));
        // reset values
        value = 0.0;
      }
      value += fForward[i];
    }
  }

  CalcNoOfFitBins();

  return true;
}

//--------------------------------------------------------------------------
// PrepareRawViewData (protected)
//--------------------------------------------------------------------------
/**
 * <p>Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing
 * without any life time correction.
 * <p>The following steps are preformed:
 * -# check if view packing is whished.
 * -# check that 'first good data bin', 'last good data bin', and 't0' makes any sense
 * -# packing (i.e. rebinnig)
 * -# calculate theory
 *
 * <b>return:</b>
 * - true, if everything went smooth
 * - false, otherwise.
 *
 * \param runData raw run data handler
 * \param histoNo forward histogram number
 */
Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo)
{
  // check if view_packing is wished
  Int_t packing = fPacking;
  if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) {
    packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking;
  }

  // calculate necessary norms
  Double_t dataNorm = 1.0, theoryNorm = 1.0;
  if (fScaleN0AndBkg) {
    dataNorm = 1.0/ (packing * (fTimeResolution * 1.0e3)); // fTimeResolution us->ns
  } else if (!fScaleN0AndBkg && (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0)) {
    theoryNorm = static_cast<Double_t>(fMsrInfo->GetMsrPlotList()->at(0).fViewPacking)/static_cast<Double_t>(fPacking);
  }

  // raw data, since PMusrCanvas is doing ranging etc.
  // start = the first bin which is a multiple of packing backward from first good data bin
  Int_t start = fGoodBins[0] - (fGoodBins[0]/packing)*packing;
  // end = last bin starting from start which is a multiple of packing and still within the data
  Int_t end   = start + ((fForward.size()-start)/packing)*packing;
  // 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) - ((static_cast<Int_t>(fT0s[0])+offset)/packing)*packing;
    end = start + ((fForward.size()-start)/packing)*packing;
    std::cerr << std::endl << ">> PRunSingleHisto::PrepareRawViewData(): **WARNING** data range was not provided, will try data range start = " << start << ".";
    std::cerr << std::endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE.";
    std::cerr << std::endl;
  }
  // check if start, end, and t0 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 << ">> PRunSingleHisto::PrepareRawViewData(): **ERROR** start data bin doesn't make any sense!";
    std::cerr << std::endl;
    return false;
  }
  // 3rd check if end is within proper bounds
  if ((end < 0) || (end > static_cast<Int_t>(fForward.size()))) {
    std::cerr << std::endl << ">> PRunSingleHisto::PrepareRawViewData(): **ERROR** end data bin doesn't make any sense!";
    std::cerr << std::endl;
    return false;
  }

  // everything looks fine, hence fill data set
  Int_t t0 = static_cast<Int_t>(fT0s[0]);
  Double_t value = 0.0;
  // data start time = (binStart - 0.5) + pack/2 - t0, with pack and binStart used as double
  fData.SetDataTimeStart(fTimeResolution*((static_cast<Double_t>(start)-0.5) + static_cast<Double_t>(packing)/2.0 - static_cast<Double_t>(t0)));
  fData.SetDataTimeStep(fTimeResolution*packing);

  for (Int_t i=start; i<end; i++) {
    if (((i-start) % packing == 0) && (i != start)) { // fill data
      value *= dataNorm;
      fData.AppendValue(value);
      if (value == 0.0)
        fData.AppendErrorValue(1.0);
      else
        fData.AppendErrorValue(TMath::Sqrt(value*dataNorm));
      // reset values
      value = 0.0;
    }
    value += fForward[i];
  }

  CalcNoOfFitBins();

  // fill theory vector for kView
  // 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 asymmetry
  Double_t N0;
  // check if norm is a parameter or a function
  if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
    N0 = par[fRunInfo->GetNormParamNo()-1];
  } else { // norm is a function
    // get function number
    UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
    // evaluate function
    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
  }
  N0 *= theoryNorm;

  // get tau
  Double_t tau;
  if (fRunInfo->GetLifetimeParamNo() != -1)
    tau = par[fRunInfo->GetLifetimeParamNo()-1];
  else
    tau = PMUON_LIFETIME;

  // get background
  Double_t bkg;
  if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted
    if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
      if (fRunInfo->GetBkgRange(0) >= 0) { // background range given
        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 << ">> PRunSingleHisto::PrepareRawViewData(): **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;
      }
      bkg = fBackground;
    } else { // fixed bkg given
      bkg = fRunInfo->GetBkgFix(0);
    }
  } else { // bkg fitted
    bkg = par[fRunInfo->GetBkgFitParamNo()-1];
  }
  bkg *= theoryNorm;

  // 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 = 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);
  }

  Double_t time;
  Double_t theoryValue;
  for (UInt_t i=0; i<size; i++) {
    time = fData.GetTheoryTimeStart() + i*fData.GetTheoryTimeStep();
    theoryValue = fTheory->Func(time, par, fFuncValues);
    if (fabs(theoryValue) > 1.0e10) {  // dirty hack needs to be fixed!!
      theoryValue = 0.0;
    }
    fData.AppendTheoryValue(N0*TMath::Exp(-time/tau)*(1.0+theoryValue)+bkg);
  }

  // clean up
  par.clear();

  return true;
}

//--------------------------------------------------------------------------
// PrepareViewData (protected)
//--------------------------------------------------------------------------
/**
 * <p>Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing
 * with life time correction, i.e. the exponential decay is removed.
 * <p>The following steps are preformed:
 * -# check if view packing is whished.
 * -# check that 'first good data bin', 'last good data bin', and 't0' makes any sense
 * -# transform data sets (see below).
 * -# calculate theory
 *
 * <p> Muon life time corrected data: Starting from
 * \f[ N(t) = N_0 e^{-t/\tau} [ 1 + A(t) ] + \mathrm{Bkg} \f]
 * it follows that
 * \f[ A(t) = (-1) + e^{+t/\tau}\, \frac{N(t)-\mathrm{Bkg}}{N_0}. \f]
 * For the error estimate only the statistical error of \f$ N(t) \f$ is used, and hence
 * \f[ \Delta A(t) = \frac{e^{t/\tau}}{N_0}\,\sqrt{\frac{N(t)}{p}} \f]
 * where \f$ p \f$ is the packing, and \f$ N(t) \f$ are the packed data, i.e.
 * \f[ N(t_i) = \frac{1}{p}\, \sum_{j=i}^{i+p} n_j \f]
 * with \f$ n_j \f$ the raw histogram data bins.
 *
 * <b>return:</b>
 * - true, if everything went smooth
 * - false, otherwise
 *
 * \param runData raw run data handler
 * \param histoNo forward histogram number
 */
Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histoNo)
{
  // check if view_packing is wished. This is a global option for all PLOT blocks!
  Int_t packing = fPacking;
  if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) {
    packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking;
  }
  // check if rrf_packing is present. This is a global option for all PLOT blocks, since operated on a single set of data.
  if (fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking > 0) {
    packing = fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking;
  }

  // calculate necessary norms
  Double_t dataNorm = 1.0, theoryNorm = 1.0;
  if (fScaleN0AndBkg) {
    dataNorm = 1.0/ (packing * (fTimeResolution * 1.0e3)); // fTimeResolution us->ns
  } else if (!fScaleN0AndBkg && (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0)) {
    theoryNorm = static_cast<Double_t>(fMsrInfo->GetMsrPlotList()->at(0).fViewPacking)/static_cast<Double_t>(fPacking);
  }

  // transform raw histo data. This is done the following way (for details see the manual):
  // for the single histo fit, just the rebinned raw data are copied
  // first get start data, end data, and t0
  Int_t t0 = static_cast<Int_t>(fT0s[0]);

  // start = the first bin which is a multiple of packing backward from first good data bin
  Int_t start = fGoodBins[0] - (fGoodBins[0]/packing)*packing;
  // end = last bin starting from start which is a multiple of packing and still within the data
  Int_t end   = start + ((fForward.size()-start)/packing)*packing;

  // 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) - ((static_cast<Int_t>(fT0s[0])+offset)/packing)*packing;
    end = start + ((fForward.size()-start)/packing)*packing;
    std::cerr << std::endl << ">> PRunSingleHisto::PrepareViewData(): **WARNING** data range was not provided, will try data range start = " << start << ".";
    std::cerr << std::endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE.";
    std::cerr << std::endl;
  }

  // check if start, end, and t0 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 << ">> PRunSingleHisto::PrepareViewData(): **ERROR** start data bin doesn't make any sense!";
    std::cerr << std::endl;
    return false;
  }
  // 3rd check if end is within proper bounds
  if ((end < 0) || (end > static_cast<Int_t>(fForward.size()))) {
    std::cerr << std::endl << ">> PRunSingleHisto::PrepareViewData(): **ERROR** end data bin doesn't make any sense!";
    std::cerr << std::endl;
    return false;
  }

  // everything looks fine, hence fill data set

  // 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 asymmetry
  Double_t N0;
  // check if norm is a parameter or a function
  if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter
    N0 = par[fRunInfo->GetNormParamNo()-1];
  } else { // norm is a function
    // get function number
    UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
    // evaluate function
    N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
  }
  N0 *= theoryNorm;

  // get tau
  Double_t tau;
  if (fRunInfo->GetLifetimeParamNo() != -1)
    tau = par[fRunInfo->GetLifetimeParamNo()-1];
  else
    tau = PMUON_LIFETIME;

  // get background
  Double_t bkg;
  if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted
    if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval)
      if (fRunInfo->GetBkgRange(0) >= 0) { // background range given
        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 << ">> PRunSingleHisto::PrepareViewData(): **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;
      }
      bkg = fBackground;
    } else { // fixed bkg given
      bkg = fRunInfo->GetBkgFix(0);
    }
  } else { // bkg fitted
    bkg = par[fRunInfo->GetBkgFitParamNo()-1];
  }
  bkg *= theoryNorm;

  Double_t value = 0.0;
  Double_t expval = 0.0;
  Double_t rrf_val = 0.0;
  Double_t time = 0.0;

  // data start time = (binStart - 0.5) + pack/2 - t0, with pack and binStart used as double
  fData.SetDataTimeStart(fTimeResolution*((static_cast<Double_t>(start)-0.5) + static_cast<Double_t>(packing)/2.0 - static_cast<Double_t>(t0)));
  fData.SetDataTimeStep(fTimeResolution*packing);

  // data is always normalized to (per nsec!!)
  Double_t gammaRRF = 0.0, wRRF = 0.0, phaseRRF = 0.0;
  if (fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq == 0.0) { // normal Data representation
    for (Int_t i=start; i<end; i++) {
      if (((i-start) % packing == 0) && (i != start)) { // fill data
        value *= dataNorm;
        // since the packing counter is already at the end of the bin, the time needs be shifted back by pack*time_resolution
        time = (((static_cast<Double_t>(i)-0.5) + static_cast<Double_t>(packing)/2.0 - static_cast<Double_t>(t0)))*fTimeResolution - static_cast<Double_t>(packing)*fTimeResolution;
        expval = TMath::Exp(+time/tau)/N0;
        fData.AppendValue(-1.0+expval*(value-bkg));
        fData.AppendErrorValue(expval*TMath::Sqrt(value*dataNorm));
        value = 0.0;
      }
      value += fForward[i];
    }
  } else { // RRF representation
    // check which units shall be used
    switch (fMsrInfo->GetMsrPlotList()->at(0).fRRFUnit) {
      case RRF_UNIT_kHz:
        gammaRRF = TMath::TwoPi()*1.0e-3;
        break;
      case RRF_UNIT_MHz:
        gammaRRF = TMath::TwoPi();
        break;
      case RRF_UNIT_Mcs:
        gammaRRF = 1.0;
        break;
      case RRF_UNIT_G:
        gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi();
        break;
      case RRF_UNIT_T:
        gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi()*1.0e4;
        break;
      default:
        gammaRRF = TMath::TwoPi();
        break;
    }
    wRRF = gammaRRF * fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq;
    phaseRRF = fMsrInfo->GetMsrPlotList()->at(0).fRRFPhase / 180.0 * TMath::Pi();

    Double_t error = 0.0;
    for (Int_t i=start; i<end; i++) {
      if (((i-start) % packing == 0) && (i != start)) { // fill data
        fData.AppendValue(2.0*value/packing); // factor 2 needed because cos(a)cos(b) = 1/2(cos(a+b)+cos(a-b))
        fData.AppendErrorValue(expval*TMath::Sqrt(error/packing));
        value = 0.0;
        error = 0.0;
      }
      time = (static_cast<Double_t>(i)-t0)*fTimeResolution;
      expval = TMath::Exp(+time/tau)/N0;
      rrf_val = (-1.0+expval*(fForward[i]/(fTimeResolution*1.0e3)-bkg))*TMath::Cos(wRRF * time + phaseRRF);
      value += rrf_val;
      error += fForward[i]*dataNorm;
    }
  }

  CalcNoOfFitBins();

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

  // calculate theory
  Double_t theoryValue;
  UInt_t size = fForward.size()/packing;
  const Int_t factor = 8; // 8 times more points for the theory (if fTheoAsData == false)
  UInt_t rebinRRF = 0;

  if (wRRF == 0) { // no RRF
    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);
    }
  } else { // RRF
    rebinRRF = static_cast<UInt_t>((TMath::Pi()/2.0/wRRF)/fTimeResolution); // RRF time resolution / data time resolution
    fData.SetTheoryTimeStart(fData.GetDataTimeStart());
    fData.SetTheoryTimeStep(TMath::Pi()/2.0/wRRF/rebinRRF); // = theory time resolution as close as possible to the data time resolution compatible with wRRF
  }

  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 (wRRF != 0.0) {
      theoryValue *= 2.0*TMath::Cos(wRRF * time + phaseRRF);
    }
    if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!!
      theoryValue = 0.0;
    }
    fData.AppendTheoryValue(theoryValue);
  }

  // if RRF filter the theory with a FIR Kaiser low pass filter
  if (wRRF != 0.0) {
    // rebin theory to the RRF frequency
    if (rebinRRF != 0) {
      Double_t dval = 0.0;
      PDoubleVector theo;
      for (UInt_t i=0; i<fData.GetTheory()->size(); i++) {
        if ((i % rebinRRF == 0) && (i != 0)) {
          theo.push_back(dval/rebinRRF);
          dval = 0.0;
        }
        dval += fData.GetTheory()->at(i);
      }
      fData.SetTheoryTimeStart(fData.GetTheoryTimeStart()+static_cast<Double_t>(rebinRRF-1)*fData.GetTheoryTimeStep()/2.0);
      fData.SetTheoryTimeStep(rebinRRF*fData.GetTheoryTimeStep());
      fData.ReplaceTheory(theo);
      theo.clear();
    }

    // filter theory
    CalculateKaiserFilterCoeff(wRRF, 60.0, 0.2); // w_c = wRRF, A = -20 log_10(delta), Delta w / w_c = (w_s - w_p) / (2 w_c)
    FilterTheo();
  }

  // clean up
  par.clear();

  return true;
}

//--------------------------------------------------------------------------
// GetProperT0 (private)
//--------------------------------------------------------------------------
/**
 * \brief Determines time-zero (t0) values for all histograms using hierarchical fallback.
 *
 * Time-zero (t0) marks the muon arrival time in each detector histogram, the reference
 * point from which decay time is measured. This method uses a priority system to find
 * t0 values:
 *
 * <b>Priority hierarchy (highest to lowest):</b>
 * -# <b>RUN block t0:</b> Explicitly specified in the RUN block (highest priority)
 * -# <b>GLOBAL block t0:</b> Default t0 for all runs in the GLOBAL block
 * -# <b>Data file t0:</b> Stored in the raw data file (from previous analysis)
 * -# <b>Estimated t0:</b> Automatic estimation (UNRELIABLE, prints warning)
 *
 * <b>For ADDRUN support:</b>
 * If multiple runs are added (fRunInfo->GetRunNameSize() > 1), this method also
 * determines t0 values for each added run (fAddT0s) using the same hierarchy.
 * Proper t0 alignment is essential for correct ADDRUN operation.
 *
 * <b>Algorithm:</b>
 * -# Resize fT0s vector to histogram count (number of grouped detectors)
 * -# Initialize all t0 values to -1.0 (sentinel for "not set")
 * -# Fill from RUN block (if specified)
 * -# Fill from GLOBAL block where still -1.0
 * -# Fill from data file where still -1.0
 * -# Fill from estimation where still -1.0 (prints **WARNING**)
 * -# Validate all t0 values are within histogram bounds
 * -# If ADDRUN present: repeat steps 2-6 for each added run
 *
 * <b>Validation:</b>
 * After fallback, checks that each t0 satisfies:
 * \f[
 * 0 \leq t_0 \leq N_{\rm bins}
 * \f]
 * If validation fails, returns false with error message.
 *
 * \param runData Pointer to raw run data handler for main run
 * \param globalBlock Pointer to GLOBAL block from MSR file
 * \param histoNo Vector of histogram indices (zero-based, after redGreen offset correction)
 *
 * \return true if all t0 values found and validated, false if any t0 is out of bounds
 *
 * \warning Estimated t0 values (fallback option #4) are often UNRELIABLE, especially
 *          for low-energy muons (LEM). Manual specification in MSR file is strongly
 *          recommended. A warning is printed to std::cerr when estimation is used.
 *
 * \note This method updates fT0s (main run) and fAddT0s (ADDRUN) member variables.
 *       It also updates the MSR file handler with found t0 values for persistence.
 *
 * \see PrepareData(), fT0s, fAddT0s
 */
Bool_t PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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.0) || (fAddT0s[i-1][j] > static_cast<Int_t>(addRunData->GetDataBin(histoNo[j])->size()))) {
          std::cerr << std::endl << ">> PRunSingleHisto::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 the data range (first good bin / last good bin).
 *
 * Establishes which histogram bins contain valid muon decay data by
 * finding the "first good bin" (fgb) and "last good bin" (lgb). This
 * range excludes:
 * - Pre-t0 bins (before muon arrival)
 * - Early bins affected by detector dead time or pileup
 * - Late bins with insufficient statistics
 *
 * <b>Priority hierarchy (highest to lowest):</b>
 * -# <b>RUN block:</b> Explicitly specified fgb/lgb in RUN block
 * -# <b>GLOBAL block:</b> Default fgb/lgb from GLOBAL block
 * -# <b>Auto-estimation:</b> Fallback estimates with warning
 *
 * <b>Auto-estimation (if not specified):</b>
 * - <b>fgb:</b> t0 + 10 ns (to avoid dead time issues)
 * - <b>lgb:</b> End of histogram (all bins)
 *
 * <b>Validation:</b>
 * -# Check fgb < lgb (swap if necessary)
 * -# Check 0 ≤ fgb < histogram length
 * -# Check 0 ≤ lgb ≤ histogram length
 * -# If lgb > histogram length: clamp to (length - 1) and print warning
 *
 * <b>Storage:</b>
 * Results are stored in:
 * - fGoodBins[0] = fgb (first good bin index)
 * - fGoodBins[1] = lgb (last good bin index)
 *
 * These values are used by:
 * - PrepareFitData() to determine packing range
 * - GetProperFitRange() as fallback for fit range
 *
 * \return true if data range is valid and within bounds, false if validation fails
 *
 * \warning Auto-estimated ranges may not be appropriate for all detectors.
 *          Explicit specification in MSR file is strongly recommended.
 *
 * \note This method is called by PrepareData() after histogram grouping
 *       but before packing and fit range determination.
 *
 * \see PrepareData(), GetProperFitRange(), fGoodBins
 */
Bool_t PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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 the fit range (start and end times for χ² calculation).
 *
 * Establishes the time window [t_start, t_end] over which the fit will be
 * performed. The fit range can be specified in two ways:
 *
 * <b>Specification methods:</b>
 * -# <b>Time-based:</b> `fit <start> <end>` in microseconds
 *    - Example: `fit 0.1 10.0` (fit from 0.1 to 10.0 μs after t0)
 * -# <b>Bin-based:</b> `fit fgb[+offset0] lgb[-offset1]` in bins
 *    - Example: `fit fgb+10 lgb-20` (fit from 10 bins after fgb to 20 bins before lgb)
 *
 * <b>Priority hierarchy (highest to lowest):</b>
 * -# <b>RUN block time-based:</b> `fit <start> <end>` in RUN block
 * -# <b>RUN block bin-based:</b> `fit fgb+n0 lgb-n1` in RUN block
 * -# <b>GLOBAL block time-based:</b> `fit <start> <end>` in GLOBAL block
 * -# <b>GLOBAL block bin-based:</b> `fit fgb+n0 lgb-n1` in GLOBAL block
 * -# <b>Auto-fallback:</b> Use entire data range [fgb, lgb]
 *
 * <b>Bin-based conversion to time:</b>
 * When fit range is given in bins, it's converted to time:
 * \f[
 * t_{\rm start} = (\text{fgb} + n_0 - t_0) \times \Delta t
 * \f]
 * \f[
 * t_{\rm end} = (\text{lgb} - n_1 - t_0) \times \Delta t
 * \f]
 *
 * where:
 * - fgb/lgb = first/last good bin from GetProperDataRange()
 * - n₀/n₁ = offsets (can be positive or negative)
 * - t₀ = time-zero bin
 * - Δt = time resolution (fTimeResolution in μs)
 *
 * <b>Storage and updates:</b>
 * - fFitStartTime, fFitEndTime are set to the determined range
 * - If bin-based, the converted time values are written back to the MSR
 *   data structure for log file reporting
 *
 * <b>Fallback behavior:</b>
 * If no fit range is specified anywhere, uses the entire data range:
 * \f[
 * t_{\rm start} = (\text{fgb} - t_0) \times \Delta t
 * \f]
 * \f[
 * t_{\rm end} = (\text{lgb} - t_0) \times \Delta t
 * \f]
 * and prints a warning to std::cerr.
 *
 * \param globalBlock Pointer to GLOBAL block from MSR file
 *
 * \note This method is called by PrepareData() after GetProperDataRange()
 *       has established fGoodBins[0] and fGoodBins[1].
 *
 * \see PrepareData(), GetProperDataRange(), SetFitRangeBin(), CalcNoOfFitBins()
 */
void PRunSingleHisto::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 << ">> PRunSingleHisto::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;
  }
}

//--------------------------------------------------------------------------
// EstimateN0 (private)
//--------------------------------------------------------------------------
/**
 * \brief Automatically estimates the normalization parameter N₀ from data.
 *
 * Provides an intelligent initial guess for N₀ to help MINUIT convergence.
 * The estimate is based on the maximum count rate in the fit range, accounting
 * for muon decay and background.
 *
 * <b>When estimation is performed:</b>
 * - MSR file requests estimation (estimate_n0 flag in GLOBAL block)
 * - Norm parameter is a fit parameter (not fixed, not a function)
 * - Parameter step size ≠ 0 (i.e., not fixed)
 *
 * <b>When estimation is skipped:</b>
 * - Norm is a function (paramNo > MSR_PARAM_FUN_OFFSET)
 * - Norm parameter is fixed (step = 0)
 * - Invalid parameter number
 *
 * <b>Estimation algorithm:</b>
 * -# Find maximum value in fit range: max_data = max(N(t) in fit range)
 * -# Find corresponding time t_max
 * -# Extract or estimate background B
 * -# Correct for exponential decay: N₀_est = (max_data - B) / exp(-t_max/τ_μ)
 * -# Adjust for scaling if fScaleN0AndBkg is true
 * -# Update parameter value and step size in MSR parameter list
 *
 * <b>Background handling:</b>
 * - If background is fitted: extract current background parameter value
 * - If fixed background given: use fixed value
 * - If background range given: use fBackground estimate
 * - Otherwise: assume B = 0
 *
 * <b>Scaling adjustment:</b>
 * If fScaleN0AndBkg is true (normalizing to 1/ns), the estimate is divided by:
 * \f[
 * \text{scale factor} = \text{packing} \times (t_{\rm res} \times 1000)
 * \f]
 *
 * \note This method modifies the MSR parameter list in place, updating both
 *       the parameter value and the step size (for MINUIT error estimation).
 *
 * \see IsScaleN0AndBkg(), EstimateBkg(), PrepareFitData()
 */
void PRunSingleHisto::EstimateN0()
{
  // check that 'norm' in the msr-file run block is indeed a parameter number.
  // in case it is a function, nothing will be done.
  UInt_t paramNo = fRunInfo->GetNormParamNo();
  if (paramNo > 10000) // i.e. fun or map
    return;

  // get the parameters
  PMsrParamList *param = fMsrInfo->GetMsrParamList();
  assert(param);

  if (paramNo > param->size()) {
    std::cerr << std::endl << ">> PRunSingleHisto::EstimateN0: **ERROR** found parameter number " << paramNo << ", which is larger than the number of parameters = " << param->size() << std::endl;
    return;
  }

  // check if N0 is fixed. If this is the case, do NOT estimate N0
  if (param->at(paramNo-1).fStep == 0.0) // N0 parameter fixed
    return;


  // check that 'backgr.fit' in the msr-file run block is indeed a parameter number.
  // in case it is a function, nothing will be done.
  Int_t paramNoBkg = fRunInfo->GetBkgFitParamNo();
  Bool_t scaleBkg = true;
  Double_t bkg=0.0, errBkg=1.0;
  if ((paramNoBkg > 10000) || (paramNoBkg == -1)) { // i.e. fun or map
    scaleBkg = false;
  } else {
    if (paramNoBkg-1 < static_cast<Int_t>(param->size())) {
      bkg = param->at(paramNoBkg-1).fValue;
      errBkg = param->at(paramNoBkg-1).fStep;
    }
  }

  // estimate N0
  Double_t dt = fTimeResolution;
  Double_t tau = PMUON_LIFETIME;

  UInt_t t0 = static_cast<UInt_t>(round(fT0s[0]));
  Double_t dval = 0.0;
  Double_t nom = 0.0;
  Double_t denom = 0.0;
  Double_t xx = 0.0;

  // calc nominator
  for (UInt_t i=t0; i<fForward.size(); i++) {
    xx = exp(-dt*static_cast<Double_t>(i-t0)/tau);
    nom += xx;
  }

  // calc denominator
  for (UInt_t i=t0; i<fForward.size(); i++) {
    xx = exp(-dt*static_cast<Double_t>(i-t0)/tau);
    dval = fForward[i];
    if (dval > 0)
      denom += xx*xx/dval;
  }
  Double_t N0 = nom/denom;  

  if (fScaleN0AndBkg) {
    N0 /= fTimeResolution*1.0e3;
  } else {
    N0 *= fPacking;
  }

  Double_t rescale = 1;
  if ((param->at(paramNo-1).fValue != 0.0) && scaleBkg) {
    rescale = N0 / param->at(paramNo-1).fValue;
    bkg *= rescale;
    errBkg *= rescale;
  }

  std::cout << ">> PRunSingleHisto::EstimateN0: found N0=" << param->at(paramNo-1).fValue << ", will set it to N0=" << N0 << std::endl;
  if (scaleBkg)
    std::cout << ">> PRunSingleHisto::EstimateN0: found Bkg=" << param->at(paramNoBkg-1).fValue << ", will set it to Bkg=" << bkg << std::endl;
  fMsrInfo->SetMsrParamValue(paramNo-1, N0);
  fMsrInfo->SetMsrParamStep(paramNo-1, sqrt(fabs(N0)));
  if (scaleBkg) {
    fMsrInfo->SetMsrParamValue(paramNoBkg-1, bkg);
    fMsrInfo->SetMsrParamStep(paramNoBkg-1, errBkg);
  }
}

//--------------------------------------------------------------------------
// EstimateBkg (private)
//--------------------------------------------------------------------------
/**
 * \brief Estimates background count rate from pre-t0 bins.
 *
 * Calculates the average background rate from bins before the muon pulse
 * arrives. For pulsed beam facilities (PSI, RAL, TRIUMF), adjusts the
 * background interval to be a multiple of the beam period to avoid
 * systematic biases from beam structure.
 *
 * <b>Algorithm:</b>
 * -# Extract background range [start, end] from MSR file (in bins)
 * -# Validate start < end (swap if necessary)
 * -# If pulsed beam (PSI/RAL/TRIUMF):
 *    - Calculate interval duration in time: t_bkg = (end - start) × t_res × packing
 *    - Find number of complete beam cycles: N_cycles = floor(t_bkg / T_beam)
 *    - Adjust end bin to match N_cycles × T_beam exactly
 * -# Validate start and end are within histogram bounds
 * -# Sum counts in [start, end]: Σ fForward[i]
 * -# Calculate average: fBackground = Σ counts / (end - start)
 *
 * <b>Beam periods:</b>
 * - <b>PSI:</b> 19.75 ns (50.63 MHz cyclotron)
 * - <b>RAL (ISIS):</b> 320 ns (3.125 MHz target)
 * - <b>TRIUMF:</b> 43.0 ns (23.26 MHz cyclotron)
 * - <b>Other facilities:</b> No period correction applied
 *
 * <b>Why adjust to beam period?</b>
 * Pulsed beams have time-dependent backgrounds from:
 * - Flash (instantaneous background from beam pulse)
 * - Prompt particles
 * - Pion background
 *
 * Averaging over complete beam cycles ensures unbiased background estimates
 * by including all phases of the pulsed structure.
 *
 * <b>Edge cases:</b>
 * - If interval < 1 beam period: uses original end bin (no correction)
 * - If start ≥ histogram length: returns false with error
 * - If end ≥ histogram length: returns false with error
 *
 * \param histoNo Forward histogram number (for error messages, currently not directly used)
 *
 * \return true if background estimated successfully, false if bins out of bounds
 *
 * \note The estimated background is stored in fBackground member variable
 *       and subtracted from data in PrepareFitData() if not fitted.
 *
 * \see PrepareFitData(), ACCEL_PERIOD_PSI, ACCEL_PERIOD_RAL, ACCEL_PERIOD_TRIUMF
 */
Bool_t PRunSingleHisto::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
  Int_t start = fRunInfo->GetBkgRange(0);
  Int_t end   = fRunInfo->GetBkgRange(1);
  if (end < start) {
    std::cout << std::endl << "PRunSingleHisto::EstimatBkg(): end = " << end << " > start = " << start << "! Will swap them!";
    Int_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*fPacking); // 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*fPacking));
    std::cout << std::endl << "PRunSingleHisto::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 << ">> PRunSingleHisto::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 << ">> PRunSingleHisto::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);

  if (fScaleN0AndBkg)
    fBackground = bkg / (fTimeResolution * 1e3); // keep background (per 1 nsec) for chisq, max.log.likelihood, fTimeResolution us->ns
  else
    fBackground = bkg * fPacking;  // keep background (per bin)

  fRunInfo->SetBkgEstimated(fBackground, 0);

  return true;
}

//--------------------------------------------------------------------------
// IsScaleN0AndBkg (private)
//--------------------------------------------------------------------------
/**
 * \brief Determines if N₀ and background should be normalized to 1/ns.
 *
 * Checks whether N₀ and background parameters should be scaled to represent
 * count rates per nanosecond (1/ns) rather than counts per packed bin.
 *
 * <b>Default behavior:</b> Scaling is ENABLED (true)
 *
 * This makes fitted parameters physically meaningful and independent of packing:
 * - N₀ represents the initial muon decay rate at t=0 in counts/ns
 * - Background B represents constant background rate in counts/ns
 *
 * <b>To disable scaling:</b> Add to MSR file COMMAND block:
 * \code
 * SCALE_N0_BKG FALSE
 * \endcode
 *
 * <b>When to disable scaling:</b>
 * - When N₀ and B should represent total counts per packed bin
 * - When comparing with older analysis that didn't use scaling
 * - When packing is 1 (no difference between modes)
 *
 * <b>Effect on fit parameters:</b>
 * - <b>Scaled (default):</b> N₀ and B independent of packing choice
 * - <b>Unscaled:</b> N₀ and B depend on packing value
 *
 * <b>Implementation details:</b>
 * Scaling is applied in:
 * - PrepareFitData(): Data is divided by (packing × t_res × 1000)
 * - CalcChiSquare(): χ² is multiplied by (packing × t_res × 1000)
 * - CalcMaxLikelihood(): -2ln(L) is multiplied by normalizer
 * - EstimateBkg(): Background estimate is divided by (t_res × 1000)
 *
 * These operations cancel out mathematically but keep parameters in 1/ns units.
 *
 * \return true if N₀ and background should be scaled to 1/ns (default),
 *         false if they should represent counts per packed bin
 *
 * \note This method is called during construction to set fScaleN0AndBkg.
 *
 * \see CalcChiSquare(), CalcMaxLikelihood(), PrepareFitData(), EstimateBkg()
 */
Bool_t PRunSingleHisto::IsScaleN0AndBkg()
{
  Bool_t willScale = true;

  PMsrLines *cmd = fMsrInfo->GetMsrCommands();
  for (UInt_t i=0; i<cmd->size(); i++) {
    if (cmd->at(i).fLine.Contains("SCALE_N0_BKG", TString::kIgnoreCase)) {
      TObjArray *tokens = nullptr;
      TObjString *ostr = nullptr;
      TString str;
      tokens = cmd->at(i).fLine.Tokenize(" \t");
      if (tokens->GetEntries() != 2) {
        std::cerr << std::endl << ">> PRunSingleHisto::IsScaleN0AndBkg(): **WARNING** Found uncorrect 'SCALE_N0_BKG' command, will ignore it.";
        std::cerr << std::endl << ">> Allowed commands: SCALE_N0_BKG TRUE | FALSE" << std::endl;
        return willScale;
      }
      ostr = dynamic_cast<TObjString*>(tokens->At(1));
      str = ostr->GetString();
      if (!str.CompareTo("FALSE", TString::kIgnoreCase)) {
        willScale = false;
      }
      // clean up
      if (tokens)
        delete tokens;
    }
  }

  return willScale;
}