Q-PYTHIA: a medium-modified implementation of final state radiation

We present a Monte Carlo implementation, within PYTHIA, of medium-induced gluon radiation in the final state branching process. Medium effects are introduced through an additive term in the splitting functions computed in the multiple-soft scattering approximation. The observable effects of this modification are studied for different quantities as fragmentation functions and the hump-backed plateau, and transverse momentum and angular distributions. The anticipated increase of intra-jet multiplicities, energy loss of the leading particle and jet broadening are observed as well as modifications of naive expectations based solely on analytical calculations. This shows the adequacy of a Monte Carlo simulator for jet analyses. Effects of hadronization are found to wash out medium effects in the soft region, while the main features remain. To show the performance of the implementation and the feasibility of our approach in realistic experimental situations we provide some examples: fragmentation functions, nuclear suppression factors, jet shapes and jet multiplicities. The package containing the modified routines is available for public use. This code, which is not an official PYTHIA release, is called Q-PYTHIA. We also include a short manual to perform the simulations of jet quenching.


Introduction
Two of the most striking experimental observations at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory [1] are: The suppression of particles with large transverse momentum produced in nucleus-nucleus collisions compared with the expectations from an incoherent superposition of nucleon-nucleon ones, and; The disappearance or suppression of the yield of high transverse momentum particles in the region in azimuth opposite to a high transverse momentum one. These two phenomena are manifestations of medium effects which are usually known under the generic name of jet quenching.
The most common explanation for jet quenching is given in terms of radiative energy loss of a high-energy parton traversing a dense partonic medium, see the reviews [2]. Within this framework, the amount of energy loss measured allows to map medium properties as energy density, temperature, etc. Models implementing these ideas have been successful in describing the mentioned phenomena for light hadrons, while problems persist for heavy flavors, see e.g. [3].
In spite of these phenomenological successes, both the measurements and their theoretical description in terms of radiative energy loss are subject to several limitations which restrict their applicability to more general situations. On the experimental side, these measurements at RHIC suffer from a trigger bias -due to the requirement of a high transverse momentum particle in the event -which affects the production mechanism. This bias is very difficult to consider in analytical models as it comes strongly related to energy-momentum constraints, while most models have been developed within high-energy approximations.
On the purely theoretical side, the existing developments have been mainly designed for dealing with one-particle inclusive distributions, with multi-gluon emission treated by simple ansatzs [4]. Although this prescription is reasonable for the inclusive quantities characterizing RHIC measurements until recently, it clearly lacks of the needed degree of refinement for the more general situations that we are facing. Additionally, although radiative energy loss still provides the most successful phenomenology of jet quenching, other effects, as collisional energy loss (e.g. [5]) could also play a role -in fact, collisional energy loss would correspond to non-eikonal (finite energy) corrections to the computed gluon radiation spectra. Also, several implementations of radiative energy loss computed, in particular, using different prescriptions to treat the medium averages exist [6].
To better characterize the produced medium and to distinguish among these different possibilities, both more differential observables at large transverse momentum and unbiased measurements like jets, are required [7]. While RHIC is starting to look at some of these new opportunities [8], the Large Hadron Collider [9] will be the ideal place due both to the higher collision energy and to the characteristics of the detectors.
Radiative energy loss implies a modification of the standard QCD radiation pattern, whose generic expectations are a decrease of the energy of the leading particle in the jet (jet quenching), an enhancement of the intra-jet multiplicity and an increase of the typical emission angle with respect to the jet axis (jet broadening) [2,7]. While analytical attempts exist to study such modification and expectations within several approximation to the QCD cascading process (e.g. [10,11]), the proper tool for considering the QCD branching process in the final state, with full energy-momentum conservation, is a Monte Carlo simulator.
In spite of the fact that a probabilistic interpretation of radiation in a medium requires phenomenological assumptions, the practical advantages of a Monte Carlo are numerous. First, it allows the access to observables other than the limited single inclusive measurements, such as different jet shapes, jet multiplicities, multiparticle intra-jet correlations, etc. Moreover, such an implementation makes it possible to explore new physical mechanisms in jet development, as the interplay of the multi-gluon radiation with the medium length, effects of the color flow and reconnections, effects of recoil with the medium and others. A Monte Carlo tool is also of utmost importance for the experimental analysis. Codes as PYTHIA [12] or HERWIG [13] are in the core of the software that the different experimental collaborations apply for correcting the data, e.g. in the calibration of the jet energy. Clearly, under these conditions, a code including medium effects in the jet development provides an essential contribution to the corresponding analyses in the heavy-ion program. Several implementations of radiative energy loss in Monte Carlo codes exist [14]. Each of them focus on specific aspects and relies on different simplifications.
In this paper we present a Monte Carlo with medium-modified final-state radiation, based on the ideas described in [15] (a preliminary presentation of the Monte Carlo can be found in [16]). There, medium effects enter as an additive correction to the standard, vacuum splitting functions. The modification is taken from the full medium-induced gluon radiation spectrum computed in the multiple-soft scattering approximation (also known as Baier-Dokshitzer-Mueller-Peigné-Schiff, BDMPS) and used in previous phenomenology of jet quenching [2]. In this case, the spectra depends on two quantities, which, for a static scenario are the length, L, of the medium and the transport coefficient,q. The second quantity encodes all the information about the medium properties and is taken as a fitting parameter in practical applications. It also relates, through a single parameter, the inelastic energy loss and the angular broadening. This medium-modification is implemented at the level of the splitting functions in the standard final state showering routine PYSHOW in PYTHIA.
A medium-modified parton shower should also include the interplay between the finite size of the medium and the space-time evolution of the jet. This interplay is known to play an important role to suppress radiation at the inclusive one-gluon level due to destructive interference from different scattering centers. The corresponding case for multi-gluon radiation, in particular the role of the different ordering variables, has not yet been solved from first principles. Here, we will adopt a hybrid approach in which we consider virtuality as a good ordering variable, but correct for the finite formation time of the gluons at each step during evolution. In its present form, our implementation does not consider effects as the recoil of the scattering centers, and consequently elastic energy loss is not taken into account, or the modification of the color structure of the cascade by exchanges with the medium. These, and other mechanisms, can be implemented by further modification of the shower routine PYSHOW which we plan to do in the future.
The paper is organized as follows: In the next Section, we focus on the description of the implementation of medium-induced gluon radiation in the t-ordered final-state-radiation routine in PYTHIA. For a discussion of the basis of the modification of the splitting functions, we refer the reader to [15]. Some results are then presented in Section 3. We end with some conclusions and outlook in Section 4. Besides, some practical instructions to run the code are given in Appendix A. A publicly available version of our routine, called Q-PYTHIA, can be found in [17].

Medium-modified splittings
In this Subsection we briefly recall the treatment of the medium modification of parton splitting. As discussed in detail in [15], the spectrum for medium-induced gluon radiation [2] results in the sum of a vacuum term plus a medium-induced term: dI tot dz dp 2 The former can be related with the small-x limit of the standard Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) splitting functions P vac (z) [18], dI vac dz dp 2 In this two Eqs., x = 1 − z is the momentum fraction of the parent parton taken by the emitted gluon, p T the transverse momentum of the emitted gluon with respect to the direction of its parent, and C R the quadratic Casimir of the color representation of the parent parton.
The modification of the DGLAP splitting function is now performed by matching the medium and the vacuum cases, Eqs. (2.1) and (2.2), reading P tot (z) = P vac (z) → P tot (z) = P vac (z) + ∆P (z, t,q, L, E), ∆P (z, t,q, L, The correction ∆P (z, t,q, L, E) depends not only on the energy fraction z but also on the virtuality t = z(1 − z)p 2 T of the radiating parton and its energy E, and on the medium characteristics relevant for radiative energy loss: transport coefficientq (which corresponds to the mean squared transverse momentum transferred from the medium to the parton traversing it per mean free path) and medium length L. For practical reasons, we will use as variables characterizing the medium the accumulated transverse momentumqL and the gluon characteristic frequency ω c =qL 2 /2. The medium-induced gluon radiation spectrum dI med /dz dp 2 T will be taken in the multiple-soft scattering (or BDMPS) approximation [2].
This ansatz constrains only the small-x part of the correction ∆P (z, t,q, L, E). Different large-x extensions have been essayed, yielding similar results as discussed in [15]. We will use the default extension discussed there, consisting in multiplying the medium-induced gluon spectrum by (1 + z 2 )/2 if the parent parton is a quark and by z if the parent parton is a gluon (in this case we also impose a symmetrization around z = 1/2).
The additivity of a vacuum and a medium contribution is a common feature of the spectra when computed for the one-gluon inclusive case. Medium-modified splitting functions in a higher-twist formalism have been found to present the same factorization under the assumption of exponentiation of the elementary probabilities [19]. Nevertheless, a formal proof or disproof of this separation and of the factorization of no-splitting and splitting probabilities (respectively given by the Sudakov form factor, see below, and by the DGLAP splitting functions), like those available in the vacuum -see e.g. [20] and references therein -, is still missing for a radiating parton traveling through a colored medium.

Basics steps for a Monte Carlo implementation
Basically, a branching algorithm must solve the following problem: given a parton coming from a branching (or production) point with coordinates (t 1 ,x 1 ), with t 1 the virtuality and x 1 its energy fraction, which are the coordinates (t 2 ,x 2 ) for the next branching?
Ignoring, for simplicity, parton labels, the Sudakov form factor gives the probability for a parton not to branch (in a resolvable manner) while evolving from an initial scale t 0 to another scale t 1 . Consequently, ∆(t 2 )/∆(t 1 ) stands for the probability of evolving from t 1 to t 2 without branching. Thus t 2 can be generated by solving the equation

5)
R being a random number between 0 and 1. The energy fraction kept by the parton in the next branching z 2 can be diced down by solving the equation: with R ′ another random number between 0 and 1. Eqs. (2.5) and (2.6) are the two basic steps of a branching Monte Carlo algorithm. Now let us address in more detail the procedure followed by PYTHIA. The scale t of the evolution in PYTHIA is chosen to be the invariant mass squared of the decaying particles. The showering is performed by the creation of pairs of partons, conserving four-momentum at each step. If, in the course of the evolution, a parton a reaches a scale t a that is smaller than any minimal mass of any allowed branching mode, then its mass and its scale will be set to its current mass m 0 a and it will leave the parton shower without further branching. The evolution has a scale cut-off t 0 = Q 2 0 . In this setup, it is now conceptually simple to introduce our medium-modified splitting functions (2.3) which enter at steps (2.5) and (2.6) described above. A straightforward implementation is, however, not possible because of technical reasons: PYTHIA approximates the splitting functions by their z → 1 form to simplify the dicing procedure; this approximation is later corrected by a rejection method. The first change we make is then to solve all these integrals numerically at every step. So, this procedure modifies standard PYTHIA even for the vacuum case. In Fig. 1 we check that our modification of PYTHIA, in which we implement the exact vacuum splitting function without the z → 1 approximation, matches very well the results of usual PYTHIA. Once this modification is performed, we add the medium term to the vacuum splitting functions at every step of the evolution. Note also that our medium modifications are for g → gg and q(q) → q(q)g splittings [15]. We do not consider any modification of the g → qq branching because its splitting probability is not singular at z → 1.
In Eqs. (2.4) and (2.6) the value of the maximum initial virtuality, t max ∝ E 2 jet , of the lower virtuality limit, t 0 = 1 GeV 2 , the lower and upper limits of the z-integrals, z ± ≡ z ± (t), and the scale at which α s runs, as well as all other aspects of the t-ordered evolution (for example, the implementation of angular ordering which is a default), are the PYTHIAv6.4.18 defaults, see [12]. They can be changed by the user.
Finally, let us note that we have implemented the medium modifications only in the t-ordered routine for final state radiation. Actually, three variables have be chosen in the literature as evolution variables: virtuality t, squared transverse momentum p 2 T and angle θ. The fortran implementation PYTHIAv6.4 sets t as the default evolution variable, while it offers the possibility of using p 2 T which has become the default in the C++ version PYTHIAv8 [12]. On the other hand, angle θ is used in HERWIG++ [21] and t (plus angular ordering) in SHERPA [22]. Our medium modification of the splitting functions, Eq. (2.3), is completely general and applicable to any evolution variable. We plan its implementation in different variables in future works [23].

Length and energy evolution
The presence of a medium makes it necessary to consider the interplay between the momentum evolution variable (virtuality in our case) and the space-time dimensions of the medium. For single gluon emission the presence of a finite-size medium is known [2] to produce an interference effect known as the Landau-Pomeranchuk-Migdal effect. A full treatment of such problem for multigluon emission is still unsolved theoretically and, in fact, it is linked to the factorization problem mentioned previously. Here, we will use an assumption which relies on the factorization underlying our approach: we will treat the ordering as in usual PYTHIA, but introducing a correction procedure at every splitting to take into account the space-time evolution of the shower. We consider, as well, the energy degradation of the shower at each splitting -notice that the medium-modified splitting functions (2.3) depend now on the energy of the parent parton. More specifically, the length traveled by a parton before a gluon decoheres from its wave function and is radiated, can be estimated [2] by the gluon formation length l coh = 2ω i /k 2 T,i , where ω i and k T,i are the energy and transverse momentum (with respect to its parent parton) of the i-th emitted gluon, respectively.
The shower begins with a parton that faces the full length of the medium L, so the medium effects on the probability of the first branching are evaluated at L. The coherence length of the emitted gluon is then computed being its next branching evaluated at L − l coh . The process is iterated. Possible effects of dynamical expansion of the medium are handled analogously by considering that gluons are formed at times given by their formation lengths.
Also the energy degradation is considered. For a process a( , the medium effects in the branching process of a is considered at energy E b + E c , while the subsequent branchings of b and c, if any, are considered at E b and E c respectively. In our default results for the medium, both the evolution in length and the energy degradation are considered. The separate effect of these aspects of evolution will be discussed below.
All these modifications have been implemented in the t-ordered final state radiation routine PYSHOW, available as a code named Q-PYTHIA in [17]. A brief manual is provided in Appendix A. Note that neither the modification of the medium length nor the energy degradation are considered in usual analytical approaches to radiative energy loss [2,4,15]: all splittings are considered for a parton with the same energy and medium length or density.

Results
As mentioned in the Introduction, the generic expectations of medium-induced gluon radiation formalisms are: • A softening of the spectra: jet quenching.
• An increase of the multiplicity.
• An angular broadening of the jet: jet broadening.
These are the main features of parton energy loss that should be reflected in our implementation. However, the energy-momentum conservation within the parton shower may play a distorting role on these naive expectations based on the one-gluon inclusive spectra, as we will see later. Notice, also, that an energy-momentum flow between the parton shower and the medium (not taken into account in our present implementation) could further affect the cascade, especially for the softer particles.
To see the results, we run PYTHIA with our medium modifications. For this model description, we first focus on representative results which illustrate the effects of medium-induced gluon radiation in the branching process, to provide in the last Subsection some examples of practical interest. The statistics we have used is 10 5 generated events for each case. Unless otherwise indicated, the maximum initial virtuality has been taken as t max = 2E 2 , with E the energy of the parton which initiates the shower.

Results at parton level
To see the results at parton level, we run our medium-modified PYSHOW on a gluon of energy E jet = 10 and 100 GeV moving along the positive z-axis, and we study the intra-jet distribution of final partons in energy fraction (actually in ξ = ln (E jet /p) with p = | p| the modulus of the momentum of the final particle -the hump-backed plateau plot), transverse momentum p T = p 2 x + p 2 y of the final particle, and polar angle θ = acos(p z /p) of the final particle. The medium length is fixed to L = 2 and 5 fm, considering an homogeneous piece of matter defined by the planes (x, y, 0) and (x, y, L).
The results are shown in Fig. 2. We observe a suppression of low-ξ particles and a large enhancement of particles with large ξ-values, as expected due to the energy degradation of the leading particle (energy loss). We also observe a suppression of high-p T particles and the corresponding  enhancement of intermediate-p T particles. The p T -spectrum should be softer than vacuum at low transverse momentum since low-p T particles should be kicked towards higher values of the transverse momentum. However we find a clear enhancement of low-p T particles. Here the energy-momentum conservation in PYTHIA is causing a large effect as any increase of the multiplicity means that the total jet energy is shared among a larger number of partons. The angular distribution broadens also with increasing transport coefficient, as expected.

Effects of hadronization
To introduce hadronization effects, two back-to-back gluons of the same energy E jet running along the z-axis with opposite directions, in color singlet, are considered. PYTHIA hadronization routines are run after final state radiation. All final hadrons as produced by standard PYTHIA are taken into account. Now we consider an homogeneous piece of matter defined by the planes (x, y, −L) and (x, y, L), with the back-to-back pair initially produced at (0, 0, 0).
The corresponding results are shown in Fig. 3. The multiplicity enhancement at large ξ and small p T is less pronounced, particularly at small energies. We conclude that, as expected, hadronization tends to wash out medium effects in the soft region. The broadening of the angular distribution remains, and the region perpendicular to the initial parton directions is partially filled.

Different effects on evolution: medium length and energy degradation
In Fig. 4 we study the different effects discussed in Subsection 2.3. More concretely, we show three cases: (1) a shower with neither energy degradation nor evolution in length, so medium effects are computed at the initial values in all branching ('no evolution in energy or length'); (2) considering the energy degradation but not the evolution in length ('evolution in energy but not in length'); (3) our default: energy degradation and evolution in length ('evolution in energy and length'). We see that the medium effects are the largest in case (2). This is due to the fact that radiative energy loss is roughly independent of the parton energy [2] so, as the energy of the parton is reduced at each branching, they become more and more important. We also see that the evolution in length, case (3), causes a decrease of the medium effects with respect to the previous case, as expected.

Four examples: fragmentation functions, R AA , jet shapes and jet multiplicities
In Subsections 3.1 and 3.2 we have shown several results which illustrate the impact of our modification of the final-state parton shower on the hump-backed plateau, transverse momentum and angular distribution of partons and hadrons. In this Subsection we will indicate how some other observables are modified. We do not aim here to any quantitative comparison with present or future experimental situations, a work which is left for future studies.
First, we show in Fig. 5 the medium modification of the fragmentation function for π 0 's and hadrons expressed through the ratio of the fragmentation function in medium over that in vacuum. These fragmentation functions have been obtained through hadronization using the same setup as in Subsection 3.2: two back-to-back gluons in color singlet. Actually the information is exactly the same as that contained in the hump-backed plateau, though plotted in different variables. The effects already observed in the hump-backed plateau: enhancement of small-z and suppression of high-z hadrons with increasingq or L, are evident here. Second, we show results for the nuclear suppression factor, defined as the ratio between the number of particles measured over the expectation from a superposition of independent nucleon-  Figure 4: Upper plots: Intra-jet parton distributions in ξ (left), p T (middle) and θ (right) for a gluon of initial energy E jet = 10 GeV in a medium of length L = 5 fm and for different transport coefficientsq = 0 (black) and 10 (colored lines) GeV 2 /fm. In red we show the results for no evolution in energy or length (1), in green those for evolution in energy but not in length (2), and in blue those for evolution in energy and length (3), see the text. Lower plots: Id. but for E jet = 100 GeV and L = 2 fm.
nucleon collisions. Its deviation from 1 is a reflection of the nuclear effects. The measurements at RHIC [1] of its value being much smaller than 1 in central AuAu collisions, was the indication of jet quenching. In Fig. 6 we present results for the nuclear suppression factor for charged particles defined as For the latter the geometry is that of a 0 − 10% central PbPb collisions and the treatment of the production points and of the quenching parameters is done like in the PQM model [24]. In short, the production points of the hard scatterings are distributed in the nuclear overlapping area according to the probability of binary nucleon-nucleon collisions. Then,q and in-medium path length L are computed locally through two integrals of the density of binary nucleon-nucleon collisions along parton trajectories isotropically distributed in azimuth, see the Appendix A. In this model the only free parameter is the scale of the transport coefficient k (in fm). Using the value reported in Ref. [24] to reproduce single-inclusive RHIC R AA , k = 6 · 10 6 fm (which corresponds to an average q = 14 GeV 2 /fm in this case) we get a suppression factor of order 5 at p T > 5 GeV in semi-quantitative agreement with RHIC experimental data [1]. While no comparison to experimental data is aimed here, we note that these results are in good agreement with those from Ref. [24] in the PQM model which considers the energy loss of the leading parton through the simple ansatz [4] usually assumed in previous jet quenching phenomenology at RHIC. This implies that the introduction of evolution in virtuality and of energy-momentum conservation does not modify strongly, as also observed in [15], the conclusions about medium characteristics at RHIC extracted in previously existing frameworks.
Third, we present some results on jet shapes. In Fig. 7 we compute the jet shape in pp collisions, with and without medium corrections, at √ s N N = 5.5 TeV using the anti-k T recombination algorithm [25] within the FastJet package [26]. For the produced particles considered for jet definition, we follow the Les Houches accord [27] (considering all decays within 10 mm from the production point). We proceed by finding the hardest jet in the event with radius R j = R max = 1.5 and transverse momentum p j T . Its constituents are recursively reclustered with smaller resolution R sj < R j into subjets of transverse momentum p sj T . The jet shape φ(r) = p T sj / p T j is the fraction of the jet transverse momentum inside the hardest subjet. For the generation of jets, at minimum p min T = 100 GeV is required in PYTHIA. Quenching is implemented for the geometry 0 − 10% central PbPb collisions as it was done for the previous example of the computation of R AA .
We find that in a vacuum jet, ≃ 80% of the initial jet transverse energy is kept inside a subjet of radius R sj = 0.3 (corresponding to r = R sj /R j = 0.2 in the figure). Quenching leads to a redistribution and broadening of the jet transverse energy: according to our plot, ≃ 60% of the initial jet energy is resolved with a subcone of radius R sj = 0.3 in the case of q = 14 GeV 2 /fm. Finally, we show in Fig. 8 the distributions in jet multiplicity (i.e. the number of reconstructed jets in an event) for the partonic decay products of a gluon jet of energy 100 GeV, using the anti-k T recombination algorithm, versus the radius of the algorithm. As expected, for fixed radius the distributions with a higher number of jets are enhanced by medium effects. Implementing cuts on the minimum p T of the reconstructed jets shifts all the distributions towards smaller radii. subtraction Note that all the extracted conclusions about jets can only be taken as illustrative of the effects of quenching in our model. Further ongoing work on jet reconstruction within a heavy ion environment using different jet-definition algorithms and state-of-the-art background subtraction techniques [28] is required to further substantiate these conclusions.

Conclusions and outlook
We present a medium modification of the final-state branching process implemented through an additive term in the splitting functions following the ideas in Ref. [15]. In order to develop a tool for public use we have modified the parton shower routines within PYTHIA [12] in two ways: first we have included the exact leading order vacuum splitting functions through numerical integrations -this modification was checked to agree with normal PYTHIA; second we added a term to these splitting functions given by the medium-induced gluon radiation spectra within the BDMPS approximation. This implementation, which is not an official PYTHIA release, is nicknamed Q-PYTHIA and is publicly available at the site [17].
At large jet energies, radiative processes are expected to dominate the modification of the parton shower evolution in a medium. This leading contribution is the one included here. Other contributions which could reveal to be relevant in different regions of phase space and are not yet taken into account in our formalism, will be considered in the future. These include: i) non-eikonal corrections to the splitting functions; ii) medium-modified color reconnections; iii) energy flow from and to the medium; iv) different hadronization scenarios; v) role of the ordering variable, etc.
The main aim of the present work is to provide the tools for the upcoming data on reconstructed jets in heavy-ion collisions. These data are starting to appear at RHIC energies and they are expected to be extremely relevant at LHC energies. For this reason we do not seek here any comparison with present experimental data. We have presented, instead, a set of observables to check  Figure 8: Distributions in jet multiplicity for the partonic decay products of a gluon jet of energy 100 GeV, using the anti-k T recombination algorithm, versus the radius of the algorithm. Solid lines refer to results without medium, dashed lines for a medium with L = 5 fm andq = 10 GeV 2 /fm. Black/red/blue/magenta lines refer to 1-,2-,4-and n ≥ 6-jet distributions respectively. The upper plot contains all reconstructed jets without any cut on their p T , while the lower plot shows the results for jets with p j T > 2 GeV. See the text for explanations. the feasibility of our approach, in particular to reproduce qualitative or semi-quantitative features expected for in-medium jets. We have shown the suppression of leading particles in fragmentation functions as well as the corresponding enhancement of the hump-backed plateau; the broadening of the angular distributions; and the enhancement of intra-jet multiplicities which corresponds, as expected, to a larger number of splittings. Hadronization tends to diminish the medium effects on soft particles but the expected features remain. Energy-momentum conservation is seen to play a relevant role with respect to naive expectations for some of these quantities. This fact provides further support to the use of Monte Carlo techniques for studying the present problem.
Additionally, we have presented some examples to illustrate the ability of our formalism in simulating experimentally accessible quantities. In particular, we have presented the broadening of the jet energy and the enhancement of the n-jet distributions. Both of them are computed in a very idealized setup, in particular without any background, and should not be considered, in any respect, predictions to be compared with experimental data before a rigorous feasibility analysis is performed. In order to check the degree of convergence of our approach to previous ones, we have also computed the R AA using the same geometry and parameters as in Ref. [24] obtaining a reasonable agreement with the findings at RHIC.
Compilation and running have been tested in different platforms and compilers (g77, gfortran) and with different operative systems (MacOS and several versions of linux). Note that in version 1.0.1 linking with CERNLIB is no longer required. The performance is considerably slower than standard PYTHIA, noticeably for large energies of the showering partons and large mediumqL and ω c . This is due in part to the complexity of the generation of virtualities and momentum fractions from the medium-induced splitting functions, but mainly to the larger amount of branchings that occur. Work to improve the performance, and to include more complex examples of medium models, is under progress.