Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 tutorial requires 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.094486
Number of function calls = 69
42201.6649413923442
Time using TF1::Integral : 0.37215
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 * = nullptr){
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 exact
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 );
cum10->SetBinError(i, std::numeric_limits<double>::epsilon() * cum1->GetBinContent(i) );
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(false);
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(false);
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 a(i)
Definition RSha256.hxx:99
#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 x2
Option_t Option_t TPoint TPoint const char x1
Functor1D class for one-dimensional functions.
Definition Functor.h:95
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:233
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition TF1.cxx:2531
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:670
void SetTitle(const char *title) override
Change/set the title.
Definition TH1.cxx:6709
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:9197
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3066
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:9213
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition TH1.cxx:3113
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5052
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition TH1.cxx:8981
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...
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void Stop()
Stop the stopwatch.
@ 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:640
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
Authors
M. Slawinska, L. Moneta

Definition in file mathmoreIntegration.C.