#if __cplusplus >= 201103L
#define ROOT_CPLUSPLUS11 1
#endif
#include "TROOT.h"
#include "TClass.h"
#include "TMethod.h"
#include "TMath.h"
#include "TF1.h"
#include "TMethodCall.h"
#include <TBenchmark.h>
#include "TError.h"
#include "TInterpreter.h"
#include "TFormula.h"
#include <cassert>
#include <iostream>
#include <unordered_map>
using namespace std;
#ifdef WIN32
#pragma optimize("",off)
#endif
#include "v5/TFormula.h"
ClassImp(TFormula)
static const TString gNamePrefix = "TFormula__";
static std::unordered_map<std::string, void *> gClingFunctions = std::unordered_map<std::string, void * >();
Bool_t TFormula::IsOperator(const char c)
{
char ops[] = { '+','^', '-','/','*','<','>','|','&','!','=','?'};
Int_t opsLen = sizeof(ops)/sizeof(char);
for(Int_t i = 0; i < opsLen; ++i)
if(ops[i] == c)
return true;
return false;
}
Bool_t TFormula::IsBracket(const char c)
{
char brackets[] = { ')','(','{','}'};
Int_t bracketsLen = sizeof(brackets)/sizeof(char);
for(Int_t i = 0; i < bracketsLen; ++i)
if(brackets[i] == c)
return true;
return false;
}
Bool_t TFormula::IsFunctionNameChar(const char c)
{
return !IsBracket(c) && !IsOperator(c) && c != ',' && c != ' ';
}
Bool_t TFormula::IsDefaultVariableName(const TString &name)
{
return name == "x" || name == "z" || name == "y" || name == "t";
}
Bool_t TFormula::IsScientificNotation(const TString & formula, int i)
{
if ( (formula[i] == 'e' || formula[i] == 'E') && (i > 0 && i < formula.Length()-1) ) {
if ( (isdigit(formula[i-1]) || formula[i-1] == '.') && ( isdigit(formula[i+1]) || formula[i+1] == '+' || formula[i+1] == '-' ) )
return true;
}
return false;
}
Bool_t TFormula::IsHexadecimal(const TString & formula, int i)
{
if ( (formula[i] == 'x' || formula[i] == 'X') && (i > 0 && i < formula.Length()-1) && formula[i-1] == '0') {
if (isdigit(formula[i+1]) )
return true;
static char hex_values[12] = { 'a','A', 'b','B','c','C','d','D','e','E','f','F'};
for (int jjj = 0; jjj < 12; ++jjj) {
if (formula[i+1] == hex_values[jjj])
return true;
}
}
return false;
}
bool TFormulaParamOrder::operator() (const TString& a, const TString& b) const {
if ( a[0] == 'p' && a.Length() > 1) {
if ( b[0] == 'p' && b.Length() > 1) {
TString lhs = a(1,a.Length()-1);
TString rhs = b(1,b.Length()-1);
if (lhs.IsDigit() && rhs.IsDigit() )
return (lhs.Atoi() < rhs.Atoi() );
}
else {
return true;
}
}
else {
if ( b[0] == 'p' && b.Length() > 1)
return false;
if (a.IsDigit() && b.IsDigit() )
return (a.Atoi() < b.Atoi() );
}
return a < b;
}
TFormula::TFormula()
{
fName = "";
fTitle = "";
fClingInput = "";
fReadyToExecute = false;
fClingInitialized = false;
fAllParametersSetted = false;
fMethod = 0;
fNdim = 0;
fNpar = 0;
fNumber = 0;
fClingName = "";
fFormula = "";
}
static bool IsReservedName(const char* name){
if (strlen(name)!=1) return false;
for (auto const & specialName : {"x","y","z","t"}){
if (strcmp(name,specialName)==0) return true;
}
return false;
}
TFormula::~TFormula()
{
if (!TestBit(TFormula::kNotGlobal) && gROOT ) {
R__LOCKGUARD2(gROOTMutex);
gROOT->GetListOfFunctions()->Remove(this);
}
if(fMethod)
{
fMethod->Delete();
}
int nLinParts = fLinearParts.size();
if (nLinParts > 0) {
for (int i = 0; i < nLinParts; ++i) delete fLinearParts[i];
}
}
#ifdef OLD_VERSION
TFormula::TFormula(const char *name, Int_t nparams, Int_t ndims)
{
fName = name;
fTitle = "";
fClingInput = "";
fReadyToExecute = false;
fClingInitialized = false;
fAllParametersSetted = false;
fMethod = 0;
fNdim = ndims;
fNpar = 0;
fNumber = 0;
fClingName = "";
fFormula = "";
for(Int_t i = 0; i < nparams; ++i)
{
TString parName = TString::Format("%d",i);
DoAddParameter(parName,0,false);
}
}
#endif
TFormula::TFormula(const char *name, const char *formula, bool addToGlobList) :
TNamed(name,formula),
fClingInput(formula),fFormula(formula)
{
fReadyToExecute = false;
fClingInitialized = false;
fMethod = 0;
fNdim = 0;
fNpar = 0;
fNumber = 0;
FillDefaults();
if (addToGlobList && gROOT) {
TFormula *old = 0;
R__LOCKGUARD2(gROOTMutex);
old = dynamic_cast<TFormula*> ( gROOT->GetListOfFunctions()->FindObject(name) );
if (old)
gROOT->GetListOfFunctions()->Remove(old);
if (IsReservedName(name))
Error("TFormula","The name %s is reserved as a TFormula variable name.\n",name);
else
gROOT->GetListOfFunctions()->Add(this);
}
SetBit(kNotGlobal,!addToGlobList);
if (!fFormula.IsNull() ) {
PreProcessFormula(fFormula);
PrepareFormula(fFormula);
}
}
TFormula::TFormula(const TFormula &formula) : TNamed(formula.GetName(),formula.GetTitle())
{
fReadyToExecute = false;
fClingInitialized = false;
fMethod = 0;
fNdim = formula.GetNdim();
fNpar = formula.GetNpar();
fNumber = formula.GetNumber();
fFormula = formula.GetExpFormula();
FillDefaults();
if (!TestBit(TFormula::kNotGlobal) && gROOT ) {
R__LOCKGUARD2(gROOTMutex);
TFormula *old = (TFormula*)gROOT->GetListOfFunctions()->FindObject(formula.GetName());
if (old)
gROOT->GetListOfFunctions()->Remove(old);
if (IsReservedName(formula.GetName())) {
Error("TFormula","The name %s is reserved as a TFormula variable name.\n",formula.GetName());
} else
gROOT->GetListOfFunctions()->Add(this);
}
PreProcessFormula(fFormula);
PrepareFormula(fFormula);
}
TFormula& TFormula::operator=(const TFormula &rhs)
{
if (this != &rhs) {
rhs.Copy(*this);
}
return *this;
}
Int_t TFormula::Compile(const char *expression)
{
TString formula = expression;
if (formula.IsNull() ) {
formula = fFormula;
if (formula.IsNull() ) formula = GetTitle();
}
if (formula.IsNull() ) return -1;
if (IsValid() && formula == fFormula ) return 0;
if (!fFormula.IsNull() ) Clear();
fFormula = formula;
if (fVars.empty() ) FillDefaults();
PreProcessFormula(fFormula);
bool ret = PrepareFormula(fFormula);
return (ret) ? 0 : 1;
}
void TFormula::Copy(TObject &obj) const
{
TNamed::Copy(obj);
TFormula & fnew = dynamic_cast<TFormula&>(obj);
fnew.fClingParameters = fClingParameters;
fnew.fClingVariables = fClingVariables;
fnew.fFuncs = fFuncs;
fnew.fVars = fVars;
fnew.fParams = fParams;
fnew.fConsts = fConsts;
fnew.fFunctionsShortcuts = fFunctionsShortcuts;
fnew.fFormula = fFormula;
fnew.fNdim = fNdim;
fnew.fNpar = fNpar;
fnew.fNumber = fNumber;
fnew.SetParameters(GetParameters());
int nLinParts = fnew.fLinearParts.size();
if (nLinParts > 0) {
for (int i = 0; i < nLinParts; ++i) delete fnew.fLinearParts[i];
fnew.fLinearParts.clear();
}
nLinParts = fLinearParts.size();
if (nLinParts > 0) {
fnew.fLinearParts.reserve(nLinParts);
for (int i = 0; i < nLinParts; ++i) {
TFormula * linearNew = new TFormula();
TFormula * linearOld = (TFormula*) fLinearParts[i];
if (linearOld) {
linearOld->Copy(*linearNew);
fnew.fLinearParts.push_back(linearNew);
}
else
Warning("Copy","Function %s - expr %s has a dummy linear part %d",GetName(),GetExpFormula().Data(),i);
}
}
fnew.fClingInput = fClingInput;
fnew.fReadyToExecute = fReadyToExecute;
fnew.fClingInitialized = fClingInitialized;
fnew.fAllParametersSetted = fAllParametersSetted;
fnew.fClingName = fClingName;
if (fMethod) {
if (fnew.fMethod) delete fnew.fMethod;
TMethodCall *m = new TMethodCall(*fMethod);
fnew.fMethod = m;
}
fnew.fFuncPtr = fFuncPtr;
}
void TFormula::Clear(Option_t * )
{
fNdim = 0;
fNpar = 0;
fNumber = 0;
fFormula = "";
fClingName = "";
if(fMethod) fMethod->Delete();
fMethod = nullptr;
fClingVariables.clear();
fClingParameters.clear();
fReadyToExecute = false;
fClingInitialized = false;
fAllParametersSetted = false;
fFuncs.clear();
fVars.clear();
fParams.clear();
fConsts.clear();
fFunctionsShortcuts.clear();
int nLinParts = fLinearParts.size();
if (nLinParts > 0) {
for (int i = 0; i < nLinParts; ++i) delete fLinearParts[i];
}
fLinearParts.clear();
}
bool TFormula::PrepareEvalMethod()
{
if(!fMethod)
{
fMethod = new TMethodCall();
Bool_t hasParameters = (fNpar > 0);
Bool_t hasVariables = (fNdim > 0);
TString prototypeArguments = "";
if(hasVariables)
{
prototypeArguments.Append("Double_t*");
}
if(hasVariables && hasParameters)
{
prototypeArguments.Append(",");
}
if(hasParameters)
{
prototypeArguments.Append("Double_t*");
}
fMethod->InitWithPrototype(fClingName,prototypeArguments);
if(!fMethod->IsValid())
{
Error("Eval","Can't find %s function prototype with arguments %s",fClingName.Data(),prototypeArguments.Data());
return false;
}
CallFunc_t * callfunc = fMethod->GetCallFunc();
TInterpreter::CallFuncIFacePtr_t faceptr = gCling->CallFunc_IFacePtr(callfunc);
fFuncPtr = faceptr.fGeneric;
}
return true;
}
void TFormula::InputFormulaIntoCling()
{
if(!fClingInitialized && fReadyToExecute && fClingInput.Length() > 0)
{
gCling->Declare(fClingInput);
fClingInitialized = PrepareEvalMethod();
}
}
void TFormula::FillDefaults()
{
const TString defvars[] = { "x","y","z","t"};
const pair<TString,Double_t> defconsts[] = { {"pi",TMath::Pi()}, {"sqrt2",TMath::Sqrt2()},
{"infinity",TMath::Infinity()}, {"e",TMath::E()}, {"ln10",TMath::Ln10()},
{"loge",TMath::LogE()}, {"c",TMath::C()}, {"g",TMath::G()},
{"h",TMath::H()}, {"k",TMath::K()},{"sigma",TMath::Sigma()},
{"r",TMath::R()}, {"eg",TMath::EulerGamma()},{"true",1},{"false",0} };
const pair<TString,TString> funShortcuts[] =
{ {"sin","TMath::Sin" },
{"cos","TMath::Cos" }, {"exp","TMath::Exp"}, {"log","TMath::Log"}, {"log10","TMath::Log10"},
{"tan","TMath::Tan"}, {"sinh","TMath::SinH"}, {"cosh","TMath::CosH"},
{"tanh","TMath::TanH"}, {"asin","TMath::ASin"}, {"acos","TMath::ACos"},
{"atan","TMath::ATan"}, {"atan2","TMath::ATan2"}, {"sqrt","TMath::Sqrt"},
{"ceil","TMath::Ceil"}, {"floor","TMath::Floor"}, {"pow","TMath::Power"},
{"binomial","TMath::Binomial"},{"abs","TMath::Abs"},
{"min","TMath::Min"},{"max","TMath::Max"},{"sign","TMath::Sign" },
{"sq","TMath::Sq"}
};
std::vector<TString> defvars2(10);
for (int i = 0; i < 9; ++i)
defvars2[i] = TString::Format("x[%d]",i);
for(auto var : defvars)
{
int pos = fVars.size();
fVars[var] = TFormulaVariable(var,0,pos);
fClingVariables.push_back(0);
}
for(auto con : defconsts)
{
fConsts[con.first] = con.second;
}
for(auto fun : funShortcuts)
{
fFunctionsShortcuts[fun.first] = fun.second;
}
}
void TFormula::HandlePolN(TString &formula)
{
Int_t polPos = formula.Index("pol");
while(polPos != kNPOS)
{
SetBit(kLinear,1);
Bool_t defaultVariable = false;
TString variable;
Int_t openingBracketPos = formula.Index('(',polPos);
Bool_t defaultCounter = openingBracketPos == kNPOS;
Bool_t defaultDegree = true;
Int_t degree,counter;
TString sdegree;
if(!defaultCounter)
{
sdegree = formula(polPos + 3,openingBracketPos - polPos - 3);
if (!sdegree.IsDigit() ) defaultCounter = true;
}
if (!defaultCounter) {
degree = sdegree.Atoi();
counter = TString(formula(openingBracketPos+1,formula.Index(')',polPos) - openingBracketPos)).Atoi();
}
else
{
Int_t temp = polPos+3;
while(temp < formula.Length() && isdigit(formula[temp]))
{
defaultDegree = false;
temp++;
}
degree = TString(formula(polPos+3,temp - polPos - 3)).Atoi();
counter = 0;
}
fNumber = 300 + degree;
TString replacement = TString::Format("[%d]",counter);
if(polPos - 1 < 0 || !IsFunctionNameChar(formula[polPos-1]) || formula[polPos-1] == ':' )
{
variable = "x";
defaultVariable = true;
}
else
{
Int_t tmp = polPos - 1;
while(tmp >= 0 && IsFunctionNameChar(formula[tmp]) && formula[tmp] != ':')
{
tmp--;
}
variable = formula(tmp + 1, polPos - (tmp+1));
}
Int_t param = counter + 1;
Int_t tmp = 1;
while(tmp <= degree)
{
replacement.Append(TString::Format("+[%d]*%s^%d",param,variable.Data(),tmp));
param++;
tmp++;
}
TString pattern;
if(defaultCounter && !defaultDegree)
{
pattern = TString::Format("%spol%d",(defaultVariable ? "" : variable.Data()),degree);
}
else if(defaultCounter && defaultDegree)
{
pattern = TString::Format("%spol",(defaultVariable ? "" : variable.Data()));
}
else
{
pattern = TString::Format("%spol%d(%d)",(defaultVariable ? "" : variable.Data()),degree,counter);
}
if (!formula.Contains(pattern)) {
Error("HandlePolN","Error handling polynomial function - expression is %s - trying to replace %s with %s ", formula.Data(), pattern.Data(), replacement.Data() );
break;
}
formula.ReplaceAll(pattern,replacement);
polPos = formula.Index("pol");
}
}
void TFormula::HandleParametrizedFunctions(TString &formula)
{
map< pair<TString,Int_t> ,pair<TString,TString> > functions;
functions.insert(make_pair(make_pair("gaus",1),make_pair("[0]*exp(-0.5*(({V0}-[1])/[2])*(({V0}-[1])/[2]))","[0]*exp(-0.5*(({V0}-[1])/[2])*(({V0}-[1])/[2]))/(sqrt(2*pi)*[2])")));
functions.insert(make_pair(make_pair("landau",1),make_pair("[0]*TMath::Landau({V0},[1],[2],false)","[0]*TMath::Landau({V0},[1],[2],true)")));
functions.insert(make_pair(make_pair("expo",1),make_pair("exp([0]+[1]*{V0})","")));
functions.insert(make_pair(make_pair("crystalball",1),make_pair("[0]*ROOT::Math::crystalball_function({V0},[3],[4],[2],[1])","[0]*ROOT::Math::crystalball_pdf({V0},[3],[4],[2],[1])")));
functions.insert(make_pair(make_pair("breitwigner",1),make_pair("[0]*ROOT::Math::breitwigner_pdf({V0},[2],[1])","[0]*ROOT::Math::breitwigner_pdf({V0},[2],[4],[1])")));
functions.insert(make_pair(make_pair("cheb0" ,1),make_pair("ROOT::Math::Chebyshev0({V0},[0])","")));
functions.insert(make_pair(make_pair("cheb1" ,1),make_pair("ROOT::Math::Chebyshev1({V0},[0],[1])","")));
functions.insert(make_pair(make_pair("cheb2" ,1),make_pair("ROOT::Math::Chebyshev2({V0},[0],[1],[2])","")));
functions.insert(make_pair(make_pair("cheb3" ,1),make_pair("ROOT::Math::Chebyshev3({V0},[0],[1],[2],[3])","")));
functions.insert(make_pair(make_pair("cheb4" ,1),make_pair("ROOT::Math::Chebyshev4({V0},[0],[1],[2],[3],[4])","")));
functions.insert(make_pair(make_pair("cheb5" ,1),make_pair("ROOT::Math::Chebyshev5({V0},[0],[1],[2],[3],[4],[5])","")));
functions.insert(make_pair(make_pair("cheb6" ,1),make_pair("ROOT::Math::Chebyshev6({V0},[0],[1],[2],[3],[4],[5],[6])","")));
functions.insert(make_pair(make_pair("cheb7" ,1),make_pair("ROOT::Math::Chebyshev7({V0},[0],[1],[2],[3],[4],[5],[6],[7])","")));
functions.insert(make_pair(make_pair("cheb8" ,1),make_pair("ROOT::Math::Chebyshev8({V0},[0],[1],[2],[3],[4],[5],[6],[7],[8])","")));
functions.insert(make_pair(make_pair("cheb9" ,1),make_pair("ROOT::Math::Chebyshev9({V0},[0],[1],[2],[3],[4],[5],[6],[7],[8],[9])","")));
functions.insert(make_pair(make_pair("cheb10",1),make_pair("ROOT::Math::Chebyshev10({V0},[0],[1],[2],[3],[4],[5],[6],[7],[8],[9],[10])","")));
functions.insert(make_pair(make_pair("gaus",2),make_pair("[0]*exp(-0.5*(({V0}-[1])/[2])^2 - 0.5*(({V1}-[3])/[4])^2)","")));
functions.insert(make_pair(make_pair("landau",2),make_pair("[0]*TMath::Landau({V0},[1],[2],false)*TMath::Landau({V1},[3],[4],false)","")));
functions.insert(make_pair(make_pair("expo",2),make_pair("exp([0]+[1]*{V0})","exp([0]+[1]*{V0}+[2]*{V1})")));
map<TString,Int_t> functionsNumbers;
functionsNumbers["gaus"] = 100;
functionsNumbers["landau"] = 400;
functionsNumbers["expo"] = 200;
functionsNumbers["crystalball"] = 500;
formula.ReplaceAll("xygaus","gaus[x,y]");
formula.ReplaceAll("xylandau","landau[x,y]");
formula.ReplaceAll("xyexpo","expo[x,y]");
for(map<pair<TString,Int_t>,pair<TString,TString> >::iterator it = functions.begin(); it != functions.end(); ++it)
{
TString funName = it->first.first;
Int_t funPos = formula.Index(funName);
while(funPos != kNPOS)
{
Int_t lastFunPos = funPos + funName.Length();
Int_t iposBefore = funPos - 1;
if (iposBefore >= 0) {
assert( iposBefore < formula.Length() );
if (isalpha(formula[iposBefore] ) ) {
break;
}
}
Bool_t isNormalized = false;
if (lastFunPos < formula.Length() ) {
isNormalized = (formula[lastFunPos] == 'n');
if (isNormalized) lastFunPos += 1;
if (lastFunPos < formula.Length() ) {
if (isalnum(formula[lastFunPos] ) ) break;
if (formula[lastFunPos] != '[' && formula[lastFunPos] != '(' && ! IsOperator(formula[lastFunPos] ) ) {
funPos = formula.Index(funName,lastFunPos);
continue;
}
}
}
if(isNormalized)
{
SetBit(kNormalized,1);
}
std::vector<TString> variables;
Int_t dim = 0;
TString varList = "";
Bool_t defaultVariable = false;
Int_t openingBracketPos = funPos + funName.Length() + (isNormalized ? 1 : 0);
Int_t closingBracketPos = kNPOS;
if(openingBracketPos > formula.Length() || formula[openingBracketPos] != '[')
{
dim = 1;
variables.resize(dim);
variables[0] = "x";
defaultVariable = true;
}
else
{
closingBracketPos = formula.Index(']',openingBracketPos);
varList = formula(openingBracketPos+1,closingBracketPos - openingBracketPos - 1);
dim = varList.CountChar(',') + 1;
variables.resize(dim);
Int_t Nvar = 0;
TString varName = "";
for(Int_t i = 0 ; i < varList.Length(); ++i)
{
if(IsFunctionNameChar(varList[i]))
{
varName.Append(varList[i]);
}
if(varList[i] == ',')
{
variables[Nvar] = varName;
varName = "";
Nvar++;
}
}
if(varName != "")
{
variables[Nvar] = varName;
}
}
if(dim != it->first.second)
{
pair<TString,Int_t> key = make_pair(funName,dim);
if(functions.find(key) == functions.end())
{
Error("PreProcessFormula","%d dimension function %s is not defined as parametrized function.",dim,funName.Data());
return;
}
break;
}
Int_t openingParenthesisPos = (closingBracketPos == kNPOS) ? openingBracketPos : closingBracketPos + 1;
bool defaultCounter = (openingParenthesisPos > formula.Length() || formula[openingParenthesisPos] != '(');
Int_t counter;
if(defaultCounter)
{
counter = 0;
}
else
{
counter = TString(formula(openingParenthesisPos+1,formula.Index(')',funPos) - openingParenthesisPos -1)).Atoi();
}
TString body = (isNormalized ? it->second.second : it->second.first);
if(isNormalized && body == "")
{
Error("PreprocessFormula","%d dimension function %s has no normalized form.",it->first.second,funName.Data());
break;
}
for(int i = 0 ; i < body.Length() ; ++i)
{
if(body[i] == '{')
{
i += 2;
Int_t num = TString(body(i,body.Index('}',i) - i)).Atoi();
TString variable = variables[num];
TString pattern = TString::Format("{V%d}",num);
i -= 2;
body.Replace(i, pattern.Length(),variable,variable.Length());
i += variable.Length()-1;
}
else if(body[i] == '[')
{
Int_t tmp = i;
while(tmp < body.Length() && body[tmp] != ']')
{
tmp++;
}
Int_t num = TString(body(i+1,tmp - 1 - i)).Atoi();
num += counter;
TString replacement = TString::Format("%d",num);
body.Replace(i+1,tmp - 1 - i,replacement,replacement.Length());
i += replacement.Length() + 1;
}
}
TString pattern;
if(defaultCounter && defaultVariable)
{
pattern = TString::Format("%s%s",
funName.Data(),
(isNormalized ? "n" : ""));
}
if(!defaultCounter && defaultVariable)
{
pattern = TString::Format("%s%s(%d)",
funName.Data(),
(isNormalized ? "n" : ""),
counter);
}
if(defaultCounter && !defaultVariable)
{
pattern = TString::Format("%s%s[%s]",
funName.Data(),
(isNormalized ? "n":""),
varList.Data());
}
if(!defaultCounter && !defaultVariable)
{
pattern = TString::Format("%s%s[%s](%d)",
funName.Data(),
(isNormalized ? "n" : ""),
varList.Data(),
counter);
}
TString replacement = body;
if (fNumber == 0 && formula.Length() <= (pattern.Length()-funPos) +1 ) {
fNumber = functionsNumbers[funName] + 10*(dim-1);
}
formula.Replace(funPos,pattern.Length(),replacement,replacement.Length());
funPos = formula.Index(funName);
}
}
}
void TFormula::HandleExponentiation(TString &formula)
{
Int_t caretPos = formula.Last('^');
while(caretPos != kNPOS)
{
TString right,left;
Int_t temp = caretPos;
temp--;
if(formula[temp] == ')')
{
Int_t depth = 1;
temp--;
while(depth != 0 && temp > 0)
{
if(formula[temp] == ')')
depth++;
if(formula[temp] == '(')
depth--;
temp--;
}
if (depth == 0) temp++;
}
do {
temp--;
if (temp>=2 && IsScientificNotation(formula, temp-1) ) temp-=3;
}
while(temp >= 0 && !IsOperator(formula[temp]) && !IsBracket(formula[temp]) );
assert(temp+1 >= 0);
Int_t leftPos = temp+1;
left = formula(leftPos, caretPos - leftPos);
temp = caretPos;
temp++;
if (temp >= formula.Length() ) {
Error("HandleExponentiation","Invalid position of operator ^");
return;
}
if(formula[temp] == '(')
{
Int_t depth = 1;
temp++;
while(depth != 0 && temp < formula.Length())
{
if(formula[temp] == ')')
depth--;
if(formula[temp] == '(')
depth++;
temp++;
}
temp--;
}
else {
if (formula[temp] == '-' || formula[temp] == '+' ) temp++;
Int_t depth = 0;
while(temp < formula.Length() && ( (depth > 0) || !IsOperator(formula[temp]) ) )
{
temp++;
if (temp>=2 && IsScientificNotation(formula, temp) ) temp+=2;
if (formula[temp] == '(') depth++;
if (depth > 0 && formula[temp] == ')') depth--;
}
}
right = formula(caretPos + 1, (temp - 1) - caretPos );
TString pattern = TString::Format("%s^%s",left.Data(),right.Data());
TString replacement = TString::Format("pow(%s,%s)",left.Data(),right.Data());
formula.Replace(leftPos,pattern.Length(),replacement,replacement.Length());
caretPos = formula.Last('^');
}
}
void TFormula::HandleLinear(TString &formula)
{
Int_t linPos = formula.Index("@");
if (linPos == kNPOS ) return;
Int_t NofLinParts = formula.CountChar((int)'@');
assert(NofLinParts > 0);
fLinearParts.reserve(NofLinParts + 1);
Int_t Nlinear = 0;
bool first = true;
while(linPos != kNPOS)
{
SetBit(kLinear,1);
Int_t temp = 0;
TString left;
if (first) {
temp = linPos - 1;
while(temp >= 0 && formula[temp] != '@')
{
temp--;
}
left = formula(temp+1,linPos - (temp +1));
}
temp = linPos + 1;
while(temp < formula.Length() && formula[temp] != '@')
{
temp++;
}
TString right = formula(linPos+1,temp - (linPos+1));
TString pattern = (first) ? TString::Format("%s@%s",left.Data(),right.Data()) : TString::Format("@%s",right.Data());
TString replacement = (first) ? TString::Format("([%d]*(%s))+([%d]*(%s))",Nlinear,left.Data(),Nlinear+1,right.Data()) : TString::Format("+([%d]*(%s))",Nlinear,right.Data());
Nlinear += (first) ? 2 : 1;
formula.ReplaceAll(pattern,replacement);
if (first) {
TFormula *lin1 = new TFormula("__linear1",left,false);
fLinearParts.push_back(lin1);
}
TFormula *lin2 = new TFormula("__linear2",right,false);
fLinearParts.push_back(lin2);
linPos = formula.Index("@");
first = false;
}
}
void TFormula::PreProcessFormula(TString &formula)
{
formula.ReplaceAll("**","^");
formula.ReplaceAll("++","@");
formula.ReplaceAll(" ","");
HandlePolN(formula);
HandleParametrizedFunctions(formula);
HandleExponentiation(formula);
HandleLinear(formula);
formula.ReplaceAll("--","- -");
formula.ReplaceAll("++","+ +");
}
Bool_t TFormula::PrepareFormula(TString &formula)
{
fFuncs.clear();
fReadyToExecute = false;
ExtractFunctors(formula);
fFormula = formula;
fClingInput = formula;
fFormula.ReplaceAll("{","");
fFormula.ReplaceAll("}","");
fFuncs.sort();
fFuncs.unique();
ProcessFormula(fClingInput);
if (fNumber != 0) SetPredefinedParamNames();
return fReadyToExecute && fClingInitialized;
}
void TFormula::ExtractFunctors(TString &formula)
{
TString name = "";
TString body = "";
for(Int_t i = 0 ; i < formula.Length(); ++i )
{
if(formula[i] == '[')
{
Int_t tmp = i;
i++;
TString param = "";
while(formula[i] != ']' && i < formula.Length())
{
param.Append(formula[i++]);
}
i++;
if (param.IsDigit() ) param.Insert(0,'p');
param.ReplaceAll("\\s"," ");
DoAddParameter(param,0,false);
TString replacement = TString::Format("{[%s]}",param.Data());
formula.Replace(tmp,i - tmp, replacement,replacement.Length());
fFuncs.push_back(TFormulaFunction(param));
continue;
}
if (formula[i] == '\"') {
do {
i++;
} while(formula[i] != '\"');
}
if (IsScientificNotation(formula, i) )
continue;
if (IsHexadecimal(formula, i) ) {
while ( !IsOperator(formula[i]) && i < formula.Length() ) {
i++;
}
continue;
}
if(isalpha(formula[i]) && !IsOperator(formula[i]))
{
while( IsFunctionNameChar(formula[i]) && i < formula.Length())
{
if (formula[i] == ':' && ( (i+1) < formula.Length() ) ) {
if ( formula[i+1] == ':' ) {
name.Append("::");
i+=2;
continue;
}
else
break;
}
name.Append(formula[i++]);
}
if(formula[i] == '(')
{
i++;
if(formula[i] == ')')
{
fFuncs.push_back(TFormulaFunction(name,body,0));
name = body = "";
continue;
}
Int_t depth = 1;
Int_t args = 1;
while(depth != 0 && i < formula.Length())
{
switch(formula[i])
{
case '(': depth++; break;
case ')': depth--; break;
case ',': if(depth == 1) args++; break;
}
if(depth != 0)
{
body.Append(formula[i++]);
}
}
Int_t originalBodyLen = body.Length();
ExtractFunctors(body);
formula.Replace(i-originalBodyLen,originalBodyLen,body,body.Length());
i += body.Length() - originalBodyLen;
fFuncs.push_back(TFormulaFunction(name,body,args));
}
else
{
TObject *obj = gROOT->GetListOfFunctions()->FindObject(name);
TFormula * f = dynamic_cast<TFormula*> (obj);
if (!f) {
TF1 * f1 = dynamic_cast<TF1*> (obj);
if (f1) f = f1->GetFormula();
}
if (f) {
TString replacementFormula = f->GetExpFormula();
PreProcessFormula(replacementFormula);
int nparOffset = 0;
std::vector<TString> newNames;
if (fNpar > 0) {
nparOffset = fNpar;
newNames.resize(f->GetNpar() );
for (int jpar = f->GetNpar()-1; jpar >= 0; --jpar ) {
TString pj = TString(f->GetParName(jpar));
if ( pj[0] == 'p' && TString(pj(1,pj.Length())).IsDigit() ) {
TString oldName = TString::Format("[%s]",f->GetParName(jpar));
TString newName = TString::Format("[p%d]",nparOffset+jpar);
replacementFormula.ReplaceAll(oldName,newName);
newNames[jpar] = newName;
}
else
newNames[jpar] = f->GetParName(jpar);
}
}
ExtractFunctors(replacementFormula);
for (int jpar = 0; jpar < f->GetNpar(); ++jpar) {
if (nparOffset> 0) {
assert((int) newNames.size() == f->GetNpar() );
SetParameter(newNames[jpar], f->GetParameter(jpar) );
}
else
SetParameter(f->GetParName(jpar), f->GetParameter(jpar) );
}
replacementFormula.Insert(0,'(');
replacementFormula.Insert(replacementFormula.Length(),')');
formula.Replace(i-name.Length(),name.Length(), replacementFormula, replacementFormula.Length());
i += replacementFormula.Length()-name.Length();
name = "";
continue;
}
TString replacement = TString::Format("{%s}",name.Data());
formula.Replace(i-name.Length(),name.Length(),replacement,replacement.Length());
i += 2;
fFuncs.push_back(TFormulaFunction(name));
}
}
name = body = "";
}
}
void TFormula::ProcessFormula(TString &formula)
{
for(list<TFormulaFunction>::iterator funcsIt = fFuncs.begin(); funcsIt != fFuncs.end(); ++funcsIt)
{
TFormulaFunction & fun = *funcsIt;
if(fun.fFound)
continue;
if(fun.IsFuncCall())
{
map<TString,TString>::iterator it = fFunctionsShortcuts.find(fun.GetName());
if(it != fFunctionsShortcuts.end())
{
TString shortcut = it->first;
TString full = it->second;
Ssiz_t index = formula.Index(shortcut,0);
while ( index != kNPOS) {
Ssiz_t i2 = index + shortcut.Length();
if ( (index > 0) && (isalpha( formula[index-1] ) || formula[index-1] == ':' )) {
index = formula.Index(shortcut,i2);
continue;
}
if (i2 < formula.Length() && formula[i2] != '(') {
index = formula.Index(shortcut,i2);
continue;
}
formula.Replace(index, shortcut.Length(), full);
Ssiz_t inext = index + full.Length();
index = formula.Index(shortcut,inext);
fun.fFound = true;
}
}
if(fun.fName.Contains("::"))
{
std::string name(fun.fName);
size_t index = name.rfind("::");
assert(index != std::string::npos);
TString className = fun.fName(0,fun.fName(0,index).Length());
TString functionName = fun.fName(index + 2, fun.fName.Length());
Bool_t silent = true;
TClass *tclass = new TClass(className,silent);
const TList *methodList = tclass->GetListOfAllPublicMethods();
TIter next(methodList);
TMethod *p;
while ((p = (TMethod*) next()))
{
if (strcmp(p->GetName(),functionName.Data()) == 0 &&
(fun.GetNargs() <= p->GetNargs() && fun.GetNargs() >= p->GetNargs() - p->GetNargsOpt() ) )
{
fun.fFound = true;
break;
}
}
}
if(!fun.fFound)
{
TFunction * f = (TFunction*) gROOT->GetListOfGlobalFunctions(true)->FindObject(fun.fName);
if (f && fun.GetNargs() <= f->GetNargs() && fun.GetNargs() >= f->GetNargs() - f->GetNargsOpt() )
{
fun.fFound = true;
}
}
if(!fun.fFound)
{
if (gDebug)
Info("TFormula","Could not find %s function with %d argument(s)",fun.GetName(),fun.GetNargs());
fun.fFound = false;
}
}
else
{
TFormula *old = (TFormula*)gROOT->GetListOfFunctions()->FindObject(gNamePrefix + fun.fName);
if(old)
{
assert(false);
fun.fFound = true;
TString pattern = TString::Format("{%s}",fun.GetName());
TString replacement = old->GetExpFormula();
PreProcessFormula(replacement);
ExtractFunctors(replacement);
formula.ReplaceAll(pattern,replacement);
continue;
}
map<TString,TFormulaVariable>::iterator varsIt = fVars.find(fun.GetName());
if(varsIt!= fVars.end())
{
TString name = (*varsIt).second.GetName();
Double_t value = (*varsIt).second.fValue;
AddVariable(name,value);
if(!fVars[name].fFound)
{
fVars[name].fFound = true;
int varDim = (*varsIt).second.fArrayPos;
if (varDim >= fNdim) {
fNdim = varDim+1;
for ( auto &v : fVars) {
if (v.second.fArrayPos < varDim && !v.second.fFound ) {
AddVariable(v.first, v.second.fValue);
v.second.fFound = true;
}
}
}
}
TString pattern = TString::Format("{%s}",name.Data());
TString replacement = TString::Format("x[%d]",(*varsIt).second.fArrayPos);
formula.ReplaceAll(pattern,replacement);
fun.fFound = true;
continue;
}
TString funname = fun.GetName();
if (funname.Contains("x[") && funname.Contains("]") ) {
TString sdigit = funname(2,funname.Index("]") );
int digit = sdigit.Atoi();
if (digit >= fNdim) {
fNdim = digit+1;
for (int j = 0; j < fNdim; ++j) {
TString vname = TString::Format("x[%d]",j);
if (fVars.find(vname) == fVars.end() ) {
fVars[vname] = TFormulaVariable(vname,0,j);
fVars[vname].fFound = true;
AddVariable(vname,0.);
}
}
}
fun.fFound = true;
TString pattern = TString::Format("{%s}",funname.Data());
formula.ReplaceAll(pattern,funname);
continue;
}
auto paramsIt = fParams.find(fun.GetName());
if(paramsIt != fParams.end())
{
TString pattern = TString::Format("{[%s]}",fun.GetName());
if(formula.Index(pattern) != kNPOS)
{
TString replacement = TString::Format("p[%d]",(*paramsIt).second);
formula.ReplaceAll(pattern,replacement);
}
fun.fFound = true;
continue;
}
else {
}
map<TString,Double_t>::iterator constIt = fConsts.find(fun.GetName());
if(constIt != fConsts.end())
{
TString pattern = TString::Format("{%s}",fun.GetName());
TString value = TString::Format("%lf",(*constIt).second);
formula.ReplaceAll(pattern,value);
fun.fFound = true;
continue;
}
fun.fFound = false;
}
}
if(!fReadyToExecute)
{
fReadyToExecute = true;
Bool_t hasVariables = (fNdim > 0);
Bool_t hasParameters = (fNpar > 0);
if(!hasParameters)
{
fAllParametersSetted = true;
}
if (hasParameters && ! hasVariables) {
fNdim = 1;
AddVariable("x",0);
hasVariables = true;
}
Bool_t hasBoth = hasVariables && hasParameters;
Bool_t inputIntoCling = (formula.Length() > 0);
if (inputIntoCling) {
std::string inputFormula = std::string(formula);
TString argumentsPrototype =
TString::Format("%s%s%s",(hasVariables ? "Double_t *x" : ""), (hasBoth ? "," : ""),
(hasParameters ? "Double_t *p" : ""));
fClingName = gNamePrefix;
R__LOCKGUARD2(gROOTMutex);
auto funcit = gClingFunctions.find(inputFormula);
if (funcit != gClingFunctions.end() ) {
fFuncPtr = ( TInterpreter::CallFuncIFacePtr_t::Generic_t) funcit->second;
fClingInitialized = true;
inputIntoCling = false;
}
auto hasher = gClingFunctions.hash_function();
fClingName = TString::Format("%s__id%zu",gNamePrefix.Data(),(unsigned long) hasher(inputFormula) );
fClingInput = TString::Format("Double_t %s(%s){ return %s ; }", fClingName.Data(),argumentsPrototype.Data(),inputFormula.c_str());
if(inputIntoCling) {
InputFormulaIntoCling();
if (fClingInitialized) {
R__LOCKGUARD2(gROOTMutex);
gClingFunctions.insert ( std::make_pair ( inputFormula, (void*) fFuncPtr) );
}
}
else {
fAllParametersSetted = true;
fClingInitialized = true;
}
}
}
if (!fClingInitialized) {
Bool_t allFunctorsMatched = true;
for(list<TFormulaFunction>::iterator it = fFuncs.begin(); it != fFuncs.end(); it++)
{
if(!it->fFound)
{
allFunctorsMatched = false;
if (it->GetNargs() == 0)
Error("ProcessFormula","\"%s\" has not been matched in the formula expression",it->GetName() );
else
Error("ProcessFormula","Could not find %s function with %d argument(s)",it->GetName(),it->GetNargs());
}
}
if (!allFunctorsMatched) {
Error("ProcessFormula","Formula \"%s\" is invalid !", GetExpFormula().Data() );
fReadyToExecute = false;
}
}
if (fClingInitialized && fReadyToExecute) {
auto itvar = fVars.begin();
do
{
if ( ! itvar->second.fFound ) {
itvar = fVars.erase(itvar);
}
else
itvar++;
}
while( itvar != fVars.end() );
}
}
void TFormula::SetPredefinedParamNames() {
if (fNumber == 0) return;
if (fNumber == 100) {
SetParName(0,"Constant");
SetParName(1,"Mean");
SetParName(2,"Sigma");
return;
}
if (fNumber == 110) {
SetParName(0,"Constant");
SetParName(1,"MeanX");
SetParName(2,"SigmaX");
SetParName(3,"MeanY");
SetParName(4,"SigmaY");
return;
}
if (fNumber == 200) {
SetParName(0,"Constant");
SetParName(1,"Slope");
return;
}
if (fNumber == 400) {
SetParName(0,"Constant");
SetParName(1,"MPV");
SetParName(2,"Sigma");
return;
}
if (fNumber == 500) {
SetParName(0,"Constant");
SetParName(1,"Mean");
SetParName(2,"Sigma");
SetParName(3,"Alpha");
SetParName(4,"N");
return;
}
if (fNumber == 600) {
SetParName(0,"Constant");
SetParName(1,"Mean");
SetParName(2,"Gamma");
return;
}
return;
}
const TObject* TFormula::GetLinearPart(Int_t i) const
{
if (!fLinearParts.empty()) {
int n = fLinearParts.size();
if (i < 0 || i >= n ) {
Error("GetLinearPart","Formula %s has only %d linear parts - requested %d",GetName(),n,i);
return nullptr;
}
return fLinearParts[i];
}
return nullptr;
}
void TFormula::AddVariable(const TString &name, double value)
{
if(fVars.find(name) != fVars.end() )
{
TFormulaVariable & var = fVars[name];
var.fValue = value;
if(var.fArrayPos < 0)
{
var.fArrayPos = fVars.size();
}
if(var.fArrayPos >= (int)fClingVariables.size())
{
fClingVariables.resize(var.fArrayPos+1);
}
fClingVariables[var.fArrayPos] = value;
}
else
{
TFormulaVariable var(name,value,fVars.size());
fVars[name] = var;
fClingVariables.push_back(value);
if (!fFormula.IsNull() ) {
ProcessFormula(fClingInput);
}
}
}
void TFormula::AddVariables(const TString *vars, const Int_t size)
{
Bool_t anyNewVar = false;
for(Int_t i = 0 ; i < size; ++i)
{
const TString & vname = vars[i];
TFormulaVariable &var = fVars[vname];
if(var.fArrayPos < 0)
{
var.fName = vname;
var.fArrayPos = fVars.size();
anyNewVar = true;
var.fValue = 0;
if(var.fArrayPos >= (int)fClingVariables.capacity())
{
Int_t multiplier = 2;
if(fFuncs.size() > 100)
{
multiplier = TMath::Floor(TMath::Log10(fFuncs.size()) * 10);
}
fClingVariables.reserve(multiplier * fClingVariables.capacity());
}
fClingVariables.push_back(0.0);
}
}
if(anyNewVar && !fFormula.IsNull())
{
ProcessFormula(fClingInput);
}
}
void TFormula::SetName(const char* name)
{
if (IsReservedName(name)) {
Error("SetName","The name \'%s\' is reserved as a TFormula variable name.\n"
"\tThis function will not be renamed.",name);
} else {
auto listOfFunctions = gROOT->GetListOfFunctions();
TObject* thisAsFunctionInList = nullptr;
R__LOCKGUARD2(gROOTMutex);
if (listOfFunctions){
thisAsFunctionInList = listOfFunctions->FindObject(this);
if (thisAsFunctionInList) listOfFunctions->Remove(thisAsFunctionInList);
}
TNamed::SetName(name);
if (thisAsFunctionInList) listOfFunctions->Add(thisAsFunctionInList);
}
}
void TFormula::SetVariables(const pair<TString,Double_t> *vars, const Int_t size)
{
for(Int_t i = 0; i < size; ++i)
{
pair<TString,Double_t> v = vars[i];
if(fVars.find(v.first) != fVars.end())
{
fVars[v.first].fValue = v.second;
fClingVariables[fVars[v.first].fArrayPos] = v.second;
}
else
{
Error("SetVariables","Variable %s is not defined.",v.first.Data());
}
}
}
Double_t TFormula::GetVariable(const char *name) const
{
TString sname(name);
if(fVars.find(sname) == fVars.end())
{
Error("GetVariable","Variable %s is not defined.",sname.Data());
return -1;
}
return fVars.find(sname)->second.fValue;
}
Int_t TFormula::GetVarNumber(const char *name) const
{
TString sname(name);
if(fVars.find(sname) == fVars.end())
{
Error("GetVarNumber","Variable %s is not defined.",sname.Data());
return -1;
}
return fVars.find(sname)->second.fArrayPos;
}
TString TFormula::GetVarName(Int_t ivar) const
{
if (ivar < 0 || ivar >= fNdim) return "";
for ( auto & v : fVars) {
if (v.second.fArrayPos == ivar) return v.first;
}
Error("GetVarName","Variable with index %d not found !!",ivar);
return TString();
}
void TFormula::SetVariable(const TString &name, Double_t value)
{
if(fVars.find(name) == fVars.end())
{
Error("SetVariable","Variable %s is not defined.",name.Data());
return;
}
fVars[name].fValue = value;
fClingVariables[fVars[name].fArrayPos] = value;
}
void TFormula::DoAddParameter(const TString &name, Double_t value, Bool_t processFormula)
{
if(fParams.find(name) != fParams.end() )
{
int ipos = fParams[name];
if (ipos < 0)
{
ipos = fParams.size();
fParams[name] = ipos;
}
if(ipos >= (int)fClingParameters.size())
{
if(ipos >= (int)fClingParameters.capacity())
fClingParameters.reserve( TMath::Max(int(fParams.size()), ipos+1));
fClingParameters.insert(fClingParameters.end(),ipos+1-fClingParameters.size(),0.0);
}
fClingParameters[ipos] = value;
}
else
{
fNpar++;
int pos = fParams.size();
auto ret = fParams.insert(std::make_pair(name,pos));
if (ret.second) {
if (ret.first == fParams.begin() )
pos = 0;
else {
auto previous = (ret.first);
--previous;
pos = previous->second + 1;
}
if (pos < (int)fClingParameters.size() )
fClingParameters.insert(fClingParameters.begin()+pos,value);
else {
if (pos > (int)fClingParameters.size() )
Warning("inserting parameter %s at pos %d when vector size is %d \n",name.Data(),pos,(int)fClingParameters.size() );
if(pos >= (int)fClingParameters.capacity())
fClingParameters.reserve( TMath::Max(int(fParams.size()), pos+1));
fClingParameters.insert(fClingParameters.end(),pos+1-fClingParameters.size(),0.0);
fClingParameters[pos] = value;
}
for ( auto it = ret.first; it != fParams.end(); ++it ) {
it->second = pos;
pos++;
}
}
if (processFormula) {
fClingInput.ReplaceAll(name,TString::Format("[%s]",name.Data() ) );
ProcessFormula(fClingInput);
}
}
}
Int_t TFormula::GetParNumber(const char * name) const {
auto it = fParams.find(name);
if(it == fParams.end())
{
return -1;
}
return it->second;
}
Double_t TFormula::GetParameter(const char * name) const
{
int i = GetParNumber(name);
if (i == -1) {
Error("GetParameter","Parameter %s is not defined.",name);
return TMath::QuietNaN();
}
return GetParameter( GetParNumber(name) );
}
Double_t TFormula::GetParameter(Int_t param) const
{
if(param >=0 && param < (int) fClingParameters.size())
return fClingParameters[param];
Error("GetParameter","wrong index used - use GetParameter(name)");
return TMath::QuietNaN();
}
const char * TFormula::GetParName(Int_t ipar) const
{
if (ipar < 0 || ipar >= fNpar) return "";
for ( auto & p : fParams) {
if (p.second == ipar) return p.first.Data();
}
Error("GetParName","Parameter with index %d not found !!",ipar);
return TString();
}
Double_t* TFormula::GetParameters() const
{
if(!fClingParameters.empty())
return const_cast<Double_t*>(&fClingParameters[0]);
return 0;
}
void TFormula::GetParameters(Double_t *params) const
{
for(Int_t i = 0; i < fNpar; ++i)
{
if (Int_t(fClingParameters.size()) > i)
params[i] = fClingParameters[i];
else
params[i] = -1;
}
}
void TFormula::SetParameter(const char *name, Double_t value)
{
SetParameter( GetParNumber(name), value);
#ifdef OLDPARAMS
if(fParams.find(name) == fParams.end())
{
Error("SetParameter","Parameter %s is not defined.",name.Data());
return;
}
fParams[name].fValue = value;
fParams[name].fFound = true;
fClingParameters[fParams[name].fArrayPos] = value;
fAllParametersSetted = true;
for(map<TString,TFormulaVariable>::iterator it = fParams.begin(); it != fParams.end(); it++)
{
if(!it->second.fFound)
{
fAllParametersSetted = false;
break;
}
}
#endif
}
#ifdef OLDPARAMS
void TFormula::SetParameters(const pair<TString,Double_t> *params,const Int_t size)
{
for(Int_t i = 0 ; i < size ; ++i)
{
pair<TString,Double_t> p = params[i];
if(fParams.find(p.first) == fParams.end())
{
Error("SetParameters","Parameter %s is not defined",p.first.Data());
continue;
}
fParams[p.first].fValue = p.second;
fParams[p.first].fFound = true;
fClingParameters[fParams[p.first].fArrayPos] = p.second;
}
fAllParametersSetted = true;
for(map<TString,TFormulaVariable>::iterator it = fParams.begin(); it != fParams.end(); it++)
{
if(!it->second.fFound)
{
fAllParametersSetted = false;
break;
}
}
}
#endif
void TFormula::DoSetParameters(const Double_t *params, Int_t size)
{
if(!params || size < 0 || size > fNpar) return;
if (size != (int) fClingParameters.size() ) {
Warning("SetParameters","size is not same of cling parameter size %d - %d",size,int(fClingParameters.size()) );
for(Int_t i = 0; i < size; ++i)
{
TString name = TString::Format("%d",i);
SetParameter(name,params[i]);
}
return;
}
fAllParametersSetted = true;
std::copy(params, params+size, fClingParameters.begin() );
}
void TFormula::SetParameters(const Double_t *params)
{
DoSetParameters(params,fNpar);
}
void TFormula::SetParameters(Double_t p0,Double_t p1,Double_t p2,Double_t p3,Double_t p4,
Double_t p5,Double_t p6,Double_t p7,Double_t p8,
Double_t p9,Double_t p10)
{
if(fNpar >= 1) SetParameter(0,p0);
if(fNpar >= 2) SetParameter(1,p1);
if(fNpar >= 3) SetParameter(2,p2);
if(fNpar >= 4) SetParameter(3,p3);
if(fNpar >= 5) SetParameter(4,p4);
if(fNpar >= 6) SetParameter(5,p5);
if(fNpar >= 7) SetParameter(6,p6);
if(fNpar >= 8) SetParameter(7,p7);
if(fNpar >= 9) SetParameter(8,p8);
if(fNpar >= 10) SetParameter(9,p9);
if(fNpar >= 11) SetParameter(10,p10);
}
void TFormula::SetParameter(Int_t param, Double_t value)
{
if (param < 0 || param >= fNpar) return;
assert(int(fClingParameters.size()) == fNpar);
fClingParameters[param] = value;
}
void TFormula::SetParNames(const char *name0,const char *name1,const char *name2,const char *name3,
const char *name4, const char *name5,const char *name6,const char *name7,
const char *name8,const char *name9,const char *name10)
{
if(fNpar >= 1) SetParName(0,name0);
if(fNpar >= 2) SetParName(1,name1);
if(fNpar >= 3) SetParName(2,name2);
if(fNpar >= 4) SetParName(3,name3);
if(fNpar >= 5) SetParName(4,name4);
if(fNpar >= 6) SetParName(5,name5);
if(fNpar >= 7) SetParName(6,name6);
if(fNpar >= 8) SetParName(7,name7);
if(fNpar >= 9) SetParName(8,name8);
if(fNpar >= 10) SetParName(9,name9);
if(fNpar >= 11) SetParName(10,name10);
}
void TFormula::SetParName(Int_t ipar, const char * name)
{
if (ipar < 0 || ipar > fNpar) {
Error("SetParName","Wrong Parameter index %d ",ipar);
return;
}
TString oldName;
for ( auto &it : fParams) {
if (it.second == ipar) {
oldName = it.first;
fParams.erase(oldName);
fParams.insert(std::make_pair(name, ipar) );
break;
}
}
if (oldName.IsNull() ) {
Error("SetParName","Parameter %d is not existing.",ipar);
return;
}
ReplaceParamName(fFormula, oldName, name);
}
void TFormula::ReplaceParamName(TString & formula, const TString & oldName, const TString & name){
if (!formula.IsNull() ) {
bool found = false;
for(list<TFormulaFunction>::iterator it = fFuncs.begin(); it != fFuncs.end(); ++it)
{
if(oldName == it->GetName())
{
found = true;
it->fName = name;
break;
}
}
if(!found)
{
Error("SetParName","Parameter %s is not defined.",oldName.Data());
return;
}
TString newName = name;
newName.ReplaceAll(" ","\\s");
TString pattern = TString::Format("[%s]",oldName.Data());
TString replacement = TString::Format("[%s]",newName.Data());
formula.ReplaceAll(pattern,replacement);
}
}
Double_t TFormula::EvalPar(const Double_t *x,const Double_t *params) const
{
return DoEval(x, params);
}
Double_t TFormula::Eval(Double_t x, Double_t y, Double_t z, Double_t t) const
{
double xxx[4] = {x,y,z,t};
return DoEval(xxx);
}
Double_t TFormula::Eval(Double_t x, Double_t y , Double_t z) const
{
double xxx[3] = {x,y,z};
return DoEval(xxx);
}
Double_t TFormula::Eval(Double_t x, Double_t y) const
{
double xxx[2] = {x,y};
return DoEval(xxx);
}
Double_t TFormula::Eval(Double_t x) const
{
double * xxx = &x;
return DoEval(xxx);
}
Double_t TFormula::DoEval(const double * x, const double * params) const
{
if(!fReadyToExecute)
{
Error("Eval","Formula is invalid and not ready to execute ");
for(auto it = fFuncs.begin(); it != fFuncs.end(); ++it)
{
TFormulaFunction fun = *it;
if(!fun.fFound)
{
printf("%s is unknown.\n",fun.GetName());
}
}
return TMath::QuietNaN();
}
if (!fClingInitialized) {
Error("Eval","Formula is invalid or not properly initialized - try calling TFormula::Compile");
return TMath::QuietNaN();
#ifdef EVAL_IS_NOT_CONST
TString oldClingName = fClingName;
fClingName.Replace(fClingName.Index("_0x")+1,fClingName.Length(), TString::Format("%p",this) );
fClingInput.ReplaceAll(oldClingName, fClingName);
InputFormulaIntoCling();
#endif
}
Double_t result = 0;
void* args[2];
double * vars = (x) ? const_cast<double*>(x) : const_cast<double*>(fClingVariables.data());
args[0] = &vars;
if (fNpar <= 0)
(*fFuncPtr)(0, 1, args, &result);
else {
double * pars = (params) ? const_cast<double*>(params) : const_cast<double*>(fClingParameters.data());
args[1] = &pars;
(*fFuncPtr)(0, 2, args, &result);
}
return result;
}
TString TFormula::GetExpFormula(Option_t *option) const
{
TString opt(option);
if (opt.IsNull() ) return fFormula;
opt.ToUpper();
if (opt.Contains("CLING") ) {
std::string clingFunc = fClingInput.Data();
std::size_t found = clingFunc.find("return");
std::size_t found2 = clingFunc.rfind(";");
if (found == std::string::npos || found2 == std::string::npos) {
Error("GetExpFormula","Invalid Cling expression - return default formula expression");
return fFormula;
}
TString clingFormula = fClingInput(found+7,found2-found-7);
if (!opt.Contains("P")) return clingFormula;
int i = 0;
while (i < clingFormula.Length()-2 ) {
if (clingFormula[i] == 'p' && clingFormula[i+1] == '[' && isdigit(clingFormula[i+2]) ) {
int j = i+3;
while ( isdigit(clingFormula[j]) ) { j++;}
if (clingFormula[j] != ']') {
Error("GetExpFormula","Parameters not found - invalid expression - return default cling formula");
return clingFormula;
}
TString parNumbName = clingFormula(i+2,j-i-2);
int parNumber = parNumbName.Atoi();
assert(parNumber < fNpar);
TString replacement = TString::Format("%f",GetParameter(parNumber));
clingFormula.Replace(i,j-i+1, replacement );
i += replacement.Length();
}
i++;
}
return clingFormula;
}
if (opt.Contains("P") ) {
TString expFormula = fFormula;
int i = 0;
while (i < expFormula.Length()-2 ) {
if (expFormula[i] == '[') {
int j = i+1;
while ( expFormula[j] != ']' ) { j++;}
if (expFormula[j] != ']') {
Error("GetExpFormula","Parameter names not found - invalid expression - return default formula");
return expFormula;
}
TString parName = expFormula(i+1,j-i-1);
TString replacement = TString::Format("%g",GetParameter(parName));
expFormula.Replace(i,j-i+1, replacement );
i += replacement.Length();
}
i++;
}
return expFormula;
}
Warning("GetExpFormula","Invalid option - return defult formula expression");
return fFormula;
}
void TFormula::Print(Option_t *option) const
{
printf(" %20s : %s Ndim= %d, Npar= %d, Number= %d \n",GetName(),GetTitle(), fNdim,fNpar,fNumber);
printf(" Formula expression: \n");
printf("\t%s \n",fFormula.Data() );
TString opt(option);
opt.ToUpper();
if (opt.Contains("V") ) {
if (fNdim > 0) {
printf("List of Variables: \n");
assert(int(fClingVariables.size()) >= fNdim);
for ( int ivar = 0; ivar < fNdim ; ++ivar) {
printf("Var%4d %20s = %10f \n",ivar,GetVarName(ivar).Data(), fClingVariables[ivar]);
}
}
if (fNpar > 0) {
printf("List of Parameters: \n");
if ( int(fClingParameters.size()) < fNpar)
Error("Print","Number of stored parameters in vector %lu in map %lu is different than fNpar %d",fClingParameters.size(), fParams.size(), fNpar);
assert(int(fClingParameters.size()) >= fNpar);
for ( int ipar = 0; ipar < fNpar ; ++ipar) {
printf("Par%4d %20s = %10f \n",ipar,GetParName(ipar), fClingParameters[ipar] );
}
}
printf("Expression passed to Cling:\n");
printf("\t%s\n",fClingInput.Data() );
}
if(!fReadyToExecute)
{
Warning("Print","Formula is not ready to execute. Missing parameters/variables");
for(list<TFormulaFunction>::const_iterator it = fFuncs.begin(); it != fFuncs.end(); ++it)
{
TFormulaFunction fun = *it;
if(!fun.fFound)
{
printf("%s is unknown.\n",fun.GetName());
}
}
}
if(!fAllParametersSetted)
{
}
}
void TFormula::Streamer(TBuffer &b)
{
if (b.IsReading() ) {
UInt_t R__s, R__c;
Version_t v = b.ReadVersion(&R__s, &R__c);
if (v <= 8 && v > 3 && v != 6) {
ROOT::v5::TFormula * fold = new ROOT::v5::TFormula();
fold->Streamer(b, v, R__s, R__c, TFormula::Class());
TFormula fnew(fold->GetName(), fold->GetExpFormula() );
*this = fnew;
SetParameters(fold->GetParameters() );
if (!fReadyToExecute ) {
Error("Streamer","Old formula read from file is NOT valid");
Print("v");
}
delete fold;
return;
}
else if (v > 8) {
b.ReadClassBuffer(TFormula::Class(), this, v, R__s, R__c);
if (fFormula.IsNull() ) return;
std::vector<double> parValues = fClingParameters;
auto paramMap = fParams;
fNpar = fParams.size();
fClingParameters.clear();
FillDefaults();
PreProcessFormula(fFormula);
PrepareFormula(fFormula);
if (fNpar != (int) parValues.size() ) {
Error("Streamer","number of parameters computed (%d) is not same as the stored parameters (%d)",fNpar,int(parValues.size()) );
Print("v");
}
assert(fNpar == (int) parValues.size() );
std::copy( parValues.begin(), parValues.end(), fClingParameters.begin() );
if (fParams.size() != paramMap.size() ) {
Warning("Streamer","number of parameters list found (%lu) is not same as the stored one (%lu) - use re-created list",fParams.size(),paramMap.size()) ;
}
else
fParams = paramMap;
if (!TestBit(kNotGlobal)) {
R__LOCKGUARD2(gROOTMutex);
gROOT->GetListOfFunctions()->Add(this);
}
if (!fReadyToExecute ) {
Error("Streamer","Formula read from file is NOT ready to execute");
Print("v");
}
return;
}
else {
Error("Streamer","Reading version %d is not supported",v);
return;
}
}
else {
b.WriteClassBuffer(TFormula::Class(),this);
}
}