Introduction to the PanScales framework, version 0.1

In this article, we document version 0.1 of the PanScales code for parton shower simulations. With the help of a few examples, we discuss basic usage of the code, including tests of logarithmic accuracy of parton showers. We expose some of the numerical techniques underlying the logarithmic tests and include a description of how users can implement their own showers within the framework. Some of the simpler logarithmic tests can be performed in a few minutes on a modern laptop. As an early step towards phenomenology, we also outline some aspects of a preliminary interface to Pythia, for access to its hard matrix elements and its hadronisation modules.


Introduction
Parton showers lie at the core of the majority of experimental and phenomenological studies in collider physics.At the LHC, they connect the electroweak and TeV momentum scales of hard-scattering processes, where the relevant degrees of freedom are perturbative quarks and gluons, with the non-perturbative physics of hadrons at scales of a few hundred MeV.As such, parton showers account for physics across several orders of magnitude in momentum scales.In QCD, large logarithms typically appear in the presence of large momentum scale hierarchies, which have to be resummed to all orders in the strong coupling to obtain physically sensible results.One of the frontiers of the development of parton showers is to understand, demonstrate and improve their logarithmic accuracy, with analytic resummations providing crucial inputs, as well as reference results for comparison.
This PANSCALES release includes two main NLL-accurate parton showers, PanGlobal and PanLocal.They have had their next-to-leading-logarithmic (NLL) accuracy tested for e + e − [1] and colour-singlet production in pp [5,6] collisions, as well as Deep Inelastic Scattering (DIS) and Vector Boson Fusion (VBF) processes [8].They include state-of-the-art handling of subleading colour corrections, which for many processes and observables allows for full colour accuracy at LL (and often beyond) [2,5].They also include the treatment of both collinear and soft spin correlations [3][4][5], and first steps towards matching [7], as well as elements towards NNLL accuracy [9] (the latter two just for e + e − collisions).Finally, the codebase contains early versions of features that are yet to be discussed in physics research papers, in particular an interface with the PYTHIA8.3 event generator [26], which can be used to provide hardprocess generation and hadronisation.Despite the inclusion of the interface to PYTHIA8.3, the code is not yet at a stage of maturity that is suitable for extensive comparisons to experimental data.This is notably because of the absence of finite quark-mass effects, the need for further work on matching with higher (fixed) order effects, as well as tuning of the nonperturbative parameters of the shower and of any hadronisation model with which it is used (see e.g.Refs.[27][28][29]).
This manuscript is structured as follows.Section 2 focuses on basic usage of the code, illustrating: the build procedure (section 2.1); stand-alone event generation (section 2.2); the use of the code for carrying out basic logarithmic tests of parton showers (section 2.3); usage with a preliminary interface to PYTHIA8.3 (section 2.4); and details for carrying out validation of the code and building it with higher numerical precision (section 2.5).Section 3 illustrates some of the techniques that underlie the logarithmic tests, discussing both double-logarithmic global event shape observables (section 3.1) and single-logarithmic non-global observables (section 3.2).Section 4 gives a brief discussion of how to use the PanScales framework to implement a new shower, which provides a relatively straightforward way to gain access to the colour, spin-handling and logarithmic-accuracy testing facilities.We close in section 5 with an outlook.

Basic usage
The PANSCALES code requires a C++14 compiler and a Fortran 95 compiler, the GSL library, CMake (≥ 3.7) and, for some scripts, Python (≥ 3.6) with matplotlib installed.Some features (higher-precision builds) require the MPFR [30] and QD [31] libraries (see section 2.5 for details).It includes several third party codes, notably fjcore [32] for jet finding and hoppet [33] for PDF handling and the Catch2 library for unit-testing (see the 3rdPartyCode.mdfile for further details).The code can also be linked to LHAPDF [34] and to PYTHIA8.3.The PANSCALES code is released under the GNU GPL v3 license.

Downloading and building the code
PANSCALES can be obtained from the git repository git clone --recursive https://gitlab.com/panscales/panscales-0.X This comes with the following directory structure: • shower-code/: main shower code and some simple examples • pythia-interface/: code for the PANSCALES interface to PYTHIA8.3 [26] • analyses/nll-validation/: code for the full set of NLL validation tests • helpers/: contains the double_exp numerical type and a copy of fjcore [32] • submodules/: contains HOPPET [33], CmdLine (for command-line processing and help) and AnalysisTools (for run steering and histogramming) • scripts/: a variety of helper scripts • unit-tests/: code for low-level unit tests To build the code and examples in double-precision, do the following cd panscales-0.X/shower-code ../scripts/build.py-j The scripts/build.py script uses CMake to organise the build, which by default is placed in the build-double/ subdirectory.Advanced use of CMake is described in the BUILD.mdfile.

Standalone event generation
To run showering for the e + e − → qq process and analyse the events, do build-double/example-ee -shower panglobal -beta 0 -process ee2qq \ -physical-coupling -rts 91.1876 -nev 100000 \ -out example-results/example-ee.dat This will take about 5−10 seconds, and produce an output file with histograms for a range of event shapes.Runs are generally configured with command-line options, for example with the first line indicating the use of the PanGlobal shower.The shower ordering variable is v = k t e −β PS |η| , with k t and η = − ln tan θ /2 respectively the transverse momentum and pseudorapidity of the emission with respect to its parent.The choice β PS = 0 (set with the -beta 0 command-line argument) corresponds to transverse-momentum ordering.
The available command-line options can be explored with the -h flag, and command-line options can also be placed in a card-file and read with the -argfile card-file.txtoption.For convenience, the main options are also listed in an OPTIONS.mdfile (any executable can be made to generate such a file by adding -markdown-help to the command line).More details on the PANSCALES interface can be found by examining the example-ee.ccfile.Much of the code also has doxygen documentation, which can be obtained by running doxygen from the shower-code/ directory.Note that in standalone mode, as given above, currently all events have unit weight, i.e. in order to recover a physical cross-section from the above run the histograms should be multiplied by the appropriate cross section for the hard process.Furthermore, in most cases the Born event is a fixed configuration rather than being sampled over.
A final comment is that the event record (in the class panscales::Event, Appendix A.1) holds only the particles as they appear after all showering, rather than containing all intermediate steps as in some other codes.The reason for this is that for showers with global recoil, storing all intermediate steps would take memory of order n 2 for an n-particle event. 1 Access to the event record at intermediate steps can however be enabled, cf.Appendix A.2.By default the code prints the first event, but this can be changed to print N events with the -max-print Similar examples for pp → Z and DIS are given respectively in the example-pp.ccand example-dis.cc,with further explanations in an EXAMPLES.mdfile.Note that the above examples do not by default integrate over the hard-process kinematics.That functionality is instead available through the interface with PYTHIA8.3 (cf.section 2.4).

Logarithmic tests
Within the shower-code/ directory, the example-global-nll-ee.pyscript, and associated example-global-nll-ee.ccprogram, illustrate how to carry out a basic test of the NLL accuracy of a shower for global event-shape observables.It can be used as ./example-global-nll-ee.py --njobs NJOBS --shower panglobal and will test the NLL accuracy of the e + e − PanGlobal shower with β PS = 0 (i.e.k t ordered), for event shapes like the Cambridge-algorithm y 23 [35] and Lund observables [1].The Lund observables measure either the sum or the maximum of the k t i e −β obs |η i | for primary Lundplane declusterings i [36].The NLL test works with β obs = 0, i.e. examining just the relative transverse momenta (k t ) of the declusterings.
The --njobs flag takes an integer which should be equal to the number of cores that one wishes to run on.Depending on the machine, the script takes around 15 to 50 CPU minutes, i.e. a few minutes of wall-time on a modern multi-core machine.On completion, for each observable, it will produce plots such as those in Fig. 1, showing the α s → 0 limit of the ratio Σ shower (λ)/Σ NLL (λ), where Σ(λ) is the cross section for α s ln O to be smaller than λ, with O being the value of the observable.The plot is given as a function of λ, as is standard for NLL tests [1]: for an NLL-correct shower, the α s → 0 limit of the ratio will be equal to 1.The user can examine the results for a non-NLL shower by replacing panglobal on the command line with dipole-kt, which provides a standard k t -ordered dipole shower [5], much like those [37][38][39] in standard public tools.
The shower-code directory also includes an example-nonglobal-nll-ee.py script for testing non-global logarithms in the context of energy flow into an angular slice.It can be used as follows ./example-nonglobal-nll-ee.py --njobs NJOBS --shower panglobal The script again generates a plot with the ratio to the NLL (single-logarithmic) result.
The above command lines serve mainly to illustrate the more straightforward logarithmic accuracy tests and provide only a subset of the functionality required for a full set of tests.
For example, the global event-shape tests in the example-global-nll-ee.pyscript are limited to β PS = β obs = 0. Further discussion of considerations for logarithmic tests is given in section 3.

Usage with Pythia8
While we do not yet recommend the PANSCALES code for phenomenological production purposes, those wishing to start exploring such applications may try the interface with the PYTHIA8.3 [26] generator code, which we have tested with version 8.3.10.This enables the use of PYTHIA8.3 to generate the hard process,2 as well as for hadronisation, while PANSCALES carries out the parton showering at the scales in between.It also enables access to PYTHIA8.3 I/O, e.g. for outputting HepMC files [40].
To compile and run the code, enter the pythia-interface/ directory and run 3./get-pythia.sh../scripts/build.py--build-lib -j --with-lhapdf The main example program simulates Drell-Yan production.The default invocation requires the CT14lo set [41] from LHAPDF [34], which needs to be installed on the system e.g. by running lhapdf install CT14lo. 4Running build-double/main-dy -physical-coupling -lhapdf-set CT14lo \ -shower panglobal -nev 1e5 -out main-dy.datwill generate 10 5 Drell-Yan events and histogram the Z rapidity, mass and transverse momentum.It does not (yet) include matching to higher-order matrix elements, so kinematic distributions such as the transverse momentum of the colour singlet p t Z are sensible only in the resummation region, i.e. at low p t Z .Note that the above run will produce warnings concerning x regions where the PDF set is badly behaved.Most other PDF sets have issues with negative or zero parton distribution functions, especially at large x, that cause the PANSCALES code to throw an exception after some number of events.Ultimately, we intend to make PANSCALES more tolerant of ill-behaved PDFs, but at this stage of development we have taken the approach that it is safer to abort the run than to continue generating events when a clear problem has arisen, even if only a rare occurrence.
The same directory contains a range of other examples that can be run with PYTHIA8.3, and the header of each example illustrates how to use it.We also provide an example of interfacing to RIVET [42], in which case one needs to make sure RIVET is installed and the code is compiled with it (through -cmake-options="-DWITH_RIVET=on").This allows for an easy comparison with data, although we should note that none of the PANSCALES showers are currently tuned, nor do they include the effects of quark masses.We comment further at the end of this subsection.
The PANSCALES-PYTHIA8.3 interface transfers the event into the PYTHIA8.3 event record after each shower emission.This is done to have access to the full PYTHIA8.3functionality in the future, i.e. interleaving multi-parton interactions (MPI) with showering of the hard process.Note that at this moment, we do not have the functionality to run with MPI, but hadronisation can be added through the flag -hadron-level.
The above examples all use the same PANSCALES event-loop framework as in the main shower-code/ directory.We also distribute examples with a standard PYTHIA8.3structure.These are to be found in the main-pythia02.cc and main-pythia06.ccfiles, which are based on the corresponding main02.ccand main06.ccexamples from the PYTHIA8.3 distribution.Note that only a restricted set of options is supported in this form, which can be found at the end of PanScalesPythiaModule.cc.
In thinking about the interface between the PanScales showers and Pythia's hadronisation, it is perhaps useful to comment briefly on the interplay between hadronisation and the parton shower's logarithmic resummation.Consider an event with a hard scattering scale Q and some infrared and collinear (IRC) safe observable (e.g. a jet rate) that is sensitive to a transversemomentum scale p t ≪ Q.Any discussion in terms of quarks and gluons implicitly comes with a requirement that the transverse momentum scale p t be much larger than Λ QCD ∼ 1 GeV.We expect that hadronisation (and associated tuning parameters, such as those of the Lund string model) will act primarily at scales of the order of Λ QCD .Thus for the observable at scale p t , the relative impact of hadronisation should be suppressed by a power of Λ QCD /p t .Consequently hadronisation will not modify the resummation properties of the shower.Of course in practice one may also be interested in physics at a few GeV or in quantities that are not IRC safe.In those cases, the interplay between shower resummation and hadronisation may be nontrivial: for example more accurate shower resummation may provide a better initial condition for hadronisation models.A more detailed study of these questions would be interesting, but is beyond the scope of this manual.

Code validation and more advanced builds
To validate that the code is generating expected results, enter the shower-code/validation/ directory and run ./validate-showers.py -j which runs a range of validation tests in parallel across all available cores.This takes a total of about 300−1000 CPU seconds on a typical laptop, carrying out of the order of 100 separate short runs, each with different settings, verifying that they give histograms that are identical to those in a set of reference files.It is possible to carry out validation runs with larger numbers of events, but one should be aware that there can be differences due to varying floating-point behaviours across different hardware.Additionally, lower-level unit tests can be found in the unit-tests/ directory.These are based on the Catch2 library and are mainly used for checks of foundational classes, such as the main Momentum class (see the README.mdfile in that directory for further details).The unit tests and a subset of the validation tests can also be run with the continuous integration script scripts/CI-build-and-test.py.
It is sometimes useful to build the main code in different numerical precisions, e.g. for logarithmic tests that probe very disparate energy scales and angles.For this, the general build script has an option -builds X which essentially invokes cmake with a suitable set of configuration options.Specifically, one runs

../scripts/build.py --builds X [-j]
where X is a space-separated list that contains one or more of the following options: double, ddreal, qdreal, doubleexp, mpfr4096.
The ddreal and qdreal options require at least version 2.3.23 of the QD library [31] and have precisions of about twice and four times that of a double type, with speeds that are about 10 and 200−300 times slower.
The doubleexp type was developed specifically for the PANSCALES project and has the same relative precision as double, but a much larger range of exponents (stored in an additional 64-bit integer), which is useful when exploring finite values of α s ln v with very small values of α s and correspondingly large values of ln v.It is about 3−10 times slower than double.
The mpfr4096 type is based on the MPFR library and has 4096 bits for the mantissa, i.e.
about 75 times higher than double precision, or equivalently a little over 1000 decimal digits of precision.It is about 1500−2500 times slower than double (using version 4.2.1 of the library).
For most purposes, the double and doubleexp types are sufficient, notably when used with methods that track the differences between directions of dipole ends in addition to the actual momenta of the corresponding partons [2].That tracking can be enabled at run time with the -use-diffs option and has only a modest speed penalty, of the order of 10%.The higher precision types are, however, important when developing new showers and testing the correctness of any parts of the code that carry out the dedicated calculations with direction differences.

Further details for logarithmic tests
In this section we provide some insight into features of the code that facilitate shower logarithmic accuracy tests, together with a more detailed discussion of some of the underlying methodology than has been given in previous work.Some of the discussion below concerns tests that go beyond the simple ones illustrated in section 2.3.Code for these more advanced tests is to be found in the directory analyses/nll-validation/ with usage explanations in the corresponding README.mdfile.

Global event shapes
We start by considering the global event shape tests of section 2.3 and specifically examine the commands that are run by the example-global-nll-ee.pyscript.For each of several α s values, that script executes one or more commands of the following kind build-double/example-global-nll-ee -Q 1.0 -shower panglobal -beta 0.0 \ -alphas 0.0016 -nloops 2 \ -lambda-obs-min -0.5 -lnkt-cutoff -327.5 -dynamic-lncutoff -15 \ -weighted-generation -nev-cache 75000.0\ -spin-corr off -use-diffs -nev 750000 -rseq 11 \ -out example-results/lambda-0.5-alphas0.0016-rseq11.res The second line indicates the value of α s (Q) in the MS scheme (working in units of s ≡ Q = 1) and that a two-loop running coupling is to be used.The next two lines contain two separate critical elements, which we discuss below in more detail.Further options are -spin-corr off, which turns off spin correlations (leaving them on adds about 50% to the run time).The -use-diffs option turns on tracking of direction differences for higher numerical precision.It is not strictly necessary for this example, but the speed cost is relatively limited, at about 10%.The -rseq 11 argument specifies the random sequence that is used.
The choice of α s values (cf.also Fig. 1) involves a balance between several considerations.
1. CPU usage differs across values of α s .For example taking too small a value of α s (and correspondingly large logarithms) requires the use of the doubleexp type, which is slower.Furthermore, in some cases, smaller α s values bring slower Monte Carlo convergence (notably when one needs to use the methods of Section 3.1.3).

At order N p LL one is looking for residual effects of O (α p−1 s
).For NLL tests, we have p = 1 and so we are measuring an O (1) effect, i.e. independent of α s .But for, say, an NNLL test, one would be looking to measure O (α s ) effects, and taking values of α s that are too small would make it more difficult to identify those effects.
3. The ability to take the shower cutoff sufficiently small that cutoff effects are negligible, while also not encountering the Landau pole in the coupling.This complicates the use of larger values of α s , because one approaches the Landau pole more rapidly.CPU usage then increases significantly, notably because of the higher particle multiplicities in each event.
4. The size of residual contributions that are beyond the order of extrapolation for the α s → 0 limit.This favours taking smaller values of α s .
As concerns the last of these points, one should estimate a systematic error for the extrapolation, for example by using two different sets of α s points.While this is not done for Fig. 1, it is a part of all tests in analyses/nll-validation/ that rely on extrapolation of α s → 0.
Two further critical elements that we now need to discuss in more detail are (1) a dynamic generation cutoff, which is necessary to prevent event particle multiplicities from becoming intractable and (2) weighted event generation, which is necessary in order to explore regions with large Sudakov suppression.

Dynamic generation cutoff
Let us discuss the following part of the command line: -lambda-obs-min -0.5 -lnkt-cutoff -327.5 -dynamic-lncutoff -15 This indicates that the minimal value of λ = α s ln v is −0.5, which corresponds to ln k t /Q = −312.5 with α s = 0.0016.The choice -lnkt-cutoff -327.5 allows the shower to run somewhat below the k t scale associated with the minimal value of λ.This is important, because multiple emissions break the immediate relation between the shower scale and the observable, e.g. because of shower recoil effects, and because the observable sums contributions from multiple emissions.These effects are properly accounted for only if the shower is allowed to evolve sufficiently far below the single-emission scale that is equivalent to the smallest value of the observable that is probed.
If we were to straightforwardly run with the very small cutoff indicated above, then the average event particle multiplicity, n, would be prohibitively large.Specifically, it scales as [43], where λ c = α s ln k t,cutoff /Q.For the parameters above, this would give n ∼ 10 8 .To resolve this issue, in addition to the fixed k t cutoff we also use a dynamic cutoff.We track the shower scale v 1 of the first shower emission.Knowing that we are restricting our attention to β obs = β PS observables and that standard global event shapes are recursively infrared-and collinear safe [44], we can be sure that emissions with much lower values of v will not modify the event shapes.Accordingly, when the shower reaches an ln v scale that is sufficiently far below that of the first emission, we simply stop showering.That choice is specified by the -dynamic-lncutoff -15 argument, which stops the showering when ln v < ln v 1 + δ ln v = ln v 1 − 15.This ensures that the multiplicity is kept limited, e.g. for the above run it averages to about 20.
To validate this approach, the dependence of the normalised cross section Σ(λ) on the size of the dynamic cutoff is shown in Fig. 2a (red line, corresponding to a k t -like observable, i.e. β obs = 0).It is illustrated for the pp → Z process, as carried out for the studies of Ref. [6], and one observes good convergence of Σ(λ) as the dynamic cutoff δ ln v is taken below −9 (for the β PS = β obs case discussed, δ ln v ≡ δ ln O in the plot).The value of −15 that we use above leaves a comfortable margin.
The above procedure is sufficient when the angular power β PS in the definition of the shower ordering variable coincides with the angular power β obs in the observable.When β PS ̸ = β obs a more elaborate procedure is needed.As a starting point, we need to be able to evaluate the order of magnitude of the contribution of each emission i to the observable, O approx,i .Given knowledge of β PS and β obs , this can be determined from the value of ln v and the (approximate) rapidity of the emission.Specifically we take it to have the form where ln b i is a rapidity-like variable (cf.Section 4), the c v and c b are coefficients that are independent of the dipole kinematics, and the remaining constant term does depend on the dipole kinematics. 5In general Eq.(1) will agree with the true contribution to the observable from emission i to within a factor of order 1.We then determine if the following condition holds ln O approx,i < max where on the right-hand side the max operation runs over prior accepted emissions j and δ ln O generalises the δ ln v discussed above.If Eq. (2) holds, then emission i is discarded and showering continues.Showering is subsequently stopped when ln v is sufficiently small such that for all emission rapidities one can be sure that Eq. (2) will always hold.
Note that for β PS ̸ = β obs this procedure is guaranteed to be safe only for showers that respect the PANSCALES condition that a given emission does not impact other emissions far in the Lund plane. 6Results are illustrated in Fig. 2, for various combinations of β PS and β obs , showing the relative change in Σ between successive pairs of values of δ ln O , corresponding to the extremities of each horizontal bar.One sees a behaviour that is consistent with an exponentially vanishing effect as δ ln O becomes more negative.Again, a choice of −15 should be more than adequate for NLL logarithmic tests.Note that the need for a more sophisticated dynamic veto and cutoff is not the only challenge that arises with β PS ̸ = β obs .Further considerations are discussed below in section 3.1.3.

Weighted event generation for β PS = β obs
Now we turn to the following part of the command line at the beginning of section 3.1 -weighted-generation -nev-cache 75000.0This is needed to address the fact that Σ becomes infinitesimally small in the limit α s → 0 for fixed λ = α s ln v, specifically ln Σ ∼ λ/α s .To get a more concrete sense of the challenge, consider that Σ ∼ 10 −103 for λ = −0.5 and α s = 0.0016 as in the command line at the start of section 3.1.With unweighted event generation, it would take orders of magnitude longer than the age of the universe to explore that region.
We address this challenge by greatly enhancing the number of events whose first emission has an extremely small value of ln v, assigning a suitable weight to those events so as to reproduce the correct final Σ distribution.
We divide the full evolution range [ln v max , ln v min ] into a set of n consecutive bins, each defined by their upper boundaries, ln v + i .Writing the shower Sudakov form factor between two scales v and v ′ as ∆(v, v ′ ) (for the Born event), we precompute the part of the Sudakov form factor associated with each bin. 7The precomputed Sudakov can be inspected in the output file.
For each event, we choose a generation bin i randomly, with a probability p i that we take proportional to ln v + i − ln v + i+1 .The shower then starts from scale v + i .If the first emission scale is above v + i+1 , the shower continues down to the dynamic cutoff as explained in section 3.1.1.If the first emission scale is below v + i+1 , the emission is discarded and showering is restarted from scale v + i , repeating this until one generates an emission above v + i+1 .The weight assigned to the event is where the factor 1/p i compensates for the likelihood of choosing the window and the second factor accounts for the Sudakov form-factor probability of having the first emission between v + i and v + i+1 .In practice we choose the bins such that ln v + i scales as − i, which in a fixed-coupling approximation ensures that the Sudakov form factor decreases by a similar factor from one bin to the next.We take the total number of bins to be n = ln ∆(v + , v − )/ ln(u) where u is a number of O(1/2), so that the probability of each successive bin is about half that of the previous bin.This ensures that each attempt at starting the showering has a ∼ 50% chance of generating an emission within the given window, while also ensuring that the event weight tracks the actual first-emission Sudakov probability to within about a factor of two.Fig. 3a shows a validation of the correctness of the procedure for a physical value of the coupling, α s = 0.1, where it is straightforward to obtain high accuracy with both weighted and unweighted approaches.The plot shows the distribution of ln v 1 for the first emission when the generation is unweighted (red solid lines) and weighted (blue dashed lines).The two approaches agree to within statistical errors.Fig. 3b shows the size of the relative statistical error in the two approaches, illustrating the superiority of weighted generation for small ln v 1 values, and also illustrating that its statistical error is fairly independent of ln v 1 , with a mild sawtooth structure that lines up with the generation bin edges (indicated as vertical dotted lines).

Weighted event generation for β PS ̸ = β obs
The above weighted event-generation procedure is very powerful when the Lund contour of the observable lines up with that of the shower evolution variable, i.e. β obs = β PS .If this is not the case, a given value of the observable receives contributions from different evolution windows with often widely differing weights, significantly worsening statistical convergence.
Let us start by explaining the structure of our approach, which is illustrated also with the help of Fig. 4. Ultimately, we want to generate all possible events such that a given observable O is below some threshold, e.g.α s ln O = λ < −0.5.As a first step will consider how to generate events with a condition on an approximation to the observable, O max ≡ max i {O approx,i }, with the O approx,i as defined in Eq. ( 1).Specifically, we will consider how to generate events such that O max is in some range O − < O max < O + , i.e. the green band in Fig. 4, labelled "approximate observable window".Then we will show how to use event-samples in different approximate-observable windows so as to obtain a set of events with a constraint on the actual observable, exploiting the fact that O and O max differ at most by a factor of order one.

C enh
Figure 4: Illustration of some of the main steps in the weighted event-generation approach that is used for β PS ̸ = β obs , for a specific target observable window.Beyond what is shown in the figure, one important element is that if, after showering, there are no emissions in the approximate observable window, the event is discarded.The approach is additionally supplemented with the dynamic generation cutoff of section 3.1.1.
The complication that we face in generating an event sample with a constraint on O max is that the shower contours (blue lines of Fig. 4) do not align with observable contours (e.g.edges of the green band).However, knowing the scaling of the observable, it is possible to analytically work out the largest shower scale, v a , that can generate an emission with approximate observable-value O + .We therefore start the showering from that scale v a , with an initial weight equal to w = ∆(v max , v a ), where v max is the largest kinematically accessible scale.We then need to ensure that the showering does not generate any emissions with O approx,i > O + .The simplest approach would be to veto every event that has any of the emissions i contributing such that O approx,i > O + .This would be correct, but would in general lead to a very small fraction of surviving events.That survival fraction corresponds roughly to the Sudakov form factor associated with the area between the top of the green band (O + ) and the thick blue v a line in Fig. 4.
To work around this problem, we employ a procedure that relates to an approach first introduced in Ref. [45] in the context of multi-jet merging.Specifically, we use an enhancement C enh > 1 for the probability of generating individual emissions.If an emission i has O approx,i > O + , the emission is discarded, the event is kept, but the event weight w is multiplied by a factor such that The final step is that for a given value of the constraint on the actual observable (e.g.α s ln O < λ = −0.5),we need to add together several contiguous O ± windows above the constraint, plus one final window without a lower bound (i.e.O − = 0).Note that Fig. 4 shows only one of the windows.The need for several windows is because the actual value of the observable can be smaller than the approximate value of the observable by some O (1) factor.The number of windows that is needed depends on both the shower and the observable.Typically if the actual observable can be substantially below the approximate observable, one needs a greater number of windows above the value of the true observable constraint.One should always check that the highest window contributes minimally to the result, i.e. less than the total statistical uncertainty after summing all windows.
Let us close with a comment on event weights.Within a given approximate observable window, from one event to the next Eq.( 4) will have been applied a different number of times, reflecting the number of emission that have been discarded in that event.There will therefore be a spread of event weights.The width of the event weight distribution gets smaller as C enh is made larger, but there is an associated slow-down of the showering, because of the increased emission probability.The optimal choice for C enh involves some balance between these two aspects.In practice we use C enh = 20 for showers with β PS = 0, and C enh = 100 for showers with β PS = 0.5.
In its current form, the code (analyses/nll-validation/shower-global-obs.cc)performs separate runs for each of several windows and the python script that manages the calculations (run-nll-tests.py)adds the results together.A command line that sets up the use of the above algorithm is, for instance build-double/shower-global-obs -Q 1.0 -no-spin-corr -nloops 2 -alphas 0.008 \ -shower panglobal -beta 0.0 -colour NODS \ -lambda-obs -0.5 -beta-obs 0.5 -lnkt-cutoff -77.5 \ -strategy RouletteEnhanced -enhancement-factor 20.0 \ -ln-obs-buffer 3.5 -nln-obs-div 7 -veto-buffer -15.0 \ -use-diffs -rseq 11 -nev 70000 \ -out nll-tests-tmp/panglobal00-fapx0.5-lambda0.5-as0.008-NODS-rseq11.dat The second and third lines set up a shower and an observable class with different β PS and β obs values.The fourth line sets up the strategy described above and the associated value of the enhancement factor, corresponding to C enh .The fifth line indicates that 7 approximate observable windows are explored, extending to e 3.5 times above the target value of the observable.The code performs runs in each of the 7 separate approximate-observable windows, producing one output file for each (Fig. 4 corresponds to a single window).The -veto-buffer -15 argument on the fifth line sets the size of the dynamic veto and cutoff, as explained for β PS ̸ = β obs at the end of section 3.1.1.If a user wishes to explore new observables or new showers, they are strongly advised to manually verify that the contribution to the cross section from the highest window is sufficiently suppressed relative to the total result.This can be done by examining the output files from the above command, there being one file per observable window. 8One way of understanding this is to think that one starts with a number of replicas of the shower that will evolve in parallel.The number of replicas is equal to C enh , and the enhancement of the emission probability can be interpreted as being equivalent to the total probability of emissions occurring in any of the replicas.Concentrating specifically on emissions with O approx,i > O + , the first time any of the replicas generates an emission with O approx,i > O + , that replica is simply discarded, leaving C enh − 1 replicas.The factor 4) is simply the ratio of surviving to original replicas.As the shower continues, the remaining weight C enh is then shared back out across C enh replicas again, and so the procedure continues.

Single-logarithmic observables, e.g. non-global logarithms
An example of a single-logarithmic test, such as the transverse energy flow in a rapidity slice in e + e − → qq collisions, can be performed with the example-nonglobal-nll-ee.py script in the shower-code/ directory.This will execute the following command build-doubleexp/example-nonglobal-nll-ee -Q 1.0 -shower panglobal -beta 0.0 \ -nloops 2 -colour CATwoCF \ -slice-maxrap 1.0 -lambda-obs-min -0.5 \ -alphas 1e-09 -lnkt-cutoff -501000000.0 -ln-obs-margin -11 \ -strategy CentralRap -half-central-rap-window 10 \ -spin-corr off -nev 10000 -rseq 11 \ -out example-results/lambda-0.5-alphas1e-09-rseq11.res The second line indicates that a 2-loop running coupling is to be used and that the colour scheme is a large-N C scheme in which C A = 2C F = 8/3. 9The third line indicates the size of the rapidity slice in which the in-slice energy flow E t will be measured, and sets the minimum value of λ = α s ln E t /Q = −0.5.
The fourth line indicates that we run at an infinitesimal value of the coupling, α s (Q) = 10 −9 (still working in units of s ≡ Q = 1), with a sufficiently small k t cutoff, e −501000000 , so as to cover the region down to λ = −0.5.The point of using such a small value of α s is to avoid a substantial extrapolation to α s = 0, alleviating the need for runs at multiple values of α s .It is physically possible to use such an infinitesimal value because the observable Σ(λ, α s ) is independent of α s for α s → 0. This is in contrast with the case of global event shapes where ln Σ(λ, α s ) ∼ λ/α s .The huge range of orders of magnitude of momenta requires the use of the doubleexp type.
As with the case of global event-shape tests, a straightforward run with the above parameters would give multiplicities that are much too high to be managed (in fact, the situation is even worse, because of the much smaller value of α s ).The mitigation procedure is different in this case: the combination of arguments on the fifth line causes the shower to only generate emissions within a window where the absolute rapidity is less than 11 (with respect to the emitter or spectator).For showers that satisfy the PANSCALES conditions, as long as λ is not too large and the window is large enough, the results should be (and are) independent of the window size.Corresponding arguments also exist for an analogous modification of the shower to generate only hard-collinear emissions, as is relevant for many spin-correlation, fragmentation function and PDF evolution tests.
Note that for verifying next-to-single-logarithmic accuracy [12,13,48] for non-global logarithms [9] it is necessary to enable the double-soft matrix element corrections and associated virtual corrections (with the -double-soft command-line argument, available only for the PanGlobal showers with β = 0 and β = 0.5, for e + e − collisions).One also has to modify the above command line to run at multiple small but finite values of α s , so as to be able to determine the first derivative of Σ(λ, α s ) with respect to α s in the α s → 0 limit.Doing so accurately requires considerable computer resources.Users who wish to explore this are advised, in a first instance, to run with the -split-dipole-frame option for the PanGlobal showers with β PS = 0, which is the most computationally efficient setup. 9Recall that the PanScales showers are logarithmically accurate for non-global logarithms only in the large-N C limit, though the subleading-colour schemes [2] are numerically close to the full-colour results [46,47].The schemes of Ref. [2] can be obtained by replacing CATwoCF with NODS or Segment.The leading-colour test can also be performed with 2 , using the option CFHalfCA.When running with small but finite values of 2 is that the former has a greater cancellation between the C A and n f terms in the running-coupling β-function, which ensures that the Landau pole is further in the infrared for any given value of α s .

Implementing a new shower
The PANSCALES framework allows for relatively straightforward addition of new dipole and antenna showers with alternative kinematic maps or ordering variables.This allows the user to leverage the existing code for colour and spin handling, as well as the infrastructure for logarithmic accuracy tests and the interface to PYTHIA8.3.As an example, our own toy implementation of the PYTHIA8.3 final-state shower involves about 250 lines of header (ShowerToyPythia8.hh),most of which are boilerplate code, and a further 250 lines in ShowerToyPythia8.cc,more than half of which are comments or blank lines.For a slightly more elaborate example that handles also initial-state radiation, the user may wish to look at the ShowerDipoleKt class.
Here we outline some of the aspects that a user should keep in mind in implementing their own new shower.Firstly, the code is in the panscales namespace.Inspecting the code, the user will see that rather than double, many variables are in precision_type: this corresponds to the precision that was chosen in the build step.Generically, double is used for logarithmic variables and acceptance probabilities, while precision_type is to be used for (non-logarithmic) kinematic variables and associated matrix-element calculations, where rounding errors and large exponents may be encountered.
Much of the core work of showering is carried out by the ShowerRunner class, which is essentially agnostic as to the specific details of any given shower.The basic work flow of ShowerRunner is depicted in Fig. 5.The ShowerRunner class is constructed by passing a pointer to any class that derives from ShowerBase, whose role is to handle a small, welldefined set of shower-specific tasks, such as providing the acceptance probability and implementing the shower's kinematic map.The functions that the user should provide to the shower are also shown in Fig. 5.For concreteness, with this distribution, we have supplied two files, ShowerUserDefined.hh and ShowerUserDefined.cc, derived from ShowerBase, which one can modify.The functions that need to be touched are marked with "USER-TODO" in the code.The user-defined shower can then be run with the -shower user-defined command line argument.Note that the template is meant for e + e − showers.For additional features like initial-state radiation, matching or the implementation of double-soft currents, the reader should examine the structure of other existing showers.
The class ShowerUserDefined implements: • Member functions that provide textual descriptive information about the shower.
• Member functions that provide structural information about the shower, notably whether the shower is a dipole shower (only the emitter splits), or an antenna shower (both the emitter and spectator can split) via the function only_emitter_splits(), and n_elements_per_dipoles() (two for a dipole shower, one for an antenna).For showers with a global kinematic map, one where a dipole branching also impacts the kinematics of particles not in the dipole, the function is_global() should return true (more on this later).
The two sub-classes encode the data associated with an emission and most of the implementation of the shower splitting: • ShowerUserDefined::EmissionInfo derives from ShowerBase::EmissionInfo.
A pointer to its base class is used by ShowerRunner as the main structure to keep track of the emission through the various steps of constructing that emission, and is generally  passed to all shower-specific functions.Much of the required functionality is already present in the base class.Users may extend it to cache quantities computed at specific steps of a single emission that is trialled, so as to not reproduce them again at a later stage.An example would be the absolute k t of an emission, or the longitudinal momentum fraction z, which might first be calculated in the early stages of Fig. 5 and then used later when generating the new dipole kinematics.
• ShowerUserDefined::Element derives from ShowerBase::Element. Let us focus here on the Element class.For a dipole shower, where each dipole end is associated with a distinct kinematic map, there will be two Elements per dipole, one for each end.
For an antenna shower, there will be just one Element per dipole.The branching kinematics are governed by three variables: • ln v, the logarithm of the (dimensionful) ordering variable.
• ln b, a logarithmic variable for the longitudinal degree of freedom of the emission.In the soft-collinear region of a back-to-back dipole it might, for example, map directly to a rapidity, or to ln z where z is the radiated particle's momentum fraction.
• φ, the azimuthal angle for the emission.
The main non-trivial member functions that need to be implemented are the following: • lnb_extent_const() and lnb_extent_lnv_coeff(): for any given ln v, ShowerRunner will take the ln b generation range to be taken equal to lnb_extent_const() where typically the first term would be a function of the dipole kinematics (as encoded in the corresponding Element) and the second term would not.The use of a simple analytic form for the extent facilitates the generation of the ln v distribution.
• lnb_generation_range(lnv): for a given value of the evolution scale lnv, returns a Range object containing the limits of ln b generation.The difference between the upper and lower limits should coincide with the extent expected from the previous bullet point.Not all ln b values in the range need be kinematically valid.
• lnv_lnb_max_density(): returns the maximal emission density (for a dln v dln b dφ 2π measure).Typically this would be the soft-collinear limit of the emission density where α max s is the maximal value that the strong coupling can take. 10  • lnkt_approx(lnv,lnb): returns the logarithm of the transverse momentum of an emitted parton with that ln v, ln b combination.The result should be exact in the softcollinear limit, but does not need to be exact in the soft large-angle or in the hardcollinear limits.It is typically used for the evaluation of α s (k t ), as well as for some of the dynamic emission vetoing used for logarithmic accuracy tests.
• eta_approx(lnv,lnb): similar, but returns the rapidity of the emitted parton in the soft-collinear limit, used in the computation of colour transition points and dynamic emission vetoing.
• acceptance_probability(emission_info): sets information subsequently used by ShowerRunner in order to calculate the probability that the dipole splits.It makes use of the generation variables lnv and lnb, as cached in emission_info.Specifically it sets the following member variables in emission_info: -emitter_weight_rad_gluon : weight for the emitter to radiate a gluon; -emitter_weight_rad_quark : weight for the emitter to radiate a (anti-)quark; -spectator_weight_rad_gluon : weight for the spectator to radiate a gluon; -spectator_weight_rad_quark : weight for the spectator to radiate a (anti-)quark. 10All current shower implementations have a private member pointer Element::_shower, and the base shower class has a max_alphas() member, as well as a qcd() member that supply access respectively to α max s and QCD constants.Together these provide the information needed to calculate the maximum density.
The weights to radiate a gluon/quark depend on the splitting functions.Only for an antenna shower will the spectator weights be non-zero.They are to include a factor of the maximal allowed value of the strong coupling α max s , with the ShowerRunner class then accounting for an α s (k t )/α max s factor.In initial-state branchings, ShowerRunner also accounts for an appropriate PDF ratio factor.The user may choose to use the splitting functions implemented in QCD.hh.This can be done by calling fill_dglap_splitting_weights, which needs the longitudinal radiated momentum fraction z (where z → 0 indicates the soft limit).In addition, when using spin correlations, the user should set -z_radiation_wrt_emitter : collinear momentum fraction with respect to the emitter; -z_radiation_wrt_spectator: collinear momentum fraction with respect to the spectator.
The acceptance_probability(...) function returns a bool, where false indicates that the generation variables are definitely outside the kinematic limit.For most showers, if it returns true then the generation variables would normally be inside the kinematic limit.The emitter and spectator splitting probabilities are then used later in ShowerRunner to accept/reject a splitting, and choose the splitting channel.
• do_kinematics(emission_info, rp): constructs the post-branching momenta of the emitter, the spectator and the newly radiated particle.For this, one may again use the cached generation variables lnv, lnb and phi, alongside any other quantity that the user stored in emission_info.The post-branching momenta should be stored in emission_info under the names emitter_out, spectator_out and radiation.This function returns a bool, where true indicates that the generation variables were inside the physical kinematic limit. 11Note that the pre-branching emitter and spectator momenta should be taken from a variable rp of type RotatedPieces.This class is part of the framework for handling directional differences.It provides a rotated version of the dipole with one or other of its particles aligned along the z axis, which allows the user to retain high precision in the branching kinematics (specifically, small components along the x or y axis) without explicit knowledge of the underlying direction-difference structures.The direction-difference infrastructure then takes care of deducing the correct momenta and direction differences in the original frame.
• update_event(...): In the Element base class, the update_event(...) member function takes care of replacing the pre-branching emitter and spectator particles with the post-branching ones, adds the radiated particle to the event, and takes care of some of the bookkeeping associated with colour handling.However, other particles in the event are by default not modified.Therefore, for showers with a global kinematic map, this function needs to implement additional operations on the rest of the event (e.g. a boost or rescaling).These would typically be preceded by an explicit call to ShowerBase::Element::update_event(...).The event is then further processed 11 It is possible for the acceptance_probability(...) function to return true even outside the kinematic limits.In this case there are two possible avenues for imposing the kinematic limit.One is for the shower's Element class to overload the base-class member function check_after_acceptance_probability(...), which gets called after α s and PDF factors have been incorporated into the branching probabilities and those have been used to decide to continue with the emission generation.It is called before φ and the splitting channel are known.Alternatively do_kinematics(...) is called with knowledge of the φ value and channel, and can return false if the emission is outside the kinematics.The former is the only route that is currently valid in order for spin correlations to be correctly accounted for.by ShowerRunner, updating the event dipoles, colour and spin-density structure in addition to any caching associated with the event generation.Note that cached information associated with dipoles other than the newly-created dipole and the splitting dipole are by default not updated, unless the ShowerUserDefined::is_global() function returns true.
Once implemented, the new shower can be run by using the flag -shower user-defined.
If the user would like their shower to work with direction differences, they should inspect how this is implemented in existing showers.Aside from the RotatedPieces discussion above, calculations of dot products in determining dipole invariants should make explicit use of knowledge of direction differences (available from the dirdiff_3_minus_3bar member variable of element.dipole()),and do_kinematics(...) should use 3-vectors in its internal calculations to avoid triggering off-mass-shell errors.Furthermore if the shower carries out any global boosts, these need to be performed in a way that correctly boosts also the full set of dipole direction differences.The user is invited to inspect the code of existing showers for further details.
It is important also to test the correctness of the direction-differences implementation.Typically we do this by first running a double-precision build with a physically sensible range and verifying that results are identical with/without the -use-diffs argument.Then we create a build in the doubleexp type, running with -use-diffs and a logarithmic range of about 1000 (and correspondingly low α s , so that the multiplicity stays of the order of 10−100) and compare the output to a build with the much slower mpfr4096 type, running without -use-diffs.Again the results should be identical, though typically a few iterations are likely to be necessary to identify all sources of potential loss of precision.A final comment is that, by default, the logarithmic-accuracy tests of section 2.3 run in double precision, with a ln v range reaching about 300.However, for this to work, the shower should not ever do more than take the square of a momentum, otherwise the result will exceed the exponent that can be represented in double precision.If this is a problem, the user should modify a configuration flag in example-global-nll-ee.pyso as to use doubleexp (which is somewhat slower).

Conclusions
This 0.1 series release of the PANSCALES code allows users to start exploring its features and techniques, notably for tests of logarithmic accuracy of parton showers.It also demonstrates an early version of the interface to PYTHIA8.3.While we do not yet recommend its use for phenomenological applications, we hope that this early release of the code will provide a foundation for exploring connections with other projects, so as to enable the wide ecosystem of collider physics tools to benefit from the validated logarithmic accuracy of the PANSCALES showers.LS, and RP\R\231001, GPS) and by the Science and Technology Facilities Council (STFC) under grants ST/T000864/1 (MvB, GPS), ST/X000761/1 (GPS), ST/T000856/1 (KH) and ST/X000516/1 (KH), ST/T001038/1 (MD) and ST/00077X/1 (MD).LS is supported by the Australian Research Council through a Discovery Early Career Researcher Award (project number DE230100867).The work of PM is funded by the European Union (ERC, grant agreement No. 101044599).Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency.Neither the European Union nor the granting authority can be held responsible for them.We also thank each others' institutes for hospitality during the course of this work.

A Event record A.1 Basic event information
Here we give an example event record, printed by default for the first event of any run (the first N events can be printed with the command-line option -max-print N).The printout includes the hard event before showering, and the final fully showered event.
The ID corresponds to PDG code for the particle, a non-zero value of IS indicates that the particle is in the initial state and HS is non-zero for particles that are considered to be part of the hard system (here, particle 2, which is the Z boson).The DIS variable is relevant only to DIS and VBF processes and in particular for the latter it indicates which of the two sides of the VBF process a given particle is associated with.From the point of view of an analysis, all particles with IS equal to zero are valid final-state particles.Within an analysis, the particles are available through f_event.particles(),which returns a vector of Particle objects.That same f_event object also can be used to access the event dipoles.

A.2 Access to step-by-step event information
As discussed in Section 2.2, by default PANSCALES does not store the event structure separately after each branching.However it is possible to ask for the corresponding information to be stored.This is achieved by calling f_shower_runner.set_step_by_step_event_caching(true) in the user_startup phase of an analysis.Once this is done, after each event the user can access f_shower_runner->cached_step_by_step_event() which has various enquiry functions, for example for retrieving the branching variables or the event at any given step of the shower.Note that this interface may still evolve.For debugging purposes, there is also a compile-time option (-DPSVERBOSE=on) that can be passed to cmake.When enabled, running a shower provides very extensive output as to the shower's individual steps.

Figure 1 :
Figure 1: Example results from tests of NLL accuracy of the PanGlobal (β PS = 0) shower for three global event-shape observables, the Cambridge algorithm y 23 and the Lund β obs = 0 sum and maximum observables, corresponding respectively to the three panels.Each panel shows the ratio of the shower cumulative-distribution result to the NLL calculation.Each coloured line corresponds to a specific value of α s , while the black line gives the extrapolation to α s = 0.The results are shown as a function of the maximum allowed value of λ = α s ln O, where O is the observable.

5 obs 5 obsFigure 2 :
Figure 2: Dependence of Σ(λ) on the size, δ ln O , of the parton-shower dynamic cutoff.The results are shown for the S j,β obs class of observables in Drell-Yan production, which sum p t i e −β obs |η i | over all final-state jets i (with p t defined as the transverse momentum with respect to the beam).The dependence on the size of the dynamic cutoff size is shown normalised to the NLL prediction for Σ, for each of three values of β obs : 0 (red), 0.5 (green) and 1.0 (blue).Two showers are shown (a) PanGlobal with β PS = 0 and (b) PanLocal (dipole) with β PS = 0.5.The upper panels show the ratio of Σ to the NLL result, showing the convergence as δ ln O is made more negative (note that the ratio is not expected to go to 1, because of the finite value of α s ).The lower panels show the difference between Σ/Σ NLL results for successive δ ln O values, again normalised to the NLL result, giving a clearer view of the degree of convergence.The dashed lines help illustrate that the behaviour is consistent with exponential dependence on δ ln O .

Figure 3 :
Figure 3: Validation of the weighted event generation.We show the results for Dipole-k t for pp → H collisions with s/m H = 100 and y H = 0. We take α s = 0.1 and target λ = −0.5.(a) The distribution of ln v 1 values for the first emission, with and without weighted generation, illustrating that they agree.The bottom panel shows the ratio between the weighted and unweighted results minus 1, where the statistical uncertainty is indicated with the red band.(b) The relative statistical uncertainty in the ln v 1 distribution for both unweighted and weighted shower generation.
̸ = β obs S t a r t s h o w e r a t l n ( v a ) Sudakov ∆(v max , v a ) C e n h = 1 b e l o w l n ( v b ) , w → w C enh −1

Figure 5 :
Figure 5: Simplified structure of a single trial for an emission in the shower.Steps on the left are performed by the centralised ShowerRunner class.Shower-specific steps, shown on the right, are the responsibility of the ShowerUserDefined implementation, as well as its sub-classes.The EmissionInfo sub-class stores all required information needed for the current branching, such as the splitting variables and constructed kinematics.The emitting dipole is accessible in the Element sub-class, which contains the majority of the shower-specific code.The matching, implementation of the spin correlations throughout the event evolution, and the double-soft splitting functionalities are not illustrated here.
Choose lnv and elementReturn lnb rangeChoose lnbReturn lnkt for α s It is responsible for carrying out almost all of the shower-specific work.In case of multiple types of elements (i.e.II, IF and FF dipoles), one might choose to derive distinct element classes for each, typically each from an intermediate ShowerUserDefined::Element base class (see e.g.ShowerDipoleKt).