RECAPPHRI2014005
Three photon production to NLO+PS accuracy at the LHC
M. K. Mandal ^{1}^{1}1mandal@hri.res.in , Prakash Mathews ^{2}^{2}2prakash.mathews@saha.ac.in , V. Ravindran^{3}^{3}3ravindra@imsc.res.in , Satyajit Seth ^{4}^{4}4satyajit.seth@saha.ac.in
[.5cm]
Regional Centre for Acceleratorbased Particle Physics
HarishChandra Research Institute, Chhatnag Road, Jhunsi,
Allahabad 211 019, India
[.5cm]
Saha Institute of Nuclear Physics,
1/AF Bidhan Nagar,
Kolkata 700 064, India
[.5cm]
The Institute of Mathematical Sciences
C.I.T Campus, 4th Cross St, Tharamani Chennai,
Tamil Nadu 600 113, India
[.5cm]
Abstract
In this paper, we present the nexttoleading order predictions for three photon production in the Standard Model, matched to the parton shower using the MC@NLO formalism. We have studied the role of parton shower on various observables and shown a selection of results for the TeV LHC.
1 Introduction
Wealth of data from the Large Hadron Collider (LHC) and the Tevatron involving large number of leptons, gauge bosons and hadrons in the final state not only provides ample opportunity to test the predictions of the Standard Model (SM), but also constrains various physics scenarios in the beyond standard model (BSM). Signatures of BSM are often plagued by the large SM background and hence careful study of wide variety of SM processes has been underway [1]. Precise predictions for such SM processes are important as the quantum corrections are often comparable to the BSM effects. In addition, they are essential to reduce the theoretical uncertainties of the leading order (LO) predictions, that arise from the missing higher order quantum corrections through the renormalisation and factorisation scales. This necessitates the calculation of the nexttoleading order (NLO) quantum effects through Quantum Chromodynamics (QCD) radiative corrections to these SM observables at the hadron colliders. Presently, the phenomenological results for almost all physical processes of interest at the LHC are available at this accuracy due to the tremendous advancement in automatising the calculation of the virtual and real emission contributions. However, the situation becomes more involved in case of nexttonextto leading order (NNLO) calculation due to a lot of technical difficulties. In order to improve the theoretical predictions in a consistent way, it is customary to take into account all the higher order contributions, which are important in the complementary kinematic regions of phase space corresponding to the phase space relevant for the fixed order evaluation. However, to act correctly on these effects invoke a lot of technical problems. In practice, it can be approximated via the parton shower (PS) algorithm, which not only gives a reasonable estimate of these effects in the collinear kinematic regions of the phase space, but also provides a very realistic final state configuration. In other words, parton level predictions have to be gone through such showering of multi partons and recombination of these partons into hadrons through a hadronisation mechanism in order to compare them against the experimental data. Such predictions require careful matching of results at various orders to avoid double counting. Thus, NLO SM results supplemented with parton showering can provide a more reliable as well as realistic predictions that can serve in testing various BSM scenarios. Till now, there exists mainly two different algorithms incorporating the matching of a NLO calculation to parton showers, namely MC@NLO [2] and POWHEG [3]. We shall adopt the MC@NLO algorithm here, which has already been implemented and completely automated in aMC@NLO [4].
In this article, we revisit the three photon production process at the LHC at NLO in QCD and present study the consequences of matching it with the parton shower. Triplephoton production provides a background to technipion production in association with a photon, where the technipion decays into a photon pair [5]. This process has already been studied at LO [6], as well as at NLO level [7] in QCD. We extend the analysis including the effect of parton shower to get a realistic estimate of various kinematical distributions. We quantify the improvement in the predictions at small transverse momentum regions of the final state particles and the stabilisation of cross section against the variation of the factorisation and renormalisation scales.
This paper is organized as follows: in section 2, we have described the details of the calculation, mainly the virtual as well as the real emission contribution. The numerical results of the fixed order calculation together with the NLO+PS accurate results have been discussed in section 3 and finally, we conclude in section 4.
2 Calculational Details
LO ()) contributions to the production of three photons at the LHC come from quark antiquark annihilation processes. At NLO in QCD, we encounter virtual as well as real emission contributions resulting from an additional parton, namely quark or antiquark or gluon. Virtual amplitudes are already at , hence only the interference of them with the LO Born amplitudes will contribute to the NLO level. The real emission processes at NLO level come from two types of processes namely gluon emissions from the LO processes and scattering of a quark (antiquark) and a gluon producing three photons along with a quark (antiquark). The ultraviolet (UV) divergences coming from the virtual contributions and the infrared (IR) divergences originated from the virtual as well as real emission contributions, need to be removed through the addition of proper counter terms. The resulting IRsafe parton level cross section up to NLO can be written as,
(1)  
The first term is the Born contribution; is the phase space measure of the three photon final states and is the observable function which depends on the kinematic variables through the momenta of the external particles i.e., . The second term corresponds to virtual corrections to the Born process. They are often divergent when the loop momentum becomes very large and these UV divergences are first regularised and then renormalised using the counter terms given in the third term. The fourth term represents the real emission contributions at the NLO level come from parton emissions from the initial and/or final state partons. Due to massless quarks, antiquarks and gluons participating in the hard processes, both virtual and real emission contributions encounter soft and collinear divergences. The divergences coming from soft gluons and from collinear partons in the final state of the real emission processes get cancelled with those coming from the virtual processes. The remaining collinear divergences from the initial states are removed by adding mass counter terms given in the last term of eq. 1. The details of obtaining UV renormalised virtual contributions are discussed in the next section. The real emission contributions and the corresponding mass counter terms are obtained with the help of MadFKS [8], a set of routines available in the aMC@NLO [4], which along with our inhouse FORTRAN routines for calculating virtual contributions, can provide results on an eventbyevent basis in terms of four momenta of all the particles involved in the scattering process and we use them to obtain the observables that we require to study. In the following subsections, we sketch a systematic outline of the complete computational procedure.
2.1 Virtual contribution
Virtual contribution comes from the interference between the Born diagrams and the one loop corrected virtual diagrams. The number of virtual diagrams to order for the three photon production is forty eight. Up to permutations of the final state photons, we find 1 pentagon diagram, 2 box diagrams, 3 triangle diagrams and 2 bubble diagrams. We have used QGRAF [9] to generate both LO and NLO amplitudes. It generates the symbolic description of the Feynman diagrams in terms of propagators and vertices. We have written a FORM [10] code, which translates the output of QGRAF into a suitable format, that can be used for further symbolic manipulations. We have supplied Feynman rules, identities for Dirac gamma matrices, equations of motion through this code and performed various simplifications at the amplitude level. The loop integrals are regulated using dimensional regularisation. Both Lorentz contractions and Dirac gamma matrix simplifications are done in spacetime dimensions. Both UV and IR divergences appear as poles in and they have been calculated using scheme. Writing the virtual contribution in the following way:
(2) 
where is the born amplitude and s’ are the distinct topologies of virtual diagrams, we compute only one particular topology and then the permutations of photon momenta and their polarisations gave us the remaining contributions.
The reduction of tensor integrals to scalar ones in dimensions is done using the standard procedure a` la PassarinoVeltman [11]. The tensor integrals that appear at one loop level are of the form
(3) 
where
(4) 
One can decompose the above tensor integral in terms of scalar coefficients as follows:
(5) 
where the square bracket implies the nonequivalent symmetrisation by giving the full set of nonequivalent permutations. We have written a FORM code for the purpose of doing this tensor reduction and expressed the virtual contributions in terms of these scalar coefficients. As described in [12], these coefficients are related to the scalar integrals in different spacetime dimensions in the following way,
(6) 
where is a generalized scalar integral in shifted spacetime dimension. These integrals in the shifted dimensions can be expressed in terms of integrals in dimensions using the dimensional recurrence relations discussed in [13]. In this approach, inverse Gram determinants that result from the recurrence relations, often spoil the numerical stability of the integral. There exists a handful of solutions to this problem in the literature [14, 15, 16, 17, 18, 19, 20, 21, 22]. Recently, an elegant approach has been put forward in [23], where the authors have found signed minor algebraic relation, which avoids the appearance of inverse Gram determinants and thereby introducing a set of higher dimensional scalar integrals to cope with the small Gram determinants. These higher dimensional scalar integrals have been evaluated numerically after employing a series expansion in the small Gram region. This whole algorithm has been implemented in the numerical package, named PJFry [24, 25], which we use to evaluate numerically the scalar coefficients of the tensor integral for every phase space point in dimensions. PJFry reduction library uses QCDLoop [26] and OneLOop [27] to evaluate the scalar integrals in dimensions. In order to validate our FORM codes, namely those ones that perform conversion of output of QGRAF to FORM readable symbolic expressions, reduction of tensor integrals to scalar coefficients and also to validate FORTRAN routines, which evaluate the virtual contributions numerically using PJFry, we recalculated the virtual corrections of the diphoton production process in both SM and BSM to order . We compared our results thoroughly against the results presented in [28, 29] and found an excellent agreement between these two. Using our FORM codes and FORTRAN routines along with the publicly available packages, viz. QGRAF, PJFry, QCDLoop and OneLOop, we have evaluated the virtual contributions to the three photon production process at level. We find that after UV renormalisation, the IR poles namely double and single poles in are in accordance with the expectation. We express the virtual contribution of the three photon production in a form suitable for further analysis as follows:
(7) 
where is the strong coupling evaluated at the the renormalisation scale , is the partonic center of mass energy and the colour factor is: for . comes from the colourlinked Born amplitude , whereas denotes the finite virtual contribution that has been computed numerically. Note that the IR poles in are in agreement with the universal behaviour of soft and collinear partons.
2.2 Real emission contribution
Real emission contributions come from gluon emission from the Born processes as well as from the scattering of a quark/antiquark and a gluon producing a quark/antiquark and three photons. We use aMC@NLO [4] framework not only to compute these contributions along with the mass factorisation terms required to remove the initial state collinear singularities, but also to obtain the NLO results matched with PS. Within aMC@NLO, the standalone package MadGraph [30] generates all the required matrix elements both at LO as well as at NLO level. As already discussed in the previous subsection 2.1, we have prepared a set of external codes to deal with the virtual correction part and made an interface to implement it within MadFKS [8], which separates out the soft and collinear configurations in the real emission processes using the FKS subtraction scheme [31] and provides IRdivergent and IRsafe contributions separately along with the mass factorisation terms, that take part in removing the initial state collinear singularities coming from the virtual and the real emission processes. In the FKS subtraction scheme, the phase space is partitioned in such a way that each partition contains at most one soft and one collinear divergences. This is done by introducing a set of positivedefinite functions, where the s’ are chosen in such a way that they vanish in all singular limits not related to: (i) a particle becoming soft, (ii) particles and becoming collinear, obeying the restriction that the sum over all such pairs must be equal to identity. This ensures that each term of the sum is finite throughout the phase space, except when the energy of particle goes to zero or particles and become collinear. Now, after finding out the exact position of the divergences for a given partition, the generalized plus distribution is used to regulate them. All these steps are systematically automated in MadFKS within the MadGraph5 environment. We have explicitly checked the cancellation of the soft and collinear divergences among the virtual, real and mass factorisation terms at different regions of the phase space thereby confirming the perfect implementation of all the above mentioned external inputs within the aMC@NLO framework. The events, that are generated using aMC@NLO, also include the Monte Carlo counter terms to take care of the MC@NLO matching and thereby preventing the occurrence of any double counting at the time of matching to PS. These events are then showered by HERWIG [32], PYTHIA [33] parton shower to get the realistic events.
Photons are produced not only at the partonic level, but also through the fragmentation of partons into photons and a jet of hadrons can often be collinear to them. This necessitates the inclusion of nonperturbative fragmentation functions. At NLO level, the QED collinear divergence can arise when one of the final state parton becomes collinear to a photon. This can be factorised in a universal manner and then removed by adding counter terms, which renormalise the fragmentation functions, thereby bringing in a scale dependence at the partonic cross sections through the fragmentation functions, which is known as fragmentation scale. An alternate isolation criteria has been proposed in [34], using which one can obtain an observable in which fragmentation contribution is minimised and at the same time, the IR safety of that observable is guaranteed. We call it Frixione isolation here after and use this isolation for our analysis. It works in the following way: define a cone centered around each photon with a radius in the rapidityazimuthal angle () plane, where . Now, it is demanded that the sum of hadronic transverse energy inside any concentric circle of radius would be less than an amount given by the function . This function can be chosen in such a way that lesser and lesser hadronic energy is allowed as we move closer to a given photon. Because of the fact that goes to zero as , the partons that are collinear to photon are removed while the soft partons are kept intact thereby guaranteeing the QCD IR safety. For our analysis, we have taken the following canonical choice for , i.e.,
(8) 
where is the transverse energy of the photon and , , are three parameters that are to be set while applying this isolation criteria.^{5}^{5}5Effects of photon frangmentation and different isolation prescriptions have very recently been studied [35]
3 Numerical Results
In this section, we present the results for various kinematic distributions relevant to the production of three photon in SM at the LHC with the centerofmass energy TeV. Here we list the input parameters used for the whole computation:
(9) 
These values of , and ensure that the mass of the Wboson ( GeV) and the value of () remain closer to the experimental values. We have considered massless quarks with five flavours () throughout our calculation.
In our present study, we have used MSTW2008(N)LO parton distribution function with errors estimated at 68 CL for the (N)LO and it also sets the value of the strong coupling at (N)LO in QCD. The factorisation scale () and the renormalisation scale () are set equal to a central scale, which is the invariant mass of the three photon final states i.e., .
For the fixed order (N)LO calculation, we have taken the following choices of cuts: rapidity of each photon , separation between any two photons in the () plane , where . In addition, we have studied a variety of differential distributions applying two types of cuts on the transverse momentum of each photon i.e., GeV and GeV in the fixed order analysis. Unless stated otherwise, we consider GeV as the generic choice of cut on photons transverse momenta. Parameters involved in the Frixione isolation are set as: , and .
LHC  LO [pb]  NLO [pb]  Kfactor 

GeV  2.36  
GeV  2.16 
3.1 Fixed order Analysis
In table 1, we have shown the results of total cross sections for fixed order LO and NLO using the central choice of and for two different cuts. To begin with, we present some distributions of few selective kinematical variables at fixed order LO and NLO. Photons are ordered according to their transverse momentum. The hardest photon with maximum transverse momentum is denoted by . Like wise, represents the second hardest photon and the softest photon is labelled as . In fig. 1, we have shown transverse momentum distribution of at LO and NLO in the left panel and in the right panel, distribution of invariant mass of the three photon system has been plotted. The lower insets show the binbybin distribution of the Kfactor for the corresponding observable. We find that, for low transverse momentum, the Kfactor is large as it is due to the fact that the recoil against the extra parton helps to fulfil the transverse momentum cut, which was not possible at LO. The left panel of the fig. 2 shows the distribution of the rapidity of the hardest photon, whereas the rapidity of the three photon system is shown in the right panel. The distribution of the Kfactor is shown for the corresponding variables in the lower insets. Unlike fig. 1, the Kfactors in fig. 2 appear to be mostly steady indicating the affinity of these observables towards the photons having fairly high transverse momenta.
n  [pb]  



0.4 
1  
2  
1  1  
2 
All the above distributions show a substantial effect of radiative corrections on this process. This is mainly
because of the inclusion of new subprocesses at the NLO, as quarkgluon subprocesses begin to contribute at this order and
due to the enhancement in the phase space. In fig. 3, we have plotted the separation between the ordered
photons in the () plane obeying the selection cut: , where .
We have checked the rapidity differences between these photons are quite small. Therefore, the peaks arising in these distributions
near the angle (180), suggest that the emitted photons are mostly backtoback. The hardest photon is separated form
the softest one i.e., , by at least at LO, whereas at NLO they can be very
close as permitted by the selection cut due to the emission of an extra radiation at this level.
Besides, we have checked the effect of variation of Frixione isolation parameters i.e., , and n. Though Frixione isolation has no effect on the LO cross section, the dependency of the NLO cross section on these isolation parameters is shown in table 2. From eq. (8), it is evident that the NLO crosssection increases when decreases and it also increases with increasing . In fig. 4, we have shown the transverse momentum distribution of the hardest photon by varying the value of (left panel) from to for a fixed value of and =1, where the Kfactors vary from to . The right panel shows the same distribution for a fixed value of and and in this case, Kfactors vary from to . It is evident from table 2, as well as from fig. 4 that, for a fixed choice of value, the NLO crosssection is large for and , whereas it becomes much smaller for and indicating the fact that smaller increases the crosssection when is larger [36]. It is also clear from table 2, that the effect of varying , keeping n and fixed, is quite minimal. Similar studies with changing the value , provide same kind of distributions analogous to fig. 4.
3.2 Discussion on NLO+PS
In this section, we compare the fixed order NLO result with the NLO results matched with PS (NLO+PS) with two different showering algorithm, namely HW6 and PY6. For the showering purpose, parton level events are generated using very loose cuts: GeV, , with the following Frixione isolation parameters: , and . We have explicitly checked that the events thus produced, remain unbiased in total rates and differential distributions after showering and hadronisation for this choice of kinematical cuts and Frixione isolation parameters. These events are then showered with HERWIG6 (HW6) and PYTHIA6 (PY6) and we have imposed the same set of analysis cuts that we used in the fixed order analysis along with the generic cut on the transverse momentum of the photon at the time of showering.
The scale dependencies of the results are calculated by varying and independently around the central value via the following assignment: and , where and are varied between the range [1/2,2] independently. Various ratios of , and that appear as arguments of logarithms in the perturbative expansion to NLO are within the range [1/2,2]. The scale uncertainty band is the envelope of the results obtained by varying this and within this range [29]. The PDF uncertainties are estimated with the Hessian method, as given by the MSTW [37] collaboration. We have plotted fractional uncertainty, which is defined as the ratio of the variation about the central value divided by the central value, being a good indicator of the uncertainties. These uncertainty bands can be generated automatically at the time of parton level event generation by storing additional information, sufficient to determine via a reweighting technique, at no extra CPU cost within the aMC@NLO framework as described in [38].
We have shown distribution for HW6 and PY6 together the fixed order NLO result, in fig. 5. It is clear that at low values, NLO+PS (for both HW6 and PY6) result shows the effect of all order resummation of the large logarithms, hereby suppressing the cross section leading to a meaningful value, while the fixed order NLO result diverges for . At low , PY6 result is different from the HW6 result as the soft and collinear emissions constituting the parton shower are treated differently. PYTHIA generates more softer spectra than HERWIG in this region and as a result of this, these two showers show different behaviour as expected [39]. At high , the NLO fixed order and NLO+PS (for both HW6 and PY6) results are in agreement as in this region, the hard emissions are dominant and they are correctly described by the NLO hard cross section. In the middle and lower insets of fig. 5, we have presented the fractional scale and PDF uncertainties of the NLO+PS result for HW6 and PY6 respectively which increase with increasing [39]. We do not find any significant differences in case of studying fractional uncertainties using these two different showers. Therefore, in the rest of the figures, we present the fractional uncertainty plots only for HW6.
We now present the results for various kinematical distributions to NLO accuracy, matched with PS (labelled as NLO+PS), for both HW6 and PY6 with the specified analysis cuts. We have adopted a consistent pattern for all the rest of the distributions. In each case, within the main frame, three curves corresponding to the distributions in fixed order NLO (solid red) and NLO+PS using HW6 (dashed blue) and NLO+PS using PY6 (dotted black) are shown. The middle inset shows fractional scale uncertainty (dashed cyan) and fractional pdf uncertainty (solid violet), while the lower inset shows the ratio between NLO+PS and NLO for HW6 (dashed blue) and for PY6 (solid black).
In the left panel of fig. 6, we have shown the plots for transverse momentum distribution of the hardest photon and the right panel shows the distribution of the invariant mass distribution of the three photons. We do not find much difference in the results of two showers HW6 and PY6. In both distributions, NLO results are very little larger than the NLO+PS results and this is due to the fact that the QCD radiation becomes softer when we demand all the photons to satisfy a high cut (i.e. GeV) and damping of PDFs at large Bjorken values further subdue its effect at the parton level, whereas aMC@NLO produces more events with hard and central jets resulting in the suppression of these distributions after showering. In fig 7, we have depicted the plot showing rapidity distribution of the three photon system and as expected, we observe that the NLO result is slightly harder than the NLO+PS result. However, the ratio in the lower inset shows that PY6 generated events give larger contribution than HW6 in the large rapidity region indicating that PY6 produces significantly large number of radiations than HW6 in the full kinematically available phase space.
4 Conclusion
Precise and realistic predictions of both signal and background processes at hadron colliders are now possible due to tremendous developments in the computational methods and the availability of the state of the art computational tools. We have used packages, namely QGRAF, PJFry, aMC@NLO to study the three photon production process at the NLO level in QCD for the LHC taking into account the parton shower effects and realistic experimental cuts. In addition, we have developed some codes that build the interfaces among these different analytical and numerical tools. We have plotted different kinematic observables and discussed the consequences of showering the fixed order NLO results with two different showering algorithm HERWIG and PYTHIA. We have also discussed the effects of scanning over the Frixione isolation parameters on the NLO cross section. We find our predictions are less sensitive to scale uncertainties and choice of PDFs and hence more suited for direct comparison with the data from the experiments.
Acknowledgments
We would like to thank Rikkert Frederix for many helpful discussions and acknowledge Valery Yundin for his help in running PJFry package. The work of MKM and VR has been partially supported by funding from Regional Center for Acceleratorbased Particle Physics (RECAPP), Department of Atomic Energy, Govt. of India. PM and SS would like to acknowledge HPC cluster computing facility of Theory Division, SINP.
References
 [1] J.M. Campbell, J.W. Huston, W.J. Stirling, Rep. Prog. Phys. 70 (2007) 89, The SM and NLO Multileg Working Group: Summary report, arXiv:1003.1241, The SM and NLO Multileg and SM MC Working Groups,: Summary Report, arXiv:1203.6803.
 [2] S. Frixione and B. R. Webber, JHEP 0206, 029 (2002) [hepph/0204244].

[3]
P. Nason,
JHEP 0411, 040 (2004)
[hepph/0409146];
S. Frixione, P. Nason and C. Oleari, JHEP 0711, 070 (2007) [arXiv:0709.2092 [hepph]].  [4] R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, R. Pittau and P.Torrielli, Phys. Lett. B701 (2011) 427; R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, R.Pittau and P. Torrielli, JHEP 09 (2011) 061.
 [5] A. Zerwekh, C. Dib, R. Rosenfeld, Phys. Lett. B549, 154158 (2002). [hepph/0207270].
 [6] M. Golden, S.R. Sharpe, Nucl. Phys. B261 (1985) 217, V. Barger, T. Han, Phys. Lett. B212 (1988) 117.
 [7] G. Bozzi, F. Campanario, M. Rauch, D. Zeppenfeld, Phys. Rev. D84 (2011) 074028.
 [8] R. Frederix, S. Frixione, F. Maltoni and T. Stelzer, JHEP 10 (2009) 003.
 [9] P. Nogueira, Journal of Computational Physics 105 (1993) 279.
 [10] M.Tentyukov and J.A.M. Vermaseren, hepph/0702279.
 [11] L. M. Brown and R. P. Feynman, Phys. Rev. 85 (1952) 231â244, G. ât Hooft and M. J. G. Veltman, Nucl. Phys. B153 (1979) 365â401, G. Passarino and M. J. G. Veltman,Nucl. Phys. B160 (1979) 151.
 [12] A. I. Davydychev, Phys. Lett. B263 (1991) 107.
 [13] O.V. Tarasov, Phys. Rev. D54 (1996) 6479, J. Fleischer, F.Jegerlehner, and O.V. Tarasov, Nucl. Phys. B566 (2000) 423.
 [14] W. L. van Neerven and J. A. M. Vermaseren, Phys. Lett. B137 (1984) 241
 [15] G. J. van Oldenborgh and J. A. M. Vermaseren, Z. Phys. C46 (1990) 425.
 [16] J. M. Campbell, E. W. N. Glover, and D. J. Miller, Nucl. Phys. B498 (1997) 397.
 [17] T. Binoth, J. P. Guillet, and G. Heinrich, Nucl. Phys. B572 (2000) 361.
 [18] F. Campanario, JHEP 1110, 070 (2011), F. Campanario, Q. Li, M. Rauch, M. Spira, JHEP 1306, 069 (2013), F. Campanario, H. Czy Ì , J. Gluza, M. Gunia, T. Riemann, G. Rodrigo, V. Yundin, JHEP 1402, 114 (2014)
 [19] W. T. Giele and E. W. N. Glover, JHEP 04 (2004)029.
 [20] A. Denner and S. Dittmaier, Nucl. Phys. B734 (2006) 62â115.
 [21] T. Binoth, J. P. Guillet, G. Heinrich, E. Pilon, and C. Schubert, JHEP 10 (2005) 015.
 [22] T. Diakonidis, J. Fleischer, J. Gluza, K. Kajda, T. Riemann, and J. B. Tausk, Phys. Rev. D80 (2009) 036003.
 [23] J. Fleischer and T. Riemann, Phys. Rev. D83 (2011) 073004.
 [24] V. Yundin, Ph.D thesis, HumboldtUniversitat zu Berlin, 2012, http://edoc.huberlin.de/docviews/abstract.php?id=39163.
 [25] https://github.com/Vayu/PJFry/.
 [26] R.K. Ellis and G. Zanderighi, JHEP 02 (2008) 002.
 [27] A. van Hameren, Comput. Phys. Commun. 182 (2011) 2427.
 [28] M. C. Kumar, P. Mathews, V. Ravindran and A. Tripathi, Nucl. Phys. B 818 (2009) 28. M. C. Kumar, P. Mathews, V. Ravindran and A. Tripathi, Phys. Lett. B 672 (2009) 45
 [29] R. Frederix, M. K. Mandal, P. Mathews, V. Ravindran, S. Seth, P. Torrielli, M. Zaro, JHEP 1212 (2012) 102;
 [30] T. Stelzer and W.F. Long, Comput. Phys. Commun. 81 (1994) 357, J. Alwall et al., JHEP 0709 (2007) 028.
 [31] S. Frixione, Z. Kunszt and A. Signer, Nucl. Phys. B467 (1996) 399, S. Frixione, Nucl. Phys. B507 (1997) 295
 [32] G. Marchenini et al., HERWIG, Comput. Phys. Commun. 67 (1992) 465, G. Corcella et al., HERWIG 6.5, JHEP 01 (2001) 010, G. Corcella et al., HERWIG 6.5 release note, hepph/0210213.
 [33] T. Sjostrand, S. Mrenna and P.Z. Skands, JHEP 05 (2006) 026.
 [34] S. Frixione, Phys.Lett. B429 (1998) 369.
 [35] J. M. Campbell, C. Williams, arXiv:1403.2641.
 [36] Vittorio Del Duca, Fabio Maltoni, Zoltan Nagy, Zoltan Trocsanyi, JHEP 04 (2003) 059
 [37] A. D. Martin, W. J. Stirling, R. S. Thorne, and G. Watt, Eur. Phys. J. C63 (2009) 189â285.
 [38] R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, R. Pittau and P. Torrielli, JHEP 02 (2012) 099.
 [39] Paolo Torrielli, Stefano Frixione, JHEP 04 (2010)110.