improve the doxygen docu of PRunMuMinus.*
This commit is contained in:
@@ -47,7 +47,19 @@
|
|||||||
// Constructor
|
// Constructor
|
||||||
//--------------------------------------------------------------------------
|
//--------------------------------------------------------------------------
|
||||||
/**
|
/**
|
||||||
* <p>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()
|
PRunMuMinus::PRunMuMinus() : PRunBase()
|
||||||
{
|
{
|
||||||
@@ -70,12 +82,43 @@ PRunMuMinus::PRunMuMinus() : PRunBase()
|
|||||||
// Constructor
|
// Constructor
|
||||||
//--------------------------------------------------------------------------
|
//--------------------------------------------------------------------------
|
||||||
/**
|
/**
|
||||||
* <p>Constructor
|
* \brief Main constructor initializing a μ⁻ run from MSR file and raw data.
|
||||||
*
|
*
|
||||||
* \param msrInfo pointer to the msr-file handler
|
* Performs comprehensive initialization for negative muon analysis:
|
||||||
* \param rawData raw run data
|
*
|
||||||
* \param runNo number of the run within the msr-file
|
* 1. <b>Base Class Initialization:</b>
|
||||||
* \param tag tag showing what shall be done: kFit == fitting, kView == viewing
|
* - Calls PRunBase constructor with MSR/data handlers
|
||||||
|
* - Initializes theory engine and parameter mappings
|
||||||
|
*
|
||||||
|
* 2. <b>Packing Validation (CRITICAL for μ⁻):</b>
|
||||||
|
* - 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. <b>Member Initialization:</b>
|
||||||
|
* - 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. <b>Data Preparation:</b>
|
||||||
|
* - 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) :
|
PRunMuMinus::PRunMuMinus(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) :
|
||||||
PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData)
|
PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData)
|
||||||
@@ -114,7 +157,10 @@ PRunMuMinus::PRunMuMinus(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t
|
|||||||
// Destructor
|
// Destructor
|
||||||
//--------------------------------------------------------------------------
|
//--------------------------------------------------------------------------
|
||||||
/**
|
/**
|
||||||
* <p>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()
|
PRunMuMinus::~PRunMuMinus()
|
||||||
{
|
{
|
||||||
@@ -125,12 +171,38 @@ PRunMuMinus::~PRunMuMinus()
|
|||||||
// CalcChiSquare
|
// CalcChiSquare
|
||||||
//--------------------------------------------------------------------------
|
//--------------------------------------------------------------------------
|
||||||
/**
|
/**
|
||||||
* <p>Calculate chi-square.
|
* \brief Calculates χ² between μ⁻ data and theory (least-squares fit metric).
|
||||||
*
|
*
|
||||||
* <b>return:</b>
|
* Computes the chi-squared statistic using the standard formula:
|
||||||
* - chisq value
|
* \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)
|
||||||
|
*
|
||||||
|
* <b>OpenMP Parallelization:</b>
|
||||||
|
* - 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
|
||||||
|
*
|
||||||
|
* <b>Theory Pre-calculation:</b>
|
||||||
|
* 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<Double_t>& par)
|
Double_t PRunMuMinus::CalcChiSquare(const std::vector<Double_t>& par)
|
||||||
{
|
{
|
||||||
@@ -172,12 +244,33 @@ Double_t PRunMuMinus::CalcChiSquare(const std::vector<Double_t>& par)
|
|||||||
// CalcChiSquareExpected (public)
|
// CalcChiSquareExpected (public)
|
||||||
//--------------------------------------------------------------------------
|
//--------------------------------------------------------------------------
|
||||||
/**
|
/**
|
||||||
* <p>Calculate expected chi-square.
|
* \brief Calculates expected χ² assuming theory is the true distribution (diagnostic).
|
||||||
*
|
*
|
||||||
* <b>return:</b>
|
* Computes the expected chi-squared using theory values as the expected counts:
|
||||||
* - chisq value == 0.0
|
* \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)
|
||||||
|
*
|
||||||
|
* <b>OpenMP Parallelization:</b> 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<Double_t>& par)
|
Double_t PRunMuMinus::CalcChiSquareExpected(const std::vector<Double_t>& par)
|
||||||
{
|
{
|
||||||
@@ -221,12 +314,48 @@ Double_t PRunMuMinus::CalcChiSquareExpected(const std::vector<Double_t>& par)
|
|||||||
// CalcMaxLikelihood
|
// CalcMaxLikelihood
|
||||||
//--------------------------------------------------------------------------
|
//--------------------------------------------------------------------------
|
||||||
/**
|
/**
|
||||||
* <p>Calculate log max-likelihood. See http://pdg.lbl.gov/index.html
|
* \brief Calculates negative log-likelihood for Poisson statistics (low-count fit metric).
|
||||||
*
|
*
|
||||||
* <b>return:</b>
|
* Computes the maximum likelihood estimator assuming Poisson-distributed histogram counts:
|
||||||
* - log max-likelihood value
|
* \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)
|
||||||
|
*
|
||||||
|
* <b>When to use likelihood vs. χ²:</b>
|
||||||
|
* - <b>Low counts (< 10-20 per bin):</b> Use likelihood (this method)
|
||||||
|
* - <b>High counts (> 20 per bin):</b> Either method works, χ² is faster
|
||||||
|
* - <b>Uncertain errors:</b> 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 χ²)
|
||||||
|
*
|
||||||
|
* <b>Edge Case Handling:</b>
|
||||||
|
* - 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<Double_t>& par)
|
Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector<Double_t>& par)
|
||||||
{
|
{
|
||||||
@@ -282,9 +411,19 @@ Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector<Double_t>& par)
|
|||||||
// GetNoOfFitBins (public)
|
// GetNoOfFitBins (public)
|
||||||
//--------------------------------------------------------------------------
|
//--------------------------------------------------------------------------
|
||||||
/**
|
/**
|
||||||
* <p>Calculate the number of fitted bins for the current fit range.
|
* \brief Returns the number of bins included in the fit range.
|
||||||
*
|
*
|
||||||
* <b>return:</b> 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()
|
UInt_t PRunMuMinus::GetNoOfFitBins()
|
||||||
{
|
{
|
||||||
|
|||||||
@@ -33,49 +33,455 @@
|
|||||||
#include "PRunBase.h"
|
#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
|
class PRunMuMinus : public PRunBase
|
||||||
{
|
{
|
||||||
public:
|
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();
|
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);
|
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();
|
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);
|
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);
|
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);
|
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();
|
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();
|
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);
|
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; }
|
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; }
|
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; }
|
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();
|
virtual void CalcNoOfFitBins();
|
||||||
|
|
||||||
protected:
|
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();
|
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);
|
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);
|
virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo);
|
||||||
|
|
||||||
private:
|
private:
|
||||||
UInt_t fNoOfFitBins; ///< number of bins to be fitted
|
UInt_t fNoOfFitBins; ///< Number of bins within fit range (between fStartTimeBin and fEndTimeBin)
|
||||||
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
|
|
||||||
|
|
||||||
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);
|
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();
|
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);
|
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|||||||
Reference in New Issue
Block a user