format: text

text:
# Contents

## 1 Introduction

This work describe a new version of an algorithm for producing random points according to an arbitrary, user-defined distribution in the $`n`$-dimensional space – much improved with respect to the original version of Ref. . The new implementation is realized in the C++ programming language in a fully object-oriented manner111 The early C++ version of Foam was coded by M. Ciesla and M. Slusarczyk . It was a translation of a version 1.x from Fortran77 to C++. Analogous translation to the JAVA language was also done.. Since the changes, with respect to Ref. , both in the algorithm and in the implementation are quite essential, a complete description of the method and the new code is provided, instead of only an update of Ref. .

For the problem of function minimization, integration (quadrature), there are plenty of general-purpose programs, which can be applied to an arbitrary user-defined function. “General-purpose” means that all these tools work, in principle, for a very wide range of user-defined functions. For the multidimensional Monte Carlo simulation problem, that is for the problem of randomly generating points according to a given $`n`$-dimensional distribution, there is precious few examples of the General-Purpose Monte Carlo Simulators (GPMCS), that is programs that work (in principle) for arbitrary distribution . As in these works, here we are concerned mainly with the MC applications to high energy physics, that is to simulation of the differential distributions in the multi-particle Lorentz invariant phase space provided the Quantum Field Theory, see also classic Ref.  on this subject. An example of the work on GPMCS applied to other fields see the interesting works222In these works there is far more emphasis on the parallel computing aspects of the integration (not simulation) than on the cell geometry, as compared with our work. of Refs. .

Let us also note that the two-dimensional cellular MC sampler VESKO2 with the primitive binary split was already included in the program LESKO-F of Ref.  long time ago. Very similar programs were also described in Refs. 333We thank S. Kawabata for drawing our attention to these works.. Still another very interesing class of GPMCSs for the high energy physics, based of the Metropolis algorith , is described in Ref. .

Two essential reasons for a realtive scarcity of GPMCSs, as compared to other programming tools, are the lack of novel ideas for an efficient algorithm and the need of much CPU power and memory – only recently available or affordable.

GPMCS is essentially a random-number generator of points in multidimensional space with a non-uniform user-defined probability distribution. Inevitably the GPMCS works in two stages: exploration and generation444Exploration and generation could be done simultaneously, at the expense of complications in the algorithm and the code.. During an exploration phase, the GPMCS is “digesting” the entire shape of the $`n`$-dimensional distribution $`\rho (\stackrel{}{x})`$ to be generated, memorizing it as efficiently as possible, using all available CPU processing power and memory555The procedure of memorizing multidimensional distribution $`\rho (\stackrel{}{x})0`$ is a kind of interpolation, in which the grid of cells is denser in places where the distribution peaks and/or varies strongly.. In Foam, the exploration phase is the phase of the build-up of the system of cells covering entirely the integration space, which will be called “foam”, produced in the process of the binary split of the cells. In the generation phase, GPMCS provides a method of the MC generation of the random points $`\stackrel{}{x}`$ exactly according to $`\rho (\stackrel{}{x})`$. The vector $`x=\stackrel{}{x}=(x_1,x_2,\mathrm{},x_n)`$ will also be called in the following a Monte Carlo event. In Foam, the MC generation is very simple: a cell is chosen randomly, and next, a point is generated within the cell with uniform distribution; see below for more details. The value of the integrand is already estimated in the exploration; it can be calculated with an arbitrary precision in the generation phase.

During the exploration Foam constructs a distribution $`\rho ^{}(x)`$, which is uniform within each cell, and is used for the MC generation. Events are weighted with the weight $`w=\rho /\rho ^{}`$. The quality of the distribution of this weight, measured in terms of the weight distribution parameters, such as variance and ratio of maximum to average, is determined by the quality of the exploration. The basic principle of the Foam algorithm is that the parameters of the anticipated “target weight distribution” in the MC generation phase are used as a driving force guiding the cell build-up (exploration). In the case of a successful exploration, weighted MC events can be turned efficiently into unweighted ones with the usual rejection method, that is with a small rejection rate.

Since the exploration phase may be CPU-time consuming, it is natural to expect that GPMCS has a built-in mechanism of persistency, i.e., there is a mechanism of writing into a mass storage (computer disk) the whole information on the memorized shape of the distribution obtained from the exploration phase, such that the generation of the MC events can be (re)started at any later time, without any need of repeating the time-consuming exploration. One small step further is to require that the generation of events with GPMCS can be stopped at any time, the entire status of the GPMCS can be written on the disk, and the generation of the next event can be resumed at any later time upon reading the stored information; the next generated event will be such as if there had been no break in the generation process. In fact, this is what we shall really mean in the following as a persistency mechanism for GPMCS, and what is actually implemented in the Foam. There, the persistency is realized using the ROOT666 The use of ROOT is optional in Foam. However, version of Foam without ROOT does not feature any kind of persistency. package .

The GPMCS programs will always be limited to “small dimensions”. With currently available computers, “small” means in practice $`n10`$, up to $`n16`$ for “mildly” singular distributions. This is already quite satisfactory, especially if we remember that this limit will pushed higher, as the available hardware gets more powerful, without any need of modifying the existing code777The present implementation of Foam is fully based on the dynamic allocation of the memory and the space dimension is a user-defined parameter.. Twelve years from now, with portable computer featuring a 100 GHz processor and 1 TByte disk the same version of Foam will work efficiently for even higher dimensions.

Foam has been developed having in mind that it will be used as part of a bigger MC program, typically, to generate a subset of variables in which a model distribution is the most singular (has strong peaks). This is why we are not so much concerned by the fact that the cellular algorithm of Foam is inefficient for, say, 150 variables. The user is supposed to select $`n16`$ “wild variables”  and apply Foam to them. For the remaining “mild variables”, Foam may merely serve as a provider of the uniform random numbers, if the user of Foam wishes to exploit that option. On the other hand, for smaller MC problems, Foam may play the role of a “stand-alone MC generator” or “stand-alone MC integrator”. Also, from the following description of the various modes of the use of Foam it will be clear that the subprogram providing the model distribution to Foam can have quite a complicated structure. Nevertheless, this user-provided part of the program will be smaller than a solution without Foam, because Foam provides for essential functionalities concerning the optimization of the MC weight distribution. This remark is especially true for the case of implementation of multibranching with the help of Foam.

It is worth mentioning that Foam is not based on the “principle of factorizability” of the integrand distribution, $`\rho (\stackrel{}{x})=_1^n\rho _i(x_i)`$, on which VEGAS-family programs are built .

The outline of the paper is the following: in Section 2 we describe the cellular Foam algorithm, delegating the description of the cell division procedure to Section 3. Section 4 is devoted to the description of the Foam code in C++. The use of Foam is described in Section 5 and examples of the numerical results (MC efficiency) are given in Section 6. Conclusions and Appendices on the variance minimization finalize the paper.

## 2 The Foam algorithm

As already mentioned, the execution of the Foam algorithm is clearly separated into the first stage of the “distribution exploration” consisting of the “build-up of the foam of cells”, which in a sense memorizes the $`n`$-dimensional shape of the distribution, and the second stage of the actual “MC generation”, see Fig. 1. The most essential part of the present Foam algorithm is the procedure of the binary split of the cell, in which it is decided which cell is picked up for the next split and the necessary parameters of the geometry of the cell split are determined. This part of the Foam algorithm description is delegated to the next section. In the present section we describe other, more general, aspects of the Foam algorithm.

### 2.1 Cellular exploration of the distribution

The most obvious method to minimize the variance or maximum weight of the target weight distribution in the MC generation, already considered some 40 years ago, is to split the integration domain into many cells, so that the distribution $`\rho (\stackrel{}{x})`$ is approximated by $`\rho ^{}(\stackrel{}{x})`$, which is constant within each cell888 For the MC simulation, our main aim, a more sophisticated interpolation of $`\rho (\stackrel{}{x})`$ within a cell does not seem to be worth the effort – it would be interesting if our main aim was the integration of $`\rho (\stackrel{}{x})`$.. This is a cellular class of the general-purpose MC algorithms999The term “stratified sampling”, used in the literature, has in our opinion a narrower meaning than “cellular class”..

The immediate questions are: What kind or shapes of the cells to use and how to cover the integration domain with cells? The reader may find in Ref.  an example of a rather general discussion of these questions. In the Foam program the user may opt for one of the three geometries of the cells: (1) simplices, (2) hyperrectangles and (3) Cartesian products of simplices and hyperrectangles. For these particular types of cells there exists an efficient method of parametrizing them in the computer memory and handling their geometry.

The system of many cells can be created and reorganized all at once, as in VEGAS-type programs , or in a more evolutionary way, as the cell split process of this work. In the Foam algorithm we rely on the binary split of cells. Starting from the entire integration domain (unit hyperrectangle or simplex) cells are split into two daughter cells, step by step, until the user-defined memory limit is reached. The choice of next cell to be split and the geometry of the split in the exploration phase are driven by the “target weight distribution” of the generation process, see Section 3. The important advantage of any cell split algorithm is that it assures automatically the full coverage of the integration domain – simply because the primary root cell is identical to the entire integration domain and the two daughter cells always cover the parent cell entirely. The problem of blind spots discussed in Ref.  is avoided by construction.

In the early version of Foam , there was a possibility in the algorithm that the “unsuccessful” branch in the tree of all cells can be erased and rebuild. This was called “collapse” and “rebuild”. In the present version this option was removed101010A “flush method”, which erases the entire foam of cells from the computer memory and allows for its reinitialization is, however, available., because the experience with many testing functions has shown that the algorithm of the cell build-up is rather “deterministic” and the “rebuild” procedure was usually leading to a new branch of foam with about the same features as the old one.

Let us finally remark that the version of the cellular algorithm presented in this paper is, in fact, result of many experiments with different variants of the cellular algorithm. The presented version is the best one out of several development versions. In the code one may still see some “hooks” and unused features (class members or methods) related to these alternative variants. We have left them just in the case that some new idea of improving the algorithm emerges, or for certain kinds of debugging/testing.

### 2.2 Variance reduction versus maximum weight reduction

In the construction of the Foam algorithm, most effort was invested into a minimization of the ratio of the maximum weight to the average weight $`w_{\mathrm{max}}/w`$. This parameter is essential, if we want to transform variable-weight events into $`w=1`$ events, at the latter stage of the MC generation111111We provide optionally in Foam for the rejection leading to $`w=1`$ events..

Minimizing the maximum weight $`w_{\mathrm{max}}`$ is not the same as minimizing the variance $`\sigma =\sqrt{w^2w^2}`$. Minimizing $`w_{\mathrm{max}}`$ may be more difficult – but it is worth an effort because Foam is really meant to be a part of a bigger MC program, where it is usually essential that the “inner part” of the program provides events with an excellent weight distribution, or even $`w=1`$ events. Nevertheless, minimizing the variance is also implemented in Foam and available optionally. It can be useful if one is satisfied with the variable-weight events, and/or if the main aim is the evaluation of the integral and not the MC simulation.

The difference between the above two options is well illustrated in Fig. 2, which shows two examples of the evolution of the MC weight distribution due to a gradual increase of the number of cells. For the default configuration, Foam optimises the ratio $`w_{\mathrm{max}}/w`$. This case is shown in plots (a–c) in Fig. 2. Here, the weight distribution features a sharper and sharper drop of the weight distribution at $`w=1`$, while increasing the number of cells. Also, the average weight increases gradually and the weight distribution gets narrower. The optional case of the optimization of $`\sigma /w`$ is shown in plots (d–f) of Fig. 2. In this case the variance is decreasing with the growing numbers of cells. On the other hand, the maximum weight is noticably higher than before. All weight distributions were obtained for the same 2-dimensional testing function $`\rho _b(x)`$, used also in Section 6.

### 2.3 Hyperrectangles or simplices?

In Ref.  simplical cells have been chosen instead of simpler hyperrectangles, mainly because of the author’s “prejudice” that simplices may adapt more efficiently to complicated singularities in the distribution $`\rho (x)`$ spanned along subspaces, not necessarily parallel to the axes of the global reference frame. Hyperrectangles tend to remember the orientation of the parent hyperrectangle, while simplices feature, in principle, a kind of “angular mobility”, i.e. they can forget the orientation of great-grand-parents, and adapt to the orientation of the singularity in $`\rho (x)`$. However, an experience with tens of testing functions has shown that in many cases hyperrectangles provide the same or even better final MC efficiency than simplices, for the same number of cells. Moreover, simplices have certain additional disadvantages. At present, Foam with simplices is practically limited to rather low dimensions $`n5`$, because in most cases the entire integration domain is a unit hyperrectangle, which has to be divided into $`n!`$ simplices, where $`n!`$ quickly becomes a large number121212 Mapping of the hyperrectangle into a simplex is possible, however, one should use transformation with $`|x/y|0`$, in order to avoid producing nasty singularities located at the vertices, edges and walls of the simplex.. This limitation is, of course, not valid, if the entire integration domain is actually a simplex of high dimensionality instead of a hyperrectangle. (Foam can be configured to start cell evolution from a simplex or Cartesian product of a simplex and a hyperrectangle.) Furthermore, geometry manipulations in the simplical case require calculation of many determinants – this slows down the program execution at higher dimensions. In addition, in the present implementation, the memory consumption in a simplical foam build-up is $`16\times n`$ Bytes/cell, while for hyperrectangles we have found a method that limits memory consumption to below $`80`$ Bytes/cell independently of $`n`$, see Section 2.6. We can therefore reach easily the level of $`10^6`$ hyperrectangular cells at any dimension (in practice $`n16`$) and about $`50,000`$ simplical cells, for $`n5`$. As we see, hyperrectangular foam seems to win on many fronts. Nevertheless, we keep simplical foam as an option, because in certain applications one will definitely encounter distributions for which it turns out to be more efficient to use simplices, in spite of all their limitations, at least for a subset of the integration variables.

### 2.4 Build up of the foam and data organization

The foam of cells is built-up by starting from the root cell, which is the entire integration domain, through a process of binary split of a parent cell into two daughter cells. The root cell is either a unit hypercube $`0x_i1`$ (default) or a simplex $`0x_1x_2x_3\mathrm{}x_n1`$. Also a Cartesian product of these two shapes available on option. Any cell being a product of the cell split can be also a hyperrectangle, a simplex or Cartesian product of the $`k`$-dimensional hyperrectangle and $`n`$-dimensional simplex, with the total dimensionality $`k+n`$. If the starting root cell is a hypercube and cells are simplical (or of mixed type) then the root cell is immediately divided into $`n!`$ simplical (or mixed type) cells.

Each cell is explored immediately after its creation. In the exploration of the cell, about $`100`$-$`1000`$ MC events (the user may define this number) are generated inside the cell with flat (uniform) distribution and using MC weight equal to $`\rho (x)`$; certain averages and certain integrals over the cell are estimated. Also, the best geometry of the binary split of the cell is established and recorded for future use. In this way, every created cell is ready for an immediate split. The determination of the best split is described in fine detail in Section 3. In the exploration the estimate of the integral $`R_I=_{\omega _I}\rho (x)𝑑x^n`$ is calculated for each cell $`\omega _I`$. Far more important is another functional $`R_{loss}|_I=_{\omega _I}\rho _{loss}(x)𝑑x^n`$, see Section 3 for its definition, which determines the evolution of the foam and the split of the cell. Next cell to be divided into two is a cell chosen randomly, according to a probability proportional to $`R_{loss}|_I`$ or, optionally, a cell with the largest $`R_{loss}|_I`$.

The process of cell division continues until the user-defined maximum number $`N_c`$ of cells is reached. $`N_c`$ includes, in addition to the normal cells called active cell, also all cells that were split, i.e. all parents, grand-parents and great-grand-parents inluding root cell, which we shall call inactive cells. Usually, when we refer to cells, we mean both active and inactive ones. Keeping inactive cells in the record may look like a waste of the memory, but owing to the binary character of the cell split, the loss is only a mere factor of 2. There are many reasons, which make it profitable to keep all inactive cells (including the root cell) in the memory; in particular, as we shall see in Section 2.6, keeping all cells in the record will help us encode all cells in the memory in an economic way, that at higher dimensions we finally gain in terms of total consumption of the CPU memory. Furthermore, for certain quantities that are the integrals over the cell, such as $`R_I`$, we do the following: just after the split, when a new, more precise value of $`R_I`$ is known for the daughter cells – the value of the $`R_I`$ of the parent cell is updated with the sum of the contributions from two daughter cells. This correcting procedure is repeated for all ancestor cells up to the root cell. In this way, the root cell (and any other inactive cell) always keeps track of the actual value of the total $`R_I`$ during the whole foam build-up process. This can be done for any other integral quantity as well, and can be exploited for various purposes.

Since the maximum number of cells $`N_c`$ is defined at the beginning of the foam build-up, all the cell objects and/or other related objects (vertices) are allocated in the computer memory at once, at the very beginning of the cell build-up. On the other hand, the cell objects are organized as a multiply-linked list, with pointers towards parents and daughters. In addition, an array of pointers to all active cells is created at the end of the foam build-up, in order to speed-up the MC generation.

Let us now briefly explain how the geometry of an individual cell is parametrized and stored in the memory. It is relatively easy to parametrize an $`n`$-dimensional hyperrectangle or simplex in a way that does not require much computer memory. An $`n`$-dimensional simplex is fully determined by its $`n+1`$ vertices. Since most of the vertices are shared by the adjacent simplices, the most efficient method is to build an array of all vertices131313In the foam build-up every new vertex is appended at the end of the array. $`\stackrel{}{V}_K,K=1,2,\mathrm{},N_V`$, each of them being an $`n`$-component vector, and to define every simplex as a list of $`n+1`$ vertex indices (integers or pointers) $`K_1,K_2,\mathrm{},K_{n+1}`$. For $`N_c`$ simplical cells resulting from the binary split of a single “root” simplex cell, the number of vertices is $`n+1+N_c`$, because each binary split adds one new vertex. (We include in $`N_c`$ also cells that have got split). The interior points of the simplex are parametrized as follows

$$\stackrel{}{x}=\underset{ip}{\overset{n}{}}\lambda _i(\stackrel{}{V}_{K_i}\stackrel{}{V}_{K_p}),\lambda _i>0,\underset{ip}{}\lambda _i<1,i=1,2,\mathrm{},n,$$
(1)
using basis vectors relative to the $`p`$-th vertex. The above method would be inefficient for $`n`$-dimensional hyperrectangles, because memorizing all $`2^n`$ vertices would require too much memory at higher dimensions. Instead, we use another way of parametrization: each hyperrectangle is defined by the $`n`$-dimensional vector $`\stackrel{}{q}`$ defining the origin of the cell and another vector $`\stackrel{}{h}=(h_1,h_2,\mathrm{},h_n)`$, where each component $`h_i`$ is the length of the hyperrectangle along the $`i`$-th direction. This is even clearer from the explicit parametrization of the interior of the hyperrectangle:

$$x_i=q_i+\lambda _ih_i,0<\lambda _i<1,i=1,2,\mathrm{},k.$$
(2)
For cells with mixed topology, we apply Eq. (2) for $`i=1,2,\mathrm{},k`$ and Eq. (1) for $`i=k+1,k+2,\mathrm{},k+n`$. In Section 2.6 we describe an optional method of storing hyperrectangular cells, in which just two integer numbers are recorded instead of two vectors $`\stackrel{}{q}`$ and $`\stackrel{}{h}`$ (two 2-Byte integers instead of $`2n`$ of 8-Byte floating-point numbers). This method is implemented for the hyperrectangular part of the space only.

### 2.5 Monte Carlo generation

Once the build-up of the cells is finished, the Monte Carlo generation takes place. There is no need for any reorganization of the cells. MC generation can be started immediately. The only thing, that is done at the very end of the foam build-up is the preparation of the list of pointers to all active cells and the array of corresponding $`R_I^{}`$.

The MC point is generated in two steps. First, a cell is chosen with a probability proportional to $`R_I^{}=_{xCell_I}\rho ^{}(x)`$ and next a MC point $`x`$ is chosen with uniform probability inside the cell. The MC weight $`w=\rho ^{}(x)/\rho (x)`$ is associated with the event. For a successful foam of cells, the MC weight is close to 1 and the user may turn weighted events into $`w=1`$ event by means of the rejection method (with the acceptance rate $`w/w_{\mathrm{max}}`$). Foam can do this also for the user. However, the user can sometimes organize the calculation of the $`w`$ and bookkeeping of other parameters of the weight better, in a way that best fits his own aims. This is why the mode of variable weights MC events is also available. The total integral, usually necessary for the proper normalization of the MC sample, is calculated using $`R=R^{}w`$. Foam program provides both the exact value of the $`R^{}`$ and the MC estimate of the integral $`R`$.

### 2.6 Economic use of the computer memory

The actual implementation of the single cell object occupies about 80 Bytes (it could be shrinked to about 40 Bytes, if necessary) of various integer and double precision attributes, plus the dimension-dependent part. In the case of a simplical cell, each new cell adds one $`n`$-component double-precision vector (vertex) and the total memory consumption is therefore of order $`(80+8\times n)`$ Bytes/cell. For $`n=5`$ and 100K cells it is therefore $``$ 15 MBytes of the memory, still an affordable amount. For the hyperrectangle cells we have to count two $`n`$-component double-precision vectors per cell, that is about $`(80+16\times n)`$ Bytes/cell. For the $`10^6`$ cells and $`n=15`$ that would mean $``$ 340 MBytes for the entire foam of cells, which could be annoying. Fortunately, we have found a method of substantially reducing the memory consumption for a hyperrectangular foam. As discussed in Section 3 the geometry of the cell division is fully determined in terms of two integers: one of them is the index of an edge to which the division plane is perpendicular and the other one defines the position of a division plane. The position parameter is a rational number, and only the integer numerator has to be remembered, while the denominator is common to all cells. The above two integers define uniquely the position of the two daughter cells relative to a parent cell. With this method the memory consumption is down to about 80 Bytes/cell, independently of $`n`$ in the present implementation141414In fact, it can be reduced below 40 Bytes/cell, if really necessary.. There is, however, a price to be payed in terms of CPU time. For the generation of the point inside a cell, or even the evaluation of the weight, we need the “absolute” components of $`x`$, that is in the reference frame of the root cell, not relative to vertices of the cell. It is therefore necessary to use a procedure (a method in the class of cells) to construct the absolute position of a given cell “in flight”. This is done by means of tracing all ancestors of a given cell up to the root cell and translating position and size with respect to its parent into absolute ones, step by step, finally relative to the root cell. It is implemented exploiting the organization of cell objects into a linked binary tree. The average number of ancestor cells to be traced back from a given active cell up to the root cell for $`N_c=10^6`$ cells is on the average about $`\mathrm{ln}_2N_c20`$. This may cause about $`20\%`$ increase in the CPU time of the MC generation – an affordable price, if we remember that the MC efficiency improves mainly with the increase of $`N_c`$. In principle, this kind of memory-saving arrangement is also possible for simplical cells; however, in this case the CPU time overhead would be larger, because of the necessity of the full linear transformation at each step, on the way from a given cell up to the root cell. In the case of hyperrectangular cells the transformation is much simpler (and faster); it is the translation and/or dilatation along a single spatial direction at each step.

### 2.7 CPU time saving solution

The final MC efficiency is improved mainly by increasing the number of cells $`N_c`$. The CPU time of the cell build-up is $`Tn\times N_c\times N_{samp}`$, where $`N_{samp}`$ is the number of MC events used in the exploration of each newly created cell. The important practical question is: Can one somehow reduce $`N_{samp}`$ without much loss of final MC efficiency, in order to be able to increase $`N_c`$, within the same CPU time budget?

A simple solution is the following: during the MC exploration of a new cell we continuously monitor an accumulated “number of effective events with $`w=1`$” defined  as $`N_{eff}=(w_i)^2/w_i^2`$, and terminate cell exploration when151515The actual limit of equivalent events per bin is the user-defined parameter, not necessarily equal to 25. $`N_{eff}/n_{bin}>25`$, where $`n_{bin}`$ is the number of bins in each histogram, which is used to estimate the best division direction/edge parameters. This method helps to cut on the total CPU time, because the increase of $`N_{samp}`$ is not wasted for cells in which the distribution $`\rho (x)`$ is already varying very little. At the later stage of the foam evolution this happens quite often. In this method the user may set $`N_{samp}`$ to a very high value and the program will distribute the total CPU time (in terms of $`N_{samp}`$) among all cells economically, giving more CPU time to those cells really need it, i.e. to cells with the stronger variation of $`\rho (x)`$.

### 2.8 Inhibited variables – flat dependence

In some cases the user may not want Foam to intervene into certain variables in the distribution $`\rho (x)`$, simply because there is little or no dependence on them in $`\rho (x)`$. The user may draw, of course, these variables directly from any uniform random number generator. He may, however, find it more convenient to get them from the Foam program. This is easily implemented in Foam: any variable $`x_i`$ may be “inhibited” for the purpose of cell-splitting procedure. In the Foam code it is actually done in such a way that Foam excludes this variable (edge) from the procedure of determining the best binary division of the cell. This provision makes practical sense mainly for the hyperrectangular part of the variable subspace.

In Fig. 3 we show two 2-dimensional foams (250 cells) for the same testing distribution $`\rho (x)`$ (two Gaussian peaks on the diagonal). In one of them (right plot) we have inhibited a split in the first variable, that is for $`x_1`$.

### 2.9 Predefined split points – provision for very narrow peaks

In the practical applications, see Refs. , one may encounter in certain variables extremely narrow spikes (narrow resonances). The Foam exploration algorithm may find it difficult to locate these spikes with the usual method of MC sampling in the cells, at the early stage of the foam build-up. For very narrow spikes, or too low number of requested cells, it may not find them at all! The user usually knows in advance the position of these spikes and the Foam should have a built-in mechanism to exploit this knowledge.

The solution is simple. (It applies for the hyperrectangular subspace of the parameter space only.) The user has a possibility to provide Foam, for each variable, with the list of a number of predefined values: the first splitting positions of the root cell. In the Foam algorithm, it is checked if the list of predefined division points is not empty. If it is the case, then instead of adopting the division parameter from the usual procedure described in Section 3, Foam takes the division parameter from the list, and removes it from the list. In this way the first few division points are taken from the “user-defined menu”, if available, and the next ones are chosen with the usual methods. For narrow spikes this method helps Foam to locate and surround them with as a dense group of cells as necessary.

In Fig. 4 we show an example with two Gaussian peaks in which we requested the Foam program to use the three predefined division points for the $`x_1`$ variable. They are clearly seen as three vertical division lines dividing the entire root cell. In the present case, peaks are not so narrow and there is no real need for a predefined division. The example is just to illustrate the principle of the method.

### 2.10 Mapping of variables

If the structure of the singularities is known and/or Foam is unable to get a reasonable weight distribution for a reasonable number of cells, it is then worth performing an additional change of variables, so that the transformation Jacobian compensates for the singularities, at least partly. In such a case the user subprogram provides Foam with the distribution

$$\rho ^{}(y)=\frac{d\rho }{dy_1\mathrm{}dy_n}=\rho (x_1(y),x_2(y),\mathrm{},_n(y))\left|\frac{x^{(j)}(y)}{y}\right|,$$
(3)
instead of the original $`\rho (x)=d^n\rho /dx^n`$. For each vector $`y`$ generated by Foam, the image vector $`x`$ is well known in the subprogram calculating $`\rho ^{}(y)`$. A mechanism for exporting $`x`$ to the outside world usually has to be provided by the user, because Foam cannot help – it does not know anything about $`x`$; it only knows $`y`$.

Note that in the limiting case of the “ideal mapping” we have

$$\left|\frac{x(y)}{y}\right|\frac{R}{\rho (x)};$$
(4)
consequently, $`\rho ^{}(y)=R`$ and in this case Foam would play merely the role of a provider of the random numbers for $`y`$.

The user of Foam may also need to apply mapping in the case of a “weak” integrable singularity in the distribution $`\rho `$ like $`\mathrm{log}(x)`$ or $`\sqrt{x}`$. Foam can deal with them by brute force, at the expense of a larger number of cells. However, a wiser approach is to apply mapping, in order to remove such singularities from the distribution.

In the next section we shall describe how to combine the mapping method with the multibranching. Such a mixture is well known as the most powerful method of improving the efficiency of the Monte Carlo method.

### 2.11 Provisions for the multibranching

In the following we elaborate on the various methods of implementing multibranching MC method  with the help of Foam.

#### 2.11.1 Single discrete variable

As a warm-up exercise, let us consider the question: Is Foam capable to generate (and sum up) a discrete variable $`i=1,2,\mathrm{},N,`$ according to the (unnormalized) distribution $`r_1,\mathrm{},r_N`$? Of course it can. The simplest way is to define an auxiliary 1-dimensional distribution

$$\rho (x)=r_i,\text{for}\frac{i1}{N}x\frac{i}{N},i=1,2,\mathrm{},N.$$
(5)
The user subprogram providing the above $`\rho (x)`$ is trivial. If plotted, this $`\rho (x)`$ would look like a histogram with $`N`$ equal-width bins. Foam will build up its own grid of cells (intervals), and if we request a high enough number of cells (that is $`N_c>N`$), it will approximate the above $`\rho (x)`$ very well, with its own “histogram-like” distribution $`\rho ^{}(x)`$. However, the Foam approximation will never be ideal, because Foam is not able to detect the exact position of the discontinuities in $`\rho (x)`$. (Nevertheless, this will be a workable solution with a very good weight distribution.) The present Foam algorithm provides for an essential improvement: one may predefine the division points as $`x^{(i)}=i/N`$, $`i=1,\mathrm{},N`$, and set the number of cells to be $`N_cN`$. In such a case Foam will define its cells matching exactly the shape of $`\rho (x)`$. It will generate points with $`w1`$ and provide the exact sum $`R=R^{}=r_i`$, already at the end of the foam build-up. In the MC generation, Foam will randomly generate variable $`x(0,1)`$, which can be translated easily into discrete $`i`$ $`1iN`$.

#### 2.11.2 Discrete and continuous variables

How does the above extend to the case of the distribution $`\rho `$ depending on one discrete variable and the usual $`n`$ continuous variables? For such a distribution $`\rho (y_1,\mathrm{},y_n,i)`$ we define

$$\rho (x_1,x_2,\mathrm{},x_{n+1})=\rho (x_1,\mathrm{},x_n,i),\text{for}\frac{i1}{N}x_{n+1}\frac{i}{N},i=1,2,\mathrm{},N,$$
(6)
in completely analogy with Eq. (5). As previously, we provide for the variable $`x_{n+1}`$ a list of predefined division points $`x_{n+1}^{(i)}=i/N,i=1,2,\mathrm{},N`$, and, of course, we request $`N_c>>N`$. There is still one small problem: Foam may “by mistake” perform an unnecessary cell division for the variable $`x_{n+1}`$, simply due to statistical errors in the “projection histogram” described in Section 3.4.1. This problem is solved in Foam in an elegant way: in addition to providing for $`x_{n+1}`$ predefined division points, the user of Foam may declare $`x_{n+1}`$ as an “inhibited variable” in the sense of Section 2.8. In this case Foam will still split cells according to a list of predefined division points for $`x_{n+1}`$, but will not perform any additional division in this variable! (The translation of the continuous $`x_{n+1}`$ to the discrete $`i`$ is done as usual.) The above method is the basic method of implementation of the “multibranching” (or “multichannel”) MC method using Foam. Let us call it “predefined and inhibited division”, for short PAID method. We shall also describe below how to combine the PAID method with mapping, etc.

In order to appreciate more fully the advantages of PAID, let us consider a more straightforward implementation of the multibranching. In the object-oriented environment one may easily construct $`N`$ instances of the Foam object, each of them for the $`n`$-dimensional function $`\rho (x_1,\mathrm{},x_n,i)`$, initialize them (creating a separate foam of cells) and generate event $`(x,i)`$ with the associated weight $`w_i(x)`$. Index $`i`$ can be chosen according to probability $`p_i=R_i^{}/_jR_j^{}`$, where $`R_i^{}`$ are provided by the $`i`$-th object of the Foam class (at the end of its initialization). The total weight of the event is $`w(x,i)=w_i(x)/p_i`$. Let us call this scenario an “externally organized multibranching”, for short EOM.

Both methods have certain advantages and disadvantages. In PAID the user does not need to organize the optimal/efficient generation of the branching index $`i`$. The root cell is divided into $`N`$ equal-size sub-root cells, which then evolve separately into an independent system of cells, adapting individually to the singularities in the $`i`$-th component of $`\rho `$. Foam adjusts the relative importance of the sub-root cells and their descendants, and finds automatically the optimal number of the division cells in the $`N`$ subfoams within the requested total memory limit. In the EOM scheme these adjustments for the individual branches has to be done by the user. On the other hand, in some cases, the user may want to configure the Foam objects for each branch individually. In the EOM scheme, this can be done for each Foam object separately. In the PAID scheme this cannot be done, because all cells have the same properties, the cell-split algorithm is the same, the cell geometry is common, etc. In most cases, the PAID method will be preferred, because it is easier to organize.

In the following we shall concentrate on the PAID scheme. In this case, the normalization integral is provided by the Foam at the end of the exploration phase, and it includes the sum over a discrete variable:

$$R^{}=\underset{i=1}{\overset{N}{}}\rho _i^{}(x)𝑑x^n.$$
(7)
We also have the usual relation between the average weight and the integral

$$R=\underset{i=1}{\overset{N}{}}\rho (x)_i𝑑x^n=R^{}w.$$
(8)
The above method extends trivially to the case of several discrete variables. As already stressed, the relative probabilities of the discrete components $`p_iR_i^{}=\rho _i^{}(x)𝑑x^n`$ in the MC generation are automatically adjusted by the Foam algorithm, so that the maximum weight or the total invariance is minimized. In the EOM scheme, arranging this in the user program would require an extra programming effort, while in the PAID method this comes for free, by exploiting standard Foam functionality.

#### 2.11.3 Multilayer method

There is an alternative PAID-type method of dealing with the problem of the discrete variable, which generates points according to $`\rho (x_1,\mathrm{},x_n,i),i=1,2,\mathrm{},N`$. It will produce the same distribution but will differ from PAID in the MC efficiency, in terms of the maximum weight or variance. One may simply generate, with the help of Foam, the $`n`$-dimensional auxiliary distribution

$$\overline{\rho }(x_1,\mathrm{},x_n)=\underset{j=1}{\overset{N}{}}\rho (x_1,\mathrm{},x_n,j)$$
(9)
and next, for each generated $`x`$, choose randomly a discrete variable $`i`$ according to the probability

$$p_i(x)=\rho (x,i)/\underset{j=1}{\overset{n}{}}\rho (x,j).$$
(10)
Let us call it PAID, or the multilayer method. This method is slightly less convenient to implement, as is clearly seen for $`n=0`$, where the user effectively has to generate by himself the discrete variable $`i=1,2,\mathrm{},N`$ according to the above probability, by means of creating an inverse cumulative distribution, mapping a random number into $`i`$, etc., while in the standard PAID scenario all this job is done by the Foam program161616The mapping $`x_{n+1}i`$ is a simple arithmetic operation.. Furthermore, in the PAID method each component distribution $`\rho (x,j)`$ may have a “cleaner” structure of the singularities than the sum. Consequently, in the PAID method Foam will probably find it easier to learn the shape of each component distribution than that of the sum in PAID. These two kinds of equivalent multibranching algorithms, PAID and PAID, are described and analysed in Ref. . The PAID method is used in the KKMC event generator of Ref.  to generate index $`i`$ numbering the type (flavour) of the produced quark or lepton.

#### 2.11.4 Multibranching and mapping

However, the most important reason for setting up Foam according to the PAID scenario, with the separate foam build-up for each component distributions $`\rho (x,j)`$, is that for each component one may apply individually adjusted mapping of variables, which makes every component distribution much less singular. The combination of the mapping and multibranching is probably the most powerful known methods of improving MC efficiency . How it can be actually realized with the help of Foam depends on the properties of the distribution $`\rho (x)`$ to be generated. In the case where we have an explicit sum over many components

$$\rho (x_1,\mathrm{},x_n)=\underset{j=1}{\overset{N}{}}\rho (x_1,\mathrm{},x_n,j),$$
(11)
each of the components being positive, with a distinctly different and well known structure of the singularities, we would recommend the use of Foam in the PAID scheme. Knowing the structure of singularities, we may be able to introduce mapping in each component separately, which compensates for these singularities with the Jacobian factor. In such a case Foam is provided with the following distribution:

$$\rho (y_1,\mathrm{},y_n,j)=\rho (x_1^{(j)}(y),\mathrm{},x_n^{(j)}(y),j)\left|\frac{x^{(j)}(y)}{y}\right|,j=1,2,\mathrm{},N,$$
(12)
understanding that the translation of the discrete index $`j`$ into a continuous variable $`y_{n+1}`$, is done in the usual way. The foam of cells is, of course, build-up in the $`y`$-variables, different for each $`j`$-th branch. The user is fully responsible for the proper mapping $`x^{(j)}(y),j=1,2,\mathrm{},N`$, and the calculation of the Jacobian factor in every component (branch). In the user subprogram providing the $`\rho `$-distribution, the variable $`y_{n+1}`$ will first be translated into index $`j`$ and then, depending on the value of $`j`$, a given type of mapping will be applied. For the outside part of the code, the index $`j`$ can be made available, or it may be hidden (erased from the record), depending on the needs of a specific application.

In some cases, however, one does not know any unique way of splitting the $`\rho (x)`$ into well defined positive components as Eq. (12), but one may have a rough idea about the leading singularities. This means, that one is able to construct the distribution

$$\rho (x_1,\mathrm{},x_n)\overline{\rho }(x_1,\mathrm{},x_n)=\underset{j=1}{\overset{N}{}}\overline{\rho }(x_1,\mathrm{},x_n,j),$$
(13)
where $`\overline{\rho }(x,j)`$ have the same type of the leading singularities as $`\rho (x)`$, and one controls the normalization of singularities in $`\rho (x)`$ up to a constant factor; that is for $`x`$ in the neighbourhood of the $`j`$-th singular point, a line or a (hyper)plane, only $`\overline{\rho }(x,j)`$ really matters, that is $`\rho (x)C_j\overline{\rho }(x,j)`$, where $`C_j`$ is not known a priori171717This definition is not very precise; it roughly means that each component is approximately a product of the singular factors and cannot be reduced into the sum of these..

In addition, let us assume that we are able to compensate for the singularities in each $`\overline{\rho }(x,j)`$ exactly by dedicated mapping specific to singularities in the $`j`$-th branch. In other words, the mapping is ideal in each branch:

$$\left|\frac{x^{(j)}(y)}{y}\right|=\frac{\overline{R}_j}{\overline{\rho }(x,j)}.$$
(14)
The above also means, that we know analytically the exact values of the integrals181818In fact, we could normalize $`\overline{\rho }(x,j)`$ to unity, $`\overline{R}_j=1`$, if we wanted.: $`\overline{R}_j=\overline{\rho }(x,j)𝑑x^n`$.

In such a case we may employ the algorithm of Foam successfully by means of defining the “branching ratio”

$$b_j(y_1,\mathrm{},y_n)=\overline{\rho }(y_1,\mathrm{},y_n,j)/\overline{\rho }(y_1,\mathrm{},y_n),b_i(y)=1,$$
(15)
constructing the distribution to be digested by Foam as

$$\rho (y_1,\mathrm{},y_n,j)=b_j(x)\rho (x_1^{(j)}(y),\mathrm{},x_n^{(j)}(y))\left|\frac{x^{(j)}(y)}{y}\right|=\frac{\overline{R}_j\rho \left(x(y)\right)}{_l\overline{\rho }(x(y),l)},$$
(16)
and proceeding as in the PAID scheme described previously.

The role of the function $`b_j(x)`$ is to isolate out from $`\rho (x)`$ “a layer” including just one known type of singularity. In order to see how this method works, let us consider a $`\delta ^{(n)}(xa)`$ shape singularity in the $`j`$-th component (narrow Gaussian peak etc.) of size $`ϵ`$. Then, in the vicinity of that singularity we have $`b_j(x)=1`$ and $`\rho (x)C_j\overline{\rho }(x,j)`$, while further away we have $`\rho (x)𝒪(ϵ^n)`$. The Foam program will, of course, include the $`C_j`$ factor properly in the normalization, and build up the foam of cells everywhere, close to a singularity and far away. It will do it really not in the $`x`$ variables but in the $`y`$ variables. Now, the mapping $`xy`$ (specific to the $`j`$-th branch) will expand the singularity neighbourhood of $`𝒪(ϵ)`$ into a region of size of $`𝒪(1)`$, while the $`y`$-image of the remaining $`x`$-space will be of the size $`𝒪(ϵ)`$. This can be the source of the following drawback: in the small $`y`$-domain of $`𝒪(ϵ)`$ described above, at the positions of the singularities from the other components $`ij`$, $`\rho (y,j)`$ may get narrow spikes or dips of height of $`𝒪(1)`$, such that their integral contribution will be negligible, of $`𝒪(ϵ^n)`$, and this may not look as a problem. However, the Foam algorithm may find it difficult to locate these structures, and this may lead to a small but finite bias of the generated distributions and calculated integrals. This should be kept in mind, and special tests (MC runs with a maximum number of cells, and high MC statistics) should be performed, in order to check that this effect is absent.

The above method is quite similar to that of Ref. . One difference is that in the latter there are several iterations with an aim af adjusting the relative normalization of the components $`\overline{\rho }(x_1,\mathrm{},x_n,j)`$ to $`\rho (x)`$. Our scheme could effectively be regarded as the method of Ref.  with just one iteration; that is the first step being the foam build-up, and the second step (1st iteration) being the MC simulation. One iteration is sufficient in the limit of vanishing overlap of the components $`\overline{\rho }(x_1,\mathrm{},x_n,j)`$ in the entire $`\overline{\rho }(x)`$. While in the method of Ref.  a better adjustment is provided by the next iterations, in Foam the cellular adaptive method provides an extra mileage. One cannot therefore say which one is better in general – it depends on the distribution $`\rho (x)`$.

In fact, in the PAID scheme with the mapping, extra iterations are also possible. It can be done as follows: (a) read $`R_j`$ from all $`N`$ “leading cells” after the foam build-up, (b) rescale $`\overline{\rho }(y_1,\mathrm{},y_n,j)(R_j/\overline{R}_j)\overline{\rho }(y_1,\mathrm{},y_n,j)`$, and (c) repeat the foam build-up191919When programming with Foam it is possible to erase all cells from memory and rebuild them. for the new branching ratios $`b_j(x)`$ in Eq. (16). The above procedure can be repeated. Whether such an iteration is profitable depends on the particular distribution – we expect that in most cases it is not necessary, thanks to the adaptive capabilities of Foam.

Last, but not least, let us also consider the case of a sum of integrals with different dimensionality or, in other words, the distribution in which the number $`n_i`$ of the continuous variables $`x_1,\mathrm{},x_{n_i}`$ depends on a certain discrete “master variable” $`i=1,\mathrm{},N`$ (for example $`n_i=i`$):

$$R=\underset{i=1}{\overset{N}{}}\rho _i(x_1,\mathrm{},x_{n_i}).$$
(17)
Foam can deal with this case too. The simplest solution is to find the maximum dimension $`n_{\mathrm{max}}`$ and add extra dummy variables on which some of $`\rho _i`$ does not depend, such that formally all sub-distributions have the same dimension $`n_{\mathrm{max}}`$. In this way one is back to a situation described earlier, and may apply the PAID method, with or without the additional mapping. The slight drawback of this solution is that in the present implementation of Foam we cannot inhibit the unnecessary cell divisions across the directions of the newly introduced dummy variables – simply because they are not the same in all branches202020A more sophisticated procedure of inhibiting the division could be implemented in Foam, if there is a strong demand.. This is the reason why this kind of problem can, in some cases, be dealt with more efficiently using the EOM scenario, with a separate Foam object for each branch.

For additional practical examples on how to realize multibranching with Foam, see Section 6.

## 3 Cell split algorithm and geometry

As already indicated, our algorithm of the cell split covers two strategies: (A) minimization of the maximum weight $`w_{\mathrm{max}}`$ and (B) minimization of the variance $`\sigma `$, where both $`w_{\mathrm{max}}`$ and $`\sigma `$ are calculated in the Monte Carlo generation, using the MC weight $`w=\rho /\rho ^{}`$. The distribution $`\rho ^{}`$ is the result of the exploration (it is constant over each cell) and is frozen at the end of exploration. The subsequent MC event generation of events is done according to $`\rho ^{}(x)`$. Its integral $`R^{}=\rho ^{}(x)d^nx`$ has to be known exactly before the start of the MC generation. On the other hand, the best estimate of the integral $`R=\rho (x)d^nx`$ is obtained, up to a statistical error, at the end of the MC event generation with the usual relation: $`R=R^{}w_\rho ^{}`$. The average $`\mathrm{}_\rho ^{}`$ is over events generated according to $`\rho ^{}`$.

There is another important ingredient in the algorithm of the cell split: in addition to the auxiliary distributions $`\rho ^{}(x)`$ we also define another distribution $`\rho _{loss}(x)`$ related to the integrand $`\rho (x)`$. The important role of the distribution $`\rho _{loss}(x)`$ is to guide the build-up of the foam of cells; the function $`R_{loss}=\rho _{loss}d^nx`$ is minimized in the process – its value is decreasing step by step, at each cell split. Obviously, both $`\rho _{loss}(x)`$ and $`\rho ^{}(x)`$ are evolving step by step during the foam build-up. Once the division process is finished, the distribution $`\rho _{loss}(x)`$ is not used anymore; MC events are generated with $`\rho ^{}(x)`$. Nevertheless, $`\rho _{loss}(x)`$ is strongly related to the properties of the weight distribution in the MC generation phase.

(A) In the case when our ultimate aim is to minimize $`w_{\mathrm{max}}`$, we define

$$\begin{array}{cc}& \rho ^{}(x)\underset{yCell_I}{\mathrm{max}}\rho (y),\text{for}xCell_I,\hfill \\ & R_{loss}=d^nx[\rho ^{}(x)\rho (x)]=d^nx\rho _{loss}(x).\hfill \end{array}$$
(18)
The distribution $`\rho _{loss}`$ is the difference between the “ceiling distribution” $`\rho ^{}`$ and the actual distribution $`\rho `$ from which it is derived. The rejection rate in the final MC run is just proportional to the integral over the loss distribution $`\rho _{loss}(x)`$ by construction, i.e. the rejection rate $`=R_{loss}/R`$. (This justifies the name “loss”.) The distribution $`\rho _{loss}(x)`$ also has a clear geometrical meaning, see below.

(B) In the case when we do not care so much about the maximum weight and the rejection rate, but we want rather to minimize the ratio of the variance to the average of the weight, $`\sigma /w`$, in the final MC generation, we are led to the following definition:

$$\begin{array}{cc}& \rho ^{}(x)\sqrt{\rho ^2_I},\text{for}xCell_I,\hfill \\ & \rho _{loss}(x)\sqrt{\rho ^2_I}\rho _I,\text{for}xCell_I.\hfill \end{array}$$
(19)
The average $`\mathrm{}_I`$ is over the $`I`$-th cell; see Appendix A for a detailed derivation of the above prescription. The ratio $`\sigma /w`$ in the final MC generation is a monotonous ascending function of $`R_{loss}=\rho _{loss}(x)𝑑x^n`$, see Appendix A. Consequently, minimization of $`R_{loss}`$ is equivalent to minimization of $`\sigma /w`$.

### 3.1 Rules governing the binary split of a cell

The basic rules governing the development of the foam of cells are the following:

* For the next cell to be split we chose a cell with the largest212121In the Foam code there is also an option of choosing the next cell to be split randomly, according to probability proportional to $`R_{loss}`$, instead of a cell with the largest $`R_{loss}`$. $`R_{loss}`$.
* The position/direction of a plane dividing a parent cell into two daughter cells $`\omega \omega ^{}+\omega ^{\prime \prime }`$ is chosen to get the largest possible decrease $`\mathrm{\Delta }R_{loss}=R_{loss}^\omega R_{loss}^\omega ^{}R_{loss}^{\omega ^{\prime \prime }}`$.

How is the split of a given cell into two daughter cells in step (b) done in practice? The method relies upon a small MC exercise within a cell, in which a few hundreds of events are generated with a flat distribution. They are weighted with $`\rho (x)`$ and projected onto $`k`$ (hyperrectangular case) or $`n(n+1)/2`$ (simplical case) edges of the cells. In the mixed case of a cell being the Cartesian product of a $`k`$-dimensional hyperrectangle and an $`n`$-dimensional simplex, there are $`k+n(n+1)/2`$ projections/edges. Resulting histograms are analysed and the best “division edge” and “division hyperplane position” are found – the one for which the estimate (forecast) of the $`\mathrm{\Delta }R_{loss}`$ is the largest. In the actual Foam algorithm, each new-born cell is immediately explored, its $`R_{loss}`$, $`R`$ and $`R^{}`$ are calculated, and the best candidate as to the direction and position of the dividing plane are established and memorized, as the attributes of the cell; see below for details. In this way, every newly created cell is ready for an immediate binary division.

### 3.2 Geometry of the binary split of a cell

In Fig. 5 we show a 3-dimensional cell being a simplex or a hyperrectangle and we visualize the geometry of their split.

Let us first describe the split of the $`n`$-dimensional simplical parent cell into two daughter cells. In this case  we need to know which of $`n(n+1)/2`$ edges defined by any pair of vertices of a given simplex is the best for the split. Suppose that it is an edge defined by a pair of indices $`(i,j),ij,`$ where $`i,j=1,2,\mathrm{},N_V`$, of the vertices $`(\stackrel{}{V}_{K_i},\stackrel{}{V}_{K_j})`$, see Sect. 2.4 for the method of numbering the vertices. A new vertex $`V_{N_V+1}`$ is put somewhere on the line (edge) in between the two vertices:

$$V_{N_V+1}=\lambda \stackrel{}{V}_{K_i}+(1\lambda _{div})\stackrel{}{V}_{K_j},\mathrm{\hspace{0.33em}\hspace{0.33em}0}<\lambda _{div}<1,$$
(20)
where the division parameter $`\lambda _{div}`$ is determined by using an elaborate procedure described later in this section, and the number of vertices is updated $`N_VN_V+1`$. With the new vertex two daughter simplices are formed with the following two lists of vertices (their pointers):

$$\begin{array}{c}\hfill (K_1,K_2,\mathrm{},K_{i1},(N_V+1),K_{i+1},\mathrm{},K_{j1},K_j,K_{j+1},\mathrm{},K_n,K_{n+1}),\\ \hfill (K_1,K_2,\mathrm{},K_{i1},K_i,K_{i+1},\mathrm{},K_{j1},(N_V+1),K_{j+1},\mathrm{},K_n,K_{n+1}).\end{array}$$
(21)
For the $`k`$-dimensional hyperrectangular cell defined with a pair of vectors $`(\stackrel{}{q},\stackrel{}{h})`$ we first decide on the direction of the division hyperplane. Assuming that this hyperplane is perpendicular to the $`i`$-th direction, the two daughter cells (a) and (b) are defined with the two pairs of new vectors as follows:

$$\begin{array}{cc}& \stackrel{}{q}^{(a)}=(q_1,q_2,\mathrm{},q_k),\stackrel{}{h}^{(a)}=(h_1,h_2,\mathrm{},h_{i1},h_i\lambda _{div},h_{i+1},\mathrm{},h_k),\hfill \\ & \stackrel{}{q}^{(b)}=(q_1,\mathrm{},q_{i1},q_i+h_i\lambda _{div},h_{i+1},\mathrm{},q_k),\stackrel{}{h}^{(b)}=(h_1,\mathrm{},h_{i1},h_i(1\lambda _{div}),h_{i+1},\mathrm{},h_k).\hfill \end{array}$$
(22)
The 3-dimensional case of the simplical and hyperrectangular cell split made in this way is illustrated in Fig. 5.

### 3.3 Projecting points into an edge

Before we describe the determination of the division edge and of the division parameter $`\lambda _{div}`$, let us still discuss certain geometric aspects of the Foam algorithm – that is how we project a point $`\stackrel{}{x}`$ inside a cell onto one of the edges of the cell. In the case of a simplex the edges are numbered by the pair of indices $`(i,j)`$, $`i>j`$, which number edges spanned by a pair of vertices222222As explained in Sect. 2.4, the numbering of vertices is done using pointers $`K_i`$ to the elements of the array of vertices. $`(\stackrel{}{V}_{K_i},\stackrel{}{V}_{K_j})`$, while in the case of the hyperrectangle the $`i`$-th edge is spanned by the pair of vectors $`\stackrel{}{q}`$ and $`\stackrel{}{q}_+=(q_1,\mathrm{},q_{i1},q_i+h_i,\mathrm{},q_n)`$. The point $`x`$ inside a cell is projected onto the edge and parametrized using the parameter $`\lambda (0,1)`$. The parameter $`\lambda `$ will be used to define an auxiliary projection $`d\rho /d\lambda `$ for each edge in the following subsection. In particular we have to know how to evaluate $`\lambda `$ in an efficient way. In the case of a simplex cell we have:

$$\stackrel{}{x}_{proj}=\lambda \stackrel{}{V}_{K_i}+(1\lambda _{ij})\stackrel{}{V}_{K_j},\mathrm{\hspace{0.33em}\hspace{0.33em}0}<\lambda _{ij}<1,$$
(23)
where

$$\begin{array}{cc}& \lambda _{ij}(\stackrel{}{x})=\frac{|\mathrm{Det}_i|}{|\mathrm{Det}_i|+|\mathrm{Det}_j|},\hfill \\ & \mathrm{Det}_i=\mathrm{Det}(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_{i1},\stackrel{}{r}_{i+1},\mathrm{},\stackrel{}{r}_n,\stackrel{}{r}_{n+1}),\hfill \\ & \mathrm{Det}_j=\mathrm{Det}(\stackrel{}{r}_1,\mathrm{},\stackrel{}{r}_{j1},\stackrel{}{r}_{j+1},\mathrm{},\stackrel{}{r}_n,\stackrel{}{r}_{n+1}),\stackrel{}{r}_j=\stackrel{}{V}_{K_j}\stackrel{}{x},\hfill \end{array}$$
(24)
and $`\mathrm{Det}(x_1,x_2,\mathrm{},x_n)`$ is determinant. The case of a hyperrectangular cell is much simpler:

$$\lambda _i=(x_iq_i)/h_i.$$
(25)
Obviously, owing to the time-consuming evaluation of the determinants at higher dimensions, the above projection procedure will be much slower for simplices than for hyperrectangles.

In Fig. 6 we illustrate a projection procedure into six edges for 3-dimensional simplex and in Fig. 7 the case of the three edges of the 3-dimensional hyperrectangular cell.

### 3.4 Determination of an optimal division edge and of $`\lambda _{div}`$

Our aim is to find out which division hyperplane, cutting through which edge, provides the best gain of the total integral $`R_{loss}`$, summed over two daughters, as compared to the parent cell. In order to do that, we first analyse all possible division planes, for all edges, and find out the best one, in terms of the gain in $`R_{loss}`$. In other words, we go through all edges ($`k`$ edges for a hyperrectangle and/or $`n(n+1)/2`$ for a simplex) and find, for each edge, the optimal parameter $`\lambda _{div}`$ and the corresponding best gain in $`R_{loss}`$. We then compare between the gains in $`R_{loss}`$ for all edges and define the optimal edge as the one with the best gain in $`R_{loss}`$. The procedure of finding the best $`\lambda _{div}`$ is essentially the same for simplical and hyperrectangular cells; on the other hand, in the algorithm for finding the best $`\lambda _{div}`$ there is a difference between the cases of optimization of the maximum weight and of the variance; see the following discussion.

#### 3.4.1 Optimization of the maximum weight – choosing $`\lambda _{div}`$

Let us consider first the case of finding the best $`\lambda _{div}`$ for $`R_{loss}`$ corresponding to the optimization of the maximum weight. The 1-dimensional case is a good starting point. The cell in this case is just an interval $`(q,q+h)`$ and $`\lambda =(xq)/h`$. In the left-hand part of Fig. 8, we see a histogram with $`N_b`$ bins, made of 1000 events generated inside a cell (interval) using the weight $`w=\rho (x)`$, so that the histogram represents approximately the distribution $`d\rho /d\lambda `$, $`\lambda (0,1)`$. This distribution (histogram) peaks close to the lower edge. The function $`\rho ^{}(x)=\mathrm{max}_{xCell}\rho (x)`$ is constant over the cell and is depicted as an upper horizontal line marked “Parent $`\rho ^{}`$”. The contribution of this particular (parent) cell to $`R_{loss}=_{cell}\rho _{loss}𝑑x=_{cell}(\rho ^{}(x)\rho (x))𝑑x`$ is easily recognized as an area between the line marked “Parent $`\rho ^{}`$” and the histogram line. If we have stopped the exploration at this stage, with this parent cell, then in the MC run points would be generated with the flat “Parent $`\rho ^{}`$” and the weight would be $`w=\rho (x)/\rho `$. Turning weighted events into unweighted ones by accepting $`r<w`$ events and rejecting $`rw`$, where $`0r1`$ is a uniform random number, would correspond to generating points $`(\lambda ,r)`$ within the rectangle below the “Parent $`\rho ^{}`$” line, accepting all points below the histogram line and rejecting all those above, in the area marked “Parent $`R_{loss}`$”. This justifies the subscript “loss”.

The best cell division is found by examining all $`N_b1`$ end-points $`\lambda =q+ih/N_b`$, $`i=1,2,\mathrm{},N_b1`$ of the bins in the histogram, as a possible candidate for the division point (plane in two and more dimensions) between the two daughter cells, and choosing the best one. In the right-hand part of Fig. 8 we have marked such a candidate division point with a star. For a given division point, we determine for two daughter cells the new “ceiling function” $`\rho ^{}`$; in Fig. 8 it is the line marked “New $`\rho ^{}`$”. For each daughter cell we evaluate $`R_{loss}`$. The summary $`R_{loss}`$ for both daughter cells is easily recognized as an area between the line marked “New $`\rho ^{}`$” and the histogram. Of course, we automatically get the new total $`R_{loss}`$ smaller than for the original parent cell! This procedure is repeated for all possible $`j=1,2,\mathrm{},N_b1`$ division points and each time we record the net gain in $`\mathrm{\Delta }_jR_{loss}=R_{loss,parent}R_{loss,daughter1}R_{loss,daughter2}`$. For the actual best division point we chose the division point with the largest gain $`\mathrm{\Delta }_jR_{loss}`$. In Fig. 8 the star marks the best division point.

Let us now consider the 2-dimensional distribution $`\rho (x)`$ depicted in the left part of Fig. 9, which is non-zero within the narrow strip along the four edges of the rectangle. In the simplical mode Foam divides starting $`n`$-dimensional hyperrectangle into $`n!`$ simplices – in this case into two triangles. We concentrate on the division procedure of the lower triangle; see the right part of Fig. 9. In the exploration of this triangular cell we use 1000 MC points and project them onto three edges. The corresponding three histograms are shown in Fig. 10, where the middle histogram represents the projection onto the diagonal – this is why it features two distinct peaks at the ends. In all three plots we have also drawn the curve for $`\rho ^{}(x)`$ for the best hypothetical split. The most promising split in terms of the gain in $`R_{loss}`$ turn outs to be related to the middle plot and is marked in the right part of Fig. 9.

The reader may notice that the $`\rho ^{}(x)`$ in the middle plot of Fig. 10 is not of the type discussed above, because it has two discontinuities instead of one. This is because in Foam we have introduced an important refinement of the algorithm to find an optimal $`\lambda _{div}`$. One may easily notice that the algorithm described above could not correctly locate the drop in the distribution $`d\rho /d\lambda `$ of the middle plot, because there are two equally strong peaks at the ends of this distribution232323 The algorithm would pick up $`\lambda _{div}`$ at a random point between the two peaks.. In the improved version of the algorithm the search of the optimal $`\lambda _{div}`$ uses all pairs of the bin edges $`(\lambda _i,\lambda _j)=(q+ih/N_b,q+jh/N_b)`$, $`0i<jN_b`$. For every pair $`(i,j)`$ a new ceiling function $`\rho ^{}(x)`$ is determined, such that it is unchanged outside the subinterval $`(\lambda _i,\lambda _j)`$ and is “majorizing” the histogram bins inside this subinterval. Once we find out the best pair $`(i,j)`$ in terms of $`R_{loss}`$, then we take either $`\lambda _i`$ or $`\lambda _j`$ as a division point $`\lambda _{div}`$ (at least one of them is not equal to 0 or 1). In the case of two or more peaks in $`d\rho /d\lambda `$ the resulting division point $`\lambda _{div}`$ happens to be close to one of the edges of the gap between the two peaks. This feature prevents the Foam algorithm (at least partly) from placing a new division plane across a void in the multidimensional distribution $`\rho (x)`$. In other words such a void will “repel” the division hyperplanes. In the case of the double peak structure in the middle plot of Fig. 10, the improved algorithm will of course allocate the big value of $`R_{loss}`$ to a new cell (interval), which includes the gap. In the next step of the foam build-up this cell (interval) will have a big chance to be split, and for this split the position of the split point will be located at the second edge of the gap. This is exactly what we need for an efficient foam evolution.

#### 3.4.2 Optimization of the variance – choosing $`\lambda _{div}`$

Let us consider now the case of finding out the best $`\lambda _{div}`$ for $`R_{loss}`$ corresponding to optimization of the variance. The strategy is again to choose $`\lambda _{div}`$ by minimizing $`R_{loss}`$. In Fig. 11 we illustrate our algorithm on the example of the three projections of the triangular cell. The three projections correspond to a triangular cell in two dimensions. (We do not specify $`\rho (x)`$, as it is irrelevant for the purpose of our explanation.) In the upper row of the three plots in Fig. 11 we show as a horizontal line the value of the distribution $`\rho _{loss}=\sqrt{\rho ^2}`$ (it is the same for all three projections). The solid histogram is the distribution $`d\rho /d\lambda `$ and the dashed histogram is the distribution $`d\rho _{loss}/d\lambda `$ calculated bin by bin using $`\sqrt{\rho ^2}`$, treating every bin as a separate cell. The properly normalized difference of the above two distributions $`\sqrt{\rho ^2}\rho `$ is plotted separately as the dashed histograms in the lower row of the three plots in Fig. 11. The horizontal line for the total $`\rho _{loss}`$ is shown once again there. The histogram of $`d\rho _{loss}/d\lambda `$ gives us an idea, where the biggest source of the variance is located and our aim is to “trap” it properly with an intelligent choice of $`\lambda _{div}`$. We follow an algorithm similar to that of the maximum weight minimization, namely we loop over pairs of the bin edges $`(\lambda _i,\lambda _j)=(q+ih/N_b,q+jh/N_b)`$, $`0i<jN_b`$, and for every pair we calculate $`R_{loss}`$ inside the interval $`(\lambda _i,\lambda _j)`$ and outside it. We find out which $`(i,j)`$ provides the biggest gain $`\mathrm{\Delta }_{ij}R_{loss}=R_{loss,parent}R_{loss,Inside}R_{loss,Outside}`$. In the lower row of the plots in Fig. 11 we show as a solid line the distribution of $`R_{loss}`$ for the best pair $`(i,j)`$. Depending on the peak structure, at least one of the division points of the optimal pair $`(\lambda _i,\lambda _j)`$ is different from 0 or 1, and we take this one as $`\lambda _{div}`$. In Fig. 11 the chosen $`\lambda _{div}`$ are marked with black triangles. The above procedure is done for each edge, and the best $`\mathrm{\Delta }_{ij}R_{loss}`$ is used as a guide to define an edge for which the next cell division will be executed. The information about the best edge and the best division point $`\lambda _{div}`$ is recorded in the cell object. As seen in Fig. 11, $`\lambda _{div}`$ tends to fall at the position where $`d\rho _{loss}/d\lambda `$ drops or increases sharply. Note that since the division point is always at the edge of the bin, it is therefore a rational number, $`\lambda _{div}=j/N_b`$. This has interesting consequences, since the number of the bins $`N_b`$ is fixed, it is therefore sufficient to memorize this integer index $`j`$ (2 Bytes) together with the integer index of the division cell edge (also 2 Bytes) as an attribute of the cell, in order to define fully and uniquely the geometry of the division of the cell! See Section 2.6 for more details on how this is exploited to save the computer memory needed to encode the entire foam of cells.

#### 3.4.3 Concluding remarks on the cell division algorithm

The algorithm of the split of the cell is the important and most sophisticated part of the new Foam. Let us therefore add a couple of final remarks:

* The new, much improved procedure of the choice of the division plane is the most significant difference242424 I would like to thank A. Para for a discussion which ignited this new development. with respect to Foam 1.x of Ref. .
* The choice of the edge based on the histograms for each edge makes sense if we use histograms with at least 2–4 bins and at least 100 MC events per cell. This might be a serious limitation for these $`\rho (x)`$, which require a lot of CPU time per function call.

Finally, in Fig. 12 we show examples of the evolution of the foam of cells as the number of the cells grows gradually. The 2-dimensional case is easily visualized and we show it in Fig. 12 for triangular and rectangular cells. In the upper six plots, $`\rho (x)`$ features a circular ridge, in the bottom two plots it is concentrated along the antidiagonal $`x_1+x_2=1`$, the last one corresponding to the $`\rho (x)`$ of Fig 9.

### 3.5 Limitations

We are aware that the present procedure of selecting the next cell for the split and the procedure of defining the division plane, although quite sophisticated, is not a perfect one and has certain shortcomings. Some of them can probably be corrected for, but some of them are inherent to the method. In Fig. 13 we show a surprisingly simple example of a function for which our method of finding a good division point $`\lambda _{div}`$ fails. It fails simply because both projections of $`\rho (x)`$ onto two edges of the rectangle are just flat and our procedure will pick up some $`\lambda _{div}`$ randomly within $`(0,1)`$, while the most economic division point is in the middle $`\lambda _{div}=1/2`$. On the other hand, although the Foam algorithm gets disoriented for the first division, it will recover and correct for the false start in the next divisions rather quickly. It will eliminate the two voids from its area of interest.

Let us notice that the distribution of Fig. 13 violates maximally the “principle of factorizability” $`\rho (\stackrel{}{x})=_1^n\rho _i(x_i)`$, the principle on which the VEGAS family of the programs is built . Contrary to VEGAS the problem with factorizability in Foam is not a general one, but is limited to a single cell and usually goes away after the cell split. Nevertheless, the algorithm of Foam analysing projections on all edges in a single cell is relying on the “principle of factorizability”.

## 4 The Foam code

At present, the C++ version of the Foam code is more advanced than the Fortran77 version. (We do not plan to develop F77 code any further.) In this section we shall therefore describe mainly the C++ code.

The code of the Foam version 1.x was originally written in Fortran77 with popular language extensions, such as long variable names etc. It was already written in an object-oriented style, as much as it had been possible. In particular the main classes TFOAM and TFCELL of the present C++ version were already present under a certain form. The important shortcoming of the F77 version is the lack of dynamic allocation of the memory. Otherwise, it has most of the functionality of the C++ version; see later in this section for a list of limitations.

### 4.1 Description of C++ classes

In Table 1 we list all classes of the Foam package. The main class is TFOAM, which is the MC simulator itself. It is served by the class TFCELL of the cell objects, and three auxiliary classes TFVECT, TFMATRIX and TFPARTITION. The other classes are not related directly to the Foam algorithm – they are utilities used by Foam: random-number generator class TPSEMAR  and the histogramming class TFHIST. The class TFDIST provides a menu of the distributions for testing Foam. In the following we shall describe in more detail the key classes TFOAM and TFCELL.

### 4.2  TFOAM class

TFOAM is the main class. Every new instance of this class (properly initialized) is another independent Foam event generator. In Tables 2 and 3 we provide a full list of the data members of the class TFOAM and their short description. As seen in these tables, we have added the prefix “m\_” to all names of the data members, so that they are visually different, in the code, from the other variables (which is recommended practice in C++ coding).

Generally, one may notice that many data members could be declared (and allocated) as the local variables in the procedures, instead of being data members. For example, vector m\_MCvect, which transpors random numbers out from the random number generator m\_PseMar could be declared locally at every place where m\_PseMar is called. We opted for a more “static” structure of the data, with more than necessary of the data members in the class, at the expense of the human readability of the code, in order to: (a) facilitate the implementation persistency with ROOT, (b) gain in execution speed and (c) facilitate the translation to other languages.

Most of the methods (procedures) of the class TFOAM are listed in Table 4. We omitted in this table “setters” and “getters”, which provide access to some data members, and simple inline functions, such as sqr for squaring a double variable. Data members which are served by the setters and getters are marked in Tables 2 and 3 by the superscripts “$`s`$” or/and “$`g`$”.

Let us now characterize briefly the role of most important methods of the class TFOAM in the Foam algorithm.

#### 4.2.1 Procedures for Foam initialization and foam build-up

The constructor TFOAM(const char\*) is for creating an object of the class TFOAM. Its parameter is the name given by the user to an object. The principal role of this constructor is to initialize data members to their default values – no memory allocation is done at this stage. After resetting all kinds of steering parameters of the Foam to preferred values (using setters) the user calls the Initialize method, which builds up the foam of cells. The two methods InitVertices and InitCells allocate arrays of vertices and cells (pointers) with empty cells. The empty cells are allocated/filled using CellFill. Next comes the procedure Grow which loops over cells, picking up the most promising cell for the split, either randomly using PeekRand or deterministically using Peekmax. The chosen cell is split using Divide. It is, however, the procedure Explore called by Divide (and by InitCells for the root cell) which does the most important job in the Foam build-up: it performs a small MC run for each newly allocated daughter cell. It calculates how profitable the future split of the cell will be and defines the optimal cell division geometry with the help of Carver or Varedu procedures, for maximum weight or variance optimization respectively. All essential results of the exploration are written into the explored cell object. At the very end of the foam build-up, MakeActiveList is invoked to create a list of pointers to all active cells, for the purpose of the quick access during the MC generation. The procedure Explore uses two procedures MakeLambda and MakeAlpha, which generate randomly (uniformly) coordinates of the MC points inside a given cell. The above sequence of the procedure calls is depicted in Fig. 14.

#### 4.2.2 Procedures for MC generation

The MC generation of a single MC event is done by invoking MakeEvent, which chooses randomly a cell with the help of procedure252525Method GenerCell1 exists, but it is not used. GenerCell2 and, next, the internal coordinates of the point within the cell using MakeLambda and/or MakeAlpha. The absolute coordinates of the MC event are calculated and stored in the data member double-precision vector m\_MCvect. The MC weight is calculated using an external procedure that provides the density distribution $`\rho (x)`$, which is represented by the pointer m\_Rho. A class to which the object m\_Rho belongs must inherit from the abstract class TFOAM\_INTEGRAND. The MC event (double-precision vector) and its weight are available through getters GetMCvect and GetMCwt. Note that the variables of the hyperrectangular subspace come first in the m\_MCvect, before variables of the simplical subspace.

The user may alternatively call MCgenerate, which invokes MakeEvent and provides a MC event and its weight, all at the same time.

#### 4.2.3 Procedures for finalization and debugging

The use of the method Finalize is not mandatory. It prints statistics and calculates the estimate of the integral using the average weight from the MC run. The amount of printed information depends on the values of m\_chat. For the normalization of the plots and integrals, the user needs to know the exact value of $`R^{}=\rho ^{}(x)𝑑x`$, which is provided by the method GetIntNorm or Finalize. The actual value of the integrand from the MC series is provided by GetIntegMC. Note that for the convenience of the user GetIntNorm provides $`R^{}`$ or MC estimate of $`R=\rho (x)𝑑x`$, depending on whether the MC run was with weighted or $`w=1`$ events.

Another useful finalization frocedure,

```
GetWtParams(const double eps, double &AveWt, double &WtMax, double &Sigma)

```

provides three parameters that characterize the MC weight distribution: the average weight AveWt, the “intelligent” maximum weight WtMax $`=w_{\mathrm{max}}^\epsilon `$ for a given value of eps $`=\epsilon `$ (see Sect. 6 for its definition) and the variance sigma $`=\sigma `$. In particular, in the case of unweighted events, $`w_{\mathrm{max}}^\epsilon `$ can be used as an input for the next MC run.

The Foam program includes procedure CheckAll, which checks the correctness of the pointers in the doubly linked tree of cells (this can take time for large $`N_c`$). It can sometimes be useful for debugging purposes. Another two methods PrintVertices and PrintCells can be used at any stage of the calculation in order to print the list of all cells and vertices. In the case of the two-dimensions there is a possibility to view the geometry of the cells with a 2-dimensional plot, which is either a LaTeX file produced by LaTexPlot2dim, or a ROOT file produced by RootPlot2dim.

### 4.3 TFCELL class

TFCELL is an important class of objects representing a single cell of the foam. Data members of the class are listed in Table 5.

Most of the methods of the TFCELL class are setters and getters. The non-trivial methods are GetHcub and GetHSize, which calculate the absolute position and size of hyperrectangles in the algorithm of Section 2.6, and MakeVolume, which calculates the Cartesian volume of the cell. In the simplical subspace, the volume is a determinant of a square matrix of the class TFMATRIX.

### 4.4 Persistency with the help of ROOT

The C++ language does not provide any built-in mechanism of the persistency of the classes. For this purpose we use the ROOT package , with the help of its “automatic streamers”. ROOT is a useful C++ library for histogramming, organizing large databases of identical objects of the type used in high energy physics experiments. It also provides an efficient input/output, with compressing capabilities.

Providing full persistency for C++ classes preserving all the structures of the pointers is not a trivial task in general. ROOT can do it, even for pointers. For this, our code had to be reorganized slightly: puting additional directives for ROOT in the comment lines, remooving static variable, explicit integer indices instead of pointers in some places. As a whole, this solution is relatively simple, and works correctly. In Tables 2 and 3 we have at the beginning of the description certain characteristic marks which are directives for the persistency mechanism of ROOT, see the ROOT manual  for more details.

One has to remember, when reading a TFOAM class object from the disk, that the method LinkCells() has to be invoked in order to fully reconstruct all pointers in the doubly linked tree of cells. Moreover, any object of the class TFOAM restored from the disk file will have its internal object for the random number generator and distribution function. There is a method that provides access (pointer) to these objects, if necessary. The relevant fragment of the code may look as follows:

```
 TPSEMAR *RNGen= FoamX->GetPseRan();        //get pointer of RN generator
 TFDISTR *RHO  = (TFDISTR*)FoamX->GetRho(); //get pointer of distribution

```

It might be useful if, for instance, one wants to reinitialize the random number generator used by the TFOAM class object, which has been read from the disk file, with a new “random seed”.

Foam can be used with or without ROOT. In the code all parts of the code dependent on ROOT enclosed in the pair of preprocessor commands #ifdef ROOT\_DEF ... #endif, where the ROOT\_DEF environmental variable is defined centrally in the header file ROOT\_DEF.h. Eliminating ROOT requires removing this variable and modifying makefile accordingly (the TFHST class has to be linked). The version without ROOT does not feature persistency, and is employing its own simple histogramming class TFHST instead of the ROOT class TH1D. ROOT also helps to create documentation of the Foam in the html format. We recommend a version tied up with ROOT.

### 4.5 Fortran77 version and its limitations

We also provide users with the Fortran77 versions of Foam. There are two of them at the development level 2.02 (May 2001) of the algorithm. In the first one cells can be simplical, hyperrectangular and the Cartesian product of the two. This version is limited to dimension five for the simplical subspace and not very useful for large dimensions ($`n5`$) in the hyperrectangular space, because it does not feature the memory-saving algorithm of Section 2.6. Another version, named MCell (standing for Mega-Cell), features only hyperrectangular cells; on the other hand, it includes the memory-saving algorithm described in Section 2.6. We recommend the reader to use the version MCell.

Neither of these versions can have dynamic memory allocation; they have a maximum dimensions of the integration/simulation subspaces (simplical and hyperrectangular) hard-coded in the source code. Any change of these maximum dimensions requires recompilation of the code.

Present versions in Fortran77 are substantially improved with respect to the original version of Ref. . For the option of minimizing the maximum weight, they have exactly the same algorithm (of the cell split) as the C++ version. They feature, however, an older, more primitive version of the algorithm of finding the best cell division for the variance reduction.

The structure of the programs, naming of procedures and variables, configuration parameters and their meaning are very similar in the F77 and C++ versions. Some differences in the usage will be indicated in the next subsections.

### 4.6 Future development

In the following we indicate some of the possible future developments of the Foam package. As already indicated we do not plan to develop the Fortran77 version any further. On the other hand, it would be interesting to upgrade the existing Foam version 1.x in the JAVA language to the level of the present version 2.x.

As for the C++ version, it would be a logical development to inherit the class of pseudo-random number generators TPSEMAR from a unique abstract class, and in this way to define a universal interface for a library of the number generators. We intend to collect a library of a few random number generators with a universal interface (or find one) for the use in Foam and applications based on it.

Concerning the version of Foam adapted to parallel processing, as in Refs. , we do not have plans in this direction in the immediate future. Here, we have of course in mind the use of the true CPU parallelism in the foam of cells build-up. One has to remember that in the high energy physics applications, which are our main objective, the foam of cell build-up will always be a tiny fraction of the total CPU time. The main fraction will be the subsequent MC simulation in which, as the vast experience with PC-farms in CERN and FNAL shows, one may organize the MC simulation with a low level of parallelism, with many simulators started with different random seeds running in parallel but not communicating – another specialized job combines all the results at the very end of the run. However, the first practical examples of true parallelism in the massive MC simulation for the purpose of the high energy physics experiments has appeared recently  and is used in the LEP2 data analysis.

The main algorithm of Foam is quite stable at this stage, and the main emphasis in its future development will be on making it more user friendly and better adapted to the use as a part of bigger MC projects. In particular its provisions for multibranching will become more sophisticated, as more feedback comes from real-life applications.

## 5 Usage of Foam

In the following we provide basic information on how to use the Foam program in the user application.

### 5.1 Foam distribution directory of the C++ version

The Foam package is distributed together with the demonstration main programs and some utilities in the form of about 20 files in a single UNIX directory FOAM-export-v2.05. Demonstration runs can be executed using standard make  commands as follows:

```
       make Demo-run
       make DemoPers
       make Demo-map
       make DemoNR-run

```

The essential fragments of the output from running command “make Demo-run” are shown in Appendix B. The compilation and linking procedure is encoded in the file Makefile, and it has to be checked by the user whether it conforms to the local operating system. In particular, if ROOT is used, then certain paths and environmental variables in Makefile have to be adjusted. The use of ROOT is decided by the presence of the variable ROOT\_DEF in the file ROOT\_DEF.h. Without ROOT, the user should execute command “make DemoNR-run”.

The essential part of Foam, that is class TFOAM, TFCELL and a few auxiliary classes, is located in the files TFOAM.cxx and TFOAM.h. This is the “core” of the Foam source. The source code of the other utility classes TFHST, TFMAXWT, TPSEMAR and TFDISTR are in separate files. The main programs are in files Demo.cxx and DemoPers.cxx. They should serve as a useful template for the user’s own application based on Foam.

There is also one Fortran77 source code, circe2.f262626We thank Thorsten Ohl for providing us with a preliminary version of this code ., which contains a testing density distribution (2-dimensional beamstrahlung spectrum  of the electron linear collider ) used in TFDISTR. The Makefile provides, therefore, also a useful example of linking C++ and F77 codes.

There are also two output files output-Demo.linux and output-DemoNR.linux, which the reader may use to check whether he is able to reproduce these benchmark output results.

### 5.2 Simple example of an application

A very simple example of the use of Foam may look as follows:

```
//   *** Initialization ***
double MCwt;
TFDISTR *Density1 = new TFDISTR(FunType); // Create integrand function
TPSEMAR *PseRan   = new TPSEMAR();        // Create random numb. generator
TFOAM   *FoamX    = new TFOAM("FoamX");   // Create Simulator
FoamX->SetkDim( 3);                       // Set dimension, h-rect.
FoamX->Initialize(PseRan,  Density1 );    // Initialize simulator
//   *** MC Generation ***
TFHST *hst_Wt = new TFHST(0.0,1.25, 25);  // Create weight histogram
double *MCvect =new double[3];            // Monte Carlo event
for(long loop=0; loop<1000000; loop++){
   MCwt = FoamX->MCgenerate(double *MCvect); // Generate MC event
   hst_Wt->Fill(MCwt,1.0);                   // Fill weight histogram
}
//   *** Finalization ***
double IntNorm, Errel;
FoamX->Finalize( IntNorm, Errel);  // Print statistics, get normalization
double MCresult, MCerror, AveWt, WtMax, Sigma;
FoamX->GetIntegMC( MCresult, MCerror);        // get MC integral
double eps = 0.0005;
FoamX->GetWtParams(eps, AveWt, WtMax, Sigma); // get MC wt parameters
hst_Wt->Print();                   // Print weight histogram

```

The user normally provides his own density distribution function belonging to the class that has to inherit from the following abstract class:

```
class TFOAM_INTEGRAND{ // Abstract class of distributios for Foam
 public:
  TFOAM_INTEGRAND() { };
  virtual ~TFOAM_INTEGRAND() { };
  virtual double Density(int ndim, double*) = 0;
};

```

In the above example the distribution \*Density1 belongs to the class TFDISTR, which is provided in the Foam distribution directory.

### 5.3 Configuring Foam

Foam has fourteen principal configuration parameters plus parameters inhibiting and/or predefining the division geometry in the cell split.

#### 5.3.1 Principal configuration parameters

All of the principal parameters listed in Table 6 are set to meaningful default values, hence the beginner may stay ignorant of their role for some time, and learn gradually how to exploit them in order to improve the efficiency of Foam in the actual application. All these parameters are data members of the TFOAM class, see Table. 2. If the user wants to redefine all of them, then the relevant piece of code will look as follows:

```
  FoamX->SetnDim(     nDim);
  FoamX->SetkDim(     kDim);
  FoamX->SetnCells(   nCells);
  FoamX->SetnSampl(   nSampl);
  FoamX->SetnBin(     nBin);
  FoamX->SetOptRej(   OptRej);
  FoamX->SetOptDrive( OptDrive);
  FoamX->SetOptPeek(  OptPeek);
  FoamX->SetOptEdge(  OptEdge);
  FoamX->SetOptOrd(   OptOrd);
  FoamX->SetOptMCell( OptMCell);
  FoamX->SetEvPerBin( EvPerBin);
  FoamX->SetMaxWtRej( MaxWtRej);
  FoamX->SetChat(     Chat);

```

In practical applications one will redefine only some of them. The minimum requirement is that the user sets a non-zero value of nDim or kDim such that the total dimension nDim+kDim is a non-zero positive integer.

#### 5.3.2 Inhibiting cell division in certain directions

If the user of Foam decides to inhibit the division in certain variable in the hyperrectangular subspace, this can be done with the method SetInhiDiv(int iDim, int InhiDiv) of the class TFOAM, where iDim is the dimension index for which inhibition is done and InhiDiv is the inhibition tag. This method should be used before invoking Initialize, after setting nDim and/or kDim. The relevant code may look as follows:

```
  FoamX->SetInhiDiv(0, 1); //Inhibit division of x_1
  FoamX->SetInhiDiv(1, 1); //Inhibit division of x_2

```

The allowed values are InhiDiv=0,1 and the default value is InhiDiv=0. Note that the numbering of dimensions with iDim starts from zero and variables of the hyperrectangular subspace always come first, before the simplical ones.

#### 5.3.3 Setting predefined cell division geometry

The user may predefine divisions of the root cell in certain variables in the hyperrectangular subspace using the method SetXdivPRD(int iDim, int len, double xDiv\[\]). The relevant piece of the user code may look as follows:

```
  double xDiv[3];
  xDiv[0]=0.30; xDiv[1]=0.40; xDiv[2]=0.65;
  FoamX->SetXdivPRD(0, 3, xDiv);

```

Again, this should be done before invoking Initialize, after setting nDim and/or kDim.

### 5.4 Persistency

Persistency of the Foam classes is arranged using “default streamers” of the ROOT  package. Writing a TFOAM class object into the disk file rmain.root can be done with the single Write as follows:

```
  TFile RootFile("rmain.root","RECREATE","histograms");
  ...
  FoamX->Write("FoamX"); //Writing Foam on the disk
  ...
  RootFile.Write();
  RootFile.Close();

```

The instruction FoamX-\>Write("FoamX") can be put at any place of the code after the instruction FoamX-\>Initialize(...), see example of the user code shown in Section 5.2.

Next, in another program, the TFOAM class object can be read from the disk file rmain.root as follows:

```
  TFile fileA("rmain.root");  // connect disk file
  fileA.cd();
  fileA.ls();                        // optional printout
  fileA.Map();                       // optional printout
  fileA.ShowStreamerInfo();          // optional printout
  fileA.GetListOfKeys()->Print();    // optional printout
  TFOAM  *FoamX = (TFOAM*)fileA.Get("FoamX"); // find object
  FoamX->LinkCells();  // restore pointers of the binary tree of cells
  FoamX->CheckAll(1);                // optional x-check of pointers

```

and at this point the FoamX object is ready to generate MC events, as in the MC generation part of the code shown in Section 5.2.

### 5.5 Fortran77 versions

The distribution directory FoamF77-2.02-export contains the README file, two demonstration main programs DemoFoam.f and DemoMCell.f to be compiled and run with the help of commands

```
        make DemoFoam
        make DemoMCell

```

encoded in the Makefile. The outputs from the above runs can be compared with the benchmark outputs output-DemoFoam-linux and output-DemoMCell-linux.

The basic Foam source files are: FoamA.f with header file FoamA.h and MCellA.f with header file MCellA.h.

For the description of the input (configuration) parameters, see comments in FoamA.f and MCellA.f respectively. The names of the configuration variables are the same as in the C++ version, except nCells, which is here renamed to nBuff. Their values and the meaning are the same.

Demonstration main programs DemoFoam.f and DemoMCell.f can serve as templates for the user application programs.

The testing main program uses the histogramming package GLK of the KKMC program , which the user may replace with any other histogramming package.

## 6 Numerical studies and example applications

In the following subsection we examine MC efficiency of the Foam in a series of numerical exercises. In some of them we shall also show examples of the Foam application with the distributions relevant to everyday practice in high energy physics.

### 6.1 Dependence of the Foam efficiency on the configuration parameters

As a first numerical exercise, we examine the dependence of the Foam efficiency on the most important (input) configuration parameters, including the dimension of the space.

In Table 7 we collect results from many MC runs for various dimensions, numbers of cells and numbers of MC events in single cell exploration, varying also the type of the cells. We always use the same test distribution $`\rho _{camel}(x)`$ as in Ref. , which features two relatively narrow gaussian peaks placed on the diagonal. The 2-dimensional version of this distribution is shown in Fig. 15. We have used the non-default values nBin=4 and EvPerBin=50 and the default values for the other configuration parameters. In this table all tests were done for the default option of reduction of the maximum weight, OptDrive=2 (for results with OptDrive=1, see next table). The efficiency $`w_{\mathrm{max}}^\epsilon /w`$ is calculated using maximum weight $`w_{\mathrm{max}}^\epsilon `$ defined as in272727The $`\epsilon `$-dependent maximum weight is defined such that events with $`w>w_{\mathrm{max}}^\epsilon `$ contribute an $`\epsilon `$-fraction to the total integral. It is numerically more stable in the numerical evaluation than the one defined as the largest weight in the MC run. Ref.  for $`\epsilon =0.0005`$. The maximum weight $`w_{\mathrm{max}}^\epsilon `$ is calculated with the help of the small auxiliary class TFMAXWT.

In Table 7 the efficiency of the MC run measured in terms of $`w_{\mathrm{max}}/w`$ and $`\sigma /w`$. The value of the integral $`R\pm \mathrm{\Delta }_{statist.}R`$, shown in the last four columns, was obtained from the MC run in which the total number of MC events was $`2\times 10^6`$. The value of the integral $`R`$ is well known, it is equal to 1, within $`10^5`$.

The following observations based, on the results of Table 7, can be made:

* Looking at the results, for total dimension $`n=4`$ we see that the hyperrectangular cells clearly provide better MC efficiency than simplical ones. All other results are for hyperrectangular cells.
* All of the results support the observation that the MC efficiency depends critically on the number of cells. In particular, see results for $`n=6`$, the increase of nSamp (No. of MC events in cell exploration) beyond a certain value does not help at all.
* In the case of the very inefficient Foam, see $`n=12`$ with $`\sigma /w6`$, the estimate of the MC statistical error can be misleading. We see an indication that one should not trust runs with $`\sigma /w>3`$.
* For this particular testing function the dimension $`n=12`$ requires a minimum of 400k cells and resulting MC efficiency of order of 1% is barely acceptable282828However, we still get the correct value of the integral within 0.2%..

In Table 8 we repeat the exercise of Table 7 for the option of the variance reduction OptDrive=1 at four dimensions. As compared to Table 7 we see a net improvement in the variance and deterioration of the $`w_{\mathrm{max}}^\epsilon `$. This agrees with the expectations.

### 6.2 Comparison with Foam 1.x and classic VEGAS

In Table 9 we update the comparison of Foam and VEGAS of Ref. , adding results for the new hyperrectangular option. The simplical results are now clearly improved with respect to Ref. , because of the better cell division algorithm. Generally, the hyperrectangular cell mode provides as good an efficiency as the simplical one. However, one should keep in mind that Foam with hyperrectangular cells is a factor of 2 or more faster in the execution.

### 6.3 Example of sharply peaked distribution

In Table 10 we examine the dependence of the MC efficiency/error on the various input configuration parameters of Foam. All these numerical results are for the distribution

$$\rho _g(x)=\frac{\mu x_2}{(x_1+x_21)^2+\mu ^2},$$
(26)
which, for $`\mu =10^6`$, has a very sharp ridge across the diagonal $`x_1+x_2=1`$. This distribution is taken from Refs. and is related to the photon distribution at high energy electron–positron colliders.

What can we learn from the results in Table 10? First of all, in the first line, we see a spectacular failure of Foam with rectangular cells292929Setting the number of the MC events in a single rectangle exploration to $`10^4`$ cures the problem partly.. The value of the integral is wrong by a factor of 2 and the statistical error is underestimated. This illustrates the problem of the the lack of “angular mobility” of the rectangular cells discussed in Section 2.3; rectangles are unable to align with the singularity along the diagonal.

This we illustrate in the left plot of Fig. 16, for rectangular 1000 cells, where we clearly see “blind spots”. In the right-hand plot the triangular cells in the foam are clearly aligned with the diagonal ridge. In the rows 2 and 3 of Table 10 we see the reasonable numerical results for the triangular foam, which are for the maximum weight reduction and variance reduction options respectively; the other configuration parameters are rather close to the default ones. In rows 4 and 5 we are playing with the increase of the cell number and in the rows 6 and 7 with the number of MC events used in the cell exploration. Finally in rows 8 and 9 we change the binning of the histograms used in the MC cell exploration. As we see, the most profitable in terms of MC efficiency/precision is the increase of the number of cells however, adjusting other parameters can also help. In all cases we display variable nCalls, the number of calls of the density distribution function during the foam build-up. In the best result of line 5 with 10000 triangular cells we have obtained a 5-digit precision for about $`10^7`$ function calls303030In Ref.  the same precision for the same function was attained for about $`10^8`$ function calls..

Summarizing, we see from the above exercise that the user of Foam has a possibility to adjust several configuration parameters, so that the MC efficiency for a given distribution is improved quite significantly.

### 6.4 Decay of the $`\tau `$ lepton into 3 pions

In Table 11 we collect numerical results for an example of the application of Foam to the very practical problem of the MC simulation of the decay $`\tau \nu \pi ^{}\pi ^+\pi ^{}`$, according to the distribution (matrix element squared) taken from the MC program TAUOLA .

The amplitude of the decay process contains two distinct parts due to two Feynman diagrams, see Fig. 17, which have peaks due to the $`a_1`$ resonance and the $`\rho `$ resonance. There are two peaks due to $`\rho `$ resonances that partly overlap in the Lorentz invariant phase-space, such that the actual shape of the differential distribution is rather complicated. We took for this exercise subroutine DPHTRE of TAUOLA, in which nine random numbers313131Including two random numbers of the subroutine SPHERA and two (Euler) angles corresponding to the overall rotation of the entire event. are replaced by the nine variables of Foam. The 4-particle phase-space is 8-dimensional. The ninth variable is due to two branches in the phase-space parametrization of TAUOLA, and in case of Foam as well (the method is similar to that of Section 2.11). In cases of “no mapping and no multibranching” we are back to eight dimensions. The variables $`x_1`$ and $`x_2`$ of Foam represent (up to a linear transformation) the two effective masses of the $`3\pi `$ and $`2\pi `$ systems. The next four variables $`x_i,i=3,4,5,6`$ are the polar variables $`\mathrm{cos}\theta `$ and $`\varphi `$ of the pions in the rest frame of the $`3\pi `$ and $`2\pi `$ systems – this is a completely standard phase-space parametrization, see Ref. , and also Ref. . The variables $`x_7`$ and $`x_8`$ are reserved for the overall rotation, and the last one, $`x_9`$, is mapped into branch index $`j=1,2`$, see Section 2.11. Variables $`x_7`$ and $`x_8`$ are inhibited (no cell-split in them), because the distribution does not depend on them. Variable $`x_9`$ (if present) has a predefined division value equal to $`0.5`$ and is inhibited for the division (see Section 2.11). In Fig. 18 we show the decay distribution as a function of $`x_1`$ and $`x_2`$ (without mapping), fixing the other variables to some values. The distribution is clearly a non-trivial one.

In Table 11 we show the MC efficiency of Foam for a gradually increasing number of hyperrectangular cells, with the mapping compensating for the Breit–Wigner peaks of the $`a_1`$ and $`\rho `$ resonances (as in Ref. ) and without. One example with simplical cells is also included.

The most striking result in Table 11 is the comparison of lines (a2) and (b1): the Foam algorithm with only 20 cells is performing equally well as the doubly-branched mapping compensating for the resonance peaks of the $`a_1`$ and $`\rho `$. When going to higher number of cells, the MC efficiency in the cases with and without mapping becomes almost the same. This is expected, since Foam also does the mapping compensating for the resonances on its own. From row (c0) we see also that the simplical mode of Foam is clearly underperforming. We think that Foam configured with 1000 cells, see row (c1) in Table 11, is an economic solution to this particular MC problem323232The net profit with respect to TAUOLA  would be a three times faster program, and more importantly, a significantly simpler code.. (Mapping due to resonances is not really necessary in this case.)

### 6.5 Beamstrahlung spectrum

In Fig. 19 we show the foam of rectangular cells for the 2-dimensional beamstrahlung spectrum $`D(z_1,z_2)`$ of the electron–positron collider at 500 GeV, encoded in the program circe2 of Refs. . It should be stressed that this spectrum is not known analytically but rather through a numerical fit to results of the machine simulation or (in the future) from an experiment. In order to avoid (integrable) infinite singularities at $`z_i=1`$ in $`D(z_1,z_2)`$ we use the variables $`t_i=(1z_i)^{0.1},i=1,2`$ in Fig. 19. For this exercise we used foam of 500 cells getting the MC efficiency $`\sigma /w=0.41`$ and $`w/w_{\mathrm{max}}^\epsilon =0.64`$ (for $`\epsilon =0.0005`$); enough for practical application (can be easily improved by adding more cells).

Generating $`D(z_1,z_2)`$ is not really such a very much important and difficult problem. A more interesting problem is to generate the distribution $`D(z_1,z_2)\sigma (sz_1z_2)`$, where $`\sigma (s)`$ is the cross section of some physics process, which may have a strong singularity of its own, such as a resonance or threshold factor. Such a problem was already treated with the help of the Foam program in KKMC  event generator and the study of Ref..

## 7 Conclusions

The author hopes that this new adaptive tool for constructing efficient MC programs will find its way to many applications in high energy physics and beyond. The main points on the new Foam algorithm and the program are the following:

* Foam is a versatile adaptive, general-purpose Monte Carlo simulator.
* The Foam algorithm is based on the cellular division of the integration domain.
* The geometry of the “foam of cells” is rather simple, cells of simplical and/or hyperrectangular shape are constructed in the process of a binary split.
* It works in principle for arbitrary distribution – no assumption of factorizability as in VEGAS of Ref..
* Foam is reducing maximum weight of the weight distribution, it can therefore provide unweighted events. The variance reduction, useful for the integration and generating weighted events, is also available.
* Encoding cells with memory-efficient methods allows up to $`10^6`$ cells to be built in the computer memory of a typical desktop computer.
* The rules for choosing the next cell for division and the division geometry are based on rather sophisticated analysis of the projection on the cell edges. This costs CPU time, which becomes the main barrier towards higher MC efficiency.
* Foam can deal efficiently with rather strongly peaked distributions, up to relatively high dimensions, $`12`$ dimensions, with today’s desktop computers.

Acknowledgements

I would like to thank T. Ohl, A. Para, W. Płaczek, E. Richter-Wa̧s, F. Tkachov and Z. Wa̧s for interesting and useful discussions and R. Brun for help in getting persistency using ROOT. The support and warm hospitality of the CERN Theory Division and DESY Zeuthen is kindly acknowledged. The useful assistance of the Parasoft Company in debugging the C++ code (with the Insure++ tool) is also acknowledged.

## Appendix A Variance optimization

Suppose we have already constructed the cells $`\omega _1,\omega _1,\mathrm{},\omega _N`$ and within each cell we have defined the function $`\rho ^{}(x)`$ constant over the cell, $`\rho ^{}(x)=\rho _I^{}`$. The integral $`R_I^{}=_{\omega _i}\rho ^{}(x)𝑑x^n=\rho _I^{}V_I`$ is known, because the volume of the cell $`V_I`$ is known. The function $`\rho (x)`$ is not, in general, constant over the cell and the weight $`w=\rho (x)/\rho ^{}(x)`$ is used to determine its integral in the usual way: $`R=R^{}w_\rho ^{}=_Iw_{\rho _I^{}}R_I`$, where the average

$$a_{\rho _I^{}}=\frac{1}{R_I^{}}_{\omega _I}\rho _I^{}(x)a(x)𝑑x^n$$

is defined for the $`I`$-th cell alone.

The question is now the following: preserving the geometry of the cells, can we get smaller variance by simply changing the probabilities of the generation of the cells? Rescaling $`\rho _I^{}\rho _I^{\prime \prime }=\lambda _I\rho _I^{}`$ does not affect the integral $`R`$, because the change of the normalization of $`R^{}`$ and $`w_\rho ^{}`$ is cancelling. It is convenient to assume that the above rescaling preserves $`R^{}=R^{\prime \prime }`$ and the total average weight $`w_\rho ^{}=w_{\rho ^{\prime \prime }}`$, that is $`\lambda _I`$ obey the constraint $`_IR_I^{}=_IR_I^{}\lambda _I=`$ const. With the above constraint in mind, we now ask: for what values of $`\lambda _I`$ is the dispersion of the total weight $`\sigma ^2=w^2_{\rho ^{\prime \prime }}w_{\rho ^{\prime \prime }}^2`$ minimal?333333 This problem was, of course, often considered in the past, see for instance Refs. . We outline here the solution for the sake of completeness and convenience of the reader. Since by construction $`w_{\rho ^{\prime \prime }}`$ is independent of $`\lambda _I`$, we may look only for a minimum of $`w^2_{\rho ^{\prime \prime }}`$. With the standard methods we get a (local) minimum condition:

$$\frac{}{\lambda _I}\left\{w^2_{\rho ^{\prime \prime }}+\mathrm{\Lambda }R^{\prime \prime }\right\}=\frac{}{\lambda _I}\left\{\frac{1}{R^{}}\underset{I}{}_{\omega _I}\left(\frac{\rho (x)}{\rho _I^{}\lambda _I}\right)^2\rho _I^{}\lambda _I𝑑x^n+\mathrm{\Lambda }\left(\underset{I}{}R_I^{}\lambda _I\right)\right\}=0,$$
(27)
where $`\mathrm{\Lambda }`$ is the Lagrange multiplier. The solution of the minimum condition

$$\frac{R_I^{}}{R^{}}w^2_{\rho _I^{}}\frac{1}{\lambda _I^2}\mathrm{\Lambda }R_I^{}=0$$
(28)
is simply $`\lambda _I\mathrm{const}\times \sqrt{w^2_{\rho _I^{}}}`$, or more precisely

$$\lambda _I=\frac{\sqrt{w^2_{\rho _I^{}}}}{_J\frac{R_J^{}}{R^{}}\sqrt{w^2_{\rho _J^{}}}}.$$
(29)
The value of $`w^2_{\mathrm{min}}w^2_{\rho ^{\prime \prime }}`$ calculated at the minimum, that is for $`\lambda _I`$ of Eq. (29), is also rather simple

$$w^2_{\mathrm{min}}=\underset{I}{}\frac{R_I^{}}{R^{}}w^2_{\rho _I^{}}\frac{1}{\lambda _I^2}=\left(\underset{J}{}\frac{R_J^{}}{R^{}}\sqrt{w^2_{\rho _J^{}}}\right)^2=\left(\sqrt{w^2_{\rho _J^{}}}_\rho ^{}\right)^2.$$
(30)
Let us note that the functional $`\sqrt{w^2_{\mathrm{min}}}`$, which we intend to minimize in the process of the evolution of the foam (cell split) is a simple sum of contributions from all cells. Consequently, when working out details of the split of a given cell we may calculate the gain in terms of the total dispersion independently of other cells, see also below. This very convenient feature is exploited in the algorithm of the foam build-up.

Adopting (temporarily) the following normalization conventions: $`\rho _I^{}=1`$, $`w=\rho (x)`$ and $`R_I^{}=V_I`$, $`\lambda _I\sqrt{\rho ^2}_I`$, we get

$$\rho _I^{\prime \prime }=\sqrt{\rho ^2_I},R_I^{\prime \prime }=V_I\sqrt{\rho ^2_I},$$
(31)
where the average $`\mathrm{}_I`$ is understood as defined for points uniformly distributed within the $`I`$-th cell.

In the Foam algorithm we do not go, of course, through a judicious adjustment of the relative importance of the cells, which minimizes the variance as described above; but instead, we simply declare that the distribution $`\rho ^{}(x)`$ is defined by the optimum solution of Eq. (31). Once we have done that, for the MC weight $`w`$ defined with respect to such a new $`\rho ^{}(x)`$, we find out that $`w^2_I=1`$ for each cell and, also for all cells, $`w^2=1`$. What is then minimized in the process of the cell division is the ratio of the variance to the average weight343434 We exploit here the relation $`w=R/R^{}`$, and we should keep in mind that $`R`$ is constant during the variance minimization.

$$\frac{\sigma ^2}{w^2}=\left(\frac{R^{}}{R}\right)^21.$$
(32)
This quantity is not so convenient to optimize in the process of the cell split (foam evolution) and we rather chose to minimize a closely related ‘linearized” quantity

$$R_{loss}=R\left(\sqrt{\frac{\sigma ^2}{w^2}+1}1\right)=\underset{I}{}V_I(\sqrt{w^2_I}w_I)=\rho _{loss}(x)𝑑x^n=R^{}R,$$
(33)
which is a sum of contributions from all cells and is a monotonous ascending function of $`\sigma /w`$. In the process of the cell division $`\omega _I\omega _{Ia}\omega _{Ib}`$, we decrease $`\sigma /w`$ step by step, by playing with the geometry of the cell split, so that the gain in the total $`R_{loss}=R^{}R`$ due to a given cell split is as big as possible. It is convenient that the contributions from the other cells to $`R_{loss}`$ are unchanged. In this way every cell split will lead to a smaller and smaller $`\sigma /w`$ in the final MC run.

## Appendix B Output of the demonstration program in C++

```
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
F                                                                              F
F                   ****************************************                   F
F                   ******     TFOAM::Initialize      ******                   F
F                   ****************************************                   F
F                                                      FoamX                   F
F    Version =            2.05   =               Release date:  2002.02.13     F
F       nDim =          0 =  Dimension of the simplical sub-space              F
F       kDim =          2 =  Dimension of the hyper-cubical sub-space          F
F     nCells =       1000 =  Requested number of Cells (half of them active)   F
F     nSampl =        500 =  No of MC events in exploration of a cell cell     F
F       nBin =          8 =  No of bins in histograms, MC exploration of cell  F
F   EvPerBin =         25 =  Maximum No effective_events/bin, MC exploration   F
F   OptDrive =          2 =  Type of Driver   =1,2 for Sigma,WtMax             F
F    OptEdge =          0 =  Decides whether vertices are included in the MC   F
F    OptPeek =          0 =  Type of the cell Peek =0,1 for maximum, random    F
F     OptOrd =          0 =  Root cell hyp-cub. or simplex, =0,1               F
F   OptMCell =          1 =  MegaCell option, slim memory for hyp-cubes        F
F   OptDebug =          1 =  Additional debug histogram, SetDirectory(1)       F
F   OptCu1st =          1 =  Numbering of dimensions  starts with h-cubic      F
F     OptRej =          0 =  MC rejection on/off for OptRej=0,1                F
F   MaxWtRej =               2   =       Maximum wt in rejection for wt=1 evts F
F                                                                              F
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

```
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
F                                                                              F
F                   ***  TFOAM::Initialize FINISHED!!!   ***                   F
F     nCalls =     246799 =            Total number of function calls          F
F     XPrime =       1.2748942   =     Primary total integral                  F
F     XDiver =      0.27441356   =     Driver  total integral                  F
F   MCresult =       1.0004806   =     Estimate of the true MC Integral        F
F                                                                              F
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```

```
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
F                                                                              F
F                   ****************************************                   F
F                   ******      TFOAM::Finalize       ******                   F
F                   ****************************************                   F
F     NevGen =       2000 = Number of generated events in the MC generation    F
F     LastVe =         -1 = Number of vertices (only for simplical option)     F
F     nCalls =     248799 = Total number of function calls                     F
F                   ----------------------------------------                   F
F      AveWt =      0.78836949   =     Average MC weight                       F
F      WtMin =   1.3022085e-06   =     Minimum MC weight (absolute)            F
F      WtMax =       1.0227719   =     Maximum MC weight (absolute)            F
F                   ----------------------------------------                   F
F     XPrime =       1.2748942   =     Primary total integral, R_prime         F
F     XDiver =      0.27441356   =     Driver  total integral, R_loss          F
F                   ----------------------------------------                   F
F      IntMC =       1.0050877 +-    0.0052357675  = Result of the MC Integral F
F   MCerelat =    0.0052092644   =     Relative error of the MC intgral        F
F  <w>/WtMax =      0.80037512   =              MC efficiency, acceptance rate F
F  Sigma/<w> =      0.23296539   =              MC efficiency, variance/ave_wt F
F      WtMax =           0.985   =              WtMax(esp= 0.0005)             F
F      Sigma =       0.1836628   =              variance of MC weight          F
F                                                                              F
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

```