LMU/musrfit
5
0
Fork 0
Code Issues 5 Pull Requests Actions Packages Projects Releases 4 Activity
Files
8a863516740a3830f4dcf85cbceef6393866914e
musrfit/src/classes/PRunAsymmetryBNMR.cpp
Andreas Suter 4d5ad0a00c
All checks were successful
Build and Deploy Documentation / build-and-deploy (push) Successful in 17s
Details
improve doxygen documentation of PRunAsymmetryBNMR.*
2025-11-14 10:21:07 +01:00

1987 lines
82 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
/***************************************************************************
PRunAsymmetryBNMR.cpp
Author: Zaher Salman
Based on PRunAsymmetry.cpp by Andreas Suter
e-mail: zaher.salman@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2018-2025 by Zaher Salman *
* zaher.salman@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. *
***************************************************************************/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_GOMP
#include <omp.h>
#endif
#include <stdio.h>
#include <iostream>
#include <TString.h>
#include <TObjArray.h>
#include <TObjString.h>
#include "PMusr.h"
#include "PRunAsymmetryBNMR.h"
//--------------------------------------------------------------------------
// 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()
{
fNoOfFitBins = 0;
fPacking = -1;
fTheoAsData = false;
// the 2 following variables are need in case fit range is given in bins, and since
// the fit range can be changed in the command block, these variables need to be accessible
fGoodBins[0] = -1;
fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
}
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* \brief Main constructor that initializes β-NMR asymmetry fitting.
*
* 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)
{
// the 2 following variables are need in case fit range is given in bins, and since
// the fit range can be changed in the command block, these variables need to be accessible
fGoodBins[0] = -1;
fGoodBins[1] = -1;
fStartTimeBin = -1;
fEndTimeBin = -1;
fPacking = fRunInfo->GetPacking();
if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block
fPacking = fMsrInfo->GetMsrGlobal()->GetPacking();
}
if (fPacking == -1) { // this should NOT happen, somethin is severely wrong
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **SEVERE ERROR**: Couldn't find any packing information!";
std::cerr << std::endl << ">> This is very bad :-(, will quit ...";
std::cerr << std::endl;
fValid = false;
return;
}
// check if alpha and/or beta is fixed --------------------
PMsrParamList *param = msrInfo->GetMsrParamList();
// check if alpha is given
Bool_t alphaFixedToOne = false;
Bool_t alphaSetDefault = false;
if (fRunInfo->GetAlphaParamNo() == -1) { // no alpha given
// std::cerr << std::endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** no alpha parameter given! This is needed for an asymmetry fit!";
// std::cerr << std::endl;
// fValid = false;
// return;
alphaSetDefault = true;
} else if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > (Int_t)param->size())) { // check if alpha parameter is within proper bounds
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** alpha parameter no = " << fRunInfo->GetAlphaParamNo();
std::cerr << std::endl << ">> This is out of bound, since there are only " << param->size() << " parameters.";
std::cerr << std::endl;
fValid = false;
return;
} else { // check if alpha is fixed
if (((*param)[fRunInfo->GetAlphaParamNo()-1].fStep == 0.0) &&
((*param)[fRunInfo->GetAlphaParamNo()-1].fValue == 1.0))
alphaFixedToOne = true;
}
// check if beta is given
Bool_t betaFixedToOne = false;
if (fRunInfo->GetBetaParamNo() == -1) { // no beta given hence assuming beta == 1
betaFixedToOne = true;
} else if ((fRunInfo->GetBetaParamNo() < 0) || (fRunInfo->GetBetaParamNo() > (Int_t)param->size())) { // check if beta parameter is within proper bounds
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** beta parameter no = " << fRunInfo->GetBetaParamNo();
std::cerr << std::endl << ">> This is out of bound, since there are only " << param->size() << " parameters.";
std::cerr << std::endl;
fValid = false;
return;
} else { // check if beta is fixed
if (((*param)[fRunInfo->GetBetaParamNo()-1].fStep == 0.0) &&
((*param)[fRunInfo->GetBetaParamNo()-1].fValue == 1.0))
betaFixedToOne = true;
}
// set fAlphaBetaTag
if (alphaFixedToOne && betaFixedToOne) // alpha == 1, beta == 1
fAlphaBetaTag = 1;
else if (!alphaFixedToOne && betaFixedToOne & !alphaSetDefault) // alpha != 1, beta == 1
fAlphaBetaTag = 2;
else if (alphaFixedToOne && !betaFixedToOne) // alpha == 1, beta != 1
fAlphaBetaTag = 3;
else if (!alphaFixedToOne && betaFixedToOne & alphaSetDefault) // alpha ??, beta == 1
fAlphaBetaTag = 5;
else if (!alphaFixedToOne && !betaFixedToOne & alphaSetDefault) // alpha ??, beta != 1
fAlphaBetaTag = 6;
else
fAlphaBetaTag = 4;
// calculate fData
if (!PrepareData())
fValid = false;
}
//--------------------------------------------------------------------------
// 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()
{
fForwardp.clear();
fForwardpErr.clear();
fBackwardp.clear();
fBackwardpErr.clear();
fForwardm.clear();
fForwardmErr.clear();
fBackwardm.clear();
fBackwardmErr.clear();
}
//--------------------------------------------------------------------------
// CalcChiSquare (public)
//--------------------------------------------------------------------------
/**
* \brief Calculates chi-square for β-NMR asymmetry fit.
*
* Computes χ² by comparing the asymmetry function with the theory:
* χ² = Σ[(A_data - A_theory)²/σ²]
*
* 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<Double_t>& par)
{
Double_t chisq = 0.0;
Double_t diff = 0.0;
Double_t asymFcnValue = 0.0;
Double_t a, b, f;
// calculate functions
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
}
// calculate chi square
Double_t time(1.0),alphaest;
Int_t i;
// determine alpha/beta
alphaest = fRunInfo->GetEstimatedAlpha();
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
a = 1.0;
b = 1.0;
break;
case 2: // alpha != 1, beta == 1
if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
a = par[fRunInfo->GetAlphaParamNo()-1];
} else { // alpha is function
// get function number
UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
b = 1.0;
break;
case 3: // alpha == 1, beta != 1
a = 1.0;
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
b = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
break;
case 4: // alpha != 1, beta != 1
if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
a = par[fRunInfo->GetAlphaParamNo()-1];
} else { // alpha is function
// get function number
UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
b = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
break;
case 5: // alpha ?? , beta == 1
a = alphaest;
b = 1.0;
break;
case 6: // alpha ??, beta != 1
a = alphaest;
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
b = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
break;
default:
a = 1.0;
b = 1.0;
break;
}
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
// for a given set of parameters---which should be done outside of the parallelized loop.
// For all other functions it means a tiny and acceptable overhead.
asymFcnValue = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,f) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
f = fTheory->Func(time, par, fFuncValues)/2.0;
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0));
diff = fData.GetValue()->at(i) - asymFcnValue;
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
}
return chisq;
}
//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
/**
* \brief Calculates expected chi-square (not yet implemented).
*
* 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 from MINUIT (unused)
* \return Always returns 0.0 (placeholder)
*/
Double_t PRunAsymmetryBNMR::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
return 0.0;
}
//--------------------------------------------------------------------------
// CalcMaxLikelihood (public)
//--------------------------------------------------------------------------
/**
* \brief Calculates maximum likelihood estimator (not yet implemented).
*
* 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<Double_t>& par)
{
std::cout << std::endl << "PRunAsymmetryBNMR::CalcMaxLikelihood(): not implemented yet ..." << std::endl;
return 1.0;
}
//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
* \brief Returns the number of bins used in the fit.
*
* Calls CalcNoOfFitBins() to update fNoOfFitBins before returning the value.
*
* \return Number of bins included in the fit range
*/
UInt_t PRunAsymmetryBNMR::GetNoOfFitBins()
{
CalcNoOfFitBins();
return fNoOfFitBins;
}
//--------------------------------------------------------------------------
// SetFitRangeBin (public)
//--------------------------------------------------------------------------
/**
* \brief Dynamically changes the fit range (used in COMMAND block).
*
* 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] ...
*
* 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)
{
TObjArray *tok = nullptr;
TObjString *ostr = nullptr;
TString str;
Ssiz_t idx = -1;
Int_t offset = 0;
tok = fitRange.Tokenize(" \t");
if (tok->GetEntries() == 3) { // structure FIT_RANGE fgb+n0 lgb-n1
// handle fgb+n0 entry
ostr = dynamic_cast<TObjString*>(tok->At(1));
str = ostr->GetString();
// check if there is an offset present
idx = str.First("+");
if (idx != -1) { // offset present
str.Remove(0, idx+1);
if (str.IsFloat()) // if str is a valid number, convert is to an integer
offset = str.Atoi();
}
fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution;
// handle lgb-n1 entry
ostr = dynamic_cast<TObjString*>(tok->At(2));
str = ostr->GetString();
// check if there is an offset present
idx = str.First("-");
if (idx != -1) { // offset present
str.Remove(0, idx+1);
if (str.IsFloat()) // if str is a valid number, convert is to an integer
offset = str.Atoi();
}
fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution;
} else if ((tok->GetEntries() > 3) && (tok->GetEntries() % 2 == 1)) { // structure FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]]
Int_t pos = 2*(fRunNo+1)-1;
if (pos + 1 >= tok->GetEntries()) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
std::cerr << std::endl << ">> will ignore it. Sorry ..." << std::endl;
} else {
// handle fgb+n0 entry
ostr = static_cast<TObjString*>(tok->At(pos));
str = ostr->GetString();
// check if there is an offset present
idx = str.First("+");
if (idx != -1) { // offset present
str.Remove(0, idx+1);
if (str.IsFloat()) // if str is a valid number, convert is to an integer
offset = str.Atoi();
}
fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution;
// handle lgb-n1 entry
ostr = static_cast<TObjString*>(tok->At(pos+1));
str = ostr->GetString();
// check if there is an offset present
idx = str.First("-");
if (idx != -1) { // offset present
str.Remove(0, idx+1);
if (str.IsFloat()) // if str is a valid number, convert is to an integer
offset = str.Atoi();
}
fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution;
}
} else { // error
std::cerr << std::endl << ">> PRunAsymmetryBNMR::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
std::cerr << std::endl << ">> will ignore it. Sorry ..." << std::endl;
}
// clean up
if (tok) {
delete tok;
}
}
//--------------------------------------------------------------------------
// CalcNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
* \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()
{
// In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly
fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (fStartTimeBin < 0)
fStartTimeBin = 0;
fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
fEndTimeBin = fData.GetValue()->size();
if (fEndTimeBin > fStartTimeBin)
fNoOfFitBins = static_cast<UInt_t>(fEndTimeBin - fStartTimeBin);
else
fNoOfFitBins = 0;
}
//--------------------------------------------------------------------------
// CalcTheory (protected)
//--------------------------------------------------------------------------
/**
* \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()
{
// feed the parameter vector
std::vector<Double_t> par;
PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
for (UInt_t i=0; i<paramList->size(); i++)
par.push_back((*paramList)[i].fValue);
// calculate functions
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
}
// calculate asymmetry
Double_t asymFcnValue = 0.0;
Double_t a, b, f, alphaest;
Double_t time;
// Get estimated alpha
alphaest = fRunInfo->GetEstimatedAlpha();
for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(time, par, fFuncValues);
break;
case 2: // alpha != 1, beta == 1
if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
a = par[fRunInfo->GetAlphaParamNo()-1];
} else { // alpha is function
// get function number
UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0));
break;
case 3: // alpha == 1, beta != 1
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
b = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-f*(b+1.0)/(2.0+f*(b-1.0));
break;
case 4: // alpha != 1, beta != 1
if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
a = par[fRunInfo->GetAlphaParamNo()-1];
} else { // alpha is function
// get function number
UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
b = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0));
break;
case 5: // alpha ?? , beta == 1
a = alphaest;
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0));
break;
case 6: // alpha ??, beta != 1
a = alphaest;
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
b = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0));
break;
default:
asymFcnValue = 0.0;
break;
}
fData.AppendTheoryValue(asymFcnValue);
}
// clean up
par.clear();
}
//--------------------------------------------------------------------------
// PrepareData (protected)
//--------------------------------------------------------------------------
/**
* \brief Prepares β-NMR asymmetry data for fitting or viewing.
*
* 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
*
* 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()
{
if (!fValid)
return false;
// keep the Global block info
PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal();
// get forward/backward histo from PRunDataHandler object ------------------------
// get the correct run
PRawRunData *runData = fRawData->GetRunData(*(fRunInfo->GetRunName()));
if (!runData) { // run not found
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareData(): **ERROR** Couldn't get run " << fRunInfo->GetRunName()->Data() << "!";
std::cerr << std::endl;
return false;
}
// keep the field from the meta-data from the data-file
fMetaData.fField = runData->GetField();
// keep the energy from the meta-data from the data-file
fMetaData.fEnergy = runData->GetEnergy();
// keep the temperature(s) from the meta-data from the data-file
for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
fMetaData.fTemp.push_back(runData->GetTemperature(i));
// collect histogram numbers
PUIntVector forwardHistoNo;
PUIntVector backwardHistoNo;
for (UInt_t i=0; i<fRunInfo->GetForwardHistoNoSize(); i++) {
forwardHistoNo.push_back(fRunInfo->GetForwardHistoNo(i));
if (!runData->IsPresent(forwardHistoNo[i])) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareData(): **PANIC ERROR**:";
std::cerr << std::endl << ">> forwardHistoNo found = " << forwardHistoNo[i] << ", which is NOT present in the data file!?!?";
std::cerr << std::endl << ">> Will quit :-(";
std::cerr << std::endl;
// clean up
forwardHistoNo.clear();
backwardHistoNo.clear();
return false;
}
}
for (UInt_t i=0; i<fRunInfo->GetBackwardHistoNoSize(); i++) {
backwardHistoNo.push_back(fRunInfo->GetBackwardHistoNo(i));
if (!runData->IsPresent(backwardHistoNo[i])) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareData(): **PANIC ERROR**:";
std::cerr << std::endl << ">> backwardHistoNo found = " << backwardHistoNo[i] << ", which is NOT present in the data file!?!?";
std::cerr << std::endl << ">> Will quit :-(";
std::cerr << std::endl;
// clean up
forwardHistoNo.clear();
backwardHistoNo.clear();
return false;
}
}
/* //as35
if (forwardHistoNo.size() != backwardHistoNo.size()) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareData(): **PANIC ERROR**:";
std::cerr << std::endl << ">> # of forward histograms different from # of backward histograms.";
std::cerr << std::endl << ">> Will quit :-(";
std::cerr << std::endl;
// clean up
forwardHistoNo.clear();
backwardHistoNo.clear();
return false;
}
*/ //as35
// keep the time resolution in (s)
fTimeResolution = runData->GetTimeResolution()/1.0e3;
std::cout.precision(10);
std::cout << std::endl << ">> PRunAsymmetryBNMR::PrepareData(): time resolution=" << std::fixed << runData->GetTimeResolution() << "(ms)" << std::endl;
// get all the proper t0's and addt0's for the current RUN block
if (!GetProperT0(runData, globalBlock, forwardHistoNo, backwardHistoNo)) {
return false;
}
// keep the histo of each group at this point (addruns handled below)
std::vector<PDoubleVector> forward, backward;
forward.resize(forwardHistoNo.size()); // resize to number of groups
for (UInt_t i=0; i<forwardHistoNo.size(); i++) {
forward[i].resize(runData->GetDataBin(forwardHistoNo[i])->size());
forward[i] = *runData->GetDataBin(forwardHistoNo[i]);
}
backward.resize(backwardHistoNo.size()); // resize to number of groups
for (UInt_t i=0; i<backwardHistoNo.size(); i++) {
backward[i].resize(runData->GetDataBin(backwardHistoNo[i])->size());
backward[i] = *runData->GetDataBin(backwardHistoNo[i]);
}
// check if addrun's are present, and if yes add data
// check if there are runs to be added to the current one
if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present
PRawRunData *addRunData;
for (UInt_t i=1; i<fRunInfo->GetRunNameSize(); i++) {
// get run to be added to the main one
addRunData = fRawData->GetRunData(*(fRunInfo->GetRunName(i)));
if (addRunData == nullptr) { // couldn't get run
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareData(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!";
std::cerr << std::endl;
return false;
}
// add forward run
UInt_t addRunSize;
for (UInt_t k=0; k<forwardHistoNo.size(); k++) { // fill each group
addRunSize = addRunData->GetDataBin(forwardHistoNo[k])->size();
for (UInt_t j=0; j<addRunData->GetDataBin(forwardHistoNo[k])->size(); j++) { // loop over the bin indices
// make sure that the index stays in the proper range
if ((static_cast<Int_t>(j)+static_cast<Int_t>(fAddT0s[i-1][2*k])-static_cast<Int_t>(fT0s[2*k]) >= 0) &&
(j+static_cast<Int_t>(fAddT0s[i-1][2*k])-static_cast<Int_t>(fT0s[2*k]) < addRunSize)) {
forward[k][j] += addRunData->GetDataBin(forwardHistoNo[k])->at(j+static_cast<Int_t>(fAddT0s[i-1][2*k])-static_cast<Int_t>(fT0s[2*k]));
}
}
}
// add backward run
for (UInt_t k=0; k<backwardHistoNo.size(); k++) { // fill each group
addRunSize = addRunData->GetDataBin(backwardHistoNo[k])->size();
for (UInt_t j=0; j<addRunData->GetDataBin(backwardHistoNo[k])->size(); j++) { // loop over the bin indices
// make sure that the index stays in the proper range
if ((static_cast<Int_t>(j)+static_cast<Int_t>(fAddT0s[i-1][2*k+1])-static_cast<Int_t>(fT0s[2*k+1]) >= 0) &&
(j+static_cast<Int_t>(fAddT0s[i-1][2*k+1])-static_cast<Int_t>(fT0s[2*k+1]) < addRunSize)) {
backward[k][j] += addRunData->GetDataBin(backwardHistoNo[k])->at(j+static_cast<Int_t>(fAddT0s[i-1][2*k+1])-static_cast<Int_t>(fT0s[2*k+1]));
}
}
}
}
}
// set forward histo data of the first group
fForwardp.resize(forward[0].size());
fForwardm.resize(forward[0].size());
for (UInt_t i=0; i<fForwardp.size(); i++) {
fForwardp[i] = forward[0][i];
fForwardm[i] = forward[1][i];
}
// set backward histo data of the first group
fBackwardp.resize(backward[0].size());
fBackwardm.resize(backward[0].size());
for (UInt_t i=0; i<fBackwardp.size(); i++) {
fBackwardp[i] = backward[0][i];
fBackwardm[i] = backward[1][i];
}
// subtract background from histogramms ------------------------------------------
if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given
if (fRunInfo->GetBkgRange(0) >= 0) { // background range given
if (!SubtractEstimatedBkg())
return false;
} else { // no background given to do the job, try to estimate it
fRunInfo->SetBkgRange(static_cast<Int_t>(fT0s[0]*0.1), 0);
fRunInfo->SetBkgRange(static_cast<Int_t>(fT0s[0]*0.6), 1);
fRunInfo->SetBkgRange(static_cast<Int_t>(fT0s[1]*0.1), 2);
fRunInfo->SetBkgRange(static_cast<Int_t>(fT0s[1]*0.6), 3);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareData(): **WARNING** Neither fix background nor background bins are given!";
std::cerr << std::endl << ">> Will try the following:";
std::cerr << std::endl << ">> forward: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1);
std::cerr << std::endl << ">> backward: bkg start = " << fRunInfo->GetBkgRange(2) << ", bkg end = " << fRunInfo->GetBkgRange(3);
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS MAKES ANY SENSE! Better check ...";
std::cerr << std::endl;
if (!SubtractEstimatedBkg())
return false;
}
} else { // fixed background given
if (!SubtractFixBkg())
return false;
}
UInt_t histoNo[2] = {forwardHistoNo[0], backwardHistoNo[0]};
// get the data range (fgb/lgb) for the current RUN block
if (!GetProperDataRange(runData, histoNo)) {
return false;
}
// get the fit range for the current RUN block
GetProperFitRange(globalBlock);
// everything looks fine, hence fill data set
Bool_t status;
switch(fHandleTag) {
case kFit:
status = PrepareFitData();
break;
case kView:
status = PrepareViewData(runData, histoNo);
break;
default:
status = false;
break;
}
// clean up
forwardHistoNo.clear();
backwardHistoNo.clear();
return status;
}
//--------------------------------------------------------------------------
// SubtractFixBkg (private)
//--------------------------------------------------------------------------
/**
* <p>Subtracts a fixed background from the raw data. The background is given
* in units of (1/bin); for the Asymmetry representation (1/ns) doesn't make too much sense.
* The error propagation is done the following way: it is assumed that the error of the background
* is Poisson like, i.e. \f$\Delta\mathrm{bkg} = \sqrt{\mathrm{bkg}}\f$.
*
* Error propagation:
* \f[ \Delta f_i^{\rm c} = \pm\left[ (\Delta f_i)^2 + (\Delta \mathrm{bkg})^2 \right]^{1/2} =
* \pm\left[ f_i + \mathrm{bkg} \right]^{1/2}, \f]
* where \f$ f_i^{\rm c} \f$ is the background corrected histogram, \f$ f_i \f$ the raw histogram
* and \f$ \mathrm{bkg} \f$ the fix given background.
*
* <b>return:</b>
* - true
*
*/
Bool_t PRunAsymmetryBNMR::SubtractFixBkg()
{
Double_t dval;
// Order in RunInfo structure Fp, Fm, Bp, Bm
for (UInt_t i=0; i<fForwardp.size(); i++) {
// keep the error, and make sure that the bin is NOT empty
if (fForwardp[i] != 0.0)
dval = TMath::Sqrt(fForwardp[i]);
else
dval = 1.0;
fForwardpErr.push_back(dval);
fForwardp[i] -= fRunInfo->GetBkgFix(0);
// keep the error, and make sure that the bin is NOT empty
if (fForwardm[i] != 0.0)
dval = TMath::Sqrt(fForwardm[i]);
else
dval = 1.0;
fForwardmErr.push_back(dval);
fForwardm[i] -= fRunInfo->GetBkgFix(1);
// keep the error, and make sure that the bin is NOT empty
if (fBackwardp[i] != 0.0)
dval = TMath::Sqrt(fBackwardp[i]);
else
dval = 1.0;
fBackwardpErr.push_back(dval);
fBackwardp[i] -= fRunInfo->GetBkgFix(2);
// keep the error, and make sure that the bin is NOT empty
if (fBackwardm[i] != 0.0)
dval = TMath::Sqrt(fBackwardm[i]);
else
dval = 1.0;
fBackwardmErr.push_back(dval);
fBackwardm[i] -= fRunInfo->GetBkgFix(3);
}
return true;
}
//--------------------------------------------------------------------------
// SubtractEstimatedBkg (private)
//--------------------------------------------------------------------------
/**
* <p>Subtracts the background which is estimated from a given interval (typically before t0).
*
* The background corrected histogramms are:
* \f$ f_i^{\rm c} = f_i - \mathrm{bkg} \f$, where \f$ f_i \f$ is the raw data histogram,
* \f$ \mathrm{bkg} \f$ the background estimate, and \f$ f_i^{\rm c} \f$ background corrected
* histogram. The error on \f$ f_i^{\rm c} \f$ is
* \f[ \Delta f_i^{\rm c} = \pm \sqrt{ (\Delta f_i)^2 + (\Delta \mathrm{bkg})^2 } =
* \pm \sqrt{f_i + (\Delta \mathrm{bkg})^2} \f]
* The background error \f$ \Delta \mathrm{bkg} \f$ is
* \f[ \Delta \mathrm{bkg} = \pm\frac{1}{N}\left[\sum_{i=0}^N (\Delta f_i)^2\right]^{1/2} =
* \pm\frac{1}{N}\left[\sum_{i=0}^N f_i \right]^{1/2},\f]
* where \f$N\f$ is the number of bins over which the background is formed.
*
* <b>return:</b>
* - true
*/
Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg()
{
Double_t beamPeriod = 0.0;
// check if data are from PSI, RAL, or TRIUMF
if (fRunInfo->GetInstitute()->Contains("psi"))
beamPeriod = ACCEL_PERIOD_PSI;
else if (fRunInfo->GetInstitute()->Contains("ral"))
beamPeriod = ACCEL_PERIOD_RAL;
else if (fRunInfo->GetInstitute()->Contains("triumf"))
beamPeriod = ACCEL_PERIOD_TRIUMF;
else
beamPeriod = 0.0;
// check if start and end are in proper order
Int_t start[2] = {fRunInfo->GetBkgRange(0), fRunInfo->GetBkgRange(2)};
Int_t end[2] = {fRunInfo->GetBkgRange(1), fRunInfo->GetBkgRange(3)};
for (UInt_t i=0; i<2; i++) {
if (end[i] < start[i]) {
std::cout << std::endl << "PRunAsymmetryBNMR::SubtractEstimatedBkg(): end = " << end[i] << " > start = " << start[i] << "! Will swap them!";
UInt_t keep = end[i];
end[i] = start[i];
start[i] = keep;
}
}
// calculate proper background range
for (UInt_t i=0; i<2; i++) {
if (beamPeriod != 0.0) {
Double_t timeBkg = static_cast<Double_t>(end[i]-start[i])*(fTimeResolution*fPacking); // length of the background intervall in time
UInt_t fullCycles = static_cast<UInt_t>(timeBkg/beamPeriod); // how many proton beam cylces can be placed within the proposed background intervall
// correct the end of the background intervall such that the background is as close as possible to a multiple of the proton cylce
end[i] = start[i] + static_cast<UInt_t>((fullCycles*beamPeriod)/(fTimeResolution*fPacking));
std::cout << "PRunAsymmetryBNMR::SubtractEstimatedBkg(): Background " << start[i] << ", " << end[i] << std::endl;
if (end[i] == start[i])
end[i] = fRunInfo->GetBkgRange(2*i+1);
}
}
// check if start is within histogram bounds
if ((start[0] < 0) || (start[0] >= fForwardp.size()) ||
(start[1] < 0) || (start[1] >= fBackwardp.size())) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::SubtractEstimatedBkg(): **ERROR** background bin values out of bound!";
std::cerr << std::endl << ">> histo lengths (f/b) = (" << fForwardp.size() << "/" << fBackwardp.size() << ").";
std::cerr << std::endl << ">> background start (f/b) = (" << start[0] << "/" << start[1] << ").";
return false;
}
// check if end is within histogram bounds
if ((end[0] < 0) || (end[0] >= fForwardp.size()) ||
(end[1] < 0) || (end[1] >= fBackwardp.size())) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::SubtractEstimatedBkg(): **ERROR** background bin values out of bound!";
std::cerr << std::endl << ">> histo lengths (f/b) = (" << fForwardp.size() << "/" << fBackwardp.size() << ").";
std::cerr << std::endl << ">> background end (f/b) = (" << end[0] << "/" << end[1] << ").";
return false;
}
// calculate background
Double_t bkgp[2] = {0.0, 0.0};
Double_t errBkgp[2] = {0.0, 0.0};
Double_t bkgm[2] = {0.0, 0.0};
Double_t errBkgm[2] = {0.0, 0.0};
// forward
for (UInt_t i=start[0]; i<=end[0]; i++) {
bkgp[0] += fForwardp[i];
bkgm[0] += fForwardm[i];
}
errBkgp[0] = TMath::Sqrt(bkgp[0])/(end[0] - start[0] + 1);
bkgp[0] /= static_cast<Double_t>(end[0] - start[0] + 1);
std::cout << std::endl << ">> estimated pos hel forward histo background: " << bkgp[0];
errBkgm[0] = TMath::Sqrt(bkgp[0])/(end[0] - start[0] + 1);
bkgm[0] /= static_cast<Double_t>(end[0] - start[0] + 1);
std::cout << std::endl << ">> estimated neg hel forward histo background: " << bkgm[0];
// backward
for (UInt_t i=start[1]; i<=end[1]; i++) {
bkgp[1] += fBackwardp[i];
bkgm[1] += fBackwardm[i];
}
errBkgp[1] = TMath::Sqrt(bkgp[1])/(end[1] - start[1] + 1);
bkgp[1] /= static_cast<Double_t>(end[1] - start[1] + 1);
std::cout << std::endl << ">> estimated pos hel backward histo background: " << bkgp[1];
errBkgm[1] = TMath::Sqrt(bkgm[1])/(end[1] - start[1] + 1);
bkgm[1] /= static_cast<Double_t>(end[1] - start[1] + 1);
std::cout << std::endl << ">> estimated neg hel backward histo background: " << bkgm[1];
// correct error for forward, backward
Double_t errVal = 0.0;
for (UInt_t i=0; i<fForwardp.size(); i++) {
if (fForwardp[i] > 0.0)
errVal = TMath::Sqrt(fForwardp[i]+errBkgp[0]*errBkgp[0]);
else
errVal = 1.0;
fForwardpErr.push_back(errVal);
if (fBackwardp[i] > 0.0)
errVal = TMath::Sqrt(fBackwardp[i]+errBkgp[1]*errBkgp[1]);
else
errVal = 1.0;
fBackwardpErr.push_back(errVal);
if (fForwardm[i] > 0.0)
errVal = TMath::Sqrt(fForwardm[i]+errBkgm[0]*errBkgm[0]);
else
errVal = 1.0;
fForwardmErr.push_back(errVal);
if (fBackwardm[i] > 0.0)
errVal = TMath::Sqrt(fBackwardm[i]+errBkgm[1]*errBkgm[1]);
else
errVal = 1.0;
fBackwardmErr.push_back(errVal);
}
// subtract background from data
for (UInt_t i=0; i<fForwardp.size(); i++) {
fForwardp[i] -= bkgp[0];
fBackwardp[i] -= bkgp[1];
fForwardm[i] -= bkgm[0];
fBackwardm[i] -= bkgm[1];
}
fRunInfo->SetBkgEstimated(bkgp[0], 0);
fRunInfo->SetBkgEstimated(bkgp[1], 1);
fRunInfo->SetBkgEstimated(bkgm[0], 3);
fRunInfo->SetBkgEstimated(bkgm[1], 4);
// Get estimate for alpha once
Double_t alpha;
alpha = EstimateAlpha();
fRunInfo->SetEstimatedAlpha(alpha);
return true;
}
//--------------------------------------------------------------------------
// PrepareFitData (protected)
//--------------------------------------------------------------------------
/**
* <p>Take the pre-processed data (i.e. grouping and addrun are preformed) and form the asymmetry for fitting.
* Before forming the asymmetry, the following checks will be performed:
* -# check if data range is given, if not try to estimate one.
* -# check that if a data range is present, that it makes any sense.
* -# check that 'first good bin'-'t0' is the same for forward and backward histogram. If not adjust it.
* -# pack data (rebin).
* -# if packed forward size != backward size, truncate the longer one such that an asymmetry can be formed.
* -# calculate the asymmetry: \f$ A_i = (f_i^c-b_i^c)/(f_i^c+b_i^c) \f$
* -# calculate the asymmetry errors: \f$ \delta A_i = 2 \sqrt{(b_i^c)^2 (\delta f_i^c)^2 + (\delta b_i^c)^2 (f_i^c)^2}/(f_i^c+b_i^c)^2\f$
*/
Bool_t PRunAsymmetryBNMR::PrepareFitData()
{
// transform raw histo data. This is done the following way (for details see the manual):
// first rebin the data, than calculate the asymmetry
// everything looks fine, hence fill packed forward and backward histo
PRunData forwardpPacked;
PRunData backwardpPacked;
PRunData forwardmPacked;
PRunData backwardmPacked;
Double_t valuep = 0.0;
Double_t errorp = 0.0;
Double_t valuem = 0.0;
Double_t errorm = 0.0;
// forward
for (Int_t i=fGoodBins[0]; i<fGoodBins[1]; i++) {
if (fPacking == 1) {
forwardpPacked.AppendValue(fForwardp[i]);
forwardpPacked.AppendErrorValue(fForwardpErr[i]);
forwardmPacked.AppendValue(fForwardm[i]);
forwardmPacked.AppendErrorValue(fForwardmErr[i]);
} else { // packed data, i.e. fPacking > 1
if (((i-fGoodBins[0]) % fPacking == 0) && (i != fGoodBins[0])) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per bin
valuep /= fPacking;
valuem /= fPacking;
forwardpPacked.AppendValue(valuep);
forwardmPacked.AppendValue(valuem);
if (valuep == 0.0) {
forwardpPacked.AppendErrorValue(1.0);
} else {
forwardpPacked.AppendErrorValue(TMath::Sqrt(errorp)/fPacking);
}
if (valuem == 0.0) {
forwardmPacked.AppendErrorValue(1.0);
} else {
forwardmPacked.AppendErrorValue(TMath::Sqrt(errorm)/fPacking);
}
valuep = 0.0;
errorp = 0.0;
valuem = 0.0;
errorm = 0.0;
}
valuep += fForwardp[i];
errorp += fForwardpErr[i]*fForwardpErr[i];
valuem += fForwardm[i];
errorm += fForwardmErr[i]*fForwardmErr[i];
}
}
// backward
for (Int_t i=fGoodBins[2]; i<fGoodBins[3]; i++) {
if (fPacking == 1) {
backwardpPacked.AppendValue(fBackwardp[i]);
backwardpPacked.AppendErrorValue(fBackwardpErr[i]);
backwardmPacked.AppendValue(fBackwardm[i]);
backwardmPacked.AppendErrorValue(fBackwardmErr[i]);
} else { // packed data, i.e. fPacking > 1
if (((i-fGoodBins[2]) % fPacking == 0) && (i != fGoodBins[2])) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per bin
valuep /= fPacking;
valuem /= fPacking;
backwardpPacked.AppendValue(valuep);
backwardmPacked.AppendValue(valuem);
if (valuep == 0.0) {
backwardpPacked.AppendErrorValue(1.0);
} else {
backwardpPacked.AppendErrorValue(TMath::Sqrt(errorp)/fPacking);
}
if (valuem == 0.0) {
backwardmPacked.AppendErrorValue(1.0);
} else {
backwardmPacked.AppendErrorValue(TMath::Sqrt(errorm)/fPacking);
}
valuep = 0.0;
errorp = 0.0;
valuem = 0.0;
errorm = 0.0;
}
valuep += fBackwardp[i];
errorp += fBackwardpErr[i]*fBackwardpErr[i];
valuem += fBackwardm[i];
errorm += fBackwardmErr[i]*fBackwardmErr[i];
}
}
// check if packed forward and backward hist have the same size, otherwise take the minimum size
UInt_t noOfBins = forwardpPacked.GetValue()->size();
if (forwardpPacked.GetValue()->size() != backwardpPacked.GetValue()->size()) {
if (forwardpPacked.GetValue()->size() > backwardpPacked.GetValue()->size())
noOfBins = backwardpPacked.GetValue()->size();
}
// form asymmetry including error propagation
Double_t asym,error;
Double_t fp, bp, efp, ebp;
Double_t fm, bm, efm, ebm;
// fill data time start, and step
// data start at data_start-t0 shifted by (pack-1)/2
fData.SetDataTimeStart(fTimeResolution*(static_cast<Double_t>(fGoodBins[0])-fT0s[0]+static_cast<Double_t>(fPacking-1)/2.0));
fData.SetDataTimeStep(fTimeResolution*static_cast<Double_t>(fPacking));
for (UInt_t i=0; i<noOfBins; i++) {
// to make the formulae more readable
fp = forwardpPacked.GetValue()->at(i);
bp = backwardpPacked.GetValue()->at(i);
efp = forwardpPacked.GetError()->at(i);
ebp = backwardpPacked.GetError()->at(i);
fm = forwardmPacked.GetValue()->at(i);
bm = backwardmPacked.GetValue()->at(i);
efm = forwardmPacked.GetError()->at(i);
ebm = backwardmPacked.GetError()->at(i);
// check that there are indeed bins
if (fp+bp != 0.0)
asym = (fp-bp) / (fp+bp) - (fm-bm) / (fm+bm);
else
asym = 0.0;
fData.AppendValue(asym);
// calculate the error
if (fp+bp != 0.0)
errorp = 2.0/((fp+bp)*(fp+bp))*TMath::Sqrt(bp*bp*efp*efp+ebp*ebp*fp*fp);
else
errorp = 1.0;
if (fm+bm != 0.0)
errorm = 2.0/((fm+bm)*(fm+bm))*TMath::Sqrt(bm*bm*efm*efm+ebm*ebm*fm*fm);
else
errorp = 1.0;
error = TMath::Sqrt(errorp*errorp+errorm*errorm);
fData.AppendErrorValue(error);
}
CalcNoOfFitBins();
// clean up
fForwardp.clear();
fForwardpErr.clear();
fBackwardp.clear();
fBackwardpErr.clear();
fForwardm.clear();
fForwardmErr.clear();
fBackwardm.clear();
fBackwardmErr.clear();
return true;
}
//--------------------------------------------------------------------------
/**
* <p>Take the pre-processed data (i.e. grouping and addrun are preformed) and form the asymmetry for view representation.
* -# check that data range is present, that it makes any sense.
* -# check that 'first good bin'-'t0' is the same for forward and backward histogram. If not adjust it.
* -# pack data (rebin).
* -# if packed forward size != backward size, truncate the longer one such that an asymmetry can be formed.
* -# calculate the asymmetry: \f$ A_i = (\alpha f_i^c-b_i^c)/(\alpha \beta f_i^c+b_i^c) \f$
* -# calculate the asymmetry errors: \f$ \delta A_i = 2 \sqrt{(b_i^c)^2 (\delta f_i^c)^2 + (\delta b_i^c)^2 (f_i^c)^2}/(f_i^c+b_i^c)^2\f$
* -# calculate the theory vector.
*
* \param runData raw run data needed to perform some crosschecks
* \param histoNo histogram number (within a run). histoNo[0]: forward histogram number, histNo[1]: backward histogram number
*/
Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
{
// check if view_packing is wished
Int_t packing = fPacking;
if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) {
packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking;
}
// feed the parameter vector
std::vector<Double_t> par;
PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
for (UInt_t i=0; i<paramList->size(); i++)
par.push_back((*paramList)[i].fValue);
// transform raw histo data. This is done the following way (for details see the manual):
// first rebin the data, than calculate the asymmetry
// first get start data, end data, and t0
Int_t start[2] = {fGoodBins[0], fGoodBins[2]};
Int_t end[2] = {fGoodBins[1], fGoodBins[3]};
Int_t t0[2] = {static_cast<Int_t>(fT0s[0]), static_cast<Int_t>(fT0s[1])};
// check if the data ranges and t0's between forward/backward are compatible
Int_t fgb[2];
if (start[0]-t0[0] != start[1]-t0[1]) { // wrong fgb aligning
if (abs(start[0]-t0[0]) > abs(start[1]-t0[1])) {
fgb[0] = start[0];
fgb[1] = t0[1] + start[0]-t0[0];
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **WARNING** needed to shift backward fgb from ";
std::cerr << start[1] << " to " << fgb[1] << std::endl;
} else {
fgb[0] = t0[0] + start[1]-t0[1];
fgb[1] = start[1];
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **WARNING** needed to shift forward fgb from ";
std::cerr << start[0] << " to " << fgb[0] << std::endl;
}
} else { // fgb aligning is correct
fgb[0] = start[0];
fgb[1] = start[1];
}
Int_t val = fgb[0]-packing*(fgb[0]/packing);
do {
if (fgb[1] - fgb[0] < 0)
val += packing;
} while (val + fgb[1] - fgb[0] < 0);
start[0] = val;
start[1] = val + fgb[1] - fgb[0];
// make sure that there are equal number of rebinned bins in forward and backward
UInt_t noOfBins0 = (runData->GetDataBin(histoNo[0])->size()-start[0])/packing;
UInt_t noOfBins1 = (runData->GetDataBin(histoNo[1])->size()-start[1])/packing;
if (noOfBins0 > noOfBins1)
noOfBins0 = noOfBins1;
end[0] = start[0] + noOfBins0 * packing;
end[1] = start[1] + noOfBins0 * packing;
// check if start, end, and t0 make any sense
// 1st check if start and end are in proper order
for (UInt_t i=0; i<2; i++) {
if (end[i] < start[i]) { // need to swap them
Int_t keep = end[i];
end[i] = start[i];
start[i] = keep;
}
// 2nd check if start is within proper bounds
if ((start[i] < 0) || (start[i] > static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size()))) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **ERROR** start data bin doesn't make any sense!";
std::cerr << std::endl;
return false;
}
// 3rd check if end is within proper bounds
if ((end[i] < 0) || (end[i] > static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size()))) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **ERROR** end data bin doesn't make any sense!";
std::cerr << std::endl;
return false;
}
// 4th check if t0 is within proper bounds
if ((t0[i] < 0) || (t0[i] > static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size()))) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **ERROR** t0 data bin doesn't make any sense!";
std::cerr << std::endl;
return false;
}
}
// everything looks fine, hence fill packed forward and backward histo
PRunData forwardpPacked;
PRunData backwardpPacked;
PRunData forwardmPacked;
PRunData backwardmPacked;
Double_t valuep = 0.0;
Double_t errorp = 0.0;
Double_t valuem = 0.0;
Double_t errorm = 0.0;
Double_t value = 0.0;
Double_t error = 0.0;
// forward
for (Int_t i=start[0]; i<end[0]; i++) {
if (packing == 1) {
forwardpPacked.AppendValue(fForwardp[i]);
forwardpPacked.AppendErrorValue(fForwardpErr[i]);
forwardmPacked.AppendValue(fForwardm[i]);
forwardmPacked.AppendErrorValue(fForwardmErr[i]);
} else { // packed data, i.e. packing > 1
if (((i-start[0]) % packing == 0) && (i != start[0])) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per bin
valuep /= packing;
forwardpPacked.AppendValue(valuep);
valuem /= packing;
forwardmPacked.AppendValue(valuem);
if (valuep == 0.0) {
forwardpPacked.AppendErrorValue(1.0);
} else {
forwardpPacked.AppendErrorValue(TMath::Sqrt(errorp)/packing);
}
if (valuem == 0.0) {
forwardmPacked.AppendErrorValue(1.0);
} else {
forwardmPacked.AppendErrorValue(TMath::Sqrt(errorm)/packing);
}
valuep = 0.0;
errorp = 0.0;
valuem = 0.0;
errorm = 0.0;
}
valuep += fForwardp[i];
errorp += fForwardpErr[i]*fForwardpErr[i];
valuem += fForwardm[i];
errorm += fForwardmErr[i]*fForwardmErr[i];
}
}
// backward
for (Int_t i=start[1]; i<end[1]; i++) {
if (packing == 1) {
backwardpPacked.AppendValue(fBackwardp[i]);
backwardpPacked.AppendErrorValue(fBackwardpErr[i]);
backwardmPacked.AppendValue(fBackwardm[i]);
backwardmPacked.AppendErrorValue(fBackwardmErr[i]);
} else { // packed data, i.e. packing > 1
if (((i-start[1]) % packing == 0) && (i != start[1])) { // fill data
// in order that after rebinning the fit does not need to be redone (important for plots)
// the value is normalize to per bin
valuep /= packing;
valuem /= packing;
backwardpPacked.AppendValue(valuep);
backwardmPacked.AppendValue(valuem);
if (valuep == 0.0) {
backwardpPacked.AppendErrorValue(1.0);
} else {
backwardpPacked.AppendErrorValue(TMath::Sqrt(errorp)/packing);
}
if (valuem == 0.0) {
backwardmPacked.AppendErrorValue(1.0);
} else {
backwardmPacked.AppendErrorValue(TMath::Sqrt(errorm)/packing);
}
valuep = 0.0;
errorp = 0.0;
valuem = 0.0;
errorm = 0.0;
}
valuep += fBackwardp[i];
errorp += fBackwardpErr[i]*fBackwardpErr[i];
valuem += fBackwardm[i];
errorm += fBackwardmErr[i]*fBackwardmErr[i];
}
}
// check if packed forward and backward hist have the same size, otherwise take the minimum size
UInt_t noOfBins = forwardpPacked.GetValue()->size();
if (forwardpPacked.GetValue()->size() != backwardpPacked.GetValue()->size()) {
if (forwardpPacked.GetValue()->size() > backwardpPacked.GetValue()->size())
noOfBins = backwardpPacked.GetValue()->size();
}
// form asymmetry including error propagation
Double_t asym;
Double_t fp, bp, efp, ebp, alpha, beta = 1.0;
Double_t fm, bm, efm, ebm;
// set data time start, and step
// data start at data_start-t0
fData.SetDataTimeStart(fTimeResolution*(static_cast<Double_t>(start[0])-t0[0]+static_cast<Double_t>(packing-1)/2.0));
fData.SetDataTimeStep(fTimeResolution*static_cast<Double_t>(packing));
// Get estimate for alpha once
alpha = EstimateAlpha();
// get the proper alpha and beta
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
alpha = 1.0;
beta = 1.0;
break;
case 2: // alpha != 1, beta == 1
if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
alpha = par[fRunInfo->GetAlphaParamNo()-1];
} else { // alpha is function
// get function number
UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
beta = 1.0;
break;
case 3: // alpha == 1, beta != 1
alpha = 1.0;
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
beta = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
break;
case 4: // alpha != 1, beta != 1
if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter
alpha = par[fRunInfo->GetAlphaParamNo()-1];
} else { // alpha is function
// get function number
UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
beta = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
break;
case 5: // alpha ?? , beta == 1
// use estimated value
beta = 1.0;
break;
case 6: // alpha ??, beta != 1
// use estimated value
if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter
beta = par[fRunInfo->GetBetaParamNo()-1];
} else { // beta is a function
// get function number
UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
// evaluate function
beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par, fMetaData);
}
break;
default:
break;
}
for (UInt_t i=0; i<forwardpPacked.GetValue()->size(); i++) {
// to make the formulae more readable
fp = forwardpPacked.GetValue()->at(i);
bp = backwardpPacked.GetValue()->at(i);
efp = forwardpPacked.GetError()->at(i);
ebp = backwardpPacked.GetError()->at(i);
fm = forwardmPacked.GetValue()->at(i);
bm = backwardmPacked.GetValue()->at(i);
efm = forwardmPacked.GetError()->at(i);
ebm = backwardmPacked.GetError()->at(i);
// check that there are indeed bins
if (fp+bp != 0.0)
asym = (alpha*fp-bp) / (alpha*beta*fp+bp) - (alpha*fm-bm) / (alpha*beta*fm+bm);
else
asym = 0.0;
fData.AppendValue(asym);
// calculate the error
if (fp+bp != 0.0)
errorp = 2.0/((fp+bp)*(fp+bp))*TMath::Sqrt(bp*bp*efp*efp+ebp*ebp*fp*fp);
else
errorp = 1.0;
if (fm+bm != 0.0)
errorm = 2.0/((fm+bm)*(fm+bm))*TMath::Sqrt(bm*bm*efm*efm+ebm*ebm*fm*fm);
else
errorm = 1.0;
error = TMath::Sqrt(errorp*errorp+errorm*errorm);
fData.AppendErrorValue(error);
}
CalcNoOfFitBins();
// clean up
fForwardp.clear();
fForwardpErr.clear();
fBackwardp.clear();
fBackwardpErr.clear();
fForwardm.clear();
fForwardmErr.clear();
fBackwardm.clear();
fBackwardmErr.clear();
// fill theory vector for kView
// calculate functions
for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
}
// calculate theory
UInt_t size = runData->GetDataBin(histoNo[0])->size();
Int_t factor = 8; // 8 times more points for the theory (if fTheoAsData == false)
fData.SetTheoryTimeStart(fData.GetDataTimeStart());
if (fTheoAsData) { // calculate theory only at the data points
fData.SetTheoryTimeStep(fData.GetDataTimeStep());
} else {
// finer binning for the theory (8 times as many points = factor)
size *= factor;
fData.SetTheoryTimeStep(fData.GetDataTimeStep()/(Double_t)factor);
}
Double_t time;
for (UInt_t i=0; i<size; i++) {
time = fData.GetTheoryTimeStart() + static_cast<Double_t>(i)*fData.GetTheoryTimeStep();
value = fTheory->Func(time, par, fFuncValues);
if (fabs(value) > 10.0) { // dirty hack needs to be fixed!!
value = 0.0;
}
fData.AppendTheoryValue(value);
}
// clean up
par.clear();
return true;
}
//--------------------------------------------------------------------------
// 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
*
* \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 forwardHistoNo histogram number vector of forward; forwardHistoNo = msr-file forward + redGreen_offset - 1
* \param backwardHistoNo histogram number vector of backwardward; backwardHistoNo = msr-file backward + redGreen_offset - 1
*
* <b>return:</b>
* - true if everthing went smooth
* - false, otherwise.
*/
Bool_t PRunAsymmetryBNMR::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &forwardHistoNo, PUIntVector &backwardHistoNo)
{
// feed all T0's
// first init T0's, T0's are stored as (forward T0, backward T0, etc.)
fT0s.clear();
// this strange fT0 size estimate is needed in case #forw histos != #back histos
size_t size = 2*forwardHistoNo.size();
if (backwardHistoNo.size() > forwardHistoNo.size())
size = 2*backwardHistoNo.size();
fT0s.resize(size);
for (UInt_t i=0; i<fT0s.size(); i++) {
fT0s[i] = -1.0;
}
// fill in the T0's from the msr-file (if present)
for (UInt_t i=0; i<fRunInfo->GetT0BinSize(); i++) {
fT0s[i] = fRunInfo->GetT0Bin(i);
}
// fill in the missing T0's from the GLOBAL block section (if present)
for (UInt_t i=0; i<globalBlock->GetT0BinSize(); i++) {
if (fT0s[i] == -1) { // i.e. not given in the RUN block section
fT0s[i] = globalBlock->GetT0Bin(i);
}
}
// fill in the missing T0's from the data file, if not already present in the msr-file
for (UInt_t i=0; i<forwardHistoNo.size(); i++) {
if (fT0s[2*i] == -1.0) // i.e. not present in the msr-file, try the data file
if (runData->GetT0Bin(forwardHistoNo[i]) > 0.0) {
fT0s[2*i] = runData->GetT0Bin(forwardHistoNo[i]);
fRunInfo->SetT0Bin(fT0s[2*i], 2*i);
}
}
for (UInt_t i=0; i<backwardHistoNo.size(); i++) {
if (fT0s[2*i+1] == -1.0) // i.e. not present in the msr-file, try the data file
if (runData->GetT0Bin(backwardHistoNo[i]) > 0.0) {
fT0s[2*i+1] = runData->GetT0Bin(backwardHistoNo[i]);
fRunInfo->SetT0Bin(fT0s[2*i+1], 2*i+1);
}
}
// fill in the T0's gaps, i.e. in case the T0's are NEITHER in the msr-file and NOR in the data file
for (UInt_t i=0; i<forwardHistoNo.size(); i++) {
if (fT0s[2*i] == -1.0) { // i.e. not present in the msr-file and data file, use the estimated T0
fT0s[2*i] = runData->GetT0BinEstimated(forwardHistoNo[i]);
fRunInfo->SetT0Bin(fT0s[2*i], 2*i);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
std::cerr << std::endl << ">> run: " << fRunInfo->GetRunName()->Data();
std::cerr << std::endl << ">> will try the estimated one: forward t0 = " << runData->GetT0BinEstimated(forwardHistoNo[i]);
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
std::cerr << std::endl;
}
}
for (UInt_t i=0; i<backwardHistoNo.size(); i++) {
if (fT0s[2*i+1] == -1.0) { // i.e. not present in the msr-file and data file, use the estimated T0
fT0s[2*i+1] = runData->GetT0BinEstimated(backwardHistoNo[i]);
fRunInfo->SetT0Bin(fT0s[2*i+1], 2*i+1);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
std::cerr << std::endl << ">> run: " << fRunInfo->GetRunName()->Data();
std::cerr << std::endl << ">> will try the estimated one: backward t0 = " << runData->GetT0BinEstimated(backwardHistoNo[i]);
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
std::cerr << std::endl;
}
}
// check if t0 is within proper bounds
for (UInt_t i=0; i<forwardHistoNo.size(); i++) {
if ((fT0s[2*i] < 0) || (fT0s[2*i] > static_cast<Int_t>(runData->GetDataBin(forwardHistoNo[i])->size()))) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperT0(): **ERROR** t0 data bin (" << fT0s[2*i] << ") doesn't make any sense!";
std::cerr << std::endl << ">> forwardHistoNo " << forwardHistoNo[i];
std::cerr << std::endl;
return false;
}
}
for (UInt_t i=0; i<backwardHistoNo.size(); i++) {
if ((fT0s[2*i+1] < 0) || (fT0s[2*i+1] > static_cast<Int_t>(runData->GetDataBin(backwardHistoNo[i])->size()))) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::PrepareData(): **ERROR** t0 data bin (" << fT0s[2*i+1] << ") doesn't make any sense!";
std::cerr << std::endl << ">> backwardHistoNo " << backwardHistoNo[i];
std::cerr << std::endl;
return false;
}
}
// check if addrun's are present, and if yes add the necessary t0's
if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present
PRawRunData *addRunData;
fAddT0s.resize(fRunInfo->GetRunNameSize()-1); // resize to the number of addruns
for (UInt_t i=1; i<fRunInfo->GetRunNameSize(); i++) {
// get run to be added to the main one
addRunData = fRawData->GetRunData(*(fRunInfo->GetRunName(i)));
if (addRunData == 0) { // couldn't get run
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperT0(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!";
std::cerr << std::endl;
return false;
}
// feed all T0's
// first init T0's, T0's are stored as (forward T0, backward T0, etc.)
fAddT0s[i-1].clear();
fAddT0s[i-1].resize(2*forwardHistoNo.size());
for (UInt_t j=0; j<fAddT0s[i-1].size(); j++) {
fAddT0s[i-1][j] = -1.0;
}
// fill in the T0's from the msr-file (if present)
for (Int_t j=0; j<fRunInfo->GetAddT0BinSize(i-1); j++) {
fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i-1, j);
}
// fill in the T0's from the data file, if not already present in the msr-file
for (UInt_t j=0; j<forwardHistoNo.size(); j++) {
if (fAddT0s[i-1][2*j] == -1.0) // i.e. not present in the msr-file, try the data file
if (addRunData->GetT0Bin(forwardHistoNo[j]) > 0.0) {
fAddT0s[i-1][2*j] = addRunData->GetT0Bin(forwardHistoNo[j]);
fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j], i-1, 2*j);
}
}
for (UInt_t j=0; j<backwardHistoNo.size(); j++) {
if (fAddT0s[i-1][2*j+1] == -1.0) // i.e. not present in the msr-file, try the data file
if (addRunData->GetT0Bin(backwardHistoNo[j]) > 0.0) {
fAddT0s[i-1][2*j+1] = addRunData->GetT0Bin(backwardHistoNo[j]);
fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j+1], i-1, 2*j+1);
}
}
// fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file
for (UInt_t j=0; j<forwardHistoNo.size(); j++) {
if (fAddT0s[i-1][2*j] == -1.0) { // i.e. not present in the msr-file and data file, use the estimated T0
fAddT0s[i-1][2*j] = addRunData->GetT0BinEstimated(forwardHistoNo[j]);
fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j], i-1, 2*j);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
std::cerr << std::endl << ">> run: " << fRunInfo->GetRunName(i)->Data();
std::cerr << std::endl << ">> will try the estimated one: forward t0 = " << addRunData->GetT0BinEstimated(forwardHistoNo[j]);
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
std::cerr << std::endl;
}
}
for (UInt_t j=0; j<backwardHistoNo.size(); j++) {
if (fAddT0s[i-1][2*j+1] == -1.0) { // i.e. not present in the msr-file and data file, use the estimated T0
fAddT0s[i-1][2*j+1] = addRunData->GetT0BinEstimated(backwardHistoNo[j]);
fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j+1], i-1, 2*j+1);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
std::cerr << std::endl << ">> run: " << fRunInfo->GetRunName(i)->Data();
std::cerr << std::endl << ">> will try the estimated one: backward t0 = " << runData->GetT0BinEstimated(backwardHistoNo[j]);
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
std::cerr << std::endl;
}
}
}
}
return true;
}
//--------------------------------------------------------------------------
// 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.
*
* \param runData raw run data needed to perform some crosschecks
* \param histoNo histogram number (within a run). histoNo[0]: forward histogram number, histNo[1]: backward histogram number
*
* <b>return:</b>
* - true if everthing went smooth
* - false, otherwise.
*/
Bool_t PRunAsymmetryBNMR::GetProperDataRange(PRawRunData* runData, UInt_t histoNo[2])
{
// first get start/end data
Int_t start[2] = {fRunInfo->GetDataRange(0), fRunInfo->GetDataRange(2)};
Int_t end[2] = {fRunInfo->GetDataRange(1), fRunInfo->GetDataRange(3)};
// check if data range has been provided in the RUN block. If not, try the GLOBAL block
if (start[0] == -1) {
start[0] = fMsrInfo->GetMsrGlobal()->GetDataRange(0);
}
if (start[1] == -1) {
start[1] = fMsrInfo->GetMsrGlobal()->GetDataRange(2);
}
if (end[0] == -1) {
end[0] = fMsrInfo->GetMsrGlobal()->GetDataRange(1);
}
if (end[1] == -1) {
end[1] = fMsrInfo->GetMsrGlobal()->GetDataRange(3);
}
Double_t t0[2] = {fT0s[0], fT0s[1]};
Int_t offset = static_cast<Int_t>((10.0e-3/fTimeResolution)); // needed in case first good bin is not given, default = 10ns
// check if data range has been provided, and if not try to estimate them
if (start[0] < 0) {
start[0] = static_cast<Int_t>(t0[0])+offset;
fRunInfo->SetDataRange(start[0], 0);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** data range (forward) was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start[0] << ".";
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS DOES MAKE ANY SENSE.";
std::cerr << std::endl;
}
if (start[1] < 0) {
start[1] = static_cast<Int_t>(t0[1])+offset;
fRunInfo->SetDataRange(start[1], 2);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** data range (backward) was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start[1] << ".";
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS DOES MAKE ANY SENSE.";
std::cerr << std::endl;
}
if (end[0] < 0) {
end[0] = runData->GetDataBin(histoNo[0])->size();
fRunInfo->SetDataRange(end[0], 1);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** data range (forward) was not provided, will try data range end = " << end[0] << ".";
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS DOES MAKE ANY SENSE.";
std::cerr << std::endl;
}
if (end[1] < 0) {
end[1] = runData->GetDataBin(histoNo[1])->size();
fRunInfo->SetDataRange(end[1], 3);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** data range (backward) was not provided, will try data range end = " << end[1] << ".";
std::cerr << std::endl << ">> NO GUARANTEE THAT THIS DOES MAKE ANY SENSE.";
std::cerr << std::endl;
}
// check if start, end, and t0 make any sense
// 1st check if start and end are in proper order
for (UInt_t i=0; i<2; i++) {
if (end[i] < start[i]) { // need to swap them
Int_t keep = end[i];
end[i] = start[i];
start[i] = keep;
}
// 2nd check if start is within proper bounds
if ((start[i] < 0) || (start[i] > static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size()))) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **ERROR** start data bin doesn't make any sense!";
std::cerr << std::endl;
return false;
}
// 3rd check if end is within proper bounds
if (end[i] < 0) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **ERROR** end data bin (" << end[i] << ") doesn't make any sense!";
std::cerr << std::endl;
return false;
}
if (end[i] > static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size())) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** end data bin (" << end[i] << ") > histo length (" << static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size()) << ").";
std::cerr << std::endl << ">> Will set end = (histo length - 1). Consider to change it in the msr-file." << std::endl;
std::cerr << std::endl;
end[i] = static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size())-1;
}
// 4th check if t0 is within proper bounds
if ((t0[i] < 0) || (t0[i] > static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size()))) {
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **ERROR** t0 data bin doesn't make any sense!";
std::cerr << std::endl;
return false;
}
}
// check that start-t0 is the same for forward as for backward, otherwise take max(start[i]-t0[i])
if (fabs(static_cast<Double_t>(start[0])-t0[0]) > fabs(static_cast<Double_t>(start[1])-t0[1])){
start[1] = static_cast<Int_t>(t0[1] + static_cast<Double_t>(start[0]) - t0[0]);
end[1] = static_cast<Int_t>(t0[1] + static_cast<Double_t>(end[0]) - t0[0]);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange **WARNING** needed to shift backward data range.";
std::cerr << std::endl << ">> given: " << fRunInfo->GetDataRange(2) << ", " << fRunInfo->GetDataRange(3);
std::cerr << std::endl << ">> used : " << start[1] << ", " << end[1];
std::cerr << std::endl;
}
if (fabs(static_cast<Double_t>(start[0])-t0[0]) < fabs(static_cast<Double_t>(start[1])-t0[1])){
start[0] = static_cast<Int_t>(t0[0] + static_cast<Double_t>(start[1]) - t0[1]);
end[0] = static_cast<Int_t>(t0[0] + static_cast<Double_t>(end[1]) - t0[1]);
std::cerr << std::endl << ">> PRunAsymmetryBNMR::GetProperDataRange **WARNING** needed to shift forward data range.";
std::cerr << std::endl << ">> given: " << fRunInfo->GetDataRange(0) << ", " << fRunInfo->GetDataRange(1);
std::cerr << std::endl << ">> used : " << start[0] << ", " << end[0];
std::cerr << std::endl;
}
// keep good bins for potential latter use
fGoodBins[0] = start[0];
fGoodBins[1] = end[0];
fGoodBins[2] = start[1];
fGoodBins[3] = end[1];
return true;
}
//--------------------------------------------------------------------------
// GetProperFitRange (private)
//--------------------------------------------------------------------------
/**
* \brief Determines the proper fit range for the run.
*
* Fit ranges can be specified in two ways:
* - Time format: fit <start> <end> 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)
{
// set fit start/end time; first check RUN Block
fFitStartTime = fRunInfo->GetFitRange(0);
fFitEndTime = fRunInfo->GetFitRange(1);
// if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now
if (fRunInfo->IsFitRangeInBin()) {
fFitStartTime = (fGoodBins[0] + fRunInfo->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt
fFitEndTime = (fGoodBins[1] - fRunInfo->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution; // (lgb-n1-t0)*dt
// write these times back into the data structure. This way it is available when writting the log-file
fRunInfo->SetFitRange(fFitStartTime, 0);
fRunInfo->SetFitRange(fFitEndTime, 1);
}
if (fFitStartTime == PMUSR_UNDEFINED) { // fit start/end NOT found in the RUN block, check GLOBAL block
fFitStartTime = globalBlock->GetFitRange(0);
fFitEndTime = globalBlock->GetFitRange(1);
// if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now
if (globalBlock->IsFitRangeInBin()) {
fFitStartTime = (fGoodBins[0] + globalBlock->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt
fFitEndTime = (fGoodBins[1] - globalBlock->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution; // (lgb-n1-t0)*dt
// write these times back into the data structure. This way it is available when writting the log-file
globalBlock->SetFitRange(fFitStartTime, 0);
globalBlock->SetFitRange(fFitEndTime, 1);
}
}
if ((fFitStartTime == PMUSR_UNDEFINED) || (fFitEndTime == PMUSR_UNDEFINED)) {
fFitStartTime = (fGoodBins[0] - fT0s[0]) * fTimeResolution; // (fgb-t0)*dt
fFitEndTime = (fGoodBins[1] - fT0s[0]) * fTimeResolution; // (lgb-t0)*dt
std::cerr << ">> PRunSingleHisto::GetProperFitRange(): **WARNING** Couldn't get fit start/end time!" << std::endl;
std::cerr << ">> Will set it to fgb/lgb which given in time is: " << fFitStartTime << "..." << fFitEndTime << " (usec)" << std::endl;
}
}
//--------------------------------------------------------------------------
// EstimateAlpha (private)
//--------------------------------------------------------------------------
/**
* \brief Estimates the α parameter from histogram data.
*
* 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()
{
Double_t SumF[2] = {0.0, 0.0};
Double_t SumB[2] = {0.0, 0.0};
Double_t alpha = 1;
// forward
for (Int_t i=0; i<fForwardp.size(); i++) {
SumF[0] += fForwardp[i];
SumF[1] += fForwardm[i];
}
// backward
for (Int_t i=0; i<fBackwardp.size(); i++) {
SumB[0] += fBackwardp[i];
SumB[1] += fBackwardm[i];
}
// Spit out estimated alpha value from total counts (Bp+Bm)/(Fp+Fm)
if ( (SumF[0]+SumF[1]) == 0) {
alpha = 1;
} else {
alpha = (SumB[0]+SumB[1])/(SumF[0]+SumF[1]);
}
std::cout << std::endl << ">> PRunAsymmetryBNMR::EstimateAlpha(): alpha estimate=" << alpha << std::endl;
return alpha;
}
Reference in New Issue View Git Blame Copy Permalink