Three-jet production in POWHEG

We present an implementation of the production of three jets at NLO plus parton-shower effects in the POWHEG BOX. Using the recently introduced MiNLO procedure for setting the renormalization and factorization scales, we are able to obtain a generator that is also well behaved when the third jet becomes unresolved. We compare key distributions computed at the NLO level, at the level of the POWHEG hard emission and after full shower by PYTHIA, Pythia8 and HERWIG6. We also compare our three-jet generator with the already available dijet POWHEG generator.


Introduction
The production of hadronic jets is an ever present phenomenon in hadronic collisions. Jets are the manifestation of the production of coloured particles with large transverse momentum, and in hadronic collisions this phenomenon is very frequent, due to the relatively large size of the strong coupling constant, and to the presence of coloured incoming partons.
Electroweak processes with associated production of QCD jets are an ever present background both to Standard Model studies and to searches of new physics. It is therefore mandatory to understand these phenomena to our best. For this reason, basic QCD jetproduction processes can constitute a framework where we can test our ability to simulate jet phenomena.
We stress that QCD jet production is more difficult to understand and simulate with respect to associated jet-production phenomena. In fact, in the latter case, our initial process already involves the production of a massive object, with a relatively well defined mass. This sets the relevant scale and momentum fractions for the parton distribution functions (pdfs), and the associated jet production probes the values of pdfs around this point. On the other hand, in the basic QCD jet-production processes, these scales and momentum fractions are instead determined by the jet system, that is not as well known. Small errors in the determination of the jet energy induce larger uncertainties due to the steep fall of the luminosity as a function of the mass of the produced system. Furthermore, at relatively low total transverse energies, we are approaching the high-energy regime, that is not usually dealt rigorously by current Monte Carlo implementations. 1 Thus, a reasonable understanding of basic QCD jet simulation can increase our confidence that we can also model associated jet-production phenomena in a reliable way.
A NLO-accurate generator for dijet production that can be interfaced to parton-shower generators (i.e. a NLO+PS generator), using the POWHEG method [2,3], was implemented in ref. [4] (the dijet generator from now on), within the POWHEG BOX framework [5].
In the present work, we implement a NLO+PS generator, built in the POWHEG BOX framework, for the production of three jets (the trijet generator from now on). Basically, we include the 2 → 3 parton scattering processes and all the QCD corrections to them, that include, besides the virtual corrections, all the 2 → 4 parton scattering processes at leading order. We neglect parton masses throughout. The trijet implementation is carried out within the POWHEG BOX V2 framework. 2 The NLO virtual matrix elements for three-jet production were computed for the first time in refs. [6,7,8]. Compact expressions for the real contributions are also available from refs. [9,10,11,12]. We have used these results as coded in the c++ program NLOJET++ [13]. The other missing ingredients, needed to set up a POWHEG BOX generator, are the colourand spin-correlated Born amplitudes. These are easily obtained using the MadGraph4 [14] interface to the POWHEG BOX developed in ref. [15]. The Born phase space, due to the complex singularity structure of the Born amplitude, has been built using a multi-channel technique.
The Born process in the trijet generator has several singular kinematic regions, associated to pairs of final-state partons becoming collinear, or one parton acquiring small transverse momentum. A further, overall, singular configuration is the one where all partons have small transverse momentum. We find that the MiNLO [16] procedure for setting the scales and assigning Sudakov form factors is particularly helpful here, since it tames the divergences in all kinematic regions but the overall one. As we will discuss in the following, we also find that, when using MiNLO, two-jet inclusive observables are fairly well described, so that we do not need to worry about the impact of configurations close to the limit where one jet becomes unresolved, and furthermore we only need to deal with the problems of the overall singularity of the Born-level cross section, since the others are regulated by MiNLO.
The paper is organized as follows. In section 2 we describe the construction of the Born phase space, and the multi-channel technique that we have used in order to probe adequately all the singular regions, and other technical details about the trijet implementation. In section 3 we discuss the checks that we have carried out in order to validate our cross-section formulae. In section 4 we compare the output of our generator at different levels. After the discussion of the common settings for the comparisons in section 4.1, we compare among each other the NLO, the Les Houches Event (LHE) level and the shower results in section 4.2. The LHE level is the stage where POWHEG has already generated the hardest radiation, but no other radiation has been added by the subsequent shower programs. The purpose of this comparison is to determine how the final result is built up. In section 4.3 we compare the NLO and the MiNLO-improved NLO results, in order to show at what level the MiNLO procedure differs from the NLO results obtained with a standard scale choice. In section 4.4 we compare among each other the NLO, LHE and shower results, when MiNLO is turned on. Since MiNLO regulates the divergences related to the third jet becoming soft or collinear, but not those related to the whole event having small total transverse energy, we discuss how to enforce some physicality requests on the small transverse-energy region in section 4.5. In section 4.6, we compare the MiNLO trijet results with the dijet ones, when considering quantities inclusive in the third jet. We give our conclusions in section 5. A short discussion on the choice and setting of the dynamical scales in the POWHEG BOX is presented in appendix A.

Technical details
In this section we discuss a few technical details of the trijet implementation: the multichannel Born phase space, the generation cuts and the production of weighted events, with and without MiNLO.

Phase-space generation with multi-channel technique
The Born cross section for the production of three partons has several sharp peaks, in correspondence to the singular regions where soft and collinear singularities are approached. A standard way to integrate a many-peak function is by using a multi-channel technique. In the following we illustrate how we implemented it on trijet production.
We label the particles as follows: with 1 and 2 we indicate the two incoming partons, using the label 0 when we refer to both, and with 3, 4, 5 the final-state ones. Momentum conservation in the center-of-mass frame at the Born level is then given by (2.1) The Born cross section has 9 singular regions, according to the 9 possible choices of emitteremitted couples. We label them with two indexes: the first index identifies the emitting particle and the second the emitted one: -6 final-state regions: {35, 53, 45, 54, 34, 43}, where, for example, 35 is the singular region associated with parton 5 being emitted by parton 3.
-3 initial-state regions: {03, 04, 05}, where we treat as one the singular region associated with either of the incoming partons.
In order to perform an efficient importance sampling in each singular region, we introduce a function of the kinematic variables that approaches 1 only in one singular region and goes to zero fast enough in all the others. To do so, we define the following quantities where, for final-state partons, with E i = p 0 i , and θ ij the angle between parton i and parton j, in the center-of-mass frame, and where θ 1j is the angle between the direction of the first incoming beam and the outgoing parton j. It is clear from their definition that when two final-state partons become collinear or when one parton becomes soft, the corresponding S ij diverges. The same can be said about S 0j when the jth parton becomes soft or collinear with respect to the incoming beam. Defining we can then write the following identity where we have introduced a self-explanatory notation in eq. (2.7). Each term on the right-hand-side of eqs. (2.6) or (2.7) approaches 1 only in one particular singular region. For example, the termsS 0j approach 1 when the jth parton is either soft or collinear to any of the two incoming beams, and go to 0 when other singular regions are approached.
Similarly, terms of the formS ij approach 1 when the jth parton is either soft or collinear to the final-state parton i, and go to 0 when other singular regions are approached.
We insert then 1 written as eq. (2.7) in the formula of the invariant phase-space element We can now choose the best parametrization of the kinematic variables (i.e. the momenta p i ) in terms of the Monte Carlo integration variables, in order to do an importance sampling for each of the terms of the sum. Each dΦ B in the sum has then a different parametrization in terms of the Monte Carlo integration variables, so that each of them can be seen to depend on the summation indexes. We will indicate this by adding the subscript ij to each phase-space element volume In the trijet generator, each phase-space volume (dΦ B ) kj is computed using a replica of the automatic machinery (see refs. [3,5]) for the generation of the real phase space, starting from the Born phase space for dijet production, i.e. starting from the 2 → 2 phase space and attaching an extra parton, with different importance sampling according to the singular region where it has been adapted to. When the subroutine for the computation of the Born phase space is invoked, one extra random number is used to choose, with equal probability, one of the 9 different parametrization of eq. (2.9). The returned Jacobian is then the product of (dΦ B ) kj and of the corresponding suppression functionS kj , multiplied by 9, in order to compensate for the 1/9 factor introduced by choosing to evaluate randomly only one single term of the sum.
The real phase space for trijet production is then built as usual by the POWHEG BOX automatic machinery on top of the Born kinematics.

Generation cuts or weighted events
As already stated, the Born cross section for three-jet production has several singular regions, associated to a pair of final-state partons becoming collinear among each other, or to a final-state parton becoming collinear to an initial-state parton, or becoming soft. Because of these singularities, an unweighted generator would end up generating all events in these singular regions. This problem is usually handled by requiring that the final-state partons satisfy some generation cuts, such as to avoid the singular configurations. In this case, of course, one should make sure that the contributions arising from the neglected regions of phase space do not end up affecting observables of interest. In the case at hand, one may require that the three partons form well separated jets. If final-state observables do require at least three jets, it is unlikely that the neglected regions would contribute to them. However, one should always check that the results are independent upon these cuts.
An alternative method is to generate weighted events. One chooses a weight function that diverges when approaching the singular regions. This is done as follows. One introduces a function of the Born phase-space kinematics, F (Φ B ), that vanishes in the singular regions, such that the following integral of the differential cross section σ is finite. One then generates the phase-space points with a probability proportional to (dσ/dΦ B )F (Φ B ), assigning to each point a weight 1/F (Φ B ). In this way one should not worry about the independence of the result upon generation cuts. Observables that do depend upon the singular regions will typically receive rare contributions with large weights, yielding large errors. This method was used in the POWHEG BOX since ref. [17], where the function F was dubbed "Born suppression factor". In the POWHEG BOX framework the factor F was applied to all cross-section contributions (i.e. not only the Born term), and, in the case of real terms and collinear remnants, it was computed as a function of the underlying-Born kinematics. Here we stress that, in spite of the name, no phase-space regions are really suppressed. In fact, the effect of the F factor at the generation level is exactly compensated by the 1/F weight of the event.
In the present case, we have considered the following F function (2.12) where S 1 and S 2 are suitably chosen scales, q 1 , q 2 and q 3 are the square of the transverse momenta of the three final-state particles with respect to the beam axis, and q ij is the relative transverse momentum squared of particle i and j, defined in the partonic centerof-mass frame as (2.14) Furthermore The role of F 2 is to handle the singular region associated to all partons having small transverse momentum. It also plays the role of increasing the importance sampling in the region of large transverse-momentum jets, a feature that is needed in order to properly cover a large range of transverse energy.

MiNLO
If we apply the MiNLO procedure to the trijet generator, the factor F 1 discussed above is no longer needed. All singular regions, except for the overall one, are regulated by the MiNLO Sudakov form factors. The MiNLO form factor is exactly as described in ref. [16], with the only freedom of choosing the scale of the basic process, 3 that in the case of ref. [16] (dealing with Higgs boson production in association with jets) was taken equal to the Higgs boson virtuality. In the present case, we have chosen the scale of the basic process to be equal to the sum of the transverse momenta of the two final pseudoparticles (after the MiNLO clustering has taken place).

Checks of the code
We have performed several checks on our code. We have generated subroutines for the virtual corrections using NJET [18], GoSam [19,20,21,22,23,24,25] and HELAC-NLO [26] and compared them to the virtual contributions obtained with the routines in our program. We have found that the NJET, GoSam and HELAC-NLO results were in agreement among each other for all subprocesses. We also found agreement with the NLOJET++ routines, except for the qqggg amplitude, and all its crossings, where there was a problem in the colour sum. After fixing this, we found perfect agreement. The Born contribution, together with the colour-and spin-correlated Born amplitude, were generated using the MadGraph4 POWHEG BOX interface. This allowed also the generation of the real contribution according to MadGraph4, that was thus checked against the (much faster) one that we have implemented in our code.
As a further check, the full NLO calculation was compared with the one of ref. [27] by first comparing the plots displayed there with the results of our code, that we ran using the same parton distribution functions and the same scales. We found agreement within statistical errors. In order to carry out this comparison we needed to run the POWHEG BOX with the same dynamical scale choice of ref. [27]. Although the default scale choice in the POWHEG BOX is computed as a function of the underlying-Born kinematics, it is possible to set it up in such a way that any scale choice can be implemented. See appendix A for more details.

Comparing trijet results at various levels
The trijet generator can be used to compute three-jet observables at the NLO level, at the POWHEG Les Houches Event level (i.e. with the hardest emission generated according to the shower technique implemented in POWHEG), and after the full shower. The MiNLO feature can be turned on already at the NLO stage. In the present section we compare the output of our generator at these levels.

General settings for the forthcoming comparisons
We consider jet production at the 8 TeV LHC. We use the CT10nlo pdf set [28]. We remind the reader that any of the modern pdf sets can be used [29,30], and our choice is only a matter of definiteness. We consider the shower output at the parton level, without the inclusion of hadronization effects and multiple interactions, since, at this stage, we are not comparing our result with jet data. 4 We interface our generator to PYTHIA (version 6.4.25), Pythia8 (version 8.183) and HERWIG (version 6.510). When using PYTHIA we use the Perugia 0 tune [31] (pytune(320)). We use Pythia8 and HERWIG with their default tunes.
In the POWHEG settings we have included the doublefsr 1 option, and the modification of the scalup prescription obtained by setting changescalup 1. Both these features are illustrated in ref. [32].
We encountered severe stability problems when using HERWIG, showing up as spikes in our final histograms. We investigated them, and found that they were related to events having very small transverse energy (of the order of few GeV) at the Les Houches level, and developing very high transverse momentum jets (above 50 GeV) after shower. The cause of these problems was photon emission from quarks, that apparently does not comply with the scalup veto in HERWIG. These problems disappeared completely by setting vpcut=1D30, that switches off photon radiation. Thus, all our HERWIG results were obtained with this setting. We verified that it has no visible effect on our results, but for the disappearance of the spikes.
The use of PYTHIA with the trijet generator in combination with MiNLO requires particular care. In fact, in this case, PYTHIA is unable to shower a sensible fraction of very small H T events. These events are not going to contribute to physically interesting distributions. Thus, they should be treated as events that did not pass the cuts, and should be counted when dividing the total weight entering a histogram bin by the total number of events. It turns out that, when PYTHIA finds such events, it silently discards them and loads a new one. In the analysis setup of the POWHEG BOX, the total number of events is usually computed as the number of times that the analysis routine is called, and, because of the aforementioned behaviour of PYTHIA, the fraction of discarded events is thus not counted. We coded a workaround for this problem in our analysis routines. A user adopting different analysis frameworks must make sure not to incur this problem.
We adopt the following values for the parameters entering the F function in eqs. (2.11) and (2.13): In addition, in order to avoid uninteresting regions with very small transverse momenta that may cause numerical problems, we reject all Born configurations such that min{q i , q ij } < (0.3 GeV) 2 . The value of the factorization and renormalization scale is taken equal to H T /2, computed on the partonic configuration of the underlying-Born kinematics. Finally, jets are reconstructed using the anti-k T algorithm [33] as implemented in the Fastjet package [34,35], with jet radius R.
The results we will show in the next sections have been obtained by generating 2.4 M events in two runs, with and without MiNLO. The runs have been performed on a 48 core machine, and they took roughly 37 and 100 hours, respectively.

NLO, LHE and shower-level comparisons
We begin by showing in fig. 1 the comparison of the fixed-NLO and the LHE-level results for the transverse-momentum distribution of the third jet. The PYTHIA and HERWIG showered output, compared to the NLO result, are displayed in figs. 2 and 3.
We first remark that the small transverse-momentum suppression of the LHE, PYTHIA and HERWIG showered results is simply due to the fact that, because of our F function, events with transverse momentum smaller than 50 GeV are very rarely generated.
We observe that when R = 1 all results are in better agreement. Differences arise for smaller values of R, and can be understood as follows. First of all, it can be easily checked  that the LHE level result has very mild dependence upon R. This is due to the fact that the splitting of the third parton into two collinear partons has a very strong Sudakov suppression. In fact, in this case the POWHEG Sudakov form factor is the product of the form factors for vetoing harder final-state splittings of all final-state partons, times the form factor for vetoing harder initial-state radiation. Because of this suppression, partons are relatively well separated, and a small R dependence is observed. 5 When completing the shower, further splitting processes can take place at a suitable rate, and the R dependence is reinstated. Notice that no Sudakov suppression for radiation is included in the fixedorder calculation, yielding a visible R dependence. It is clear, however, that in the trijet generator the R dependence will mostly arise at the shower stage. This is a desirable feature. In fact we do not expect the NLO result to be reliable in this region, since, among other things, it lacks the Sudakov form factor and the appropriate scale choice in the coupling constant. Because of this, the shower algorithm will acquire the responsibility to  reliably describe the R dependence.
In figs. 2 and 3 we observe a disturbing difference between PYTHIA and HERWIG. In the latter, the R dependence is much stronger. Our expectation is that the shower result should be determined by two elements: on one side, the introduction of the correct Sudakov form factor for the splitting process (that will tend to reduce collinear splitting and thus increase the third jet cross section at smaller R), and multiple emissions, that will tend to increase collinear splitting processes, and thus reduce the jet cross section at smaller R. The net effect is an increase of the shower cross section at R = 0.5 (with respect to the R = 1 value) that is around 10% in PYTHIA, but is more of the order of 20% in HERWIG. Furthermore, the HERWIG result shows an increasing discrepancy with the fixed order result at large transverse momenta of the third jet in the R = 1 case. We have no good understanding of why this is the case. On the other hand, the Pythia8 result is compatible with the PYTHIA one.
When going from a fixed-NLO result to an LHE one, the most striking differences are  The 100 GeV cut on the third jet is imposed for the following reason. If no cuts on the remaining jets are imposed, as p j 4 T increases, also p j 1 T , p j 2 T and p j 3 T must increase, and we are thus probing an overall property of the cross section. With the cut on the third jet, by studying the p j 4 T spectrum below the third-jet cut, we are studying the soft-collinear radiation dynamics from the three-jet Born configuration. We thus expect, for example, that the fixed-order result (that for this quantity has only leading order accuracy) diverges at very small transverse momenta, in contrast with the LHE result, where the soft-collinear region for the emission of the fourth jet is strongly Sudakov suppressed. We also remark that, in this case, the R dependence of the NLO result is not at all reliable, since no further partons are emitted in this framework beyond the fourth.

Comparison between the NLO results with a standard choice of scales and NLO+MiNLO
Before turning to the discussion of the trijet results when MiNLO is active, we perform a comparison of the bare NLO calculation with or without MiNLO, whose purpose is to show that only minor differences are seen in the region where the three jets are resolved. Here we report in fig. 6 such comparison for the transverse-momentum distribution of the third jet. Our standard scale choice, when MiNLO is not used, is to set µ R and µ F equal to H T /2, where H T is computed on the kinematics of the underlying-Born configuration. We see that the two results are fairly compatible, except in the very small transverse-momentum region. Here, the NLO with the standard scale choice grows much faster in magnitude, and becomes large and negative in the first bin. With MiNLO, the small transverse-momentum region is better behaved, as expected. However, we remind the reader that, in fig. 6, this region is divergent also with MiNLO, since it is dominated by the production of low transverse-momentum jets.

NLO, LHE and shower-level comparisons with MiNLO
It is interesting to present the results for the transverse momentum of the third jet, when the MiNLO feature is turned on. Also in this case, for R = 1 we see good agreement among the NLO, LHE and the showered results. We illustrate these results in figs. 7, 8 and 9. It turns out, however, that for R = 0.5 the NLO result is below the LHE one by a larger amount with respect to the no-MiNLO case, which leads to a slightly better agreement of the showered and NLO results.

The low transverse-energy region
When showering a POWHEG-generated LHE configuration, we usually assume that the jet structure of the event is only marginally affected by the shower. In jet production, in particular, we assume that the configurations having very small transverse energy (i.e. H T ) at the Les Houches level should not contribute significantly to events that pass the jet cuts.  It turns out, however, that such configurations have diverging cross section, and thus one may worry that the very small probability that the shower has for building up relatively hard jets starting from LHE configurations with small H T , may end up being amplified by an unphysically large cross section. A similar problem arises when we consider associated jets in a hard phenomenon. In this case, however, the MiNLO procedure ensures that the cross section for the Les Houches event is physically well behaved. On the other hand, in the case of jet production, the Sudakov resummation is not enough to guarantee a physical behaviour at small transverse energies.
In order to study this potential problem, we have taken a very simple approach. We have determined the cross section for events with transverse energy above a given cut at the Les Houches level. We have then found that for H T > 10 GeV such cross section is about 60 mb, going down to around 30 mb for H T > 20 GeV. Without imposing this cut, the cross section reaches 1000 mb. The diverging behaviour is, in fact, limited by the tiny cut-off that we impose upon the kinematics to avoid divergences. Since the inelastic cross section at the 7 TeV LHC is around 70 mb [36], it seems reasonable to impose a cut on the transverse energy of our events at the Les Houches level, in the range between 10 and 20 GeV. We have found that, for the shower Monte Carlo programs that we have considered, and with our settings, the Les Houches level cut has visible impact only on events with very small transverse momenta. For example, we find sensible differences in the distribution of the transverse momentum of the third jet only for p j 3 T 5 GeV, for the cut H T > 10 GeV, and for p j 3 T 10 GeV for the cut H T > 20 GeV. This is reasonable to expect, since the H T of the event is at least three times the transverse momentum of the third jet. In spite of this, we have preferred to maintain the Les Houches level H T > 20 GeV cut as our standard, since it seems, in all cases, unreasonable to have events with a cross section that becomes of the same order of the total inelastic cross section.

Comparison of the trijet+MiNLO results with the dijet results for inclusive quantities
When using MiNLO, the trijet generator becomes predictive also for inclusive quantities, i.e. for observables that do not necessarily require the presence of a third jet. It is hard to quantify theoretically the accuracy of such predictions. We have shown that, for inclusive quantities, the VJ generators (i.e. generators for Higgs or W/Z boson production in association with a jet) improved with the MiNLO procedure yield an accuracy that is better then LO, but not quite at the NLO level, unless one makes a careful tuning of the procedure [37].
In the case of jets, the argument of ref. [37] cannot be applied as is, since the soft-singularity structure of the two-parton production process is quite involved. Rather than trying to understand theoretically what is the accuracy of the trijet+MiNLO generator for inclusive quantities, here we simply compare its inclusive distributions to those obtained with the NLO-accurate dijet generator. We begin by showing in fig. 10 the transverse-momentum distribution of the third jet, p j 3 T , in events where p j 2 T > 20 GeV. The aim of the figure is to compare the Sudakov effects affecting the production of the third jet, introduced by the POWHEG machinery, in the case of the dijet generator, with those introduced in the trijet generator by the MiNLO procedure. We see that, in both cases, the small transverse-momentum region is properly regulated by the Sudakov form factor. We remind the reader that, as far as the large transverse-momentum region is concerned, the trijet generator has NLO accuracy for this observable, while the dijet one is only LO accurate. The slightly strange features at very low transverse momentum, that we observe with the Pythia8 shower, are shared for the same observable by the PYTHIA result, while, in the HERWIG case, we see a smoother behaviour, as shown in figs. 11 and 12.
We now turn to the comparison of the trijet+MiNLO and dijet generators for inclusive-jet distributions. We first compare the rapidity distribution of an inclusive jet in fig. 13. No interesting differences are seen between the trijet+MiNLO and dijet curves, other than the obvious difference in normalization, due to the cut in transverse momen-  tum, that could also be evinced from the transverse-momentum distribution. We show this distribution only for the case of Pythia8 shower, since we obtain similar results when using HERWIG or PYTHIA.
Another interesting distribution is the transverse-momentum of an inclusive jet, displayed in figs. 14, 15 and 16 for Pythia8, PYTHIA and HERWIG respectively. We notice the remarkable agreement between the trijet+MiNLO and dijet generators.
In conclusion, we find that the use of MiNLO considerably improves the trijet generator also in the region where one jet becomes unresolved. Observe that, while the dijet generator is NLO accurate for these distributions, the trijet+MiNLO generator is at most LO accurate. Thus, we do not advocate the use of the trijet+MiNLO generator for one or two jet inclusive distributions. However, the fact that also these regions are treated consistently gives us confidence that we should not be afraid to make predictions for three-jet observables also in the region where one jet is relatively close to being soft, or relatively  close to a collinear configuration. We thus recommend that our generator is used with the MiNLO feature turned on.

Scale-variation bands
In this section we show the scale-variation bands for some key distributions. The purpose of this section is threesome: to show that scale uncertainties can be easily computed, to give an idea of the uncertainty involved, and to show that the remarkable agreement of the trijet+MiNLO and dijet results for inclusive quantities is not accidental. A more thorough uncertainty study will be carried out in a forecoming publication where comparisons with available data will be considered. We illustrate in figs. 17, 18 and 19 the scale-variation bands for the transverse-momentum distribution of the third hardest jet, and the rapidity distribution and transverse momentum of the inclusive jet, for trijet+MiNLO at the LHE level. We compare the trijet+MiNLO scale variation with the same scale variation in dijet production.
The scale variation can be performed in a fast way, without regenerating the event sample, using the reweighting tool in the POWHEG BOX V2, that allows for a very fast re-  evaluation of the weight associated to each event. 6 The scale-variation band is obtained by taking the envelope of the 7 differential cross sections computed by multiplying the reference factorization and renormalization scales by the factors K F and K R , respectively, where In all the three distributions we are showing, we get a comparable scale-band size in the trijet+MiNLO and dijet results, of the order of 20%, and there is a very good degree of overlapping for the inclusive-jet rapidity and transverse momentum.

Conclusions
In this paper, we have presented an implementation of a NLO plus parton-shower generator for three-jet production, built in the framework of the POWHEG BOX V2.
We have compared key kinematic distributions at different levels: NLO, Les Houches event level, and after the shower performed by PYTHIA, Pythia8 and HERWIG6. We found very good agreement between the NLO and the PYTHIA and Pythia8 results, for variables that are correctly described by a fixed-order calculation. Slightly worse agreement is found between the NLO and the HERWIG-showered results.
We have also applied the recently-proposed MiNLO procedure, for the scale assignment in the NLO calculation, to our generator. We have found that MiNLO considerably improves the trijet generator also in the regions where one jet becomes unresolved. The fact that also these regions are treated consistently gives us confidence that we can make predictions for three-jet observables also in the region where one jet is relatively close to being soft, or relatively close to a collinear configuration.
We have seen that the trijet+MiNLO and dijet results display remarkable consistency among each other. On the other hand, we have also evidence that the kind of shower generator that is used for the final shower has a non-negligible impact on the result, especially for relatively-small jet cone sizes.
The code can be downloaded following the instructions in the POWHEG BOX web site http://powhegbox.mib.infn.it.
Fund grant K-101482 and is thankful to the Aspen Center for Physics for warm hospitality where part of this work was carried out. The computations shown in this paper were partially performed on the HPC facility of the University of Debrecen (NIIF Institute, Hungary).

A. Scale options in the POWHEG BOX
In the POWHEG BOX, the factorization and renormalization scales are usually set as a function of the underlying-Born kinematics. Since the POWHEG BOX can also be used as a parton level, fixed-order generator (by setting testplots 1 in the powheg.input file), it is convenient, at times, to remove this restriction (for example, in order to compare the fixed order NLO output to other codes). This is done as follows. If one sets the variable btlscalereal 1 in the powheg.input file, the internal flag flg btildepart is used to distinguish the Born, virtual and subtraction-term contributions from the real one. When flg btildepart equals 'b', the program is computing the Born or virtual. When it is set to 'r' it is computing the real contribution. The user can then modify the set fac ren scales subroutine, so that, on the basis of the value of flg btildepart, the program uses the Born or real kinematics to compute the scales.
Another ambiguity in the scale choice has to do with the computation of the subtraction terms. It is acceptable to use for them the same scales used for the real contributions. On the other hand, it is also acceptable to use for them the scales of the corresponding underlying-Born configuration. In order to implement also this option, one sets the variable btlscalect 1 in the powheg.input file. If this variable is set, the set fac ren scales subroutine is called with flg btildepart equals to 'c' when the subtraction terms are computed.