LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

allow alpha, beta in the RUN block being expressed via a function.

Browse Source
This commit is contained in:
suter_a 2020-05-25 16:21:16 +02:00
parent 9af6b1fb8e
commit c930298972
5 changed files with 499 additions and 143 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -5,7 +5,7 @@ if (CMAKE_VERSION GREATER_EQUAL 3.12)
cmake_policy(SET CMP0075 NEW)
endif (CMAKE_VERSION GREATER_EQUAL 3.12)
project(musrfit VERSION 1.6.0 LANGUAGES C CXX)
project(musrfit VERSION 1.6.1 LANGUAGES C CXX)
#--- musrfit specific options -------------------------------------------------
option(nexus "build optional NeXus support. Needed for ISIS" OFF)

74
src/classes/PMsrHandler.cpp
View File

@ -819,15 +819,24 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
} else if (sstr.BeginsWith("alpha ")) {
fout.width(16);
fout << std::left << "alpha";
// check of alpha is given as a function
if (fRuns[runNo].GetAlphaParamNo() >= MSR_PARAM_FUN_OFFSET)
fout << "fun" << fRuns[runNo].GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
else
fout << fRuns[runNo].GetAlphaParamNo() << std::endl;
fout << std::endl;
} else if (sstr.BeginsWith("beta ")) {
fout.width(16);
fout << std::left << "beta";
if (fRuns[runNo].GetBetaParamNo() >= MSR_PARAM_FUN_OFFSET)
fout << "fun" << fRuns[runNo].GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
else
fout << fRuns[runNo].GetBetaParamNo() << std::endl;
fout << std::endl;
} else if (sstr.BeginsWith("norm")) {
fout.width(16);
fout << std::left << "norm";
// check if norm is give as a function
// check if norm is given as a function
if (fRuns[runNo].GetNormParamNo() >= MSR_PARAM_FUN_OFFSET)
fout << "fun" << fRuns[runNo].GetNormParamNo()-MSR_PARAM_FUN_OFFSET;
else
@ -1939,14 +1948,24 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, std::map<UInt_t, TString
if (fRuns[i].GetAlphaParamNo() != -1) {
fout.width(16);
fout << std::left << "alpha";
fout << fRuns[i].GetAlphaParamNo() << std::endl;
// check if alpha is give as a function
if (fRuns[i].GetAlphaParamNo() >= MSR_PARAM_FUN_OFFSET)
fout << "fun" << fRuns[i].GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET;
else
fout << fRuns[i].GetAlphaParamNo();
fout << std::endl;
}
// beta
if (fRuns[i].GetBetaParamNo() != -1) {
fout.width(16);
fout << std::left << "beta";
fout << fRuns[i].GetBetaParamNo() << std::endl;
// check if beta is give as a function
if (fRuns[i].GetBetaParamNo() >= MSR_PARAM_FUN_OFFSET)
fout << "fun" << fRuns[i].GetBetaParamNo()-MSR_PARAM_FUN_OFFSET;
else
fout << fRuns[i].GetBetaParamNo();
fout << std::endl;
}
// norm
@ -3379,6 +3398,12 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
param.SetAlphaParamNo(ival);
else
error = true;
} else if (str.Contains("fun")) {
Int_t no;
if (FilterNumber(str, "fun", MSR_PARAM_FUN_OFFSET, no))
param.SetAlphaParamNo(no);
else
error = true;
} else {
error = true;
}
@ -3401,6 +3426,12 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
param.SetBetaParamNo(ival);
else
error = true;
} else if (str.Contains("fun")) {
Int_t no;
if (FilterNumber(str, "fun", MSR_PARAM_FUN_OFFSET, no))
param.SetBetaParamNo(no);
else
error = true;
} else {
error = true;
}
@ -4582,9 +4613,12 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
// handle a single PLOT block
while ((iter1 != iter2) && !error) {
TString line = iter1->fLine;
if (line.First('#') != -1) // remove trailing comment before proceed
line.Resize(line.First('#'));
if (iter1->fLine.Contains("PLOT")) { // handle plot header
tokens = iter1->fLine.Tokenize(" \t");
if (line.Contains("PLOT")) { // handle plot header
tokens = line.Tokenize(" \t");
if (!tokens) {
std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
std::cerr << std::endl << std::endl;
@ -4605,9 +4639,9 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
delete tokens;
tokens = nullptr;
}
} else if (iter1->fLine.Contains("lifetimecorrection", TString::kIgnoreCase)) {
} else if (line.Contains("lifetimecorrection", TString::kIgnoreCase)) {
param.fLifeTimeCorrection = true;
} else if (iter1->fLine.Contains("runs", TString::kIgnoreCase)) { // handle plot runs
} else if (line.Contains("runs", TString::kIgnoreCase)) { // handle plot runs
TComplex run;
PStringNumberList *rl;
std::string errorMsg;
@ -4623,7 +4657,7 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
case MSR_PLOT_ASYM_RRF:
case MSR_PLOT_NON_MUSR:
case MSR_PLOT_MU_MINUS:
rl = new PStringNumberList(iter1->fLine.Data());
rl = new PStringNumberList(line.Data());
if (!rl->Parse(errorMsg, true)) {
std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
std::cerr << std::endl << ">> Error Message: " << errorMsg;
@ -4643,14 +4677,14 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
error = true;
break;
}
} else if (iter1->fLine.Contains("range ", TString::kIgnoreCase)) { // handle plot range
} else if (line.Contains("range ", TString::kIgnoreCase)) { // handle plot range
// remove previous entries
param.fTmin.clear();
param.fTmax.clear();
param.fYmin.clear();
param.fYmax.clear();
tokens = iter1->fLine.Tokenize(" \t");
tokens = line.Tokenize(" \t");
if (!tokens) {
std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
std::cerr << std::endl << std::endl;
@ -4700,14 +4734,14 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
delete tokens;
tokens = nullptr;
}
} else if (iter1->fLine.Contains("sub_ranges", TString::kIgnoreCase)) {
} else if (line.Contains("sub_ranges", TString::kIgnoreCase)) {
// remove previous entries
param.fTmin.clear();
param.fTmax.clear();
param.fYmin.clear();
param.fYmax.clear();
tokens = iter1->fLine.Tokenize(" \t");
tokens = line.Tokenize(" \t");
if (!tokens) {
std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
std::cerr << std::endl << std::endl;
@ -4762,11 +4796,11 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
delete tokens;
tokens = nullptr;
}
} else if (iter1->fLine.Contains("use_fit_ranges", TString::kIgnoreCase)) {
} else if (line.Contains("use_fit_ranges", TString::kIgnoreCase)) {
param.fUseFitRanges = true;
// check if y-ranges are given
tokens = iter1->fLine.Tokenize(" \t");
tokens = line.Tokenize(" \t");
if (!tokens) {
std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo;
std::cerr << std::endl << std::endl;
@ -5207,12 +5241,20 @@ UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx)
paramVector.push_back(fRuns[idx].GetBkgFitParamNo());
// get alpha parameter if present (asymmetry fit)
if (fRuns[idx].GetAlphaParamNo() != -1)
if (fRuns[idx].GetAlphaParamNo() != -1) {
if (fRuns[idx].GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) // parameter
paramVector.push_back(fRuns[idx].GetAlphaParamNo());
else // function
funVector.push_back(fRuns[idx].GetAlphaParamNo() - MSR_PARAM_FUN_OFFSET);
}
// get beta parameter if present (asymmetry fit)
if (fRuns[idx].GetBetaParamNo() != -1)
if (fRuns[idx].GetBetaParamNo() != -1) {
if (fRuns[idx].GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) // parameter
paramVector.push_back(fRuns[idx].GetBetaParamNo());
else // function
funVector.push_back(fRuns[idx].GetBetaParamNo() - MSR_PARAM_FUN_OFFSET);
}
// go through the theory block and collect parameters
// possible entries: number -> parameter, fun<number> -> function, map<number> -> maps

174
src/classes/PRunAsymmetry.cpp
View File

@ -111,7 +111,16 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UIn
return;
}
// check if alpha parameter is within proper bounds
if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > static_cast<Int_t>(param->size()))) {
if (fRunInfo->GetAlphaParamNo() >= MSR_PARAM_FUN_OFFSET) { // alpha seems to be a function
if ((fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET < 0) ||
(fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET > msrInfo->GetNoOfFuncs())) {
std::cerr << std::endl << ">> PRunAsymmetry::PRunAsymmetry(): **ERROR** alpha parameter is a function with no = " << fRunInfo->GetAlphaParamNo();
std::cerr << std::endl << ">> This is out of bound, since there are only " << msrInfo->GetNoOfFuncs() << " functions.";
std::cerr << std::endl;
fValid = false;
return;
}
} else if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > static_cast<Int_t>(param->size()))) {
std::cerr << std::endl << ">> PRunAsymmetry::PRunAsymmetry(): **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;
@ -120,10 +129,11 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UIn
}
// check if alpha is fixed
Bool_t alphaFixedToOne = false;
if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is parameter
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
@ -196,6 +206,58 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
Double_t time(1.0);
Int_t i;
// determine alpha/beta
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);
}
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);
}
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);
}
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);
}
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.
@ -206,34 +268,12 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
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,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
#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();
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(time, par, fFuncValues);
break;
case 2: // alpha != 1, beta == 1
a = par[fRunInfo->GetAlphaParamNo()-1];
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
break;
case 3: // alpha == 1, beta != 1
b = par[fRunInfo->GetBetaParamNo()-1];
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
break;
case 4: // alpha != 1, beta != 1
a = par[fRunInfo->GetAlphaParamNo()-1];
b = par[fRunInfo->GetBetaParamNo()-1];
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
break;
default:
asymFcnValue = 0.0;
break;
}
diff = fData.GetValue()->at(i) - asymFcnValue;
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
}
@ -429,18 +469,46 @@ void PRunAsymmetry::CalcTheory()
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);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (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);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = 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);
}
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);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
break;
@ -1138,16 +1206,44 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
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);
}
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);
}
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);
}
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);
}
break;
default:
break;
@ -1349,16 +1445,44 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]
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);
}
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);
}
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);
}
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);
}
break;
default:
break;

174
src/classes/PRunAsymmetryBNMR.cpp
View File

@ -206,56 +206,90 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector<Double_t>& par)
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);
}
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);
}
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);
}
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);
}
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);
}
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);
alphaest = fRunInfo->GetEstimatedAlpha();
#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,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
#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();
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(time, par, fFuncValues);
break;
case 2: // alpha != 1, beta == 1
a = par[fRunInfo->GetAlphaParamNo()-1];
f = fTheory->Func(time, par, fFuncValues)/2.0;
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
b = par[fRunInfo->GetBetaParamNo()-1];
f = fTheory->Func(time, par, fFuncValues)/2.0;
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
a = par[fRunInfo->GetAlphaParamNo()-1];
b = par[fRunInfo->GetBetaParamNo()-1];
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));
break;
case 5: // alpha ?? , beta == 1
a = alphaest;
f = fTheory->Func(time, par, fFuncValues)/2.0;
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;
b = par[fRunInfo->GetBetaParamNo()-1];
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));
break;
default:
asymFcnValue = 0.0;
break;
}
diff = fData.GetValue()->at(i) - asymFcnValue;
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
}
@ -455,18 +489,46 @@ void PRunAsymmetryBNMR::CalcTheory()
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);
}
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);
}
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);
}
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);
}
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;
@ -477,7 +539,14 @@ void PRunAsymmetryBNMR::CalcTheory()
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);
}
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;
@ -1287,16 +1356,44 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
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);
}
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);
}
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);
}
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);
}
break;
case 5: // alpha ?? , beta == 1
// use estimated value
@ -1304,7 +1401,14 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2
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);
}
break;
default:
break;

132
src/classes/PRunAsymmetryRRF.cpp
View File

@ -188,6 +188,58 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
Double_t time(1.0);
Int_t i;
// determine alpha/beta
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);
}
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);
}
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);
}
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);
}
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.
@ -198,34 +250,12 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
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,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
#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();
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(time, par, fFuncValues);
break;
case 2: // alpha != 1, beta == 1
a = par[fRunInfo->GetAlphaParamNo()-1];
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0));
break;
case 3: // alpha == 1, beta != 1
b = par[fRunInfo->GetBetaParamNo()-1];
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0));
break;
case 4: // alpha != 1, beta != 1
a = par[fRunInfo->GetAlphaParamNo()-1];
b = par[fRunInfo->GetBetaParamNo()-1];
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
break;
default:
asymFcnValue = 0.0;
break;
}
diff = fData.GetValue()->at(i) - asymFcnValue;
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
}
@ -421,18 +451,46 @@ void PRunAsymmetryRRF::CalcTheory()
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);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (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);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = 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);
}
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);
}
f = fTheory->Func(time, par, fFuncValues);
asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0));
break;
@ -1034,16 +1092,44 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]
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);
}
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);
}
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);
}
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);
}
break;
default:
break;