ROOT   Reference Guide
mathmoreIntegration.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_math
3/// \notebook -nodraw
4/// Example on the usage of the adaptive 1D integration algorithm of MathMore
5/// it calculates the numerically cumulative integral of a distribution (like in this case the BreitWigner)
6/// to execute the macro type it (you need to compile with AClic)
7///
8/// ~~~{.cpp}
9/// root[0] .x mathmoreIntegration.C+
10/// ~~~
11///
12/// This tutorials require having libMathMore built with ROOT.
13///
14/// To build mathmore you need to have a version of GSL >= 1.8 installed in your system
15/// The ROOT configure will automatically find GSL if the script gsl-config (from GSL) is in your PATH,.
16/// otherwise you need to configure root with the options --gsl-incdir and --gsl-libdir.
17///
18/// \macro_image
19/// \macro_output
20/// \macro_code
21///
22/// \authors M. Slawinska, L. Moneta
23
24#include "TMath.h"
25#include "TH1.h"
26#include "TCanvas.h"
27#include "TLegend.h"
28
29/*#include "TLabel.h"*/
30#include "Math/Functor.h"
32#include "Math/IFunction.h"
33#include "Math/Integrator.h"
34#include <iostream>
35
36#include "TStopwatch.h"
37#include "TF1.h"
38
39#include <limits>
40
41//!calculates exact integral of Breit Wigner distribution
42//!and compares with existing methods
43
44int nc = 0;
45double exactIntegral( double a, double b) {
46
47 return (TMath::ATan(2*b)- TMath::ATan(2*a))/ TMath::Pi();
48}
49double func( double x){
50 nc++;
51 return TMath::BreitWigner(x);
52}
53
54// TF1 requires the function to have the ( )( double *, double *) signature
55double func2(const double *x, const double * = 0){
56 nc++;
57 return TMath::BreitWigner(x[0]);
58}
59
60
61
62
63void testIntegPerf(double x1, double x2, int n = 100000){
64
65
66 std::cout << "\n\n***************************************************************\n";
67 std::cout << "Test integration performances in interval [ " << x1 << " , " << x2 << " ]\n\n";
68
69 TStopwatch timer;
70
71 double dx = (x2-x1)/double(n);
72
73 //ROOT::Math::Functor1D<ROOT::Math::IGenFunction> f1(& TMath::BreitWigner);
75
76 timer.Start();
78 double s1 = 0.0;
79 nc = 0;
80 for (int i = 0; i < n; ++i) {
81 double x = x1 + dx*i;
82 s1+= ig.Integral(x1,x);
83 }
84 timer.Stop();
85 std::cout << "Time using ROOT::Math::Integrator :\t" << timer.RealTime() << std::endl;
86 std::cout << "Number of function calls = " << nc/n << std::endl;
87 int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
88
89
90
91 //TF1 *fBW = new TF1("fBW","TMath::BreitWigner(x)",x1, x2); // this is faster but cannot measure number of function calls
92 TF1 *fBW = new TF1("fBW",func2,x1, x2,0);
93
94 timer.Start();
95 nc = 0;
96 double s2 = 0;
97 for (int i = 0; i < n; ++i) {
98 double x = x1 + dx*i;
99 s2+= fBW->Integral(x1,x );
100 }
101 timer.Stop();
102 std::cout << "Time using TF1::Integral :\t\t\t" << timer.RealTime() << std::endl;
103 std::cout << "Number of function calls = " << nc/n << std::endl;
104 pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);
105
106
107}
108
109void DrawCumulative(double x1, double x2, int n = 100){
110
111 std::cout << "\n\n***************************************************************\n";
112 std::cout << "Drawing cumulatives of BreitWigner in interval [ " << x1 << " , " << x2 << " ]\n\n";
113
114
115 double dx = (x2-x1)/double(n);
116
117 TH1D *cum0 = new TH1D("cum0", "", n, x1, x2); //exact cumulative
118 for (int i = 1; i <= n; ++i) {
119 double x = x1 + dx*i;
120 cum0->SetBinContent(i, exactIntegral(x1, x));
121
122 }
123
124 // alternative method using ROOT::Math::Functor class
126
127
129
130 TH1D *cum1 = new TH1D("cum1", "", n, x1, x2);
131 for (int i = 1; i <= n; ++i) {
132 double x = x1 + dx*i;
133 cum1->SetBinContent(i, ig.Integral(x1,x));
134 }
135
136
137 TF1 *fBW = new TF1("fBW","TMath::BreitWigner(x, 0, 1)",x1, x2);
138
139
140 TH1D *cum2 = new TH1D("cum2", "", n, x1, x2);
141 for (int i = 1; i <= n; ++i) {
142 double x = x1 + dx*i;
143 cum2->SetBinContent(i, fBW->Integral(x1,x));
144 }
145
146 TH1D *cum10 = new TH1D("cum10", "", n, x1, x2); //difference between 1 and exact
147 TH1D *cum20 = new TH1D("cum23", "", n, x1, x2); //difference between 2 and excact
148 for (int i = 1; i <= n; ++i) {
149 double delta = cum1->GetBinContent(i) - cum0->GetBinContent(i);
150 double delta2 = cum2->GetBinContent(i) - cum0->GetBinContent(i);
151 //std::cout << " diff for " << x << " is " << delta << " " << cum1->GetBinContent(i) << std::endl;
152 cum10->SetBinContent(i, delta );
154 cum20->SetBinContent(i, delta2 );
155 }
156
157
158 TCanvas *c1 = new TCanvas("c1","Integration example",20,10,800,500);
159 c1->Divide(2,1);
160 c1->Draw();
161
162 cum0->SetLineColor(kBlack);
163 cum0->SetTitle("BreitWigner - the cumulative");
164 cum0->SetStats(0);
165 cum1->SetLineStyle(2);
166 cum2->SetLineStyle(3);
167 cum1->SetLineColor(kBlue);
168 cum2->SetLineColor(kRed);
169 c1->cd(1);
170 cum0->DrawCopy("h");
171 cum1->DrawCopy("same");
172 //cum2->DrawCopy("same");
173 cum2->DrawCopy("same");
174
175 c1->cd(2);
176 cum10->SetTitle("Difference");
177 cum10->SetStats(0);
178 cum10->SetLineColor(kBlue);
179 cum10->Draw("e0");
180 cum20->SetLineColor(kRed);
181 cum20->Draw("hsame");
182
183 TLegend * l = new TLegend(0.11, 0.8, 0.7 ,0.89);
184 l->AddEntry(cum10, "GSL integration - analytical ");
185 l->AddEntry(cum20, "TF1::Integral - analytical ");
186 l->Draw();
187
188
189 c1->Update();
190 std::cout << "\n***************************************************************\n";
191
192
193}
194
195
196
197void mathmoreIntegration(double a = -2, double b = 2)
198{
199 DrawCumulative(a, b);
200 testIntegPerf(a, b);
201}
202
#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:2544
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
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