Constructor + * \brief Default constructor for RRF single histogram fitting class. + * + * Initializes all member variables to safe default values: + * - fNoOfFitBins = 0 (no bins to fit) + * - fBackground = 0 (will be estimated or set from MSR file) + * - fBkgErr = 1.0 (default error estimate) + * - fRRFPacking = -1 (invalid until set from GLOBAL block) + * - fTheoAsData = false (high-resolution theory grid) + * - fGoodBins[0,1] = -1 (calculated from data range) + * - fN0EstimateEndTime = 1.0 μs (default N₀ estimation window) + * + * \warning This constructor creates an invalid object that cannot be used + * until properly initialized with MSR file data. Use the full + * constructor for normal operation. + * + * \see PRunSingleHistoRRF(PMsrHandler*, PRunDataHandler*, UInt_t, EPMusrHandleTag, Bool_t) */ PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase() { @@ -75,12 +90,44 @@ PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase() // Constructor //-------------------------------------------------------------------------- /** - *
Constructor + * \brief Main constructor for RRF 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 RRF single histogram run object by: + * -# Validating GLOBAL block presence (mandatory for RRF analysis) + * -# Validating RRF frequency specification (rrf_freq in GLOBAL block) + * -# Validating RRF packing specification (rrf_packing in GLOBAL block) + * -# Calling PrepareData() to load histogram and perform RRF transformation + * + * GLOBAL Block Requirements: + * The RRF fit type requires the following entries in the GLOBAL block: + * - \c rrf_freq: Rotation frequency with unit (e.g., "13.554 T", "183.7 MHz") + * - \c rrf_packing: Number of bins to average (integer) + * - \c rrf_phase: (optional) Initial phase in degrees + * + * Error Handling: + * If any validation fails, the constructor: + * - Outputs detailed error message to stderr + * - Sets fValid = false + * - Returns immediately (PrepareData() is not called) + * + * \param msrInfo Pointer to MSR file handler (NOT owned, must outlive this object) + * \param rawData Pointer to raw run data handler (NOT owned, must outlive this object) + * \param runNo Zero-based index of the RUN block in the MSR file + * \param tag Operation mode: kFit (fitting) or kView (viewing/plotting) + * \param theoAsData If true, theory calculated only at data points (for viewing); + * if false, theory uses finer time grid (8× data resolution) + * + * \warning GLOBAL block with RRF parameters is MANDATORY for this fit type. + * Always check IsValid() after construction. + * + * \note After construction, check IsValid() to ensure initialization succeeded. + * Common failure modes: + * - Missing GLOBAL block + * - Missing rrf_freq specification + * - Missing rrf_packing specification + * - Data file not found or histogram missing + * + * \see PrepareData(), PrepareFitData(), PrepareViewData() */ PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) : PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData) @@ -131,7 +178,17 @@ PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *ra // Destructor //-------------------------------------------------------------------------- /** - *
Destructor + * \brief Destructor for RRF single histogram fitting class. + * + * Cleans up dynamically allocated memory: + * - Clears the forward histogram data vector (fForward) + * - Other vectors (fM, fMerr, fW, fAerr) are local to PrepareFitData + * and cleared automatically + * + * Base class destructor (PRunBase) handles cleanup of: + * - Theory objects + * - Function value arrays + * - Other shared resources */ PRunSingleHistoRRF::~PRunSingleHistoRRF() { @@ -142,12 +199,42 @@ PRunSingleHistoRRF::~PRunSingleHistoRRF() // CalcChiSquare (public) //-------------------------------------------------------------------------- /** - *
Calculate chi-square.
+ * \brief Calculates χ² between RRF-transformed data and theory (least-squares fit metric).
*
- * return:
- * - chisq value
+ * Computes the standard chi-square goodness-of-fit statistic for RRF asymmetry:
+ * \f[
+ * \chi^2 = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{\sigma_i^2}
+ * \f]
*
- * \param par parameter vector iterated by minuit2
+ * Unlike standard single histogram fitting, no explicit N₀ or exponential decay
+ * factors appear since the RRF transformation already produces dimensionless
+ * asymmetry with properly propagated errors.
+ *
+ * Algorithm:
+ * -# Evaluate all user-defined functions from FUNCTIONS block
+ * -# Pre-evaluate theory at t=1.0 to initialize any stateful functions
+ * (e.g., LF relaxation, user functions with internal state)
+ * -# Loop over fit range bins [fStartTimeBin, fEndTimeBin)
+ * -# For each bin: calculate time, evaluate theory, accumulate χ²
+ *
+ * OpenMP Parallelization:
+ * When compiled with OpenMP (HAVE_GOMP defined):
+ * - Dynamic scheduling with chunk size = max(10, N_bins / N_processors)
+ * - Private variables per thread: i, time, diff
+ * - Reduction performed on chisq accumulator
+ * - Thread-safe due to pre-evaluation of theory at t=1.0
+ *
+ * \param par Parameter vector from MINUIT minimizer, containing current
+ * estimates of all fit parameters
+ *
+ * \return χ² value (sum over all bins in fit range). Minimize this value
+ * during fitting to find optimal parameters.
+ *
+ * \note The theory function is evaluated in the RRF frame. The THEORY block
+ * should describe the low-frequency RRF signal, not the laboratory frame
+ * precession.
+ *
+ * \see CalcChiSquareExpected(), CalcMaxLikelihood()
*/
Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector Calculate expected chi-square.
+ * \brief Calculates expected χ² using theory variance instead of data variance.
*
- * return:
- * - chisq value
+ * Computes the expected chi-square where the error estimate in the denominator
+ * comes from the theory prediction rather than the data:
+ * \f[
+ * \chi^2_{\rm exp} = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{A_{\rm RRF}^{\rm theory}(t_i)}
+ * \f]
*
- * \param par parameter vector iterated by minuit2
+ * This metric is useful for:
+ * - Diagnostic purposes to assess fit quality
+ * - Detecting systematic deviations from the model
+ * - Comparing with standard χ² to identify error estimation issues
+ *
+ * Algorithm:
+ * Same as CalcChiSquare() but uses theory value as variance estimate instead
+ * of measured error bars. OpenMP parallelization is applied when available.
+ *
+ * \param par Parameter vector from MINUIT minimizer
+ *
+ * \return Expected χ² value. For a good fit, this should be approximately
+ * equal to the number of degrees of freedom (N_bins - N_params).
+ *
+ * \warning Theory values must be positive for valid variance estimate.
+ * Negative theory values can lead to incorrect χ² calculation.
+ *
+ * \see CalcChiSquare()
*/
Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector Calculate log maximum-likelihood. See http://pdg.lbl.gov/index.html
+ * \brief Calculates maximum likelihood for RRF data (NOT YET IMPLEMENTED).
*
- * return:
- * - log maximum-likelihood value
+ * Maximum likelihood estimation for RRF single histogram data is more complex
+ * than for raw histograms due to the non-linear transformation from
+ * Poisson-distributed counts to RRF asymmetry.
*
- * \param par parameter vector iterated by minuit2
+ * Theoretical Background:
+ * For raw histogram data, the likelihood is:
+ * \f[
+ * \mathcal{L} = \prod_i \frac{\mu_i^{n_i} e^{-\mu_i}}{n_i!}
+ * \f]
+ * where \f$\mu_i\f$ is the expected count and \f$n_i\f$ is the observed count.
+ *
+ * For RRF-transformed data, the error propagation through the transformation
+ * must be properly accounted for in the likelihood function.
+ *
+ * Current Implementation:
+ * Returns 0.0 (not implemented). Use χ² minimization (CalcChiSquare) instead.
+ *
+ * \param par Parameter vector from MINUIT minimizer (unused)
+ *
+ * \return 0.0 (not implemented)
+ *
+ * \todo Implement proper maximum likelihood for RRF data by:
+ * -# Deriving the likelihood function for transformed asymmetry
+ * -# Accounting for error propagation through RRF transformation
+ * -# Including correlations introduced by packing
+ *
+ * \see CalcChiSquare() for currently supported fit metric
*/
Double_t PRunSingleHistoRRF::CalcMaxLikelihood(const std::vector Calculate theory for a given set of fit-parameters.
+ * \brief Evaluates theory function at all data points for viewing/plotting.
+ *
+ * Calculates the theoretical RRF asymmetry using the current MSR parameter
+ * values and stores results in fData for display. This method is called
+ * after fitting to generate the theory curve overlay.
+ *
+ * Algorithm:
+ * -# Extract parameter values from MSR parameter list
+ * -# Evaluate all user-defined functions from FUNCTIONS block
+ * -# Loop over data points (size matches RRF-packed data)
+ * -# Calculate time: t = dataTimeStart + i × dataTimeStep
+ * -# Evaluate theory: P(t) = Func(t, par, funcValues)
+ * -# Store results via fData.AppendTheoryValue()
+ *
+ * Theory Function:
+ * The theory is evaluated directly in the RRF frame. The THEORY block should
+ * specify the low-frequency RRF signal (after transformation), not the
+ * laboratory-frame high-frequency precession.
+ *
+ * Example THEORY block for RRF analysis:
+ * \code
+ * THEORY
+ * asymmetry 1 (amplitude)
+ * simpleGss 2 (Gaussian relaxation)
+ * TFieldCos map1 fun1 (frequency shift, phase)
+ * \endcode
+ *
+ * where the frequency in TFieldCos is the difference frequency (ω - ω_RRF).
+ *
+ * \note Unlike CalcChiSquare(), this method does not return a value.
+ * Results are stored internally in fData.fTheory vector.
+ *
+ * \see CalcChiSquare(), PrepareViewData()
*/
void PRunSingleHistoRRF::CalcTheory()
{
@@ -289,9 +451,23 @@ void PRunSingleHistoRRF::CalcTheory()
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
- * Calculate the number of fitted bins for the current fit range.
+ * \brief Returns the number of bins included in the current fit range.
*
- * return: number of fitted bins.
+ * Triggers CalcNoOfFitBins() to ensure the bin count is current based on
+ * the latest fit range settings, then returns the cached value.
+ *
+ * The number of fit bins is needed for:
+ * - Calculating degrees of freedom: ν = N_bins - N_params
+ * - Reduced χ²: χ²_red = χ² / ν
+ * - Statistical diagnostics and fit quality assessment
+ *
+ * \return Number of RRF-packed bins within [fFitStartTime, fFitEndTime].
+ * This accounts for RRF packing: fewer bins than raw data.
+ *
+ * \note The fit range may be modified during fitting by COMMANDS block
+ * instructions. Always call this method to get the current count.
+ *
+ * \see CalcNoOfFitBins(), SetFitRangeBin()
*/
UInt_t PRunSingleHistoRRF::GetNoOfFitBins()
{
@@ -304,15 +480,51 @@ UInt_t PRunSingleHistoRRF::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 Sets fit range using bin-offset syntax from COMMANDS block.
*
- * nXY are offsets which can be used to shift, limit the fit range.
+ * Allows dynamic modification of the fit range during fitting, as specified
+ * in the COMMANDS block. This is used for systematic studies where the fit
+ * range needs to be varied.
*
- * \param fitRange string containing the necessary information.
+ * Syntax:
+ * \code
+ * FIT_RANGE fgb[+n0] lgb[-n1] # single range for all runs
+ * FIT_RANGE fgb+n00 lgb-n01 fgb+n10 lgb-n11 ... # per-run ranges
+ * \endcode
+ *
+ * where:
+ * - \c fgb = first good bin (symbolic, replaced by actual value)
+ * - \c lgb = last good bin (symbolic, replaced by actual value)
+ * - \c +n0 = positive offset added to fgb
+ * - \c -n1 = positive offset subtracted from lgb
+ *
+ * Conversion to Time:
+ * \f[
+ * t_{\rm start} = ({\rm fgb} + n_0 - t_0) \times \Delta t
+ * \f]
+ * \f[
+ * t_{\rm end} = ({\rm lgb} - n_1 - t_0) \times \Delta t
+ * \f]
+ *
+ * where Δt is the raw time resolution (before RRF packing).
+ *
+ * Multiple Run Handling:
+ * - Single pair: Applied to all runs
+ * - Multiple pairs: Must match number of RUN blocks; each run uses its own range
+ * - Run selection: Uses (2 × runNo + 1) to index into token array
+ *
+ * Example:
+ * \code
+ * COMMANDS
+ * FIT_RANGE fgb+100 lgb-500 # start 100 bins after fgb, end 500 bins before lgb
+ * \endcode
+ *
+ * \param fitRange String containing FIT_RANGE specification from COMMANDS block
+ *
+ * \warning Invalid syntax produces error message but does not throw exception.
+ * The previous fit range values are retained on error.
+ *
+ * \see GetProperFitRange(), CalcNoOfFitBins()
*/
void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange)
{
@@ -394,7 +606,37 @@ void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange)
// CalcNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
- * Calculate the number of fitted bins for the current fit range.
+ * \brief Calculates start/end bin indices from fit time range.
+ *
+ * Converts the fit range from time (μs) to RRF-packed bin indices.
+ * This method is called whenever the fit range may have changed.
+ *
+ * Conversion Formulas:
+ * \f[
+ * i_{\rm start} = \lceil \frac{t_{\rm start} - t_{\rm data,0}}{\Delta t_{\rm data}} \rceil
+ * \f]
+ * \f[
+ * i_{\rm end} = \lfloor \frac{t_{\rm end} - t_{\rm data,0}}{\Delta t_{\rm data}} \rfloor + 1
+ * \f]
+ *
+ * where:
+ * - \f$t_{\rm data,0}\f$ = fData.GetDataTimeStart() (first RRF-packed bin center)
+ * - \f$\Delta t_{\rm data}\f$ = fData.GetDataTimeStep() (RRF-packed bin width)
+ *
+ * Bounds Checking:
+ * - fStartTimeBin clamped to [0, data size)
+ * - fEndTimeBin clamped to [0, data size]
+ * - fNoOfFitBins = 0 if fEndTimeBin <= fStartTimeBin
+ *
+ * Side Effects:
+ * Updates member variables:
+ * - fStartTimeBin: First bin index in fit range (inclusive)
+ * - fEndTimeBin: Last bin index in fit range (exclusive)
+ * - fNoOfFitBins: Number of bins = fEndTimeBin - fStartTimeBin
+ *
+ * \note Time step includes RRF packing: dataTimeStep = rawTimeRes × fRRFPacking
+ *
+ * \see GetNoOfFitBins(), SetFitRangeBin()
*/
void PRunSingleHistoRRF::CalcNoOfFitBins()
{
@@ -416,17 +658,50 @@ void PRunSingleHistoRRF::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 preparation orchestrator for RRF single histogram analysis.
*
- * return:
- * - true if everthing went smooth
- * - false, otherwise.
+ * Coordinates the loading and preprocessing of histogram data before RRF
+ * transformation. This method handles all operations common to both fitting
+ * and viewing modes.
+ *
+ * Processing Steps:
+ * -# Validate object state (return false if already marked invalid)
+ * -# Get GLOBAL block reference for settings
+ * -# Load raw run data from data handler using run name from MSR file
+ * -# Extract metadata from data file:
+ * - Magnetic field (fMetaData.fField)
+ * - Beam energy (fMetaData.fEnergy)
+ * - Sample temperature(s) (fMetaData.fTemp)
+ * -# Collect histogram numbers from RUN block forward specification
+ * -# Validate all specified histograms exist in data file
+ * -# Determine time resolution: ns → μs conversion
+ * -# Determine t0 values via GetProperT0()
+ * -# Initialize forward histogram from first group
+ * -# Handle addrun (co-add multiple runs with t0 alignment)
+ * -# Handle grouping (combine multiple detectors with t0 alignment)
+ * -# Determine data range via GetProperDataRange()
+ * -# Determine fit range via GetProperFitRange()
+ * -# Dispatch to PrepareFitData() or PrepareViewData() based on tag
+ *
+ * Run Addition (addrun):
+ * When multiple runs are specified in the RUN block, histograms are co-added
+ * with t0 alignment:
+ * \code
+ * fForward[i] += addRunData[i + addT0 - mainT0]
+ * \endcode
+ *
+ * Detector Grouping:
+ * When multiple forward histograms are specified, they are summed with
+ * t0 alignment to form a single combined histogram.
+ *
+ * \return true if data preparation succeeds, false on any error:
+ * - Run data not found
+ * - Histogram not present in data file
+ * - Invalid t0 determination
+ * - Data range validation failure
+ * - Fit/view data preparation failure
+ *
+ * \see PrepareFitData(), PrepareViewData(), GetProperT0(), GetProperDataRange()
*/
Bool_t PRunSingleHistoRRF::PrepareData()
{
@@ -560,21 +835,73 @@ Bool_t PRunSingleHistoRRF::PrepareData()
// PrepareFitData (protected)
//--------------------------------------------------------------------------
/**
- * Take the pre-processed data (i.e. grouping and addrun are preformed) and form the RRF 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.
- * -# estimate N0
- * -# RRF transformation
- * -# packing (i.e rebinning)
+ * \brief Performs full RRF transformation on histogram data for fitting.
*
- * return:
- * - true, if everything went smooth
- * - false, otherwise
+ * Takes the pre-processed histogram (grouping and addrun already applied) and
+ * transforms it into RRF asymmetry suitable for fitting. This is the core
+ * method implementing the RRF analysis technique.
*
- * \param runData raw run data handler
- * \param histoNo forward histogram number
+ * Processing Steps:
+ *
+ * 1. Frequency Analysis:
+ * - Calls GetMainFrequency() on raw data to find dominant precession frequency
+ * - Used to determine optimal N₀ estimation window (integer oscillation cycles)
+ * - Prints "optimal packing" suggestion for user reference
+ *
+ * 2. Background Handling:
+ * - If fixed background specified: use directly from RUN block
+ * - If background range given: estimate via EstimateBkg()
+ * - If neither: estimate range as [0.1×t0, 0.6×t0] with warning
+ * - Subtract background: N(t) → N(t) - B
+ *
+ * 3. Lifetime Correction:
+ * \f[
+ * M(t) = [N(t) - B] \cdot e^{+t/\tau_\mu}
+ * \f]
+ * - Removes exponential decay from signal
+ * - Error: \f$\sigma_M = e^{+t/\tau_\mu} \sqrt{N(t) + \sigma_B^2}\f$
+ * - Weights: \f$w = 1/\sigma_M^2\f$
+ *
+ * 4. N₀ Estimation:
+ * - Call EstimateN0() with frequency information
+ * - Uses average of M(t) over integer number of oscillation cycles
+ * - Returns N₀ and its error σ_N₀
+ *
+ * 5. Asymmetry Extraction:
+ * \f[
+ * A(t) = \frac{M(t)}{N_0} - 1
+ * \f]
+ * Error propagation:
+ * \f[
+ * \sigma_A(t) = \frac{e^{+t/\tau_\mu}}{N_0} \sqrt{N(t) + \left(\frac{N(t)-B}{N_0}\right)^2 \sigma_{N_0}^2}
+ * \f]
+ *
+ * 6. RRF Rotation:
+ * \f[
+ * A_{\rm RRF}(t) = 2 \cdot A(t) \cdot \cos(\omega_{\rm RRF} t + \phi_{\rm RRF})
+ * \f]
+ * - ω_RRF from GLOBAL block (rrf_freq converted to angular frequency)
+ * - φ_RRF from GLOBAL block (rrf_phase in radians)
+ * - Factor 2 compensates for discarded counter-rotating component
+ *
+ * 7. RRF Packing:
+ * - Average over fRRFPacking consecutive bins
+ * - Data: \f$A_{\rm packed} = \frac{1}{n}\sum_{i=1}^{n} A_{\rm RRF}(t_i)\f$
+ * - Error: \f$\sigma_{\rm packed} = \frac{\sqrt{2}}{n}\sqrt{\sum_{i=1}^{n} \sigma_A^2(t_i)}\f$
+ * - √2 factor accounts for doubling in RRF rotation
+ *
+ * 8. Time Grid Setup:
+ * - Data time start: accounts for packing offset (center of first packed bin)
+ * - Data time step: rawTimeResolution × fRRFPacking
+ *
+ * \param runData Raw run data handler (for background estimation)
+ * \param histoNo Forward histogram number (0-based, for background error messages)
+ *
+ * \return true on success, false on error:
+ * - Background estimation failure
+ * - Invalid bin ranges
+ *
+ * \see PrepareViewData(), GetMainFrequency(), EstimateN0(), EstimateBkg()
*/
Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t histoNo)
{
@@ -712,20 +1039,55 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his
// PrepareViewData (protected)
//--------------------------------------------------------------------------
/**
- * Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing
- * with life time correction, i.e. the exponential decay is removed.
- * The following steps are preformed:
- * -# check if view packing is whished.
- * -# check that 'first good data bin', 'last good data bin', and 't0' makes any sense
- * -# transform data sets (see below).
- * -# calculate theory
+ * \brief Prepares RRF-transformed data for viewing/plotting.
*
- * return:
- * - true, if everything went smooth
- * - false, otherwise
+ * Takes the pre-processed histogram and prepares both data and theory curves
+ * for display. This method is used when the operation mode is kView rather
+ * than kFit.
*
- * \param runData raw run data handler
- * \param histoNo forward histogram number
+ * Processing Steps:
+ *
+ * 1. Data Preparation:
+ * - Calls PrepareFitData() to perform full RRF transformation
+ * - Data is ready for display after this step
+ *
+ * 2. View Packing Check:
+ * - Checks if additional view packing is specified in PLOT block
+ * - If viewPacking > fRRFPacking: additional packing would be applied (TODO)
+ * - If viewPacking < fRRFPacking: warning issued and view packing ignored
+ * (cannot unpack already-packed RRF data)
+ *
+ * 3. Theory Grid Setup:
+ * Determines theory evaluation resolution:
+ * - If fTheoAsData = true: theory calculated only at data points
+ * - If fTheoAsData = false: theory calculated on 8× finer grid for smooth curves
+ *
+ * Time grid:
+ * - Theory time start = data time start
+ * - Theory time step = data time step / factor (where factor = 1 or 8)
+ *
+ * 4. Theory Evaluation:
+ * - Extract parameter values from MSR parameter list
+ * - Evaluate all user-defined functions from FUNCTIONS block
+ * - Loop over theory grid points
+ * - Evaluate theory: P(t) = Func(t, par, funcValues)
+ * - Apply sanity check: |P(t)| > 10 → set to 0 (dirty hack for edge cases)
+ * - Store results via fData.AppendTheoryValue()
+ *
+ * Theory Function:
+ * The theory is evaluated directly in the RRF frame. The THEORY block should
+ * specify the low-frequency signal after RRF transformation, not the
+ * laboratory-frame high-frequency precession.
+ *
+ * \param runData Raw run data handler (passed to PrepareFitData)
+ * \param histoNo Forward histogram number (passed to PrepareFitData)
+ *
+ * \return true on success, false on error (typically from PrepareFitData)
+ *
+ * \note The 8× theory resolution provides smoother curves for visualization
+ * and better Fourier transforms without affecting fit performance.
+ *
+ * \see PrepareFitData(), CalcTheory()
*/
Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t histoNo)
{
@@ -793,21 +1155,52 @@ Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t hi
// 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 and validates t0 values for all histograms in the run.
*
- * \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
+ * Searches for time-zero (t0) values through a priority cascade and performs
+ * validation. t0 marks the muon arrival time and is critical for proper
+ * timing alignment.
*
- * return:
- * - true if everthing went smooth
- * - false, otherwise.
+ * t0 Search Priority:
+ * -# RUN block: t0 values specified in the MSR file RUN block
+ * -# GLOBAL block: fallback t0 values from GLOBAL block
+ * -# Data file: t0 values stored in the raw data file header
+ * -# Automatic estimation: estimated from histogram shape (with warning)
+ *
+ * Algorithm:
+ * -# Initialize fT0s vector with size = number of forward histograms (grouping)
+ * -# Set all t0 values to -1.0 (unset marker)
+ * -# Fill from RUN block if present
+ * -# Fill remaining -1.0 values from GLOBAL block
+ * -# Fill remaining -1.0 values from data file
+ * -# Fill remaining -1.0 values from automatic estimation (with warning)
+ * -# Validate all t0 values are within histogram bounds
+ *
+ * Addrun Handling:
+ * When multiple runs are co-added (addrun), each additional run needs its own
+ * t0 value for proper time alignment:
+ * -# Initialize fAddT0s[runIdx] for each addrun
+ * -# Apply same priority cascade for each addrun
+ * -# Validate addrun t0 values
+ *
+ * Validation:
+ * - t0 must be ≥ 0
+ * - t0 must be < histogram length
+ * - Error message and return false on validation failure
+ *
+ * \param runData Pointer to raw run data for histogram access and data file t0
+ * \param globalBlock Pointer to GLOBAL block for fallback t0 values
+ * \param histoNo Vector of histogram indices (0-based) for forward grouping
+ *
+ * \return true if valid t0 found for all histograms, false on:
+ * - t0 out of histogram bounds
+ * - addrun data not accessible
+ * - addrun t0 validation failure
+ *
+ * \warning Automatic t0 estimation is unreliable for LEM and other specialized
+ * setups. Always specify t0 explicitly for best results.
+ *
+ * \see GetProperDataRange(), PrepareData()
*/
Bool_t PRunSingleHistoRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo)
{
@@ -931,14 +1324,40 @@ Bool_t PRunSingleHistoRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *gl
// 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 valid data range (first/last good bins) for analysis.
*
- * return:
- * - true if everthing went smooth
- * - false, otherwise.
+ * Establishes the boundaries of usable histogram data through a priority
+ * cascade. The data range defines which bins contain valid signal (after t0,
+ * before noise/overflow).
+ *
+ * Data Range Search Priority:
+ * -# RUN block: data range specified in MSR file RUN block
+ * -# GLOBAL block: fallback data range from GLOBAL block
+ * -# Automatic estimation: estimate from t0 (with warning)
+ * - Start: t0 + 10 ns offset
+ * - End: histogram length
+ *
+ * Validation:
+ * -# Check start < end (swap if reversed)
+ * -# Check start ≥ 0 and start < histogram size
+ * -# Check end ≥ 0 (adjust if > histogram size with warning)
+ *
+ * Result:
+ * Sets fGoodBins[0] (first good bin) and fGoodBins[1] (last good bin).
+ * These are used for:
+ * - RRF transformation range
+ * - Fit range specification in bins (fgb+n0 lgb-n1)
+ * - Data validity checking
+ *
+ * \return true if valid data range determined, false on:
+ * - Start bin out of bounds
+ * - End bin negative
+ * - Other range validation failures
+ *
+ * \note The data range is typically wider than the fit range. Data range
+ * defines where valid data exists; fit range defines what is fitted.
+ *
+ * \see GetProperFitRange(), GetProperT0()
*/
Bool_t PRunSingleHistoRRF::GetProperDataRange()
{
@@ -1019,16 +1438,45 @@ Bool_t PRunSingleHistoRRF::GetProperDataRange()
// GetProperFitRange (private)
//--------------------------------------------------------------------------
/**
- * Get the proper fit range. There are two possible fit range commands:
- * fit gets the maximum frequency of the given data.
+ * \brief Finds the dominant precession frequency in raw histogram data.
*
- * \param raw data set.
+ * Performs Fourier analysis on the raw histogram to identify the main muon
+ * spin precession frequency. This is essential for:
+ * - Determining optimal N₀ estimation window (integer oscillation cycles)
+ * - Calculating suggested "optimal packing" for user information
+ * - Validating that the RRF frequency is appropriate
+ *
+ * Algorithm:
+ * -# Create ROOT TH1F histogram from data in good bin range [fGoodBins[0], fGoodBins[1]]
+ * -# Set histogram binning to match time structure
+ * -# Apply Fourier transform using PFourier class
+ * -# Use strong apodization (windowing) to reduce spectral leakage
+ * -# Search power spectrum for maximum above 10 MHz (ignores DC component)
+ * -# Return frequency at maximum power
+ *
+ * Frequency Search Constraints:
+ * - Ignores frequencies below 10 MHz (DC and low-frequency noise)
+ * - Searches for local maximum (rising then falling power)
+ * - Returns 0 if no clear maximum found
+ *
+ * Output Information:
+ * The method prints diagnostic information:
+ * - Detected maximum frequency (MHz)
+ * - Suggested optimal packing for ~5-8 points per cycle
+ *
+ * \param data Raw histogram data vector (counts per bin)
+ *
+ * \return Maximum frequency in MHz, or 0 if no maximum found above 10 MHz.
+ *
+ * \note The frequency is used to determine the N₀ estimation window:
+ * window = ceil(fN0EstimateEndTime × freqMax / 2π) × (2π / freqMax)
+ * This ensures an integer number of complete oscillation cycles.
+ *
+ * \see EstimateN0(), PrepareFitData()
*/
Double_t PRunSingleHistoRRF::GetMainFrequency(PDoubleVector &data)
{
@@ -1117,9 +1597,50 @@ Double_t PRunSingleHistoRRF::GetMainFrequency(PDoubleVector &data)
// EstimateN0 (private)
//--------------------------------------------------------------------------
/**
- * Estimate the N0 for the given run.
+ * \brief Estimates initial normalization N₀ from lifetime-corrected histogram data.
*
- * \param errN0
+ * Calculates N₀ by averaging the lifetime-corrected signal M(t) over an initial
+ * time window. The window is chosen to span an integer number of complete
+ * precession cycles to minimize bias from oscillations.
+ *
+ * N₀ Estimation Formula:
+ * \f[
+ * N_0 = \frac{1}{n} \sum_{i=0}^{n-1} M(t_i)
+ * \f]
+ *
+ * where M(t) = [N(t) - B] × exp(+t/τ_μ) is the lifetime-corrected histogram
+ * stored in fM vector.
+ *
+ * Window Determination:
+ * The end bin is calculated to include an integer number of oscillation cycles:
+ * \f[
+ * n_{\rm end} = {\rm round}\left( \left\lceil \frac{T \cdot f_{\rm max}}{2\pi} \right\rceil \cdot \frac{2\pi}{f_{\rm max} \cdot \Delta t} \right)
+ * \f]
+ *
+ * where:
+ * - T = fN0EstimateEndTime (default 1.0 μs)
+ * - f_max = main precession frequency from GetMainFrequency()
+ * - Δt = fTimeResolution
+ *
+ * Error Estimation:
+ * \f[
+ * \sigma_{N_0} = \frac{\sqrt{\sum_{i=0}^{n-1} w_i^2 \sigma_{M_i}^2}}{\sum_{i=0}^{n-1} w_i}
+ * \f]
+ *
+ * where w_i = 1/σ_M_i² are the weights stored in fW vector.
+ *
+ * Output Information:
+ * Prints diagnostic message: "N₀ = value (error)"
+ *
+ * \param errN0 [out] Reference to store the estimated error on N₀
+ * \param freqMax Main precession frequency (MHz) from GetMainFrequency()
+ *
+ * \return Estimated N₀ value (counts × exp(+t/τ))
+ *
+ * \note The simple average (rather than weighted average) is used for N₀ itself,
+ * while the error uses the full weight information.
+ *
+ * \see GetMainFrequency(), PrepareFitData()
*/
Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax)
{
@@ -1152,13 +1673,47 @@ Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax)
// EstimatBkg (private)
//--------------------------------------------------------------------------
/**
- * Estimate the background for a given interval.
+ * \brief Estimates background level and error from pre-t0 histogram bins.
*
- * return:
- * - true, if everything went smooth
- * - false, otherwise
+ * Calculates the background (random coincidences, accidentals) from a
+ * specified range of histogram bins, typically before t0. The background
+ * range is adjusted to span an integer number of accelerator RF cycles
+ * for proper averaging.
*
- * \param histoNo forward histogram number of the run
+ * Accelerator RF Periods:
+ * - PSI: ACCEL_PERIOD_PSI ns
+ * - RAL: ACCEL_PERIOD_RAL ns
+ * - TRIUMF: ACCEL_PERIOD_TRIUMF ns
+ * - Unknown: no RF adjustment
+ *
+ * Algorithm:
+ * -# Get background range [start, end] from RUN block
+ * -# Swap start/end if in wrong order (with message)
+ * -# Determine accelerator RF period from institute name
+ * -# Adjust end to span integer number of RF cycles
+ * -# Validate range is within histogram bounds
+ * -# Calculate mean background:
+ * \f$ B = \frac{1}{n}\sum_{i={\rm start}}^{{\rm end}} N(i) \f$
+ * -# Calculate background error (standard deviation):
+ * \f$ \sigma_B = \sqrt{\frac{1}{n-1}\sum_{i={\rm start}}^{{\rm end}} (N(i) - B)^2} \f$
+ * -# Store results in fBackground and fBkgErr
+ * -# Update RUN block with estimated background value
+ *
+ * Output Information:
+ * Prints diagnostic messages:
+ * - Adjusted background range (if RF adjustment applied)
+ * - Background value and error: "fBackground = value (error)"
+ *
+ * \param histoNo Forward histogram number (for error message context)
+ *
+ * \return true on success, false on error:
+ * - Background start out of histogram bounds
+ * - Background end out of histogram bounds
+ *
+ * \note The RF cycle adjustment ensures proper averaging over the
+ * accelerator bunch structure, reducing systematic bias.
+ *
+ * \see PrepareFitData()
*/
Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo)
{
diff --git a/src/include/PRunSingleHistoRRF.h b/src/include/PRunSingleHistoRRF.h
index 91c3aad57..947beba22 100644
--- a/src/include/PRunSingleHistoRRF.h
+++ b/src/include/PRunSingleHistoRRF.h
@@ -33,59 +33,502 @@
#include "PRunBase.h"
/**
- * Class handling single histogram fit type.
+ * \brief Class for fitting single histogram data in a Rotating Reference Frame (RRF).
+ *
+ * PRunSingleHistoRRF implements single histogram analysis with Rotating Reference Frame
+ * transformation. This technique is essential for high transverse field (TF) μSR
+ * measurements where the precession frequency is too high to resolve directly.
+ *
+ * \section rrf_physics Physics Background
+ *
+ * The Problem with High-Field μSR:
+ * - At fields B > ~0.5 T, muon precession frequency ω = γ_μ B exceeds ~70 MHz
+ * - Standard time-domain analysis struggles: requires sub-ns binning → low statistics
+ * - Frequency-domain (Fourier) analysis has limited resolution
+ *
+ * The RRF Solution:
+ * - Transform to a reference frame rotating at frequency ω_RRF ≈ ω_signal
+ * - In rotating frame: effective frequency Δω = ω_signal - ω_RRF << ω_signal
+ * - Low-frequency signal can be fit with standard time-domain techniques
+ * - Preserves full information: relaxation rates, frequency distributions
+ *
+ * \section rrf_transform Mathematical Transformation
+ *
+ * Data processing steps:
+ * -# Background subtraction: N(t) → N(t) - B
+ * -# Lifetime correction: N(t) - B → [N(t) - B] × exp(+t/τ_μ) = M(t)
+ * -# N₀ estimation: Fit M(t) over initial time window to extract N₀
+ * -# Asymmetry extraction: A(t) = M(t)/N₀ - 1
+ * -# RRF rotation: A_RRF(t) = 2 × A(t) × cos(ω_RRF t + φ_RRF)
+ * -# RRF packing: Average over multiple bins to filter high-frequency components
+ *
+ * The factor of 2 in step 5 compensates for the loss of the counter-rotating component
+ * when the high-frequency part is suppressed by packing.
+ *
+ * Theory function in RRF:
+ * The polarization function P(t) transforms as:
+ * \f[
+ * P_{\rm RRF}(t) = 2 \cdot P(t) \cdot \cos(\omega_{\rm RRF} t + \phi_{\rm RRF})
+ * \f]
+ *
+ * For a simple oscillation P(t) = A·cos(ωt + φ)·exp(-λt), this becomes:
+ * \f[
+ * P_{\rm RRF}(t) = A \cdot [\cos((\omega - \omega_{\rm RRF})t + \phi - \phi_{\rm RRF}) + \text{high freq.}] \cdot e^{-\lambda t}
+ * \f]
+ *
+ * After RRF packing, only the low-frequency component Δω = ω - ω_RRF remains.
+ *
+ * \section rrf_requirements GLOBAL Block Requirements
+ *
+ * RRF analysis REQUIRES the following entries in the GLOBAL block:
+ * - rrf_freq: RRF rotation frequency (value + unit: MHz, Mc, T, G, kG)
+ * - rrf_phase: RRF initial phase in degrees
+ * - rrf_packing: Number of bins to average (filters high frequencies)
+ *
+ * \section rrf_msr MSR File Example
+ *
+ * \code
+ * GLOBAL
+ * rrf_freq 13.554 T # RRF frequency equivalent to 13.554 Tesla
+ * rrf_phase 45.0 # RRF phase in degrees
+ * rrf_packing 25 # Average 25 bins to suppress high-freq component
+ * data 200 60000 # Data range (fgb lgb)
+ * fit 0.1 8.0 # Fit range in μs
+ *
+ * RUN data/run2425 PSI MUE4 PSI MUSR-ROOT
+ * fittype 8 (SingleHistoRRF)
+ * map 1
+ * forward 1
+ * background 50 150
+ * t0 210.5
+ * \endcode
+ *
+ * \section rrf_workflow Data Processing Workflow
+ *
+ * 1. PrepareData(): Load raw histogram, determine t0, group/add histograms
+ * 2. PrepareFitData():
+ * - Subtract background (estimated or fixed)
+ * - Apply lifetime correction: exp(+t/τ_μ)
+ * - Estimate N₀ using Fourier analysis to find main frequency
+ * - Extract asymmetry: A(t) = M(t)/N₀ - 1
+ * - Apply RRF rotation: A_RRF(t) = 2·A(t)·cos(ω_RRF t + φ_RRF)
+ * - Pack data (average over rrf_packing bins)
+ * - Calculate errors with proper propagation
+ *
+ * \section rrf_errors Error Propagation
+ *
+ * RRF asymmetry error (unpacked):
+ * \f[
+ * \sigma_{A}(t) = \frac{e^{t/\tau_\mu}}{N_0} \sqrt{N(t) + \left(\frac{N(t)-B}{N_0}\right)^2 \sigma_{N_0}^2}
+ * \f]
+ *
+ * After RRF packing (n bins):
+ * \f[
+ * \sigma_{A_{\rm RRF}}^{\rm packed} = \frac{\sqrt{2}}{n} \sqrt{\sum_{i=1}^{n} \sigma_{A}^2(t_i)}
+ * \f]
+ *
+ * The √2 factor accounts for the doubling in the RRF transformation.
+ *
+ * \section rrf_applications Applications
+ *
+ * High-TF μSR experiments:
+ * - Knight shift measurements (NMR-like)
+ * - Vortex state studies in type-II superconductors
+ * - Diamagnetic/paramagnetic shift measurements
+ * - Internal field distributions in magnets
+ *
+ * \see PRunSingleHisto for standard (non-RRF) single histogram analysis
+ * \see PRunAsymmetryRRF for RRF asymmetry (forward-backward) analysis
+ * \see PRunBase for base class interface
*/
class PRunSingleHistoRRF : public PRunBase
{
public:
+ /**
+ * \brief Default constructor creating an empty, invalid RRF single histogram run object.
+ *
+ * Initializes member variables to default/invalid states:
+ * - fNoOfFitBins = 0
+ * - fBackground = 0
+ * - fBkgErr = 1.0 (default error estimate)
+ * - fRRFPacking = -1 (invalid, must be set from GLOBAL block)
+ * - fTheoAsData = false
+ * - fGoodBins[0,1] = -1 (unset)
+ * - fN0EstimateEndTime = 1.0 μs
+ *
+ * This constructor exists for container compatibility. The resulting object
+ * cannot be used until properly initialized.
+ */
PRunSingleHistoRRF();
+
+ /**
+ * \brief Main constructor initializing RRF single histogram run from MSR file and data.
+ *
+ * Performs comprehensive validation and initialization:
+ *
+ * 1. GLOBAL Block Validation (MANDATORY):
+ * - Checks GLOBAL block is present (required for RRF)
+ * - Validates RRF frequency is specified (rrf_freq)
+ * - Validates RRF packing is specified (rrf_packing)
+ * - Sets fValid=false and returns on any validation failure
+ *
+ * 2. Member Initialization:
+ * - Extracts fRRFPacking from GLOBAL block
+ * - Sets fN0EstimateEndTime = 1.0 μs
+ * - Initializes fGoodBins to -1 (determined in PrepareData)
+ *
+ * 3. Data Preparation:
+ * - Calls PrepareData() to load histogram and perform RRF transformation
+ * - Sets fValid=false if data preparation fails
+ *
+ * \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 resolution: true = at data points, false = 8× finer grid
+ *
+ * \warning Always check IsValid() after construction. GLOBAL block with RRF
+ * parameters is MANDATORY for this fit type.
+ *
+ * \see PrepareData() for RRF transformation details
+ */
PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData);
+
+ /**
+ * \brief Virtual destructor releasing allocated resources.
+ *
+ * Clears the forward histogram vector. Base class handles theory objects.
+ */
virtual ~PRunSingleHistoRRF();
+ /**
+ * \brief Calculates χ² between RRF-transformed data and theory.
+ *
+ * Computes chi-squared for RRF single histogram fitting:
+ * \f[
+ * \chi^2 = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{\sigma_i^2}
+ * \f]
+ *
+ * Algorithm:
+ * -# Evaluate user-defined FUNCTIONS with current parameters
+ * -# Pre-evaluate theory at t=1.0 to initialize LF/user functions (thread-safety)
+ * -# Loop over fit range with OpenMP parallelization
+ * -# Sum squared differences weighted by inverse variance
+ *
+ * OpenMP Parallelization:
+ * - Dynamic scheduling with chunk = max(10, N_bins / N_processors)
+ * - Reduction on chisq accumulator
+ *
+ * \param par Parameter vector from MINUIT minimizer
+ * \return χ² value (sum over all bins in fit range)
+ *
+ * \note Unlike standard single histogram, no N₀/exp(-t/τ) factors needed since
+ * RRF transformation already produces dimensionless asymmetry.
+ */
virtual Double_t CalcChiSquare(const std::vector