diff --git a/src/classes/PRunMuMinus.cpp b/src/classes/PRunMuMinus.cpp index 1c6b7c2b9..337fa2c72 100644 --- a/src/classes/PRunMuMinus.cpp +++ b/src/classes/PRunMuMinus.cpp @@ -47,7 +47,19 @@ // Constructor //-------------------------------------------------------------------------- /** - *
Constructor + * \brief Default constructor creating an empty, invalid μ⁻ run object. + * + * Initializes all member variables to default/invalid states: + * - Bin counts and indices set to -1 (invalid) + * - Packing set to -1 (unspecified - will cause error if used) + * - Theory mode set to false (high-resolution grid) + * - Handle tag set to kEmpty (uninitialized) + * + * This constructor is needed for creating vectors of PRunMuMinus objects, + * but the resulting object cannot be used for fitting until properly + * initialized via the main constructor. + * + * \see PRunMuMinus(PMsrHandler*, PRunDataHandler*, UInt_t, EPMusrHandleTag, Bool_t) */ PRunMuMinus::PRunMuMinus() : PRunBase() { @@ -70,12 +82,43 @@ PRunMuMinus::PRunMuMinus() : PRunBase() // Constructor //-------------------------------------------------------------------------- /** - *
Constructor + * \brief Main constructor initializing a μ⁻ run from MSR file and raw data. * - * \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 + * Performs comprehensive initialization for negative muon analysis: + * + * 1. Base Class Initialization: + * - Calls PRunBase constructor with MSR/data handlers + * - Initializes theory engine and parameter mappings + * + * 2. Packing Validation (CRITICAL for μ⁻): + * - Attempts to read packing from RUN block + * - Falls back to GLOBAL block if not in RUN block + * - SEVERE ERROR if packing == -1 (unspecified) + * - Packing is MANDATORY for μ⁻ data (unlike some other fit types) + * + * 3. Member Initialization: + * - Good bin markers set to -1 (determined later) + * - Fit range bins initialized to -1 (invalid until set) + * - Number of fit bins set to 0 + * + * 4. Data Preparation: + * - Calls PrepareData() to load and preprocess histogram + * - If PrepareData() fails, run is marked invalid + * + * The object is marked as invalid (fValid=false) if: + * - Packing parameter is missing + * - PrepareData() fails (file not found, invalid t0, etc.) + * + * \param msrInfo Pointer to MSR file handler (must remain valid for object lifetime) + * \param rawData Pointer to raw data handler for loading histogram files + * \param runNo Run number (0-based index in MSR file RUN blocks) + * \param tag Operation mode: kFit (fitting), kView (display/plotting) + * \param theoAsData Theory calculation mode: true = at data points, false = high-resolution + * + * \warning Invalid objects should not be used for fitting. Check IsValid() after construction. + * + * \see PrepareData() for data preprocessing details + * \see PRunBase constructor for base initialization */ PRunMuMinus::PRunMuMinus(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) : PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData) @@ -114,7 +157,10 @@ PRunMuMinus::PRunMuMinus(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t // Destructor //-------------------------------------------------------------------------- /** - *
Destructor + * \brief Destructor cleaning up allocated resources. + * + * Releases memory used by the forward histogram vector. Other cleanup + * (theory objects, base class resources) is handled by the PRunBase destructor. */ PRunMuMinus::~PRunMuMinus() { @@ -125,12 +171,38 @@ PRunMuMinus::~PRunMuMinus() // CalcChiSquare //-------------------------------------------------------------------------- /** - *
Calculate chi-square.
+ * \brief Calculates χ² between μ⁻ data and theory (least-squares fit metric).
*
- * return:
- * - chisq value
+ * Computes the chi-squared statistic using the standard formula:
+ * \f[ \chi^2 = \sum_{i={\rm start}}^{\rm end} \frac{(N_i^{\rm data} - N_i^{\rm theory})^2}{\sigma_i^2} \f]
*
- * \param par parameter vector iterated by minuit2
+ * Algorithm:
+ * 1. Evaluate FUNCTIONS block for user-defined functions
+ * 2. Pre-calculate theory once at t=1.0 (thread-safe initialization for LF/user functions)
+ * 3. Loop over fit range bins (fStartTimeBin to fEndTimeBin)
+ * 4. For each bin:
+ * - Calculate time from bin index
+ * - Evaluate theory function at that time
+ * - Compute squared difference weighted by error
+ * 5. Sum contributions (with OpenMP reduction if available)
+ *
+ * OpenMP Parallelization:
+ * - When compiled with GOMP support, χ² calculation is parallelized
+ * - Loop divided into dynamic chunks for load balancing
+ * - Chunk size: (N_bins / N_processors), minimum 10 bins per chunk
+ * - Private variables per thread: i, time, diff
+ * - Reduction performed on final chisq sum
+ *
+ * Theory Pre-calculation:
+ * The initial call to fTheory->Func(time=1.0, ...) ensures thread-safe
+ * initialization for London field (LF) and user functions that cache
+ * computationally expensive calculations per parameter set.
+ *
+ * \param par Parameter vector from MINUIT with current parameter values
+ * \return Chi-squared value (sum of weighted squared residuals)
+ *
+ * \see CalcMaxLikelihood() for Poisson-based alternative (better for low counts)
+ * \see CalcTheory() for theory evaluation without χ² calculation
*/
Double_t PRunMuMinus::CalcChiSquare(const std::vector Calculate expected chi-square.
+ * \brief Calculates expected χ² assuming theory is the true distribution (diagnostic).
*
- * return:
- * - chisq value == 0.0
+ * Computes the expected chi-squared using theory values as the expected counts:
+ * \f[ \chi^2_{\rm expected} = \sum_{i} \frac{(N_i^{\rm data} - N_i^{\rm theory})^2}{N_i^{\rm theory}} \f]
*
- * \param par parameter vector iterated by minuit2
+ * This is a statistical diagnostic for evaluating:
+ * - Quality of error estimates (if errors are correct, χ²/ν ≈ 1)
+ * - Validity of Gaussian approximation (breaks down for low counts)
+ * - Over/under-dispersion in data relative to Poisson expectations
+ *
+ * For Poisson-distributed data with large counts:
+ * - χ²_expected ≈ number of degrees of freedom
+ * - χ²_expected / ν ≈ 1 indicates proper error estimation
+ *
+ * Algorithm is identical to CalcChiSquare() except:
+ * - Error denominator is σ² = N_theory (Poisson variance)
+ * - Instead of using fData.GetError()->at(i)
+ *
+ * OpenMP Parallelization: Same as CalcChiSquare()
+ *
+ * \param par Parameter vector from MINUIT
+ * \return Expected χ² value (currently returns 0.0 - calculation done but not returned)
+ *
+ * \note Current implementation performs the calculation but returns 0.0.
+ * The calculated value could be returned for diagnostic purposes.
+ *
+ * \see CalcChiSquare() for standard χ² calculation
*/
Double_t PRunMuMinus::CalcChiSquareExpected(const std::vector Calculate log max-likelihood. See http://pdg.lbl.gov/index.html
+ * \brief Calculates negative log-likelihood for Poisson statistics (low-count fit metric).
*
- * return:
- * - log max-likelihood value
+ * Computes the maximum likelihood estimator assuming Poisson-distributed histogram counts:
+ * \f[ -2\ln L = 2\sum_{i} \left[N_i^{\rm theory} - N_i^{\rm data}\ln(N_i^{\rm theory})\right] \f]
*
- * \param par parameter vector iterated by minuit2
+ * This is derived from the Poisson probability:
+ * \f[ P(n|\\lambda) = \frac{\\lambda^n e^{-\\lambda}}{n!} \f]
+ *
+ * Taking negative log-likelihood and multiplying by 2 gives a metric that:
+ * - Is minimized at the best-fit parameters (like χ²)
+ * - Approaches χ² in the high-count (Gaussian) limit
+ * - Handles low counts correctly (no Gaussian approximation needed)
+ *
+ * When to use likelihood vs. χ²:
+ * - Low counts (< 10-20 per bin): Use likelihood (this method)
+ * - High counts (> 20 per bin): Either method works, χ² is faster
+ * - Uncertain errors: Use likelihood (no error estimates needed)
+ *
+ * Algorithm:
+ * 1. Evaluate FUNCTIONS block
+ * 2. Pre-calculate theory at t=1.0 (thread-safe initialization)
+ * 3. Loop over fit range bins
+ * 4. For each bin:
+ * - Calculate theory prediction
+ * - Check for negative theory (warning + skip)
+ * - Add Poisson likelihood term:
+ * - If data > 10^-9: (theo - data) + data·ln(data/theo)
+ * - If data ≈ 0: (theo - data) only (avoid log(0))
+ * 5. Multiply sum by 2 (convention for comparison with χ²)
+ *
+ * Edge Case Handling:
+ * - Negative theory: Print warning, skip bin (should not happen with proper model)
+ * - Zero data: Use simplified formula without log term
+ * - Zero theory: Would cause log(0), but checked via negative theory guard
+ *
+ * Reference: Particle Data Group (PDG), Statistics Review
+ * http://pdg.lbl.gov/index.html
+ *
+ * \param par Parameter vector from MINUIT
+ * \return -2×ln(L) value (minimize during fitting, comparable to χ²)
+ *
+ * \see CalcChiSquare() for standard least-squares metric
*/
Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector Calculate the number of fitted bins for the current fit range.
+ * \brief Returns the number of bins included in the fit range.
*
- * return: number of fitted bins.
+ * Recalculates and returns the count of histogram bins between fStartTimeBin
+ * and fEndTimeBin. This count is used for:
+ * - Degrees of freedom calculation: ν = N_bins - N_params
+ * - Reduced χ²: χ²_red = χ² / ν
+ * - Statistical quality assessment
+ *
+ * Internally calls CalcNoOfFitBins() to ensure fNoOfFitBins is up-to-date.
+ *
+ * \return Number of bins within fit range (fEndTimeBin - fStartTimeBin)
+ *
+ * \see CalcNoOfFitBins() for bin range calculation from time range
*/
UInt_t PRunMuMinus::GetNoOfFitBins()
{
diff --git a/src/include/PRunMuMinus.h b/src/include/PRunMuMinus.h
index ed58c57dd..75731ef71 100644
--- a/src/include/PRunMuMinus.h
+++ b/src/include/PRunMuMinus.h
@@ -33,49 +33,455 @@
#include "PRunBase.h"
/**
- * Will eventually handle the \f$\mu^{-}\f$ fitting etc.
+ * \brief Class for handling negative muon (μ⁻) histogram fits.
+ *
+ * PRunMuMinus implements single-histogram fitting specialized for negative muon
+ * measurements. Negative muons (μ⁻) have significantly different properties compared
+ * to positive muons (μ⁺), requiring dedicated handling:
+ *
+ * \section muminus_physics Physics of Negative Muons
+ *
+ * Key differences from μ⁺:
+ * - Charge: Negatively charged (antimuon has positive charge)
+ * - Behavior in matter: μ⁻ are captured by nuclei rather than stopping interstitially
+ * - Cascade process: μ⁻ cascade down atomic orbitals before nuclear capture
+ * - Lifetime: Effective lifetime is reduced by nuclear capture (typically ~100-1000 ns vs. 2.2 μs free lifetime)
+ * - Decay signature: Different detector response due to capture vs. decay
+ * - Applications: Elemental analysis, nuclear structure, muonic atom spectroscopy
+ *
+ * \section muminus_data Data Structure
+ *
+ * Unlike asymmetry measurements, μ⁻ fits use a single forward histogram:
+ * - One detector histogram (typically forward detector)
+ * - No asymmetry calculation (no forward/backward pair)
+ * - Direct histogram fitting with background subtraction
+ * - Bin packing for improved statistics
+ *
+ * \section muminus_analysis Analysis Workflow
+ *
+ * 1. Data Loading: Load single histogram from data file
+ * 2. Background Subtraction: Remove constant or estimated background
+ * 3. Time-Zero Determination: Identify μ⁻ arrival time (t0)
+ * 4. Bin Packing: Rebin data if specified (required parameter)
+ * 5. Theory Calculation: Evaluate capture-modified exponential decay
+ * 6. Fit Metric: χ² or maximum likelihood minimization
+ *
+ * \section muminus_theory Typical Theory Functions
+ *
+ * Common models for μ⁻ data include:
+ * - Simple exponential: N(t) = N₀ · exp(-t/τ_eff)
+ * - With relaxation: N(t) = N₀ · exp(-t/τ_eff) · [1 + A·cos(ωt + φ)]
+ * - Multi-component: Sum of exponentials for different capture sites
+ *
+ * where τ_eff combines free decay and nuclear capture rates:
+ * \f[ \frac{1}{\tau_{\rm eff}} = \frac{1}{\tau_{\mu}} + \lambda_{\rm capture} \f]
+ *
+ * \section muminus_usage Example Usage
+ *
+ * \code
+ * // In MSR file RUN block:
+ * RUN data/muminus_run2425 MUD PSI MUSR-ROOT (name beamline)
+ * fittype 10 (MuMinus)
+ * map 1 (forward histogram)
+ * forward 1
+ * packing 10 (required!)
+ * background 50 150 (pre-t0 background)
+ * data 200 2000
+ * t0 210.5
+ * fit 0.1 10.0
+ * \endcode
+ *
+ * \see PRunSingleHisto for standard positive muon single histogram fits
+ * \see PRunBase for base class interface and common functionality
*/
class PRunMuMinus : public PRunBase
{
public:
+ /**
+ * \brief Default constructor creating an empty, invalid μ⁻ run object.
+ *
+ * Creates an uninitialized run object with all values set to defaults.
+ * This constructor is needed for creating vectors of PRunMuMinus objects
+ * but the object is not usable until properly initialized via the main constructor.
+ */
PRunMuMinus();
+
+ /**
+ * \brief Main constructor initializing a μ⁻ run from MSR file and raw data.
+ *
+ * Performs comprehensive initialization:
+ * 1. Validates packing parameter (REQUIRED for μ⁻ - cannot be -1)
+ * 2. Initializes base class (PRunBase) with MSR and data handlers
+ * 3. Calls PrepareData() to load and preprocess the histogram
+ * 4. Validates successful initialization
+ *
+ * The packing parameter is mandatory for μ⁻ data and is obtained from:
+ * - RUN block "packing" entry (if specified)
+ * - Falls back to GLOBAL block "packing" entry
+ * - SEVERE ERROR if neither is specified → run marked invalid
+ *
+ * \param msrInfo Pointer to MSR file handler (must remain valid)
+ * \param rawData Pointer to raw data handler for histogram loading
+ * \param runNo Run number (0-based index in MSR file RUN blocks)
+ * \param tag Operation mode: kFit (fitting), kView (display/plotting)
+ * \param theoAsData Theory mode: true = at data points, false = high-resolution
+ *
+ * \warning If PrepareData() fails, the object is marked invalid (fValid=false)
+ */
PRunMuMinus(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData);
+
+ /**
+ * \brief Virtual destructor cleaning up allocated resources.
+ *
+ * Releases memory used by the forward histogram vector (fForward).
+ * Base class destructor handles cleanup of theory objects and other shared resources.
+ */
virtual ~PRunMuMinus();
+ /**
+ * \brief Calculates χ² between μ⁻ data and theory (least-squares fit metric).
+ *
+ * Computes the chi-squared statistic for the forward histogram:
+ * \f[ \chi^2 = \sum_{i={\rm start}}^{\rm end} \frac{(N_i^{\rm data} - N_i^{\rm theory})^2}{\sigma_i^2} \f]
+ *
+ * where:
+ * - N_i^data is the background-corrected, packed histogram count in bin i
+ * - N_i^theory is the theory prediction from THEORY block
+ * - σ_i is the propagated error (including background uncertainty)
+ * - Sum runs from fStartTimeBin to fEndTimeBin (fit range)
+ *
+ * Implementation uses OpenMP parallelization when available for performance.
+ * Each thread processes a chunk of bins independently, with reduction for final sum.
+ *
+ * \param par Parameter vector from MINUIT with current parameter values
+ * \return Chi-squared value (lower is better; minimize during fitting)
+ *
+ * \see CalcMaxLikelihood() for alternative fit metric (better for low counts)
+ * \see CalcTheory() for theory function evaluation details
+ */
virtual Double_t CalcChiSquare(const std::vector