Files
musrfit/src/include/PRunAsymmetry.h
T
2026-01-25 07:45:44 +01:00

285 lines
12 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.
/***************************************************************************
PRunAsymmetry.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 _PRUNASYMMETRY_H_
#define _PRUNASYMMETRY_H_
#include "PRunBase.h"
//---------------------------------------------------------------------------
/**
* \brief Class for handling standard μSR asymmetry fits.
*
* PRunAsymmetry implements asymmetry fitting for conventional μSR experiments,
* where the asymmetry is calculated from forward and backward detector histograms:
*
* \f[ A(t) = \frac{F(t) - \alpha B(t)}{F(t) + \alpha B(t)} \f]
*
* The class supports various configurations:
* - α and β parameters (detector efficiency corrections)
* - Background subtraction (fixed or estimated)
* - Bin packing for improved statistics
* - Rotating reference frame (RRF) data handling
* - Multiple histogram grouping and run adding
*
* The α/β parameter configurations are:
* - Tag 1: α = β = 1 (no corrections)
* - Tag 2: α ≠ 1, β = 1 (forward/backward asymmetry correction)
* - Tag 3: α = 1, β ≠ 1 (alternative correction)
* - Tag 4: α ≠ 1, β ≠ 1 (both corrections active)
*
* \see PRunBase for the base class providing common functionality
* \see PRunAsymmetryBNMR for β-NMR asymmetry (helicity-dependent)
* \see PRunSingleHisto for single histogram fits
*/
class PRunAsymmetry : public PRunBase
{
public:
/// Default constructor
PRunAsymmetry();
/**
* \brief Main constructor for μSR asymmetry fitting.
* \param msrInfo Pointer to MSR file handler
* \param rawData Pointer to raw run data handler
* \param runNo Run number within the MSR file
* \param tag Operation mode (kFit for fitting, kView for viewing)
* \param theoAsData If true, calculate theory only at data points; if false, calculate additional points for Fourier
*/
PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData);
/// Destructor
virtual ~PRunAsymmetry();
/**
* \brief Calculates chi-square for the current parameter set.
* \param par Parameter vector from MINUIT
* \return Chi-square value
*/
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
/**
* \brief Calculates expected chi-square (for statistical analysis).
* \param par Parameter vector from MINUIT
* \return Expected chi-square value
*/
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
/**
* \brief Calculates maximum likelihood estimator.
* \param par Parameter vector from MINUIT
* \return Maximum likelihood value
*/
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
/**
* \brief Calculates theoretical asymmetry function.
*
* Computes the theory values for the μSR asymmetry based on the
* current parameters and fit function.
*/
virtual void CalcTheory();
/**
* \brief Returns the number of bins used in the fit.
* \return Number of fit bins
*/
virtual UInt_t GetNoOfFitBins();
/**
* \brief Sets the fit range in bins.
* \param fitRange Fit range string (format depends on configuration)
*/
virtual void SetFitRangeBin(const TString fitRange);
/**
* \brief Returns the first bin used in the fit.
* \return Start time bin index
*/
virtual Int_t GetStartTimeBin() { return fStartTimeBin; }
/**
* \brief Returns the last bin used in the fit.
* \return End time bin index
*/
virtual Int_t GetEndTimeBin() { return fEndTimeBin; }
/**
* \brief Returns the packing factor.
* \return Number of bins combined (1 = no packing)
*/
virtual Int_t GetPacking() { return fPacking; }
/**
* \brief Calculates the number of bins to be fitted.
*
* Determines fNoOfFitBins based on the fit range and data availability.
*/
virtual void CalcNoOfFitBins();
protected:
/**
* \brief Prepares all data for fitting or viewing.
* \return True on success, false on error
*
* Main data preparation routine that handles background subtraction,
* packing, and asymmetry calculation from forward/backward histograms.
*/
virtual Bool_t PrepareData();
/**
* \brief Prepares data specifically for fitting.
* \return True on success, false on error
*
* Sets up data structures for the fitting process, including determining
* fit ranges and calculating the number of fit bins.
*/
virtual Bool_t PrepareFitData();
/**
* \brief Prepares data for viewing/plotting.
* \param runData Pointer to raw run data
* \param histoNo Array of histogram numbers [0]=forward, [1]=backward
* \return True on success, false on error
*/
virtual Bool_t PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]);
/**
* \brief Prepares rotating reference frame (RRF) data for viewing.
* \param runData Pointer to raw run data
* \param histoNo Array of histogram numbers [0]=forward, [1]=backward
* \return True on success, false on error
*
* Special handling for RRF data where asymmetry is calculated in a
* rotating coordinate frame.
*/
virtual Bool_t PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]);
private:
UInt_t fAlphaBetaTag; ///< Tag indicating α/β configuration: 1=both unity, 2=α free/β unity, 3=α unity/β free, 4=both free
UInt_t fNoOfFitBins; ///< Number of bins included in the fit
Int_t fPacking; ///< Bin packing factor from RUN or GLOBAL block
Bool_t fTheoAsData; ///< If true, theory calculated only at data points; if false, extra points for nicer Fourier transforms
PDoubleVector fForward; ///< Forward detector histogram data
PDoubleVector fForwardErr; ///< Forward detector histogram errors
PDoubleVector fBackward; ///< Backward detector histogram data
PDoubleVector fBackwardErr; ///< Backward detector histogram errors
Int_t fGoodBins[4]; ///< Good bin boundaries: [0]=forward first, [1]=forward last, [2]=backward first, [3]=backward last
Int_t fStartTimeBin; ///< First bin index for fitting
Int_t fEndTimeBin; ///< Last bin index for fitting
/**
* \brief Subtracts fixed background from histograms.
*
* Subtracts user-specified fixed background values from forward and backward histograms.
* Background values are read from the MSR file (e.g., "backgr.fix 2 3" for forward/backward).
*
* Error propagation:
* \f[ \Delta f_i^{\rm c} = \pm\sqrt{(\Delta f_i)^2 + (\Delta \mathrm{bkg})^2} = \pm\sqrt{f_i + \mathrm{bkg}} \f]
*
* where \f$ f_i^{\rm c} \f$ is the background-corrected histogram, \f$ f_i \f$ is the raw histogram,
* and \f$ \mathrm{bkg} \f$ is the fixed background value.
*
* \return True on success, false if background values are missing
*/
Bool_t SubtractFixBkg();
/**
* \brief Estimates and subtracts background from histograms.
*
* Calculates background from a specified bin range (typically before t0) and subtracts it
* from forward and backward histograms. The background range is adjusted to align with
* accelerator beam cycles when applicable (PSI, RAL, TRIUMF).
*
* Background calculation:
* \f[ \mathrm{bkg} = \frac{1}{N}\sum_{i=0}^N f_i \f]
*
* Error propagation:
* \f[ \Delta f_i^{\rm c} = \pm\sqrt{(\Delta f_i)^2 + (\Delta \mathrm{bkg})^2} = \pm\sqrt{f_i + (\Delta \mathrm{bkg})^2} \f]
*
* where \f$ N \f$ is the number of background bins and
* \f[ \Delta \mathrm{bkg} = \pm\frac{1}{N}\sqrt{\sum_{i=0}^N f_i} \f]
*
* \return True on success, false if background range is out of bounds
*/
Bool_t SubtractEstimatedBkg();
/**
* \brief Retrieves proper t0 values for all histograms.
*
* Determines t0 (time zero) values for forward and backward histograms from:
* 1. Individual RUN block t0 entries
* 2. GLOBAL block t0 values (fallback)
* 3. Data file t0 values (final fallback)
*
* Also handles addT0 corrections for runs to be added together. The t0 values
* are critical for proper alignment of histograms in time.
*
* \param runData Pointer to raw run data containing histogram information
* \param globalBlock Pointer to global MSR block with default t0 values
* \param forwardHisto Vector of forward histogram indices
* \param backwardHistoNo Vector of backward histogram indices
* \return True on success, false if t0 values cannot be determined
*/
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &forwardHisto, PUIntVector &backwardHistoNo);
/**
* \brief Retrieves proper data range for histograms.
*
* Determines the "good bins" range for data analysis from:
* 1. Individual RUN block data range settings
* 2. GLOBAL block data range (fallback)
* 3. Full histogram range (final fallback)
*
* The good bins define which portion of the histograms contains valid data,
* excluding noisy or problematic regions at the start/end.
*
* \param runData Pointer to raw run data
* \param histoNo Array of histogram numbers [0]=forward, [1]=backward
* \return True on success, false on error
*/
virtual Bool_t GetProperDataRange(PRawRunData* runData, UInt_t histoNo[2]);
/**
* \brief Determines the proper fit range from global block.
*
* Extracts fit range settings from the GLOBAL block if not specified in the RUN block.
* The fit range defines the time window used for χ² minimization, specified either in
* time units (μs) or bin numbers.
*
* \param globalBlock Pointer to global MSR block containing default fit range
*/
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);
};
#endif // _PRUNASYMMETRY_H_