Production of tau lepton pairs with high pT jets at the LHC and the TauSpinner reweighting algorithm

The TauSpinner algorithm allows to modify the physics of the Monte Carlo generated samples due to the changed assumptions of event production dynamics, without re-generating events. To each event it attributes weights: the spin effects of tau-lepton production or decay, or the production mechanism are modified. There is no need to repeat the detector response simulation. We document the extension to 2 to 4 processes in which the matrix elements for the parton-parton scattering amplitudes into a tau-lepton pair and two outgoing partons are used. Tree-level matrix elements for the Standard Model processes, including the Higgs boson production are used. Automatically generated codes by MadGraph5 have been adapted. Tests of the matrix elements, reweighting algorithm and numerical results are presented. For averaged tau lepton polarisation, we perform comparison of 2 to 2 and 2 to 4 matrix elements used to calculate the spin weight in pp to tau tau j j events. We show, that for events with tau-lepton pair close to the Z-boson peak, the tau-lepton polarisation calculated using 2 to 4 matrix elements is very close to the one calculated using 2 to 2 Born process only. For the m_(tautau) masses above the Z-boson peak, the effect from including 2 to 4 matrix elements is also marginal, however when restricting into subprocesses qq,q bar q to tau tau j j only, it can lead to a 10% difference on the predicted tau-lepton polarisation. Choice of electroweak scheme can have significant impact. The modification of the electroweak or strong interaction can be performed with the re-weighting technique. TauSpinner v.2.0.0, allows to introduce non-standard couplings for the Higgs boson and study their effects in the vector-boson-fusion. The discussion is relegated to forthcoming publications.


Introduction
With the data collected so far by LHC experiments, there was not much interest to explore physics of τ-lepton decays, with the exception of exploiting τ leptons in searches for rare or Standard-Model-forbidden decay channels, see eg. [1]. However, τ-lepton signatures can provide a powerful tools in many areas, like studies of hard processes characteristics, measurements of properties of Higgs boson(s) [2,3], or in searches for New Physics [4,5,6].
The τ leptons cannot be observed directly due to their short life-time. All decay products are observed, with the exception of ν's. There are more than 20 different τ decay channels, each of them leading to a somewhat distinct signature. This makes a preparation of observables involving τ decays laborious. However, such efforts can be rewarding, because τ-lepton spin polarization can be measured directly, contrary to the case of electron or muon signatures, giving better insight into the nature of its production mechanism, e.g. the properties of resonances decaying to τ leptons. This is the main motivation for developing TauSpinner, an algorithm to simplify the task of exploring the τ physics potential, which could be used for evaluation/modification of event samples including τ decays.
In the first release, the program algorithms were focused on longitudinal spin effects only [7]. Already TauSpinner ver.1.1 handled these effects with the help of the appropriate spin weight attributed to each event. In this way, spin effects could be introduced, or removed, from the sample. With time, variety of extensions were introduced. Since Ref. [8], a second weight was introduced which allows to manipulate the production process by adding additional contributions or completely replacing the production process with an alternative one, including for example an exchange of a new intermediate particle.
Ref. [9] brought a possibility of modifying transverse spin effects in the cascade τ decays of intermediate Higgs boson. Later, Ref. [10] enabled the transverse spin effects for the case of τ leptons produced in Drell-Yan processes to be studied as well.
With time, technical options or important precision improvements were introduced too. In [9], an option to attribute helicity states to τ-leptons was introduced. One should keep in mind, that because of quantum entanglement, the assignment of a definite helicity state to intermediate τ's is necessarily subject to an approximation. However, for spin weight calculation, the complete spin density matrix is taken into account and in general, approximation is not used. With later publication [10], one-loop electroweak (EW) corrections also became available for the Drell-Yan parton-process qq → Z/γ * → ττ.
Let us mention another technical option. Initially the program was expected to work for samples, where spin effects are either taken into account in full, or are absent. One can however configure TauSpinner algorithm to work on generated samples where only part of spin effects is taken into account (only some components of the density matrix used) and to correct them to full spin effects.
Until now, for calculations of spin weights, TauSpinner algorithm was always using the Born-level (2 → 2) scattering amplitudes convoluted with the corresponding parton distribution functions (PDFs). Kinematic configurations of the incoming/outgoing partons were reconstructed from the four-momenta of outgoing τ leptons and incoming protons (using c.m. collision energy), and somewhat elaborated kinematical transformations were used for calculating an effective scattering angle of the assumed Born process.
The validity and precision of this approximation became of a concern, especially for configurations with high momentum transfers in the t-channel and for outgoing particles with high transverse momentum (p T ) that accompany decay products of the electroweak bosons. In such cases, more elaborated description of the production process dynamics is needed. The aim of the present paper is to describe an improved version of TauSpinner 2.0.0 which now includes hard processes featuring tree-level parton matrix elements for production of a τ-lepton pair and two jets. Numerical test, and some results of physics interest, will be also presented.
The paper is organized as follows: In Section 2 we recall assumptions used for the Monte Carlo reweighting techniques, in particular for the modeling of kinematic distributions in the multi-dimensional phase-space. We then define the master formula used by TauSpinner for modeling spin correlations of τ-lepton decay products in events with different topologies in proton-proton collisions. Section 3 documents details of the tree-level matrix elements used for the calculation of weights in pp → ττ j j events. The implemented functionality is based on automatically produced FORTRAN code from MadGraph5 package [11] for processes of the Drell-Yan-type and of the Standard Model Higgs boson production in vector boson fusion (VBF) processes, which have been later manually modified and adapted. Numerical effects of different choices for electroweak and QCD interactions initialization are presented in the last two subsections. We classify parton level processes into groups, which are then used in the following Section 4 for technical tests. We explain details of the modification which we have introduced to the initialization of MadGraph5 generated amplitudes and emphasize the necessity of using the effective sin 2 θ e f f W for the calculation of the coupling constants to correctly model the measured spin asymmetries in the Drell-Yan process. This is even more important for a correct generation of angular distributions of leptons in the decay frame of intermediate Z bosons. Then we discuss combinatorial and CP symmetries that allow us to reduce the number of parton subprocesses for which distinct codes of spin amplitudes are needed. (Appendix A is devoted to describe technical details of the introduced extension of TauSpinner.) Numerical results shown in Section 5 are divided into three parts. The first one is devoted to the evaluation of systematic biases present if the (2 → 2) variant of TauSpinner is used for spin effects (or for the matrix element weights) in pp → τ τ j j processes. Next, we present numerical consequences of the choice of the electroweak scheme, in particular: (i) in the τ τ j j production, (ii) in the calculation of the spin correlation matrix used for the generation of τ decays, for the observable distributions. Section 6, closes the paper. Somewhat lengthy collection of tests are relegated to Appendices B and C.
In the present paper we concentrate on physics oriented aspects of new implementations. All technical details and a description of available options, resulting not only from the present work but also from the previous publications on TauSpinner, will be collected in a forthcoming publication. The most important points for technical aspects of the program use are nonetheless presented in Appendix A. Benchmark outputs from the programs are relegated to the project web page [12].

Theoretical basis
Before we start the discussion of new implementations in the TauSpinner and present numerical results, let us shortly recall the basis of the approach being used. For the Monte Carlo techniques of calculating integrals or simulating series of events, to be well established in the mathematical formalism, one has to define the phase-space and the function one is going to integrate. One can parametrize the integral in the following form where on the right-hand side, the sum runs over n dimensional vectorsx i j of random numbers (eachx i j in the [0, 1] range) which define the point in the hypercube of coordinates. The N denotes number of events used.
The function g consists of several components: the phase-space Jacobian resulting from the use ofx j coordinates for the phase-space parametrization; the matrix element squared calculated for a given process at prepared phase-space-point; and finally the acceptance function which is zero outside the desired integration region. Uniformly distributed random numbersx j are used as Monte Carlo integration variables in formula (1). The average value of g, calculated over the event sample, gives the value of integral G. In practical applications a lot of refinements are necessary to assure acceptable speed of calculation and numerical stability. From a single sample of events several observables can be obtained simultaneously, e.g. in the form of differential distributions (histograms).
For the convenience of calculating multi-dimensional observables one introduces rejection techniques. The event i (constructed from random-number variablesx i j ) is accepted if an additional randomly generated number is smaller than g(x i 1 ,x i 2 , ...,x i n )/g max ; otherwise the event is rejected. The result of the integral is then equal to g max × n accepted n generated . Statistical error of this estimate can be calculated using standard textbook Monte Carlo methods. In such a method, one has to assure that for the allowedx j range, the condition 0 ≤ g ≤ g max holds. The accepted events are distributed according to dG and can be used as a starting sample for the next step of the generation of weighted (or weight 1) events.
The principle goal of the TauSpinner program is to un-do, modify or supersede the discussed above rejection. Let us assume that the sample of events, for which the program will be used, are distributed accordingly, with all details, to the known production mechanism described by the formula In Eq. (2), the ∑ i, j,k,l extends over all possible configurations of incoming and outgoing partons for the processes of i(p 1 ) j(p 2 ) → k(p 3 ) l(p 4 ) τ + τ − . The p i stand for the 4-momenta of incoming/outgoing partons, x 1 and x 2 stand for energy fractions of the beams carried by the incoming partons, parton distribution functions are denoted as f i (x 1 ), f j (x 2 ) respectively for the first and the second incoming proton. The parton-level flux factor is denoted as Φ f lux and the phase-space volume element as dΩ.
Finally the parton-level matrix element M i, j,k,l completes the formula. Obviously parton distributions (PDFs) are dependent on parton flavour configurations. In Eq. (2) the τ decay phase-space and the corresponding matrix elements are omitted. Even though it amounts to semi-factorization, exploited by TauSpinner algorithms, we omit for now also the discussion of τ-spin correlation matrix. They are not essential for the clarification of requirements needed for TauSpinner algorithms. For the calculation of TauSpinner weights in the case of replacing one production mechanism A with another one B, one has to take into account not only differences in the matrix elements and PDFs but also, potentially, in Φ f lux and dΩ. Thus the respective weight 1 is calculated as follows Although the factors Φ f lux and dΩ may cancel between the numerator and denominator in the case when all incoming and outgoing partons are considered to be massless, they still may differ due to symmetry factors which are different for identical or distinct flavours of partons.
The physics processes of interest are the Standard Model processes in pp collision with two opposite-sign τ leptons and 2 jets (quarks or gluons) in the final state 2 . Such processes are described at the tree level by (2 → 4) matrix elements, with intermediate states being single or double Z,W, γ * , H or fermion exchange in the s-or t-channel. Depending on the initial state, tree-level matrix elements are of the order of α S α EW or α 2 EW , involving sometimes triple WW Z couplings. More details are given in Table 1. We will limit our implementation to the tree-level only, but with the emphasis on controlling the spin configurations.

Incorporating MadGraph generated code into TauSpinner
There are automated programs for generating codes of spin amplitudes calculation. In the development of TauSpinner we have used MadGraph5 [13]. Let us recall some details of this step of the program development to explain the adopted procedure, which may be useful in future for introducing anomalous couplings or new physics models.
The FORTRAN code for calculating matrix elements squared (ME 2 ) is generated using MadGraph5 with the following commands: a) import model sm-ckm b) with default definition of "multiparticles" p = g u c d s u˜c˜d˜sj = g u c d s u˜c˜d˜sc ) for the Higgs signal processes generate p p > j j h, h > ta+ tad) for the Drell-Yan-type SM background processes generate p p > j j ta+ ta-/ h QED=4 e) and print the output using output standalone "directory name".
Setting the parameter QED=4 enforces generation of diagrams up to 4th order in the electroweak couplings. Other settings are initialized as in the default version of the MadGraph5 setup. The generated codes for the individual subprocesses are then grouped together into subroutines, depending on the flavour of initial state partons, and named accordingly. For example, SUBROUTINE UDX (P,I3,I4,H1,H2,KEY,ANS) corresponds to processes initiated by ud partons. X after the letter U,D,S and C means the antiquark, i.e. UXCX corresponds to processes initiated byūc, while GUX -processes initiated by gū. The input variables are: real matrix P(0:3,6) for fourmomenta of incoming and outgoing particles, integers I3,I4 for the Particle Data Group (PDG) identifiers for final parton flavours, integers H1,H2 stand for outgoing τ helicity states; integer KEY selects the requested matrix element for the SM background (KEY=0), the SM Higgs boson 3 (KEY=1), ANS returns the calculated value of the matrix-element squared. According to the value of I3,I4,KEY the corresponding subroutine generated by MadGraph5 is called 4 . The TauSpinner user usually will not access to the KEY variable. Before integrating these subroutines into the TauSpinner program, a number of modifications have been done for the following reasons: a) Since MadGraph5 by default sums and averages over spins of incoming and outgoing particles, while we are interested in τ spin states, the generated codes have to be modified to keep track of the τ polarization; b) Moreover, since the subroutines and internal functions generated by MadGraph5 have the same names for all subprocesses SMATRIX(P,ANS), the names had to be changed to be unique for each subprocess. To be more specific, for the Higgs signal subprocess ud → cd h, h → τ + τ − the generated subroutine name is changed to UDX_CDX_H(P,H1,H2,ANS), while for the background ud → cd τ + τ − process, the generated subroutine name is changed to UDX_CDX_noH(P,H1,H2,ANS). Table 1: List of implemented processes for calculating matrix element squared grouped into categories, which differ by flavours of incoming partons. For each category, FORTRAN files with implemented subroutines for calculating the matrix element square, grouped by the flavour of incoming partons, are given in the second column. Examples of processes in each category are given in the last column. Partially redundant codes for some of the processes are used for tests only, this is the case of amplitudes stored in files UCX.f and CUX.f, the amplitudes of these two files can be obtained from each other by CP symmetry.

Category of Corresponding FORTRAN files Processes Matrix Elements
For other processes and internal functions similar convention is used, see Table 1. Note that for example for the processes with cs quarks in the initial state, exchange of W ′ s is allowed, the final states cannot include gluons and the only allowed final states are: cs, cd, us, ud. After taking into account permutation of incoming and outgoing partons and CP symmetric states this gives in total 4 × 4 × 2 = 32 non-zero contributions to the sum of Eq. (2). This is the case both for Drell-Yan-type background and Higgs-boson production processes. For the remaining processes the codes listed in Table 1 are also used with the help of CP symmetry or re-ordering of partons.
At the parton level each of the incoming or outgoing parton can be one of flavours:bcsūd g d u s c b, with Particle Data Group (PDG) identifiers: -5, -4, -3, -2, -1, 21, 1, 2, 3, 4, 5 respectively. For processes with two incoming partons, two outgoing τ leptons and two outgoing patrons that gives 11 4 possibilities, most of them with the zero contribution, and many available one from another by relations following from CP symmetries and/or permutations of incoming and/or outgoing partons. Grouped by the type of initial state partons, the subroutines listed in Table 1 are currently limited to the first two flavour families. The matrix elements for processes involving b-quarks are not yet implemented 5 .
Also, for practical purposes, for a pair of final-state parton flavours k = l, the MadGraph5 generated codes have been obtained for a definite ordering (k, l), but not for (l, k), to reduce the number of generated configurations. When TauSpinner is invoked, the flavour configuration of outgoing partons is unknown and it takes into account both possibilities: thus a compensating factor 1+δ i j 2 has to be introduced. This is because of the organization of the sum in Eq. (3).

Topologies and the dynamical structure of subprocesses
The number of contributing subprocesses is very large. For the case of the non-Higgs Drell-Yan-type background processes, in which the τ-pair originates either from the vector boson decay (including also cascade decays) or from multi-peripheral vector-boson fusion processes, MadGraph5 generates 82 subprocesses with partons belonging to the first two generations of quarks, or gluons. Subprocesses in which all partons are of the same flavour (like uū → uūτ + τ − ) receive contributions from 64 Feynman diagrams, subprocesses with two pairs of flavours -either 43 diagrams (if one pair is of up-type and the other down-type, like uū → ssτ + τ − ) or 32 diagrams (if both pairs are either down-or up-type, like uū → ccτ + τ − ), subprocesses with three or four different flavours -11 diagrams (like us → udτ + τ − ), and subprocesses with two quarks and two gluons -16 diagrams. As far as the dynamical structure of the amplitudes is concerned, there are all together seven different topologies of Feynman diagrams, with representatives shown in Fig. 1. Which of them contribute to a given subprocess depends on flavours of incoming and outgoing partons. Irrespectively of their origin, in all processes the polarizations of τ leptons are strongly correlated due to the helicity-conserving couplings to the vector bosons. The spin correlations of the produced τ pair depend on the relative size of the subprocesses with vector and pseudo-vector couplings contributing to the given final state configuration. For example, in the case of qq → τ + τ − qq , see Fig. 1, diagram (d) contributes with 100% polarised τ's since  they couple directly to W ± . In diagram (g), the polarisation of Z/γ * is different than in the Born-like production because Z/γ * decaying to τ + τ − originates from the WW Z/γ * vertex. This leads to a distinct polarisation of τ leptons. For the Higgs signal processes the τ pairs originate from the Higgs boson decay, as imposed at the generation level, and the number of subprocesses is reduced to 67. Each subprocess receives contributions from at most two Feynman diagrams, since with massless quarks of the first two generations, the Higgs boson can originate either from the vector boson fusion or from Higgs-strahlung diagrams, as illustrated in Fig. 2. Depending on the flavour configuration of incoming partons, mediating boson is W or Z, which leads to almost 10 GeV shift between resonance invariant mass of the outgoing pair of jets in case of Higgs-strahlung process. The helicity-flipping scalar coupling to the Higgs boson results in the opposite spin correlation as compared to the case of the Drell-Yan-process. The individual τ polarization is absent.
Concerning the analytic structure of the differential cross sections, it is determined by topologies of contributing diagrams to a particular subprocess. For example, s-channel propagators will result in a resonance enhancement, while the t-channel ones may lead to collinear or soft singularities (in the limit m 2 W /s ≪ 1, m 2 Z /s ≪ 1) regulated either by the phase space cuts or by the virtuality of the attached boson line. Understanding differences in analytic structures of subprocesses will turn important when discussing tests of reweighing technique of TauSpinner in Subsection 4.2.
Technically speaking, the sums in Eq. (1) or (3) defining the production weights used in TauSpinner consist of 9 4 (11 4 if b-quarks are allowed) elements, which are potentially distinct and require their own subroutines for the matrix element calculation. Since most of the elements are equal zero, or some matrix elements are related to others by permutation of partons and/or CP symmetries, special interfacing procedure is prepared to exploit those relations. It reduces significantly the computation time and size of the program code. Details are given in Appendix A.

EW scheme and parameters
In early versions of TauSpinner the electroweak interactions were embedded into an effective (2 → 2) Born process for qq → τ + τ − . Its analytic form is given by Eqs. (3)-(5) and Table 2 of Ref. [14]. The adopted scheme is fully compatible with the one of Tauola universal interface [15]. It is using the lowest order ME for the qq → Z/γ * → ττ process, however with the effective value for the sin 2 θ e f f W and running Z-boson width. Such a choice corresponds to a partial resummation of higher order electroweak effects, exactly as it was adopted at the time of precision tests of the Standard Model at LEP [16], with the remaining loop weak corrections at the per mille level.
However, since the effects of WW boxes can be numerically significant for τ lepton pairs of large virtuality or large invariant mass, there is an option to include genuine weak loop effects into TauSpinner effective Born, already since its version 1.4.0 of June 2014 as well. It can be done for TauSpinner in a manner similar to Tauola universal interface [15] because this process is implemented in both codes in the same way. Table 2 compares numerical values of the input parameters for the (2 → 2) and (2 → 4) processes. Variants of initialization for (2 → 4) processes are explained in Table 3. It is worth to point out that by using over-constrained set of parameters essential effects of the loop corrections are taken into account, providing the results for τ-lepton polarisation close to LEP measurements [16]. As the parameters are not independent, this can lead to problems if the input values are not consistent, especially when applied to processes other than (2 → 2), which is the main focus of our paper.
Let us now turn to the details of electroweak schemes used for the matrix elements of the (2 → 4) hard subprocesses entering the pp → ττ j j. The code from MadGraph5 has its own initialisation module consistent with the so called G F scheme, which uses G F , α QED and m Z as input parameters, see Table 2 (and EWSH=1 scheme in Table 3). As such, it uses tree-level (equivalent to on-shell) definition of the weak mixing angle sin 2 222246. This is far from the measured value from the Z boson couplings to fermions; also the constant width of Z-boson is used. Since the τ-lepton polarization is very sensitive to the value of the mixing angle, and for both Tauola and TauSpinner the τ physics is important target, such a LO implementation in the G F scheme is not sufficiently realistic. This is even more serious issue for the angular distributions of leptons themselves, making such a scheme phenomenologically inadequate to any observable that relies on directions of leptons. Alternatively, one could adopt the scheme with G F , m Z and sin 2 θ e f f W as input parameters (and EWSH=2 scheme in Table 3), but then the predicted tree-level W -boson mass is away from the measured value which would result in distorted spectra of jets coming from W decays (and shift in the resonance structure of the matrix element). In some regions of the phase space the distortion can reach 40%. One can also use scheme with G F , m Z and m W , as input parameters (and EWSH=3 scheme in Table 3), but the on-shell definition sin 2 θ W = 1 − M 2 W /M 2 Z = 0.222246 will lead back to far from measured value of sin 2 θ W . There are two options: either include EW loop corrections simultaneously with QCD corrections, or adopt an effective scheme which would allow at tree-level to account correctly for the τ-lepton polarization at the Z-boson peak and physical W -boson mass. Since the former is beyond the scope of the present paper, we take the second option.
To this end, we define an effective scheme with θ  Table 3). Although being in principle flavour dependent, the value sin 2 θ e f f W is flavour universal with an accuracy of order 0.1%. Effectively such a procedure amounts to the inclusion of some of higher order EW corrections to the Zτ + τ − vertex. 6 This value is used in all vertices, also in the triple gauge-boson coupling since the WW Z coupling is essential for the gauge cancellation and it must match the couplings in other Feynman diagrams, forming together the gauge invariant part of the whole amplitude. In our case we are not aiming at a careful theoretical study of higher order corrections; instead we checked numerically that the introduction of dominant loop corrections to Zτ + τ − vertex through the effective sin 2 θ e f f W does not lead to numerically important consequences for the WW Z vertex. For example, the effect of the mismatch of WW Z and Z ff couplings for the case of qq → qqττ subprocess is small, see Fig. 7, in Section 4.4. Thus, we gain consistency with observables, such as τ-polarization or τ-directions, which would otherwise be off by ∼ 40% at the expense of breaking EW relations in higher order of perturbation theory. Moreover, since TauSpinner is used to reweigh events, as given in Eq. (3), the uncertainties of our procedure should to a large degree cancel out.
For the purpose of comparison of the predicted τ-lepton polarisation at the Z-boson peak, we provide four initialisation options for the (2 → 4) matrix elements, the first three motivated by the schemes used in [18], and the fourth one corresponding to θ e f f W . They are specified in Table 3. Scheme labeled EWSH = 4 is the only one numerically appropriate for use when predicting τ polarisation, and taking into account configurations with two additional jets, as shown in Section 5. For technical testing purposes we also introduce a scheme like EWSH=4 but with modified WW Z coupling by 5% which we label as EWSH=5. Table 2: Input parameters for initialising couplings calculations for (2 → 2) in Tauola code and (2 → 4) in MadGraph5 code. Note that in Tauola code α QED = α QED (Q 2 = 0) is used as an input for calculation of the Z couplings as well. This leads, in principle, to an over-all missing factor of ( α QED (0) ) 2 . It can be thus dropped off, as long as it cancels out in calculation of weights, the ratios of differential cross-sections. The numerical values of CKM matrix are taken from Ref. [11].

Type
Tauola code Input/Calculated Table 3: Implemented EW schemes, the recommended EW scheme is EWSH=4 which gives the τ lepton polarisation on the Z-boson mass peak, in agreement with the measurement at LEP1 [19], and physical W boson mass.

QCD scales and parton density functions
The distribution version of TauSpinner is interfaced with LHAPDF v6 library [20]. User has the freedom of choosing renormalization and factorization scales, within the constraint that µ F = µ R , otherwise minor re-coding is necessary. To this end we have implemented four predefined choices for the scale µ 2 as should be expected for our processes: where sums are taken over final state particles of hard scattering process. For the α s (µ 2 ) we provide, as a default, a simple choice of the µ 2 dependence, following the leading logarithmic formula, with the starting point α s (M 2 Z ) = 0.118. The same value of α s is used for the case of the fixed coupling constant, that is for scalePDFOpt=0.
The reweighting procedure of TauSpinner itself may be used to study numerically the effects of different scale choices, as well as for the electroweak schemes, see the discussion later in Section 5 and Appendix A.2.

Tests of matrix elements using fixed kinematical configurations
For the purpose of testing the consistency of implemented codes, generated with MadGraph5 and modified as explained in Sect. 3.1, we have chosen a fixed kinematic configuration at the parton level 7 . For such kinematics we have calculated the matrix element squared for all possible helicity configurations of all subprocesses using the codes implemented in TauSpinner and checked against the numerical values obtained directly from MadGraph5. The agreement of at least 6 significant digits has been confirmed.

Tests of matrix elements using series of generated events
As further tests of the internal consistency of matrix element implementation in TauSpinner we have used the reweighting procedure by comparing a number of kinematic distributions obtained in two different ways: the first one obtained directly from events generated for a specified parton level process REF (a reference distribution REF), and the second one (GEN reweighted) obtained by reweighting with TauSpinner events generated for a different process GEN. These tests have been performed in a few steps as follows.
• Series of 10 million events each for a number of different processes in pp → ττ j j (with specified flavours of final state jets, or for subprocesses with selected flavours of incoming partons) with MadGraph5_aMC@NLO [11] v2.3.3 at LO have been generated. Samples were generated for pp collisions at the c.m. energy of 13 TeV using CTEQ6L1 PDFs [21] linked through LHAPDF v6 interface. Renormalization and factorization scales were fixed to µ R = µ F = m Z . Only very loose selection criteria at the generation level were applied: invariant mass of the ττ pair was required to be in the range 8 m ττ = 60 − 130 GeV, and jets to be separated by ∆R j j > 0.1 and with transverse momenta p j T > 1 GeV. A complete configuration file used for events generation is given in the file MadgraphCards.txt which is included for reference in TAUOLA/TauSpinner/examples/example-VBF/benchfiles directory.
• The testing program was reading generated events stored in the LesHouches Event File format [22] filtering the ones of a given ID1, ID2, ID3, ID4 configuration of flavour of incoming/outgoing partons corresponding to the process GEN. The weight wt ME allowing to transform this subset of events into the equivalent of reference REF one, was calculated as and kinematic distributions of reweighted events (GEN reweighted) were compared to distributions of the reference process REF (ID1, ID2, ID3', ID4'). Note that for this test to be meaningful one has to select processes with the same initial state partons, so that the dependence on the structure functions cancels out. A very good agreement between the REF and GEN reweighted distributions was found for 10 different kinematic distributions for several configurations of (ID1, ID2, ID3, ID4, ID3', ID4'). It has shown a very good numerical stability, which was not obvious from the beginning as events corresponding to the REF and GEN processes may have very different kinematic distributions due to their specific topologies and resonance structures of Feynman diagrams.
• In the next step, the tests were repeated, but now reweighting the matrix elements convoluted with the structure functions of the incoming patrons and summing over final states restricted to the selected sub-groups (named respectively C and D of parton level processes). In this case the weight is calculated as where the notation as for Eq. (3) is used, except that now the ∑ C,D mean that summation is restricted to processes belonging to sub-groups C, D, respectively. For testing the code implementation for the Drell-Yan process the groups, listed in the first column of Table 1, were reweighted, one to another.
The reweighting tests performed between sub-groups of processes, and later, between groups of processes listed in Table 1, allowed to check relative normalization of amplitudes. Again, a good agreement has been found 9 . For the tests, the following kinematical distributions were used: − Pseudorapidity of an outgoing parton j. − Pseudorapidity gap of outgoing partons. − Rapidity of the ττ and j j systems. − Transverse momentum of the ττ and j j system. − Invariant mass of the ττ and j j system. − Longitudinal momentum of the ττ and ττ j j. − Cosine of the azimuthal angle of τ lepton in the ττ rest frame.
Let us discuss some of these results, shown in Fig. 3 and Fig. 4 (the complete set of distributions is shown in Appendices B and C). In each plot the distribution REF for the reference process is shown as a black histogram, while the red histogram shows the distribution for a different process GEN. Both histograms are obtained directly from the MadGraph5 generated samples of REF and GEN processes, respectively. Now the histogram GEN is reweighted using TauSpinner and the resulting reweighted histogram is represented by the red points with error bars. For the test to be successful the red points should follow the black histogram; the ratio of the REF and GEN reweighted distributions is shown in the bottom panel of each figure.
Let us note, that in our tests, we reweight events of substantially different dynamical structures over the multi-dimensional phase-space. This may be not evident from the histograms shown in figures, which can be both for the REF and GEN reweighted distributions rather regular and similar. Nevertheless, several bins of GEN reweighted distributions with small errors can be found to lie below the REF distribution, whereas a few above with large errors. This second category of bins is populated by a few events, which originate from the flat distribution of the GEN process, receiving high weight due to some resonance/collinear configuration of the REF process. This is a technical difficulty for the testing, but is not an issue of the actual use of TauSpinner when all subprocesses are used together. To confirm that the observed deviations are not significant statistically we have reproduced plots from Fig. 3 and Fig. 4 for four independent series of events. We observed that bins with large error or sequences of few bins with large deviations were randomly distributed between these series strongly indicating that observed deviations are of statistical origin. As primarily we are not interested in use of implemented code to reweight between the groups of parton level processes, for checking general correctness of its implementation it was sufficient to use four statistically independent samples only. In practical applications, contributions from all processes will be merged together and weights will become less dispersed. Similar tests have been performed for the Higgs boson production. Fig. 5 shows the comparison of generated and reweighed distributions for the jet pseudorapidity and for the pseudorapidity gap between jets in the case of qq and qq processes. As the resonant structure in the m j j distribution coming from Z → qq and W → qq is different in REF and GEN processes, results of some bins feature unexpectedly large statistical fluctuations.
Finally, let us stress that simple, but nonetheless, necessary check have been done as well: from the inspection of the control outputs we confirmed that the dominant contributions to cross sections are distinct for Drell-Yan and Higgs production processes, and that the slopes of energy spectra of τ-decay products are of a proper sign. That confirms that our installation is free of possible trivial errors in spin implementation. Drell-Yan processes j j τ τ → p p Figure 6: Impact on the matrix element calculation of parton shower smearing, as explained in the text. On the left, the difference of spin weights calculated with and without ISR parton shower kinematic smearing is shown. On the right, the ratio of matrix element weights calculated for the two cases is shown. Sample of 10000 events was used.

Tests with hard process + parton shower events
After technical tests at the hard process level (convoluted with structure functions), we turn to check the algorithm on events where the incoming parton momenta can not be assumed to be along the beam direction due to the presence of the parton shower in the initial state (ISR). For that purpose, we have taken events generated with MadGraph5 and added ISR with the default version of Pythia 8.2 (as described in Ref. [23]). The two statistically correlated samples were constructed and used by TauSpinner for calculation of spin weights (wt spin ) and production weights (wt prod ). Fig. 6 shows the number of events as a function of differences for the spin weights calculated for each event from configurations with and without ISR parton shower. Similarly, shown is the ratio of wt prod weights calculated for configurations with and without ISR parton shower. One can see from Fig. 6 (left plot), that the spin weights for the cases with and without ISR are strongly correlated. Majority of events reside in central bins of the distribution and the difference in weights is smaller than the bin width. Also the matrix element weights for the two cases are strongly correlated, see Fig. 6 (right plot). Majority of events reside in central bins. We can conclude that, similarly as in the past [7] for the (2 → 2) process, the algorithm which is applied to kinematics of the hard process particles effectively removes impact of the initial state transverse momentum and leads to results which are stable with respect to the presence of extra showering. This test is of more physical nature, since in such a case Eq. (2) does not hold for the distribution of reweighted events and, as a consequence, reweighting with Eq. (3) is featuring an approximation, which we have validated with this test. Note that adding ISR means only that the system of partons and τ leptons outgoing from the hard process underwent (as a whole) a boost and rotation before calculating matrix elements and PDF's. This justifies the evaluation of x 1 , x 2 , fraction of proton energies carried by the incoming partons in collinear approximation.

Tests on the EW schemes and WWZ coupling
For (2 → 2) process, resummation of higher order effects into effective couplings is well established. In (2 → 4) case, care is necessary, one may destroy gauge cancellations where matching of Z emissions from quark lines with the ones of the t-channel W must be preserved. In Fig. 7 we demonstrate results, where effective sin 2 θ W is used in otherwise G F scheme, one can see that varying arbitrarily of WW Z coupling by ±0.05 bring marginal effects only, even for the q q,qq processes, chosen to maximize the relative effect of WW Z coupling mismatch. The effect is negligible for the shown, most sensitive kinematical distribution studied. The estimate of the average polarisation remains unchanged. This is an expected result as for our amplitudes condition sin 2 θ W = 1 − M 2 W /M 2 Z is in principle not needed for gauge cancelation.

Numerical results
Once we have completed our technical tests, and gained confidence in the functioning of the 2 → 4 extension of TauSpinner algorithms, let us turn to presentation of numerical results. In spite of a limited scope of the present version, like lack of the loop-induced gluon coupling to the Higgs boson, or subprocesses with b-quarks, TauSpinner can already be used as a tool to obtain numerical results of interest for phenomenology. Note that b-quarks as final jets can be tagged, and should thus be treated separately, while the contribution from the b-quark PDFs is rather small. Possible applications of TauSpinner are presented below.  Figure 7: Distribution of ∆ j j, reweighted to the one corresponding to WW Z coupling (internal MadGraph5 notation GC_53) multiplied by factor 0.95 (right) and 1.05 (left), shown for q q,qq (top) and qq (bottom) Drell-Yan processes. Table 4: Comparison of the τ-lepton polarisation in ττ j j events, calculated using TauSpinner weight wt spin of (2 → 2) and (2 → 4) processes and G F EW schemes with sin 2 θ W = 0.22222. Required invariant mass of the τ pair of m Z ± 10 GeV and low threshold on outgoing partons transverse momenta, p T > 1 GeV. Rows of the Table correspond to different subsets of events generated with MadGraph5, selected accordingly to flavours of incoming partons. TauSpinner algorithm is not using this information and the average of all possible configurations is used. In case of the last collumn, we restrict the average to the ones actually used for the selected subset of events.

Average τ lepton polarisation
For calculation of weights earlier versions of TauSpinner used the elementary (2 → 2) parton level qq(gg) → Z/γ/(H) → τ − τ + amplitudes factorized out from the complex event processes. This approach can now be verified with the explicitly implemented (2 → 4) matrix elements when two hard jets are present in the calculation of the amplitudes. The physics of interest is the measurement of the Standard Model Higgs boson properties in decays to the τ leptons and its separation from the Drell-Yan background of τ-pair production. We start by confirming the overall consistency of the calculations, comparing results from (2 → 2) and (2 → 4) calculations on inclusive ττ j j events, with τ-pair around the Z-boson mass peak, but with very loose requirements on the accompanying jets, p jet T > 1 GeV. In Tables 4 and 5, the estimated average polarisation is shown using matrix elements for (2 → 2) and (2 → 4) for four categories of hard processes and for cuts selecting events at the Z peak or above. For the (2 → 4) implementation shown is also the difference when estimating polarisation using an average of all hard processes, or for only specific category.
To verify that not only the calculation of spin averaged amplitudes, but the contributions from specific helicity configurations are properly matched between (2 → 2) and (2 → 4), we have checked the E π /E τ spectra in the τ ± → π ± ν decays. This variable is sensitive to the polarisation of the ττ system and longitudinal spin correlations. To introduce spin effects to the sample, otherwise featuring non-polarized τ decays, we have used weights calculated by TauSpinner. The spin weight distribution, the visible mass of τ's decay products combined and the energy fraction carried by the π ± in τ → πν decays are compared for two different EW schemes in Fig. 8.
To emphasize possible differences between using the (2 → 2) or (2 → 4) matrix elements for calculating spin weights for ττ j j events, we have applied simplified kinematic selection inspired by the analysis of [2], called in the following VBFlike selection: transverse momenta of outgoing jets above 50 GeV; pseudorapidity gap between jets, |∆η j j | > 3.0; Transverse momenta of outgoing τ leptons of 35 GeV and 30 GeV, respectively and pseudorapidity |η τ | < 2.5. It is also required that the invariant mass of the τ-lepton pairs and jj pair is above the Z-boson peak. Results for the average polarization, are shown in Table 6. and the energy fraction of the decaying τ lepton carried by π ± in τ → π ± ν (bottom-right), weighted with (2 → 2) and (2 → 4) matrix elements and for different EW schemes. Table 6: Polarisation of the τ-lepton in ττ j j events, calculated using TauSpinner weight wt spin and (2 → 2) and (2 → 4) processes and EWSH=1 scheme with sin 2 θ W = 0.22222. For this comparison the initialisation of (2 → 2) process was also adopted to EWSH=1 scheme. Required is the invariant mass of the τ-pair and j j-pair above 120 GeV and VBF-like selection (see text).
The VBF-like selection enhances contributions from qq,qq → ττ j j processes to about 25% of the total cross section. The highest discrepancy found between the predicted τ lepton polarisation with (2 → 2) and (2 → 4) matrix element is at the level of 4% in absolute value, being relative 10% of the polarisation. Using the average (2 → 4) matrix element i.e. assuming that the initial state is not known, reduces the discrepancy by factor 2. For the polarisation averaged over all production processes the difference between (2 → 2) and (2 → 4) matrix element is at the level of 1.0 -1.5% in absolute value, which is only 2% relative effect.
The above results indicate strongly that TauSpinner in the (2 → 2) mode is sufficient for the evaluation of spin effects observable in τ decays. The (2 → 4) mode is useful mainly for validations or systematic studies.
Please note, that results were obtained with the G F -on-shell scheme EWSH=1, thus are different from physically expected values. Let us continue now with the discussion of typical initializations used in calculations of matrix elements for (2 → 4) processes.

EW scheme dependence
In initialization of programs like MadGraph5, the tree level formula for weak mixing angle, sin 2 θ W = 1 − M 2 W /M 2 Z = 0.222246, is often used following the EWSH=1 or EWSH=3 schemes described previously. This theoretically motivated choice is quite distant from the sin 2 θ e f f W = 0.23147 describing the ratio of vector to axial vector couplings of Z-boson to fermions and which is used in the EWSH=2 scheme. The LO approximation used in MadGraph5 initialization and in our tests so far, can not be used for the program default initialization. We are constrained by the measured values of the M W , M Z and sin 2 θ e f f W and that is why the EWSH=4 scheme is chosen as a default.
One must keep in mind, that the τ-lepton polarisation in Z-boson decays is very sensitive to the scheme used for the electroweak sector. Table 7 gives numbers for the average polarisation in the case of G F and effective EW schemes. Results for using (2 → 2) and (2 → 4) matrix elements coincide within statistical error, for event sample with rather loose kinematical cuts. On the contrary, results of calculations strictly following the G F scheme are off by 50 % with respect to experimentally measured value, -0.1415 ± 0.0059, see Table 7. This must be taken into account if results are compared with the data, as it was done in the LEP times [24].

Summary and outlook
In this paper, new developments of the TauSpinner program for calculation of spin and matrix-element weights for the previously generated events have been presented. The extension of the program enables the calculation of spin and matrix-element weights with the help of (2 → 4) amplitudes convoluted with parton distribution functions. Required only are kinematical configurations of the outgoing τ leptons, their decay products and two accompanying jets.
The comparisons of results of the new version of TauSpinner, where matrix elements feature additional jets, and the previous one where the Born-level (2 → 2) matrix element is used, offer the possibility to evaluate systematic errors due to the neglect of transverse momentum of jets in calculating spin weights. We have found that for observables sensitive to spin, the bias was not exceeding 0.01 for sufficiently inclusive observables with tagged jets.
Numerical tests and technical details on how the new option of the program can be used were discussed. Special emphasis was put on spin effects sensitive to variants for SM electroweak schemes used in the generation of samples and available in initialization of TauSpinner. The effect of using different electroweak schemes can be as big as 50% of the spin effect and can be even larger for angular distribution of outgoing τ leptons. For the configurations of final states with a pair of jets close to the W mass the effect can be also high, up to 40%. For applications of TauSpinner, we recommend the effective scheme leading to results on τ polarisation, m Z and m W close to measurements.
The phenomena of τ decay and production are separated by the τ lifetime. This simplifying feature is used in organizing the programs. As a consequence for the generated Monte Carlo sample, different variant of the electroweak initialization may be used for generation of τ lepton momenta and later, for implementation of spin effects in the τ's decays. Such a flexibility of the code may be a desirable feature: TauSpinner weight calculation can be also adjusted to situation when the matrix element weight and spin weight have to be calculated with distinct initializations.
Numerical results in the paper were obtained with the help of weights. Not only spin weight, but also the production weight has been used to effectively replace the matrix element of the generation. This feature, introduced and explained in Ref. [8], was targeting an implementation of anomalous contributions. However, its use can easily be adopted for studies of the electroweak sector initialization. This helps to get results quicker thanks to correlated sample method. It provides technical advantage for the future, namely the possibility for the use of externally provided matrix elements or initialisation of EW schemes.
In tests discussed in this paper we have used MadGraph5 generated events for the pp → τ + τ − j j process. The incoming partons were distributed according to PDFs, but in most cases neither the transverse momentum of incoming state nor additional initial state jets were allowed. We will return to this point in the future, with greater attention. We may also be able to extend, with the help of the program developed for the present paper, the work on factorization of the effective Born of (2 → 2) configuration. This, in turn, will help to check the factorization of additional p T activity of our parton parton → τ + τ − j j hard processes from the pp collision. The special case of processes with a single hard jet in the final state and its corresponding matrix elements will be also useful for such tests and we plan to return to such topic in the near future.
The program is now ready for studies with matrix elements featuring extensions of SM amplitudes. Systematic errors have been discussed. Let us stress, that results of such studies depend on the definition of observable and need to be repeated whenever new observables or selection cuts are introduced. In the evaluation of impact of new physics or variation of SM parameters on experimentally accessible distribution it is necessary to compare results of calculations which differ by such changes. The Monte Carlo simulations are used, whenever detector acceptance and other effects are to be taken into account.
A Comments on the code organization and how to use it.
In this Section, we collect information on how to use the (2 → 4) option of the TauSpinner program. We will concentrate on aspects, which are important to demonstrate the general scheme and organization of the new functionality of the code. We assume that the reader is already familiar with previous versions of TauSpinner or Refs. [9,10]. 9. Several scenarios (models) of the hard process for calculating corresponding spin weight (note that at the same time the weight for the production matrix elements is calculated), are possible. At the initialisation step the choice is made and the technical internal parameter KEY is set. We give below some details: • KEY=1 for the Standard Model Higgs process, matrix elements explained in Section 3.1 is used.
• KEY=0 for non-Higgs Drell-Yan-like processes, matrix elements explained in Section 3.1 is used.
• KEY> 1 is reserved for non-standard calculations, that is when nonSM=true option is used 11 , and the matrix element calculations are modified with the routines provided by the user. Provisions with KEY=3 have been prepared for the non-standard Higgs-like production process and KEY=2 for the Drell-Yan-like. The required choice is made implicitly, at the initialisation step when setting the pointer to the user provided function vbfdistrModif and selecting initialization variable nonSM2=1 which will set internal global variable nonSM=true. The result of default calculations will be passed to vbfdistrModif function to be overwritten with user-driven modifications to the matrix elements, without the need of re-coding and recompiling standard TauSpinner library. • The internal parameter KEY in general does not need explanation. However, as it is passed latter to the methods for calculating α s or matrix elements, which may be replaced from user main program by re-setting the pointer, documentation was necessary.
10. The method double getTauSpin() returns helicities attributed to τ leptons on the statistical basis. It is the same method as already implemented for the (2 → 2) case. 11. The value which is returned by the method double getWtNonSM() depends on the configuration of two flags: nonSM and relWTnonSM. The above discussed weight features matrix elements squared and summed over spin degrees of freedom. It has similar functionality as already implemented for the (2 → 2) case. It is supposed to supplement the spin weight WT of TauSpinner. In general, spin weight differs for the SM and nonSM calculation. The ratio of these two has to be used for modifying decay product kinematic distributions. Finally let us point out, that also the helicity of τ's will be attributed at the nonSM step of the calculation, corresponding to the chosen nonSM model.
Let us bring some further points on the details of the use of the example.
• This example has been prepared to read events in HepMC format. An additional tool lhe-to-hepmc.exe convertin LHE event to HepMC has been provided as well.
• To read events form the data file, the method read_particles_for_VBF stored in file TAUOLA/TauSpinner/examples/example-VBF/read_particles_for_VBF.cxx is used. It is invoked as follows: int status = read_particles_for_VBF(input_file,p1,p2,X,p3,p4,tau1,tau2,tau1_daughters,tau2_daughters); which reads consecutive event and retrieves the following information: four-momenta of incoming and outgoing partons, denoted as p1,p2, and p3,p4, respectively; four-momenta of outgoing τ leptons, tau1, tau2 and tau1_daughters,tau2_daughters which stand for lists of decay products (their four-momenta and PDG-id's) and, if available, also the four-momentum of an intermediate resonance X and its PDG-id. It returns status=1 if no event to read was found (this is specific to the method chosen for reading the events and is used in the user program only).
Let us stress, that the above interface to read event record is not a part of TauSpinner library. It is used in the demonstration program and is adopted to the particular conditions. It is expected to be replaced by the user with the customized one. Implementation of such a method must match conventions for the format and information on different particles in the stored event. For example, for formats (eg. HepMC [25] or lhe [22]) distinct conventions are used in Monte Carlo generators: i.e. for relations among particles and intermediate states (resonances) which may be explicitly written into event record or omitted. The same is true for production/decay vertices, status codes, etc. With a variety of conventions used, it can be highly error-prone and lead to necessity of non trivial implementations in the code, like can be seen in read_particles_for_VBF method provided in the distribution tar-ball. See eg. [15,26] for discussion of similar difficulties in other projects. In fact, for TauSpinner algorithms this is less an issue, as the information on intermediate and incoming states is not used, even though it can be very useful for testing purposes.
The information required by TauSpinner algorithms for (2 → 4) processes which must be read in, is limited to: fourmomenta of outgoing jets and of τ ± including all their decay products and PDG identifiers of τ ± and their decay products. The four momenta of incoming partons and the intermediate resonance are not needed but can be used for tests. The imminent next step is to exploit also, if available, information on the four-momentum of the intermediate Z/γ/H state (or any other non-standard resonance) which decays to τ ± lepton pair as well, as they can be used to tackle the effect of the QED bremsstrahlung in its decay, similar as it was done for the TauSpinner algorithms in (2 → 2) case.

A.2 Initialisation methods
• Matrix elements: The TauSpinner library includes codes for calculation of matrix elements squared for all (first two families) parton level cross sections of (2 → 4) processes. Use of this library functions can be over-loaded, with the user's own matrix elements implementation by providing respective function vbfdistrModif(...). Its usage is activated at initialization with command TauSpinner The following arguments are passed to this function, and this must be obeyed in its declaration: -The first four arguments I1,I2,I3,I4 denote PDG identifiers [27] of incoming and outgoing partons (for gluon ID=21) . -The following two H1,H2= ± 1 denote helicities of outgoing τ + and τ − .
-Matrix P [6][4] encapsulates four-momenta of all incoming/outgoing partons and τ ± leptons. They are for massless partons and for massive τ leptons. Energy momentum conservation is required at the double precision level. The partons are not expected to be in the phase-space regions close to the collinear/soft boundaries.
-The parameter KEY=0 is reserved for the SM default processes of Drell-Yan-type (all Feynman diagrams included, but the ones with H → τ + τ − ), while KEY=1 for SM processes with the Higgs production and its decay to τ-lepton pair. In these two cases vbfdistrModif is not activated. For KEY=2,3 the SM calculation (again respectively for Drell-Yan and Higgs processes) is performed first and the result is passed into vbfdistrModif() where it can be just modified, before being used for final weight calculations. The user may choose to modify the value of the default calculations for all or only for subset of processes involved. This is why complete information of the initial and final state configurations is exposed. If a completely new calculation is to be performed using the above method, then it is advised to use options KEY=4,5, so that the Standard Model calculation will be avoided (to save CPU) and result=0 will be passed to vbfdistrModif(). The KEY=4,5 is reserved for optional use of vbfdistrModif().
• Electroweak schemes: In Table 3 options of initialization for the EW schemes implemented in TauSpinner  The EWSH_ref will set initialization as used for the default calculation, and EWSH_variant for the reweighting with modified amplitudes. The choices 1, 2, 3, 4 correspond to EWSH=1, EWSH=2, EWSH=3, EWSH=4 respectively. The default EWSH=4, as explained in the main text, leads to correct τ lepton polarisations and angular distributions. As it causes at tree-level inconsistencies in the calculation of the WW Z coupling, we provide an additional option, EWSH=5, for which parameter setting as for EWSH=4 is used, but with the WW Z coupling modified by 5%. It can be used for testing sensitivity of the analysed distributions to the missed higher order corrections to the WW Z coupling. For more discussion, see Section 3.3.
• PDFs and α s : Any PDF set from LHAPDF5 library [28] can be used for calculating spin weight. The choice can be configured by setting, string name="cteq6ll.LHpdf"; LHAPDF::initPDFSetByName(name); The choice of renormalization and factorization scales (imposed is case of µ F = µ R ) can be set with the help of the following command: int QCDdefault=1; // QCD scheme to be used for default 2 ->4 calculation. int QCDvariant=1; // QCD scheme to be used in optional matrix element reweighting (nonSM2=1). setPDFOpt(QCDdefault,QCDvariant); The choice can be different for the default (SM) calculation and the variant one (nonSM=true), see Appendix A.1, point 9. The Q 2 evolution and starting value of α s used in PDF's is internally defined by the LHAPDF5 library. For the matrix element calculations we do not impose consistent definition of α s but it can be enforced by the user, see next point. As a default, we fix starting point at α s (M Z ) = 0.1180 value and evolve it with Q 2 with a simple formula of Eq. (4).
• User own α s in matrix element calculation: User can supersede the simple, leading logarithmic function provided by us for α s (Q 2 ) used in the matrix element calculation (Eq. (4)) with his preferred one, and pass it to the program. The function calculating α s has to have the following arguments: alphasModif(double Q2,int scalePDFOpt, int KEY) In alphasModif one can also use directly a method LHAPDF::alphasPDF(sqrt(Q2)) of LHAPDF5 library [28], the same assuring consistency between value of α s in the matrix element and the structure functions. Such function can be used by executing set_alphasModif(alphasModif); . An example of such setup has been provided in example-VBF.cxx program.

A.3 Random number initialization
In most of the calculations the TauSpinner algorithms are not using random numbers. However, there are two exceptions. In both cases random generators from TAUOLA are used, see Appendix C.12 of Ref. [15].
• The helicity states attribution uses Tauola::RandomDouble. It should be replaced by the user, with the help of Tauola::setRandomGenerator(double (*gen)()) method and then properly initialized, with distinct seed for each parallel run. In our example program the actual command is Tauola::setRandomGenerator( randomik ); • If the read_particles_for_VBF.cxx code is required to generate τ decays, then a second random generator, coded in FORTRAN has to be also initialized with distinct seed for each individual parallel run: Tauola::setSeed(int ijklin, int ntotin, int ntot2n).

A.4 Main program -an example
The following files are prepared for the user prototype program in the TAUOLA/TauSpinner/examples/example-VBF directory • The user example program example-VBF.cxx.
• The prototype method read_particles_for_VBF.cxx to read in events stored in HepMC format is prepared specifically for MadGraph5 generated events.
• The README file which contains auxiliary information.
Only the program example-VBF.cxx is generic, and does not depend on the specific environment for event generation. This is why we provide an extract from this code below. For the TauSpinner library to work, the τ decay products must be present in the event. In case they are absent, like e.g. in events generated with τ's as final states in MadGraph5, we prepared settings for their decays in Tauola library using the mode of not-polarised τ decays and Tauola universal interface. Such additional processing is implemented in read_particles_for_VBF.cxx code. In our demonstration program for TauSpinner spin correlations between τ leptons are then introduced, using (2 → 4) matrix elements and calculating respective spin weight. The purpose of the example is to demonstrate the default initialisation of the TauSpinner program and a flow of the main event loop.
The necessary changes were introduced, and from now on the first argument ID passed by TauSpinner library to the userdefined function nonSM_adopt activated by the pointer: set_nonSM_born( nonSM_adopt ) of /TAUOLA/TauSpinner/examples/tau-reweight-test.cxx denotes the incoming parton flavour PDGid. In constrast, in earlier version of TauSpinner library [8], it was possible to invoke user-defined function nonSM_adopt activated by the pointer in set_nonSM_born( nonSM_adopt ) The first argument of this method was passing to the user function the information if incoming parton was up-or down-type quark only, without specifyingits family affiliation.

B Tests of reweighting the differential cross-sections for Drell-Yan-like processes
In this Appendix we show in Figs. 9 and 10 a complete set of kinematic distributions validating implementation of (2 → 4) non-Higgs Drell-Yan-like processes. We split pp → ττ j j events into four groups, depending on the initial partons, see Table 1 for definition of parton level processes. We use the implemented (2 → 4) matrix elements to calculate per event a weight, wt C→D prod = dσ D /dσ C , see Eq. (6), defined as a ratio of the cross-sections for events of groups C and D. The expression is similar to Eq. (3) except that the sum is over subprocesses which belong to the chosen groups C or W . We apply wt C→D prod to events from the group C and compare both the absolute normalisations and shapes of the re-weighted distributions with the distributions of events from the group D.
These tests were done on the large statistics samples and have been repeated between each groups of processes and within groups between subgroups. The achieved agreement between the reference and re-weighted distributions validates the correctness of the implemented matrix elements.

C Tests of reweighting differential cross-sections for Higgs boson production
Similar tests, as discussed in Appendix B, have been repeated for the pp → H(→ ττ) j j processes. Results are shown in Figs. 11 -14. Very good agreement between the reference and re-weighted distributions is observed, both for shapes and relative normalisations.

D Optimalization of interface to Standard Model matrix element calculation
The steering function for calculating (2 → 4) matrix elements squared REAL*8 FUNCTION VBFDISTR (ID1,ID2,ID3,ID4,HH1,HH2,PP,KEYIN) is coded in FORTRAN and stored in the VBF_distr.f file. Before invoking calculation of particular matrix element squared of the Standard Model, it performs several steps of filtering to speed up the CPU needed for numerical calculations by setting matrix element squared to zero without calculation for the cases when configuration of partonic PDG identifiers for incoming and outgoing partons imply that it is the case. The following conditions are consecutively checked (strictly in the given order 13 ). Each condition must be passed to go to the next one, and finally to invoke the matrix element calculation.  Figure 9: Shown generated g q → ττ j j (thin red line) after reweighing to q q (qx qx) → ττ j j (red points). Reference q q (qx qx) → ττ j j distribution shown with black line.  Figure 10: Shown generated g q → ττ j j (thin red line) after reweighting to q q (qx qx) → ττ j j (red points). Reference q q (qx qx) → ττ j j distribution shown with black line.  Figure 13: Shown example of tests distributions for generated process q qx → H(→ ττ) j j (thin red line) after reweighting to q q, qx qx → H(→ ττ) j j process (red points). Reference q q, qx qx → H(→ ττ) j j distribution shown with black line.