Description

hep-mc is a C++11 template library for Monte Carlo integration. The following integration algorithms are available:

PLAIN Integrator,
VEGAS Integrator [3] [4], and a
Multi Channel Integrator with adaptive weight optimization [2].

How to use the documentation

If you are in a hurry, have a look at a few examples which showcase each feature of this library separately. If you want to browse the documentation instead, first have a look at Integrand Functions, to see how a function that can be integrated must look like. Next, you should decide which integration algorithm (one from the list above) you want to use; go to their respective page to decide and to learn how to use them. Finally, you probably would like to use the result in the program (by default it is printed to the standard output). Read up on Checkpointing system to get the Results. For more complicated use cases, learn how to use Differential Distributions to simultaneously integrate arbitrarily many (similar) integrands.

Example

The following example illustrates how to integrate the square-function using VEGAS:

#include "hep/mc.hpp"
#include <cstddef>
#include <iostream>
#include <vector>
// the function that shall be integrated
double square(hep::mc_point<double> const& x)
{
    return 3.0 * x.point()[0] * x.point()[0];
}
int main()
{
    // print reference result
    std::cout << ">> computing integral of 3*x^2 from 0 to 1 which is 1.0\n\n";
    // perform 5 iteration with 1000 calls each; this function will also call
    // vegas_verbose_callback after each iteration which in turn prints the
    // individual iterations
    auto results = hep::vegas(
        hep::make_integrand<double>(square, 1),
        std::vector<std::size_t>(5, 1000)
    ).results();
    // results contains the estimations for each iteration. We could take the
    // result from last iteration, but here we instead choose to combine the
    // results of all iterations but the first one in a cumulative result
    auto result = hep::accumulate<hep::weighted_with_variance>(
        results.begin() + 1, results.end());
    double chi_square_dof = hep::chi_square_dof<hep::weighted_with_variance>(
        results.begin() + 1, results.end());
    // print the cumulative result
    std::cout << ">> cumulative result (excluding first iteration):\n>> N="
        << result.calls() << " I=" << result.value() << " +- " << result.error()
        << " chi^2/dof=" << chi_square_dof << "\n";
    return 0;
}

This program will produce the following output:

computing integral of x^2 from 0 to 1 which is 0.333333

iteration 0 finished.
this iteration: N=1000 E=0.33691 +- 0.00966589 (2.86898%)
all iterations: N=1000 E=0.33691 +- 0.00966589 (2.86898%) chi^2/dof=inf

iteration 1 finished.
this iteration: N=1000 E=0.336297 +- 0.00381044 (1.13306%)
all iterations: N=2000 E=0.33638 +- 0.00354493 (1.05385%) chi^2/dof=0.00347637

iteration 2 finished.
this iteration: N=1000 E=0.332836 +- 0.00176979 (0.531731%)
all iterations: N=3000 E=0.333543 +- 0.00158343 (0.47473%) chi^2/dof=0.401838

iteration 3 finished.
this iteration: N=1000 E=0.332933 +- 0.00132342 (0.397504%)
all iterations: N=4000 E=0.333184 +- 0.00101545 (0.304772%) chi^2/dof=0.296971

iteration 4 finished.
this iteration: N=1000 E=0.333122 +- 0.000828935 (0.248838%)
all iterations: N=5000 E=0.333147 +- 0.000642145 (0.192751%) chi^2/dof=0.22328

cumulative result (without first iteration):
N=4000 I=0.33313 +- 0.000643567 chi^2/dof=0.246959