ROOT  6.06/09
Reference Guide
testAnalyticalIntegrals.cxx
Go to the documentation of this file.
1 //
2 // TestAnalyticalIntegrals
3 //
4 //
5 // Created by Aurélie Flandi on 10.09.14.
6 //
7 //
8 
9 #include <TStyle.h>
10 #include <TPaveText.h>
11 #include <TCanvas.h>
12 #include <iostream>//for cout
13 #include <TROOT.h>
14 #include <TFile.h>
15 #include "TMath.h"
16 #include <TF1.h>
17 #include <TLegend.h>
18 #include <TStopwatch.h>
19 #include <TApplication.h>
20 #include <Math/PdfFuncMathCore.h>//for pdf
21 #include <Math/ProbFuncMathCore.h>//for cdf
22 #include <Math/WrappedTF1.h>
23 #include <Math/Integrator.h>
24 
25 
26 using namespace std;
27 
28 void testAnalyticalIntegrals()
29 {
30 
31 //compare analytical integration with numerical one
32 
33  ROOT::Math::Integrator ig;
34 
35  TCanvas * c1 = new TCanvas("pol3","pol3",800,1000);
36  c1->Divide(3,3);
37  int ipad = 0;
38  {
39  TF1 *f = new TF1("TRY","pol3",-5.,5.);
40  f -> SetParameters(1.,1.,1.,1.);
41  c1->cd(++ipad);
42  f->Draw();
43 
44  ROOT::Math::WrappedTF1 wf(*f);
45  ig.SetFunction(wf);
46 
47  double numInt = ig.Integral(-5.,5.);
48  double anaInt = f->Integral(-5.,5.);
49 
50  std::cout<<"analytical integral for " << f->GetTitle() << " = " << anaInt << std::endl;
51  std::cout<<"numerical integral for " << f->GetTitle() << " = " << numInt << std::endl;
52 
53  if (!TMath::AreEqualAbs(numInt, anaInt, 1.E-8))
54  Error("TestAnalyticalIntegral","Different integral value for %s num = %f ana = %f diff = %f",f->GetTitle(),numInt,anaInt,numInt-anaInt);
55 
56  }
57 
58  {
59  TF1 *f = new TF1("MyExp","expo",-5.,5.);
60  f -> SetParameters(0.2,-0.3);
61  c1->cd(++ipad);
62  f ->Draw();
63 
64  ROOT::Math::WrappedTF1 wf(*f);
65 
66  double anaInt = f->Integral(-5.,5.);
67  double numInt = ig.Integral(wf,-5.,5.);
68 
69  std::cout<<"analytical integral for " << f->GetTitle() << " = " << anaInt << std::endl;
70  std::cout<<"numerical integral for " << f->GetTitle() << " = " << numInt << std::endl;
71 
72  if (!TMath::AreEqualAbs(numInt, anaInt, 1.E-8))
73  Error("TestAnalyticalIntegral","Different integral value for %s num = %f ana = %f diff = %f",f->GetTitle(),numInt,anaInt,numInt-anaInt);
74 
75  }
76  {
77  TF1 *f = new TF1("MyCrystalBall","crystalball",-5.,5.);
78  f -> SetParameters(2,1,0.5,2.,0.9);
79  c1->cd(++ipad);
80  f ->Draw();
81 
82  ROOT::Math::WrappedTF1 wf(*f);
83 
84  double anaInt = f->Integral(-5.,5.);
85  double numInt = ig.Integral(wf,-5.,5.);
86 
87  std::cout<<"analytical integral for " << f->GetTitle() << " = " << anaInt << std::endl;
88  std::cout<<"numerical integral for " << f->GetTitle() << " = " << numInt << std::endl;
89 
90  if (!TMath::AreEqualAbs(numInt, anaInt, 1.E-8))
91  Error("TestAnalyticalIntegral","Different integral value for %s num = %f ana = %f diff = %f",f->GetTitle(),numInt,anaInt,numInt-anaInt);
92 
93  }
94  {
95  // CB with alpha < 0
96  TF1 *f = new TF1("MyCrystalBall","crystalball",-5.,5.);
97  f -> SetParameters(2,-1,0.5,-2.,0.9);
98  c1->cd(++ipad);
99  f ->Draw();
100 
101  ROOT::Math::WrappedTF1 wf(*f);
102 
103  double anaInt = f->Integral(-5.,5.);
104  double numInt = ig.Integral(wf,-5.,5.);
105 
106  std::cout<<"analytical integral for " << f->GetTitle() << " = " << anaInt << std::endl;
107  std::cout<<"numerical integral for " << f->GetTitle() << " = " << numInt << std::endl;
108 
109  if (!TMath::AreEqualAbs(numInt, anaInt, 1.E-8))
110  Error("TestAnalyticalIntegral","Different integral value for %s num = %f ana = %f diff = %f",f->GetTitle(),numInt,anaInt,numInt-anaInt);
111  }
112  {
113  TF1 *f = new TF1("MyGauss","gaus",-5.,5.);
114  f -> SetParameters(2.,0.,0.3);
115  c1->cd(++ipad);
116  f ->Draw();
117  ROOT::Math::WrappedTF1 wf(*f);
118 
119  double anaInt = f->Integral(-5.,5.);
120  double numInt = ig.Integral(wf,-5.,5.);
121 
122  std::cout<<"analytical integral for " << f->GetTitle() << " = " << anaInt << std::endl;
123  std::cout<<"numerical integral for " << f->GetTitle() << " = " << numInt << std::endl;
124 
125  if (!TMath::AreEqualAbs(numInt, anaInt, 1.E-8))
126  Error("TestAnalyticalIntegral","Different integral value for %s num = %f ana = %f diff = %f",f->GetTitle(),numInt,anaInt,numInt-anaInt);
127 }
128  {
129  TF1 *f = new TF1("MyGauss","gausn",-5.,5.);
130  f -> SetParameters(2.,0.,0.3);
131  c1->cd(++ipad);
132  f ->Draw();
133 
134  ROOT::Math::WrappedTF1 wf(*f);
135 
136  double anaInt = f->Integral(-5.,5.);
137  double numInt = ig.Integral(wf,-5.,5.);
138 
139  std::cout<<"analytical integral for " << f->GetTitle() << " = " << anaInt << std::endl;
140  std::cout<<"numerical integral for " << f->GetTitle() << " = " << numInt << std::endl;
141 
142  if (!TMath::AreEqualAbs(numInt, anaInt, 1.E-6))
143  Error("TestAnalyticalIntegral","Different integral value for %s num = %20.10f ana = %20.10f diff = %20.10f",f->GetTitle(),numInt,anaInt,numInt-anaInt);
144 
145  }
146 
147  {
148  TF1 *f = new TF1("MyExp","landau",-5.,5.);
149  f -> SetParameters(2.,2,0.3);
150  c1->cd(++ipad);
151  f ->Draw();
152 
153  ROOT::Math::WrappedTF1 wf(*f);
154 
155  double anaInt = f->Integral(-5.,5.);
156  double numInt = ig.Integral(wf,-5.,5.);
157 
158  std::cout<<"analytical integral for " << f->GetTitle() << " = " << anaInt << std::endl;
159  std::cout<<"numerical integral for " << f->GetTitle() << " = " << numInt << std::endl;
160 
161  if (!TMath::AreEqualAbs(numInt, anaInt, 1.E-8))
162  Error("TestAnalyticalIntegral","Different integral value for %s num = %f ana = %f diff = %f",f->GetTitle(),numInt,anaInt,numInt-anaInt);
163 
164  }
165  {
166  TF1 *f = new TF1("MyExp","landaun",-5.,5.);
167  f -> SetParameters(2.,-2.,0.3);
168  c1->cd(++ipad);
169  f ->Draw();
170  ROOT::Math::WrappedTF1 wf(*f);
171 
172  double anaInt = f->Integral(-5.,5.);
173  double numInt = ig.Integral(wf,-5.,5.);
174 
175  std::cout<<"analytical integral for " << f->GetTitle() << " = " << anaInt << std::endl;
176  std::cout<<"numerical integral for " << f->GetTitle() << " = " << numInt << std::endl;
177 
178  if (!TMath::AreEqualAbs(numInt, anaInt, 1.E-8))
179  Error("TestAnalyticalIntegral","Different integral value for %s num = %f ana = %f diff = %f",f->GetTitle(),numInt,anaInt,numInt-anaInt);
180 
181  }
182 }
183 
184 int main(int argc, char **argv)
185 {
186  bool showGraphics = false;
187  // Parse command line arguments
188  for (Int_t i=1 ; i<argc ; i++) {
189  std::string arg = argv[i] ;
190  if (arg == "-g") {
191  showGraphics = true;
192  }
193  if (arg == "-h") {
194  cerr << "Usage: " << argv[0] << " [-g] [-v]\n";
195  cerr << " where:\n";
196  cerr << " -g : graphics mode\n";
197  cerr << " -v : verbose mode";
198  cerr << endl;
199  return -1;
200  }
201  }
202 
203  TApplication* theApp = 0;
204  if ( showGraphics )
205  theApp = new TApplication("App",&argc,argv);
206 
207  testAnalyticalIntegrals();
208 
209  if ( showGraphics )
210  {
211  theApp->Run();
212  delete theApp;
213  theApp = 0;
214  }
215 
216  return 0;
217 }
TNamed::GetTitle
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52
c1
TCanvas * c1
Definition: legend1.C:2
TCanvas::cd
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:659
ROOT::Math::WrappedTF1
Class to Wrap a ROOT Function class (like TF1) in a IParamFunction interface of one dimensions to be ...
Definition: WrappedTF1.h:41
Int_t
int Int_t
Definition: RtypesCore.h:41
TF1::Draw
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1075
std
STL namespace.
TF1::Integral
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition: TF1.cxx:2264
SetParameters
void SetParameters(TFitEditor::FuncParams_t &pars, TF1 *func)
Restore the parameters from pars into the function.
Definition: TFitEditor.cxx:287
TApplication::Run
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
Definition: TApplication.cxx:1139
main
int main(int argc, char **argv)
Definition: testAnalyticalIntegrals.cxx:184
Error
void Error(const char *location, const char *msgfmt,...)
TMath::AreEqualAbs
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition: TMath.h:192
ROOT::Math::IntegratorOneDim::Integral
double Integral(Function &f, double a, double b)
evaluate the Integral of a function f over the defined interval (a,b)
Definition: Integrator.h:506
ROOT::Math::IntegratorOneDim
User Class for performing numerical integration of a function in one dimension.
Definition: Integrator.h:102
TMath::E
Double_t E()
Definition: TMath.h:54
TCanvas
The Canvas class.
Definition: TCanvas.h:48
f
double f(double x)
Definition: testIntegration.cxx:12
testAnalyticalIntegrals
void testAnalyticalIntegrals()
Definition: testAnalyticalIntegrals.cxx:28
ROOT::Math::IntegratorOneDim::SetFunction
void SetFunction(Function &f)
method to set the a generic integration function
Definition: Integrator.h:499
TPad::Divide
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1077
TF1
1-Dim function class
Definition: TF1.h:149
TApplication
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:45
showGraphics
bool showGraphics
Definition: GAMinTutorial.cxx:42