Particle Event Generator: A Simple-in-Use System PEGASUS version 1.0

PEGASUS is a parton-level Monte-Carlo event generator designed to calculate cross sections for a wide range of hard QCD processes at high energy $pp$ and $p\bar p$ collisions, which incorporates the dynamics of transverse momentum dependent (TMD) parton distributions in a proton. Being supplemented with off-shell production amplitudes for a number of partonic subprocesses and provided with necessary TMD gluon density functions, it produces weighted or unweighted event records which can be saved as a plain data file or a file in a commonly used Les Houches Event format. A distinctive feature of PEGASUS is an intuitive and extremely user friendly interface, allowing one to easily implement various kinematical cuts into the calculations. Results can be also presented"on the fly"with built-in tool \textsc{pegasus plotter}. A short theoretical basis is presented and detailed program description is given.


Program summary
Solution: Experimental investigations of hadron scattering processes at high energy colliders, such as the Large Hadron Collider (LHC), usually imply multiparticle final states and implement complex kinematical restrictions on the fiducial phase space. The multipurpose Monte-Carlo event generators are commonly used tools for theoretical description of the collider measurements. Most of them (for example, pythia 8.2 [6], mcfm 9.0 [7], madgraph5 amc@nlo [8] and other) use conventional (collinear) QCD factorization, which is based on the DGLAP resummation. On the other side, the hadron-level Monte-Carlo event generator cascade [9] and parton-level generator katie [10] can deal with non-zero transverse momenta of incoming off-shell partons. In particular, cascade employs the CCFM equation for initial state gluon evolution. pegasus is a newly developed parton-level Monte-Carlo event generator designed to calculate cross sections for a wide range of hard QCD processes, which also incorporates the TMD gluon dynamics in a proton. It provides all necessary components, including off-shell (dependent on the transverse momenta) production amplitudes and grid files for TMD gluon densities, interpolated automatically. pegasus offers an intuitive and extremely user friendly interface, which allows one to easily implement various kinematical cuts into the calculations. Generated events (weighted or unweighted) can be stored in the Les Houches Event file or presented "on the fly" with convenient built-in tool pegasus plotter.
The present version of pegasus is applicable for Fermilab Tevatron and CERN LHC processes.

Theoretical framework
Here we present a short review of ideas and frameworks which form the theoretical basis of pegasus. For more information we address the reader to reviews [1,2].

Kinematics and cross section formula
pegasus operates in the general framework of conventional Parton Model extended to k T -factorization approach [15,16]. It employs the factorization hypothesis to calculate the cross section of a physical process through the convolution of the (TMD or collinear) parton densities and hard scattering amplitudes. By default, the k T -factorization scheme is assumed and can be switched to the collinear QCD factorization for each of colliding particles.
pegasus refers to 2 → 1, 2 → 2 and 2 → 3 partonic subprocesses. The initial partonic state is described with Sudakov variables, namely, the light-cone momentum fractions x 1 and x 2 and non-zero parton transverse momenta k 1T and k 2T . The n-particle final state phase space is parameterized in terms of particle rapidities y i , transverse momenta p iT , and azimutal angles φ i , where i = 1 ... n [17]. The fully differential cross section reads where the 4-momenta for all particles are given in parentheses,ŝ = (k 1 + k 2 ) 2 is the subprocess invariant energy, F is the flux factor (see below) and µ 2 is the factorization scale. The longitudinal momentum fractions x 1 and x 2 are not the integration variables; they are obtained from the energy-momentum conservation laws in the light-cone projections: where m iT is the transverse mass of produced particle i. In the optional choice of collinear QCD factorization, we replace the TMD parton densities with conventional distributions, omit the integration over k 2 1T and k 2 2T and take the on-shell limit in the scattering amplitude A ab as described below.

Off-shell partonic amplitudes
The calculation of partonic amplitudes follows the standard Feynman rules, with the exception that the initial gluons are off-shell.
Off-shell gluons may have nonzero transverse momentum and an admixture of longitudinal component in the polarization vector. In accordance with the k T -factorization prescriptions, the initial gluon spin density matrix is taken in the form [16]: In the collinear limit, when k T → 0, this expression converges to the ordinary µ g * ν g = −g µν /2. This property provides continuous on-shell limit for the partonic amplitudes.
In general, off-shell initial states may cause violation of gauge invariance in the nonabelian theories. We solve this problem in a way [18]. We start with a set of "extended" diagrams, where the off-shell gluons are considered as emitted by on-shell external fields (say, quarks), so that they represent internal lines in the Feynman graphs. As all of the external lines in these graphs are on shell, the gauge invariance of the whole set is fulfilled. However, the extended gauge invariant set may contain unfactorizable diagrams that cannot be represented as a convolution of the hard scattering amplitude and gluon density functions (such as, for example, Figs. 5(e) and 5(d) in [18]). At the same time, we learn from [18] that the non-factorizable diagrams vanish in the particular gauge (3). Therefore, within this gauge, we are left with the usual diagrams of the Parton Model.
In fact, the prescription (3) is a remake of Equivalent Photon Approximation (EPA) in QED and DIS. Consider a photon emitted by an electron: e(p) → e (p )+γ(k). Then, taking trace over the electron line in the matrix element squared one obtains the polarization tensor Neglecting the second term in the right hand side in the small x limit, p k, one immediately arrives at the spin structure µ * ν ∼ L µν ∼ p µ p ν . It can be rewritten in the form (3) after using the Sudakov representation k = xp+k T and applying a gauge shift µ → µ − k µ /x. The gauge invariance of the matrix element is correct to the accuracy of O(x) 10 −4 − 10 −6 .

CCFM evolution equation
The CCFM evolution equation [4] for TMD gluon density resums both BFKL [3] large logarithms α n s ln n 1/x and α n s ln n 1/(1 − x) terms. In the limit of asymptotic energies, it is almost equivalent to BFKL, but also similar to usual DGLAP evolution [5] for large x. Moreover, it introduces angular ordering condition to treat correctly gluon coherence effects. In the leading logarthmic approximation (LLA), the CCFM equation for TMD gluon density f g (x, k 2 T , µ 2 ) can be written as where k T = q(1 − z) + k T andP gg (z, k 2 T , q 2 ) is the CCFM splitting function: The Sudakov and non-Sudakov form factors read: ln ∆ ns (z, whereᾱ s = 3α s /π. The first term in the CCFM equation is the initial TMD gluon density f (0) g (x, k 2 T , µ 2 0 ) multiplied by the Sudakov form factor, decribing the contribution of non-resolvable branchings between the starting scale µ 2 0 and scale µ 2 . The second term represents the details of the QCD evolution expressed by the convolution of the CCFM gluon splitting functionP gg (z, k 2 T , q 2 ) with the gluon density and Sudakov form factor. The angular ordering condition is introduced with the theta function, so the evolution scale µ 2 is defined by the maximum allowed angle for any gluon emission [4]. The main advantage of this approach is the implicit including of higher-order radiative corrections (namely, a part of NLO + NNLO +... terms corresponding to the initial state real gluon emissions). Some details can be found, for example, in reviews [2].
A similar equation can be also written [19] for TMD valence quark densities with replacement of the gluon splitting function by the quark one 1 . One usually takes the initial TMD gluon and valence quark distributions as where σ = µ 0 / √ 2 and q v (x, µ 2 ) is the conventional (collinear) quark density function. The parameters N , B and C can be fitted from collider data (see [19]).
The TMD (valence) quark densities for CCFM approach are not available in the current version of the pegasus.

Parton Branching approach
The Parton Branching (PB) method [21,22] allows one to solve the DGLAP equations iteratively. It gives a possibility to take into account simultaneously soft-gluon emission at z → 1 and transverse momentum q T recoils in the parton branchings along the QCD cascade. The latter leads to a natural determination of TMD density functions for both gluons and quarks. A soft-gluon resolution scale z M is introduced to separate resolvable and non-resolvable emissions, which are treated via the DGLAP splitting functions P ab (α s , z) and Sudakov form-factors, respectively. The PB equations for TMD parton densities read: where a = q or g and k T = q(1 − z) + k T . The Sudakov form factors are defined as The real-emission branching probabilities P (R) ab are obtained from splitting functions P ab by eliminating δ(1−z)-terms and substituting 1/(1−z) + → 1/(1−z). The evolution scale µ 2 can be connected with the emission angle with respect to the beam direction, that leads to the angular ordering condition µ = |q T |/(1 − z). The dependence on the z M dissapears when this angular ordering condition is applied and z M is large enough. The initial TMD parton distributions are taken in a factorized form as a product of conventional quark and gluon densities and intrinsic transverse momentum distributions (taken in gaussian form [8,9]), where all the parameters can be fitted from the collider data. Unlike the CCFM, the PB densities have the strong normalization property: where a = q or g and xa(x, µ 2 ) are the conventional parton density functions (PDFs). The PB equations can be solved numerically by an iterative Monte-Carlo method with the leading order (LO) or next-to-leading order (NLO) accuracy. The solution results in a steep decline of the parton densities at k 2 T > µ 2 . It contrasts the CCFM evolution, where the transverse momentum is allowed to be larger than the scale µ 2 , corresponding to an effective taking into account higher-order contributions 2 .

Kimber-Martin-Ryskin approach
The Kimber-Martin-Ryskin (KMR) approach [24] provides a technique to construct TMD gluon and quark densities from conventional ones by loosing the DGLAP strong ordering condition at the last evolution step, that results in k T dependence of the parton distributions. This procedure is believed to take into account effectively the major part of next-to-leading logarithmic (NLL) terms α s (α s ln µ 2 ) n−1 compared to the LLA, where terms proportional to α n s ln n µ 2 are taken into account. At the LO, the KMR method, defined for k 2 T ≥ µ 2 0 ∼ 1 GeV 2 , results in expressions for TMD quark and gluon distributions [24]: where P LO ab (z) are the usual DGLAP splitting functions at LO and µ 2 0 is the minimum scale for which DGLAP evolution is valid. The theta functions introduce the specific ordering conditions in the last evolution step, thus regulating the soft gluon singularities. The cut-off parameter ∆ usually has one of two forms, ∆ = µ/(µ + |k T |) or ∆ = |k T |/µ, that reflects the angular or strong ordering conditions, respectively. In the case of angular ordering, the parton densities are extended into the k 2 T > µ 2 region, whereas the strong ordering condition leads to a steep drop of the parton distributions beyond the scale µ 2 . At low k 2 T < µ 2 0 the behaviour of the TMD parton densities has to be modelled [24]. Usually it is assumed to be flat under strong normalization condition (9).
The Sudakov form factors allow one to include logarithmic virtual (loop) corrections, they take the form: with z max = 1 − z min = µ/(µ + |q T |). These form factors give the probability of evolving from a scale k 2 T to a scale µ 2 without parton emission. At the NLO, the TMD parton densities can be written as [25]: where p 2 . Note that both DGLAP splitting functions and conventional parton distributions should be taken with NLO accuracy. The Sudakov form factors at NLO read: However, it was demonstrated that the NLO prescription, with a good accuracy, can be significantly simplified to keep only the LO splitting functions [25] while the main effect is related to the Sudakov form factors.

Flux factor
The choice of the flux factor is another peculiarity of off-shell calculations. The definition of the flux, which is the velocity of an off-shell particle, is highly disputable and is not clear. By default, we accepted the analytic continuation of the general textbook definition [17]: Our choice is supported by a numerical experiment [26], where we have considered the production of χ cJ states in a two photon process, e + e → χ cJ + e +e and made a comparison between the prompt calculation of this 2 → 3 process and calculation based on EPA, γ + γ → χ cJ with J = 0, 1 or 2. We find that the flux taken in the form (21) provides a sensibly better fit to the exact calcualtion. The fact that the exact calculation disagrees with the factorized (collinear) calculation indicates that the conditions for the factorization theorem are yet not fulfilled. In such a situation, our choice of the flux can be regarded as a phenomenological correction for non-factorizable contributions. The same definition of the flux is adopted, for example, in [10]. However, the user can optionally choose the conventional definition F = 2x 1 x 2 s, as it is explained below.

Calculations using PEGASUS
pegasus has an intuitive and unprecedentedly user friendly interface. The calcualtions using pegasus include a few general steps common for all of the processes. So, when pegasus is running, one can select the colliding particles, pp or pp, and set their center-of-mass energy √ s. The default setting corresponds to the LHC Run II setup. Then one can select factorization scheme (TMD or collinear one) for each of the colliding particles, choose corresponding parton density function and set the parameters, important for further Monte-Carlo simulation, namely, number of iterations and number of simulated events per iteration (see Fig. 1). Next steps could be as follows (note that all these steps do not depend on each other and can be done in different order): • From the list of available processes one can select the necessary process and then (optionally) correct the default kinematical restrictions, hard scales, list of requested observables and corresponding binnings (see Figs. 2 -4). This can be done by double clicking on the requested process or via main menu (using Edit → Task option).
• For each of the observables one can manually edit the default binning according to own wishes. As another option, the binning can be uploaded immediately from the data file. Several commonly used formats (such as .yoda, .yaml, .csv or plain data format, compatible with gnuplot [27] tool) as provided by HepData repository [28] are supported.
• The user-defined setup for any process (total center-of-mass energy, selected parton densities, kinematical restrictions, binnings etc) can be saved to a configuration file in some internal format (*.pegasus). This can be done via the main menu (using File → Save or File → Save As options) or via the popup menu available on right mouse button click or via appropriate button in the button panel. Of course, the configuration file can be loaded and a user-defined setup can be used in further applications. This can be done via main menu (with File → Open option) or via popup menu or via Open button on the button panel.
• Weighted or unweighted events can be generated. This option is available via main menu Edit → Settings → Generated events or via popup menu.
• If one needs to generate the Les Houches Event file [29], one has to mark corresponding option before the calculation starts (see Fig. 1). Note that this option affects the speed of the calculations.
• The calculation will start by choosing the corresponding option in main menu (Calculation → Start), popup menu or pressing Start button on the button panel. The numerical results for requested observables will be immediately presented "on the fly" with built-in tool pegasus plotter (see Fig. 5). The calculations can be paused or even stopped (using main menu options Calculation → Pause, Calculation → Stop, corresponding buttons on the button panel or options in popup menu). Of course, during pause in the calculations, all manipulations with accumulated results in pegasus plotter are allowed (see Section 4.3).
• If there are several contributing subprocesses, there is a possibility to immediately jump to the next one (via Calculation → Next option in main menu or appropriate button on button panel or popup menu) during the calculations.

Parton density functions in a proton
Latest sets [19] of CCFM-evolved TMD gluon distributions in a proton, JH'2013 set 1 and JH'2013 set 2, are available in the pegasus. The input parameters of JH'2013 set 1 were determined from the fit to high precision HERA data on the proton structure function F 2 (x, Q 2 ), whereas the parameters of JH'2013 set 2 gluon were extracted from combined fit on both F 2 (x, Q 2 ) and F c 2 (x, Q 2 ) data (see [19]). The previous CCFM fits, namely, A0 and B0 sets [30], are available also and comparison between them can be found [19]. Technically, all these CCFM-evolved gluon densities are stored on a grid in log x, log k T and log µ and a simple linear interpolation is applied to evaluated the gluon  density for values in between the grid points. This interpolation proceeds automatically at the each event generation. The data files containing the grid points are supplied with the program package and read in at the beginning of each requested calculation 3 . The TMD gluon and quark densities can be also evaluated in the standard DGLAP scenario using the PB and/or KMR schemes (with LO or NLO accuracy). As an input for the KMR procedure, the conventional (collinear) PDFs have to be applied. Several sets of KMR and PB-based gluon densities are currently available in pegasus. So, well known NNPDF3.1 (LO) [32] and MMHT'2014 (LO) [13] parametrizations and recent analytical expressions obtained [33][34][35] in the so-called generalized double asymptotic scaling (DAS) approximation of QCD [34,35] were used as an input. The DAS approximation is connected to the asymptotic behaviour of the DGLAP evolution discovered many years ago [36]. Flat initial conditions for the DGLAP equations, applied in the generalized DAS scheme, lead to the Bessel-like behaviour for the proton PDFs at small x [34,35]. The DAS LO set 1 corresponds to "frozen" treatment of the QCD strong coupling in the infrared region: α s (µ 2 ) → α s (µ 2 + m 2 ρ ). The DAS LO set 2 is based on the idea [37] regarding the analyticity of the strong coupling at low scales. The difference between these two choices is discussed [38]. Everywhere, the cut-off parameter ∆ is taken according to angular ordering condition.
The available TMD gluon sets with the essential parameters are listed in Table 1. We note that large variety of the TMD gluon distribution functions in a proton are collected in the tmdlib package [39], which is a C++ library providing a framework and an interface to the different parametrizations.
To perform the calculations in the collinear QCD factorization (mainly at LO or treelevel NLO for some processes) several sets of conventional PDFs are available in pegasus. These are recent MMHT'2014 (LO and NLO) [13] and CT14 (LO and NLO) [40] ones and previous sets provided by CTEQ Collaboration (CTEQ'6.6) [41], MSTW'2008 (LO and NLO) [42] and GRV'94 (LO) [43]. The standalone C++ codes from MMHT [13] and CTEQ [14] groups are implemented into the pegasus code and corresponding data files containing the grid points are supplied with the program package.

Simulated processes
List of processes, currently available in the pegasus, is presented in Table 2. We would like to clarify some points, which are not mentioned in the Table 2: • Q denotes a heavy (c or b) quark, Q denotes B c , B ( * ) c , J/ψ, ψ , Υ(3S), χ cJ (1P ) or χ bJ (3P ) mesons with J = 0, 1 or 2, V stands for γ or Z/γ * and l = e or µ.
• The standard spectroscopic notation (2S+1) L (a) J is used for intermediate Fock state of heavy quark pair QQ produced with spin S, orbital angular momentum L, total angular momentum J and color representation a.
• Besides the scales, listed in the Table 2, for every process some universal scales are available, namely: total energy of the partonic subprocess √ŝ , also ŝ/4 and so-called CCFM scale ŝ + Q 2 T [19,30], where Q T is the total tranverse momentum of the final partons, see Fig. 3.
• To estimate the scale uncertainties, the hard scales above can be varied around its default value by a factor of 2 or 1/2, as it often done in the pQCD calculations. This applies to any scale, except the CCFM scale. For CCFM-based gluon densities the scale uncertainties are evaluated by a change of the default gluon distribution to so-called "+" and "-" ones of the corresponding set [19,30].
• Part of the non-logarithmic loop corrections to effective g * + g * → H 0 vertex can be absorbed in the special K-factor [51] and optionally implemented into the calculations. As default choice, this K-factor is switched off.
• According to experimental setup, an isolation criterion is applied for prompt photon production. This criterion is the following: a photon is isolated if the amount of hadronic transverse energy E had T deposited inside a cone with aperture R centered around the photon direction in the pseudo-rapidity and azimuthal angle plane, is smaller than some value E max T ("cone isolation"). The isolation requirement significantly (up to ∼ 10% of the visible cross section) reduces contribution from the so-called photon fragmentation mechanisms, not implemented into the pegasus. The isolation criterion and additional conditions which preserve our calculations from divergences have been specially discussed, for example, in [52]  Open heavy flavor production Associated gauge boson and factorization has been chosen, since in the k T -factorization approach its contribution is taken into account by the off-shell gluon-gluon fusion (see [49,50] for more details).

Quarkonium final states
Quarkonium production processes need additional explanations as they contain an important extra step: the formation of bound states.
The process starts with the production of a heavy quark-antiquark pair QQ in a hard parton interaction. The produced quark pair may be either color singlet or color octet. If the QQ state is color singlet, it can immediately convert into a meson with appropriate quantum numbers. Then the formation probability is determined by a single parameter, the radial wave function at the origin of the coordinate space, |R S (0)| 2 or |R P (0)| 2 [53][54][55][56][57]. The values of these functions can be extracted from the measured decay widths or calculated within potential models.
The situation with color-octet states is more complicated as the transition to a colorsinglet physical hadron requires the emission of extra (soft) gluons. Then, for every final state hadron h and for every intermediate QQ state n = 2S+1 L J listed in Table 2 there is a specific phenomenological long-distance matrix element (LDME) [58][59][60] responsible for such a transition O h [n] . Eventually, the production cross section for a hadron h in pp collisions is given by a sum over all possible singlet and octet QQ states: where σ pp → QQ [n] is the partial partonic cross section (1) for a QQ state n. The different states n are selected by introducing the proper projection operators in the hard scattering amplitude. The correspondence between the color singlet and color octet wave functions and respective LDMEs is given by the O h 2S+1 L J = 2N c (2J +1)|R S (0)| 2 /4π and O h 2S+1 L J = 6N c (2J + 1)|R P (0)| 2 /4π, respectively. As default choice to describe the spin structure of relevant transition amplitudes, pegasus uses the model [61], where the NRQCD emission of soft gluons is considered in terms of classical multipole radiation theory. The multipole expansion is dominated by (chromo)electric dipole transitions E1. According to [61], only a single E1 transition is needed to transform a P -wave state into an S-wave state and the structure of the respective 3 P 1 + g amplitudes is taken the same as for radiative decays of χ cJ mesons [62,63]: where p, k and l = p − k are the four-momenta of the color-octet P -wave state, emitted gluon and produced color-singlet S-wave state, µ (k), µ (l), µ (p) and µν (p) are the polarization vectors (tensor) of respective particles and e µναβ is the fully antisymmetric Levi-Civita tensor. The transformation of an S-wave state into another S-wave state (such as J/ψ or ψ meson) is treated as two successive E1 transitions 3 S J + g, 3 P 2 . For each of the two transitions the same effective coupling vertices (23) - (25) are exploited. In the case of J/ψ or Υ(3S) production, the user can optionally include feeddown from the decays of upper quarkonium states (ψ , χ cJ (1P ) or χ bJ (3P ), respectively) in addition to the direct production channels.
The absolute normalization of the transition amplitudes is not calculable within the theory; these numbers are taken as free adjustable parameters. However, the values of the LDMEs for the different partial contributions are not completely independent but are connected by the heavy quark spin symmetry (HQSS) relations [60]. All the HQSS relations between LDMEs are implemented in the pegasus default setting, taken from [46,47].

Strong coupling and masses of particles
In pegasus, the strong coupling α s can be calculated in one-loop or two-loop approximation with respect to the number of active flavors N f and Λ QCD . The choice of N f and Λ QCD is done automatically (according to the Table 1) with the choice of the TMD and/or conventional parton density in a proton. There is no possibility to change it manually since this setup is essential for determination of corresponding parton distributions.
The masses of all particles (quarks, gauge bosons, heavy quarkonia etc), their branching ratios and decay widths are fixed according to Particle Data Group (PDG) [64]. Any of these parameters can be easily changed using the convenient built-in particle data tool (see Fig. 6), which is available via main menu (Edit → Settings → Particle data) or popup menu or appropriate button on the button panel.

Generation of Les Houches Event file
pegasus is supplied with a tool to construct event files in the Les Houches Event format [29]. This is a widely accepted format to present events, which is compatible with the majority of modern general purpose Monte-Carlo generators.
The Les Houche Event (LHE) file, generated by the pegasus, consists of two main blocks: the first one contains the information about the number of the recorded events, the PDFs used, the colliding hadrons and their energies. Also the total cross section is shown. The second block represents a list of events, including the data on the interacting partons, their 4-momenta and color structure of the event. We mention the basic features of the LHE file, generated by the pegasus: • The generated events could be weighted or unweighted. In first case, the sum of all the weights is the total cross section of the subprocess.
• Polarization information is not preserved.
• A parton carries a tag according to the standard PDG numbering scheme [64].
• Conventional (collinear) parton densities in a proton are numbered according to the lhapdf scheme [65] while TMD parton distributions are numbered according to the tmdlib package [39].
The produced LHE file can then be processed with an external program to introduce some peculiar event selection, to include parton showers, to hadronize the final particles, etc. It is found to be compatible with such Monte Carlo generators as pythia [6] and cascade [9].

Program components 4.1 Random number generator
Since all the internal variables in pegasus are declared as double precision ones, double precision random numbers have to be generated in the Monte-Carlo simulations. The random number generator ranlux [66] is well suited for these purposes. It has a long period, solid theoretical foundations and is commonly used in computational physics. This random number generator is implemented into the pegasus.

Phase space integration and event generation
The multidimensional phase space integration (1) is performed with the Monte-Carlo technique and is incorporated with the routine vegas [12]. The routine vegas allows up to ten integration variables, that is enough for subprocesses considered in pegasus.
The vegas algorithm is based on a method for reducing statistical errors by using a known probability distribution function to concentrate the search in those areas of the integrand that make the greatest contribution to the final integral.
The algorithm is realised through a large number of random sample points distributed over a d-dimensional rectangular volume. The whole volume is divided into d-dimensional rectangular cells (by default, 50 divisions along each axis). The probability for a point to drop into a given cell is determined by so called sampling distribution, which is adjusted to the integrand function. The sampling distribution approximates the exact distribution by making a number of passes (iterations) over the integration region while histogramming the integrand function dσ given by the expession (1). Each iteration is used to define a sampling distribution for the next iteration. To improve the convergence in the region of high p T , the user can optionally modify the integrand function from dσ to (1 + p 4 T )dσ, see Fig. 1. The optimization of the sample grid is made automatically and needs no care from the user.
Each sample point generated by vegas represents an event in the n-particle phase space with the coordinates of the sample point responding to the values of the physical integration variables (p iT , y i , φ i ). The weight attributed to that event is given by the product of the integrand function dσ and the d-volume of the sampling cell containing that point, dw = dσdV cell . Unweighting algorithm for the generated events is provided also.
The typical time needed to generate one event, of course, depends on the requested subprocess; but, in general, is similar to time needed by other Monte-Carlo event generators like cascade or pythia. The output events can be plotted on a histogram (using built-in tool pegasus plotter ) or stored for further use in the form of an Les Houches Event file.

PEGASUS Plotter
pegasus is supplied with a built-in tool pegasus plotter, allowing one to depict easily the produced differential cross sections and immediately compare them with experimental data. As a default setting, the accumulated results for requested observables during the calculation are shown in pegasus plotter. However, it is a quite independent tool and can be used apart from any calculations made within the pegasus. The program is very simple and intuitive. Let us briefly describe main features of the tool.
As one calls the pegasus plotter from the main menu of the pegasus (using the Tool → Plotter option, or via popup menu, or by pressing the corresponding button on button panel) an empty sheet is created. The following objects, stored in a plain data files, could be added to the sheet (by choosing Edit → Add option in the main menu or popup menu available with a right mouse button click on the sheet): • Curve. The data file should consist of two columns, corresponding to the rows of x and y values.
• Histogram. The data file is the same, as for Curve. However, every y value should be mentioned twice, for the both borders of the corresponding bin.
• Filled area. A three-columns data file should contain for every x value lower and upper values of y.

Installation and running
pegasus can be downloaded from https://theory.sinp.msu.ru/dokuwiki/doku. php/pegasus/download as a precompiled executive file for Linux machines. No special installation procedure is needed. The data files containing the necessary TMD parton densities in a proton (and conventional PDFs as well) are located inside the pegasus home directory (folder data). If there are some missing data files in data folder, pegasus will inform user about that (see Fig. ?). In this case, no calculation is allowed. Location of data folder could be easily changed via main menu option Edit → Settings → Path to data folder.
For Linux machines, the executive file can be just run from a terminal as ./PEGASUS. The program demands the qwt library (version 6.1.3) [11], so the library file libqwt.so.6 should be inside the pegasus home directory. Otherwise, the path to this file should be specified with export LD LIBRARY PATH=/path/to/library.
The program was tested on ROSA Linux R8.1, ROSA Linux R11 and Ubuntu 16.04.