Logo ROOT  
Reference Guide
mathmoreIntegration.C File Reference

Detailed Description

View in nbviewer Open in SWAN Example on the usage of the adaptive 1D integration algorithm of MathMore it calculates the numerically cumulative integral of a distribution (like in this case the BreitWigner) to execute the macro type it (you need to compile with AClic)

root[0] .x mathmoreIntegration.C+

This tutorials require having libMathMore built with ROOT.

To build mathmore you need to have a version of GSL >= 1.8 installed in your system The ROOT configure will automatically find GSL if the script gsl-config (from GSL) is in your PATH,. otherwise you need to configure root with the options –gsl-incdir and –gsl-libdir.

***************************************************************
Drawing cumulatives of BreitWigner in interval [ -2 , 2 ]
***************************************************************
***************************************************************
Test integration performances in interval [ -2 , 2 ]
Time using ROOT::Math::Integrator : 0.14291
Number of function calls = 69
42201.6649413923442
Time using TF1::Integral : 0.481602
Number of function calls = 91
42201.6649413923442
#include "TMath.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TLegend.h"
/*#include "TLabel.h"*/
#include "Math/Functor.h"
#include "Math/IFunction.h"
#include <iostream>
#include "TStopwatch.h"
#include "TF1.h"
#include <limits>
//!calculates exact integral of Breit Wigner distribution
//!and compares with existing methods
int nc = 0;
double exactIntegral( double a, double b) {
return (TMath::ATan(2*b)- TMath::ATan(2*a))/ TMath::Pi();
}
double func( double x){
nc++;
}
// TF1 requires the function to have the ( )( double *, double *) signature
double func2(const double *x, const double * = 0){
nc++;
return TMath::BreitWigner(x[0]);
}
void testIntegPerf(double x1, double x2, int n = 100000){
std::cout << "\n\n***************************************************************\n";
std::cout << "Test integration performances in interval [ " << x1 << " , " << x2 << " ]\n\n";
TStopwatch timer;
double dx = (x2-x1)/double(n);
//ROOT::Math::Functor1D<ROOT::Math::IGenFunction> f1(& TMath::BreitWigner);
timer.Start();
double s1 = 0.0;
nc = 0;
for (int i = 0; i < n; ++i) {
double x = x1 + dx*i;
s1+= ig.Integral(x1,x);
}
timer.Stop();
std::cout << "Time using ROOT::Math::Integrator :\t" << timer.RealTime() << std::endl;
std::cout << "Number of function calls = " << nc/n << std::endl;
int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
//TF1 *fBW = new TF1("fBW","TMath::BreitWigner(x)",x1, x2); // this is faster but cannot measure number of function calls
TF1 *fBW = new TF1("fBW",func2,x1, x2,0);
timer.Start();
nc = 0;
double s2 = 0;
for (int i = 0; i < n; ++i) {
double x = x1 + dx*i;
s2+= fBW->Integral(x1,x );
}
timer.Stop();
std::cout << "Time using TF1::Integral :\t\t\t" << timer.RealTime() << std::endl;
std::cout << "Number of function calls = " << nc/n << std::endl;
pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
}
void DrawCumulative(double x1, double x2, int n = 100){
std::cout << "\n\n***************************************************************\n";
std::cout << "Drawing cumulatives of BreitWigner in interval [ " << x1 << " , " << x2 << " ]\n\n";
double dx = (x2-x1)/double(n);
TH1D *cum0 = new TH1D("cum0", "", n, x1, x2); //exact cumulative
for (int i = 1; i <= n; ++i) {
double x = x1 + dx*i;
cum0->SetBinContent(i, exactIntegral(x1, x));
}
// alternative method using ROOT::Math::Functor class
TH1D *cum1 = new TH1D("cum1", "", n, x1, x2);
for (int i = 1; i <= n; ++i) {
double x = x1 + dx*i;
cum1->SetBinContent(i, ig.Integral(x1,x));
}
TF1 *fBW = new TF1("fBW","TMath::BreitWigner(x, 0, 1)",x1, x2);
TH1D *cum2 = new TH1D("cum2", "", n, x1, x2);
for (int i = 1; i <= n; ++i) {
double x = x1 + dx*i;
cum2->SetBinContent(i, fBW->Integral(x1,x));
}
TH1D *cum10 = new TH1D("cum10", "", n, x1, x2); //difference between 1 and exact
TH1D *cum20 = new TH1D("cum23", "", n, x1, x2); //difference between 2 and excact
for (int i = 1; i <= n; ++i) {
double delta = cum1->GetBinContent(i) - cum0->GetBinContent(i);
double delta2 = cum2->GetBinContent(i) - cum0->GetBinContent(i);
//std::cout << " diff for " << x << " is " << delta << " " << cum1->GetBinContent(i) << std::endl;
cum10->SetBinContent(i, delta );
cum20->SetBinContent(i, delta2 );
}
TCanvas *c1 = new TCanvas("c1","Integration example",20,10,800,500);
c1->Divide(2,1);
c1->Draw();
cum0->SetTitle("BreitWigner - the cumulative");
cum0->SetStats(0);
cum1->SetLineStyle(2);
cum2->SetLineStyle(3);
c1->cd(1);
cum0->DrawCopy("h");
cum1->DrawCopy("same");
//cum2->DrawCopy("same");
cum2->DrawCopy("same");
c1->cd(2);
cum10->SetTitle("Difference");
cum10->SetStats(0);
cum10->Draw("e0");
cum20->SetLineColor(kRed);
cum20->Draw("hsame");
TLegend * l = new TLegend(0.11, 0.8, 0.7 ,0.89);
l->AddEntry(cum10, "GSL integration - analytical ");
l->AddEntry(cum20, "TF1::Integral - analytical ");
l->Draw();
c1->Update();
std::cout << "\n***************************************************************\n";
}
void mathmoreIntegration(double a = -2, double b = 2)
{
DrawCumulative(a, b);
testIntegPerf(a, b);
}
#define b(i)
Definition: RSha256.hxx:100
#define s1(x)
Definition: RSha256.hxx:91
static const double x2[5]
static const double x1[5]
@ kRed
Definition: Rtypes.h:64
@ kBlack
Definition: Rtypes.h:63
@ kBlue
Definition: Rtypes.h:64
Functor1D class for one-dimensional functions.
Definition: Functor.h:493
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:94
Template class to wrap any C++ callable object which takes one argument i.e.
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TAttLine.h:42
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
The Canvas class.
Definition: TCanvas.h:27
1-Dim function class
Definition: TF1.h:210
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition: TF1.cxx:2505
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6345
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:8662
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8678
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition: TH1.cxx:3045
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8446
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
Stopwatch class.
Definition: TStopwatch.h:28
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
return c1
Definition: legend1.C:41
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
TF1 * f1
Definition: legend1.C:11
Double_t ATan(Double_t)
Definition: TMath.h:665
Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1)
Calculate a Breit Wigner function with mean and gamma.
Definition: TMath.cxx:437
constexpr Double_t Pi()
Definition: TMath.h:38
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12
REAL epsilon
Definition: triangle.c:617
Authors
M. Slawinska, L. Moneta

Definition in file mathmoreIntegration.C.