{ "cells": [ { "cell_type": "markdown", "id": "ebbd1644", "metadata": {}, "source": [ "# mathmoreIntegration\n", "Example on the usage of the adaptive 1D integration algorithm of MathMore.\n", "\n", "It calculates the numerically cumulative integral of a distribution (like in this case the BreitWigner)\n", "to execute the macro type it (you need to compile with AClic)\n", "\n", "```cpp\n", "root[0] .x mathmoreIntegration.C+\n", "```\n", "\n", "This tutorial requires having libMathMore built with ROOT.\n", "\n", "To build mathmore you need to have a version of GSL >= 1.8 installed in your system\n", "The ROOT configure will automatically find GSL if the script gsl-config (from GSL) is in your PATH,.\n", "otherwise you need to configure root with the options --gsl-incdir and --gsl-libdir.\n", "\n", "\n", "\n", "\n", "**Author:** M. Slawinska, L. Moneta \n", "This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, May 19, 2026 at 08:25 PM." ] }, { "cell_type": "code", "execution_count": 1, "id": "e5c21a33", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:25:59.325890Z", "iopub.status.busy": "2026-05-19T20:25:59.325777Z", "iopub.status.idle": "2026-05-19T20:25:59.337398Z", "shell.execute_reply": "2026-05-19T20:25:59.336854Z" } }, "outputs": [], "source": [ "%%cpp -d\n", "#include \"TMath.h\"\n", "#include \"TH1.h\"\n", "#include \"TCanvas.h\"\n", "#include \"TLegend.h\"\n", "\n", "/*#include \"TLabel.h\"*/\n", "#include \"Math/Functor.h\"\n", "#include \"Math/WrappedFunction.h\"\n", "#include \"Math/IFunction.h\"\n", "#include \"Math/Integrator.h\"\n", "#include \n", "\n", "#include \"TStopwatch.h\"\n", "#include \"TF1.h\"\n", "\n", "#include \n", "\n", "\n", "int nc = 0;" ] }, { "cell_type": "markdown", "id": "144bb44a", "metadata": {}, "source": [ " !calculates exact integral of Breit Wigner distribution\n", "!and compares with existing methods\n", " " ] }, { "cell_type": "code", "execution_count": 2, "id": "975e76fc", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:25:59.338766Z", "iopub.status.busy": "2026-05-19T20:25:59.338650Z", "iopub.status.idle": "2026-05-19T20:25:59.341389Z", "shell.execute_reply": "2026-05-19T20:25:59.340909Z" } }, "outputs": [], "source": [ "%%cpp -d\n", "double exactIntegral( double a, double b) {\n", "\n", " return (TMath::ATan(2*b)- TMath::ATan(2*a))/ TMath::Pi();\n", "}" ] }, { "cell_type": "markdown", "id": "8171a38a", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": 3, "id": "af70b022", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:25:59.342686Z", "iopub.status.busy": "2026-05-19T20:25:59.342556Z", "iopub.status.idle": "2026-05-19T20:25:59.344838Z", "shell.execute_reply": "2026-05-19T20:25:59.344354Z" } }, "outputs": [], "source": [ "%%cpp -d\n", "double func( double x){\n", " nc++;\n", " return TMath::BreitWigner(x);\n", "}" ] }, { "cell_type": "markdown", "id": "f6d98a08", "metadata": {}, "source": [ " TF1 requires the function to have the ( )( double *, double *) signature\n", " " ] }, { "cell_type": "code", "execution_count": 4, "id": "d5a2cafb", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:25:59.346055Z", "iopub.status.busy": "2026-05-19T20:25:59.345945Z", "iopub.status.idle": "2026-05-19T20:25:59.348142Z", "shell.execute_reply": "2026-05-19T20:25:59.347689Z" } }, "outputs": [], "source": [ "%%cpp -d\n", "double func2(const double *x, const double * = nullptr){\n", " nc++;\n", " return TMath::BreitWigner(x[0]);\n", "}" ] }, { "cell_type": "markdown", "id": "13965e27", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": 5, "id": "4bbbea24", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:25:59.349413Z", "iopub.status.busy": "2026-05-19T20:25:59.349301Z", "iopub.status.idle": "2026-05-19T20:25:59.370725Z", "shell.execute_reply": "2026-05-19T20:25:59.370136Z" } }, "outputs": [], "source": [ "%%cpp -d\n", "void testIntegPerf(double x1, double x2, int n = 100000){\n", "\n", "\n", " std::cout << \"\\n\\n***************************************************************\\n\";\n", " std::cout << \"Test integration performances in interval [ \" << x1 << \" , \" << x2 << \" ]\\n\\n\";\n", "\n", " TStopwatch timer;\n", "\n", " double dx = (x2-x1)/double(n);\n", "\n", " //ROOT::Math::Functor1D f1(& TMath::BreitWigner);\n", " ROOT::Math::WrappedFunction<> f1(func);\n", "\n", " timer.Start();\n", " ROOT::Math::Integrator ig(f1 );\n", " double s1 = 0.0;\n", " nc = 0;\n", " for (int i = 0; i < n; ++i) {\n", " double x = x1 + dx*i;\n", " s1+= ig.Integral(x1,x);\n", " }\n", " timer.Stop();\n", " std::cout << \"Time using ROOT::Math::Integrator :\\t\" << timer.RealTime() << std::endl;\n", " std::cout << \"Number of function calls = \" << nc/n << std::endl;\n", " int pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);\n", "\n", "\n", "\n", " //TF1 *fBW = new TF1(\"fBW\",\"TMath::BreitWigner(x)\",x1, x2); // this is faster but cannot measure number of function calls\n", " TF1 *fBW = new TF1(\"fBW\",func2,x1, x2,0);\n", "\n", " timer.Start();\n", " nc = 0;\n", " double s2 = 0;\n", " for (int i = 0; i < n; ++i) {\n", " double x = x1 + dx*i;\n", " s2+= fBW->Integral(x1,x );\n", " }\n", " timer.Stop();\n", " std::cout << \"Time using TF1::Integral :\\t\\t\\t\" << timer.RealTime() << std::endl;\n", " std::cout << \"Number of function calls = \" << nc/n << std::endl;\n", " pr = std::cout.precision(18); std::cout << s1 << std::endl; std::cout.precision(pr);\n", "\n", "\n", "}" ] }, { "cell_type": "markdown", "id": "ab974006", "metadata": {}, "source": [ " Definition of a helper function: " ] }, { "cell_type": "code", "execution_count": 6, "id": "bbc10613", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:25:59.372155Z", "iopub.status.busy": "2026-05-19T20:25:59.372032Z", "iopub.status.idle": "2026-05-19T20:25:59.405905Z", "shell.execute_reply": "2026-05-19T20:25:59.405299Z" } }, "outputs": [], "source": [ "%%cpp -d\n", "void DrawCumulative(double x1, double x2, int n = 100){\n", "\n", " std::cout << \"\\n\\n***************************************************************\\n\";\n", " std::cout << \"Drawing cumulatives of BreitWigner in interval [ \" << x1 << \" , \" << x2 << \" ]\\n\\n\";\n", "\n", "\n", " double dx = (x2-x1)/double(n);\n", "\n", " TH1D *cum0 = new TH1D(\"cum0\", \"\", n, x1, x2); //exact cumulative\n", " for (int i = 1; i <= n; ++i) {\n", " double x = x1 + dx*i;\n", " cum0->SetBinContent(i, exactIntegral(x1, x));\n", "\n", " }\n", "\n", " // alternative method using ROOT::Math::Functor class\n", " ROOT::Math::Functor1D f1(& func);\n", "\n", "\n", " ROOT::Math::Integrator ig(f1, ROOT::Math::IntegrationOneDim::kADAPTIVE,1.E-12,1.E-12);\n", "\n", " TH1D *cum1 = new TH1D(\"cum1\", \"\", n, x1, x2);\n", " for (int i = 1; i <= n; ++i) {\n", " double x = x1 + dx*i;\n", " cum1->SetBinContent(i, ig.Integral(x1,x));\n", " }\n", "\n", "\n", " TF1 *fBW = new TF1(\"fBW\",\"TMath::BreitWigner(x, 0, 1)\",x1, x2);\n", "\n", "\n", " TH1D *cum2 = new TH1D(\"cum2\", \"\", n, x1, x2);\n", " for (int i = 1; i <= n; ++i) {\n", " double x = x1 + dx*i;\n", " cum2->SetBinContent(i, fBW->Integral(x1,x));\n", " }\n", "\n", " TH1D *cum10 = new TH1D(\"cum10\", \"\", n, x1, x2); //difference between 1 and exact\n", " TH1D *cum20 = new TH1D(\"cum23\", \"\", n, x1, x2); //difference between 2 and exact\n", " for (int i = 1; i <= n; ++i) {\n", " double delta = cum1->GetBinContent(i) - cum0->GetBinContent(i);\n", " double delta2 = cum2->GetBinContent(i) - cum0->GetBinContent(i);\n", " //std::cout << \" diff for \" << x << \" is \" << delta << \" \" << cum1->GetBinContent(i) << std::endl;\n", " cum10->SetBinContent(i, delta );\n", " cum10->SetBinError(i, std::numeric_limits::epsilon() * cum1->GetBinContent(i) );\n", " cum20->SetBinContent(i, delta2 );\n", " }\n", "\n", "\n", " TCanvas *c1 = new TCanvas(\"c1\",\"Integration example\",20,10,800,500);\n", " c1->Divide(2,1);\n", " c1->Draw();\n", "\n", " cum0->SetLineColor(kBlack);\n", " cum0->SetTitle(\"BreitWigner - the cumulative\");\n", " cum0->SetStats(false);\n", " cum1->SetLineStyle(kDashed);\n", " cum2->SetLineStyle(kDotted);\n", " cum1->SetLineColor(kBlue);\n", " cum2->SetLineColor(kRed);\n", " c1->cd(1);\n", " cum0->DrawCopy(\"h\");\n", " cum1->DrawCopy(\"same\");\n", " //cum2->DrawCopy(\"same\");\n", " cum2->DrawCopy(\"same\");\n", "\n", " c1->cd(2);\n", " cum10->SetTitle(\"Difference\");\n", " cum10->SetStats(false);\n", " cum10->SetLineColor(kBlue);\n", " cum10->Draw(\"e0\");\n", " cum20->SetLineColor(kRed);\n", " cum20->Draw(\"hsame\");\n", "\n", " TLegend * l = new TLegend(0.11, 0.8, 0.7 ,0.89);\n", " l->AddEntry(cum10, \"GSL integration - analytical \");\n", " l->AddEntry(cum20, \"TF1::Integral - analytical \");\n", " l->Draw();\n", "\n", "\n", " c1->Update();\n", " std::cout << \"\\n***************************************************************\\n\";\n", "\n", "\n", "}" ] }, { "cell_type": "markdown", "id": "b26f618d", "metadata": {}, "source": [ " Arguments are defined. " ] }, { "cell_type": "code", "execution_count": 7, "id": "e9838b3f", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:25:59.407484Z", "iopub.status.busy": "2026-05-19T20:25:59.407360Z", "iopub.status.idle": "2026-05-19T20:25:59.725887Z", "shell.execute_reply": "2026-05-19T20:25:59.725110Z" } }, "outputs": [], "source": [ "double a = -2;\n", "double b = 2;" ] }, { "cell_type": "code", "execution_count": 8, "id": "42b6ac34", "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-05-19T20:25:59.727698Z", "iopub.status.busy": "2026-05-19T20:25:59.727561Z", "iopub.status.idle": "2026-05-19T20:26:00.199743Z", "shell.execute_reply": "2026-05-19T20:26:00.199038Z" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "
\n", "
\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", "***************************************************************\n", "Drawing cumulatives of BreitWigner in interval [ -2 , 2 ]\n", "\n", "\n", "***************************************************************\n", "\n", "\n", "***************************************************************\n", "Test integration performances in interval [ -2 , 2 ]\n", "\n", "Time using ROOT::Math::Integrator :\t0.0467529\n", "Number of function calls = 69\n", "42201.6649413923442\n", "Time using TF1::Integral :\t\t\t0.143414\n", "Number of function calls = 91\n", "42201.6649413923442\n" ] } ], "source": [ "DrawCumulative(a, b);\n", "testIntegPerf(a, b);" ] } ], "metadata": { "kernelspec": { "display_name": "ROOT C++", "language": "c++", "name": "root" }, "language_info": { "codemirror_mode": "text/x-c++src", "file_extension": ".C", "mimetype": " text/x-c++src", "name": "c++" } }, "nbformat": 4, "nbformat_minor": 5 }