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.0871148
Number of function calls = 69
42201.6649413923442
Time using TF1::Integral : 0.433467
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 s1(x)
Definition: RSha256.hxx:91
@ kRed
Definition: Rtypes.h:66
@ kBlack
Definition: Rtypes.h:65
@ kBlue
Definition: Rtypes.h:66
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
Functor1D class for one-dimensional functions.
Definition: Functor.h:506
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:98
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:23
1-Dim function class
Definition: TF1.h:213
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition: TF1.cxx:2535
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:620
void SetTitle(const char *title) override
Change (i.e.
Definition: TH1.cxx:6700
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:9073
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition: TH1.cxx:3060
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:9089
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition: TH1.cxx:3107
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5025
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8857
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:274
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
@ kADAPTIVE
to be used for general functions without singularities
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)
Returns the principal value of the arc tangent of x, expressed in radians.
Definition: TMath.h:638
Double_t BreitWigner(Double_t x, Double_t mean=0, Double_t gamma=1)
Calculates a Breit Wigner function with mean and gamma.
Definition: TMath.cxx:442
constexpr Double_t Pi()
Definition: TMath.h:37
TLine l
Definition: textangle.C:4
TArc a
Definition: textangle.C:12
double epsilon
Definition: triangle.c:618
Authors
M. Slawinska, L. Moneta

Definition in file mathmoreIntegration.C.