/// \file
/// \ingroup tutorial_r
/// \notebook -nodraw
/// Numerical integration using R passing the function from ROOT
///
/// \macro_code
///
/// \author

#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;

   // infinite 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;
}