improve the doxygen docu of PRunSingleHistoRRF.*
This commit is contained in:
@@ -53,7 +53,22 @@
|
||||
// Constructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Constructor
|
||||
* \brief Default constructor for RRF single histogram fitting class.
|
||||
*
|
||||
* Initializes all member variables to safe default values:
|
||||
* - fNoOfFitBins = 0 (no bins to fit)
|
||||
* - fBackground = 0 (will be estimated or set from MSR file)
|
||||
* - fBkgErr = 1.0 (default error estimate)
|
||||
* - fRRFPacking = -1 (invalid until set from GLOBAL block)
|
||||
* - fTheoAsData = false (high-resolution theory grid)
|
||||
* - fGoodBins[0,1] = -1 (calculated from data range)
|
||||
* - fN0EstimateEndTime = 1.0 μs (default N₀ estimation window)
|
||||
*
|
||||
* \warning This constructor creates an invalid object that cannot be used
|
||||
* until properly initialized with MSR file data. Use the full
|
||||
* constructor for normal operation.
|
||||
*
|
||||
* \see PRunSingleHistoRRF(PMsrHandler*, PRunDataHandler*, UInt_t, EPMusrHandleTag, Bool_t)
|
||||
*/
|
||||
PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase()
|
||||
{
|
||||
@@ -75,12 +90,44 @@ PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase()
|
||||
// Constructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Constructor
|
||||
* \brief Main constructor for RRF single histogram fitting and viewing.
|
||||
*
|
||||
* \param msrInfo pointer to the msr-file handler
|
||||
* \param rawData raw run data
|
||||
* \param runNo number of the run within the msr-file
|
||||
* \param tag tag showing what shall be done: kFit == fitting, kView == viewing
|
||||
* Constructs a fully initialized RRF single histogram run object by:
|
||||
* -# Validating GLOBAL block presence (mandatory for RRF analysis)
|
||||
* -# Validating RRF frequency specification (rrf_freq in GLOBAL block)
|
||||
* -# Validating RRF packing specification (rrf_packing in GLOBAL block)
|
||||
* -# Calling PrepareData() to load histogram and perform RRF transformation
|
||||
*
|
||||
* <b>GLOBAL Block Requirements:</b>
|
||||
* The RRF fit type requires the following entries in the GLOBAL block:
|
||||
* - \c rrf_freq: Rotation frequency with unit (e.g., "13.554 T", "183.7 MHz")
|
||||
* - \c rrf_packing: Number of bins to average (integer)
|
||||
* - \c rrf_phase: (optional) Initial phase in degrees
|
||||
*
|
||||
* <b>Error Handling:</b>
|
||||
* If any validation fails, the constructor:
|
||||
* - Outputs detailed error message to stderr
|
||||
* - Sets fValid = false
|
||||
* - Returns immediately (PrepareData() is not called)
|
||||
*
|
||||
* \param msrInfo Pointer to MSR file handler (NOT owned, must outlive this object)
|
||||
* \param rawData Pointer to raw run data handler (NOT owned, must outlive this object)
|
||||
* \param runNo Zero-based index of the RUN block in the MSR file
|
||||
* \param tag Operation mode: kFit (fitting) or kView (viewing/plotting)
|
||||
* \param theoAsData If true, theory calculated only at data points (for viewing);
|
||||
* if false, theory uses finer time grid (8× data resolution)
|
||||
*
|
||||
* \warning GLOBAL block with RRF parameters is MANDATORY for this fit type.
|
||||
* Always check IsValid() after construction.
|
||||
*
|
||||
* \note After construction, check IsValid() to ensure initialization succeeded.
|
||||
* Common failure modes:
|
||||
* - Missing GLOBAL block
|
||||
* - Missing rrf_freq specification
|
||||
* - Missing rrf_packing specification
|
||||
* - Data file not found or histogram missing
|
||||
*
|
||||
* \see PrepareData(), PrepareFitData(), PrepareViewData()
|
||||
*/
|
||||
PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) :
|
||||
PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData)
|
||||
@@ -131,7 +178,17 @@ PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *ra
|
||||
// Destructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Destructor
|
||||
* \brief Destructor for RRF single histogram fitting class.
|
||||
*
|
||||
* Cleans up dynamically allocated memory:
|
||||
* - Clears the forward histogram data vector (fForward)
|
||||
* - Other vectors (fM, fMerr, fW, fAerr) are local to PrepareFitData
|
||||
* and cleared automatically
|
||||
*
|
||||
* Base class destructor (PRunBase) handles cleanup of:
|
||||
* - Theory objects
|
||||
* - Function value arrays
|
||||
* - Other shared resources
|
||||
*/
|
||||
PRunSingleHistoRRF::~PRunSingleHistoRRF()
|
||||
{
|
||||
@@ -142,12 +199,42 @@ PRunSingleHistoRRF::~PRunSingleHistoRRF()
|
||||
// CalcChiSquare (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Calculate chi-square.
|
||||
* \brief Calculates χ² between RRF-transformed data and theory (least-squares fit metric).
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - chisq value
|
||||
* Computes the standard chi-square goodness-of-fit statistic for RRF asymmetry:
|
||||
* \f[
|
||||
* \chi^2 = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{\sigma_i^2}
|
||||
* \f]
|
||||
*
|
||||
* \param par parameter vector iterated by minuit2
|
||||
* Unlike standard single histogram fitting, no explicit N₀ or exponential decay
|
||||
* factors appear since the RRF transformation already produces dimensionless
|
||||
* asymmetry with properly propagated errors.
|
||||
*
|
||||
* <b>Algorithm:</b>
|
||||
* -# Evaluate all user-defined functions from FUNCTIONS block
|
||||
* -# Pre-evaluate theory at t=1.0 to initialize any stateful functions
|
||||
* (e.g., LF relaxation, user functions with internal state)
|
||||
* -# Loop over fit range bins [fStartTimeBin, fEndTimeBin)
|
||||
* -# For each bin: calculate time, evaluate theory, accumulate χ²
|
||||
*
|
||||
* <b>OpenMP Parallelization:</b>
|
||||
* When compiled with OpenMP (HAVE_GOMP defined):
|
||||
* - Dynamic scheduling with chunk size = max(10, N_bins / N_processors)
|
||||
* - Private variables per thread: i, time, diff
|
||||
* - Reduction performed on chisq accumulator
|
||||
* - Thread-safe due to pre-evaluation of theory at t=1.0
|
||||
*
|
||||
* \param par Parameter vector from MINUIT minimizer, containing current
|
||||
* estimates of all fit parameters
|
||||
*
|
||||
* \return χ² value (sum over all bins in fit range). Minimize this value
|
||||
* during fitting to find optimal parameters.
|
||||
*
|
||||
* \note The theory function is evaluated in the RRF frame. The THEORY block
|
||||
* should describe the low-frequency RRF signal, not the laboratory frame
|
||||
* precession.
|
||||
*
|
||||
* \see CalcChiSquareExpected(), CalcMaxLikelihood()
|
||||
*/
|
||||
Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
|
||||
{
|
||||
@@ -189,12 +276,32 @@ Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
|
||||
// CalcChiSquareExpected (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Calculate expected chi-square.
|
||||
* \brief Calculates expected χ² using theory variance instead of data variance.
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - chisq value
|
||||
* Computes the expected chi-square where the error estimate in the denominator
|
||||
* comes from the theory prediction rather than the data:
|
||||
* \f[
|
||||
* \chi^2_{\rm exp} = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{A_{\rm RRF}^{\rm theory}(t_i)}
|
||||
* \f]
|
||||
*
|
||||
* \param par parameter vector iterated by minuit2
|
||||
* This metric is useful for:
|
||||
* - Diagnostic purposes to assess fit quality
|
||||
* - Detecting systematic deviations from the model
|
||||
* - Comparing with standard χ² to identify error estimation issues
|
||||
*
|
||||
* <b>Algorithm:</b>
|
||||
* Same as CalcChiSquare() but uses theory value as variance estimate instead
|
||||
* of measured error bars. OpenMP parallelization is applied when available.
|
||||
*
|
||||
* \param par Parameter vector from MINUIT minimizer
|
||||
*
|
||||
* \return Expected χ² value. For a good fit, this should be approximately
|
||||
* equal to the number of degrees of freedom (N_bins - N_params).
|
||||
*
|
||||
* \warning Theory values must be positive for valid variance estimate.
|
||||
* Negative theory values can lead to incorrect χ² calculation.
|
||||
*
|
||||
* \see CalcChiSquare()
|
||||
*/
|
||||
Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>& par)
|
||||
{
|
||||
@@ -238,12 +345,35 @@ Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>&
|
||||
// CalcMaxLikelihood (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Calculate log maximum-likelihood. See http://pdg.lbl.gov/index.html
|
||||
* \brief Calculates maximum likelihood for RRF data (NOT YET IMPLEMENTED).
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - log maximum-likelihood value
|
||||
* Maximum likelihood estimation for RRF single histogram data is more complex
|
||||
* than for raw histograms due to the non-linear transformation from
|
||||
* Poisson-distributed counts to RRF asymmetry.
|
||||
*
|
||||
* \param par parameter vector iterated by minuit2
|
||||
* <b>Theoretical Background:</b>
|
||||
* For raw histogram data, the likelihood is:
|
||||
* \f[
|
||||
* \mathcal{L} = \prod_i \frac{\mu_i^{n_i} e^{-\mu_i}}{n_i!}
|
||||
* \f]
|
||||
* where \f$\mu_i\f$ is the expected count and \f$n_i\f$ is the observed count.
|
||||
*
|
||||
* For RRF-transformed data, the error propagation through the transformation
|
||||
* must be properly accounted for in the likelihood function.
|
||||
*
|
||||
* <b>Current Implementation:</b>
|
||||
* Returns 0.0 (not implemented). Use χ² minimization (CalcChiSquare) instead.
|
||||
*
|
||||
* \param par Parameter vector from MINUIT minimizer (unused)
|
||||
*
|
||||
* \return 0.0 (not implemented)
|
||||
*
|
||||
* \todo Implement proper maximum likelihood for RRF data by:
|
||||
* -# Deriving the likelihood function for transformed asymmetry
|
||||
* -# Accounting for error propagation through RRF transformation
|
||||
* -# Including correlations introduced by packing
|
||||
*
|
||||
* \see CalcChiSquare() for currently supported fit metric
|
||||
*/
|
||||
Double_t PRunSingleHistoRRF::CalcMaxLikelihood(const std::vector<Double_t>& par)
|
||||
{
|
||||
@@ -256,7 +386,39 @@ Double_t PRunSingleHistoRRF::CalcMaxLikelihood(const std::vector<Double_t>& par)
|
||||
// CalcTheory (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Calculate theory for a given set of fit-parameters.
|
||||
* \brief Evaluates theory function at all data points for viewing/plotting.
|
||||
*
|
||||
* Calculates the theoretical RRF asymmetry using the current MSR parameter
|
||||
* values and stores results in fData for display. This method is called
|
||||
* after fitting to generate the theory curve overlay.
|
||||
*
|
||||
* <b>Algorithm:</b>
|
||||
* -# Extract parameter values from MSR parameter list
|
||||
* -# Evaluate all user-defined functions from FUNCTIONS block
|
||||
* -# Loop over data points (size matches RRF-packed data)
|
||||
* -# Calculate time: t = dataTimeStart + i × dataTimeStep
|
||||
* -# Evaluate theory: P(t) = Func(t, par, funcValues)
|
||||
* -# Store results via fData.AppendTheoryValue()
|
||||
*
|
||||
* <b>Theory Function:</b>
|
||||
* The theory is evaluated directly in the RRF frame. The THEORY block should
|
||||
* specify the low-frequency RRF signal (after transformation), not the
|
||||
* laboratory-frame high-frequency precession.
|
||||
*
|
||||
* Example THEORY block for RRF analysis:
|
||||
* \code
|
||||
* THEORY
|
||||
* asymmetry 1 (amplitude)
|
||||
* simpleGss 2 (Gaussian relaxation)
|
||||
* TFieldCos map1 fun1 (frequency shift, phase)
|
||||
* \endcode
|
||||
*
|
||||
* where the frequency in TFieldCos is the difference frequency (ω - ω_RRF).
|
||||
*
|
||||
* \note Unlike CalcChiSquare(), this method does not return a value.
|
||||
* Results are stored internally in fData.fTheory vector.
|
||||
*
|
||||
* \see CalcChiSquare(), PrepareViewData()
|
||||
*/
|
||||
void PRunSingleHistoRRF::CalcTheory()
|
||||
{
|
||||
@@ -289,9 +451,23 @@ void PRunSingleHistoRRF::CalcTheory()
|
||||
// GetNoOfFitBins (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Calculate the number of fitted bins for the current fit range.
|
||||
* \brief Returns the number of bins included in the current fit range.
|
||||
*
|
||||
* <b>return:</b> number of fitted bins.
|
||||
* Triggers CalcNoOfFitBins() to ensure the bin count is current based on
|
||||
* the latest fit range settings, then returns the cached value.
|
||||
*
|
||||
* The number of fit bins is needed for:
|
||||
* - Calculating degrees of freedom: ν = N_bins - N_params
|
||||
* - Reduced χ²: χ²_red = χ² / ν
|
||||
* - Statistical diagnostics and fit quality assessment
|
||||
*
|
||||
* \return Number of RRF-packed bins within [fFitStartTime, fFitEndTime].
|
||||
* This accounts for RRF packing: fewer bins than raw data.
|
||||
*
|
||||
* \note The fit range may be modified during fitting by COMMANDS block
|
||||
* instructions. Always call this method to get the current count.
|
||||
*
|
||||
* \see CalcNoOfFitBins(), SetFitRangeBin()
|
||||
*/
|
||||
UInt_t PRunSingleHistoRRF::GetNoOfFitBins()
|
||||
{
|
||||
@@ -304,15 +480,51 @@ UInt_t PRunSingleHistoRRF::GetNoOfFitBins()
|
||||
// SetFitRangeBin (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Allows to change the fit range on the fly. Used in the COMMAND block.
|
||||
* The syntax of the string is: FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]].
|
||||
* If only one pair of fgb/lgb is given, it is used for all runs in the RUN block section.
|
||||
* If multiple fgb/lgb's are given, the number N has to be the number of RUN blocks in
|
||||
* the msr-file.
|
||||
* \brief Sets fit range using bin-offset syntax from COMMANDS block.
|
||||
*
|
||||
* <p>nXY are offsets which can be used to shift, limit the fit range.
|
||||
* Allows dynamic modification of the fit range during fitting, as specified
|
||||
* in the COMMANDS block. This is used for systematic studies where the fit
|
||||
* range needs to be varied.
|
||||
*
|
||||
* \param fitRange string containing the necessary information.
|
||||
* <b>Syntax:</b>
|
||||
* \code
|
||||
* FIT_RANGE fgb[+n0] lgb[-n1] # single range for all runs
|
||||
* FIT_RANGE fgb+n00 lgb-n01 fgb+n10 lgb-n11 ... # per-run ranges
|
||||
* \endcode
|
||||
*
|
||||
* where:
|
||||
* - \c fgb = first good bin (symbolic, replaced by actual value)
|
||||
* - \c lgb = last good bin (symbolic, replaced by actual value)
|
||||
* - \c +n0 = positive offset added to fgb
|
||||
* - \c -n1 = positive offset subtracted from lgb
|
||||
*
|
||||
* <b>Conversion to Time:</b>
|
||||
* \f[
|
||||
* t_{\rm start} = ({\rm fgb} + n_0 - t_0) \times \Delta t
|
||||
* \f]
|
||||
* \f[
|
||||
* t_{\rm end} = ({\rm lgb} - n_1 - t_0) \times \Delta t
|
||||
* \f]
|
||||
*
|
||||
* where Δt is the raw time resolution (before RRF packing).
|
||||
*
|
||||
* <b>Multiple Run Handling:</b>
|
||||
* - Single pair: Applied to all runs
|
||||
* - Multiple pairs: Must match number of RUN blocks; each run uses its own range
|
||||
* - Run selection: Uses (2 × runNo + 1) to index into token array
|
||||
*
|
||||
* <b>Example:</b>
|
||||
* \code
|
||||
* COMMANDS
|
||||
* FIT_RANGE fgb+100 lgb-500 # start 100 bins after fgb, end 500 bins before lgb
|
||||
* \endcode
|
||||
*
|
||||
* \param fitRange String containing FIT_RANGE specification from COMMANDS block
|
||||
*
|
||||
* \warning Invalid syntax produces error message but does not throw exception.
|
||||
* The previous fit range values are retained on error.
|
||||
*
|
||||
* \see GetProperFitRange(), CalcNoOfFitBins()
|
||||
*/
|
||||
void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange)
|
||||
{
|
||||
@@ -394,7 +606,37 @@ void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange)
|
||||
// CalcNoOfFitBins (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Calculate the number of fitted bins for the current fit range.
|
||||
* \brief Calculates start/end bin indices from fit time range.
|
||||
*
|
||||
* Converts the fit range from time (μs) to RRF-packed bin indices.
|
||||
* This method is called whenever the fit range may have changed.
|
||||
*
|
||||
* <b>Conversion Formulas:</b>
|
||||
* \f[
|
||||
* i_{\rm start} = \lceil \frac{t_{\rm start} - t_{\rm data,0}}{\Delta t_{\rm data}} \rceil
|
||||
* \f]
|
||||
* \f[
|
||||
* i_{\rm end} = \lfloor \frac{t_{\rm end} - t_{\rm data,0}}{\Delta t_{\rm data}} \rfloor + 1
|
||||
* \f]
|
||||
*
|
||||
* where:
|
||||
* - \f$t_{\rm data,0}\f$ = fData.GetDataTimeStart() (first RRF-packed bin center)
|
||||
* - \f$\Delta t_{\rm data}\f$ = fData.GetDataTimeStep() (RRF-packed bin width)
|
||||
*
|
||||
* <b>Bounds Checking:</b>
|
||||
* - fStartTimeBin clamped to [0, data size)
|
||||
* - fEndTimeBin clamped to [0, data size]
|
||||
* - fNoOfFitBins = 0 if fEndTimeBin <= fStartTimeBin
|
||||
*
|
||||
* <b>Side Effects:</b>
|
||||
* Updates member variables:
|
||||
* - fStartTimeBin: First bin index in fit range (inclusive)
|
||||
* - fEndTimeBin: Last bin index in fit range (exclusive)
|
||||
* - fNoOfFitBins: Number of bins = fEndTimeBin - fStartTimeBin
|
||||
*
|
||||
* \note Time step includes RRF packing: dataTimeStep = rawTimeRes × fRRFPacking
|
||||
*
|
||||
* \see GetNoOfFitBins(), SetFitRangeBin()
|
||||
*/
|
||||
void PRunSingleHistoRRF::CalcNoOfFitBins()
|
||||
{
|
||||
@@ -416,17 +658,50 @@ void PRunSingleHistoRRF::CalcNoOfFitBins()
|
||||
// PrepareData (protected)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Prepare data for fitting or viewing. What is already processed at this stage:
|
||||
* -# get proper raw run data
|
||||
* -# get all needed forward histograms
|
||||
* -# get time resolution
|
||||
* -# get t0's and perform necessary cross checks (e.g. if t0 of msr-file (if present) are consistent with t0 of the data files, etc.)
|
||||
* -# add runs (if addruns are present)
|
||||
* -# group histograms (if grouping is present)
|
||||
* \brief Main data preparation orchestrator for RRF single histogram analysis.
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - true if everthing went smooth
|
||||
* - false, otherwise.
|
||||
* Coordinates the loading and preprocessing of histogram data before RRF
|
||||
* transformation. This method handles all operations common to both fitting
|
||||
* and viewing modes.
|
||||
*
|
||||
* <b>Processing Steps:</b>
|
||||
* -# Validate object state (return false if already marked invalid)
|
||||
* -# Get GLOBAL block reference for settings
|
||||
* -# Load raw run data from data handler using run name from MSR file
|
||||
* -# Extract metadata from data file:
|
||||
* - Magnetic field (fMetaData.fField)
|
||||
* - Beam energy (fMetaData.fEnergy)
|
||||
* - Sample temperature(s) (fMetaData.fTemp)
|
||||
* -# Collect histogram numbers from RUN block forward specification
|
||||
* -# Validate all specified histograms exist in data file
|
||||
* -# Determine time resolution: ns → μs conversion
|
||||
* -# Determine t0 values via GetProperT0()
|
||||
* -# Initialize forward histogram from first group
|
||||
* -# Handle addrun (co-add multiple runs with t0 alignment)
|
||||
* -# Handle grouping (combine multiple detectors with t0 alignment)
|
||||
* -# Determine data range via GetProperDataRange()
|
||||
* -# Determine fit range via GetProperFitRange()
|
||||
* -# Dispatch to PrepareFitData() or PrepareViewData() based on tag
|
||||
*
|
||||
* <b>Run Addition (addrun):</b>
|
||||
* When multiple runs are specified in the RUN block, histograms are co-added
|
||||
* with t0 alignment:
|
||||
* \code
|
||||
* fForward[i] += addRunData[i + addT0 - mainT0]
|
||||
* \endcode
|
||||
*
|
||||
* <b>Detector Grouping:</b>
|
||||
* When multiple forward histograms are specified, they are summed with
|
||||
* t0 alignment to form a single combined histogram.
|
||||
*
|
||||
* \return true if data preparation succeeds, false on any error:
|
||||
* - Run data not found
|
||||
* - Histogram not present in data file
|
||||
* - Invalid t0 determination
|
||||
* - Data range validation failure
|
||||
* - Fit/view data preparation failure
|
||||
*
|
||||
* \see PrepareFitData(), PrepareViewData(), GetProperT0(), GetProperDataRange()
|
||||
*/
|
||||
Bool_t PRunSingleHistoRRF::PrepareData()
|
||||
{
|
||||
@@ -560,21 +835,73 @@ Bool_t PRunSingleHistoRRF::PrepareData()
|
||||
// PrepareFitData (protected)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Take the pre-processed data (i.e. grouping and addrun are preformed) and form the RRF histogram for fitting.
|
||||
* The following steps are preformed:
|
||||
* -# get fit start/stop time
|
||||
* -# check that 'first good data bin', 'last good data bin', and 't0' make any sense
|
||||
* -# check how the background shall be handled, i.e. fitted, subtracted from background estimate data range, or subtacted from a given fixed background.
|
||||
* -# estimate N0
|
||||
* -# RRF transformation
|
||||
* -# packing (i.e rebinning)
|
||||
* \brief Performs full RRF transformation on histogram data for fitting.
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - true, if everything went smooth
|
||||
* - false, otherwise
|
||||
* Takes the pre-processed histogram (grouping and addrun already applied) and
|
||||
* transforms it into RRF asymmetry suitable for fitting. This is the core
|
||||
* method implementing the RRF analysis technique.
|
||||
*
|
||||
* \param runData raw run data handler
|
||||
* \param histoNo forward histogram number
|
||||
* <b>Processing Steps:</b>
|
||||
*
|
||||
* <b>1. Frequency Analysis:</b>
|
||||
* - Calls GetMainFrequency() on raw data to find dominant precession frequency
|
||||
* - Used to determine optimal N₀ estimation window (integer oscillation cycles)
|
||||
* - Prints "optimal packing" suggestion for user reference
|
||||
*
|
||||
* <b>2. Background Handling:</b>
|
||||
* - If fixed background specified: use directly from RUN block
|
||||
* - If background range given: estimate via EstimateBkg()
|
||||
* - If neither: estimate range as [0.1×t0, 0.6×t0] with warning
|
||||
* - Subtract background: N(t) → N(t) - B
|
||||
*
|
||||
* <b>3. Lifetime Correction:</b>
|
||||
* \f[
|
||||
* M(t) = [N(t) - B] \cdot e^{+t/\tau_\mu}
|
||||
* \f]
|
||||
* - Removes exponential decay from signal
|
||||
* - Error: \f$\sigma_M = e^{+t/\tau_\mu} \sqrt{N(t) + \sigma_B^2}\f$
|
||||
* - Weights: \f$w = 1/\sigma_M^2\f$
|
||||
*
|
||||
* <b>4. N₀ Estimation:</b>
|
||||
* - Call EstimateN0() with frequency information
|
||||
* - Uses average of M(t) over integer number of oscillation cycles
|
||||
* - Returns N₀ and its error σ_N₀
|
||||
*
|
||||
* <b>5. Asymmetry Extraction:</b>
|
||||
* \f[
|
||||
* A(t) = \frac{M(t)}{N_0} - 1
|
||||
* \f]
|
||||
* Error propagation:
|
||||
* \f[
|
||||
* \sigma_A(t) = \frac{e^{+t/\tau_\mu}}{N_0} \sqrt{N(t) + \left(\frac{N(t)-B}{N_0}\right)^2 \sigma_{N_0}^2}
|
||||
* \f]
|
||||
*
|
||||
* <b>6. RRF Rotation:</b>
|
||||
* \f[
|
||||
* A_{\rm RRF}(t) = 2 \cdot A(t) \cdot \cos(\omega_{\rm RRF} t + \phi_{\rm RRF})
|
||||
* \f]
|
||||
* - ω_RRF from GLOBAL block (rrf_freq converted to angular frequency)
|
||||
* - φ_RRF from GLOBAL block (rrf_phase in radians)
|
||||
* - Factor 2 compensates for discarded counter-rotating component
|
||||
*
|
||||
* <b>7. RRF Packing:</b>
|
||||
* - Average over fRRFPacking consecutive bins
|
||||
* - Data: \f$A_{\rm packed} = \frac{1}{n}\sum_{i=1}^{n} A_{\rm RRF}(t_i)\f$
|
||||
* - Error: \f$\sigma_{\rm packed} = \frac{\sqrt{2}}{n}\sqrt{\sum_{i=1}^{n} \sigma_A^2(t_i)}\f$
|
||||
* - √2 factor accounts for doubling in RRF rotation
|
||||
*
|
||||
* <b>8. Time Grid Setup:</b>
|
||||
* - Data time start: accounts for packing offset (center of first packed bin)
|
||||
* - Data time step: rawTimeResolution × fRRFPacking
|
||||
*
|
||||
* \param runData Raw run data handler (for background estimation)
|
||||
* \param histoNo Forward histogram number (0-based, for background error messages)
|
||||
*
|
||||
* \return true on success, false on error:
|
||||
* - Background estimation failure
|
||||
* - Invalid bin ranges
|
||||
*
|
||||
* \see PrepareViewData(), GetMainFrequency(), EstimateN0(), EstimateBkg()
|
||||
*/
|
||||
Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t histoNo)
|
||||
{
|
||||
@@ -712,20 +1039,55 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his
|
||||
// PrepareViewData (protected)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing
|
||||
* with life time correction, i.e. the exponential decay is removed.
|
||||
* <p>The following steps are preformed:
|
||||
* -# check if view packing is whished.
|
||||
* -# check that 'first good data bin', 'last good data bin', and 't0' makes any sense
|
||||
* -# transform data sets (see below).
|
||||
* -# calculate theory
|
||||
* \brief Prepares RRF-transformed data for viewing/plotting.
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - true, if everything went smooth
|
||||
* - false, otherwise
|
||||
* Takes the pre-processed histogram and prepares both data and theory curves
|
||||
* for display. This method is used when the operation mode is kView rather
|
||||
* than kFit.
|
||||
*
|
||||
* \param runData raw run data handler
|
||||
* \param histoNo forward histogram number
|
||||
* <b>Processing Steps:</b>
|
||||
*
|
||||
* <b>1. Data Preparation:</b>
|
||||
* - Calls PrepareFitData() to perform full RRF transformation
|
||||
* - Data is ready for display after this step
|
||||
*
|
||||
* <b>2. View Packing Check:</b>
|
||||
* - Checks if additional view packing is specified in PLOT block
|
||||
* - If viewPacking > fRRFPacking: additional packing would be applied (TODO)
|
||||
* - If viewPacking < fRRFPacking: warning issued and view packing ignored
|
||||
* (cannot unpack already-packed RRF data)
|
||||
*
|
||||
* <b>3. Theory Grid Setup:</b>
|
||||
* Determines theory evaluation resolution:
|
||||
* - If fTheoAsData = true: theory calculated only at data points
|
||||
* - If fTheoAsData = false: theory calculated on 8× finer grid for smooth curves
|
||||
*
|
||||
* Time grid:
|
||||
* - Theory time start = data time start
|
||||
* - Theory time step = data time step / factor (where factor = 1 or 8)
|
||||
*
|
||||
* <b>4. Theory Evaluation:</b>
|
||||
* - Extract parameter values from MSR parameter list
|
||||
* - Evaluate all user-defined functions from FUNCTIONS block
|
||||
* - Loop over theory grid points
|
||||
* - Evaluate theory: P(t) = Func(t, par, funcValues)
|
||||
* - Apply sanity check: |P(t)| > 10 → set to 0 (dirty hack for edge cases)
|
||||
* - Store results via fData.AppendTheoryValue()
|
||||
*
|
||||
* <b>Theory Function:</b>
|
||||
* The theory is evaluated directly in the RRF frame. The THEORY block should
|
||||
* specify the low-frequency signal after RRF transformation, not the
|
||||
* laboratory-frame high-frequency precession.
|
||||
*
|
||||
* \param runData Raw run data handler (passed to PrepareFitData)
|
||||
* \param histoNo Forward histogram number (passed to PrepareFitData)
|
||||
*
|
||||
* \return true on success, false on error (typically from PrepareFitData)
|
||||
*
|
||||
* \note The 8× theory resolution provides smoother curves for visualization
|
||||
* and better Fourier transforms without affecting fit performance.
|
||||
*
|
||||
* \see PrepareFitData(), CalcTheory()
|
||||
*/
|
||||
Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t histoNo)
|
||||
{
|
||||
@@ -793,21 +1155,52 @@ Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t hi
|
||||
// GetProperT0 (private)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Get the proper t0 for the single histogram run.
|
||||
* -# the t0 vector size = number of detectors (grouping) for forward.
|
||||
* -# initialize t0's with -1
|
||||
* -# fill t0's from RUN block
|
||||
* -# if t0's are missing (i.e. t0 == -1), try to fill from the GLOBAL block.
|
||||
* -# if t0's are missing, try t0's from the data file
|
||||
* -# if t0's are missing, try to estimate them
|
||||
* \brief Determines and validates t0 values for all histograms in the run.
|
||||
*
|
||||
* \param runData pointer to the current RUN block entry from the msr-file
|
||||
* \param globalBlock pointer to the GLOBLA block entry from the msr-file
|
||||
* \param histoNo histogram number vector of forward; histoNo = msr-file forward + redGreen_offset - 1
|
||||
* Searches for time-zero (t0) values through a priority cascade and performs
|
||||
* validation. t0 marks the muon arrival time and is critical for proper
|
||||
* timing alignment.
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - true if everthing went smooth
|
||||
* - false, otherwise.
|
||||
* <b>t0 Search Priority:</b>
|
||||
* -# RUN block: t0 values specified in the MSR file RUN block
|
||||
* -# GLOBAL block: fallback t0 values from GLOBAL block
|
||||
* -# Data file: t0 values stored in the raw data file header
|
||||
* -# Automatic estimation: estimated from histogram shape (with warning)
|
||||
*
|
||||
* <b>Algorithm:</b>
|
||||
* -# Initialize fT0s vector with size = number of forward histograms (grouping)
|
||||
* -# Set all t0 values to -1.0 (unset marker)
|
||||
* -# Fill from RUN block if present
|
||||
* -# Fill remaining -1.0 values from GLOBAL block
|
||||
* -# Fill remaining -1.0 values from data file
|
||||
* -# Fill remaining -1.0 values from automatic estimation (with warning)
|
||||
* -# Validate all t0 values are within histogram bounds
|
||||
*
|
||||
* <b>Addrun Handling:</b>
|
||||
* When multiple runs are co-added (addrun), each additional run needs its own
|
||||
* t0 value for proper time alignment:
|
||||
* -# Initialize fAddT0s[runIdx] for each addrun
|
||||
* -# Apply same priority cascade for each addrun
|
||||
* -# Validate addrun t0 values
|
||||
*
|
||||
* <b>Validation:</b>
|
||||
* - t0 must be ≥ 0
|
||||
* - t0 must be < histogram length
|
||||
* - Error message and return false on validation failure
|
||||
*
|
||||
* \param runData Pointer to raw run data for histogram access and data file t0
|
||||
* \param globalBlock Pointer to GLOBAL block for fallback t0 values
|
||||
* \param histoNo Vector of histogram indices (0-based) for forward grouping
|
||||
*
|
||||
* \return true if valid t0 found for all histograms, false on:
|
||||
* - t0 out of histogram bounds
|
||||
* - addrun data not accessible
|
||||
* - addrun t0 validation failure
|
||||
*
|
||||
* \warning Automatic t0 estimation is unreliable for LEM and other specialized
|
||||
* setups. Always specify t0 explicitly for best results.
|
||||
*
|
||||
* \see GetProperDataRange(), PrepareData()
|
||||
*/
|
||||
Bool_t PRunSingleHistoRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo)
|
||||
{
|
||||
@@ -931,14 +1324,40 @@ Bool_t PRunSingleHistoRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *gl
|
||||
// GetProperDataRange (private)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Get the proper data range, i.e. first/last good bin (fgb/lgb).
|
||||
* -# get fgb/lgb from the RUN block
|
||||
* -# if fgb/lgb still undefined, try to get it from the GLOBAL block
|
||||
* -# if fgb/lgb still undefined, try to estimate them.
|
||||
* \brief Determines valid data range (first/last good bins) for analysis.
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - true if everthing went smooth
|
||||
* - false, otherwise.
|
||||
* Establishes the boundaries of usable histogram data through a priority
|
||||
* cascade. The data range defines which bins contain valid signal (after t0,
|
||||
* before noise/overflow).
|
||||
*
|
||||
* <b>Data Range Search Priority:</b>
|
||||
* -# RUN block: data range specified in MSR file RUN block
|
||||
* -# GLOBAL block: fallback data range from GLOBAL block
|
||||
* -# Automatic estimation: estimate from t0 (with warning)
|
||||
* - Start: t0 + 10 ns offset
|
||||
* - End: histogram length
|
||||
*
|
||||
* <b>Validation:</b>
|
||||
* -# Check start < end (swap if reversed)
|
||||
* -# Check start ≥ 0 and start < histogram size
|
||||
* -# Check end ≥ 0 (adjust if > histogram size with warning)
|
||||
*
|
||||
* <b>Result:</b>
|
||||
* Sets fGoodBins[0] (first good bin) and fGoodBins[1] (last good bin).
|
||||
* These are used for:
|
||||
* - RRF transformation range
|
||||
* - Fit range specification in bins (fgb+n0 lgb-n1)
|
||||
* - Data validity checking
|
||||
*
|
||||
* \return true if valid data range determined, false on:
|
||||
* - Start bin out of bounds
|
||||
* - End bin negative
|
||||
* - Other range validation failures
|
||||
*
|
||||
* \note The data range is typically wider than the fit range. Data range
|
||||
* defines where valid data exists; fit range defines what is fitted.
|
||||
*
|
||||
* \see GetProperFitRange(), GetProperT0()
|
||||
*/
|
||||
Bool_t PRunSingleHistoRRF::GetProperDataRange()
|
||||
{
|
||||
@@ -1019,16 +1438,45 @@ Bool_t PRunSingleHistoRRF::GetProperDataRange()
|
||||
// GetProperFitRange (private)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Get the proper fit range. There are two possible fit range commands:
|
||||
* fit <start> <end> given in (usec), or
|
||||
* fit fgb+offset_0 lgb-offset_1 given in (bins), therefore it works the following way:
|
||||
* -# get fit range assuming given in time from RUN block
|
||||
* -# if fit range in RUN block is given in bins, replace start/end
|
||||
* -# if fit range is NOT given yet, try fit range assuming given in time from GLOBAL block
|
||||
* -# if fit range in GLOBAL block is given in bins, replace start/end
|
||||
* -# if still no fit range is given, use fgb/lgb.
|
||||
* \brief Determines fit time range from MSR file settings.
|
||||
*
|
||||
* \param globalBlock pointer to the GLOBAL block information form the msr-file.
|
||||
* Extracts the fit range (time window for parameter extraction) through a
|
||||
* priority cascade. Supports both time-based and bin-based specifications.
|
||||
*
|
||||
* <b>Fit Range Formats:</b>
|
||||
* - Time-based: \c fit \c 0.1 \c 8.0 (in microseconds relative to t0)
|
||||
* - Bin-based: \c fit \c fgb+100 \c lgb-500 (bins relative to good bins)
|
||||
*
|
||||
* <b>Search Priority:</b>
|
||||
* -# RUN block time-based fit range
|
||||
* -# RUN block bin-based fit range (converted to time)
|
||||
* -# GLOBAL block time-based fit range
|
||||
* -# GLOBAL block bin-based fit range (converted to time)
|
||||
* -# Default to data range (fgb to lgb) with warning
|
||||
*
|
||||
* <b>Bin-to-Time Conversion:</b>
|
||||
* \f[
|
||||
* t_{\rm start} = ({\rm fgb} + n_0 - t_0) \times \Delta t
|
||||
* \f]
|
||||
* \f[
|
||||
* t_{\rm end} = ({\rm lgb} - n_1 - t_0) \times \Delta t
|
||||
* \f]
|
||||
*
|
||||
* where Δt = fTimeResolution (raw, before RRF packing).
|
||||
*
|
||||
* <b>Result:</b>
|
||||
* Sets fFitStartTime and fFitEndTime in microseconds relative to t0.
|
||||
* These define the range used in χ² calculation.
|
||||
*
|
||||
* \param globalBlock Pointer to GLOBAL block for fallback fit range settings
|
||||
*
|
||||
* \note The converted time values are written back to the MSR data structures
|
||||
* for log file output consistency.
|
||||
*
|
||||
* \warning If no fit range is specified, the full data range is used with
|
||||
* a warning message. This may not be appropriate for all analyses.
|
||||
*
|
||||
* \see GetProperDataRange(), SetFitRangeBin(), CalcNoOfFitBins()
|
||||
*/
|
||||
void PRunSingleHistoRRF::GetProperFitRange(PMsrGlobalBlock *globalBlock)
|
||||
{
|
||||
@@ -1067,9 +1515,41 @@ void PRunSingleHistoRRF::GetProperFitRange(PMsrGlobalBlock *globalBlock)
|
||||
// GetMainFrequency (private)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>gets the maximum frequency of the given data.
|
||||
* \brief Finds the dominant precession frequency in raw histogram data.
|
||||
*
|
||||
* \param raw data set.
|
||||
* Performs Fourier analysis on the raw histogram to identify the main muon
|
||||
* spin precession frequency. This is essential for:
|
||||
* - Determining optimal N₀ estimation window (integer oscillation cycles)
|
||||
* - Calculating suggested "optimal packing" for user information
|
||||
* - Validating that the RRF frequency is appropriate
|
||||
*
|
||||
* <b>Algorithm:</b>
|
||||
* -# Create ROOT TH1F histogram from data in good bin range [fGoodBins[0], fGoodBins[1]]
|
||||
* -# Set histogram binning to match time structure
|
||||
* -# Apply Fourier transform using PFourier class
|
||||
* -# Use strong apodization (windowing) to reduce spectral leakage
|
||||
* -# Search power spectrum for maximum above 10 MHz (ignores DC component)
|
||||
* -# Return frequency at maximum power
|
||||
*
|
||||
* <b>Frequency Search Constraints:</b>
|
||||
* - Ignores frequencies below 10 MHz (DC and low-frequency noise)
|
||||
* - Searches for local maximum (rising then falling power)
|
||||
* - Returns 0 if no clear maximum found
|
||||
*
|
||||
* <b>Output Information:</b>
|
||||
* The method prints diagnostic information:
|
||||
* - Detected maximum frequency (MHz)
|
||||
* - Suggested optimal packing for ~5-8 points per cycle
|
||||
*
|
||||
* \param data Raw histogram data vector (counts per bin)
|
||||
*
|
||||
* \return Maximum frequency in MHz, or 0 if no maximum found above 10 MHz.
|
||||
*
|
||||
* \note The frequency is used to determine the N₀ estimation window:
|
||||
* window = ceil(fN0EstimateEndTime × freqMax / 2π) × (2π / freqMax)
|
||||
* This ensures an integer number of complete oscillation cycles.
|
||||
*
|
||||
* \see EstimateN0(), PrepareFitData()
|
||||
*/
|
||||
Double_t PRunSingleHistoRRF::GetMainFrequency(PDoubleVector &data)
|
||||
{
|
||||
@@ -1117,9 +1597,50 @@ Double_t PRunSingleHistoRRF::GetMainFrequency(PDoubleVector &data)
|
||||
// EstimateN0 (private)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Estimate the N0 for the given run.
|
||||
* \brief Estimates initial normalization N₀ from lifetime-corrected histogram data.
|
||||
*
|
||||
* \param errN0
|
||||
* Calculates N₀ by averaging the lifetime-corrected signal M(t) over an initial
|
||||
* time window. The window is chosen to span an integer number of complete
|
||||
* precession cycles to minimize bias from oscillations.
|
||||
*
|
||||
* <b>N₀ Estimation Formula:</b>
|
||||
* \f[
|
||||
* N_0 = \frac{1}{n} \sum_{i=0}^{n-1} M(t_i)
|
||||
* \f]
|
||||
*
|
||||
* where M(t) = [N(t) - B] × exp(+t/τ_μ) is the lifetime-corrected histogram
|
||||
* stored in fM vector.
|
||||
*
|
||||
* <b>Window Determination:</b>
|
||||
* The end bin is calculated to include an integer number of oscillation cycles:
|
||||
* \f[
|
||||
* n_{\rm end} = {\rm round}\left( \left\lceil \frac{T \cdot f_{\rm max}}{2\pi} \right\rceil \cdot \frac{2\pi}{f_{\rm max} \cdot \Delta t} \right)
|
||||
* \f]
|
||||
*
|
||||
* where:
|
||||
* - T = fN0EstimateEndTime (default 1.0 μs)
|
||||
* - f_max = main precession frequency from GetMainFrequency()
|
||||
* - Δt = fTimeResolution
|
||||
*
|
||||
* <b>Error Estimation:</b>
|
||||
* \f[
|
||||
* \sigma_{N_0} = \frac{\sqrt{\sum_{i=0}^{n-1} w_i^2 \sigma_{M_i}^2}}{\sum_{i=0}^{n-1} w_i}
|
||||
* \f]
|
||||
*
|
||||
* where w_i = 1/σ_M_i² are the weights stored in fW vector.
|
||||
*
|
||||
* <b>Output Information:</b>
|
||||
* Prints diagnostic message: "N₀ = value (error)"
|
||||
*
|
||||
* \param errN0 [out] Reference to store the estimated error on N₀
|
||||
* \param freqMax Main precession frequency (MHz) from GetMainFrequency()
|
||||
*
|
||||
* \return Estimated N₀ value (counts × exp(+t/τ))
|
||||
*
|
||||
* \note The simple average (rather than weighted average) is used for N₀ itself,
|
||||
* while the error uses the full weight information.
|
||||
*
|
||||
* \see GetMainFrequency(), PrepareFitData()
|
||||
*/
|
||||
Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax)
|
||||
{
|
||||
@@ -1152,13 +1673,47 @@ Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax)
|
||||
// EstimatBkg (private)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Estimate the background for a given interval.
|
||||
* \brief Estimates background level and error from pre-t0 histogram bins.
|
||||
*
|
||||
* <b>return:</b>
|
||||
* - true, if everything went smooth
|
||||
* - false, otherwise
|
||||
* Calculates the background (random coincidences, accidentals) from a
|
||||
* specified range of histogram bins, typically before t0. The background
|
||||
* range is adjusted to span an integer number of accelerator RF cycles
|
||||
* for proper averaging.
|
||||
*
|
||||
* \param histoNo forward histogram number of the run
|
||||
* <b>Accelerator RF Periods:</b>
|
||||
* - PSI: ACCEL_PERIOD_PSI ns
|
||||
* - RAL: ACCEL_PERIOD_RAL ns
|
||||
* - TRIUMF: ACCEL_PERIOD_TRIUMF ns
|
||||
* - Unknown: no RF adjustment
|
||||
*
|
||||
* <b>Algorithm:</b>
|
||||
* -# Get background range [start, end] from RUN block
|
||||
* -# Swap start/end if in wrong order (with message)
|
||||
* -# Determine accelerator RF period from institute name
|
||||
* -# Adjust end to span integer number of RF cycles
|
||||
* -# Validate range is within histogram bounds
|
||||
* -# Calculate mean background:
|
||||
* \f$ B = \frac{1}{n}\sum_{i={\rm start}}^{{\rm end}} N(i) \f$
|
||||
* -# Calculate background error (standard deviation):
|
||||
* \f$ \sigma_B = \sqrt{\frac{1}{n-1}\sum_{i={\rm start}}^{{\rm end}} (N(i) - B)^2} \f$
|
||||
* -# Store results in fBackground and fBkgErr
|
||||
* -# Update RUN block with estimated background value
|
||||
*
|
||||
* <b>Output Information:</b>
|
||||
* Prints diagnostic messages:
|
||||
* - Adjusted background range (if RF adjustment applied)
|
||||
* - Background value and error: "fBackground = value (error)"
|
||||
*
|
||||
* \param histoNo Forward histogram number (for error message context)
|
||||
*
|
||||
* \return true on success, false on error:
|
||||
* - Background start out of histogram bounds
|
||||
* - Background end out of histogram bounds
|
||||
*
|
||||
* \note The RF cycle adjustment ensures proper averaging over the
|
||||
* accelerator bunch structure, reducing systematic bias.
|
||||
*
|
||||
* \see PrepareFitData()
|
||||
*/
|
||||
Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo)
|
||||
{
|
||||
|
||||
@@ -33,59 +33,502 @@
|
||||
#include "PRunBase.h"
|
||||
|
||||
/**
|
||||
* <p>Class handling single histogram fit type.
|
||||
* \brief Class for fitting single histogram data in a Rotating Reference Frame (RRF).
|
||||
*
|
||||
* PRunSingleHistoRRF implements single histogram analysis with Rotating Reference Frame
|
||||
* transformation. This technique is essential for high transverse field (TF) μSR
|
||||
* measurements where the precession frequency is too high to resolve directly.
|
||||
*
|
||||
* \section rrf_physics Physics Background
|
||||
*
|
||||
* <b>The Problem with High-Field μSR:</b>
|
||||
* - At fields B > ~0.5 T, muon precession frequency ω = γ_μ B exceeds ~70 MHz
|
||||
* - Standard time-domain analysis struggles: requires sub-ns binning → low statistics
|
||||
* - Frequency-domain (Fourier) analysis has limited resolution
|
||||
*
|
||||
* <b>The RRF Solution:</b>
|
||||
* - Transform to a reference frame rotating at frequency ω_RRF ≈ ω_signal
|
||||
* - In rotating frame: effective frequency Δω = ω_signal - ω_RRF << ω_signal
|
||||
* - Low-frequency signal can be fit with standard time-domain techniques
|
||||
* - Preserves full information: relaxation rates, frequency distributions
|
||||
*
|
||||
* \section rrf_transform Mathematical Transformation
|
||||
*
|
||||
* <b>Data processing steps:</b>
|
||||
* -# <b>Background subtraction:</b> N(t) → N(t) - B
|
||||
* -# <b>Lifetime correction:</b> N(t) - B → [N(t) - B] × exp(+t/τ_μ) = M(t)
|
||||
* -# <b>N₀ estimation:</b> Fit M(t) over initial time window to extract N₀
|
||||
* -# <b>Asymmetry extraction:</b> A(t) = M(t)/N₀ - 1
|
||||
* -# <b>RRF rotation:</b> A_RRF(t) = 2 × A(t) × cos(ω_RRF t + φ_RRF)
|
||||
* -# <b>RRF packing:</b> Average over multiple bins to filter high-frequency components
|
||||
*
|
||||
* The factor of 2 in step 5 compensates for the loss of the counter-rotating component
|
||||
* when the high-frequency part is suppressed by packing.
|
||||
*
|
||||
* <b>Theory function in RRF:</b>
|
||||
* The polarization function P(t) transforms as:
|
||||
* \f[
|
||||
* P_{\rm RRF}(t) = 2 \cdot P(t) \cdot \cos(\omega_{\rm RRF} t + \phi_{\rm RRF})
|
||||
* \f]
|
||||
*
|
||||
* For a simple oscillation P(t) = A·cos(ωt + φ)·exp(-λt), this becomes:
|
||||
* \f[
|
||||
* P_{\rm RRF}(t) = A \cdot [\cos((\omega - \omega_{\rm RRF})t + \phi - \phi_{\rm RRF}) + \text{high freq.}] \cdot e^{-\lambda t}
|
||||
* \f]
|
||||
*
|
||||
* After RRF packing, only the low-frequency component Δω = ω - ω_RRF remains.
|
||||
*
|
||||
* \section rrf_requirements GLOBAL Block Requirements
|
||||
*
|
||||
* RRF analysis REQUIRES the following entries in the GLOBAL block:
|
||||
* - <b>rrf_freq:</b> RRF rotation frequency (value + unit: MHz, Mc, T, G, kG)
|
||||
* - <b>rrf_phase:</b> RRF initial phase in degrees
|
||||
* - <b>rrf_packing:</b> Number of bins to average (filters high frequencies)
|
||||
*
|
||||
* \section rrf_msr MSR File Example
|
||||
*
|
||||
* \code
|
||||
* GLOBAL
|
||||
* rrf_freq 13.554 T # RRF frequency equivalent to 13.554 Tesla
|
||||
* rrf_phase 45.0 # RRF phase in degrees
|
||||
* rrf_packing 25 # Average 25 bins to suppress high-freq component
|
||||
* data 200 60000 # Data range (fgb lgb)
|
||||
* fit 0.1 8.0 # Fit range in μs
|
||||
*
|
||||
* RUN data/run2425 PSI MUE4 PSI MUSR-ROOT
|
||||
* fittype 8 (SingleHistoRRF)
|
||||
* map 1
|
||||
* forward 1
|
||||
* background 50 150
|
||||
* t0 210.5
|
||||
* \endcode
|
||||
*
|
||||
* \section rrf_workflow Data Processing Workflow
|
||||
*
|
||||
* 1. <b>PrepareData():</b> Load raw histogram, determine t0, group/add histograms
|
||||
* 2. <b>PrepareFitData():</b>
|
||||
* - Subtract background (estimated or fixed)
|
||||
* - Apply lifetime correction: exp(+t/τ_μ)
|
||||
* - Estimate N₀ using Fourier analysis to find main frequency
|
||||
* - Extract asymmetry: A(t) = M(t)/N₀ - 1
|
||||
* - Apply RRF rotation: A_RRF(t) = 2·A(t)·cos(ω_RRF t + φ_RRF)
|
||||
* - Pack data (average over rrf_packing bins)
|
||||
* - Calculate errors with proper propagation
|
||||
*
|
||||
* \section rrf_errors Error Propagation
|
||||
*
|
||||
* RRF asymmetry error (unpacked):
|
||||
* \f[
|
||||
* \sigma_{A}(t) = \frac{e^{t/\tau_\mu}}{N_0} \sqrt{N(t) + \left(\frac{N(t)-B}{N_0}\right)^2 \sigma_{N_0}^2}
|
||||
* \f]
|
||||
*
|
||||
* After RRF packing (n bins):
|
||||
* \f[
|
||||
* \sigma_{A_{\rm RRF}}^{\rm packed} = \frac{\sqrt{2}}{n} \sqrt{\sum_{i=1}^{n} \sigma_{A}^2(t_i)}
|
||||
* \f]
|
||||
*
|
||||
* The √2 factor accounts for the doubling in the RRF transformation.
|
||||
*
|
||||
* \section rrf_applications Applications
|
||||
*
|
||||
* <b>High-TF μSR experiments:</b>
|
||||
* - Knight shift measurements (NMR-like)
|
||||
* - Vortex state studies in type-II superconductors
|
||||
* - Diamagnetic/paramagnetic shift measurements
|
||||
* - Internal field distributions in magnets
|
||||
*
|
||||
* \see PRunSingleHisto for standard (non-RRF) single histogram analysis
|
||||
* \see PRunAsymmetryRRF for RRF asymmetry (forward-backward) analysis
|
||||
* \see PRunBase for base class interface
|
||||
*/
|
||||
class PRunSingleHistoRRF : public PRunBase
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* \brief Default constructor creating an empty, invalid RRF single histogram run object.
|
||||
*
|
||||
* Initializes member variables to default/invalid states:
|
||||
* - fNoOfFitBins = 0
|
||||
* - fBackground = 0
|
||||
* - fBkgErr = 1.0 (default error estimate)
|
||||
* - fRRFPacking = -1 (invalid, must be set from GLOBAL block)
|
||||
* - fTheoAsData = false
|
||||
* - fGoodBins[0,1] = -1 (unset)
|
||||
* - fN0EstimateEndTime = 1.0 μs
|
||||
*
|
||||
* This constructor exists for container compatibility. The resulting object
|
||||
* cannot be used until properly initialized.
|
||||
*/
|
||||
PRunSingleHistoRRF();
|
||||
|
||||
/**
|
||||
* \brief Main constructor initializing RRF single histogram run from MSR file and data.
|
||||
*
|
||||
* Performs comprehensive validation and initialization:
|
||||
*
|
||||
* 1. <b>GLOBAL Block Validation (MANDATORY):</b>
|
||||
* - Checks GLOBAL block is present (required for RRF)
|
||||
* - Validates RRF frequency is specified (rrf_freq)
|
||||
* - Validates RRF packing is specified (rrf_packing)
|
||||
* - Sets fValid=false and returns on any validation failure
|
||||
*
|
||||
* 2. <b>Member Initialization:</b>
|
||||
* - Extracts fRRFPacking from GLOBAL block
|
||||
* - Sets fN0EstimateEndTime = 1.0 μs
|
||||
* - Initializes fGoodBins to -1 (determined in PrepareData)
|
||||
*
|
||||
* 3. <b>Data Preparation:</b>
|
||||
* - Calls PrepareData() to load histogram and perform RRF transformation
|
||||
* - Sets fValid=false if data preparation fails
|
||||
*
|
||||
* \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 resolution: true = at data points, false = 8× finer grid
|
||||
*
|
||||
* \warning Always check IsValid() after construction. GLOBAL block with RRF
|
||||
* parameters is MANDATORY for this fit type.
|
||||
*
|
||||
* \see PrepareData() for RRF transformation details
|
||||
*/
|
||||
PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData);
|
||||
|
||||
/**
|
||||
* \brief Virtual destructor releasing allocated resources.
|
||||
*
|
||||
* Clears the forward histogram vector. Base class handles theory objects.
|
||||
*/
|
||||
virtual ~PRunSingleHistoRRF();
|
||||
|
||||
/**
|
||||
* \brief Calculates χ² between RRF-transformed data and theory.
|
||||
*
|
||||
* Computes chi-squared for RRF single histogram fitting:
|
||||
* \f[
|
||||
* \chi^2 = \sum_{i=t_{\rm start}}^{t_{\rm end}} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{\sigma_i^2}
|
||||
* \f]
|
||||
*
|
||||
* <b>Algorithm:</b>
|
||||
* -# Evaluate user-defined FUNCTIONS with current parameters
|
||||
* -# Pre-evaluate theory at t=1.0 to initialize LF/user functions (thread-safety)
|
||||
* -# Loop over fit range with OpenMP parallelization
|
||||
* -# Sum squared differences weighted by inverse variance
|
||||
*
|
||||
* <b>OpenMP Parallelization:</b>
|
||||
* - Dynamic scheduling with chunk = max(10, N_bins / N_processors)
|
||||
* - Reduction on chisq accumulator
|
||||
*
|
||||
* \param par Parameter vector from MINUIT minimizer
|
||||
* \return χ² value (sum over all bins in fit range)
|
||||
*
|
||||
* \note Unlike standard single histogram, no N₀/exp(-t/τ) factors needed since
|
||||
* RRF transformation already produces dimensionless asymmetry.
|
||||
*/
|
||||
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
|
||||
|
||||
/**
|
||||
* \brief Calculates expected χ² using theory variance instead of data variance.
|
||||
*
|
||||
* Computes expected chi-squared where the error estimate comes from the
|
||||
* theory prediction rather than the data:
|
||||
* \f[
|
||||
* \chi^2_{\rm exp} = \sum_{i} \frac{[A_{\rm RRF}^{\rm data}(t_i) - A_{\rm RRF}^{\rm theory}(t_i)]^2}{A_{\rm RRF}^{\rm theory}(t_i)}
|
||||
* \f]
|
||||
*
|
||||
* This is useful for diagnostic purposes to check fit quality.
|
||||
*
|
||||
* \param par Parameter vector from MINUIT
|
||||
* \return Expected χ² value
|
||||
*/
|
||||
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
|
||||
|
||||
/**
|
||||
* \brief Calculates maximum likelihood (not yet implemented for RRF).
|
||||
*
|
||||
* Maximum likelihood estimation for RRF data is more complex due to the
|
||||
* transformation from Poisson-distributed counts to RRF asymmetry.
|
||||
* Currently returns 0.0.
|
||||
*
|
||||
* \param par Parameter vector from MINUIT
|
||||
* \return 0.0 (not implemented)
|
||||
*
|
||||
* \todo Implement proper ML for RRF data accounting for error propagation
|
||||
* through the RRF transformation.
|
||||
*/
|
||||
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
|
||||
|
||||
/**
|
||||
* \brief Evaluates theory function at all data points for viewing/plotting.
|
||||
*
|
||||
* Calculates RRF theory values using current fit parameters stored in
|
||||
* the MSR parameter list. The theory function P(t) is evaluated at each
|
||||
* data time point and results are stored in fData for display.
|
||||
*
|
||||
* \note Theory is evaluated directly in the RRF frame, not transformed.
|
||||
* The THEORY block should specify the low-frequency RRF signal.
|
||||
*/
|
||||
virtual void CalcTheory();
|
||||
|
||||
/**
|
||||
* \brief Returns the number of bins included in the current fit range.
|
||||
*
|
||||
* Triggers CalcNoOfFitBins() to ensure bin count is current, then returns
|
||||
* the cached value.
|
||||
*
|
||||
* \return Number of RRF-packed bins within [fFitStartTime, fFitEndTime]
|
||||
*/
|
||||
virtual UInt_t GetNoOfFitBins();
|
||||
|
||||
/**
|
||||
* \brief Sets fit range using bin-offset syntax from COMMANDS block.
|
||||
*
|
||||
* Parses fit range specification in the format:
|
||||
* - "FIT_RANGE fgb+n0 lgb-n1" (single range for all runs)
|
||||
* - "FIT_RANGE fgb+n00 lgb-n01 fgb+n10 lgb-n11 ..." (per-run ranges)
|
||||
*
|
||||
* where fgb = first good bin, lgb = last good bin, and nXY are offsets.
|
||||
*
|
||||
* Converts bin specification to time using:
|
||||
* - fFitStartTime = (fgb + offset - t0) × time_resolution
|
||||
* - fFitEndTime = (lgb - offset - t0) × time_resolution
|
||||
*
|
||||
* \param fitRange String containing FIT_RANGE specification
|
||||
*
|
||||
* \note This is called when COMMANDS block modifies fit range during fitting.
|
||||
*/
|
||||
virtual void SetFitRangeBin(const TString fitRange);
|
||||
|
||||
/**
|
||||
* \brief Returns the first bin index in the fit range.
|
||||
* \return Start bin index (0-based, after RRF packing)
|
||||
*/
|
||||
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 Calculates start/end bin indices from fit time range.
|
||||
*
|
||||
* Converts fit range times (μs) to RRF-packed bin indices:
|
||||
* - fStartTimeBin = ceil((fFitStartTime - dataTimeStart) / dataTimeStep)
|
||||
* - fEndTimeBin = floor((fFitEndTime - dataTimeStart) / dataTimeStep) + 1
|
||||
*
|
||||
* Also updates fNoOfFitBins = fEndTimeBin - fStartTimeBin.
|
||||
*
|
||||
* \note Time step accounts for RRF packing: dataTimeStep = rawTimeStep × fRRFPacking
|
||||
*/
|
||||
virtual void CalcNoOfFitBins();
|
||||
|
||||
protected:
|
||||
/**
|
||||
* \brief Main data preparation orchestrator for RRF single histogram analysis.
|
||||
*
|
||||
* Coordinates loading and preprocessing of histogram data:
|
||||
* -# Validates GLOBAL block presence and RRF parameters
|
||||
* -# Retrieves raw run data from data handler
|
||||
* -# Extracts metadata (field, energy, temperature)
|
||||
* -# Collects histogram numbers from RUN block forward specification
|
||||
* -# Validates histograms exist in data file
|
||||
* -# Determines time resolution (ns → μs conversion)
|
||||
* -# Determines t0 values via GetProperT0()
|
||||
* -# Handles addrun (co-adding multiple runs)
|
||||
* -# Handles grouping (combining multiple detectors)
|
||||
* -# Determines data range via GetProperDataRange()
|
||||
* -# Determines fit range via GetProperFitRange()
|
||||
* -# Calls PrepareFitData() or PrepareViewData() based on tag
|
||||
*
|
||||
* \return true if data preparation succeeds, false on any error
|
||||
*
|
||||
* \see PrepareFitData(), PrepareViewData()
|
||||
*/
|
||||
virtual Bool_t PrepareData();
|
||||
|
||||
/**
|
||||
* \brief Performs full RRF transformation for fitting.
|
||||
*
|
||||
* Transforms raw histogram to RRF asymmetry through these steps:
|
||||
*
|
||||
* 1. <b>Frequency Analysis:</b>
|
||||
* - Calls GetMainFrequency() to find dominant precession frequency
|
||||
* - Used to determine optimal N₀ estimation window
|
||||
*
|
||||
* 2. <b>Background Handling:</b>
|
||||
* - If background range given: estimate from data via EstimateBkg()
|
||||
* - If fixed background given: use directly
|
||||
* - Subtract background: N(t) → N(t) - B
|
||||
*
|
||||
* 3. <b>Lifetime Correction:</b>
|
||||
* - Apply exp(+t/τ_μ) to remove exponential decay
|
||||
* - Store as M(t) = [N(t) - B] × exp(+t/τ_μ)
|
||||
* - Calculate M_err = exp(+t/τ_μ) × √(N(t) + σ_B²)
|
||||
*
|
||||
* 4. <b>N₀ Estimation:</b>
|
||||
* - Call EstimateN0() to determine normalization
|
||||
* - Uses weighted average over full oscillation cycles
|
||||
*
|
||||
* 5. <b>Asymmetry Extraction:</b>
|
||||
* - A(t) = M(t)/N₀ - 1
|
||||
* - A_err(t) = exp(+t/τ)/N₀ × √(N(t) + [(N(t)-B)/N₀]² × σ_N₀²)
|
||||
*
|
||||
* 6. <b>RRF Rotation:</b>
|
||||
* - A_RRF(t) = 2 × A(t) × cos(ω_RRF × t + φ_RRF)
|
||||
* - Factor 2 compensates for discarded counter-rotating component
|
||||
*
|
||||
* 7. <b>RRF Packing:</b>
|
||||
* - Average over fRRFPacking consecutive bins
|
||||
* - Error: σ_packed = √(2 × Σσ²) / n
|
||||
*
|
||||
* 8. <b>Time Grid Setup:</b>
|
||||
* - Set data time start accounting for packing offset
|
||||
* - Set data time step = raw_resolution × fRRFPacking
|
||||
*
|
||||
* \param runData Raw run data handler
|
||||
* \param histoNo Forward histogram number (0-based)
|
||||
* \return true on success, false on error
|
||||
*/
|
||||
virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo);
|
||||
|
||||
/**
|
||||
* \brief Prepares RRF data for viewing/plotting.
|
||||
*
|
||||
* Similar to PrepareFitData() but additionally:
|
||||
* - Handles view packing (if specified and > RRF packing)
|
||||
* - Sets up theory calculation grid (8× finer than data if fTheoAsData=false)
|
||||
* - Evaluates theory function for plotting overlay
|
||||
*
|
||||
* \param runData Raw run data handler
|
||||
* \param histoNo Forward histogram number
|
||||
* \return true on success, false on error
|
||||
*
|
||||
* \note View packing < RRF packing is ignored with warning since RRF
|
||||
* packing is already applied during data transformation.
|
||||
*/
|
||||
virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo);
|
||||
|
||||
private:
|
||||
Double_t fN0EstimateEndTime; ///< end time in (us) over which N0 is estimated.
|
||||
Double_t fN0EstimateEndTime; ///< End time (μs) for N₀ estimation window. Rounded to integer number of oscillation cycles based on main frequency.
|
||||
|
||||
UInt_t fNoOfFitBins; ///< number of bins to be fitted
|
||||
Double_t fBackground; ///< needed if background range is given (units: 1/bin)
|
||||
Double_t fBkgErr; ///< estimate error on the estimated background
|
||||
Int_t fRRFPacking; ///< RRF packing for this particular run. Given in the 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 RRF-packed bins within fit range [fStartTimeBin, fEndTimeBin)
|
||||
Double_t fBackground; ///< Estimated or fixed background level in counts/bin (before packing)
|
||||
Double_t fBkgErr; ///< Statistical error on background estimate (std dev of background region)
|
||||
Int_t fRRFPacking; ///< RRF packing factor from GLOBAL block (number of raw bins averaged together)
|
||||
Bool_t fTheoAsData; ///< Theory resolution mode: true = at data points only, false = 8× finer grid for smooth Fourier transforms
|
||||
|
||||
Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb
|
||||
Int_t fGoodBins[2]; ///< Good bin range: [0] = first good bin (fgb), [1] = last good bin (lgb). Used for COMMANDS block fit range specification.
|
||||
|
||||
Int_t fStartTimeBin; ///< bin at which the fit starts
|
||||
Int_t fEndTimeBin; ///< bin at which the fit ends
|
||||
Int_t fStartTimeBin; ///< First bin index in fit range (inclusive, 0-based in RRF-packed data)
|
||||
Int_t fEndTimeBin; ///< Last bin index in fit range (exclusive: loop as i < fEndTimeBin)
|
||||
|
||||
PDoubleVector fForward; ///< forward histo data
|
||||
PDoubleVector fM; ///< vector holding M(t) = [N(t)-N_bkg] exp(+t/tau). Needed to estimate N0.
|
||||
PDoubleVector fMerr; ///< vector holding the error of M(t): M_err = exp(+t/tau) sqrt(N(t)).
|
||||
PDoubleVector fW; ///< vector holding the weight needed to estimate N0, and errN0.
|
||||
PDoubleVector fAerr; ///< vector holding the errors of estimated A(t)
|
||||
PDoubleVector fForward; ///< Forward detector histogram data (progressively transformed during preparation)
|
||||
PDoubleVector fM; ///< Lifetime-corrected histogram: M(t) = [N(t) - B] × exp(+t/τ_μ). Used for N₀ estimation.
|
||||
PDoubleVector fMerr; ///< Error on M(t): σ_M = exp(+t/τ_μ) × √(N(t) + σ_B²). Includes background error.
|
||||
PDoubleVector fW; ///< Weights for N₀ estimation: W(t) = 1/σ_M². Used in weighted average.
|
||||
PDoubleVector fAerr; ///< Asymmetry errors before RRF packing. Used for packed error calculation.
|
||||
|
||||
/**
|
||||
* \brief Determines and validates t0 values for all histograms.
|
||||
*
|
||||
* Searches for t0 in order of priority:
|
||||
* -# RUN block t0 specification
|
||||
* -# GLOBAL block t0 specification
|
||||
* -# Data file header t0 values
|
||||
* -# Automatic estimation (with warning)
|
||||
*
|
||||
* Also handles addt0 for addrun histograms.
|
||||
*
|
||||
* \param runData Raw run data for histogram access
|
||||
* \param globalBlock GLOBAL block with default t0 values
|
||||
* \param histoNo Vector of histogram indices to process
|
||||
* \return true if valid t0 found for all histograms, false on error
|
||||
*/
|
||||
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);
|
||||
|
||||
/**
|
||||
* \brief Determines valid data range (first/last good bins).
|
||||
*
|
||||
* Establishes analysis boundaries in order of priority:
|
||||
* -# RUN block data range specification
|
||||
* -# GLOBAL block data range specification
|
||||
* -# Automatic estimation (t0 + 10ns offset to end of histogram)
|
||||
*
|
||||
* Validates that range is sensible (start < end, within histogram bounds).
|
||||
*
|
||||
* \return true if valid data range determined, false on error
|
||||
*/
|
||||
virtual Bool_t GetProperDataRange();
|
||||
|
||||
/**
|
||||
* \brief Determines fit time range from MSR file settings.
|
||||
*
|
||||
* Extracts fit range in order of priority:
|
||||
* -# RUN block fit range (time-based or bin-based)
|
||||
* -# GLOBAL block fit range (time-based or bin-based)
|
||||
* -# Default to data range (fgb to lgb) with warning
|
||||
*
|
||||
* Converts bin-based specifications to time using t0 and time resolution.
|
||||
*
|
||||
* \param globalBlock GLOBAL block with default fit settings
|
||||
*/
|
||||
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);
|
||||
|
||||
/**
|
||||
* \brief Finds the dominant precession frequency in raw data.
|
||||
*
|
||||
* Performs Fourier transform on raw histogram data to identify the
|
||||
* main muon precession frequency:
|
||||
* -# Creates TH1F from data in good bin range
|
||||
* -# Applies strong apodization (windowing)
|
||||
* -# Computes power spectrum
|
||||
* -# Searches for maximum above 10 MHz (ignores DC component)
|
||||
*
|
||||
* The frequency is used to:
|
||||
* - Determine optimal N₀ estimation window (integer cycles)
|
||||
* - Calculate suggested "optimal packing" for user information
|
||||
*
|
||||
* \param data Raw histogram data vector
|
||||
* \return Maximum frequency in MHz, or 0 if not found
|
||||
*/
|
||||
virtual Double_t GetMainFrequency(PDoubleVector &data);
|
||||
|
||||
/**
|
||||
* \brief Estimates initial normalization N₀ from lifetime-corrected data.
|
||||
*
|
||||
* Calculates N₀ as weighted average of M(t) over initial time window:
|
||||
* \f[
|
||||
* N_0 = \frac{\sum_{i=0}^{n} M(t_i)}{n}
|
||||
* \f]
|
||||
*
|
||||
* The window length is chosen to include an integer number of complete
|
||||
* oscillation cycles (based on freqMax) within fN0EstimateEndTime.
|
||||
*
|
||||
* Error estimate:
|
||||
* \f[
|
||||
* \sigma_{N_0} = \frac{\sqrt{\sum w_i^2 \sigma_{M_i}^2}}{\sum w_i}
|
||||
* \f]
|
||||
*
|
||||
* \param errN0 [out] Estimated error on N₀
|
||||
* \param freqMax Main precession frequency (MHz) for window calculation
|
||||
* \return Estimated N₀ value
|
||||
*/
|
||||
virtual Double_t EstimateN0(Double_t &errN0, Double_t freqMax);
|
||||
|
||||
/**
|
||||
* \brief Estimates background from pre-t0 bins.
|
||||
*
|
||||
* Calculates background level and error from specified background range:
|
||||
* -# Validates background range is within histogram bounds
|
||||
* -# Adjusts range to integer number of accelerator RF cycles (PSI/RAL/TRIUMF)
|
||||
* -# Computes mean and standard deviation
|
||||
* -# Stores in fBackground and fBkgErr
|
||||
*
|
||||
* \param histoNo Histogram index (for error messages)
|
||||
* \return true on success, false if background range invalid
|
||||
*/
|
||||
virtual Bool_t EstimateBkg(UInt_t histoNo);
|
||||
};
|
||||
|
||||
|
||||
Reference in New Issue
Block a user