From 736c96c66e69fe6ebd85d7ea6d38bddf51ec625b Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 23 Nov 2025 17:25:45 +0100 Subject: [PATCH] improve the doxygen docu of PRunSingleHistoRRF.* --- src/classes/PRunSingleHistoRRF.cpp | 761 +++++++++++++++++++++++++---- src/include/PRunSingleHistoRRF.h | 473 +++++++++++++++++- 2 files changed, 1116 insertions(+), 118 deletions(-) diff --git a/src/classes/PRunSingleHistoRRF.cpp b/src/classes/PRunSingleHistoRRF.cpp index befaf6839..d51b15158 100644 --- a/src/classes/PRunSingleHistoRRF.cpp +++ b/src/classes/PRunSingleHistoRRF.cpp @@ -53,7 +53,22 @@ // Constructor //-------------------------------------------------------------------------- /** - *

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& par) { @@ -189,12 +276,32 @@ Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector& par) // CalcChiSquareExpected (public) //-------------------------------------------------------------------------- /** - *

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& par) { @@ -238,12 +345,35 @@ Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector& // CalcMaxLikelihood (public) //-------------------------------------------------------------------------- /** - *

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& par) { @@ -256,7 +386,39 @@ Double_t PRunSingleHistoRRF::CalcMaxLikelihood(const std::vector& par) // CalcTheory (public) //-------------------------------------------------------------------------- /** - *

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 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 fit time range from MSR file settings. * - * \param globalBlock pointer to the GLOBAL block information form the msr-file. + * Extracts the fit range (time window for parameter extraction) through a + * priority cascade. Supports both time-based and bin-based specifications. + * + * Fit Range Formats: + * - Time-based: \c fit \c 0.1 \c 8.0 (in microseconds relative to t0) + * - Bin-based: \c fit \c fgb+100 \c lgb-500 (bins relative to good bins) + * + * Search Priority: + * -# RUN block time-based fit range + * -# RUN block bin-based fit range (converted to time) + * -# GLOBAL block time-based fit range + * -# GLOBAL block bin-based fit range (converted to time) + * -# Default to data range (fgb to lgb) with warning + * + * Bin-to-Time Conversion: + * \f[ + * t_{\rm start} = ({\rm fgb} + n_0 - t_0) \times \Delta t + * \f] + * \f[ + * t_{\rm end} = ({\rm lgb} - n_1 - t_0) \times \Delta t + * \f] + * + * where Δt = fTimeResolution (raw, before RRF packing). + * + * Result: + * Sets fFitStartTime and fFitEndTime in microseconds relative to t0. + * These define the range used in χ² calculation. + * + * \param globalBlock Pointer to GLOBAL block for fallback fit range settings + * + * \note The converted time values are written back to the MSR data structures + * for log file output consistency. + * + * \warning If no fit range is specified, the full data range is used with + * a warning message. This may not be appropriate for all analyses. + * + * \see GetProperDataRange(), SetFitRangeBin(), CalcNoOfFitBins() */ void PRunSingleHistoRRF::GetProperFitRange(PMsrGlobalBlock *globalBlock) { @@ -1067,9 +1515,41 @@ void PRunSingleHistoRRF::GetProperFitRange(PMsrGlobalBlock *globalBlock) // GetMainFrequency (private) //-------------------------------------------------------------------------- /** - *

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& par); + + /** + * \brief Calculates expected χ² using theory variance instead of data variance. + * + * Computes expected chi-squared where the error estimate comes from the + * theory prediction rather than the data: + * \f[ + * \chi^2_{\rm exp} = \sum_{i} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{A_{\rm RRF}^{\rm theory}(t_i)} + * \f] + * + * This is useful for diagnostic purposes to check fit quality. + * + * \param par Parameter vector from MINUIT + * \return Expected χ² value + */ virtual Double_t CalcChiSquareExpected(const std::vector& par); + + /** + * \brief Calculates maximum likelihood (not yet implemented for RRF). + * + * Maximum likelihood estimation for RRF data is more complex due to the + * transformation from Poisson-distributed counts to RRF asymmetry. + * Currently returns 0.0. + * + * \param par Parameter vector from MINUIT + * \return 0.0 (not implemented) + * + * \todo Implement proper ML for RRF data accounting for error propagation + * through the RRF transformation. + */ virtual Double_t CalcMaxLikelihood(const std::vector& par); + + /** + * \brief Evaluates theory function at all data points for viewing/plotting. + * + * Calculates RRF theory values using current fit parameters stored in + * the MSR parameter list. The theory function P(t) is evaluated at each + * data time point and results are stored in fData for display. + * + * \note Theory is evaluated directly in the RRF frame, not transformed. + * The THEORY block should specify the low-frequency RRF signal. + */ virtual void CalcTheory(); + /** + * \brief Returns the number of bins included in the current fit range. + * + * Triggers CalcNoOfFitBins() to ensure bin count is current, then returns + * the cached value. + * + * \return Number of RRF-packed bins within [fFitStartTime, fFitEndTime] + */ virtual UInt_t GetNoOfFitBins(); + /** + * \brief Sets fit range using bin-offset syntax from COMMANDS block. + * + * Parses fit range specification in the format: + * - "FIT_RANGE fgb+n0 lgb-n1" (single range for all runs) + * - "FIT_RANGE fgb+n00 lgb-n01 fgb+n10 lgb-n11 ..." (per-run ranges) + * + * where fgb = first good bin, lgb = last good bin, and nXY are offsets. + * + * Converts bin specification to time using: + * - fFitStartTime = (fgb + offset - t0) × time_resolution + * - fFitEndTime = (lgb - offset - t0) × time_resolution + * + * \param fitRange String containing FIT_RANGE specification + * + * \note This is called when COMMANDS block modifies fit range during fitting. + */ virtual void SetFitRangeBin(const TString fitRange); + /** + * \brief Returns the first bin index in the fit range. + * \return Start bin index (0-based, after RRF 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 Calculates start/end bin indices from fit time range. + * + * Converts fit range times (μs) to RRF-packed bin indices: + * - fStartTimeBin = ceil((fFitStartTime - dataTimeStart) / dataTimeStep) + * - fEndTimeBin = floor((fFitEndTime - dataTimeStart) / dataTimeStep) + 1 + * + * Also updates fNoOfFitBins = fEndTimeBin - fStartTimeBin. + * + * \note Time step accounts for RRF packing: dataTimeStep = rawTimeStep × fRRFPacking + */ virtual void CalcNoOfFitBins(); protected: + /** + * \brief Main data preparation orchestrator for RRF single histogram analysis. + * + * Coordinates loading and preprocessing of histogram data: + * -# Validates GLOBAL block presence and RRF parameters + * -# Retrieves raw run data from data handler + * -# Extracts metadata (field, energy, temperature) + * -# Collects histogram numbers from RUN block forward specification + * -# Validates histograms exist in data file + * -# Determines time resolution (ns → μs conversion) + * -# Determines t0 values via GetProperT0() + * -# Handles addrun (co-adding multiple runs) + * -# Handles grouping (combining multiple detectors) + * -# Determines data range via GetProperDataRange() + * -# Determines fit range via GetProperFitRange() + * -# Calls PrepareFitData() or PrepareViewData() based on tag + * + * \return true if data preparation succeeds, false on any error + * + * \see PrepareFitData(), PrepareViewData() + */ virtual Bool_t PrepareData(); + + /** + * \brief Performs full RRF transformation for fitting. + * + * Transforms raw histogram to RRF asymmetry through these steps: + * + * 1. Frequency Analysis: + * - Calls GetMainFrequency() to find dominant precession frequency + * - Used to determine optimal N₀ estimation window + * + * 2. Background Handling: + * - If background range given: estimate from data via EstimateBkg() + * - If fixed background given: use directly + * - Subtract background: N(t) → N(t) - B + * + * 3. Lifetime Correction: + * - Apply exp(+t/τ_μ) to remove exponential decay + * - Store as M(t) = [N(t) - B] × exp(+t/τ_μ) + * - Calculate M_err = exp(+t/τ_μ) × √(N(t) + σ_B²) + * + * 4. N₀ Estimation: + * - Call EstimateN0() to determine normalization + * - Uses weighted average over full oscillation cycles + * + * 5. Asymmetry Extraction: + * - A(t) = M(t)/N₀ - 1 + * - A_err(t) = exp(+t/τ)/N₀ × √(N(t) + [(N(t)-B)/N₀]² × σ_N₀²) + * + * 6. RRF Rotation: + * - A_RRF(t) = 2 × A(t) × cos(ω_RRF × t + φ_RRF) + * - Factor 2 compensates for discarded counter-rotating component + * + * 7. RRF Packing: + * - Average over fRRFPacking consecutive bins + * - Error: σ_packed = √(2 × Σσ²) / n + * + * 8. Time Grid Setup: + * - Set data time start accounting for packing offset + * - Set data time step = raw_resolution × fRRFPacking + * + * \param runData Raw run data handler + * \param histoNo Forward histogram number (0-based) + * \return true on success, false on error + */ virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo); + + /** + * \brief Prepares RRF data for viewing/plotting. + * + * Similar to PrepareFitData() but additionally: + * - Handles view packing (if specified and > RRF packing) + * - Sets up theory calculation grid (8× finer than data if fTheoAsData=false) + * - Evaluates theory function for plotting overlay + * + * \param runData Raw run data handler + * \param histoNo Forward histogram number + * \return true on success, false on error + * + * \note View packing < RRF packing is ignored with warning since RRF + * packing is already applied during data transformation. + */ virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo); private: - Double_t fN0EstimateEndTime; ///< end time in (us) over which N0 is estimated. + Double_t fN0EstimateEndTime; ///< End time (μs) for N₀ estimation window. Rounded to integer number of oscillation cycles based on main frequency. - UInt_t fNoOfFitBins; ///< number of bins to be fitted - Double_t fBackground; ///< needed if background range is given (units: 1/bin) - Double_t fBkgErr; ///< estimate error on the estimated background - Int_t fRRFPacking; ///< RRF packing for this particular run. Given in the 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 + UInt_t fNoOfFitBins; ///< Number of RRF-packed bins within fit range [fStartTimeBin, fEndTimeBin) + Double_t fBackground; ///< Estimated or fixed background level in counts/bin (before packing) + Double_t fBkgErr; ///< Statistical error on background estimate (std dev of background region) + Int_t fRRFPacking; ///< RRF packing factor from GLOBAL block (number of raw bins averaged together) + Bool_t fTheoAsData; ///< Theory resolution mode: true = at data points only, false = 8× finer grid for smooth Fourier transforms - Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb + Int_t fGoodBins[2]; ///< Good bin range: [0] = first good bin (fgb), [1] = last good bin (lgb). Used for COMMANDS block fit range specification. - 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 in RRF-packed data) + Int_t fEndTimeBin; ///< Last bin index in fit range (exclusive: loop as i < fEndTimeBin) - PDoubleVector fForward; ///< forward histo data - PDoubleVector fM; ///< vector holding M(t) = [N(t)-N_bkg] exp(+t/tau). Needed to estimate N0. - PDoubleVector fMerr; ///< vector holding the error of M(t): M_err = exp(+t/tau) sqrt(N(t)). - PDoubleVector fW; ///< vector holding the weight needed to estimate N0, and errN0. - PDoubleVector fAerr; ///< vector holding the errors of estimated A(t) + PDoubleVector fForward; ///< Forward detector histogram data (progressively transformed during preparation) + PDoubleVector fM; ///< Lifetime-corrected histogram: M(t) = [N(t) - B] × exp(+t/τ_μ). Used for N₀ estimation. + PDoubleVector fMerr; ///< Error on M(t): σ_M = exp(+t/τ_μ) × √(N(t) + σ_B²). Includes background error. + PDoubleVector fW; ///< Weights for N₀ estimation: W(t) = 1/σ_M². Used in weighted average. + PDoubleVector fAerr; ///< Asymmetry errors before RRF packing. Used for packed error calculation. + /** + * \brief Determines and validates t0 values for all histograms. + * + * Searches for t0 in order of priority: + * -# RUN block t0 specification + * -# GLOBAL block t0 specification + * -# Data file header t0 values + * -# Automatic estimation (with warning) + * + * Also handles addt0 for addrun histograms. + * + * \param runData Raw run data for histogram access + * \param globalBlock GLOBAL block with default t0 values + * \param histoNo Vector of histogram indices to process + * \return true if valid t0 found for all histograms, false on error + */ virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo); + + /** + * \brief Determines valid data range (first/last good bins). + * + * Establishes analysis boundaries in order of priority: + * -# RUN block data range specification + * -# GLOBAL block data range specification + * -# Automatic estimation (t0 + 10ns offset to end of histogram) + * + * Validates that range is sensible (start < end, within histogram bounds). + * + * \return true if valid data range determined, false on error + */ virtual Bool_t GetProperDataRange(); + + /** + * \brief Determines fit time range from MSR file settings. + * + * Extracts fit range in order of priority: + * -# RUN block fit range (time-based or bin-based) + * -# GLOBAL block fit range (time-based or bin-based) + * -# Default to data range (fgb to lgb) with warning + * + * Converts bin-based specifications to time using t0 and time resolution. + * + * \param globalBlock GLOBAL block with default fit settings + */ virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock); + + /** + * \brief Finds the dominant precession frequency in raw data. + * + * Performs Fourier transform on raw histogram data to identify the + * main muon precession frequency: + * -# Creates TH1F from data in good bin range + * -# Applies strong apodization (windowing) + * -# Computes power spectrum + * -# Searches for maximum above 10 MHz (ignores DC component) + * + * The frequency is used to: + * - Determine optimal N₀ estimation window (integer cycles) + * - Calculate suggested "optimal packing" for user information + * + * \param data Raw histogram data vector + * \return Maximum frequency in MHz, or 0 if not found + */ virtual Double_t GetMainFrequency(PDoubleVector &data); + + /** + * \brief Estimates initial normalization N₀ from lifetime-corrected data. + * + * Calculates N₀ as weighted average of M(t) over initial time window: + * \f[ + * N_0 = \frac{\sum_{i=0}^{n} M(t_i)}{n} + * \f] + * + * The window length is chosen to include an integer number of complete + * oscillation cycles (based on freqMax) within fN0EstimateEndTime. + * + * Error estimate: + * \f[ + * \sigma_{N_0} = \frac{\sqrt{\sum w_i^2 \sigma_{M_i}^2}}{\sum w_i} + * \f] + * + * \param errN0 [out] Estimated error on N₀ + * \param freqMax Main precession frequency (MHz) for window calculation + * \return Estimated N₀ value + */ virtual Double_t EstimateN0(Double_t &errN0, Double_t freqMax); + + /** + * \brief Estimates background from pre-t0 bins. + * + * Calculates background level and error from specified background range: + * -# Validates background range is within histogram bounds + * -# Adjusts range to integer number of accelerator RF cycles (PSI/RAL/TRIUMF) + * -# Computes mean and standard deviation + * -# Stores in fBackground and fBkgErr + * + * \param histoNo Histogram index (for error messages) + * \return true on success, false if background range invalid + */ virtual Bool_t EstimateBkg(UInt_t histoNo); };