64 for(
int i = 0; i < p[0]; i++)
74 for(
int i = 0; i < p[0]; i++)
85 for(
int i = 0; i < p[0]; i++)
98 double integral_num(
unsigned int dim,
double*
a,
double*
b,
double* p,
double & value,
double & error)
100 if (
verbose) std::cout <<
"\nTesting IntegratorMultiDim class.." << std::endl;
106 unsigned int nmax = (
unsigned int) 1.E8;
108 unsigned int size = 0;
115 std::cout.precision(12);
116 std::cout <<
"result: \t";
117 std::cout << ig1.
Result() <<
"\t" <<
"error: \t" << ig1.
Error() << std::endl;
118 std::cout <<
"Number of function evaluations: " << ig1.
NEval() << std::endl;
119 std::cout <<
"Time using IntegratorMultiDim: \t" << timer.
RealTime() << std::endl;
120 std::cout <<
"------------------------------------" << std::endl;
122 else { std::cout <<
" . "; }
133 std::vector<double>
integral_MC(
unsigned int dim,
double*
a,
double*
b,
double* p,
double * value,
double * error)
140 std::cout <<
"\nTesting GSLMCIntegrator class.." << std::endl;
141 std::cout <<
"\t VEGAS.. " << std::endl;
142 std::cout <<
"" << std::endl;
144 else { std::cout <<
"."; }
164 error[0] = ig1.
Error();
167 std::cout <<
"result: \t";
168 std::cout << ig1.
Result() <<
"\t" <<
"error: \t" << ig1.
Error() << std::endl;
169 std::cout <<
"sigma: \t" << ig1.
Sigma();
170 std::cout <<
"\t" <<
"chi2: \t" << ig1.
ChiSqr() << std::endl;
171 std::cout << std::endl;
172 std::cout <<
"Time using GSLMCIntegrator::VEGAS :\t" << timer.
RealTime() << std::endl;
175 std::cout <<
"" <<std::endl;
176 std::cout <<
"------------------------------------" << std::endl;
177 std::cout <<
"\t MISER.. " << std::endl;
178 std::cout <<
"" << std::endl;
180 else { std::cout <<
"."; }
203 error[1] = ig2.
Error();
205 std::cout <<
"result: \t";
206 std::cout << ig2.
Result() <<
"\t" <<
"error: \t" << ig2.
Error() << std::endl;
208 std::cout <<
"Time using GSLMCIntegrator::MISER :\t" << timer.
RealTime() << std::endl;
209 std::cout <<
"" << std::endl;
210 std::cout <<
"------------------------------------" << std::endl;
211 std::cout <<
"\t PLAIN.. " << std::endl;
213 else { std::cout <<
"."; }
223 error[2] = ig3.
Error();
226 std::cout <<
"" << std::endl;
227 std::cout <<
"result: \t";
228 std::cout << ig3.
Result() <<
"\t" <<
"error: \t" << ig3.
Error() << std::endl;
229 std::cout <<
"Time using GSLMCIntegrator::PLAIN :\t" << timer.
RealTime() << std::endl;
230 std::cout <<
"" << std::endl;
232 std::vector<double>
result(3);
233 result[0] = timeVegas; result[1] = timeMiser; result[2] = timePlain;
240 unsigned int Nmax =
NMAX;
241 unsigned int size = Nmax-1;
244 TH1D *num_performance =
new TH1D(
"cubature",
"", size, 1.5, Nmax+.5);
245 TH1D *Vegas_performance =
new TH1D(
"VegasMC",
"", size, 1.5, Nmax+.5);
246 TH1D *Miser_performance =
new TH1D(
"MiserMC",
"", size, 1.5, Nmax+.5);
247 TH1D *Plain_performance =
new TH1D(
"PlainMC",
"", size, 1.5, Nmax+.5);
251 for(
unsigned int N = 2;
N <=Nmax;
N++)
254 std::cout<<
"*********************************************" << std::endl;
255 std::cout<<
"Number of dimensions: "<<
N << std::endl;
258 std::cout <<
"\n\tdim="<<
N <<
" : ";
261 double *
a =
new double[
N];
262 double *
b =
new double[
N];
265 for (
unsigned int i=0; i <
N; i++)
272 double valMC[3], errMC[3];
273 double timeNumInt =
integral_num(N, a, b, p, val0, err0);
274 std::vector<double> timeMCInt =
integral_MC(N, a, b, p, valMC, errMC);
283 for (
int j = 0; j < 3; ++j) {
285 Error(
"testMCIntegration",
"Result is not consistent for dim %d between adaptive and MC %d ",N,j);
308 num_performance->
SetTitle(
"comparison of performance");
324 legend->
AddEntry(h1,
"Cubature",
"f");
325 legend->
AddEntry(h2,
"MC Vegas",
"f");
326 legend->
AddEntry(h3,
"MC Miser",
"f");
327 legend->
AddEntry(h4,
"MC Plain",
"f");
333 std::cout <<
"\nTest Timing results\n";
334 std::cout <<
" N dim \t Adaptive MC Vegas MC Miser MC Plain \n";
335 for (
unsigned int i=1; i<=size; i++) {
337 std::cout.precision(6);
339 std::cout <<
"\t " << std::setw(12) << num_performance->
GetBinContent(i) << std::setw(12) << Vegas_performance->
GetBinContent(i)
346 int main(
int argc,
char **argv)
351 for (
Int_t i=1 ; i<argc ; i++) {
352 std::string arg = argv[i] ;
361 cout <<
"Usage: " << argv[0] <<
" [-g] [-v]\n";
363 cout <<
" -g : graphics mode\n";
364 cout <<
" -v : verbose mode";
int NEval() const
return number of function evaluations in calculating the integral
virtual void SetBarOffset(Float_t offset=0.25)
static long int sum(long int i)
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
double Result() const
return the type of the integration used
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
This class displays a legend box (TPaveText) containing several legend entries.
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
void SetFunction(const IMultiGenFunction &f)
set the integration function (must implement multi-dim function interface: IBaseFunctionMultiDim) ...
double Integral(const GSLMonteFuncPointer &f, unsigned int dim, double *a, double *b, void *p=0)
evaluate the Integral of a function f over the defined hypercube (a,b)
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Double_t Sum(const double *x, const double *p)
double Result() const
return result of integration
void SetFunction(const IMultiGenFunction &f)
method to set the a generic integration function
int main(int argc, char **argv)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Double_t SingularFun(const double *x, const double *p)
virtual void SetBarWidth(Float_t width=0.5)
void Stop()
Stop the stopwatch.
double Error() const
return the estimate of the absolute Error of the last Integral calculation
double ChiSqr()
returns chi-squared per degree of freedom for the estimate of the integral in the Vegas algorithm ...
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
double integral_num(unsigned int dim, double *a, double *b, double *p, double &value, double &error)
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
Copy this histogram and Draw in the current pad.
void Error(const char *location, const char *msgfmt,...)
double Integral(const double *xmin, const double *xmax)
evaluate the integral with the previously given function between xmin[] and xmax[] ...
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
std::vector< double > integral_MC(unsigned int dim, double *a, double *b, double *p, double *value, double *error)
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...
double Error() const
return integration error
tomato 1-D histogram with a double per channel (see TH1 documentation)}
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
void SetGenerator(GSLRandomEngine &r)
set random number generator
double Sigma()
set parameters for PLAIN method
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Double_t SimpleFun(const double *x, const double *p)
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
virtual void SetTitle(const char *title)
Change (i.e.
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
class for adaptive quadrature integration in multi-dimensions using rectangular regions.
WrappedParamFunction class to wrap any multi-dimensional function pbject implementing the operator()(...
Documentation for the Random class.