---

ROOT-R Example: Mixing ROOT and R¶

In [1]:

auto r = ROOT::R::TRInterface::Instance();

Minimisation:¶

A simple example using random points with Gaussian distribution. We use the interactive ROOT visualisation, JSROOT.

In [2]:

%jsroot on
auto c = new TCanvas("c");
TH1D h1("h1", "h1", 256, -4, 4);
h1.FillRandom("gaus");
ROOT::Math::MinimizerOptions::SetDefaultMinimizer("RMinimizer","BFGS");
h1.Fit("gaus");
c->Draw();

Calling R with command result <- optim( initialparams, minfunction,method='BFGS',control = list(ndeps=stepsizes,maxit=1000,trace=0,abstol=1.000000e-02),hessian=TRUE)
Value at minimum =247.533

****************************************
Minimizer is RMinimizer / BFGS
Chi2                      =      247.533
NDf                       =          197
NCalls                    =          110
Constant                  =      68.7203   +/-   4.08237     
Mean                      =   0.00966072   +/-   0.0118658   
Sigma                     =     0.943347   +/-   0.03857         (limited)

Fit in R and plot in ROOT¶

The next example creates an exponential fit. The idea is to create a set of numbers x,y plus noise using ROOT, pass them to R, fit the data using nls(Nonlinear Least Squares), obtain fitted coefficient and then plot data, known function and fitted function with ROOT

In [3]:

// draw a frame to define the range
TMultiGraph mg;
const Int_t n = 24;
Double_t x[n] ;
Double_t y[n] ;

Create data¶

Generate random points from $X^3$ with noise

In [4]:

TRandom rg;
rg.SetSeed(520);
for (Int_t i = 0; i < n; i++) {
      x[i] = rg.Uniform(0, 1);
      y[i] = TMath::Power(x[i], 3) + rg.Gaus() * 0.06;
}

Plot¶

In [5]:

%%cpp -d
void SetPointStyle(TGraph& g, Int_t color)
{
   g.SetMarkerColor(color);
   g.SetFillColor(kWhite);
   g.SetMarkerStyle(8);
   g.SetMarkerSize(1);    
}

In [6]:

TGraph gr1(n,x,y);
SetPointStyle(gr1, kBlue);
mg.Add(&gr1);

TF1 f_known("f_known","pow(x,3)",0,1);
TGraph gr2(&f_known);
SetPointStyle(gr2, kRed);
mg.Add(&gr2);

Passing Data to R environment¶

passing data to R for fitting

In [7]:

TVectorD vx(n, x);
TVectorD vy(n, y);
r["x"]=vx;
r["y"]=vy;
r<<"ds<-data.frame(x=x,y=y)";

Fitting with Nonlinear Least Squares¶

In [8]:

r << "m <- nls(y ~ I(x^power),data = ds, start = list(power = 1),trace = T)";

1.506598 :  1
0.3263132 :  1.885851
0.09498773 :  2.719307
0.0764783 :  3.086762
0.07622443 :  3.136312
0.07622369 :  3.139027
0.07622369 :  3.139151

In [9]:

Double_t power;
r["summary(m)$coefficients[1]"] >> power;

Plot in ROOT¶

Plot fitted function

In [10]:

TCanvas c1("c1", "Curve Fitting", 700, 500);
c1.SetGrid();

TF1 f_fitted("f_fitted","pow(x,[0])",0,1);
f_fitted.SetParameter(0,power);

TGraph gr3(&f_fitted);
SetPointStyle(gr3, kGreen);
mg.Add(&gr3);

mg.Draw("ap");

TLegend leg(0.1,0.6,0.5,0.9,"brNDC");
leg.SetHeader("Fitting x^power");
leg.AddEntry(&gr1, "Points with gaussian noise to be fitted", "P");
leg.AddEntry(&gr2, "Known function x^3", "P");
TString fmsg;
fmsg.Form(" \"Green\"  Fitted function with power=%.4lf",power);
leg.AddEntry(&gr3, fmsg, "P");
leg.Draw();   

c1.Draw();