365 lines
14 KiB
C++
365 lines
14 KiB
C++
/***************************************************************************
|
||
|
||
PRunBase.cpp
|
||
|
||
Author: Andreas Suter
|
||
e-mail: andreas.suter@psi.ch
|
||
|
||
***************************************************************************/
|
||
|
||
/***************************************************************************
|
||
* Copyright (C) 2007-2025 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. *
|
||
***************************************************************************/
|
||
|
||
#include <iostream>
|
||
|
||
#include <TROOT.h>
|
||
#include <TSystem.h>
|
||
#include <TString.h>
|
||
#include <TObjArray.h>
|
||
#include <TObjString.h>
|
||
#include <TFolder.h>
|
||
|
||
#include "PRunBase.h"
|
||
|
||
//--------------------------------------------------------------------------
|
||
// Constructor
|
||
//--------------------------------------------------------------------------
|
||
/**
|
||
* \brief Default constructor that initializes all member variables to default values.
|
||
*
|
||
* Creates an empty, invalid run object with all pointers set to nullptr and
|
||
* values set to undefined/invalid states. This constructor is needed to allow
|
||
* creation of vectors of PRunBase-derived objects.
|
||
*
|
||
* A run created with this constructor requires initialization via the main
|
||
* constructor before it can be used for fitting.
|
||
*/
|
||
PRunBase::PRunBase()
|
||
{
|
||
fRunNo = -1;
|
||
fRunInfo = nullptr;
|
||
fRawData = nullptr;
|
||
fMetaData.fField = PMUSR_UNDEFINED;
|
||
fMetaData.fEnergy = PMUSR_UNDEFINED;
|
||
fTimeResolution = -1.0;
|
||
|
||
fFitStartTime = PMUSR_UNDEFINED;
|
||
fFitEndTime = PMUSR_UNDEFINED;
|
||
|
||
fValid = true;
|
||
fHandleTag = kEmpty;
|
||
}
|
||
|
||
//--------------------------------------------------------------------------
|
||
// Constructor
|
||
//--------------------------------------------------------------------------
|
||
/**
|
||
* \brief Main constructor that initializes a run from MSR file and raw data.
|
||
*
|
||
* Performs comprehensive initialization:
|
||
* 1. Stores operation mode (fit vs. view) and pointers to MSR and data handlers
|
||
* 2. Extracts run-specific settings from the appropriate MSR RUN block
|
||
* 3. Validates function parameter mappings to ensure FUNCTIONS block is valid
|
||
* 4. Initializes metadata structures (field, energy, temperature)
|
||
* 5. Initializes function value cache for FUNCTIONS block evaluation
|
||
* 6. Creates PTheory object for evaluating the theory function
|
||
* 7. Validates that theory initialization succeeded
|
||
*
|
||
* If any initialization step fails (e.g., invalid theory, out-of-range parameters),
|
||
* the program exits with an error message. The run object is marked as valid upon
|
||
* successful completion.
|
||
*
|
||
* \param msrInfo Pointer to MSR file handler (must remain valid for object lifetime)
|
||
* \param rawData Pointer to raw data handler (must remain valid for object lifetime)
|
||
* \param runNo Run number (0-based index into MSR file RUN blocks)
|
||
* \param tag Operation mode: kFit (fitting), kView (display/plotting)
|
||
*/
|
||
PRunBase::PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag) :
|
||
fHandleTag(tag), fMsrInfo(msrInfo), fRawData(rawData)
|
||
{
|
||
fValid = true;
|
||
|
||
fRunNo = static_cast<Int_t>(runNo);
|
||
if (runNo > fMsrInfo->GetMsrRunList()->size()) {
|
||
fRunInfo = nullptr;
|
||
return;
|
||
}
|
||
|
||
// keep the run header info for this run
|
||
fRunInfo = &(*fMsrInfo->GetMsrRunList())[runNo];
|
||
|
||
// check the parameter and map range of the functions
|
||
if (!fMsrInfo->CheckMapAndParamRange(static_cast<UInt_t>(fRunInfo->GetMap()->size()), fMsrInfo->GetNoOfParams())) {
|
||
std::cerr << std::endl << "**SEVERE ERROR** PRunBase::PRunBase: map and/or parameter out of range in FUNCTIONS." << std::endl;
|
||
exit(0);
|
||
}
|
||
|
||
// init private variables
|
||
fMetaData.fField = PMUSR_UNDEFINED;
|
||
fMetaData.fEnergy = PMUSR_UNDEFINED;
|
||
fTimeResolution = -1.0;
|
||
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++)
|
||
fFuncValues.push_back(0.0);
|
||
|
||
// generate theory
|
||
fTheory = std::make_unique<PTheory>(fMsrInfo, runNo);
|
||
if (fTheory == nullptr) {
|
||
std::cerr << std::endl << "**SEVERE ERROR** PRunBase::PRunBase: Couldn't create an instance of PTheory :-(, will quit" << std::endl;
|
||
exit(0);
|
||
}
|
||
if (!fTheory->IsValid()) {
|
||
std::cerr << std::endl << "**SEVERE ERROR** PRunBase::PRunBase: Theory is not valid :-(, will quit" << std::endl;
|
||
exit(0);
|
||
}
|
||
|
||
// set fit time ranges
|
||
fFitStartTime = PMUSR_UNDEFINED;
|
||
fFitEndTime = PMUSR_UNDEFINED;
|
||
}
|
||
|
||
//--------------------------------------------------------------------------
|
||
// Destructor
|
||
//--------------------------------------------------------------------------
|
||
/**
|
||
* \brief Virtual destructor that cleans up allocated resources.
|
||
*
|
||
* Frees memory allocated for:
|
||
* - t0 value vectors (fT0s)
|
||
* - Additional run t0 vectors (fAddT0s)
|
||
* - Function value cache (fFuncValues)
|
||
*
|
||
* The PTheory object is automatically deleted via unique_ptr.
|
||
* Pointers to MSR handler and raw data handler are NOT deleted as they
|
||
* are owned by the calling code.
|
||
*/
|
||
PRunBase::~PRunBase()
|
||
{
|
||
fT0s.clear();
|
||
for (UInt_t i=0; i<fAddT0s.size(); i++)
|
||
fAddT0s[i].clear();
|
||
fAddT0s.clear();
|
||
|
||
fFuncValues.clear();
|
||
}
|
||
|
||
|
||
//--------------------------------------------------------------------------
|
||
// SetFitRange (public)
|
||
//--------------------------------------------------------------------------
|
||
/**
|
||
* \brief Sets the fit time range and recalculates the number of fitted bins.
|
||
*
|
||
* Updates the fitting window for this run. This method handles two scenarios:
|
||
* 1. Single fit range: If fitRange contains one pair, it applies to all runs
|
||
* 2. Multiple fit ranges: If fitRange contains multiple pairs, it selects the
|
||
* appropriate range based on fRunNo
|
||
*
|
||
* The method validates that:
|
||
* - The fit range vector is not empty (asserts)
|
||
* - The run number is within the fit range vector bounds
|
||
* - Start time is before end time (swaps if not)
|
||
*
|
||
* This is typically called by the FIT_RANGE command to dynamically adjust the
|
||
* fitting window during optimization or range scanning.
|
||
*
|
||
* \param fitRange Vector of (start, end) time pairs in microseconds (μs)
|
||
*
|
||
* Example: fitRange = {(0.1, 10.0), (0.2, 8.0)} applies first range to run 0,
|
||
* second range to run 1
|
||
*/
|
||
void PRunBase::SetFitRange(PDoublePairVector fitRange)
|
||
{
|
||
Double_t start=0.0, end=0.0;
|
||
|
||
assert(fitRange.size()); // make sure fitRange is not empty
|
||
|
||
if (fitRange.size()==1) { // one fit range for all
|
||
start = fitRange[0].first;
|
||
end = fitRange[0].second;
|
||
} else {
|
||
// check that fRunNo is found within fitRange
|
||
UInt_t idx=static_cast<UInt_t>(fRunNo);
|
||
if (idx < fitRange.size()) { // fRunNo present
|
||
start = fitRange[idx].first;
|
||
end = fitRange[idx].second;
|
||
} else { // fRunNo NOT present
|
||
std::cerr << std::endl << ">> PRunBase::SetFitRange(): **ERROR** msr-file run entry " << fRunNo << " not present in fit range vector.";
|
||
std::cerr << std::endl << ">> Will not do anything! Please check, this shouldn't happen." << std::endl;
|
||
return;
|
||
}
|
||
}
|
||
|
||
// check that start is before end
|
||
if (start > end) {
|
||
std::cerr << std::endl << ">> PRunBase::SetFitRange(): **WARNING** start=" << start << " is > as end=" << end;
|
||
std::cerr << std::endl << ">> Will swap them, since otherwise chisq/logLH == 0.0" << std::endl;
|
||
fFitStartTime = end;
|
||
fFitEndTime = start;
|
||
} else {
|
||
fFitStartTime = start;
|
||
fFitEndTime = end;
|
||
}
|
||
}
|
||
|
||
//--------------------------------------------------------------------------
|
||
// CleanUp (public)
|
||
//--------------------------------------------------------------------------
|
||
/**
|
||
* \brief Cleans up allocated resources.
|
||
*
|
||
* Releases memory used by the PTheory object via unique_ptr reset.
|
||
* This is called when the run processing is complete or when preparing
|
||
* for a new fit with different settings.
|
||
*
|
||
* Other data structures (vectors) are managed automatically by their destructors.
|
||
*/
|
||
void PRunBase::CleanUp()
|
||
{
|
||
fTheory.reset();
|
||
}
|
||
|
||
//--------------------------------------------------------------------------
|
||
// CalculateKaiserFilterCoeff (protected)
|
||
//--------------------------------------------------------------------------
|
||
/**
|
||
* \brief Calculates Kaiser window FIR filter coefficients for RRF smoothing.
|
||
*
|
||
* Designs a low-pass FIR filter using the Kaiser window method, which provides
|
||
* excellent control over the frequency response characteristics. The filter is
|
||
* used to smooth theory curves in rotating reference frame (RRF) fits, ensuring
|
||
* consistent comparison between filtered data and filtered theory.
|
||
*
|
||
* Algorithm (based on Oppenheim, Schafer, Buck, "Discrete-Time Signal Processing"):
|
||
* 1. Determine β parameter from attenuation requirement A
|
||
* 2. Calculate filter order m from transition width dw
|
||
* 3. Generate Kaiser window: \f$ w[n] = \frac{I_0(\beta\sqrt{1-(n/\alpha)^2})}{I_0(\beta)} \f$
|
||
* 4. Apply window to ideal sinc filter: \f$ h[n] = \frac{\sin(\omega_c n)}{\pi n} \cdot w[n] \f$
|
||
* 5. Normalize coefficients to unity gain at DC
|
||
*
|
||
* The β parameter is chosen based on attenuation A (in dB):
|
||
* - A > 50: β = 0.1102(A - 8.7)
|
||
* - 21 ≤ A ≤ 50: β = 0.5842(A - 21)^0.4 + 0.07886(A - 21)
|
||
* - A < 21: β = 0
|
||
*
|
||
* Filter order: m = (A - 8) / (2.285 × Δω × π), rounded to nearest odd integer
|
||
*
|
||
* Reference: A.V. Oppenheim, R.W. Schafer, J.R. Buck,
|
||
* "Discrete-Time Signal Processing", Pearson 2004, pp. 574ff
|
||
*
|
||
* \param wc Cutoff frequency (normalized, 0 to π rad/sample)
|
||
* \param A Attenuation in dB: A = -20 log₁₀(δ) where δ is the tolerance band
|
||
* \param dw Transition width: Δω = ω_S - ω_P (stop band - pass band frequencies)
|
||
*
|
||
* The computed coefficients are stored in fKaiserFilter and normalized for unity DC gain.
|
||
*
|
||
* \see FilterTheo() for application of these coefficients
|
||
*/
|
||
void PRunBase::CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw)
|
||
{
|
||
Double_t beta;
|
||
Double_t dt = fData.GetTheoryTimeStep();
|
||
Int_t m;
|
||
|
||
// estimate beta (see reference above, p.574ff)
|
||
if (A > 50.0) {
|
||
beta = 0.1102*(A-8.7);
|
||
} else if ((A >= 21.0) && (A <= 50.0)) {
|
||
beta = 0.5842*TMath::Power(A-21.0, 0.4) + 0.07886*(A-21.0);
|
||
} else {
|
||
beta = 0.0;
|
||
}
|
||
|
||
m = TMath::FloorNint((A-8.0)/(2.285*dw*TMath::Pi()));
|
||
// make sure m is odd
|
||
if (m % 2 == 0)
|
||
m++;
|
||
|
||
Double_t alpha = static_cast<Double_t>(m)/2.0;
|
||
Double_t dval;
|
||
Double_t dsum = 0.0;
|
||
for (Int_t i=0; i<=m; i++) {
|
||
dval = TMath::Sin(wc*(i-alpha)*dt)/(TMath::Pi()*(i-alpha)*dt);
|
||
dval *= TMath::BesselI0(beta*TMath::Sqrt(1.0-TMath::Power((i-alpha)*dt/alpha, 2.0)))/TMath::BesselI0(beta);
|
||
dsum += dval;
|
||
fKaiserFilter.push_back(dval);
|
||
}
|
||
for (UInt_t i=0; i<=static_cast<UInt_t>(m); i++) {
|
||
fKaiserFilter[i] /= dsum;
|
||
}
|
||
}
|
||
|
||
//--------------------------------------------------------------------------
|
||
// FilterTheo (protected)
|
||
//--------------------------------------------------------------------------
|
||
/**
|
||
* \brief Applies Kaiser FIR filter to theory values for RRF fits.
|
||
*
|
||
* Performs time-domain convolution of the theory function with the Kaiser filter
|
||
* coefficients computed by CalculateKaiserFilterCoeff(). This smooths the theory
|
||
* curve to match the filtering applied to RRF-transformed data, ensuring fair
|
||
* comparison during χ² calculation.
|
||
*
|
||
* The filtering operation is:
|
||
* \f[ y_{\rm filtered}[i] = \sum_{j=0}^{M-1} h[j] \cdot y_{\rm theory}[i-j] \f]
|
||
*
|
||
* where:
|
||
* - h[j] are the Kaiser filter coefficients from fKaiserFilter
|
||
* - M is the filter length
|
||
* - For i < j, the missing samples are treated as zero (causal filter)
|
||
*
|
||
* Additional processing:
|
||
* - The filtered theory replaces the original theory in fData
|
||
* - Time start is shifted backward by half the filter length to compensate
|
||
* for the group delay introduced by the symmetric FIR filter
|
||
*
|
||
* This method is only called by RRF-derived classes (PRunAsymmetryRRF, PRunSingleHistoRRF)
|
||
* after theory calculation and RRF transformation.
|
||
*
|
||
* \pre CalculateKaiserFilterCoeff() must have been called to initialize fKaiserFilter
|
||
* \pre fData must contain theory values (CalcTheory() must have been called)
|
||
*
|
||
* \see CalculateKaiserFilterCoeff() for filter coefficient calculation
|
||
*/
|
||
void PRunBase::FilterTheo()
|
||
{
|
||
PDoubleVector theoFiltered;
|
||
Double_t dval = 0.0;
|
||
const PDoubleVector *theo = fData.GetTheory();
|
||
|
||
for (UInt_t i=0; i<theo->size(); i++) {
|
||
for (UInt_t j=0; j<fKaiserFilter.size(); j++) {
|
||
if (i<j)
|
||
dval = 0.0;
|
||
else
|
||
dval += fKaiserFilter[j]*theo->at(i-j);
|
||
}
|
||
theoFiltered.push_back(dval);
|
||
dval = 0.0;
|
||
}
|
||
|
||
fData.ReplaceTheory(theoFiltered);
|
||
|
||
// shift time start by half the filter length
|
||
dval = fData.GetTheoryTimeStart() - 0.5*static_cast<Double_t>(fKaiserFilter.size())*fData.GetTheoryTimeStep();
|
||
fData.SetTheoryTimeStart(dval);
|
||
|
||
theoFiltered.clear();
|
||
}
|