improve the doxygen docu of PRunSingleHisto.*

src/classes/PRunSingleHisto.cpp

@@ -50,7 +50,17 @@
 // Constructor
 //--------------------------------------------------------------------------
 /**
- * <p>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()
 {
@@ -73,12 +83,27 @@ PRunSingleHisto::PRunSingleHisto() : PRunBase()
 // Constructor
 //--------------------------------------------------------------------------
 /**
- * <p>Constructor
+ * \brief Main constructor for single histogram fitting and viewing.
  *
- * \param msrInfo pointer to the msr-file handler
- * \param rawData raw run data
- * \param runNo number of the run within the msr-file
- * \param tag tag showing what shall be done: kFit == fitting, kView == 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)
@@ -119,7 +144,11 @@ PRunSingleHisto::PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData,
 // Destructor
 //--------------------------------------------------------------------------
 /**
- * <p>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()
 {
@@ -130,12 +159,51 @@ PRunSingleHisto::~PRunSingleHisto()
 // CalcChiSquare (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Calculate chi-square.
+ * \brief Calculates χ² between data and theory (least-squares fit metric).
  *
- * <b>return:</b>
- * - chisq value
+ * 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]
  *
- * \param par parameter vector iterated by minuit2
+ * 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)
 {
@@ -214,12 +282,43 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
 // CalcChiSquareExpected (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Calculate expected chi-square.
+ * \brief Calculates expected χ² using theory as variance (alternative fit metric).
  *
- * <b>return:</b>
- * - chisq value
+ * 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]
  *
- * \param par parameter vector iterated by minuit2
+ * <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)
 {
@@ -299,12 +398,59 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& par
 // CalcMaxLikelihood (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Calculate log maximum-likelihood. See http://pdg.lbl.gov/index.html
+ * \brief Calculates -2 log(maximum likelihood) for Poisson-distributed histogram data.
  *
- * <b>return:</b>
- * - log maximum-likelihood value
+ * 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]
  *
- * \param par parameter vector iterated by minuit2
+ * 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 (note: reduction(-:mllh) for subtraction)
+ *
+ * <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)
 {
@@ -397,12 +543,49 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
 // CalcMaxLikelihoodExpected (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Calculate expected log maximum-likelihood.
+ * \brief Calculates expected -2 log(maximum likelihood) using G-test formulation.
  *
- * <b>return:</b>
- * - log maximum-likelihood value
+ * 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]
  *
- * \param par parameter vector iterated by minuit2
+ * <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)
 {
@@ -493,7 +676,40 @@ Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector<Double_t>&
 // CalcTheory (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Calculate theory for a given set of fit-parameters.
+ * \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()
 {
@@ -557,9 +773,21 @@ void PRunSingleHisto::CalcTheory()
 // GetNoOfFitBins (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Calculate the number of fitted bins for the current fit range.
+ * \brief Returns the number of bins in the current fit range.
  *
- * <b>return:</b> number of fitted bins.
+ * 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()
 {
@@ -572,15 +800,47 @@ UInt_t PRunSingleHisto::GetNoOfFitBins()
 // SetFitRangeBin (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Allows to change the fit range on the fly. Used in the COMMAND block.
- * The syntax of the string is: FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]].
- * If only one pair of fgb/lgb is given, it is used for all runs in the RUN block section.
- * If multiple fgb/lgb's are given, the number N has to be the number of RUN blocks in
- * the msr-file.
+ * \brief Dynamically changes the fit range from COMMAND block instructions.
  *
- * <p>nXY are offsets which can be used to shift, limit the fit range.
+ * Parses and applies a FIT_RANGE command to modify the fit range on the fly,
+ * typically used during interactive fitting sessions or systematic scans.
  *
- * \param fitRange string containing the necessary information.
+ * <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)
 {
@@ -662,7 +922,31 @@ void PRunSingleHisto::SetFitRangeBin(const TString fitRange)
 // CalcNoOfFitBins (public)
 //--------------------------------------------------------------------------
 /**
- * <p>Calculate the number of fitted bins for the current fit range.
+ * \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()
 {
@@ -684,17 +968,45 @@ void PRunSingleHisto::CalcNoOfFitBins()
 // PrepareData (protected)
 //--------------------------------------------------------------------------
 /**
- * <p>Prepare data for fitting or viewing. What is already processed at this stage:
- * -# get proper raw run data
- * -# get all needed forward histograms
- * -# get time resolution
- * -# get t0's and perform necessary cross checks (e.g. if t0 of msr-file (if present) are consistent with t0 of the data files, etc.)
- * -# add runs (if addruns are present)
- * -# group histograms (if grouping is present)
+ * \brief Main data preprocessing pipeline for single histogram runs.
  *
- * <b>return:</b>
- * - true if everthing went smooth
- * - false, otherwise.
+ * 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()
 {
@@ -834,19 +1146,55 @@ Bool_t PRunSingleHisto::PrepareData()
 // PrepareFitData (protected)
 //--------------------------------------------------------------------------
 /**
- * <p>Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for fitting.
- * The following steps are preformed:
- * -# get fit start/stop time
- * -# check that 'first good data bin', 'last good data bin', and 't0' make any sense
- * -# check how the background shall be handled, i.e. fitted, subtracted from background estimate data range, or subtacted from a given fixed background.
- * -# packing (i.e rebinning)
+ * \brief Prepares histogram data for fitting (kFit mode).
  *
- * <b>return:</b>
- * - true, if everything went smooth
- * - false, otherwise
+ * 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
  *
- * \param runData raw run data handler
- * \param histoNo forward histogram number
+ * <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)
 {
@@ -1383,21 +1731,54 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo
 // GetProperT0 (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Get the proper t0 for the single histogram run.
- * -# the t0 vector size = number of detectors (grouping) for forward.
- * -# initialize t0's with -1
- * -# fill t0's from RUN block
- * -# if t0's are missing (i.e. t0 == -1), try to fill from the GLOBAL block.
- * -# if t0's are missing, try t0's from the data file
- * -# if t0's are missing, try to estimate them
+ * \brief Determines time-zero (t0) values for all histograms using hierarchical fallback.
  *
- * \param runData pointer to the current RUN block entry from the msr-file
- * \param globalBlock pointer to the GLOBLA block entry from the msr-file
- * \param histoNo histogram number vector of forward; histoNo = msr-file forward + redGreen_offset - 1
+ * 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>return:</b>
- * - true if everthing went smooth
- * - false, otherwise.
+ * <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)
 {
@@ -1521,14 +1902,48 @@ Bool_t PRunSingleHisto::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globa
 // GetProperDataRange (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Get the proper data range, i.e. first/last good bin (fgb/lgb).
- * -# get fgb/lgb from the RUN block
- * -# if fgb/lgb still undefined, try to get it from the GLOBAL block
- * -# if fgb/lgb still undefined, try to estimate them.
+ * \brief Determines the data range (first good bin / last good bin).
  *
- * <b>return:</b>
- * - true if everthing went smooth
- * - false, otherwise.
+ * 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()
 {
@@ -1609,16 +2024,60 @@ Bool_t PRunSingleHisto::GetProperDataRange()
 // GetProperFitRange (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Get the proper fit range. There are two possible fit range commands:
- * fit <start> <end> given in (usec), or
- * fit fgb+offset_0 lgb-offset_1 given in (bins), therefore it works the following way:
- * -# get fit range assuming given in time from RUN block
- * -# if fit range in RUN block is given in bins, replace start/end
- * -# if fit range is NOT given yet, try fit range assuming given in time from GLOBAL block
- * -# if fit range in GLOBAL block is given in bins, replace start/end
- * -# if still no fit range is given, use fgb/lgb.
+ * \brief Determines the fit range (start and end times for χ² calculation).
  *
- * \param globalBlock pointer to the GLOBAL block information form the msr-file.
+ * 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)
 {
@@ -1657,7 +2116,46 @@ void PRunSingleHisto::GetProperFitRange(PMsrGlobalBlock *globalBlock)
 // EstimateN0 (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Estimate the N0 for the given run.
+ * \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()
 {
@@ -1745,16 +2243,55 @@ void PRunSingleHisto::EstimateN0()
 }
 
 //--------------------------------------------------------------------------
-// EstimatBkg (private)
+// EstimateBkg (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Estimate the background for a given interval.
+ * \brief Estimates background count rate from pre-t0 bins.
  *
- * <b>return:</b>
- * - true, if everything went smooth
- * - false, otherwise
+ * 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.
  *
- * \param histoNo forward histogram number of the run
+ * <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
+ * - Slow neutron 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)
 {
@@ -1831,14 +2368,46 @@ Bool_t PRunSingleHisto::EstimateBkg(UInt_t histoNo)
 // IsScaleN0AndBkg (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Checks if N0/Bkg normalization to 1/ns is whished. The default is yes, since most of the users want to have it that way.
- * To overwrite this, one should add the line 'SCALE_N0_BKG FALSE' to the command block of the msr-file.
+ * \brief Determines if N₀ and background should be normalized to 1/ns.
  *
- * <b>return:</b>
- * - true, if scaling of N0 and Bkg to 1/ns is whished
- * - false, otherwise
+ * Checks whether N₀ and background parameters should be scaled to represent
+ * count rates per nanosecond (1/ns) rather than counts per packed bin.
  *
- * \param histoNo forward histogram number of the run
+ * <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()
 {

src/include/PRunSingleHisto.h

@@ -33,59 +33,411 @@
 #include "PRunBase.h"
 
 /**
- * <p>Class handling single histogram fit type.
+ * \brief Class for fitting single detector histograms (basic time-differential μSR).
+ *
+ * PRunSingleHisto implements the most fundamental μSR analysis: fitting a single
+ * positron detector histogram to extract relaxation parameters. This is the basis
+ * for all μSR measurements and is used when asymmetry analysis is not appropriate
+ * or desired.
+ *
+ * \section singlehisto_physics Physics and Applications
+ *
+ * <b>Single histogram measurements are used for:</b>
+ * - <b>Time-differential μSR:</b> Measure μ⁺ decay positron time spectrum
+ * - <b>Detector calibration:</b> Characterize individual detector response
+ * - <b>Transverse field (TF):</b> When forward/backward separation unnecessary
+ * - <b>High statistics:</b> Use all positrons (no F-B discrimination)
+ * - <b>Specialized geometries:</b> Non-standard detector arrangements
+ * - <b>Method development:</b> Test fitting algorithms on single detector
+ *
+ * \section singlehisto_data Data Structure and Analysis
+ *
+ * <b>Histogram content:</b>
+ * - Raw counts vs. time from a single positron detector
+ * - Time-zero (t0): Muon arrival time marking start of decay
+ * - Background: Constant or estimated from pre-t0 bins
+ * - Signal: N(t) = N₀·exp(-t/τ_μ)·P(t) + B
+ *
+ * where:
+ * - N₀ = initial count rate (normalization parameter)
+ * - τ_μ = 2.197 μs (muon lifetime)
+ * - P(t) = polarization function (contains physics: relaxation, oscillation, etc.)
+ * - B = background (random coincidences, accidentals)
+ *
+ * \section singlehisto_workflow Analysis Workflow
+ *
+ * 1. <b>Load Histogram:</b> Read raw detector counts from data file
+ * 2. <b>Determine t0:</b> Identify muon arrival time
+ * 3. <b>Background Subtraction:</b>
+ *    - Fixed: Subtract specified constant
+ *    - Estimated: Calculate from pre-t0 bins, subtract with error propagation
+ * 4. <b>Bin Packing:</b> Rebin to improve statistics (optional)
+ * 5. <b>Fit Range:</b> Select time window for parameter extraction
+ * 6. <b>Theory Evaluation:</b> N_theory(t) = N₀·exp(-t/τ_μ)·P_theory(t) + B
+ * 7. <b>Minimization:</b> χ² or maximum likelihood via MINUIT
+ *
+ * \section singlehisto_theory Theory Function
+ *
+ * The fitted function is typically:
+ * \f[ N(t) = N_0 e^{-t/\tau_\mu} P(t) + B \f]
+ *
+ * Common polarization functions P(t):
+ * - <b>Static Gaussian:</b> P(t) = exp(-σ²t²/2)
+ * - <b>Static Lorentzian:</b> P(t) = exp(-λt)
+ * - <b>Dynamic relaxation:</b> P(t) = exp(-(λt)^β)  (stretched exponential)
+ * - <b>Oscillating:</b> P(t) = cos(ωt + φ) · exp(-λt)
+ * - <b>Kubo-Toyabe:</b> Complex relaxation functions for spin systems
+ *
+ * \section singlehisto_parameters Key Parameters
+ *
+ * <b>Normalization (N₀):</b>
+ * - Can be fit parameter or fixed value
+ * - Can be derived from FUNCTIONS block
+ * - Automatically scaled to 1/ns or 1/bin (fScaleN0AndBkg)
+ *
+ * <b>Background (B):</b>
+ * - Fixed: From "background" entry in RUN block
+ * - Estimated: Calculated from pre-t0 bin range
+ * - Units: counts/bin
+ *
+ * <b>Packing:</b>
+ * - Number of consecutive bins to combine
+ * - REQUIRED parameter (RUN or GLOBAL block)
+ * - Higher packing → better statistics, worse time resolution
+ *
+ * \section singlehisto_msr MSR File Example
+ *
+ * \code
+ * RUN data/run2425 PSI MUE4 PSI MUSR-ROOT   (name beamline)
+ *   fittype         0       (SingleHisto)
+ *   map             1       (forward histogram number)
+ *   forward         1
+ *   packing         50      (combine 50 bins → one packed bin)
+ *   background      50 150  (estimate from bins 50-150)
+ *   data            200 2000 (use bins 200-2000 for analysis)
+ *   t0              210.5   (muon arrival bin)
+ *   fit             0.1 10.0 (fit from 0.1 to 10.0 μs after t0)
+ *
+ * THEORY
+ *   asymmetry      1
+ *   simpleGss      2       (Gaussian relaxation σ)
+ *   +              3       (constant background offset)
+ * \endcode
+ *
+ * \section singlehisto_vs_asymmetry Single Histo vs. Asymmetry
+ *
+ * <table>
+ * <tr><th>Feature</th><th>Single Histogram</th><th>Asymmetry</th></tr>
+ * <tr><td>Detectors</td><td>One (forward)</td><td>Two (forward + backward)</td></tr>
+ * <tr><td>Quantity fitted</td><td>N(t) counts</td><td>A(t) = (F-αB)/(F+αB)</td></tr>
+ * <tr><td>Statistics</td><td>All positrons</td><td>F and B separately</td></tr>
+ * <tr><td>Background</td><td>Additive B</td><td>Cancels in asymmetry</td></tr>
+ * <tr><td>α parameter</td><td>N/A</td><td>Required (F/B asymmetry)</td></tr>
+ * <tr><td>Use cases</td><td>TF, calibration</td><td>LF, ZF, weak TF</td></tr>
+ * </table>
+ *
+ * \see PRunAsymmetry for forward-backward asymmetry analysis
+ * \see PRunSingleHistoRRF for rotating reference frame (high-TF)
+ * \see PRunBase for base class interface and common functionality
  */
 class PRunSingleHisto : public PRunBase
 {
   public:
+    /**
+     * \brief Default constructor creating an empty, invalid single histogram run object.
+     *
+     * Initializes all member variables to default/invalid states:
+     * - fScaleN0AndBkg = true (scale to 1/ns by default)
+     * - fNoOfFitBins = 0 (no bins to fit)
+     * - fBackground = 0 (no background)
+     * - fPacking = -1 (unspecified - will cause error if not set)
+     * - fTheoAsData = false (high-resolution theory grid)
+     * - Good bins markers = -1 (unset)
+     * - Fit range bins = -1 (unset)
+     *
+     * This constructor is needed for creating vectors of PRunSingleHisto objects.
+     * The resulting object cannot be used until properly initialized via the main constructor.
+     */
     PRunSingleHisto();
+
+    /**
+     * \brief Main constructor initializing a single histogram run from MSR file and data.
+     *
+     * Performs comprehensive initialization:
+     *
+     * 1. <b>Base Class Initialization:</b>
+     *    - Calls PRunBase constructor
+     *    - Initializes theory engine, parameter mappings, FUNCTIONS block
+     *
+     * 2. <b>N₀/Background Scaling:</b>
+     *    - Calls IsScaleN0AndBkg() to determine scaling mode
+     *    - Sets fScaleN0AndBkg flag (true = scale to 1/ns, false = leave as 1/bin)
+     *
+     * 3. <b>Packing Validation (CRITICAL):</b>
+     *    - Attempts to read packing from RUN block
+     *    - Falls back to GLOBAL block if not in RUN
+     *    - SEVERE ERROR if packing == -1 (mandatory parameter)
+     *    - Marks run invalid and returns if packing not found
+     *
+     * 4. <b>Member Initialization:</b>
+     *    - Good bin markers, fit range bins set to -1 (determined later)
+     *    - Background initialized to 0 (set during PrepareData if specified)
+     *
+     * 5. <b>Data Preparation:</b>
+     *    - Calls PrepareData() to load and preprocess histogram
+     *    - If PrepareData() fails → marks run invalid
+     *
+     * The object is marked as invalid (fValid=false) if:
+     * - Packing parameter is missing from both RUN and GLOBAL blocks
+     * - PrepareData() fails (file not found, invalid t0, etc.)
+     *
+     * \param msrInfo Pointer to MSR file handler (must remain valid)
+     * \param rawData Pointer to raw data handler for histogram loading
+     * \param runNo Run number (0-based index in MSR file RUN blocks)
+     * \param tag Operation mode: kFit (fitting), kView (display/plotting)
+     * \param theoAsData Theory mode: true = at data points, false = high-resolution
+     *
+     * \warning Always check IsValid() after construction. Packing is MANDATORY.
+     *
+     * \see PrepareData() for data preprocessing details
+     * \see IsScaleN0AndBkg() for scaling determination
+     */
     PRunSingleHisto(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData);
+
+    /**
+     * \brief Virtual destructor cleaning up allocated resources.
+     *
+     * Releases memory used by the forward histogram vector (fForward).
+     * Base class destructor handles cleanup of theory objects and other shared resources.
+     */
     virtual ~PRunSingleHisto();
 
+    /**
+     * \brief Calculates χ² between histogram data and theory.
+     *
+     * Computes chi-squared for single histogram fitting:
+     * \f[ \chi^2 = \sum_{i} \frac{(N_i^{\rm data} - N_i^{\rm theory})^2}{\sigma_i^2} \f]
+     *
+     * where N_theory(t) = N₀·exp(-t/τ_μ)·P(t) + B
+     *
+     * Uses OpenMP parallelization when available. N₀ can be a fit parameter
+     * or derived from FUNCTIONS block.
+     *
+     * \param par Parameter vector from MINUIT
+     * \return χ² value (minimize during fitting)
+     */
     virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
+
+    /**
+     * \brief Calculates expected χ² based on theory predictions.
+     * \param par Parameter vector from MINUIT
+     * \return Expected χ² (for statistical diagnostics)
+     */
     virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
+
+    /**
+     * \brief Calculates maximum likelihood for Poisson-distributed histogram counts.
+     *
+     * Computes -2ln(L) for low-count data (< 10-20 counts/bin).
+     * Superior to χ² when Gaussian approximation invalid.
+     *
+     * \param par Parameter vector from MINUIT
+     * \return -2×ln(L) value
+     */
     virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
+
+    /**
+     * \brief Calculates expected maximum likelihood.
+     * \param par Parameter vector from MINUIT
+     * \return Expected -2×ln(L)
+     */
     virtual Double_t CalcMaxLikelihoodExpected(const std::vector<Double_t>& par);
+
+    /**
+     * \brief Evaluates theory function at all data points or high-resolution grid.
+     *
+     * Calculates N_theory(t) = N₀·exp(-t/τ_μ)·P_theory(t) + B
+     * using THEORY block functions. Stores results in fData for χ² calculation.
+     */
     virtual void CalcTheory();
 
+    /**
+     * \brief Returns the number of bins included in the fit range.
+     *
+     * Used for degrees of freedom: ν = N_bins - N_params
+     *
+     * \return Number of bins in fit range
+     */
     virtual UInt_t GetNoOfFitBins();
 
+    /**
+     * \brief Sets fit range using bin-offset specification (COMMANDS block syntax).
+     *
+     * Format: "fit_range fgb+n0 lgb-n1"
+     *
+     * \param fitRange String with bin offsets from good bin markers
+     */
     virtual void SetFitRangeBin(const TString fitRange);
 
+    /**
+     * \brief Returns the estimated background level.
+     * \return Background in counts/bin
+     */
     virtual Double_t GetBackground() { return fBackground; }
 
+    /**
+     * \brief Returns the first bin index in the fit range.
+     * \return Start bin index (0-based, after packing)
+     */
     virtual Int_t GetStartTimeBin() { return fStartTimeBin; }
+
+    /**
+     * \brief Returns the last bin index in the fit range (exclusive).
+     * \return End bin index (loop condition: i < fEndTimeBin)
+     */
     virtual Int_t GetEndTimeBin() { return fEndTimeBin; }
+
+    /**
+     * \brief Returns the bin packing factor.
+     * \return Number of raw bins combined into one packed bin
+     */
     virtual Int_t GetPacking() { return fPacking; }
+
+    /**
+     * \brief Returns the N₀/background scaling mode.
+     * \return true = scaled to 1/ns, false = left as 1/bin
+     */
     virtual Bool_t GetScaleN0AndBkg() { return fScaleN0AndBkg; }
 
+    /**
+     * \brief Calculates start/end bin indices from fit time range.
+     *
+     * Converts fit range (μs) to bin indices, accounting for t0, time resolution,
+     * and packing. Updates fStartTimeBin, fEndTimeBin, fNoOfFitBins.
+     */
     virtual void CalcNoOfFitBins();
 
   protected:
+    /**
+     * \brief Main data preparation orchestrator.
+     *
+     * Coordinates histogram loading and preprocessing: determines operation mode,
+     * calls PrepareFitData() or PrepareViewData(), validates success.
+     *
+     * \return True if data preparation succeeds, false on error
+     */
     virtual Bool_t PrepareData();
+
+    /**
+     * \brief Prepares histogram data for fitting.
+     *
+     * Loads forward histogram, extracts metadata, determines t0, subtracts background,
+     * packs bins, propagates errors, sets up time grid and fit ranges.
+     *
+     * \param runData Raw run data handler
+     * \param histoNo Histogram index in data file
+     * \return True on success, false if preprocessing fails
+     */
     virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo);
+
+    /**
+     * \brief Prepares raw histogram data for viewing (minimal processing).
+     *
+     * Lighter-weight preprocessing for raw histogram visualization without
+     * background subtraction or full fitting infrastructure.
+     *
+     * \param runData Raw run data handler
+     * \param histoNo Histogram index
+     * \return True on success
+     */
     virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo);
+
+    /**
+     * \brief Prepares processed histogram data for viewing/plotting.
+     *
+     * Similar to PrepareFitData() but optimized for visualization with potentially
+     * wider time range for context.
+     *
+     * \param runData Raw run data handler
+     * \param histoNo Histogram index
+     * \return True on success
+     */
     virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo);
 
   private:
-    Bool_t fScaleN0AndBkg;  ///< true=scale N0 and background to 1/ns, otherwise 1/bin
-    UInt_t fNoOfFitBins;    ///< number of bins to be fitted
-    Double_t fBackground;   ///< needed if background range is given (units: 1/bin)
-    Int_t fPacking;         ///< packing for this particular run. Either given in the RUN- or GLOBAL-block.
-    Bool_t fTheoAsData;     ///< true=only calculate the theory points at the data points, false=calculate more points for the theory as compared to data are calculated which lead to 'nicer' Fouriers
+    Bool_t fScaleN0AndBkg;  ///< Scaling mode: true = scale N₀ and B to 1/ns, false = leave as 1/bin (determined by IsScaleN0AndBkg())
+    UInt_t fNoOfFitBins;    ///< Number of bins within fit range (fStartTimeBin to fEndTimeBin)
+    Double_t fBackground;   ///< Background level in counts/bin (estimated from pre-t0 bins or fixed value from RUN block)
+    Int_t fPacking;         ///< Bin packing factor (REQUIRED: from RUN or GLOBAL block)
+    Bool_t fTheoAsData;     ///< Theory mode: true = at data points, false = high-resolution grid for smooth Fourier transforms
 
-    Int_t fGoodBins[2];     ///< keep first/last good bins. 0=fgb, 1=lgb
+    Int_t fGoodBins[2];     ///< Good bin markers for COMMANDS block: [0]=fgb (first good bin/t0), [1]=lgb (last good bin)
 
-    PDoubleVector fForward; ///< forward histo data
+    PDoubleVector fForward; ///< Forward detector histogram (background-corrected, packed)
 
-    Int_t fStartTimeBin;    ///< bin at which the fit starts
-    Int_t fEndTimeBin;      ///< bin at which the fit ends
+    Int_t fStartTimeBin;    ///< First bin index in fit range (inclusive, 0-based after packing)
+    Int_t fEndTimeBin;      ///< Last bin index in fit range (exclusive: loop as i < fEndTimeBin)
 
+    /**
+     * \brief Determines and validates t0 values for histogram.
+     *
+     * Extracts time-zero from RUN block, data file header, GLOBAL block, or
+     * automatic determination. Validates t0 is within histogram bounds.
+     *
+     * \param runData Raw run data
+     * \param globalBlock GLOBAL block settings
+     * \param histoNo Vector of histogram indices
+     * \return True if valid t0 found
+     */
     virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);
+
+    /**
+     * \brief Determines data range (region of valid histogram data).
+     *
+     * Establishes start/end bins for analysis from RUN block "data" entry.
+     * Data range is typically wider than fit range.
+     *
+     * \return True if valid data range determined
+     */
     virtual Bool_t GetProperDataRange();
+
+    /**
+     * \brief Determines fit range from MSR file settings.
+     *
+     * Extracts fit time window from RUN or GLOBAL block "fit" entry.
+     * Format: time-based (μs) or bin-based (fgb+n0 lgb-n1).
+     *
+     * \param globalBlock GLOBAL block with default fit settings
+     */
     virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);
+
+    /**
+     * \brief Estimates initial normalization N₀ from histogram data.
+     *
+     * Calculates N₀ estimate from histogram amplitude, used as starting
+     * value if N₀ is a fit parameter.
+     */
     virtual void EstimateN0();
+
+    /**
+     * \brief Estimates background from pre-t0 bins.
+     *
+     * Calculates background average and error from specified bin range
+     * (typically before t0). Sets fBackground member.
+     *
+     * \param histoNo Histogram index
+     * \return True on success, false if background range invalid
+     */
     virtual Bool_t EstimateBkg(UInt_t histoNo);
+
+    /**
+     * \brief Determines if N₀ and background should be scaled to 1/ns.
+     *
+     * Checks time resolution and fitting preferences to decide scaling mode.
+     * Returns true for standard time bins (scale to 1/ns), false otherwise.
+     *
+     * \return True if scaling should be applied
+     */
     virtual Bool_t IsScaleN0AndBkg();
 };