#include "TROOT.h"
#include "TClass.h"
#include "TMath.h"
#include "TF1NormSum.h"
#include "Math/WrappedFunction.h"
#include "Math/WrappedTF1.h"
template<class Iterator>
void FixDuplicateNames(Iterator begin, Iterator end) {
std::multimap<TString, int > parMap;
for (Iterator it = begin; it != end; ++it) {
parMap.insert( std::make_pair( *it, std::distance(begin,it) ) );
}
for ( auto & elem : parMap) {
TString name = elem.first;
int n = parMap.count( name);
if (n > 1 ) {
std::pair <std::multimap<TString,int>::iterator, std::multimap<TString,int>::iterator> ret;
ret = parMap.equal_range(name);
int i = 0;
for (std::multimap<TString,int>::iterator it=ret.first; it!=ret.second; ++it) {
*(begin+it->second) = TString::Format("%s%d",name.Data(),++i);
}
}
}
}
void TF1NormSum::InitializeDataMembers(const std::vector <std::shared_ptr < TF1 >> &functions, const std::vector <Double_t> &coeffs, Double_t scale)
{
fScale = scale;
fFunctions = functions;
fCoeffs = coeffs;
fNOfFunctions = functions.size();
fCstIndexes = std::vector < Int_t > (fNOfFunctions);
fParNames = std::vector<TString> (fNOfFunctions);
fParNames.reserve(3*fNOfFunctions);
for (unsigned int n=0; n < fNOfFunctions; n++)
{
int npar = fFunctions[n] -> GetNpar();
fCstIndexes[n] = fFunctions[n] -> GetParNumber("Constant");
fParNames[n] = TString::Format("Coeff%d",n);
if (fCstIndexes[n]!= -1)
{
fFunctions[n] -> FixParameter(fCstIndexes[n], 1.);
int k = 0;
for (int i=0; i<npar; i++)
{
if (i==fCstIndexes[n]) continue;
fParNames.push_back( fFunctions[n] -> GetParName(i) );
k++;
}
}
else {
for (int i=0; i < npar; i++)
{
fParNames.push_back( fFunctions[n] -> GetParName(i) );
}
}
if (!fFunctions[n] -> IsEvalNormalized()) fFunctions[n] -> SetNormalized(true);
}
FixDuplicateNames(fParNames.begin()+fNOfFunctions, fParNames.end());
}
TF1NormSum::TF1NormSum()
{
fNOfFunctions = 0;
fScale = 1.;
fFunctions = std::vector< std::shared_ptr < TF1 >>(0) ;
fCoeffs = std::vector < Double_t >(0) ;
fCstIndexes = std::vector < Int_t > (0);
}
TF1NormSum::TF1NormSum(const std::vector <TF1*> &functions, const std::vector <Double_t> &coeffs, Double_t scale)
{
std::vector <std::shared_ptr < TF1 > >f;
for (unsigned int i = 0; i<functions.size(); i++)
{
f[i] = std::shared_ptr < TF1 >((TF1*)functions[i]->Clone());
}
InitializeDataMembers(f,coeffs,scale);
}
TF1NormSum::TF1NormSum(TF1* function1, TF1* function2, Double_t coeff1, Double_t coeff2, Double_t scale)
{
std::vector < std::shared_ptr < TF1 > > functions(2);
std::vector < Double_t > coeffs(2);
TF1 * fnew1 = 0;
TF1 * fnew2 = 0;
if (function1) {
fnew1 = (TF1*) function1->IsA()->New();
function1->Copy(*fnew1);
}
if (function2) {
fnew2 = (TF1*) function2->IsA()->New();
function2->Copy(*fnew2);
}
if (fnew1 == nullptr || fnew2 == nullptr)
Fatal("TF1NormSum","Invalid input functions - Abort");
std::shared_ptr < TF1 > f1( fnew1);
std::shared_ptr < TF1 > f2( fnew2);
functions = {f1, f2};
coeffs = {coeff1, coeff2};
InitializeDataMembers(functions, coeffs,scale);
}
TF1NormSum::TF1NormSum(TF1* function1, TF1* function2, TF1* function3, Double_t coeff1, Double_t coeff2, Double_t coeff3, Double_t scale)
{
std::vector < std::shared_ptr < TF1 > > functions(3);
std::vector < Double_t > coeffs(3);
TF1 * fnew1 = 0;
TF1 * fnew2 = 0;
TF1 * fnew3 = 0;
if (function1) {
fnew1 = (TF1*) function1->IsA()->New();
function1->Copy(*fnew1);
}
if (function2) {
fnew2 = (TF1*) function2->IsA()->New();
function2->Copy(*fnew2);
}
if (function3) {
fnew3 = (TF1*) function3->IsA()->New();
function3->Copy(*fnew2);
}
if (!fnew1 || !fnew2 || !fnew3 )
Fatal("TF1NormSum","Invalid input functions - Abort");
std::shared_ptr < TF1 > f1( fnew1);
std::shared_ptr < TF1 > f2( fnew2);
std::shared_ptr < TF1 > f3( fnew3);
functions = {f1, f2, f3};
coeffs = {coeff1, coeff2, coeff3};
InitializeDataMembers(functions, coeffs,scale);
}
TF1NormSum::TF1NormSum(const TString &formula, Double_t xmin, Double_t xmax)
{
TF1::InitStandardFunctions();
TObjArray *arrayall = formula.Tokenize("*+");
TObjArray *arraytimes = formula.Tokenize("*") ;
Int_t noffunctions = (formula.Tokenize("+")) -> GetEntries();
Int_t nofobj = arrayall -> GetEntries();
Int_t nofcoeffs = nofobj - noffunctions;
std::vector < std::shared_ptr < TF1 > > functions(noffunctions);
std::vector < Double_t > coeffs(noffunctions);
std::vector < TString > funcstringall(nofobj);
std::vector < Int_t > indexsizetimes(nofcoeffs+1);
std::vector < Bool_t > isacoeff(nofobj);
for (int i=0; i<nofobj; i++)
{
funcstringall[i] = ((TObjString*)((*arrayall)[i])) -> GetString();
funcstringall[i].ReplaceAll(" ","");
}
Int_t j = 0;
Int_t k = 1;
for (int i=0; i<nofcoeffs+1; i++)
{
indexsizetimes[i] = ( ( ( (TObjString*)(*arraytimes)[i] ) -> GetString() ).Tokenize("+") ) -> GetEntries();
while (k < indexsizetimes[i])
{
isacoeff[k+j-1] = 0;
k++;
}
j = j+indexsizetimes[i];
if (j==nofobj) isacoeff[j-1] = 0;
else isacoeff[j-1] = 1;
k = 1;
}
k = 0;
for (int i=0; i<noffunctions; i++)
{
if (isacoeff[k]==0)
{
coeffs[i] = 1.;
TF1* f = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(funcstringall[k]));
if (!f) Error("TF1NormSum", "Function %s does not exist", funcstringall[k].Data());
functions[i] = std::shared_ptr < TF1 > ((TF1*)f->Clone(TString::Format("function_%s_%d",funcstringall[k].Data(), i)));
functions[i]->SetRange(xmin,xmax);
k++;
}
else
{
coeffs[i] = funcstringall[k].Atof();
TF1* f = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(funcstringall[k+1]));
if (!f) Error("TF1NormSum", "Function %s does not exist", funcstringall[k+1].Data());
functions[i] = std::shared_ptr < TF1 >((TF1*)f->Clone(TString::Format("function_%s_%d",funcstringall[k+1].Data(), i) ));
functions[i]->SetRange(xmin,xmax);
k=k+2;
}
}
InitializeDataMembers(functions, coeffs,1.);
}
double TF1NormSum::operator()(double* x, double* p)
{
if (p!=0) TF1NormSum::SetParameters(p);
Double_t sum = 0.;
for (unsigned int n=0; n<fNOfFunctions; n++)
{
sum += fCoeffs[n]*(fFunctions[n] -> EvalPar(x,0));
}
return fScale * sum;
}
std::vector<double> TF1NormSum::GetParameters() const {
std::vector<double> params(GetNpar() );
int offset = 0;
int nOfNonCstParams = 0;
for (unsigned int n=0; n<fNOfFunctions; n++)
{
params[n] = fCoeffs[n];
offset += nOfNonCstParams;
int k = 0;
for (int j = 0; j < fFunctions[n]->GetNpar(); ++j) {
if (j != fCstIndexes[n]) {
params[k+fNOfFunctions+offset] = fFunctions[n]->GetParameter(j);
k++;
}
}
nOfNonCstParams = k;
}
return params;
}
void TF1NormSum::SetParameters(const double* params)
{
for (unsigned int n=0; n<fNOfFunctions; n++)
{
fCoeffs[n] = params[n];
}
Int_t offset = 0;
int k = 0;
for (unsigned int n=0; n<fNOfFunctions; n++)
{
bool equalParams = true;
Double_t * funcParams = fFunctions[n]->GetParameters();
int npar = fFunctions[n]->GetNpar();
offset += k;
k = 0;
for (int i = 0; i < npar; ++i) {
if (i != fCstIndexes[n])
{
equalParams &= (funcParams[i] == params[k+fNOfFunctions+offset] );
funcParams[i] = params[k+fNOfFunctions+offset];
k++;
}
}
if (!equalParams) fFunctions[n]->Update();
}
}
void TF1NormSum::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)
{
const double params[] = {p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10};
TF1NormSum::SetParameters(params);
}
Int_t TF1NormSum::GetNpar() const
{
Int_t nofparams = 0;
for (unsigned int n=0; n<fNOfFunctions; ++n)
{
nofparams += fFunctions[n]->GetNpar();
if (fCstIndexes[n] >= 0) nofparams -= 1;
}
return nofparams + fNOfFunctions;
}