372 lines
17 KiB
C++
372 lines
17 KiB
C++
/***************************************************************************
|
|
|
|
PRunBase.h
|
|
|
|
Author: Andreas Suter
|
|
e-mail: andreas.suter@psi.ch
|
|
|
|
***************************************************************************/
|
|
|
|
/***************************************************************************
|
|
* Copyright (C) 2007-2026 by Andreas Suter *
|
|
* andreas.suter@psi.ch *
|
|
* *
|
|
* This program is free software; you can redistribute it and/or modify *
|
|
* it under the terms of the GNU General Public License as published by *
|
|
* the Free Software Foundation; either version 2 of the License, or *
|
|
* (at your option) any later version. *
|
|
* *
|
|
* This program is distributed in the hope that it will be useful, *
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
|
* GNU General Public License for more details. *
|
|
* *
|
|
* You should have received a copy of the GNU General Public License *
|
|
* along with this program; if not, write to the *
|
|
* Free Software Foundation, Inc., *
|
|
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
|
***************************************************************************/
|
|
|
|
#ifndef _PRUNBASE_H_
|
|
#define _PRUNBASE_H_
|
|
|
|
#include <vector>
|
|
#include <memory>
|
|
|
|
#include <TString.h>
|
|
|
|
#include "PMusr.h"
|
|
#include "PMsrHandler.h"
|
|
#include "PRunDataHandler.h"
|
|
#include "PTheory.h"
|
|
|
|
//------------------------------------------------------------------------------------------
|
|
/**
|
|
* \brief Abstract base class defining the interface for all μSR fit types.
|
|
*
|
|
* PRunBase establishes a common API for processing and fitting different types of μSR data
|
|
* (single histogram, asymmetry, RRF, etc.). This class serves as the foundation of the
|
|
* musrfit framework, providing core functionality and requiring derived classes to implement
|
|
* fit-type-specific algorithms.
|
|
*
|
|
* \section derived_classes Derived Classes
|
|
* Each derived class handles a specific type of μSR measurement:
|
|
* - <b>PRunSingleHisto:</b> Single detector histogram fits (basic time-differential μSR)
|
|
* - <b>PRunAsymmetry:</b> Asymmetry fits: \f$ A(t) = \frac{F(t) - \alpha B(t)}{F(t) + \alpha B(t)} \f$
|
|
* - <b>PRunSingleHistoRRF:</b> Single histogram in rotating reference frame (high-TF analysis)
|
|
* - <b>PRunAsymmetryRRF:</b> Asymmetry in rotating reference frame
|
|
* - <b>PRunAsymmetryBNMR:</b> β-NMR asymmetry (helicity-dependent measurements)
|
|
* - <b>PRunMuMinus:</b> Negative muon fits (different decay properties)
|
|
* - <b>PRunNonMusr:</b> General x-y data fits (non-μSR time series data)
|
|
*
|
|
* \section key_responsibilities Key Responsibilities
|
|
* - Loading raw histogram data from various file formats (ROOT, NeXus, WKM, etc.)
|
|
* - Extracting metadata (magnetic field, temperature, energy) from data files
|
|
* - Managing time-zero (t0) determination and validation
|
|
* - Background subtraction (fixed or estimated from pre-t0 region)
|
|
* - Time bin packing/rebinning for improved statistics
|
|
* - Theory function evaluation via PTheory interface
|
|
* - χ² or maximum likelihood calculation for fitting
|
|
* - RRF transformations with Kaiser filtering (for RRF-derived classes)
|
|
* - Run addition (combining multiple identical measurements)
|
|
* - Histogram grouping (combining multiple detectors)
|
|
*
|
|
* \section workflow Processing Workflow
|
|
* The typical processing sequence for a fit is:
|
|
* 1. <b>Construction:</b> Initialize from MSR file and raw data handler
|
|
* 2. <b>PrepareData():</b> Load and preprocess raw data (implemented by derived classes)
|
|
* - Load histograms, validate t0 values, subtract background
|
|
* - Calculate asymmetry or apply RRF transformation (if applicable)
|
|
* - Pack bins, propagate errors
|
|
* 3. <b>Fitting loop</b> (called by MINUIT):
|
|
* - <b>CalcTheory():</b> Evaluate theory function at data points
|
|
* - <b>CalcChiSquare():</b> Compute χ² between data and theory
|
|
* - MINUIT adjusts parameters to minimize χ²
|
|
* 4. <b>Visualization:</b> Theory calculated with higher resolution for plotting
|
|
*
|
|
* \section design_pattern Design Pattern
|
|
* PRunBase implements the <b>Template Method</b> design pattern:
|
|
* - Base class defines the overall workflow and common operations
|
|
* - Derived classes implement fit-type-specific algorithms (pure virtual methods)
|
|
* - Ensures consistency across different fit types while allowing specialization
|
|
*
|
|
* \section thread_safety Thread Safety
|
|
* PRunBase objects are NOT thread-safe. Each thread in parallel fitting should
|
|
* create its own PRunBase-derived object. Theory evaluation may use OpenMP
|
|
* internally (see PTheory documentation).
|
|
*
|
|
* \see PTheory for theory function evaluation
|
|
* \see PMsrHandler for MSR file parsing
|
|
* \see PRunDataHandler for raw data loading
|
|
*/
|
|
class PRunBase
|
|
{
|
|
public:
|
|
/// Default constructor
|
|
PRunBase();
|
|
|
|
/**
|
|
* \brief Constructor initializing run from MSR file and raw data.
|
|
*
|
|
* Initializes the run object by:
|
|
* - Storing pointers to MSR file handler and raw data handler
|
|
* - Extracting run-specific settings from MSR RUN block
|
|
* - Validating function parameter mappings
|
|
* - Creating PTheory object for theory evaluation
|
|
* - Initializing metadata and function value vectors
|
|
*
|
|
* \param msrInfo Pointer to MSR file handler containing all fit settings
|
|
* \param rawData Pointer to raw data handler for accessing histogram data
|
|
* \param runNo Run number (0-based index in MSR file RUN blocks)
|
|
* \param tag Operation mode: kFit (fitting), kView (display only)
|
|
*/
|
|
PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag);
|
|
|
|
/**
|
|
* \brief Virtual destructor.
|
|
*
|
|
* Cleans up allocated resources including:
|
|
* - t0 value vectors
|
|
* - Function value vectors
|
|
* - PTheory object (via unique_ptr)
|
|
*/
|
|
virtual ~PRunBase();
|
|
|
|
/**
|
|
* \brief Calculates χ² between data and theory (pure virtual).
|
|
*
|
|
* Computes the chi-squared statistic:
|
|
* \f[ \chi^2 = \sum_{i=1}^{N} \frac{(y_i^{\rm data} - y_i^{\rm theory})^2}{\sigma_i^2} \f]
|
|
*
|
|
* This is the standard least-squares metric for fitting. It assumes Gaussian
|
|
* statistics (valid for high-count data). For low-count data, consider using
|
|
* CalcMaxLikelihood() instead.
|
|
*
|
|
* Called by MINUIT during each fit iteration to evaluate parameter quality.
|
|
* Derived classes implement the specific calculation for their data type.
|
|
*
|
|
* \param par Vector of fit parameter values from MINUIT
|
|
* \return Chi-squared value (lower is better)
|
|
*
|
|
* \see CalcMaxLikelihood for alternative fit metric
|
|
*/
|
|
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par) = 0;
|
|
|
|
/**
|
|
* \brief Calculates expected chi-square for statistical analysis (pure virtual).
|
|
*
|
|
* Computes the expected χ² value based on the theory and error estimates,
|
|
* useful for evaluating the quality of error bars and fit statistics.
|
|
*
|
|
* For properly estimated errors: χ²_expected ≈ number of degrees of freedom
|
|
*
|
|
* \param par Vector of fit parameter values from MINUIT
|
|
* \return Expected chi-squared value
|
|
*
|
|
* \note Implementation is optional in derived classes; many return 0.0
|
|
*/
|
|
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par) = 0;
|
|
|
|
/**
|
|
* \brief Calculates maximum likelihood estimator (pure virtual).
|
|
*
|
|
* Computes the negative log-likelihood for Poisson statistics:
|
|
* \f[ -2\ln L = 2\sum_{i=1}^{N} \left[y_i^{\rm theory} - y_i^{\rm data} \ln(y_i^{\rm theory})\right] \f]
|
|
*
|
|
* Maximum likelihood estimation is superior to χ² for low-count data where
|
|
* the Gaussian approximation breaks down (typically when counts < 10-20 per bin).
|
|
* It naturally handles Poisson statistics without requiring error estimates.
|
|
*
|
|
* Called by MINUIT as an alternative to χ² minimization. MINUIT minimizes
|
|
* this function just like χ².
|
|
*
|
|
* \param par Vector of fit parameter values from MINUIT
|
|
* \return Negative 2 times log-likelihood (lower is better)
|
|
*
|
|
* \see CalcChiSquare for standard least-squares metric
|
|
* \note Not implemented in all derived classes
|
|
*/
|
|
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par) = 0;
|
|
|
|
/**
|
|
* \brief Sets the fit time range for this run.
|
|
*
|
|
* Updates the fitting window, useful for the FIT_RANGE command which allows
|
|
* scanning different time windows to find the optimal range for parameter
|
|
* extraction. Can be called multiple times during a fit sequence.
|
|
*
|
|
* The fit range is specified in microseconds (μs) from t0. Multiple ranges
|
|
* can be specified for different runs in a global fit.
|
|
*
|
|
* \param fitRange Vector of (start, end) time pairs in microseconds
|
|
*
|
|
* Example: fitRange[0] = (0.1, 10.0) means fit from 0.1 μs to 10.0 μs after t0
|
|
*/
|
|
virtual void SetFitRange(PDoublePairVector fitRange);
|
|
|
|
/**
|
|
* \brief Evaluates theory function at all data points (pure virtual).
|
|
*
|
|
* Calculates the expected signal at each time point using the current parameter
|
|
* values from the MSR THEORY block. This is called:
|
|
* - During each MINUIT fit iteration
|
|
* - After fitting for visualization
|
|
* - When evaluating functions from FUNCTIONS block
|
|
*
|
|
* The theory values are stored in fData for comparison with measured data.
|
|
* Derived classes implement fit-type-specific theory calculation (e.g.,
|
|
* single histogram vs. asymmetry).
|
|
*
|
|
* \see PTheory for the underlying theory evaluation engine
|
|
*/
|
|
virtual void CalcTheory() = 0;
|
|
|
|
virtual Double_t GetTimeResolution() { return fTimeResolution; } ///< returns the native raw data time resolution
|
|
|
|
/**
|
|
* \brief Returns the run number (0-based index in MSR file).
|
|
* \return Run number corresponding to position in MSR RUN blocks
|
|
*/
|
|
virtual UInt_t GetRunNo() { return fRunNo; }
|
|
|
|
/**
|
|
* \brief Returns pointer to processed data container.
|
|
*
|
|
* The PRunData object contains:
|
|
* - Background-corrected histogram data
|
|
* - Packed/rebinned data points
|
|
* - Error bars
|
|
* - Time grid information
|
|
* - Theory values (after CalcTheory() is called)
|
|
*
|
|
* \return Pointer to PRunData object with processed data
|
|
*/
|
|
virtual PRunData* GetData() { return &fData; }
|
|
|
|
virtual PMetaData GetMetaData() { return fMetaData; } ///< returns the meta data
|
|
|
|
/**
|
|
* \brief Cleans up internal data structures.
|
|
*
|
|
* Releases memory used by temporary data structures. Called when
|
|
* run processing is complete or when resetting for a new fit.
|
|
*/
|
|
virtual void CleanUp();
|
|
|
|
/**
|
|
* \brief Returns validity status of this run object.
|
|
*
|
|
* A run becomes invalid if:
|
|
* - Required data files cannot be loaded
|
|
* - MSR file settings are inconsistent
|
|
* - Theory initialization fails
|
|
* - Data preprocessing encounters errors
|
|
*
|
|
* \return True if run initialized successfully, false otherwise
|
|
*/
|
|
virtual Bool_t IsValid() { return fValid; }
|
|
|
|
protected:
|
|
Bool_t fValid; ///< Flag indicating if run object initialized successfully; false if any error occurred
|
|
|
|
EPMusrHandleTag fHandleTag; ///< Operation mode: kFit (fitting), kView (display only), kEmpty (uninitialized)
|
|
|
|
Int_t fRunNo; ///< Run number (0-based index in MSR file RUN blocks)
|
|
PMsrHandler *fMsrInfo; ///< Pointer to MSR file handler (owned externally, not deleted here)
|
|
PMsrRunBlock *fRunInfo; ///< Pointer to this run's RUN block settings within fMsrInfo
|
|
PRunDataHandler *fRawData; ///< Pointer to raw data handler (owned externally, not deleted here)
|
|
|
|
PRunData fData; ///< Processed data container: background-corrected, packed, with theory values
|
|
Double_t fTimeResolution; ///< Time resolution of raw histogram data in microseconds (μs), e.g., 0.01953125 μs for PSI GPS
|
|
PMetaData fMetaData; ///< Experimental metadata extracted from data file header (magnetic field, temperature, beam energy)
|
|
PDoubleVector fT0s; ///< Time-zero bin values for all histograms in this run (forward, backward, etc.)
|
|
std::vector<PDoubleVector> fAddT0s; ///< Time-zero bin values for additional runs to be added to main run
|
|
|
|
Double_t fFitStartTime; ///< Fit range start time in microseconds (μs) relative to t0
|
|
Double_t fFitEndTime; ///< Fit range end time in microseconds (μs) relative to t0
|
|
|
|
PDoubleVector fFuncValues; ///< Cached values of user-defined functions from FUNCTIONS block, evaluated at current parameters
|
|
std::unique_ptr<PTheory> fTheory; ///< Theory function evaluator (smart pointer, automatically deleted)
|
|
|
|
PDoubleVector fKaiserFilter; ///< Kaiser window FIR filter coefficients for smoothing RRF theory curves
|
|
|
|
/**
|
|
* \brief Prepares raw data for fitting (pure virtual).
|
|
*
|
|
* This is the main data preprocessing pipeline, implemented differently by each
|
|
* derived class according to the fit type. Common operations include:
|
|
* - Loading histogram data from raw data files
|
|
* - Validating and determining t0 values
|
|
* - Subtracting background (fixed or estimated from pre-t0 region)
|
|
* - Calculating asymmetry (for asymmetry-based fits)
|
|
* - Applying RRF transformation (for RRF fits)
|
|
* - Time bin packing/rebinning to improve statistics
|
|
* - Proper error propagation through all transformations
|
|
* - Adding multiple runs together (if specified)
|
|
* - Grouping detector histograms (if specified)
|
|
*
|
|
* Called once during object construction. If this returns false, the run
|
|
* object is marked as invalid.
|
|
*
|
|
* \return True on success, false if any preprocessing step fails
|
|
*
|
|
* \see PRunAsymmetry::PrepareData() for asymmetry-specific implementation
|
|
* \see PRunSingleHisto::PrepareData() for single histogram implementation
|
|
*/
|
|
virtual Bool_t PrepareData() = 0;
|
|
|
|
/**
|
|
* \brief carry out dead time correction
|
|
*
|
|
* \param histos histograms to be corrected
|
|
* \param histoNo histogram numbers
|
|
*/
|
|
virtual void DeadTimeCorrection(std::vector<PDoubleVector> &histos, PUIntVector &histoNo);
|
|
|
|
/**
|
|
* \brief Calculates Kaiser window FIR filter coefficients for RRF smoothing.
|
|
*
|
|
* Computes a Kaiser window finite impulse response (FIR) filter for smoothing
|
|
* theory curves in rotating reference frame (RRF) fits. The Kaiser window
|
|
* provides excellent control over the trade-off between main lobe width
|
|
* (frequency resolution) and side lobe attenuation (spectral leakage).
|
|
*
|
|
* The filter design uses the Kaiser-Bessel formula:
|
|
* \f[ w[n] = \frac{I_0\left(\beta\sqrt{1-\left(\frac{n-\alpha}{\alpha}\right)^2}\right)}{I_0(\beta)} \f]
|
|
*
|
|
* where \f$ I_0 \f$ is the modified Bessel function of the first kind,
|
|
* \f$ \alpha = (M-1)/2 \f$, and \f$ \beta \f$ is determined from the attenuation A.
|
|
*
|
|
* \param wc Cutoff frequency (normalized, 0 to π rad/sample)
|
|
* \param A Attenuation in dB (typical: 60 dB for good side lobe suppression)
|
|
* \param dw Transition width (normalized, typical: 0.1-0.3)
|
|
*
|
|
* Coefficients are stored in fKaiserFilter for use by FilterTheo().
|
|
*
|
|
* \see FilterTheo() for application of the filter
|
|
*/
|
|
virtual void CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw);
|
|
|
|
/**
|
|
* \brief Applies Kaiser FIR filter to theory values for RRF fits.
|
|
*
|
|
* Filters the theory function stored in fData using the Kaiser window
|
|
* coefficients from fKaiserFilter. This ensures the theory curve is smoothed
|
|
* in the same way as the RRF-transformed data, enabling fair comparison
|
|
* in χ² calculation.
|
|
*
|
|
* The filtering is performed via convolution in the time domain:
|
|
* \f[ y_{\rm filtered}[n] = \sum_{k=0}^{M-1} h[k] \cdot y_{\rm theory}[n-k] \f]
|
|
*
|
|
* where h[k] are the Kaiser filter coefficients and M is the filter length.
|
|
*
|
|
* Only used in RRF-derived classes (PRunAsymmetryRRF, PRunSingleHistoRRF).
|
|
*
|
|
* \pre CalculateKaiserFilterCoeff() must be called first to initialize fKaiserFilter
|
|
* \pre fData must contain theory values (CalcTheory() must have been called)
|
|
*/
|
|
virtual void FilterTheo();
|
|
};
|
|
|
|
#endif // _PRUNBASE_H_
|