#include "TF1Convolution.h"
#include "Riostream.h"
#include "TROOT.h"
#include "TObject.h"
#include "TObjString.h"
#include "TMath.h"
#include "Math/Integrator.h"
#include "Math/IntegratorMultiDim.h"
#include "Math/IntegratorOptions.h"
#include "Math/GaussIntegrator.h"
#include "Math/GaussLegendreIntegrator.h"
#include "Math/AdaptiveIntegratorMultiDim.h"
#include "Math/Functor.h"
#include "TVirtualFFT.h"
#include "TClass.h"
class TF1Convolution_EvalWrapper
{
std::shared_ptr < TF1 > fFunction1;
std::shared_ptr < TF1 > fFunction2;
Double_t fT0;
public:
TF1Convolution_EvalWrapper(std::shared_ptr<TF1> & f1 , std::shared_ptr<TF1> & f2, Double_t t)
: fFunction1(f1), fFunction2(f2), fT0(t)
{
}
Double_t operator()(Double_t x) const
{
return fFunction1->Eval(x) * fFunction2->Eval(x-fT0);
}
};
void TF1Convolution::InitializeDataMembers(TF1* function1, TF1* function2, Bool_t useFFT)
{
if (function1) {
TF1 * fnew1 = (TF1*) function1->IsA()->New();
function1->Copy(*fnew1);
fFunction1 = std::shared_ptr<TF1>(fnew1);
}
if (function2) {
TF1 * fnew2 = (TF1*) function2->IsA()->New();
function2->Copy(*fnew2);
fFunction2 = std::shared_ptr<TF1>(fnew2);
}
if (fFunction1.get() == nullptr|| fFunction2.get() == nullptr)
Fatal("InitializeDataMembers","Invalid functions - Abort");
fFunction1->GetRange(fXmin, fXmax);
Double_t range = fXmax - fXmin;
fXmin -= 0.1*range;
fXmax += 0.1*range;
fNofParams1 = fFunction1->GetNpar();
fNofParams2 = fFunction2->GetNpar();
fParams1 = std::vector<Double_t>(fNofParams1);
fParams2 = std::vector<Double_t>(fNofParams2);
fCstIndex = fFunction2-> GetParNumber("Constant");
fFlagFFT = useFFT;
fFlagGraph = false;
fNofPoints = 10000;
fParNames.reserve( fNofParams1 + fNofParams2);
for (int i=0; i<fNofParams1; i++)
{
fParams1[i] = fFunction1 -> GetParameter(i);
fParNames.push_back(fFunction1 -> GetParName(i) );
}
for (int i=0; i<fNofParams2; i++)
{
fParams2[i] = fFunction2 -> GetParameter(i);
if (i != fCstIndex) fParNames.push_back(fFunction2 -> GetParName(i) );
}
if (fCstIndex!=-1)
{
fFunction2 -> FixParameter(fCstIndex,1.);
fNofParams2 = fNofParams2-1;
fParams2.erase(fParams2.begin()+fCstIndex);
}
}
TF1Convolution::TF1Convolution(TF1* function1, TF1* function2, Bool_t useFFT)
{
InitializeDataMembers(function1,function2, useFFT);
}
TF1Convolution::TF1Convolution(TF1* function1, TF1* function2, Double_t xmin, Double_t xmax, Bool_t useFFT)
{
InitializeDataMembers(function1, function2,useFFT);
if (xmin < xmax) {
fXmin = xmin;
fXmax = xmax;
}
}
TF1Convolution::TF1Convolution(TString formula, Double_t xmin, Double_t xmax, Bool_t useFFT)
{
TF1::InitStandardFunctions();
TObjArray *objarray = formula.Tokenize("*");
std::vector < TString > stringarray(2);
std::vector < TF1* > funcarray(2);
for (int i=0; i<2; i++)
{
stringarray[i] = ((TObjString*)((*objarray)[i])) -> GetString();
stringarray[i].ReplaceAll(" ","");
funcarray[i] = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(stringarray[i]));
if (funcarray[i] == nullptr) {
TF1 * f = new TF1(TString::Format("f_conv_%d",i+1),stringarray[i]);
if (!f->GetFormula()->IsValid() )
Error("TF1Convolution","Invalid formula : %s",stringarray[i].Data() );
if (i == 0)
fFunction1 = std::shared_ptr<TF1>(f);
else
fFunction2 = std::shared_ptr<TF1>(f);
}
}
InitializeDataMembers(funcarray[0], funcarray[1],useFFT);
if (xmin < xmax) {
fXmin = xmin;
fXmax = xmax;
}
}
TF1Convolution::TF1Convolution(TString formula1, TString formula2, Double_t xmin, Double_t xmax, Bool_t useFFT)
{
TF1::InitStandardFunctions();
(TString)formula1.ReplaceAll(" ","");
(TString)formula2.ReplaceAll(" ","");
TF1* f1 = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(formula1));
TF1* f2 = (TF1*)(gROOT -> GetListOfFunctions() -> FindObject(formula2));
if (!f1) {
fFunction1 = std::shared_ptr<TF1>(new TF1("f_conv_1",formula1) );
if (!fFunction1->GetFormula()->IsValid() )
Error("TF1Convolution","Invalid formula for : %s",formula1.Data() );
}
if (!f2) {
fFunction2 = std::shared_ptr<TF1>(new TF1("f_conv_1",formula2) );
if (!fFunction2->GetFormula()->IsValid() )
Error("TF1Convolution","Invalid formula for : %s",formula2.Data() );
}
InitializeDataMembers(f1, f2,useFFT);
if (xmin < xmax) {
fXmin = xmin;
fXmax = xmax;
}
}
void TF1Convolution::MakeFFTConv()
{
if (gDebug)
Info("MakeFFTConv","Making FFT convolution using %d points in range [%g,%g]",fNofPoints,fXmin,fXmax);
std::vector < Double_t > x (fNofPoints);
std::vector < Double_t > in1(fNofPoints);
std::vector < Double_t > in2(fNofPoints);
TVirtualFFT *fft1 = TVirtualFFT::FFT(1, &fNofPoints, "R2C K");
TVirtualFFT *fft2 = TVirtualFFT::FFT(1, &fNofPoints, "R2C K");
if (fft1 == nullptr || fft2 == nullptr) {
Warning("MakeFFTConv","Cannot use FFT, probably FFTW# package is not available. Switch to numerical convolution");
fFlagFFT = false;
}
Double_t shift2 = 0.5*(fXmin+fXmax);
Double_t x2;
for (int i=0; i<fNofPoints; i++)
{
x[i] = fXmin + (fXmax-fXmin)/(fNofPoints-1)*i;
x2 = x[i] - shift2;
in1[i] = fFunction1 -> EvalPar( &x[i], nullptr);
in2[i] = fFunction2 -> EvalPar( &x2, nullptr);
fft1 -> SetPoint(i, in1[i]);
fft2 -> SetPoint(i, in2[i]);
}
fft1 -> Transform();
fft2 -> Transform();
TVirtualFFT *fftinverse = TVirtualFFT::FFT(1, &fNofPoints, "C2R K");
Double_t re1, re2, im1, im2, out_re, out_im;
for (int i=0;i<=fNofPoints/2.;i++)
{
fft1 -> GetPointComplex(i,re1,im1);
fft2 -> GetPointComplex(i,re2,im2);
out_re = re1*re2 - im1*im2;
out_im = re1*im2 + re2*im1;
fftinverse -> SetPoint(i, out_re, out_im);
}
fftinverse -> Transform();
if (!fGraphConv) fGraphConv = std::shared_ptr< TGraph >(new TGraph(fNofPoints));
for (int i=0;i<fNofPoints;i++)
{
int j = i + fNofPoints/2;
if (j >= fNofPoints) j -= fNofPoints;
fGraphConv->SetPoint(i, x[i], fftinverse->GetPointReal(j)/fNofPoints);
}
fFlagGraph = true;
}
Double_t TF1Convolution::EvalFFTConv(Double_t t)
{
if (!fFlagGraph) MakeFFTConv();
if (fGraphConv)
return fGraphConv -> Eval(t);
else
return EvalNumConv(t);
}
Double_t TF1Convolution::EvalNumConv(Double_t t)
{
TF1Convolution_EvalWrapper fconv( fFunction1, fFunction2, t);
Double_t result = 0;
if (ROOT::Math::IntegratorOneDimOptions::DefaultIntegratorType() == ROOT::Math::IntegrationOneDim::kGAUSS )
{
ROOT::Math::GaussIntegrator integrator(1e-9, 1e-9);
integrator.SetFunction(ROOT::Math::Functor1D(fconv));
if (fXmin != - TMath::Infinity() && fXmax != TMath::Infinity())
result = integrator.Integral(fXmin, fXmax);
else if (fXmin == - TMath::Infinity() && fXmax != TMath::Infinity())
result = integrator.IntegralLow(fXmax);
else if (fXmin != - TMath::Infinity() && fXmax == TMath::Infinity())
result = integrator.IntegralUp(fXmin);
else if (fXmin == - TMath::Infinity() && fXmax == TMath::Infinity())
result = integrator.Integral();
}
else
{
ROOT::Math::IntegratorOneDim integrator(fconv, ROOT::Math::IntegratorOneDimOptions::DefaultIntegratorType(), 1e-9, 1e-9);
if (fXmin != - TMath::Infinity() && fXmax != TMath::Infinity() )
result = integrator.Integral(fXmin, fXmax);
else if (fXmin == - TMath::Infinity() && fXmax != TMath::Infinity() )
result = integrator.IntegralLow(fXmax);
else if (fXmin != - TMath::Infinity() && fXmax == TMath::Infinity() )
result = integrator.IntegralUp(fXmin);
else if (fXmin == - TMath::Infinity() && fXmax == TMath::Infinity() )
result = integrator.Integral();
}
return result;
}
Double_t TF1Convolution::operator()(Double_t* t, Double_t* p)
{
if (p!=0) TF1Convolution::SetParameters(p);
Double_t result = 0.;
if (fFlagFFT) result = EvalFFTConv(t[0]);
else result = EvalNumConv(t[0]);
return result;
}
void TF1Convolution::SetNofPointsFFT(Int_t n)
{
if (n<0) return;
fNofPoints = n;
if (fGraphConv) fGraphConv -> Set(fNofPoints);
fFlagGraph = false;
}
void TF1Convolution::SetParameters(Double_t* p)
{
bool equalParams = true;
for (int i=0; i<fNofParams1; i++) {
fFunction1 -> SetParameter(i,p[i]);
equalParams &= ( fParams1[i] == p[i] );
fParams1[i] = p[i];
}
Int_t k = 0;
Int_t offset = 0;
Int_t offset2 = 0;
if (fCstIndex!=-1) offset = 1;
Int_t totalnofparams = fNofParams1+fNofParams2+offset;
for (int i=fNofParams1; i<totalnofparams; i++) {
if (k==fCstIndex)
{
k++;
offset2=1;
continue;
}
fFunction2 -> SetParameter(k,p[i-offset2]);
equalParams &= ( fParams2[k-offset2] == p[i-offset2] );
fParams2[k-offset2] = p[i-offset2];
k++;
}
if (!equalParams) fFlagGraph = false;
}
void TF1Convolution::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 params[]={p0,p1,p2,p3,p4,p5,p6,p7};
TF1Convolution::SetParameters(params);
}
void TF1Convolution::SetExtraRange(Double_t percentage)
{
if (percentage<0) return;
double range = fXmax = fXmin;
fXmin -= percentage * range;
fXmax += percentage * range;
fFlagGraph = false;
}
void TF1Convolution::SetRange(Double_t a, Double_t b)
{
if (a>=b) return;
fXmin = a;
fXmax = b;
if (fFlagFFT && ( a==-TMath::Infinity() || b==TMath::Infinity() ) )
{
Warning("TF1Convolution::SetRange()","In FFT mode, range can not be infinite. Infinity has been replaced by range of first function plus a bufferzone to avoid spillover.");
if (a ==-TMath::Infinity()) fXmin = fFunction1 -> GetXmin();
if ( b== TMath::Infinity()) fXmax = fFunction1 -> GetXmax();
SetExtraRange(0.1);
}
fFlagGraph = false;
}