LMU/musrfit
5
0
Fork 0
Code Issues 5 Pull Requests Actions Packages Projects Releases 4 Activity

improve the doxygen docu of PRunMuMinus.*

Browse Source
This commit is contained in:
suter_a
2025-11-15 08:35:24 +01:00
parent 278fd47e52
commit 3e545f1d50
2 changed files with 574 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

422
src/include/PRunMuMinus.h
View File

@@ -33,49 +33,455 @@
#include "PRunBase.h"
/**
* <p>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
*
* <b>Key differences from μ⁺:</b>
* - <b>Charge:</b> Negatively charged (antimuon has positive charge)
* - <b>Behavior in matter:</b> μ⁻ are captured by nuclei rather than stopping interstitially
* - <b>Cascade process:</b> μ⁻ cascade down atomic orbitals before nuclear capture
* - <b>Lifetime:</b> Effective lifetime is reduced by nuclear capture (typically ~100-1000 ns vs. 2.2 μs free lifetime)
* - <b>Decay signature:</b> Different detector response due to capture vs. decay
* - <b>Applications:</b> 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. <b>Data Loading:</b> Load single histogram from data file
* 2. <b>Background Subtraction:</b> Remove constant or estimated background
* 3. <b>Time-Zero Determination:</b> Identify μ⁻ arrival time (t0)
* 4. <b>Bin Packing:</b> Rebin data if specified (required parameter)
* 5. <b>Theory Calculation:</b> Evaluate capture-modified exponential decay
* 6. <b>Fit Metric:</b> χ² or maximum likelihood minimization
*
* \section muminus_theory Typical Theory Functions
*
* Common models for μ⁻ data include:
* - <b>Simple exponential:</b> N(t) = N₀ · exp(-t/τ_eff)
* - <b>With relaxation:</b> N(t) = N₀ · exp(-t/τ_eff) · [1 + A·cos(ωt + φ)]
* - <b>Multi-component:</b> 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<Double_t>& par);
/**
* \brief Calculates expected χ² based on theory predictions (statistical diagnostic).
*
* Computes the expected chi-squared assuming theory values are the "true" counts:
* \f[ \chi^2_{\rm expected} = \sum_{i} \frac{(N_i^{\rm data} - N_i^{\rm theory})^2}{N_i^{\rm theory}} \f]
*
* This diagnostic helps evaluate:
* - Quality of error estimates
* - Poisson vs. Gaussian statistics validity
* - Over/under-dispersion in data
*
* For well-estimated errors: χ²_expected ≈ number of degrees of freedom
*
* \param par Parameter vector from MINUIT
* \return Expected χ² value (currently returns 0.0 - placeholder implementation)
*
* \note Current implementation returns 0.0; full calculation performed but not returned
*/
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
/**
* \brief Calculates negative log-likelihood for Poisson statistics (low-count fit metric).
*
* Computes the maximum likelihood estimator assuming Poisson-distributed counts:
* \f[ -2\ln L = 2\sum_{i} \left[N_i^{\rm theory} - N_i^{\rm data} \ln(N_i^{\rm theory})\right] \f]
*
* Maximum likelihood is superior to χ² when:
* - Count rates are low (< 10-20 counts per bin)
* - Poisson statistics dominate (Gaussian approximation invalid)
* - Error estimation is uncertain
*
* The factor of 2 makes the likelihood comparable to χ² in the Gaussian limit.
* MINUIT minimizes this function just like χ².
*
* Implementation details:
* - Uses OpenMP parallelization for performance
* - Handles edge cases: zero data (skips log term), negative theory (warning + skip)
* - Threshold: data > 10^-9 to include log term
*
* \param par Parameter vector from MINUIT
* \return -2×ln(L) value (lower is better; minimize during fitting)
*
* \see CalcChiSquare() for standard least-squares metric
*/
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
/**
* \brief Evaluates theory function at all data points (or high-resolution grid).
*
* Calculates the expected μ⁻ decay signal using the THEORY block functions.
* The theory is evaluated either:
* - At data point times (fTheoAsData=true): For fitting and exact comparison
* - On high-resolution grid (fTheoAsData=false): For smooth plotting
*
* Theory evaluation:
* 1. Determines time grid based on fTheoAsData flag
* 2. Evaluates FUNCTIONS block (user-defined functions)
* 3. Loops over time points calling fTheory->Func(time, par, funcValues)
* 4. Stores results in fData for χ² or likelihood calculation
*
* The theory typically models:
* - Exponential decay: N₀·exp(-t/τ_eff)
* - With modulation: N₀·exp(-t/τ_eff)·[1 + A·cos(ωt + φ)]
* - Multi-component captures
*
* \see PTheory::Func() for theory function implementation
* \see PrepareData() for data grid setup
*/
virtual void CalcTheory();
/**
* \brief Returns the number of bins included in the fit range.
*
* Calculates and returns the count of histogram bins between fStartTimeBin
* and fEndTimeBin (fit range boundaries). This is used for:
* - Degrees of freedom: ν = N_bins - N_params
* - Reduced χ²: χ²_red = χ² / ν
* - Statistical analysis and diagnostics
*
* Internally calls CalcNoOfFitBins() to update fNoOfFitBins.
*
* \return Number of bins in fit range
*
* \see CalcNoOfFitBins() for bin range calculation details
*/
virtual UInt_t GetNoOfFitBins();
/**
* \brief Sets fit range using bin-offset specification (COMMANDS block syntax).
*
* Parses and applies fit range specified as bin offsets from good bin markers.
* This supports the COMMANDS block FIT_RANGE syntax for dynamic range adjustments.
*
* Format: "fit_range fgb+n0 lgb-n1"
* - fgb = first good bin (typically t0 bin)
* - lgb = last good bin (end of data range)
* - n0, n1 = offsets in bins (can be + or -)
*
* Example: "fit_range fgb+5 lgb-10"
* → Start 5 bins after t0, end 10 bins before last good bin
*
* \param fitRange String with fit range specification in bin offsets
*
* \see SetFitRange(PDoublePairVector) in PRunBase for time-based specification
*/
virtual void SetFitRangeBin(const TString fitRange);
/**
* \brief Returns the first bin index in the fit range.
* \return Start bin index (0-based, after packing applied)
*/
virtual Int_t GetStartTimeBin() { return fStartTimeBin; }
/**
* \brief Returns the last bin index in the fit range (exclusive).
* \return End bin index (loop condition: i < fEndTimeBin)
*/
virtual Int_t GetEndTimeBin() { return fEndTimeBin; }
/**
* \brief Returns the bin packing factor.
* \return Packing value (number of raw bins combined into one packed bin)
*/
virtual Int_t GetPacking() { return fPacking; }
/**
* \brief Calculates start/end bin indices from fit time range.
*
* Converts the fit range specified in time (microseconds from t0) into
* bin indices accounting for:
* - Data time start offset
* - Time resolution (bin width)
* - Packing factor
*
* Updates:
* - fStartTimeBin: First bin in fit range
* - fEndTimeBin: One past last bin (exclusive upper bound)
* - fNoOfFitBins: Count of bins in range
*
* Called automatically when fit range is changed via SetFitRange().
*
* \see SetFitRange() in PRunBase for time-based range setting
*/
virtual void CalcNoOfFitBins();
protected:
/**
* \brief Main data preparation routine for μ⁻ fitting and viewing.
*
* Orchestrates the complete data preprocessing pipeline:
* 1. Determines operation mode (fitting vs. viewing)
* 2. Calls PrepareFitData() or PrepareRawViewData() accordingly
* 3. Validates successful data loading and processing
*
* This method is called once during object construction. If it returns false,
* the run object is marked as invalid (fValid=false).
*
* \return True if data preparation succeeds, false on any error
*
* \see PrepareFitData() for fitting mode preprocessing
* \see PrepareRawViewData() for viewing mode preprocessing
*/
virtual Bool_t PrepareData();
/**
* \brief Prepares μ⁻ histogram data for fitting.
*
* Performs comprehensive data preprocessing:
* 1. Loads forward histogram from data file
* 2. Extracts metadata (field, energy, temperature)
* 3. Validates and determines t0 value
* 4. Subtracts background (fixed or estimated)
* 5. Packs bins according to packing factor
* 6. Propagates errors through all transformations
* 7. Sets up time grid for theory evaluation
* 8. Establishes data/fit ranges
*
* Background correction (if specified in RUN block):
* - Fixed background: Subtract constant value
* - Estimated background: Calculate from pre-t0 bins, subtract with error propagation
*
* \param runData Pointer to raw run data handler containing histogram
* \param histoNo Histogram index in data file (from "forward" entry in RUN block)
* \return True on success, false if data loading or preprocessing fails
*
* \see PrepareRawViewData() for alternative viewing mode processing
*/
virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo);
/**
* \brief Prepares μ⁻ histogram data for viewing/plotting (minimal processing).
*
* Lighter-weight data preprocessing for visualization:
* 1. Loads forward histogram from data file
* 2. Extracts metadata
* 3. Determines t0
* 4. Sets up data range (may extend beyond fit range)
* 5. Minimal or no background subtraction
* 6. Optional packing for display
*
* Used when tag=kView to prepare data for plotting without full fitting infrastructure.
* Typically displays wider time range than fit range for context.
*
* \param runData Pointer to raw run data handler containing histogram
* \param histoNo Histogram index in data file
* \return True on success, false if data loading fails
*
* \see PrepareFitData() for full fitting mode processing
*/
virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo);
private:
UInt_t fNoOfFitBins; ///< number of bins to be fitted
Int_t fPacking; ///< packing for this particular run. Either given in the RUN- or GLOBAL-block.
Bool_t fTheoAsData; ///< true=only calculate the theory points at the data points, false=calculate more points for the theory as compared to data are calculated which lead to 'nicer' Fouriers
UInt_t fNoOfFitBins; ///< Number of bins within fit range (between fStartTimeBin and fEndTimeBin)
Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb
/**
* \brief Bin packing factor (REQUIRED for μ⁻).
*
* Number of consecutive raw histogram bins combined into one packed bin.
* Higher packing improves statistics but reduces time resolution.
*
* Source priority:
* 1. RUN block "packing" entry
* 2. GLOBAL block "packing" entry
* 3. ERROR if neither specified (mandatory for μ⁻)
*
* Typical values: 1 (no packing), 5, 10, 20
*/
Int_t fPacking;
PDoubleVector fForward; ///< forward histo data
/**
* \brief Theory calculation mode flag.
*
* Controls theory grid resolution:
* - true: Theory evaluated only at data point times (efficient for fitting)
* - false: Theory evaluated on finer grid (smooth curves for plotting/Fourier)
*
* Set from PRunListCollection based on PLOT block settings.
*/
Bool_t fTheoAsData;
Int_t fStartTimeBin; ///< bin at which the fit starts
Int_t fEndTimeBin; ///< bin at which the fit ends
/**
* \brief Good bin markers for bin-based fit range specification.
*
* Stores reference bins used in COMMANDS block FIT_RANGE parsing:
* - fGoodBins[0]: First good bin (fgb) - typically t0 bin
* - fGoodBins[1]: Last good bin (lgb) - end of valid data range
*
* Used when fit range is specified as "fgb+n0 lgb-n1" rather than absolute times.
* Needed because COMMANDS block can dynamically change fit range during execution.
*/
Int_t fGoodBins[2];
/**
* \brief Forward detector histogram data (background-corrected, packed).
*
* Contains the processed μ⁻ decay histogram after:
* - Background subtraction
* - Bin packing
* - Conversion to appropriate units
*
* This is the data fitted against theory predictions.
* Stored in fData.GetValue() vector for χ² calculation.
*/
PDoubleVector fForward;
Int_t fStartTimeBin; ///< First bin index in fit range (inclusive, 0-based after packing)
Int_t fEndTimeBin; ///< Last bin index in fit range (exclusive: loop as i < fEndTimeBin)
/**
* \brief Determines and validates t0 values for μ⁻ histogram.
*
* Extracts time-zero (muon arrival time) from:
* 1. RUN block "t0" entry (if specified)
* 2. Data file header (if available)
* 3. GLOBAL block default (fallback)
* 4. Automatic determination from histogram structure (if enabled)
*
* For μ⁻, t0 marks the start of the decay/capture signal.
* Validates that t0 is within reasonable bounds of the histogram.
*
* \param runData Raw run data containing histogram and metadata
* \param globalBlock GLOBAL block from MSR file with default settings
* \param histoNo Vector of histogram indices to process
* \return True if valid t0 found/determined, false otherwise
*/
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);
/**
* \brief Determines data range (region of valid histogram data).
*
* Establishes the data range boundaries from RUN block specification:
* - Start bin: First bin with valid μ⁻ data (typically after t0)
* - End bin: Last bin before noise/artifacts dominate
*
* Data range is typically wider than fit range. The fit range is a
* subset of the data range optimized for parameter extraction.
*
* \return True if valid data range determined, false if specification invalid
*
* \see GetProperFitRange() for fit range (subset of data range)
*/
virtual Bool_t GetProperDataRange();
/**
* \brief Determines fit range from MSR file settings.
*
* Extracts fit range boundaries (time window for χ² calculation) from:
* 1. RUN block "fit" entry (if specified)
* 2. GLOBAL block "fit" entry (fallback)
*
* Fit range format in MSR file:
* - Time-based: "fit 0.1 10.0" (μs from t0)
* - Bin-based: "fit fgb+5 lgb-10" (offsets from markers)
*
* The fit range determines which bins contribute to χ² or likelihood.
* Choosing the optimal range is critical for accurate parameter extraction.
*
* \param globalBlock GLOBAL block from MSR file with default fit settings
*
* \see CalcNoOfFitBins() for converting time range to bin indices
*/
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);
};