Files
musrfit/src/include/PFitterFcn.h

172 lines
7.7 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
/***************************************************************************
PFitterFcn.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 _PFITTERFCN_H_
#define _PFITTERFCN_H_
#include <vector>
#include "Minuit2/FCNBase.h"
#include "PRunListCollection.h"
//--------------------------------------------------------------------------
/**
* @brief Objective function interface for ROOT Minuit2 minimization.
*
* This class implements the FCNBase interface required by ROOT's Minuit2
* minimizer. It provides the objective function (χ² or log-likelihood)
* that Minuit2 minimizes during parameter optimization.
*
* The class serves as a bridge between musrfit's data structures
* (PRunListCollection) and Minuit2's optimization algorithms, calculating
* the goodness-of-fit measure for any given parameter set.
*
* @par Fitting modes:
* - **χ² minimization:** Standard least-squares fitting for Gaussian errors
* - **Maximum likelihood:** Poisson statistics, better for low-count data
*
* @par Usage in fitting workflow:
* 1. PFitter creates a PFitterFcn instance with data and fit mode
* 2. Minuit2 calls operator()() repeatedly with trial parameter sets
* 3. operator()() calculates χ²/maxLH by evaluating theory vs. data
* 4. Minuit2 searches parameter space to minimize the returned value
* 5. Up() defines the error criterion (Δχ²=1 or ΔmaxLH=0.5 for 1σ)
*
* @see PFitter, PRunListCollection
* @see ROOT::Minuit2::FCNBase in ROOT Minuit2 documentation
*/
class PFitterFcn : public ROOT::Minuit2::FCNBase
{
public:
//----------------------------------------------------------------------
/**
* @brief Constructor for objective function.
*
* Initializes the function evaluator with preprocessed data and
* configures the error definition based on the fitting mode.
*
* @param runList Pointer to collection of preprocessed run data
* @param useChi2 If true, use χ² minimization; if false, use maximum likelihood
*
* @note The runList pointer must remain valid for the lifetime of this object.
*/
PFitterFcn(PRunListCollection *runList, Bool_t useChi2);
/**
* @brief Destructor.
*/
~PFitterFcn();
//----------------------------------------------------------------------
/**
* @brief Returns error definition for Minuit2 (Up value).
*
* The "Up" value defines what change in the objective function
* corresponds to 1σ error bars on parameters:
* - For χ² fits: Up = 1.0 (parabolic errors, Δχ²=1)
* - For max likelihood: Up = 0.5 (asymmetric errors, ΔmaxLH=0.5)
*
* This value is used by Minuit2's error analysis algorithms (HESSE, MINOS).
*
* @return Error definition value (1.0 for χ², 0.5 for likelihood)
*
* @see ROOT::Minuit2::FCNBase::Up() in Minuit2 manual
*/
Double_t Up() const { return fUp; }
//----------------------------------------------------------------------
/**
* @brief Evaluates objective function for given parameters.
*
* This is the core function called by Minuit2 during minimization.
* It computes either χ² or negative log-likelihood by:
* 1. Passing parameters to PRunListCollection
* 2. Calculating theory predictions for all runs
* 3. Comparing theory vs. data across all fitted bins
* 4. Returning the total χ²/maxLH value
*
* @param par Parameter vector with current trial values
* @return χ² value (if fUseChi2=true) or -2×log-likelihood (if fUseChi2=false)
*
* @note This function must be const as required by FCNBase interface.
* @note For likelihood fits, returns -2×ln(L) so minimization is equivalent to maximizing L.
*
* @par Performance:
* This function is called hundreds to thousands of times during
* a fit, so it's optimized for speed (parallel evaluation if OpenMP enabled).
*/
Double_t operator()(const std::vector<Double_t> &par) const;
//----------------------------------------------------------------------
/**
* @brief Returns total number of bins used in the fit across all runs.
*
* @return Total count of fitted bins (summed over all runs)
*/
UInt_t GetTotalNoOfFittedBins() { return fRunListCollection->GetTotalNoOfBinsFitted(); }
//----------------------------------------------------------------------
/**
* @brief Returns number of fitted bins for a specific run.
*
* @param idx Run index (0-based)
* @return Number of bins fitted for the specified run
*/
UInt_t GetNoOfFittedBins(const UInt_t idx) { return fRunListCollection->GetNoOfBinsFitted(idx); }
//----------------------------------------------------------------------
/**
* @brief Calculates expected χ² (or maxLH) for quality assessment.
*
* Computes the theoretical expected value of χ² assuming the model
* is correct. This is used to assess goodness-of-fit:
* - If observed χ² ≈ expected χ²: fit is consistent with data quality
* - If observed χ² >> expected χ²: systematic deviations present
* - If observed χ² << expected χ²: possible overestimated errors
*
* For single histogram fits, expected χ² = NDF. For asymmetry fits,
* the calculation is more complex due to error propagation.
*
* @param par Parameter vector for evaluation
* @param totalExpectedChisq Returns total expected χ²/maxLH (output)
* @param expectedChisqPerRun Returns expected χ²/maxLH for each run (output)
*
* @note The expectedChisqPerRun vector is resized and filled by this method.
*/
void CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_t &totalExpectedChisq, std::vector<Double_t> &expectedChisqPerRun);
private:
Double_t fUp; ///< Error definition: 1.0 for χ² (1σ = Δχ²=1), 0.5 for maxLH (1σ = ΔmaxLH=0.5)
Bool_t fUseChi2; ///< Fit mode flag: true = χ² minimization, false = max log-likelihood
PRunListCollection *fRunListCollection; ///< Pointer to preprocessed muSR data collection
};
#endif // _PFITTERFCN_H_