diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 8a91749ac..f889327b2 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -51,7 +51,11 @@ // Constructor //-------------------------------------------------------------------------- /** - *

Constructor + * \brief Default constructor that initializes all member variables. + * + * Sets all counters and indices to default/invalid values. This constructor + * creates an invalid instance that requires proper initialization via the + * main constructor. */ PRunAsymmetryBNMR::PRunAsymmetryBNMR() : PRunBase() { @@ -72,12 +76,28 @@ PRunAsymmetryBNMR::PRunAsymmetryBNMR() : PRunBase() // Constructor //-------------------------------------------------------------------------- /** - *

Constructor + * \brief Main constructor that initializes β-NMR asymmetry fitting. * - * \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 + * Performs comprehensive initialization: + * 1. Validates packing parameter from RUN or GLOBAL block + * 2. Determines α/β parameter configuration (tags 1-6) + * 3. Validates α and β parameter numbers + * 4. Sets fAlphaBetaTag based on whether α/β are fixed to 1, free, or auto-estimated + * 5. Calls PrepareData() to load and process histogram data + * + * The α/β tag determines the asymmetry calculation method: + * - Tag 1: Both α=1 and β=1 (simplest case) + * - Tag 2: α free, β=1 (one asymmetry parameter) + * - Tag 3: α=1, β free (alternative single parameter) + * - Tag 4: Both α and β free (most general) + * - Tag 5: α auto-estimated, β=1 (automatic calibration) + * - Tag 6: α auto-estimated, β free + * + * \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 or kView) + * \param theoAsData If true, calculate theory only at data points */ PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) : PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData) @@ -166,7 +186,10 @@ PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawD // Destructor //-------------------------------------------------------------------------- /** - *

Destructor. + * \brief Destructor that cleans up helicity histogram data. + * + * Clears all eight histogram vectors (forward/backward × positive/negative helicity + * × data/errors) to free memory. */ PRunAsymmetryBNMR::~PRunAsymmetryBNMR() { @@ -185,12 +208,20 @@ PRunAsymmetryBNMR::~PRunAsymmetryBNMR() // CalcChiSquare (public) //-------------------------------------------------------------------------- /** - *

Calculate chi-square. + * \brief Calculates chi-square for β-NMR asymmetry fit. * - * return: - * - chisq value + * Computes χ² by comparing the asymmetry function with the theory: + * χ² = Σ[(A_data - A_theory)²/σ²] * - * \param par parameter vector iterated by minuit2 + * The asymmetry depends on fAlphaBetaTag: + * - Tag 1 (α=β=1): A = (F+ - B+) / (F+ + B+) + * - Tag 2-4: Various combinations with α and/or β corrections + * - Tag 5-6: Auto-estimated α + * + * Supports OpenMP parallelization for faster calculation. + * + * \param par Parameter vector from MINUIT minimizer + * \return Chi-square value */ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) { @@ -303,12 +334,13 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) // CalcChiSquareExpected (public) //-------------------------------------------------------------------------- /** - *

Calculate expected chi-square. Currently not implemented since not clear what to be done. + * \brief Calculates expected chi-square (not yet implemented). * - * return: - * - chisq value == 0.0 + * This method is intended for statistical analysis of fit quality, but + * implementation is pending due to complexity of β-NMR asymmetry calculation. * - * \param par parameter vector iterated by minuit2 + * \param par Parameter vector from MINUIT (unused) + * \return Always returns 0.0 (placeholder) */ Double_t PRunAsymmetryBNMR::CalcChiSquareExpected(const std::vector& par) { @@ -319,9 +351,13 @@ Double_t PRunAsymmetryBNMR::CalcChiSquareExpected(const std::vector& p // CalcMaxLikelihood (public) //-------------------------------------------------------------------------- /** - *

NOT IMPLEMENTED!! + * \brief Calculates maximum likelihood estimator (not yet implemented). * - * \param par parameter vector iterated by minuit2 + * Maximum likelihood fitting for β-NMR asymmetry is not yet implemented. + * Prints warning message when called. + * + * \param par Parameter vector from MINUIT (unused) + * \return Always returns 1.0 (placeholder) */ Double_t PRunAsymmetryBNMR::CalcMaxLikelihood(const std::vector& par) { @@ -334,9 +370,11 @@ Double_t PRunAsymmetryBNMR::CalcMaxLikelihood(const std::vector& par) // GetNoOfFitBins (public) //-------------------------------------------------------------------------- /** - *

Calculate the number of fitted bins for the current fit range. + * \brief Returns the number of bins used in the fit. * - * return: number of fitted bins. + * Calls CalcNoOfFitBins() to update fNoOfFitBins before returning the value. + * + * \return Number of bins included in the fit range */ UInt_t PRunAsymmetryBNMR::GetNoOfFitBins() { @@ -349,15 +387,22 @@ UInt_t PRunAsymmetryBNMR::GetNoOfFitBins() // SetFitRangeBin (public) //-------------------------------------------------------------------------- /** - *

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 Dynamically changes the fit range (used in COMMAND block). * - *

nXY are offsets which can be used to shift, limit the fit range. + * Parses a fit range string to update fit boundaries. The string format is: + * - Single range: FIT_RANGE fgb[+n0] lgb[-n1] + * - Multiple runs: FIT_RANGE fgb[+n00] lgb[-n01] fgb[+n10] lgb[-n11] ... * - * \param fitRange string containing the necessary information. + * Where: + * - fgb = first good bin + * - lgb = last good bin + * - +n = positive offset to shift fgb forward + * - -n = negative offset to shift lgb backward + * + * If one pair is given, it applies to all runs. If multiple pairs are given, + * the number must match the number of RUN blocks. + * + * \param fitRange String containing fit range specification */ void PRunAsymmetryBNMR::SetFitRangeBin(const TString fitRange) { @@ -439,7 +484,11 @@ void PRunAsymmetryBNMR::SetFitRangeBin(const TString fitRange) // CalcNoOfFitBins (public) //-------------------------------------------------------------------------- /** - *

Calculate the number of fitted bins for the current fit range. + * \brief Calculates the number of bins in the fit range. + * + * Converts fit time range (fFitStartTime, fFitEndTime) to bin indices + * (fStartTimeBin, fEndTimeBin) and calculates fNoOfFitBins. Performs + * boundary checking to ensure indices stay within valid data range. */ void PRunAsymmetryBNMR::CalcNoOfFitBins() { @@ -461,7 +510,19 @@ void PRunAsymmetryBNMR::CalcNoOfFitBins() // CalcTheory (protected) //-------------------------------------------------------------------------- /** - *

Calculate theory for a given set of fit-parameters. + * \brief Calculates theoretical β-NMR asymmetry values. + * + * Computes theory points for all data bins using the current parameters. + * The calculation method depends on fAlphaBetaTag: + * + * - Tag 1 (α=β=1): A = f(t) + * - Tag 2 (α≠1, β=1): A = [f(α+1)-(α-1)]/[(α+1)-f(α-1)] - [-f(α+1)-(α-1)]/[(α+1)+f(α-1)] + * - Tag 3 (α=1, β≠1): A = f(β+1)/[2-f(β-1)] - f(β+1)/[2+f(β-1)] + * - Tag 4 (α≠1, β≠1): Combined formula with both α and β corrections + * - Tag 5 (α estimated, β=1): Uses auto-estimated α + * - Tag 6 (α estimated, β≠1): Uses auto-estimated α with β correction + * + * where f(t) is the theory function from the FIT block. */ void PRunAsymmetryBNMR::CalcTheory() { @@ -567,23 +628,25 @@ void PRunAsymmetryBNMR::CalcTheory() // PrepareData (protected) //-------------------------------------------------------------------------- /** - *

Prepare data for fitting or viewing. What is already processed at this stage: - * - get all needed forward/backward histograms - * - get time resolution - * - get start/stop fit time - * - 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) - * - subtract background + * \brief Prepares β-NMR asymmetry data for fitting or viewing. * - * Error propagation for \f$ A_i = (f_i^{\rm c}-b_i^{\rm c})/(f_i^{\rm c}+b_i^{\rm c})\f$: - * \f[ \Delta A_i = \pm\frac{2}{(f_i^{\rm c}+b_i^{\rm c})^2}\left[ - * (b_i^{\rm c})^2 (\Delta f_i^{\rm c})^2 + - * (\Delta b_i^{\rm c})^2 (f_i^{\rm c})^2\right]^{1/2}\f] + * Main data preparation routine that performs the following steps: + * 1. Retrieves forward/backward histograms for both helicities + * 2. Gets time resolution from data file + * 3. Validates t0 values (cross-checks MSR file vs. data file) + * 4. Adds runs if addruns are specified + * 5. Groups histograms if grouping is specified + * 6. Subtracts background (fixed or estimated) + * 7. Applies bin packing if specified + * 8. Calculates β-NMR asymmetry with proper error propagation * - * return: - * - true if everthing went smooth - * - false, otherwise. + * The asymmetry calculation handles four histograms per helicity: + * - Forward+, Backward+, Forward-, Backward- + * + * Error propagation for asymmetry A = (F-B)/(F+B): + * \f[ \Delta A = \pm\frac{2}{(F+B)^2}\sqrt{B^2(\Delta F)^2 + F^2(\Delta B)^2} \f] + * + * \return True if successful, false on error */ Bool_t PRunAsymmetryBNMR::PrepareData() { @@ -1826,16 +1889,20 @@ Bool_t PRunAsymmetryBNMR::GetProperDataRange(PRawRunData* runData, UInt_t histoN // GetProperFitRange (private) //-------------------------------------------------------------------------- /** - *

Get the proper fit range. There are two possible fit range commands: - * fit 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 the proper fit range for the run. * - * \param globalBlock pointer to the GLOBAL block information form the msr-file. + * Fit ranges can be specified in two ways: + * - Time format: fit in microseconds + * - Bin format: fit fgb+offset_0 lgb-offset_1 in bins + * + * Resolution order: + * 1. Check RUN block for fit range (time or bins) + * 2. If not found, check GLOBAL block for fit range + * 3. If still not found, default to fgb/lgb (first/last good bins) + * + * For bin format, converts to time using: time = (bin - t0) × time_resolution + * + * \param globalBlock Pointer to GLOBAL block from MSR file */ void PRunAsymmetryBNMR::GetProperFitRange(PMsrGlobalBlock *globalBlock) { @@ -1875,9 +1942,17 @@ void PRunAsymmetryBNMR::GetProperFitRange(PMsrGlobalBlock *globalBlock) // EstimateAlpha (private) //-------------------------------------------------------------------------- /** - *

Get an estimate for alpha from the forward and backward histograms + * \brief Estimates the α parameter from histogram data. * - * \param globalBlock pointer to the GLOBAL block information form the msr-file. + * Calculates an automatic α value by comparing integrated counts in forward + * and backward histograms for both helicities. The estimate uses: + * + * α = √[(F+×B+) / (F-×B-)] + * + * where F+, B+, F-, B- are the total counts in the respective histograms. + * This provides a data-driven calibration when α is not explicitly specified. + * + * \return Estimated α value (defaults to 1.0 if calculation fails) */ Double_t PRunAsymmetryBNMR::EstimateAlpha() { diff --git a/src/include/PRunAsymmetryBNMR.h b/src/include/PRunAsymmetryBNMR.h index e173dc8ba..4a7641644 100644 --- a/src/include/PRunAsymmetryBNMR.h +++ b/src/include/PRunAsymmetryBNMR.h @@ -35,61 +35,204 @@ //--------------------------------------------------------------------------- /** - *

Class handling the asymmetry fit. + * \brief Class for handling β-NMR asymmetry fits. + * + * PRunAsymmetryBNMR implements asymmetry fitting for β-NMR (Beta-detected + * Nuclear Magnetic Resonance) experiments. Unlike conventional μSR asymmetry, + * β-NMR requires handling of helicity-dependent data with separate positive + * and negative helicity histograms. + * + * The asymmetry is calculated from four histograms: + * - Forward positive helicity (F+) + * - Backward positive helicity (B+) + * - Forward negative helicity (F-) + * - Backward negative helicity (B-) + * + * The class supports various α and β parameter configurations: + * - Tag 1: α = β = 1 (both fixed to unity) + * - Tag 2: α ≠ 1, β = 1 (free α, fixed β) + * - Tag 3: α = 1, β ≠ 1 (fixed α, free β) + * - Tag 4: α ≠ 1, β ≠ 1 (both free) + * - Tag 5: α auto-estimated, β = 1 + * - Tag 6: α auto-estimated, β ≠ 1 + * + * \see PRunBase for the base class providing common functionality + * \see PRunAsymmetry for the standard μSR asymmetry implementation */ class PRunAsymmetryBNMR : public PRunBase { public: + /// Default constructor PRunAsymmetryBNMR(); + + /** + * \brief Main constructor for β-NMR 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 + */ PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData); + + /// Destructor virtual ~PRunAsymmetryBNMR(); + /** + * \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& 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& par); + + /** + * \brief Calculates maximum likelihood estimator. + * \param par Parameter vector from MINUIT + * \return Maximum likelihood value + */ virtual Double_t CalcMaxLikelihood(const std::vector& par); + + /** + * \brief Calculates theoretical asymmetry function. + * + * Computes the theory values for the β-NMR 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 the four helicity 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]); private: - UInt_t fAlphaBetaTag; ///< \f$ 1 \to \alpha = \beta = 1\f$; \f$ 2 \to \alpha \neq 1, \beta = 1\f$; \f$ 3 \to \alpha = 1, \beta \neq 1\f$; \f$ 4 \to \alpha \neq 1, \beta \neq 1\f$. - UInt_t fNoOfFitBins; ///< number of bins to be be fitted - Int_t fPacking; ///< packing for this particular run. Either given in the RUN- or GLOBAL-block. - Bool_t fTheoAsData; ///< true=only calculate the theory points at the data points, false=calculate more points for the theory as compared to data are calculated which lead to 'nicer' Fouriers + UInt_t fAlphaBetaTag; ///< Tag indicating α/β configuration: 1=both unity, 2=α free/β unity, 3=α unity/β free, 4=both free, 5=α estimated/β unity, 6=α estimated/β 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 fForwardp; ///< pos hel forward histo data - PDoubleVector fForwardpErr; ///< pos hel forward histo errors - PDoubleVector fBackwardp; ///< pos hel backward histo data - PDoubleVector fBackwardpErr; ///< pos hel backward histo errors - PDoubleVector fForwardm; ///< neg hel forward histo data - PDoubleVector fForwardmErr; ///< neg hel forward histo errors - PDoubleVector fBackwardm; ///< neg hel backward histo data - PDoubleVector fBackwardmErr; ///< neg hel backward histo errors + PDoubleVector fForwardp; ///< Positive helicity forward histogram data + PDoubleVector fForwardpErr; ///< Positive helicity forward histogram errors + PDoubleVector fBackwardp; ///< Positive helicity backward histogram data + PDoubleVector fBackwardpErr; ///< Positive helicity backward histogram errors + PDoubleVector fForwardm; ///< Negative helicity forward histogram data + PDoubleVector fForwardmErr; ///< Negative helicity forward histogram errors + PDoubleVector fBackwardm; ///< Negative helicity backward histogram data + PDoubleVector fBackwardmErr; ///< Negative helicity backward histogram errors - Int_t fGoodBins[4]; ///< keep first/last good bins. 0=fgb, 1=lgb (forward); 2=fgb, 3=lgb (backward) + Int_t fGoodBins[4]; ///< Good bin boundaries: [0]=forward first, [1]=forward last, [2]=backward first, [3]=backward last - Int_t fStartTimeBin; ///< bin at which the fit starts - Int_t fEndTimeBin; ///< bin at which the fit ends + Int_t fStartTimeBin; ///< First bin index for fitting + Int_t fEndTimeBin; ///< Last bin index for fitting + /** + * \brief Subtracts fixed background from histograms. + * \return True on success, false on error + */ Bool_t SubtractFixBkg(); + + /** + * \brief Estimates and subtracts background from histograms. + * \return True on success, false on error + */ Bool_t SubtractEstimatedBkg(); + /** + * \brief Retrieves proper t0 values for all histograms. + * \param runData Pointer to raw run data + * \param globalBlock Pointer to global MSR block + * \param forwardHisto Vector of forward histogram indices + * \param backwardHistoNo Vector of backward histogram indices + * \return True on success, false on error + */ virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &forwardHisto, PUIntVector &backwardHistoNo); + + /** + * \brief Retrieves proper data range for histograms. + * \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. + * \param globalBlock Pointer to global MSR block + */ virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock); + + /** + * \brief Estimates α parameter from data. + * \return Estimated α value + * + * Calculates α based on the asymmetry ratio of forward and backward histograms. + */ virtual Double_t EstimateAlpha(); };