ROOT   Reference Guide
multidimSampling.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_math
3/// \notebook
4/// Example of sampling a multi-dim distribution using the DistSampler class
5/// NOTE: This tutorial must be run with ACLIC
6///
7/// \macro_image
8/// \macro_code
9///
10/// \author Lorenzo Moneta
11
12// function (a 4d gaussian)
13#include "TMath.h"
14#include "TF2.h"
15#include "TStopwatch.h"
16#include "Math/DistSampler.h"
19#include "Math/Factory.h"
20
21#include "TKDTreeBinning.h"
22
23#include "TTree.h"
24#include "TFile.h"
25#include "TMatrixDSym.h"
26#include "TVectorD.h"
27#include "TCanvas.h"
28#include <cmath>
29
30// Gauss ND function
31// make a class in order to avoid constructing the
32// matrices for every call
33// This however requires that the code must be compiled with ACLIC
34
35bool debug = false;
36
37// Define the GausND strcture
38struct GausND {
39
40 TVectorD X;
41 TVectorD Mu;
42 TMatrixDSym CovMat;
43
44 GausND( int dim ) :
45 X(TVectorD(dim)),
46 Mu(TVectorD(dim)),
47 CovMat(TMatrixDSym(dim) )
48 {}
49
50 double operator() (double *x, double *p) {
51 // 4 parameters
52 int dim = X.GetNrows();
53 int k = 0;
54 for (int i = 0; i<dim; ++i) { X[i] = x[i] - p[k]; k++; }
55 for (int i = 0; i<dim; ++i) {
56 CovMat(i,i) = p[k]*p[k];
57 k++;
58 }
59 for (int i = 0; i<dim; ++i) {
60 for (int j = i+1; j<dim; ++j) {
61 // p now are the correlations N(N-1)/2
62 CovMat(i,j) = p[k]*sqrt(CovMat(i,i)*CovMat(j,j));
63 CovMat(j,i) = CovMat(i,j);
64 k++;
65 }
66 }
67 if (debug) {
68 X.Print();
69 CovMat.Print();
70 }
71
72 double det = CovMat.Determinant();
73 if (det <= 0) {
74 Fatal("GausND","Determinant is <= 0 det = %f",det);
75 CovMat.Print();
76 return 0;
77 }
78 double norm = std::pow( 2. * TMath::Pi(), dim/2) * sqrt(det);
79 // compute the gaussians
80 CovMat.Invert();
81 double fval = std::exp( - 0.5 * CovMat.Similarity(X) )/ norm;
82
83 if (debug) {
84 std::cout << "det " << det << std::endl;
85 std::cout << "norm " << norm << std::endl;
86 std::cout << "fval " << fval << std::endl;
87 }
88
89 return fval;
90 }
91};
92
93// Use the Math namespace
94using namespace ROOT::Math;
95
96void multidimSampling() {
97
98
99 const int N = 10000;
100 /*const int NBin = 1000;*/
101 const int DIM = 4;
102
103 double xmin[] = {-10,-10,-10, -10};
104 double xmax[] = { 10, 10, 10, 10};
105 double par0[] = { 1., -1., 2, 0, // the gaussian mu
106 1, 2, 1, 3, // the sigma
107 0.5,0.,0.,0.,0.,0.8 }; // the correlation
108
109 const int NPAR = DIM + DIM*(DIM+1)/2; // 14 in the 4d case
110 // generate the sample
111 GausND gaus4d(4);
112 TF1 * f = new TF1("functionND",gaus4d,0,1,14);
113 f->SetParameters(par0);
114
115 double x0[] = {0,0,0,0};
116 // for debugging
117 if (debug) f->EvalPar(x0,0);
118 debug = false;
119
121 for (int i = 0; i < NPAR; ++i ) {
122 if (i < DIM) f->SetParName(i, name.Format("mu_%d",i+1) );
123 else if (i < 2*DIM) f->SetParName(i, name.Format("sig_%d",i-DIM+1) );
124 else if (i < 2*DIM) f->SetParName(i, name.Format("sig_%d",i-2*DIM+1) );
125 }
126
127 /*ROOT::Math::DistSamplerOptions::SetDefaultSampler("Foam");*/
128 DistSampler * sampler = Factory::CreateDistSampler();
129 if (sampler == 0) {
130 Info("multidimSampling","Default sampler %s is not available try with Foam ",
133 }
134 sampler = Factory::CreateDistSampler();
135 if (sampler == 0) {
136 Error("multidimSampling","Foam sampler is not available - exit ");
137 return;
138 }
139
140 sampler->SetFunction(*f,DIM);
141 sampler->SetRange(xmin,xmax);
142 bool ret = sampler->Init();
143
144 std::vector<double> data1(DIM*N);
145 double v[DIM];
146 TStopwatch w;
147
148 if (!ret) {
149 Error("Sampler::Init","Error initializing unuran sampler");
150 return;
151 }
152
153 // generate the data
154 w.Start();
155 for (int i = 0; i < N; ++i) {
156 sampler->Sample(v);
157 for (int j = 0; j < DIM; ++j)
158 data1[N*j + i] = v[j];
159 }
160 w.Stop();
161 w.Print();
162
163 // fill tree with data
164 TFile * file = new TFile("multiDimSampling.root","RECREATE");
165 double x[DIM];
166 TTree * t1 = new TTree("t1","Tree from Unuran");
167 t1->Branch("x",x,"x[4]/D");
168 for (int i = 0; i < N; ++i) {
169 for (int j = 0; j < DIM; ++j) {
170 x[j] = data1[i+N*j];
171 }
172 t1->Fill();
173 }
174
175 // plot the data
176 t1->Draw("x[0]:x[1]:x[2]:x[3]","","candle");
177 TCanvas * c2 = new TCanvas();
178 c2->Divide(3,2);
179 int ic=1;
180 c2->cd(ic++);
181 t1->Draw("x[0]:x[1]");
182 c2->cd(ic++);
183 t1->Draw("x[0]:x[2]");
184 c2->cd(ic++);
185 t1->Draw("x[0]:x[3]");
186 c2->cd(ic++);
187 t1->Draw("x[1]:x[2]");
188 c2->cd(ic++);
189 t1->Draw("x[1]:x[3]");
190 c2->cd(ic++);
191 t1->Draw("x[2]:x[3]");
192
193 t1->Write();
194 file->Close();
195
196}
#define f(i)
Definition: RSha256.hxx:104
void Info(const char *location, const char *msgfmt,...)
void Error(const char *location, const char *msgfmt,...)
void Fatal(const char *location, const char *msgfmt,...)
#define N
char name[80]
Definition: TGX11.cxx:109
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
double pow(double, double)
double sqrt(double)
double exp(double)
TRObject operator()(const T1 &t1) const
static void SetDefaultSampler(const char *type)
static const std::string & DefaultSampler()
Interface class for generic sampling of a distribution, i.e.
Definition: DistSampler.h:57
const double * Sample()
sample one event and rerturning array x with coordinates
Definition: DistSampler.h:169
void SetRange(double xmin, double xmax, int icoord=0)
set range in a given dimension
Definition: DistSampler.cxx:40
void SetFunction(Function &func, unsigned int dim)
set the parent function distribution to use for sampling (generic case)
Definition: DistSampler.h:72
virtual bool Init(const char *="")
initialize the generators with the given algorithm Implemented by derived classes who needs it (like ...
Definition: DistSampler.h:101
The Canvas class.
Definition: TCanvas.h:27
1-Dim function class
Definition: TF1.h:210
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
void Print(Option_t *name="") const
Print the matrix as a table of elements.
virtual Double_t Determinant() const
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
TMatrixTSym< Element > & Similarity(const TMatrixT< Element > &n)
Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb) This is a similarity transform w...
Stopwatch class.
Definition: TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
Basic string class.
Definition: TString.h:131
A TTree represents a columnar dataset.
Definition: TTree.h:78
Int_t GetNrows() const
Definition: TVectorT.h:75
void Print(Option_t *option="") const
Print the vector as a list of elements.
Definition: TVectorT.cxx:1364
Double_t x[n]
Definition: legend1.C:17
return c2
Definition: legend2.C:14
constexpr Double_t Pi()
Definition: TMath.h:38
Definition: file.py:1
auto * t1
Definition: textangle.C:20