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 //-------------------------------------------------------------------------- /** - *
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 //-------------------------------------------------------------------------- /** - *
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) //-------------------------------------------------------------------------- /** - *
Calculate chi-square.
+ * \brief Calculates χ² between data and theory (least-squares fit metric).
*
- * return:
- * - 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]
+ *
+ * Algorithm:
+ * -# 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
+ *
+ * N₀ Parameter vs. Function Handling:
+ * - 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
+ *
+ * OpenMP Parallelization:
+ * - Dynamic scheduling with chunk size = (N_bins / N_processors), minimum 10
+ * - Private variables per thread: i, time, diff
+ * - Reduction performed on chisq sum
+ *
+ * Scaling Correction:
+ * 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 Calculate expected chi-square.
+ * \brief Calculates expected χ² using theory as variance (alternative fit metric).
*
- * return:
- * - 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
+ * Difference from Standard χ²:
+ * - 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
+ *
+ * Algorithm:
+ * -# 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
+ *
+ * OpenMP Parallelization:
+ * - 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 Calculate log maximum-likelihood. See http://pdg.lbl.gov/index.html
+ * \brief Calculates -2 log(maximum likelihood) for Poisson-distributed histogram data.
*
- * return:
- * - 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.
+ *
+ * Algorithm:
+ * -# 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
+ *
+ * Edge Cases:
+ * - Zero data (Nᵢ = 0): Uses limit: -2ln(L) → 2×N_theo
+ * - Negative theory: Skips bin and prints warning (should not occur with valid parameters)
+ * - Data threshold: Uses 10⁻⁹ to distinguish zero from non-zero data
+ *
+ * OpenMP Parallelization:
+ * - 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)
+ *
+ * When to Use Maximum Likelihood vs. χ²:
+ * - Use likelihood: Low count rates (< 100 counts/bin), asymmetric errors
+ * - Use χ²: 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 Calculate expected log maximum-likelihood.
+ * \brief Calculates expected -2 log(maximum likelihood) using G-test formulation.
*
- * return:
- * - 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
+ * Difference from CalcMaxLikelihood():
+ * - Full likelihood: includes (theo - data) + data×ln(data/theo)
+ * - Expected likelihood: 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
+ *
+ * Algorithm:
+ * -# 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.
+ *
+ * OpenMP Parallelization:
+ * - 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 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
+ *
+ * Algorithm:
+ * -# 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
+ *
+ * Time Grid:
+ * - 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)
//--------------------------------------------------------------------------
/**
- * Calculate the number of fitted bins for the current fit range.
+ * \brief Returns the number of bins in the current fit range.
*
- * return: 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)
//--------------------------------------------------------------------------
/**
- * 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.
*
- * 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.
+ * Syntax (in COMMAND block):
+ * \code
+ * FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]]
+ * \endcode
+ *
+ * where:
+ * - fgb: First good bin (start of fit range)
+ * - lgb: Last good bin (end of fit range)
+ * - +nXY / -nXY: Optional offsets to shift the range (+ extends, - contracts)
+ * - Multiple pairs: If N+1 pairs given, they apply to each of N RUN blocks
+ *
+ * Two modes:
+ * -# Single pair: `FIT_RANGE fgb lgb` applies to all runs
+ * -# Per-run pairs: `FIT_RANGE fgb₀ lgb₀ fgb₁ lgb₁ ...` applies
+ * pair i to RUN block i
+ *
+ * Algorithm:
+ * -# 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
+ *
+ * Example:
+ * \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)
//--------------------------------------------------------------------------
/**
- * 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.
+ *
+ * Algorithm:
+ * -# 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)
+ *
+ * Edge Cases:
+ * - 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)
//--------------------------------------------------------------------------
/**
- * 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.
*
- * return:
- * - true if everthing went smooth
- * - false, otherwise.
+ * Orchestrates the complete data loading and preprocessing workflow:
+ * -# Load raw data: Fetch run from PRunDataHandler using run name
+ * -# Extract metadata: Magnetic field, beam energy, temperature(s)
+ * -# Validate histograms: Check that forward histogram numbers exist in data file
+ * -# Get time resolution: Extract bin width (typically 0.1-10 ns)
+ * -# Determine t0: Call GetProperT0() for muon arrival times
+ * -# Load histogram data: Copy forward histogram bins from raw data
+ * -# Add runs (ADDRUN): If multiple runs specified, add them with t0 alignment
+ * -# Group histograms: Sum multiple detectors within a group (with t0 alignment)
+ * -# Get data range (fgb/lgb): Call GetProperDataRange() for good bin limits
+ * -# Get fit range: Call GetProperFitRange() for fit time window
+ * -# Check lifetime correction: Determine if exponential decay should be removed (for viewing)
+ * -# Dispatch to preparation:
+ * - kFit → PrepareFitData(): packing, background subtraction
+ * - kView (no lifetime corr.) → PrepareRawViewData(): packing, theory calculation
+ * - kView (with lifetime corr.) → PrepareViewData(): lifetime removal, theory
+ *
+ * ADDRUN t0 Alignment:
+ * 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.
+ *
+ * Grouping t0 Alignment:
+ * 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)
//--------------------------------------------------------------------------
/**
- * 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).
*
- * return:
- * - true, if everything went smooth
- * - false, otherwise
+ * Performs final data transformations after PrepareData() has loaded and grouped
+ * the raw histogram data:
+ * -# Estimate N₀ (optional): If MSR file requests it, call EstimateN0()
+ * -# Handle background:
+ * - 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
+ * -# Packing (rebinning): Combine consecutive bins to improve statistics:
+ * - If packing = 1: copy bins directly
+ * - If packing > 1: sum every 'packing' bins into one
+ * -# Normalization: If fScaleN0AndBkg is true, divide by (packing × t_res × 1000)
+ * to normalize counts to 1/ns
+ * -# Error calculation:
+ * - If N > 0: σ = √N (Poisson statistics)
+ * - If N = 0: σ = 1/normalizer (avoid division by zero in χ²)
+ * -# Set time grid:
+ * - Data start time: (fgb - 0.5 + pack/2 - t0) × t_res
+ * - Data time step: pack × t_res
+ * -# Calculate fit bins: Call CalcNoOfFitBins() to set fStartTimeBin, fEndTimeBin
*
- * \param runData raw run data handler
- * \param histoNo forward histogram number
+ * Packing Algorithm:
+ * \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
+ *
+ * Background Handling Priority:
+ * -# 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)
//--------------------------------------------------------------------------
/**
- * 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:
*
- * return:
- * - true if everthing went smooth
- * - false, otherwise.
+ * Priority hierarchy (highest to lowest):
+ * -# RUN block t0: Explicitly specified in the RUN block (highest priority)
+ * -# GLOBAL block t0: Default t0 for all runs in the GLOBAL block
+ * -# Data file t0: Stored in the raw data file (from previous analysis)
+ * -# Estimated t0: Automatic estimation (UNRELIABLE, prints warning)
+ *
+ * For ADDRUN support:
+ * 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.
+ *
+ * Algorithm:
+ * -# 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
+ *
+ * Validation:
+ * 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)
//--------------------------------------------------------------------------
/**
- * 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).
*
- * return:
- * - 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
+ *
+ * Priority hierarchy (highest to lowest):
+ * -# RUN block: Explicitly specified fgb/lgb in RUN block
+ * -# GLOBAL block: Default fgb/lgb from GLOBAL block
+ * -# Auto-estimation: Fallback estimates with warning
+ *
+ * Auto-estimation (if not specified):
+ * - fgb: t0 + 10 ns (to avoid dead time issues)
+ * - lgb: End of histogram (all bins)
+ *
+ * Validation:
+ * -# 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
+ *
+ * Storage:
+ * 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)
//--------------------------------------------------------------------------
/**
- * Get the proper fit range. There are two possible fit range commands:
- * fit 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.
+ *
+ * When estimation is performed:
+ * - 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)
+ *
+ * When estimation is skipped:
+ * - Norm is a function (paramNo > MSR_PARAM_FUN_OFFSET)
+ * - Norm parameter is fixed (step = 0)
+ * - Invalid parameter number
+ *
+ * Estimation algorithm:
+ * -# 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
+ *
+ * Background handling:
+ * - 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
+ *
+ * Scaling adjustment:
+ * 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)
//--------------------------------------------------------------------------
/**
- * Estimate the background for a given interval.
+ * \brief Estimates background count rate from pre-t0 bins.
*
- * return:
- * - 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
+ * Algorithm:
+ * -# 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)
+ *
+ * Beam periods:
+ * - PSI: 19.75 ns (50.63 MHz cyclotron)
+ * - RAL (ISIS): 320 ns (3.125 MHz target)
+ * - TRIUMF: 43.0 ns (23.26 MHz cyclotron)
+ * - Other facilities: No period correction applied
+ *
+ * Why adjust to beam period?
+ * 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.
+ *
+ * Edge cases:
+ * - 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)
//--------------------------------------------------------------------------
/**
- * 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.
*
- * return:
- * - 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
+ * Default behavior: 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
+ *
+ * To disable scaling: Add to MSR file COMMAND block:
+ * \code
+ * SCALE_N0_BKG FALSE
+ * \endcode
+ *
+ * When to disable scaling:
+ * - 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)
+ *
+ * Effect on fit parameters:
+ * - Scaled (default): N₀ and B independent of packing choice
+ * - Unscaled: N₀ and B depend on packing value
+ *
+ * Implementation details:
+ * 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()
{
diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h
index 5ef885b65..abe7ad46c 100644
--- a/src/include/PRunSingleHisto.h
+++ b/src/include/PRunSingleHisto.h
@@ -33,59 +33,411 @@
#include "PRunBase.h"
/**
- * 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
+ *
+ * Single histogram measurements are used for:
+ * - Time-differential μSR: Measure μ⁺ decay positron time spectrum
+ * - Detector calibration: Characterize individual detector response
+ * - Transverse field (TF): When forward/backward separation unnecessary
+ * - High statistics: Use all positrons (no F-B discrimination)
+ * - Specialized geometries: Non-standard detector arrangements
+ * - Method development: Test fitting algorithms on single detector
+ *
+ * \section singlehisto_data Data Structure and Analysis
+ *
+ * Histogram content:
+ * - 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. Load Histogram: Read raw detector counts from data file
+ * 2. Determine t0: Identify muon arrival time
+ * 3. Background Subtraction:
+ * - Fixed: Subtract specified constant
+ * - Estimated: Calculate from pre-t0 bins, subtract with error propagation
+ * 4. Bin Packing: Rebin to improve statistics (optional)
+ * 5. Fit Range: Select time window for parameter extraction
+ * 6. Theory Evaluation: N_theory(t) = N₀·exp(-t/τ_μ)·P_theory(t) + B
+ * 7. Minimization: χ² 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):
+ * - Static Gaussian: P(t) = exp(-σ²t²/2)
+ * - Static Lorentzian: P(t) = exp(-λt)
+ * - Dynamic relaxation: P(t) = exp(-(λt)^β) (stretched exponential)
+ * - Oscillating: P(t) = cos(ωt + φ) · exp(-λt)
+ * - Kubo-Toyabe: Complex relaxation functions for spin systems
+ *
+ * \section singlehisto_parameters Key Parameters
+ *
+ * Normalization (N₀):
+ * - Can be fit parameter or fixed value
+ * - Can be derived from FUNCTIONS block
+ * - Automatically scaled to 1/ns or 1/bin (fScaleN0AndBkg)
+ *
+ * Background (B):
+ * - Fixed: From "background" entry in RUN block
+ * - Estimated: Calculated from pre-t0 bin range
+ * - Units: counts/bin
+ *
+ * Packing:
+ * - 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
+ *
+ *
+ *
+ *
+ * \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. Base Class Initialization:
+ * - Calls PRunBase constructor
+ * - Initializes theory engine, parameter mappings, FUNCTIONS block
+ *
+ * 2. N₀/Background Scaling:
+ * - Calls IsScaleN0AndBkg() to determine scaling mode
+ * - Sets fScaleN0AndBkg flag (true = scale to 1/ns, false = leave as 1/bin)
+ *
+ * 3. Packing Validation (CRITICAL):
+ * - 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. Member Initialization:
+ * - Good bin markers, fit range bins set to -1 (determined later)
+ * - Background initialized to 0 (set during PrepareData if specified)
+ *
+ * 5. Data Preparation:
+ * - 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::vectorFeature Single Histogram Asymmetry Detectors One (forward) Two (forward + backward) Quantity fitted N(t) counts A(t) = (F-αB)/(F+αB) Statistics All positrons F and B separately Background Additive B Cancels in asymmetry α parameter N/A Required (F/B asymmetry) Use cases TF, calibration LF, ZF, weak TF