//Numerical integration using R
//passing the function from ROOT
#include<TMath.h>
#include<TRInterface.h>
#include<Math/Integrator.h>
#include<TF1.h>

//To integrate using R the function must be vectorized
//The idea is just to receive a vector like an argument,to evaluate 
//every element saving the result in another vector 
//and return the resultant vector.
std::vector<Double_t>  BreitWignerVectorized(std::vector<Double_t> xx)
{
  std::vector<Double_t> result(xx.size());
  for(Int_t i=0;i<xx.size();i++)
  {
    result[i]=TMath::BreitWigner(xx[i]);
  }
  return result;
}

double BreitWignerWrap( double x){ 
   return TMath::BreitWigner(x);
}


void Integration()
{

  ROOT::R::TRInterface &r=ROOT::R::TRInterface::Instance();
  
  r["BreitWigner"]=ROOT::R::TRFunctionExport(BreitWignerVectorized);
  
  Double_t value=r.Eval("integrate(BreitWigner, lower = -2, upper = 2)$value");

  std::cout.precision(18);
  std::cout<<"Integral of the BreitWigner Function in the interval [-2, 2] R        = "<<value<<std::endl;

  
  ROOT::Math::WrappedFunction<> wf(BreitWignerWrap);
  ROOT::Math::Integrator i(wf);
  value=i.Integral(-2,2);
  std::cout<<"Integral of the BreitWigner Function in the interval [-2, 2] MathMore = "<<value<<std::endl;

  
  TF1 f1("BreitWigner","BreitWignerWrap(x)");
  value=f1.Integral(-2,2);
  std::cout<<"Integral of the BreitWigner Function in the interval [-2, 2] TF1      = "<<value<<std::endl;
  
  //infinte limits
  value=r.Eval("integrate(BreitWigner, lower = -Inf, upper = Inf)$value");
  std::cout<<"Integral of BreitWigner Function in the interval [-Inf, Inf] R    = "<<value<<std::endl;

}