ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
testRandomDist.cxx
Go to the documentation of this file.
1 #include "Math/Random.h"
2 #include "Math/GSLRndmEngines.h"
3 #include "Math/DistFunc.h"
4 #include "TStopwatch.h"
5 #include "TRandom3.h"
6 #include "TRandom2.h"
7 #include "TCanvas.h"
8 #include "TH1D.h"
9 #include "TH2D.h"
10 #include "TH3D.h"
11 #include "TMath.h"
12 #include <iostream>
13 #include <cmath>
14 #include <typeinfo>
15 #ifdef HAVE_UNURAN
16 #include "UnuRanDist.h"
17 #endif
18 
19 #ifdef HAVE_CLHEP
20 #include "CLHEP/Random/RandFlat.h"
21 #include "CLHEP/Random/RandPoissonT.h"
22 #include "CLHEP/Random/RandPoisson.h"
23 #include "CLHEP/Random/RandGauss.h"
24 #include "CLHEP/Random/RandGaussQ.h"
25 #include "CLHEP/Random/RandGaussT.h"
26 #include "CLHEP/Random/RandBinomial.h"
27 #include "CLHEP/Random/JamesRandom.h"
28 #endif
29 
30 
31 #ifndef PI
32 #define PI 3.14159265358979323846264338328 /* pi */
33 #endif
34 
35 
36 
37 #ifndef NEVT
38 #define NEVT 1000000
39 #endif
40 
41 //#define TEST_TIME
42 
43 using namespace ROOT::Math;
44 #ifdef HAVE_CLHEP
45 using namespace CLHEP;
46 #endif
47 
48 static bool fillHist = false;
49 
50 
51 
52 void testDiff(TH1D & h1, TH1D & h2, const std::string & name="") {
53 
54  double chi2;
55  int ndf;
56  if (h1.GetEntries() == 0 && h2.GetEntries() == 0) return;
57  int igood = 0;
58  double prob = h1.Chi2TestX(&h2,chi2,ndf,igood,"UU");
59  std::cout << "\nTest " << name << " chi2 = " << chi2 << " ndf " << ndf << " prob = " << prob;
60  if (igood) std::cout << " \t\t" << " igood = " << igood;
61  std::cout << std::endl;
62 
63  std::string cname="c1_" + name;
64  std::string ctitle="Test of " + name;
65  TCanvas *c1 = new TCanvas(cname.c_str(), ctitle.c_str(),200,10,800,600);
66  h1.DrawCopy();
67  h2.DrawCopy("Esame");
68  c1->Update();
69 
70 
71 }
72 
73 template <class R>
74 std::string findName( const R & r) {
75 
76  std::string type = typeid(r).name();
77  if (type.find("GSL") != std::string::npos )
78  return "ROOT::Math::Random";
79  else if (type.find("TRandom") != std::string::npos )
80  return "TRandom ";
81  else if (type.find("UnuRan") != std::string::npos )
82  return "UnuRan ";
83  else if (type.find("Rand") != std::string::npos )
84  return "CLHEP ";
85 
86  return "????????? ";
87 }
88 
89 
90 template <class R>
91 void testPoisson( R & r,double mu,TH1D & h) {
92 
93  TStopwatch w;
94 
95  int n = NEVT;
96  // estimate PI
97  w.Start();
98  r.SetSeed(0);
99  for (int i = 0; i < n; ++i) {
100  int n = r.Poisson(mu );
101  if (fillHist)
102  h.Fill( double(n) );
103  }
104  w.Stop();
105  if (fillHist) { fillHist=false; return; }
106  std::cout << "Poisson - mu = " << mu << "\t\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
107  << w.CpuTime()*1.0E9/NEVT
108  << "\t(ns/call)" << std::endl;
109  // fill histogram the second pass
110  fillHist = true;
111  testPoisson(r,mu,h);
112 }
113 
114 
115 
116 // Knuth Algorith for Poisson (used also in GSL)
117 template <class R>
118 unsigned int genPoisson( R & r, double mu) {
119 
120  // algorithm described by Knuth Vol 2. 2nd edition pag. 132
121  // for generating poisson deviates when mu is large
122 
123  unsigned int k = 0;
124 
125  while (mu > 20) {
126 
127  const double alpha = 0.875; // 7/8
128  unsigned int m = static_cast<unsigned int> ((alpha*mu) );
129 
130  // generate xg according to a Gamma distribution of m
131  double sqm = std::sqrt( 2.*m -1.);
132  double pi = TMath::Pi();
133  double x,y,v;
134  do {
135  do {
136  y = std::tan( pi * r.Rndm() );
137  x = sqm * y + m - 1.;
138  }
139  while (x <= 0);
140  v = r.Rndm();
141  }
142  while ( v > (1 + y * y) * std::exp ( (m - 1) * std::log (x / (m - 1)) - sqm * y));
143 
144  // x is now distributed according to a gamma of m
145 
146  if ( x >= mu )
147  return k + r.Binomial( m-1, mu/x);
148 
149  else {
150  // continue the loop decresing mu
151  mu -= x;
152  k += m;
153  }
154  }
155  // for lower values of mu use rejection method from exponential
156  double expmu = TMath::Exp(-mu);
157  double pir = 1.0;
158  do {
159  pir *= r.Rndm();
160  k++;
161  }
162  while (pir > expmu);
163  return k -1;
164 
165 }
166 
167 
168 
169 // Numerical Receip algorithm for Poisson (used also in CLHEP)
170 template <class R>
171 unsigned int genPoisson2( R & r, double mu) {
172 
173  //double om = getOldMean();
174 
175 
176 
177 
178  if( mu < 12.0 ) {
179 
180  double expmu = TMath::Exp(-mu);
181  double pir = 1.0;
182  unsigned int k = 0;
183  do {
184  pir *= r.Rndm();
185  k++;
186  }
187  while (pir > expmu);
188  return k-1;
189  }
190  // for large mu values (should care for values larger than 2E9)
191  else {
192 
193  double em, t, y;
194  double sq, alxm, g;
195  double pi = TMath::Pi();
196 
197  sq = std::sqrt(2.0*mu);
198  alxm = std::log(mu);
199  g = mu*alxm - TMath::LnGamma(mu + 1.0);
200 
201  do {
202  do {
203  y = std::tan(pi*r.Rndm());
204  em = sq*y + mu;
205  } while( em < 0.0 );
206 
207  em = std::floor(em);
208  t = 0.9*(1.0 + y*y)* std::exp(em*alxm - TMath::LnGamma(em + 1.0) - g);
209  } while( r.Rndm() > t );
210 
211  return static_cast<unsigned int> (em);
212 
213  }
214 
215 }
216 
217 
218 template <class R>
219 void testPoisson2( R & r,double mu,TH1D & h) {
220 
221  TStopwatch w;
222 
223  int n = NEVT;
224  // estimate PI
225  w.Start();
226  r.SetSeed(0);
227  for (int i = 0; i < n; ++i) {
228  // int n = genPoisson2(r,mu);
229  int n = r.PoissonD(mu);
230  if (fillHist)
231  h.Fill( double(n) );
232  }
233  w.Stop();
234  if (fillHist) { fillHist=false; return; }
235  std::cout << "Poisson \t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
236  << w.CpuTime()*1.0E9/NEVT
237  << "\t(ns/call)" << std::endl;
238  // fill histogram the second pass
239  fillHist = true;
240  testPoisson2(r,mu,h);
241 }
242 
243 #ifdef HAVE_CLHEP
244 template<class R>
245 void testPoissonCLHEP( R & r, double mu,TH1D & h) {
246 
247  TStopwatch w;
248 
249  int n = NEVT;
250  // estimate PI
251  w.Start();
252  // r.SetSeed(0);
253  for (int i = 0; i < n; ++i) {
254  //int n = RandPoisson::shoot(mu + RandFlat::shoot());
255  int n = static_cast<unsigned int> ( r(mu) ) ;
256  if (fillHist)
257  h.Fill( double(n) );
258  }
259  w.Stop();
260  if (fillHist) { fillHist=false; return; }
261  std::cout << "Poisson - mu = " << mu << "\t\t" << findName(r) <<"\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
262  << w.CpuTime()*1.0E9/NEVT
263  << "\t(ns/call)" << std::endl;
264  // fill histogram the second pass
265  fillHist = true;
266  testPoissonCLHEP(r,mu,h);
267 }
268 
269 template<class R>
270 void testGausCLHEP( R & r,double mu,double sigma,TH1D & h) {
271 
272  TStopwatch w;
273 
274  int n = NEVT;
275  int n1 = 100;
276  int n2 = n/n1;
277  w.Start();
278  for (int i = 0; i < n2; ++i) {
279  for (int j = 0; j < n1; ++j) {
280  double x = r(mu,sigma );
281  if (fillHist)
282  h.Fill( x );
283 
284  }
285  w.Stop();
286  if (fillHist) { fillHist=false; return; }
287  std::cout << "Gaussian - mu,sigma = " << mu << " , " << sigma << "\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
288  << w.CpuTime()*1.0E9/NEVT
289  << "\t(ns/call)" << std::endl;
290  // fill histogram the second pass
291  fillHist = true;
292  testGausCLHEP(r,mu,sigma,h);
293 }
294 
295 template <class R>
296 void testFlatCLHEP( R & r,TH1D & h) {
297 
298  TStopwatch w;
299 
300  int n = NEVT;
301  w.Start();
302  //r.SetSeed(0);
303  for (int i = 0; i < n; ++i) {
304  double x = r();
305  if (fillHist)
306  h.Fill( x );
307 
308  }
309  w.Stop();
310  if (fillHist) { fillHist=false; return; }
311  std::cout << "Flat - [0,1] \t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
312  << w.CpuTime()*1.0E9/NEVT
313  << "\t(ns/call)" << std::endl;
314  // fill histogram the second pass
315  fillHist = true;
316  testFlatCLHEP(r,h);
317 }
318 
319 
320 #endif
321 
322 
323 
324 template <class R>
325 void testFlat( R & r,TH1D & h) {
326 
327  TStopwatch w;
328 
329  int n = NEVT;
330  w.Start();
331  r.SetSeed(0);
332  for (int i = 0; i < n; ++i) {
333  double x = r.Rndm();
334  if (fillHist)
335  h.Fill( x );
336 
337  }
338  w.Stop();
339  if (fillHist) { fillHist=false; return; }
340  std::cout << "Flat - [0,1] \t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
341  << w.CpuTime()*1.0E9/NEVT
342  << "\t(ns/call)" << std::endl;
343  // fill histogram the second pass
344  fillHist = true;
345  testFlat(r,h);
346 }
347 
348 
349 
350 template <class R>
351 void testGaus( R & r,double mu,double sigma,TH1D & h) {
352 
353  TStopwatch w;
354 
355  int n = NEVT;
356  // estimate PI
357  w.Start();
358  r.SetSeed(0);
359  for (int i = 0; i < n; ++i) {
360  double x = r.Gaus(mu,sigma );
361  if (fillHist)
362  h.Fill( x );
363 
364  }
365  w.Stop();
366  if (fillHist) { fillHist=false; return; }
367  std::cout << "Gaussian - mu,sigma = " << mu << " , " << sigma << "\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
368  << w.CpuTime()*1.0E9/NEVT
369  << "\t(ns/call)" << std::endl;
370  // fill histogram the second pass
371  fillHist = true;
372  testGaus(r,mu,sigma,h);
373 }
374 
375 
376 
377 template <class R>
378 void testLandau( R & r,TH1D & h) {
379 
380  TStopwatch w;
381 
382  int n = NEVT;
383  // estimate PI
384  w.Start();
385  r.SetSeed(0);
386  for (int i = 0; i < n; ++i) {
387  double x = r.Landau();
388  if (fillHist)
389  h.Fill( x );
390 
391  }
392  w.Stop();
393  if (fillHist) { fillHist=false; return; }
394  std::cout << "Landau " << "\t\t\t\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
395  << w.CpuTime()*1.0E9/NEVT
396  << "\t(ns/call)" << std::endl;
397  // fill histogram the second pass
398  fillHist = true;
399  testLandau(r,h);
400 }
401 
402 
403 
404 
405 template <class R>
406 void testBreitWigner( R & r,double mu,double gamma,TH1D & h) {
407 
408  TStopwatch w;
409 
410  int n = NEVT;
411  // estimate PI
412  w.Start();
413  r.SetSeed(0);
414  for (int i = 0; i < n; ++i) {
415  double x = r.BreitWigner(mu,gamma );
416  if (fillHist)
417  h.Fill( x );
418  }
419  w.Stop();
420  if (fillHist) { fillHist=false; return; }
421  std::cout << "Breit-Wigner - m,g = " << mu << " , " << gamma << "\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
422  << w.CpuTime()*1.0E9/NEVT
423  << "\t(ns/call)" << std::endl;
424  // fill histogram the second pass
425  fillHist = true;
426  testBreitWigner(r,mu,gamma,h);
427 }
428 
429 
430 
431 template <class R>
432 void testBinomial( R & r,int ntot,double p,TH1D & h) {
433 
434  TStopwatch w;
435 
436  int n = NEVT;
437  // estimate PI
438  w.Start();
439  r.SetSeed(0);
440  for (int i = 0; i < n; ++i) {
441  double x = double( r.Binomial(ntot,p ) );
442  if (fillHist)
443  h.Fill( x );
444  }
445  w.Stop();
446  if (fillHist) { fillHist=false; return; }
447  std::cout << "Binomial - ntot,p = " << ntot << " , " << p << "\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
448  << w.CpuTime()*1.0E9/NEVT
449  << "\t(ns/call)" << std::endl;
450  // fill histogram the second pass
451  fillHist = true;
452  testBinomial(r,ntot,p,h);
453 }
454 
455 
456 template<class R>
457 void testBinomialCLHEP( R & r,int ntot,double p,TH1D & h) {
458 
459  TStopwatch w;
460 
461  int n = NEVT;
462  // estimate PI
463  w.Start();
464  //r.SetSeed(0);
465  for (int i = 0; i < n; ++i) {
466  double x = double( r(ntot,p ) );
467  if (fillHist)
468  h.Fill( x );
469  }
470  w.Stop();
471  if (fillHist) { fillHist=false; return; }
472  std::cout << "Binomial - ntot,p = " << ntot << " , " << p << "\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
473  << w.CpuTime()*1.0E9/NEVT
474  << "\t(ns/call)" << std::endl;
475  // fill histogram the second pass
476  fillHist = true;
477  testBinomialCLHEP(r,ntot,p,h);
478 }
479 
480 
481 template <class R>
482 void testMultinomial( R & r,int ntot, TH1D & h1, TH1D & h2) {
483 
484  // generates the p distribution
485  const int nbins = h1.GetNbinsX();
486  std::vector<double> p(nbins);
487  double psum = 0;
488  for (int i = 0; i < nbins; ++i) {
489  double x1 = h1.GetBinLowEdge(i+1);
490  double x2 = x1 + h1.GetBinWidth(i+1);
491  p[i] = ROOT::Math::normal_cdf(x2) - ROOT::Math::normal_cdf(x1);
492  psum += p[i];
493  }
494  //std::cout << " psum = " << psum << std::endl;
495  // generate the multinomial
496  TStopwatch w;
497  int n = NEVT/10;
498  if (fillHist) n = 1;
499 
500  for (int i = 0; i < n; ++i) {
501  std::vector<unsigned int> multDist = r.Multinomial(ntot,p);
502  if (fillHist) {
503  for (int j = 0; j < nbins; ++j) h1.SetBinContent(j+1,multDist[j]);
504  }
505  }
506  w.Stop();
507 
508  if (!fillHist) {
509  std::cout << "Multinomial - nb, ntot = " << nbins << " , " << ntot << "\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/n << " \t"
510  << w.CpuTime()*1.0E9/n
511  << "\t(ns/call)" << std::endl;
512  }
513 
514  // check now time using poisson distribution
515 
516  w.Start();
517  for (int i = 0; i < n; ++i) {
518  for (int j = 0; j < nbins; ++j) {
519  double y = r.Poisson(ntot*p[j]);
520  if (fillHist) h2.SetBinContent(j+1,y);
521  }
522  }
523  w.Stop();
524  if (!fillHist) {
525  std::cout << "Multi-Poisson - nb, ntot = " << nbins << " , " << ntot << "\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/n << " \t"
526  << w.CpuTime()*1.0E9/n
527  << "\t(ns/call)" << std::endl;
528  }
529 
530  if (fillHist) { fillHist=false; return; }
531 
532  // fill histogram the second pass
533  fillHist = true;
534  testMultinomial(r,ntot,h1,h2);
535 
536 }
537 
538 
539 template <class R>
540 void testExp( R & r,TH1D & h) {
541 
542  TStopwatch w;
543 
544  int n = NEVT;
545  // estimate PI
546  w.Start();
547  r.SetSeed(0);
548  for (int i = 0; i < n; ++i) {
549  double x = r.Exp(1.);
550  if (fillHist)
551  h.Fill( x );
552  }
553  w.Stop();
554  if (fillHist) { fillHist=false; return; }
555  std::cout << "Exponential " << "\t\t\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
556  << w.CpuTime()*1.0E9/NEVT
557  << "\t(ns/call)" << std::endl;
558  // fill histogram the second pass
559  fillHist = true;
560  testExp(r,h);
561 }
562 
563 
564 template <class R>
565 void testCircle( R & r,TH1D & h) {
566 
567  TStopwatch w;
568 
569  int n = NEVT;
570  // estimate PI
571  w.Start();
572  r.SetSeed(0);
573  double x,y;
574  for (int i = 0; i < n; ++i) {
575  r.Circle(x,y,1.0);
576  if (fillHist)
577  h.Fill( std::atan2(x,y) );
578 
579  }
580  w.Stop();
581  if (fillHist) { fillHist=false; return; }
582 
583  std::cout << "Circle " << "\t\t\t\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
584  << w.CpuTime()*1.0E9/NEVT
585  << "\t(ns/call)" << std::endl;
586  // fill histogram the second pass
587  fillHist = true;
588  testCircle(r,h);
589 
590 }
591 
592 
593 template <class R>
594 void testSphere( R & r,TH1D & h1, TH1D & h2 ) {
595 
596 
597 #ifdef PLOT_SPHERE
598  TH2D hxy("hxy","xy",100,-1.1,1.1,100,-1.1,1.1);
599  TH3D h3d("h3d","sphere",100,-1.1,1.1,100,-1.1,1.1,100,-1.1,1.1);
600  TH1D hz("hz","z",100,-1.1,1.1);
601 #endif
602 
603  TStopwatch w;
604 
605 
606  int n = NEVT;
607  // estimate PI
608  w.Start();
609  r.SetSeed(0);
610  double x,y,z;
611  for (int i = 0; i < n; ++i) {
612  r.Sphere(x,y,z,1.0);
613  if (fillHist) {
614 
615  h1.Fill( std::atan2(x,y) );
616  h2.Fill( std::atan2( std::sqrt(x*x+y*y), z ) );
617 
618 #ifdef PLOT_SPHERE
619  hxy.Fill(x,y);
620  hz.Fill(z);
621  h3d.Fill(x,y,z);
622 #endif
623 
624  }
625 
626  }
627 
628  w.Stop();
629 
630 #ifdef PLOT_SPHERE
631  if (fillHist) {
632  TCanvas *c1 = new TCanvas("c1_xyz","sphere",220,20,800,900);
633  c1->Divide(2,2);
634  c1->cd(1);
635  hxy.DrawCopy();
636  c1->cd(2);
637  hz.DrawCopy();
638  c1->cd(3);
639  h3d.DrawCopy();
640  c1->Update();
641  }
642 #endif
643 
644  if (fillHist) { fillHist=false; return; }
645 
646  std::cout << "Sphere " << "\t\t\t\t"<< findName(r) << "\tTime = " << w.RealTime()*1.0E9/NEVT << " \t"
647  << w.CpuTime()*1.0E9/NEVT
648  << "\t(ns/call)" << std::endl;
649 
650  // fill histogram the second pass
651  fillHist = true;
652  testSphere(r,h1,h2);
653 
654 }
655 
656 
657 int testRandomDist() {
658 
659  std::cout << "***************************************************\n";
660  std::cout << " TEST RANDOM DISTRIBUTIONS NEVT = " << NEVT << std::endl;
661  std::cout << "***************************************************\n\n";
662 
663 
664 
665  Random<GSLRngMT> r;
666  TRandom3 tr;
667 #ifdef HAVE_UNURAN
668  UnuRanDist ur;
669 #else
670  TRandom2 ur;
671 #endif
672 
673  // flat
674  double xmin = 0;
675  double xmax = 1;
676  int nch = 1000;
677  TH1D hf1("hf1","FLAT ROOT",nch,xmin,xmax);
678  TH1D hf2("hf2","Flat GSL",nch,xmin,xmax);
679 
680  testFlat(r,hf1);
681  testFlat(tr,hf2);
682  testDiff(hf1,hf2,"Flat ROOT-GSL");
683 
684 #ifdef HAVE_CLHEP
685  HepJamesRandom eng;
686  RandFlat crf(eng);
687  TH1D hf3("hf3","Flat CLHEP",nch,xmin,xmax);
688  testFlatCLHEP(crf,hf3);
689  testDiff(hf3,hf1,"Flat CLHEP-GSL");
690 #endif
691 
692 
693 
694  // Poisson
695  std::cout << std::endl;
696 
697  double mu = 25;
698  xmin = std::floor(std::max(0.0,mu-5*std::sqrt(mu) ) );
699  xmax = std::floor( mu+5*std::sqrt(mu) );
700  nch = std::min( int(xmax-xmin),1000);
701  TH1D hp1("hp1","Poisson ROOT",nch,xmin,xmax);
702  TH1D hp2("hp2","Poisson GSL",nch,xmin,xmax);
703  TH1D hp3("hp3","Poisson UNR",nch,xmin,xmax);
704 
705  testPoisson(r,mu,hp1);
706  testPoisson(tr,mu,hp2);
707  testPoisson(ur,mu,hp3);
708 #ifdef HAVE_CLHEP
709  RandPoissonT crp(eng);
710  TH1D hp4("hp4","Poisson CLHEP",nch,xmin,xmax);
711  testPoissonCLHEP(crp,mu,hp4);
712 #endif
713  //testPoisson2(tr,mu,h2);
714  // test differences
715  testDiff(hp1,hp2,"Poisson ROOT-GSL");
716  testDiff(hp1,hp3,"Poisson ROOT-UNR");
717 #ifdef HAVE_CLHEP
718  testDiff(hp1,hp4,"Poisson ROOT-CLHEP");
719 #endif
720 
721  // Gaussian
722  std::cout << std::endl;
723 
724  TH1D hg1("hg1","Gaussian ROOT",nch,xmin,xmax);
725  TH1D hg2("hg2","Gaussian GSL",nch,xmin,xmax);
726  TH1D hg3("hg3","Gaussian UNURAN",nch,xmin,xmax);
727 
728 
729  testGaus(r,mu,sqrt(mu),hg1);
730  testGaus(tr,mu,sqrt(mu),hg2);
731 
732  testGaus(ur,mu,sqrt(mu),hg3);
733 #ifdef HAVE_CLHEP
734  RandGauss crg(eng);
735  TH1D hg4("hg4","Gauss CLHEP",nch,xmin,xmax);
736  testGausCLHEP(crg,mu,sqrt(mu),hg4);
737  RandGaussQ crgQ(eng);
738  testGausCLHEP(crgQ,mu,sqrt(mu),hg4);
739  RandGaussT crgT(eng);
740  testGausCLHEP(crgT,mu,sqrt(mu),hg4);
741 #endif
742 
743 
744  testDiff(hg1,hg2,"Gaussian ROOT-GSL");
745  testDiff(hg1,hg3,"Gaussian ROOT_UNR");
746 
747  // Landau
748  std::cout << std::endl;
749 
750  TH1D hl1("hl1","Landau ROOT",300,0,50);
751  TH1D hl2("hl2","Landau GSL",300,0,50);
752 
753  testLandau(r,hl1);
754  testLandau(tr,hl2);
755  testDiff(hl1,hl2,"Landau");
756 
757  // Breit Wigner
758  std::cout << std::endl;
759 
760  TH1D hbw1("hbw1","BreitWigner ROOT",nch,xmin,xmax);
761  TH1D hbw2("hbw2","BreitWigner GSL",nch,xmin,xmax);
762 
763 
764  testBreitWigner(r,mu,sqrt(mu),hbw1);
765  testBreitWigner(tr,mu,sqrt(mu),hbw2);
766  testDiff(hbw1,hbw2,"Breit-Wigner");
767 
768  // binomial
769  std::cout << std::endl;
770 
771  int ntot = 100;
772  double p =0.1;
773  xmin = 0;
774  xmax = ntot+1;
775  nch = std::min(1000,ntot+1);
776  TH1D hb1("hb1","Binomial ROOT",nch,xmin,xmax);
777  TH1D hb2("hb2","Binomial GSL",nch,xmin,xmax);
778  TH1D hb3("hb3","Binomial UNR",nch,xmin,xmax);
779 
780 
781  testBinomial(r,ntot,p,hb1);
782  testBinomial(tr,ntot,p,hb2);
783  testBinomial(ur,ntot,p,hb3);
784 #ifdef HAVE_CLHEP
785  RandBinomial crb(eng);
786  TH1D hb4("hb4","Binomial CLHEP",nch,xmin,xmax);
787  testBinomialCLHEP(crb,ntot,p,hp4);
788 #endif
789 
790 
791  testDiff(hb1,hb2,"Binomial ROOT-GSL");
792  testDiff(hb1,hb3,"Binomial ROOT-UNR");
793 
794  // multinomial
795  std::cout << std::endl;
796 
797  TH1D hmt1("hmt1","Multinomial GSL",30,-3,3);
798  TH1D hmt2("hmt2","Multi-Poisson GSL",30,-3,3);
799  TH1D hmt3("hmt3","Gaus",30,-3,3);
800  ntot = 10000;
801  testMultinomial(r,ntot,hmt1,hmt2);
802  hmt3.FillRandom("gaus",ntot);
803  testDiff(hmt1,hmt2,"Multinomial-Poisson");
804  testDiff(hmt1,hmt3,"Multinomial-gaus");
805 
806 
807  // exponential
808  std::cout << std::endl;
809 
810  TH1D he1("he1","Exp ROOT",300,0,20);
811  TH1D he2("he2","Exp GSL",300,0,20);
812 
813  testExp(r,he1);
814  testExp(tr,he2);
815  testDiff(he1,he2,"Exponential");
816 
817  // circle
818  std::cout << std::endl;
819 
820  TH1D hc1("hc1","Circle ROOT",300,-PI,PI);
821  TH1D hc2("hc2","Circle GSL",300,-PI,PI);
822 
823  testCircle(r,hc1);
824  testCircle(tr,hc2);
825  testDiff(hc1,hc2,"Circle");
826 
827 
828  // sphere
829  std::cout << std::endl;
830 
831  TH1D hs1("hs1","Sphere-Phi ROOT",300,-PI,PI);
832  TH1D hs2("hs2","Sphere-Phi GSL ",300,-PI,PI);
833  TH1D hs3("hs3","Sphere-Theta ROOT",300,0,PI);
834  TH1D hs4("hs4","Sphere-Theta GSL ",300,0,PI);
835 
836  testSphere(r,hs1,hs3);
837  testSphere(tr,hs2,hs4);
838  testDiff(hs1,hs2,"Sphere-phi");
839  testDiff(hs3,hs4,"Sphere-theta");
840 
841 
842  return 0;
843 
844 }
845 
846 int main() {
847  return testRandomDist();
848 }
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3159
xmin
float xmin
Definition: THbookFile.cxx:93
TRandom3
Random number generator class based on M.
Definition: TRandom3.h:29
ROOT::Vc::AVX::min
static Vc_ALWAYS_INLINE int_v min(const int_v &x, const int_v &y)
Definition: vector.h:433
TStopwatch::RealTime
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
TH1::Chi2TestX
virtual Double_t Chi2TestX(const TH1 *h2, Double_t &chi2, Int_t &ndf, Int_t &igood, Option_t *option="UU", Double_t *res=0) const
The computation routine of the Chisquare test.
Definition: TH1.cxx:1909
pi
const double pi
Definition: rotationApplication.cxx:282
c1
TCanvas * c1
Definition: legend1.C:2
TRandom2
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:29
testPoisson
void testPoisson(R &r, double mu, TH1D &h)
Definition: testRandomDist.cxx:91
TStopwatch::CpuTime
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
_pythonization.g
g
Definition: _pythonization.py:167
nbins
int nbins[3]
Definition: SparseDataComparer.cxx:22
TH1::GetNbinsX
virtual Int_t GetNbinsX() const
Definition: TH1.h:296
testSphere
void testSphere(R &r, TH1D &h1, TH1D &h2)
Definition: testRandomDist.cxx:594
TH1::GetEntries
virtual Double_t GetEntries() const
return the current number of entries
Definition: TH1.cxx:4051
ROOT::Math::sqrt
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
Definition: UnaryOperators.h:283
NEVT
#define NEVT
Definition: testRandomDist.cxx:38
testExp
void testExp(R &r, TH1D &h)
Definition: testRandomDist.cxx:540
TH1::GetBinWidth
virtual Double_t GetBinWidth(Int_t bin) const
return bin width for 1D historam Better to use h1.GetXaxis().GetBinWidth(bin)
Definition: TH1.cxx:8492
sqrt
double sqrt(double)
x2
static const double x2[5]
Definition: RooGaussKronrodIntegrator1D.cxx:345
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
testBinomialCLHEP
void testBinomialCLHEP(R &r, int ntot, double p, TH1D &h)
Definition: testRandomDist.cxx:457
x
Double_t x[n]
Definition: legend1.C:17
genPoisson2
unsigned int genPoisson2(R &r, double mu)
Definition: testRandomDist.cxx:171
TH1::GetBinLowEdge
virtual Double_t GetBinLowEdge(Int_t bin) const
return bin lower edge for 1D historam Better to use h1.GetXaxis().GetBinLowEdge(bin) ...
Definition: TH1.cxx:8481
testGaus
void testGaus(R &r, double mu, double sigma, TH1D &h)
Definition: testRandomDist.cxx:351
ROOT::Math::normal_cdf
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
Definition: ProbFuncMathCore.cxx:234
h2
TH2D * h2
Definition: fit2dHist.C:45
TH1::DrawCopy
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
Definition: TH1.cxx:2925
h1
TH1F * h1
Definition: legend1.C:5
z
Float_t z[5]
Definition: Ifit.C:16
testBinomial
void testBinomial(R &r, int ntot, double p, TH1D &h)
Definition: testRandomDist.cxx:432
fillHist
static bool fillHist
Definition: testRandomDist.cxx:48
ROOT::Math::Cephes::gamma
double gamma(double x)
Definition: SpecFuncCephes.cxx:339
t
TThread * t[5]
Definition: threadsh1.C:13
TH1::FillRandom
virtual void FillRandom(const char *fname, Int_t ntimes=5000)
Fill histogram following distribution in function fname.
Definition: TH1.cxx:3330
r
ROOT::R::TRInterface & r
Definition: Object.C:4
TH3D
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:309
v
SVector< double, 2 > v
Definition: Dict.h:5
TH1::SetBinContent
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:8543
testRandomDist
int testRandomDist()
Definition: testRandomDist.cxx:657
m
TMarker * m
Definition: textangle.C:8
testMultinomial
void testMultinomial(R &r, int ntot, TH1D &h1, TH1D &h2)
Definition: testRandomDist.cxx:482
floor
double floor(double)
qtexample.w
tuple w
Definition: qtexample.py:51
main
int main()
Definition: testRandomDist.cxx:846
TH3::Fill
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH3.cxx:275
xmax
float xmax
Definition: THbookFile.cxx:93
TH1D
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:613
TMath::Pi
Double_t Pi()
Definition: TMath.h:44
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:495
x1
static const double x1[5]
Definition: RooGaussKronrodIntegrator1D.cxx:327
testDiff
void testDiff(TH1D &h1, TH1D &h2, const std::string &name="")
Definition: testRandomDist.cxx:52
atan2
double atan2(double, double)
type
int type
Definition: TGX11.cxx:120
testBreitWigner
void testBreitWigner(R &r, double mu, double gamma, TH1D &h)
Definition: testRandomDist.cxx:406
y
Double_t y[n]
Definition: legend1.C:17
ROOT::Vc::AVX::max
static Vc_ALWAYS_INLINE int_v max(const int_v &x, const int_v &y)
Definition: vector.h:440
testLandau
void testLandau(R &r, TH1D &h)
Definition: testRandomDist.cxx:378
name
#define name(a, b)
Definition: linkTestLib0.cpp:5
testCircle
void testCircle(R &r, TH1D &h)
Definition: testRandomDist.cxx:565
PI
#define PI
Definition: testRandomDist.cxx:32
testPoisson2
void testPoisson2(R &r, double mu, TH1D &h)
Definition: testRandomDist.cxx:219
findName
std::string findName(const R &r)
Definition: testRandomDist.cxx:74
tan
double tan(double)
TMath::LnGamma
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:490
exp
double exp(double)
TH2::Fill
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:287
n
const Int_t n
Definition: legend1.C:16
R
TRandom3 R
a TMatrixD.
Definition: testIO.cxx:28
log
double log(double)
ROOT::Math::Random
Documentation for the Random class.
Definition: Random.h:41
TH2D
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:297
genPoisson
unsigned int genPoisson(R &r, double mu)
Definition: testRandomDist.cxx:118
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:30
testFlat
void testFlat(R &r, TH1D &h)
Definition: testRandomDist.cxx:325