445 lines
18 KiB
C++
445 lines
18 KiB
C++
/***************************************************************************
|
||
|
||
PRunSingleHisto.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 _PRUNSINGLEHISTO_H_
|
||
#define _PRUNSINGLEHISTO_H_
|
||
|
||
#include "PRunBase.h"
|
||
|
||
/**
|
||
* \brief Class for fitting single detector histograms (basic time-differential μSR).
|
||
*
|
||
* PRunSingleHisto implements the most fundamental μSR analysis: fitting a single
|
||
* positron detector histogram to extract relaxation parameters. This is the basis
|
||
* for all μSR measurements and is used when asymmetry analysis is not appropriate
|
||
* or desired.
|
||
*
|
||
* \section singlehisto_physics Physics and Applications
|
||
*
|
||
* <b>Single histogram measurements are used for:</b>
|
||
* - <b>Time-differential μSR:</b> Measure μ⁺ decay positron time spectrum
|
||
* - <b>Detector calibration:</b> Characterize individual detector response
|
||
* - <b>Transverse field (TF):</b> When forward/backward separation unnecessary
|
||
* - <b>High statistics:</b> Use all positrons (no F-B discrimination)
|
||
* - <b>Specialized geometries:</b> Non-standard detector arrangements
|
||
* - <b>Method development:</b> Test fitting algorithms on single detector
|
||
*
|
||
* \section singlehisto_data Data Structure and Analysis
|
||
*
|
||
* <b>Histogram content:</b>
|
||
* - Raw counts vs. time from a single positron detector
|
||
* - Time-zero (t0): Muon arrival time marking start of decay
|
||
* - Background: Constant or estimated from pre-t0 bins
|
||
* - Signal: N(t) = N₀·exp(-t/τ_μ)·P(t) + B
|
||
*
|
||
* where:
|
||
* - N₀ = initial count rate (normalization parameter)
|
||
* - τ_μ = 2.197 μs (muon lifetime)
|
||
* - P(t) = polarization function (contains physics: relaxation, oscillation, etc.)
|
||
* - B = background (random coincidences, accidentals)
|
||
*
|
||
* \section singlehisto_workflow Analysis Workflow
|
||
*
|
||
* 1. <b>Load Histogram:</b> Read raw detector counts from data file
|
||
* 2. <b>Determine t0:</b> Identify muon arrival time
|
||
* 3. <b>Background Subtraction:</b>
|
||
* - Fixed: Subtract specified constant
|
||
* - Estimated: Calculate from pre-t0 bins, subtract with error propagation
|
||
* 4. <b>Bin Packing:</b> Rebin to improve statistics (optional)
|
||
* 5. <b>Fit Range:</b> Select time window for parameter extraction
|
||
* 6. <b>Theory Evaluation:</b> N_theory(t) = N₀·exp(-t/τ_μ)·P_theory(t) + B
|
||
* 7. <b>Minimization:</b> χ² or maximum likelihood via MINUIT
|
||
*
|
||
* \section singlehisto_theory Theory Function
|
||
*
|
||
* The fitted function is typically:
|
||
* \f[ N(t) = N_0 e^{-t/\tau_\mu} P(t) + B \f]
|
||
*
|
||
* Common polarization functions P(t):
|
||
* - <b>Static Gaussian:</b> P(t) = exp(-σ²t²/2)
|
||
* - <b>Static Lorentzian:</b> P(t) = exp(-λt)
|
||
* - <b>Dynamic relaxation:</b> P(t) = exp(-(λt)^β) (stretched exponential)
|
||
* - <b>Oscillating:</b> P(t) = cos(ωt + φ) · exp(-λt)
|
||
* - <b>Kubo-Toyabe:</b> Complex relaxation functions for spin systems
|
||
*
|
||
* \section singlehisto_parameters Key Parameters
|
||
*
|
||
* <b>Normalization (N₀):</b>
|
||
* - Can be fit parameter or fixed value
|
||
* - Can be derived from FUNCTIONS block
|
||
* - Automatically scaled to 1/ns or 1/bin (fScaleN0AndBkg)
|
||
*
|
||
* <b>Background (B):</b>
|
||
* - Fixed: From "background" entry in RUN block
|
||
* - Estimated: Calculated from pre-t0 bin range
|
||
* - Units: counts/bin
|
||
*
|
||
* <b>Packing:</b>
|
||
* - Number of consecutive bins to combine
|
||
* - REQUIRED parameter (RUN or GLOBAL block)
|
||
* - Higher packing → better statistics, worse time resolution
|
||
*
|
||
* \section singlehisto_msr MSR File Example
|
||
*
|
||
* \code
|
||
* RUN data/run2425 PSI MUE4 PSI MUSR-ROOT (name beamline)
|
||
* fittype 0 (SingleHisto)
|
||
* map 1 (forward histogram number)
|
||
* forward 1
|
||
* packing 50 (combine 50 bins → one packed bin)
|
||
* background 50 150 (estimate from bins 50-150)
|
||
* data 200 2000 (use bins 200-2000 for analysis)
|
||
* t0 210.5 (muon arrival bin)
|
||
* fit 0.1 10.0 (fit from 0.1 to 10.0 μs after t0)
|
||
*
|
||
* THEORY
|
||
* asymmetry 1
|
||
* simpleGss 2 (Gaussian relaxation σ)
|
||
* + 3 (constant background offset)
|
||
* \endcode
|
||
*
|
||
* \section singlehisto_vs_asymmetry Single Histo vs. Asymmetry
|
||
*
|
||
* <table>
|
||
* <tr><th>Feature</th><th>Single Histogram</th><th>Asymmetry</th></tr>
|
||
* <tr><td>Detectors</td><td>One (forward)</td><td>Two (forward + backward)</td></tr>
|
||
* <tr><td>Quantity fitted</td><td>N(t) counts</td><td>A(t) = (F-αB)/(F+αB)</td></tr>
|
||
* <tr><td>Statistics</td><td>All positrons</td><td>F and B separately</td></tr>
|
||
* <tr><td>Background</td><td>Additive B</td><td>Cancels in asymmetry</td></tr>
|
||
* <tr><td>α parameter</td><td>N/A</td><td>Required (F/B asymmetry)</td></tr>
|
||
* <tr><td>Use cases</td><td>TF, calibration</td><td>LF, ZF, weak TF</td></tr>
|
||
* </table>
|
||
*
|
||
* \see PRunAsymmetry for forward-backward asymmetry analysis
|
||
* \see PRunSingleHistoRRF for rotating reference frame (high-TF)
|
||
* \see PRunBase for base class interface and common functionality
|
||
*/
|
||
class PRunSingleHisto : public PRunBase
|
||
{
|
||
public:
|
||
/**
|
||
* \brief Default constructor creating an empty, invalid single histogram run object.
|
||
*
|
||
* Initializes all member variables to default/invalid states:
|
||
* - fScaleN0AndBkg = true (scale to 1/ns by default)
|
||
* - fNoOfFitBins = 0 (no bins to fit)
|
||
* - fBackground = 0 (no background)
|
||
* - fPacking = -1 (unspecified - will cause error if not set)
|
||
* - fTheoAsData = false (high-resolution theory grid)
|
||
* - Good bins markers = -1 (unset)
|
||
* - Fit range bins = -1 (unset)
|
||
*
|
||
* This constructor is needed for creating vectors of PRunSingleHisto objects.
|
||
* The resulting object cannot be used until properly initialized via the main constructor.
|
||
*/
|
||
PRunSingleHisto();
|
||
|
||
/**
|
||
* \brief Main constructor initializing a single histogram run from MSR file and data.
|
||
*
|
||
* Performs comprehensive initialization:
|
||
*
|
||
* 1. <b>Base Class Initialization:</b>
|
||
* - Calls PRunBase constructor
|
||
* - Initializes theory engine, parameter mappings, FUNCTIONS block
|
||
*
|
||
* 2. <b>N₀/Background Scaling:</b>
|
||
* - Calls IsScaleN0AndBkg() to determine scaling mode
|
||
* - Sets fScaleN0AndBkg flag (true = scale to 1/ns, false = leave as 1/bin)
|
||
*
|
||
* 3. <b>Packing Validation (CRITICAL):</b>
|
||
* - Attempts to read packing from RUN block
|
||
* - Falls back to GLOBAL block if not in RUN
|
||
* - SEVERE ERROR if packing == -1 (mandatory parameter)
|
||
* - Marks run invalid and returns if packing not found
|
||
*
|
||
* 4. <b>Member Initialization:</b>
|
||
* - Good bin markers, fit range bins set to -1 (determined later)
|
||
* - Background initialized to 0 (set during PrepareData if specified)
|
||
*
|
||
* 5. <b>Data Preparation:</b>
|
||
* - Calls PrepareData() to load and preprocess histogram
|
||
* - If PrepareData() fails → marks run invalid
|
||
*
|
||
* The object is marked as invalid (fValid=false) if:
|
||
* - Packing parameter is missing from both RUN and GLOBAL blocks
|
||
* - PrepareData() fails (file not found, invalid t0, etc.)
|
||
*
|
||
* \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 Always check IsValid() after construction. Packing is MANDATORY.
|
||
*
|
||
* \see PrepareData() for data preprocessing details
|
||
* \see IsScaleN0AndBkg() for scaling determination
|
||
*/
|
||
PRunSingleHisto(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 ~PRunSingleHisto();
|
||
|
||
/**
|
||
* \brief Calculates χ² between histogram data and theory.
|
||
*
|
||
* Computes chi-squared for single histogram fitting:
|
||
* \f[ \chi^2 = \sum_{i} \frac{(N_i^{\rm data} - N_i^{\rm theory})^2}{\sigma_i^2} \f]
|
||
*
|
||
* where N_theory(t) = N₀·exp(-t/τ_μ)·P(t) + B
|
||
*
|
||
* Uses OpenMP parallelization when available. N₀ can be a fit parameter
|
||
* or derived from FUNCTIONS block.
|
||
*
|
||
* \param par Parameter vector from MINUIT
|
||
* \return χ² value (minimize during fitting)
|
||
*/
|
||
virtual Double_t CalcChiSquare(const std::vector<Double_t>& par);
|
||
|
||
/**
|
||
* \brief Calculates expected χ² based on theory predictions.
|
||
* \param par Parameter vector from MINUIT
|
||
* \return Expected χ² (for statistical diagnostics)
|
||
*/
|
||
virtual Double_t CalcChiSquareExpected(const std::vector<Double_t>& par);
|
||
|
||
/**
|
||
* \brief Calculates maximum likelihood for Poisson-distributed histogram counts.
|
||
*
|
||
* Computes -2ln(L) for low-count data (< 10-20 counts/bin).
|
||
* Superior to χ² when Gaussian approximation invalid.
|
||
*
|
||
* \param par Parameter vector from MINUIT
|
||
* \return -2×ln(L) value
|
||
*/
|
||
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
|
||
|
||
/**
|
||
* \brief Calculates expected maximum likelihood.
|
||
* \param par Parameter vector from MINUIT
|
||
* \return Expected -2×ln(L)
|
||
*/
|
||
virtual Double_t CalcMaxLikelihoodExpected(const std::vector<Double_t>& par);
|
||
|
||
/**
|
||
* \brief Evaluates theory function at all data points or high-resolution grid.
|
||
*
|
||
* Calculates N_theory(t) = N₀·exp(-t/τ_μ)·P_theory(t) + B
|
||
* using THEORY block functions. Stores results in fData for χ² calculation.
|
||
*/
|
||
virtual void CalcTheory();
|
||
|
||
/**
|
||
* \brief Returns the number of bins included in the fit range.
|
||
*
|
||
* Used for degrees of freedom: ν = N_bins - N_params
|
||
*
|
||
* \return Number of bins in fit range
|
||
*/
|
||
virtual UInt_t GetNoOfFitBins();
|
||
|
||
/**
|
||
* \brief Sets fit range using bin-offset specification (COMMANDS block syntax).
|
||
*
|
||
* Format: "fit_range fgb+n0 lgb-n1"
|
||
*
|
||
* \param fitRange String with bin offsets from good bin markers
|
||
*/
|
||
virtual void SetFitRangeBin(const TString fitRange);
|
||
|
||
/**
|
||
* \brief Returns the estimated background level.
|
||
* \return Background in counts/bin
|
||
*/
|
||
virtual Double_t GetBackground() { return fBackground; }
|
||
|
||
/**
|
||
* \brief Returns the first bin index in the fit range.
|
||
* \return Start bin index (0-based, after 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 Returns the bin packing factor.
|
||
* \return Number of raw bins combined into one packed bin
|
||
*/
|
||
virtual Int_t GetPacking() { return fPacking; }
|
||
|
||
/**
|
||
* \brief Returns the N₀/background scaling mode.
|
||
* \return true = scaled to 1/ns, false = left as 1/bin
|
||
*/
|
||
virtual Bool_t GetScaleN0AndBkg() { return fScaleN0AndBkg; }
|
||
|
||
/**
|
||
* \brief Calculates start/end bin indices from fit time range.
|
||
*
|
||
* Converts fit range (μs) to bin indices, accounting for t0, time resolution,
|
||
* and packing. Updates fStartTimeBin, fEndTimeBin, fNoOfFitBins.
|
||
*/
|
||
virtual void CalcNoOfFitBins();
|
||
|
||
protected:
|
||
/**
|
||
* \brief Main data preparation orchestrator.
|
||
*
|
||
* Coordinates histogram loading and preprocessing: determines operation mode,
|
||
* calls PrepareFitData() or PrepareViewData(), validates success.
|
||
*
|
||
* \return True if data preparation succeeds, false on error
|
||
*/
|
||
virtual Bool_t PrepareData();
|
||
|
||
/**
|
||
* \brief Prepares histogram data for fitting.
|
||
*
|
||
* Loads forward histogram, extracts metadata, determines t0, subtracts background,
|
||
* packs bins, propagates errors, sets up time grid and fit ranges.
|
||
*
|
||
* \param runData Raw run data handler
|
||
* \param histoNo Histogram index in data file
|
||
* \return True on success, false if preprocessing fails
|
||
*/
|
||
virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo);
|
||
|
||
/**
|
||
* \brief Prepares raw histogram data for viewing (minimal processing).
|
||
*
|
||
* Lighter-weight preprocessing for raw histogram visualization without
|
||
* background subtraction or full fitting infrastructure.
|
||
*
|
||
* \param runData Raw run data handler
|
||
* \param histoNo Histogram index
|
||
* \return True on success
|
||
*/
|
||
virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo);
|
||
|
||
/**
|
||
* \brief Prepares processed histogram data for viewing/plotting.
|
||
*
|
||
* Similar to PrepareFitData() but optimized for visualization with potentially
|
||
* wider time range for context.
|
||
*
|
||
* \param runData Raw run data handler
|
||
* \param histoNo Histogram index
|
||
* \return True on success
|
||
*/
|
||
virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo);
|
||
|
||
private:
|
||
Bool_t fScaleN0AndBkg; ///< Scaling mode: true = scale N₀ and B to 1/ns, false = leave as 1/bin (determined by IsScaleN0AndBkg())
|
||
UInt_t fNoOfFitBins; ///< Number of bins within fit range (fStartTimeBin to fEndTimeBin)
|
||
Double_t fBackground; ///< Background level in counts/bin (estimated from pre-t0 bins or fixed value from RUN block)
|
||
Int_t fPacking; ///< Bin packing factor (REQUIRED: from RUN or GLOBAL block)
|
||
Bool_t fTheoAsData; ///< Theory mode: true = at data points, false = high-resolution grid for smooth Fourier transforms
|
||
|
||
Int_t fGoodBins[2]; ///< Good bin markers for COMMANDS block: [0]=fgb (first good bin/t0), [1]=lgb (last good bin)
|
||
|
||
PDoubleVector fForward; ///< Forward detector histogram (background-corrected, packed)
|
||
|
||
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 from RUN block, data file header, GLOBAL block, or
|
||
* automatic determination. Validates t0 is within histogram bounds.
|
||
*
|
||
* \param runData Raw run data
|
||
* \param globalBlock GLOBAL block settings
|
||
* \param histoNo Vector of histogram indices
|
||
* \return True if valid t0 found
|
||
*/
|
||
virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo);
|
||
|
||
/**
|
||
* \brief Determines data range (region of valid histogram data).
|
||
*
|
||
* Establishes start/end bins for analysis from RUN block "data" entry.
|
||
* Data range is typically wider than fit range.
|
||
*
|
||
* \return True if valid data range determined
|
||
*/
|
||
virtual Bool_t GetProperDataRange();
|
||
|
||
/**
|
||
* \brief Determines fit range from MSR file settings.
|
||
*
|
||
* Extracts fit time window from RUN or GLOBAL block "fit" entry.
|
||
* Format: time-based (μs) or bin-based (fgb+n0 lgb-n1).
|
||
*
|
||
* \param globalBlock GLOBAL block with default fit settings
|
||
*/
|
||
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock);
|
||
|
||
/**
|
||
* \brief Estimates initial normalization N₀ from histogram data.
|
||
*
|
||
* Calculates N₀ estimate from histogram amplitude, used as starting
|
||
* value if N₀ is a fit parameter.
|
||
*/
|
||
virtual void EstimateN0();
|
||
|
||
/**
|
||
* \brief Estimates background from pre-t0 bins.
|
||
*
|
||
* Calculates background average and error from specified bin range
|
||
* (typically before t0). Sets fBackground member.
|
||
*
|
||
* \param histoNo Histogram index
|
||
* \return True on success, false if background range invalid
|
||
*/
|
||
virtual Bool_t EstimateBkg(UInt_t histoNo);
|
||
|
||
/**
|
||
* \brief Determines if N₀ and background should be scaled to 1/ns.
|
||
*
|
||
* Checks time resolution and fitting preferences to decide scaling mode.
|
||
* Returns true for standard time bins (scale to 1/ns), false otherwise.
|
||
*
|
||
* \return True if scaling should be applied
|
||
*/
|
||
virtual Bool_t IsScaleN0AndBkg();
|
||
};
|
||
|
||
#endif // _PRUNSINGLEHISTO_H_
|