{'fraction_non_alphanumeric': 0.04116117850953206, 'fraction_numerical': 0.012478336221837088, 'mean_word_length': 4.606266699052708, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Kaleu: a general-purpose parton-level phase space generator

Kaleu is an independent, true phase space generator. After providing it with some information about the field theory and the particular multi-particle scattering process under consideration, it returns importance sampled random phase space points. Providing it also with the total weight of each generated phase space point, it further adapts to the integration problem on the fly. It is written in Fortran, such that it can independently deal with several scattering processes in parallel.


Introduction
Efficient phase space generation is an important issue in the study of multi-particle processes at collider experiments. Parton-level defined observables are eventually formulated as phase space integrals, and the evaluation of these is the most time-consuming part of a calculation. The complicated structure of both the integrand and massive multi-particle phase space enforces the use of Monte-Carlo methods, setting the total evaluation time to the sum of the generation of phase space points and the evaluation of the integrand at the phase space points. While it is important to keep the computing time of each integrand evaluation as low as possible, it is also important to keep the number of evaluations to reach an acceptable accuracy as low as possible. The latter is commonly referred to as efficient phase space generation, and eventually is a trade off between the number of evaluations and the average computing time for a single generation of a phase space point.
Since efficient phase space generation should already be settled at leading order (LO) calculations, automatic systems for LO calculations developed over the last decade deal with this issue [2][3][4][5][6][7][8][9][10][11]. Various techniques are used in various combinations. Efficient phase space generation implies knowledge about the integrand, which is utilized both before the integral is performed, through phase space mappings, and during the integration process, using the actual value of the integrand at the phase space points. Methods using the latter are referred to as adaptive methods, and the ones used within particle physics are based on adaptive grids dividing phase space into subspaces [12][13][14][15][16], concern the optimization of the parameters in mixture distributions (the so-called adaptive multi-channel method) [17], or combine both [18]. Most important is, however, to use as much as possible information about the integrand before performing the integral, be it just when deciding which underlying phase space mapping to use in adaptive grid methods.
Eventually, the aim of the game is to have the probability density with which phase space points are generated look as much as possible as the integrand and to let it have the same peak structure. The fact that, for parton-level multi-particle calculations, this structure is essentially determined by the sum of the squares of the tree-level Feynman graphs contributing to the scattering amplitude is of great help. Since each Feynman graph encodes the phase space mapping needed for the efficient integration of its own square, the multi-channel method, in which each graph represents a channel, can be applied for efficient performance of the full integration. The problem with this approach is that the number of Feynman graphs grows factorially with the number of final-state particles. This problem was solved for the calculation of the tree-level amplitudes themselves through the use of recursive techniques [8,[19][20][21][22][23]. Regarding the phase space generation, the straightforward solution within the multi-channel method is to keep only a restricted number of channels with the highest channel weight. However, in order to determine these, the integration process still has to be started including all channels.
The problem only exists in the calculation of the weight coming with each phase space point. Only there the sum over all graphs occurs, not in the actual generation of the phase space points, for each of which only one graph encoding a phase space mapping is chosen. The reason why the weight cannot be calculated with the same recursive approach like the amplitude is that, within the multi-channel method, each individual graph has its own unique weight. And indeed, the problem is solved by departing from this approach, as was pointed out in [25], and instead giving individual weights to the vertices. Instead of choosing a whole graph encoding a mapping at once for each generation of a phase space point, a mapping is composed by choosing respective branchings. Each branching corresponds to a vertex, and has its unique weight among the list of possible branchings at a given stage in the construction of the phase space mapping. The weight associated with the generated phase space points can then be calculated following recursion relations analogous to those for the amplitude calculation.
As said, the foregoing works well assuming the peak structure of the integrand is determined by sum of the squares of the tree-level Feynman graphs contributing to the scattering amplitude, which however is not a gauge-invariant quantity. For processes involving many gluons, other phase space mappings have been designed [26][27][28].
Although many implementations of the various methods for efficient phase space generation exist, most of them are heavily integrated in automatic systems for full multi-particle calculations, making it difficult to use them separately. This write-up presents a program with which Kinematics Are Lucidly and Efficiently Utilized 1 . It makes a number of routines available which can be used for phase space generation in any given Monte Carlo program. As part of the initialization, the user has to provide masses and widths of all possible real and virtual particles in the multi-particle process under consideration, and a list of interaction vertices. KALEU uses this information to efficiently map phase space following the method of the Recursive Phase space Generator described in [25]. The generator is optimized on the fly during the integration process, with the option to also optimize the generation of invariants and polar angles further with PARNI [16]. It is a FORTRAN program, written such that several instances of the program can operate in parallel, completely independently of each other. It can be obtained from http://annapurna.ifj.edu.pl/˜hameren

The algorithms
As mentioned before, KALEU uses the algorithm of the Recursive Phase space Generator from [25]. Like most algorithms, it uses the fact that n-body phase space can be decomposed in a sequence of two-body phase spaces, referred to as splittings in the following.

Two-body phase space generation
In each splitting, two new four-momenta p 1 , p 2 are generated that sum up to an, possibly previously generated, existing four-momentum Q = p 1 + p 2 . Depending on whether 0, 1 or 2 of these momenta are on-shell, 2, 3 or 4 random variables have to be generated from which p 1 , p 2 are constructed. Most common is to generate the invariants p 2 1 , p 2 2 , the cosine z of the polar angle of p 1 in the center-of-mass frame (CMF) of Q, and the the azimuthal angel φ in the same frame. If any of the momenta is on-shell, and less than 4 random variables have to be generated, we say "the invariant is generated with a Dirac delta-distribution".
The azimuthal angle is usually not correlated with any existing kinematical variable, but the zvariable sometimes is. It may be the cosine of the polar angle w.r.t. to another, already generated, momentum in the CMF of Q. The z-variable is then linearly related to the scalar product of the two momenta. This happens, for example, in the so called open antenna generation in [28] in order to arrive at antenna densities for the efficient integration of multi-gluon amplitudes. It also happens in the the so-called t-type generation [1]. In this case, the z-variable is taken to be the cosine of the polar angle w.r.t. to one of the initial-state momenta q 1 in the CMF of Q, and is linearly related to t = (p 1 − q 1 ) 2 = (Q − q 1 − p 2 ) 2 . If this generation is formulated as 1: , these four momenta can be interpreted as the external legs of a t-channel Feynman graph. This is useful for bookkeeping purposes, as it will also be for KALEU. Both this t-type generation and the "s-type" generation with uncorrelated z variable are applied in KALEU.
The mentioned different choices of random variables are applied in order to have control over the densities following which they are distributed. They should match the behavior of the the squared matrix element one is trying to integrate. KALEU uses the densities as mentioned for example in [4].

Recursive phase space decomposition
We will illustrate how sequences of splittings are composed with the help of an explicit example, namely the process uū → ddZ. For this process, KALEU will create the list of splittings/mergings, or just vertices, in Table 1. The list contains almost all possible three-point vertices occurring in Feynman graphs for this process, i.e. all three-point vertices to be calculated in a recursive calculation of the amplitude. A few are missing, as will be explained below, and a few appear several times with differences in the last column, which will also be explained. The list only contains three-point vertices, the meaning of the particles between square brackets is explained below. There is no distinction between particles that only differ in charge since it does not matter for the kinematics. The number between parentheses encodes the momentum in the binary representation. The external particles have momenta u(1)u(16) → d(2)d(4)Z(8), and for example the gluon/photon/Z-boson attached to the d(2)d(4) pair has momentum 2 + 4 = 6.
To generate a phase space point, the list should be read upside down, and each vertex represents a splitting. The starting point is the sum of all final-state momenta 2 + 4 + 8, plus one initial-state momentum 1. The latter is necessary to include t-type kinematics. Thus a vertex is chosen among number 18 to 33. Each of these carries its own multi-channel weight, which is updated during the Monte Carlo process. Suppose vertex 29 is chosen, having W(3) , d(12) on the r.h.s. Then the same is repeated for all vertices with W(3) on the l.h.s. and with d(12) on the l.h.s.. For the latter, this involves only vertex 5. For the former, no such vertex exists in the list, because it would represent the trivial operation of subtracting the initial-state momentum 1 from 3 to obtain the final-state momentum 2. Thus the constructed Feynman graph representing the phase space mapping is given by (15) d (12) W(3) 29 5 (1) Vertex 29 represents a t-type splitting, as all vertices do with an odd momentum on the l.h.s., because they contain momentum 1. The first step in a t-type splitting is the generation of the positive invariants, in the example above the squares of momenta 12 and 2. The latter is actually external and fixed to the squared mass of the d-quark, so it is "generated following a deltadistribution". One of these invariants, in this case the square of momentum 2, is beyond the encoding of the three-point vertex, and has to be listed separately. This is the purpose of the particles between square brackets in the list of vertices. Splittings of the s-type, like vertex 12, do not have entries in this last column, as well as t-type vertices which give the possibility to be followed by another t-type splitting, like vertex 25.
For the calculation of the weight associated with a phase space pointed generated as above, the list of vertices has to be read in the order as it stands. Each "←→" should now be interpreted as "= +" i.e. each vertex represents a contribution to the "off-shell current" labeled by particle-type and momentum on the l.h.s.. These off-shell currents are probability densities, and the final density labeled by u(15), obtained after executing the whole list of vertices, is the reciprocal of the weight. Each vertex includes the density associated with the variables generated for the splitting, and the multi-channel weight for the choice of the vertex among the possibilities contributing to the same off-shell current. This way, the contribution of all Feynman graphs is included, underlining the recursive character of the algorithm.

Usage
KALEU has been designed such that the user can conveniently write his own user-defined routines suitable for a given Monte Carlo program. Some examples are included in the program. In particular, an interface to replace PHEGAS with KALEU in the HELAC/PHEGAS system [4,23,24] is included.
The first thing the user has to do is provide KALEU with all particles and possible vertices. For the example in Section 2 the following routine would be sufficient: The output model and vertx are of public derived type with private components, so the user can only declare them and pass them to routines from KALEU, in this case addparticle and addvertex. The input of these routines should be clear, it should only be mentioned that the order of the particles in addvertex does not matter. In practise, the user would of course extend the model to, for example, the full standard model, which was refrained from for brevity here. Some prepared model routines are included in the program. A routine as above would typically be called in a user-defined initialization routine to be called in the Monte Carlo program. A complete set of such user-defined routines could look as follows: module my_kaleu use avh_kaleu use avh_kaleu_model use avh_kaleu_kinem type(model_type) ,save :: mdl ! Masses and widths type(kinem_type) ,save :: kin ! Kinematics type(kaleu_type) ,save :: obj ! Instance of the phase space generator end module subroutine my_init( process ,nfinst ,ecm ,nbatch,nstep,thrs ) use my_kaleu integer ,intent(in) :: process(-2:17),nfinst,nbatch,nstep real(kind(1d0)) ,intent(in) :: ecm,thrs type(vertx_type) :: vtx call my_model( mdl,vtx ) call kaleu_put_process( mdl,vtx,kin,obj ,process ,nfinst ,ecm ) call kaleu_init_adapt( mdl,obj ,nbatch,nstep,thrs ) end subroutine subroutine my_gnrt( discard ,pkaleu ) use my_kaleu logical ,intent(out) :: discard real(kind(1d0)) ,intent(out) :: pkaleu(0:3,-2:17) call kaleu_gnrt( mdl,kin,obj ,discard ,pkaleu ) end subroutine subroutine my_wght( weight ) use my_kaleu real(kind(1d0)) ,intent(out) :: weight call kaleu_wght( mdl,kin,obj ,weight ) end subroutine subroutine my_collect( weight ) use my_kaleu real(kind(1d0)) ,intent(in) :: weight call kaleu_collect( obj ,weight ) end subroutine The module my_kaleu contains essentially what in a traditional FORTRAN program would be the common blocks. The items are of public derived type with private components, so the user can only declare them and pass them to routines from KALEU. By declaring them as arrays, several instances of the program can be created. The user-defined routines above are wrappers of KALEU-routines which need the items in my_module as arguments.
The initialisation routine my_init takes the process as input, and the particles should have the same integer labelling as in the user-defined my_model. The entries −2 and −1 of process are the initial-state particles, the (strictly) positive entries contain the final-state particles. The integer nfinst is the number of final-state particles. The input ecm is the centerof-mass energy. The numbers nbatch, nstep and thrs are related to the multi-channel optimization. This happens in nstep steps of nbatch non-zero weight events. After these nstep steps, each channel with a weight smaller than thrs times the average channel weight at the off-shell current it belongs to is discarded.
KALEU determines limits on two-particle invariants based on the masses of final-state particles, but these may be equal zero and in a practical calculation there will be phase space cuts present. The user may translate these into two-particle invariant mass cuts, and feed them to KALEU for more efficient phase space generation with call kinem_updt_smin( kin ,smin ) The entries smin(1:nfinst,1:nfinst) of the array smin(-2:17,-2:17) should be minima to final-final state two-particle invariants, and the entries smin(-2:-1,1:nfinst) should be negative, and should be maxima to final-initial state two-particle invariants. The array is considered to be symmetric.
All arguments of the generation routine my_gnrt are output. If the logical discard is true, the phase space point should get weight zero, and the array pkaleu must not be used. If discard is false, pkaleu contains the momenta, including those of the initial-state particles. The third routine my_wght returns the weight associated with the most recently generated phase space point, and the final routine my_collect takes the full weight, including the integrand, as input for optimization purposes.
The routines above would suffice to deal with e + e − scattering. In case the user would like to generate events for hadron-hadron scattering, initial-state momenta can be provided to KALEU for each event with call kinem_inst( kin ,discard ,pkaleu ) This routine should be called just before kaleu_gnrt. The array pkaleu has the same format as before, but is input here. The entries pkaleu(0:3,-2) and pkaleu(0:3,-1) are taken as the initial-state momenta. The logical discard is output, and is true if something is wrong with the momenta. Alternatively, the user can let KALEU deal with the x-variables within collinear factorization. For this option, one item should be added to the module, and a few routine calls should be added. An example is given in Appendix A. The optimization of the generation of the x-variables is performed with the help of PARNI.
An example of the use of more instances of the program is given in Appendix B.

Results
In order to assess the performance and applicability of KALEU, some results existing in literature are reproduced.  for some e + e − → 6f processes as in [29]. The numbers in the column "AMEGIC++" and "H/PHEGAS" are taken directly from that paper. The number N eval refers to the run with "H/KALEU". Table 2, Table 3, Table 4 and Table 5 show some results for e + e − → 6f processes as presented in [29]. It concerns all processes at √ s = 500 GeV, and "with QCD" if applicable. The numbers in the columns "AMEGIC++" and "H/PHEGAS" are taken from that paper. The numbers in the column "H/KALEU" were obtained by replacing PHEGAS with KALEU in the HELAC/PHEGAS system. In [29] results where obtained using 10 6 phase space points before cuts. It is, however, not clear whether this includes an optimization phase. The results with KALEU were obtained with 10 7 phase space points before cuts, including optimization. In order to compare the performances, both the standard deviation σ and √ 10 σ are presented explicitly for the latter. Some results with "H/PHEGAS" have been rounded further than in the original paper for clarity. This leads to an ambiguity for one value, which has been highlighted. The optimization was performed in 10 steps of 50 × 10 3 phase space points after cuts. These were not included in the estimation of the cross section. The last column in the tables presents the total number of phase space points passing the cuts and leading to a matrix element evaluation, including the 500 × 10 3 evaluations used for optimization. Table 6 shows some results for e + e − → 8f processes as presented in [30]. It concerns some processes with m H = 130 GeV and m t = 174.3 GeV, at √ s = 800 GeV and with QCD. The numbers in the column "CARLOMAT" are directly taken from that paper. The numbers in the column "HELAC/KALEU" were again obtained by replacing PHEGAS with KALEU in the HELAC/PHEGAS system. They were obtained with 10 8 phase space points before cuts, including optimization. The optimization was performed in 10 steps of 100 × 10 3 phase space points after cuts, which were not included in the estimation of the cross section. The last column in the tables presents the total number of phase space points passing the cuts and leading to a matrix element evaluation.
In order to present some results for hadron-hadron collisions, KALEU has been connected with ALPGEN [7]. This program deals with partonic subprocesses in hadron-hadron collisions in one Monte Carlo run, making an optimized choice of a subprocess to be generated for each event. Because KALEU can provide an independent generator for each subprocess, it can easily  for some e + e − → 6f processes as in [29]. The numbers in the column "AMEGIC++" and "H/PHEGAS" are taken directly from that paper. The number N eval refers to the run with "H/KALEU".  Table 4: Cross sections in [fb] for some e + e − → 6f processes as in [29]. The numbers in the column "AMEGIC++" and "H/PHEGAS" are taken directly from that paper. The number N eval refers to the run with "H/KALEU".  Table 5: Cross sections in [fb] for some e + e − → 6f processes as in [29]. The numbers in the column "AMEGIC++" and "H/PHEGAS" are taken directly from that paper. The number N eval refers to the run with "H/KALEU".  Table 6: Cross sections in [ab] for some e + e − → 8f processes as in [30]. The numbers in the column "CARLOMAT" are taken directly from that paper. The number N eval refers to the run with "HELAC/KALEU".
be merged into the structure of ALPGEN. This has been done such, that ALPGEN still makes the choice of a subprocess for each event, and that, given the subprocess, KALEU generates the variables x 1 , x 2 for the PDFs and the final-state momenta. It should be noted, however, that for the purpose of this write-up, the connection between the two programs has been established only to the level of cross section calculation and not, for example, to full event un-weighting. Table 7 presents results for processes of the type pp → tt+Njets and pp → tt bb+Njets. The numbers in the column "ALPGEN" are taken directly from [7]. The user of ALPGEN chooses the number of matrix element evaluations for a Monte Carlo run, not the number of generated phase space points before cuts. The last column presents the number of evaluations for the run with KALEU, including those used for optimization but omitted for the estimation of the cross section.

Summary
The program KALEU for parton-level phase space generation was presented. Given the masses and widths of all real and virtual particles, it generates importance sampled phase space points for the scattering process under consideration. Given the total weight of each event, including the value of the integrand, it optimizes further to this integrand on the fly. It is written in FORTRAN, such that it can independently deal with several scattering processes in parallel.  (8) 0.4522 (28) 10.6e+6 tt bb + 4jets 0.189 (4) 0.1851 (14) 6.2e+6 Table 7: Cross sections in [pb] for some pp collision processes at LHC as in [7]. The numbers in the column "ALPGEN" are taken directly from that paper. The number N eval refers to the run with "ALPGEN/KALEU".