Production of tau tau jj final states at the LHC and the TauSpinner algorithm: the spin-2 case

The TauSpinner algorithm is a tool that allows to modify the physics model of the Monte Carlo generated samples due to the changed assumptions of event production dynamics, but without the need of re-generating events. With the help of weights $\tau$-lepton production or decay processes can be modified accordingly to a new physics model. In a recent paper a new version TauSpinner ver.2.0.0 has been presented which includes a provision for introducing non-standard states and couplings and study their effects in the vector-boson-fusion processes by exploiting the spin correlations of $\tau$-lepton pair decay products in processes where final states include also two hard jets. In the present paper we document how this can be achieved taking as an example the non-standard spin-2 state that couples to Standard Model particles and tree-level matrix elements with complete helicity information included for the parton-parton scattering amplitudes into a $\tau$-lepton pair and two outgoing partons. This implementation is prepared as the external (user provided) routine for the TauSpinner algorithm. It exploits amplitudes generated by MadGraph5 and adopted to the TauSpinner algorithm format. Consistency tests of the implemented matrix elements, reweighting algorithm and numerical results for observables sensitive to $\tau$ polarization are presented.


Introduction
With increasing statistics collected by the LHC experiments the interests to explore final states with τ -leptons gain importance. Because of the high mass and their decays, τ -leptons may provide a sensitive window to physics beyond the Standard Model predictions. The TauSpinner algorithm, started with Ref. [1], provides a powerful tool to investigate characteristics of final states with τleptons due to modifications in underlying physics models. This is obtained with the help of weights attributed to each event from collision data or Monte Carlo generated, and thus without repeating the detector response simulation with each variant of the physics model. An approach where physics assumption variation can be introduced with weights is useful for many modern data analysis techniques. The first version of TauSpinner reweighting found its application in the domain of Standard Model measurements [2] but also for New Physics limits established for simple 2 → 2 parton level processes [3]. Later, in [4,5], TauSpinner was found useful for discussion of CP sensitive massively multi-dimensional observables in a e-mail: Z.Was@cern.ch the frame of Machine Learning techniques [6]. In a recent paper [7] an extended version of TauSpinner 2.0.0 was presented which now includes hard processes featuring tree-level parton matrix elements for production of a τ -lepton pair and two jets. It was prepared as a tool to be helpful for studying spin effects in processes of Standard Model and searches of New Physics, like in Refs. [8,9]. It found also tempting applications in the domain of implementation and discussion of variants for Standard Model electroweak calculation schemes used in simulation programs [7], and in experimental applications for Standard Model measurements [10][11][12][13].
Before discussing TauSpinner as a tool for studying observables of New Physics let us briefly recall some virtues of the TauSpinner algorithm. Since τ -leptons cannot be observed directly due to their short life-time with more than 20 different decay channels, each with somewhat distinct signature, recalculating and reanalyzing observables involving τ decays is time consuming. However the τlepton spin polarization can be inferred from their decays, contrary to the case of electron or muon signatures. Spin effects can provide a better insight into the nature of the underlying physics. Therefore efforts to explore these phenomena are worth pursuing. TauSpinner allows to greatly simplify the task of exploring the experiments' sensitivity. Evaluation of measurements significance due to different New Physics models can be performed with the help of event weights. Technical aspects of the algorithm in case of configurations with two jets accompanying a τ -lepton pair have been covered in [7].
The purpose of the present paper is to document how the user can apply the TauSpinner algorithm to his/her physics model. To this end, we take as a case study a non-standard spin-2 object coupled to SM particles. We analyze its production in proton-proton collisions and decay to a τ -lepton pair, addressing also the question to what extent τ polarization can be exploited to investigate its nature. We demonstrate how TauSpinner can facilitate such studies with the help of matrix elements for that model (or any other, provided by the user). Corresponding weight can be calculated and applied to each event of samples with full experiment simulation chains. As it is practically impossible to repeat simulations with detector response effects included for each new physics hypothesis, our procedure is beneficial and may be the only available option for the general use despite its limitations. TauSpinner algorithm can also be applied on measured data events, e.g. in the context of embedded τ lepton techniques [14].
The paper is organized as follows: Section 2 and Appendix A provide details of TauSpinner which were not discussed in Ref. [7], or treated very briefly, but are of importance for the case of New Physics models. Section 3 documents details of the tree-level matrix elements used for the calculation of spin-2 object exchange amplitudes which are later passed for weights calculation in pp → τ τ jj events. The implemented functionality is based on automatically produced FORTRAN code from Mad-Graph5 package [15] similarly as for processes of the Drell-Yan-type and of the Standard Model Higgs boson production in vector-boson-fusion processes (VBF). We explain details of the modification which we have introduced to the code of amplitudes generated with MadGraph5 (version MG5_aMC_v2.4.3). Section 4 and Appendix B are devoted to numerical results. First, tests for fixed kinematic configurations are recalled. Later, definitions of observables are given, and some distributions are presented. In Section 5 numerical results sensitive to τ polarization are presented taking the single-prong decay τ ± → π ± ν channel as a spin analyser. Section 6 summarizes and concludes the paper. 2 The TauSpinner weight wt A→B prod TauSpinner does not provide methods to generate pp collision events. Therefore, the necessary input for TauSpinner consists of a series of events, which could be of a process different than the required one but with the same outgoing final states. The events must contain information on the four momenta of (two) outgoing jets and τ -leptons with their decay products, which is necessary for the calculation of the hard process matrix elements. Flavours of incoming/outgoing partons are determined by the algorithm -there is no need to read them from the generated events. The sum over all possible configurations, weighted with PDFs, is performed. On the other hand, the information on the decay products of τ -leptons is needed for the evaluation of spin effects. Using this input, the value of the corresponding matrix elements can be calculated on the event by event basis and in particular the corresponding spin weight. As discussed in detail in Ref. [7], for each event where j stands for a quark, antiquark or a gluon, the algorithm calculates the weight which represents, for a given phase-space point ({p}) = (p 1 , p 2 , p 3 , p 4 , p τ + , p τ − ), the ratio due to the matrix element used in the generation of the sample for process (A) and the matrix elements corresponding to a New Physics model 1 (B). The evaluation of the weight in Eq. (2) requires the knowledge of contributions from all possible parton level configurations (ijkl) weighted with parton density functions f A/B i/j (x) and flux factors Φ i,j f lux . The sums run over both gluons and quark flavours alike. Note however, that the flavour is passed to the user provided matrix element routine and flavour dependence can be introduced there. For the details and explanation of a notation used in formula (2) we refer to [7]. For the purpose of calculating wt A→B prod we sum over all possible helicity configurations of outgoing τ -leptons. An event generated for the process (A), when weighted with wt A→B prod becomes an event of the process (B). Spin effects in τ decays have to be introduced separately, with TauSpinner main spin weight W T , as explained in Appendix A.
The following details need to be stressed when selecting a suitable process (A) given the target process (B). For the narrow resonance, like the Higgs state, the value of matrix elements vary greatly with the invariants built from the final state four-momenta. Therefore the numerical stability needs to be kept in mind. The TauSpinner algorithm must reconstruct invariant mass of the resonance with the precision better than 1-2 MeV from the fourmomenta of final state particles, whose energies may lie in the range of TeV. Double precision may be needed since otherwise some invariants may be inappropriately evaluated due to simple computer rounding errors. A user has to assure that the reweighting indeed works in the interesting regions of the phase-space. In particular, that the phase-space is populated for both processes (A) and (B) with not too massively distinct distributions, and that the distribution enhancements due to intermediate resonances or collinear or soft singularities have similar (matching) structure. Listed above checks require the hard processes information only. 3 New Physics model of (2 → 4) process As a case study we consider a simplified model of a massive gauge singlet spin-2 object X coupled to the SM gauge bosons. We use this model to demonstrate how to prepare and test external matrix element to be used by TauSpinner algorithm.
Scenarios with spin-2 objects have been already intensively studied in the literature in the context of LHC phenomenology [17][18][19], though none of the studies was dedicated to the analysis of X decays into τ final states. Note however that for a general study of a "Higgs"-like resonance and its parity in vector-boson-fusion processes with a τ pair as a decay product, experimental results are becoming available [20,21].
In Ref. [3] we studied a Drell-Yan-like production of τ 's through a hypothetical spin-2 object X. Building on our previous work, we study now the X production in the VBF topology, followed by X → τ + τ − decay. We start by extending the Lagrangian of Ref. [3] by a set of gauge invariant dimension 5 operators, coupling the field X to gauge boson field strength tensors B, W and G as where group indices are implicitly summed over (where appropriate). The parameter F , set to 1 TeV, is introduced to keep the coupling constants dimensionless. Note that we are agnostic on the origin of the state X, in particular we do not claim it is connected to gravity. Hence we do not couple it to the entire energy momentum tensor and couplings g X are kept as free parameters. This is in contrast to, for example Ref. [22], where the X field is coupled to the energy momentum tensor of quantized SM. After the electroweak symmetry breaking, operators in Eq. (3) generate vertices with couplings of X to photons, W ± 's, Z's and gluons; the explicit formulas for those couplings can be found in [19]. Since in this work we focus on technical aspects of incorporating the couplings of X to the EW gauge bosons, for numerical tests of the correctness of the Matrix Element implementation, we set g XBB = g Xgg = 0. Relevant diagram topologies are shown in Fig. 1: for the VBF process (Fig. 1a) and the X-strahlung process (Fig. 1b).

Generating matrix-element code using MadGraph5
The extension of the SM by spin-2 field coupled to the gauge fields as in Eq. (3), including also coupling of the X field to quarks and τ -leptons from [3], is encoded into a FeynRules [23] model. The FeynRules model file, together with its UFO output [24], is available in supplementary materials of the arXiv version of this reference. The UFO model is used to generate squared matrix elements using MadGraph5, employing the spin-2 support of the HELAS library [25]. This is done with the following set of commands (a) import model spin2_w_CKM_UFO (b) by default, "multiparticles" containers already include all massless partons p = g u c d s u~c~d~sj = g u c d s u~c~d~s( c) generate spin 2 matrix elements generate p p > j j x QED<=99 QCD<=99 NPgg<=99 NPqq<=99 NPVV<=99, x > ta+ ta-(d) write the output to disk in MadGraph's standalone mode using output standalone "directory name" NPgg, NPqq and NPVV parameters control the maximum number of g Xgg , g Xqq and g XW W , g XBB couplings, respectively. Limiting them by 99 effectively means that their number is not restricted. The model includes the CKM matrix in the Wolfenstein parametrization. As was stated above, for numerical tests we restrict ourselves setting g XBB = g Xgg = 0, though we stress again that the matrix element, coded as an example user process, contains all of them, see Appendix B for actual initialization of coupling constants.

Integrating matrix-element code into TauSpinner example
The matrix element code is based on automatically produced FORTRAN subroutines by MadGraph5 package, similarly as it has been done for processes of the Drell-Yantype and of the Standard Model Higgs boson production in vector-boson-fusion (VBF)/Higgs-strahlung processes [7]. In the spin-2 case they have been also manually modified and adapted to avoid name clashes. This technical complication is a consequence of the fact that C++ user function for the spin-2 matrix element calls FORTRAN code created by MadGraph5. We therefore can not profit from namespace functionality of C++ as a natural solution to this problem. Some name changes are necessary, as explained below. The corresponding code is stored in the directory TauSpinner/examples/example-VBF/SPIN2/ME. The generated codes for the individual sub-processes are grouped together into subroutines, depending on the flavour of initial state partons, and named accordingly. For example, Table 1. List of implemented processes contributing to the spin-2 X particle production, grouped into categories which differ by flavours of incoming partons. For each category, the names of FORTRAN files calculating squared matrix elements, for given flavour configuration of incoming partons, are given in the second column. Examples of processes in each category are given in the last column.

Category of
Corresponding FORTRAN files Processes Matrix Elements SUBROUTINE DCX_S2(P,I3,I4,H1,H2,ANS) encompasses the X production processes initiated by the dc partons. We follow our previous convention [7] where symbol X in the subroutine or internal function name after the letter U,D,S or C means that the corresponding parton is an antiquark, i.e. UXCX corresponds to processes initiated byūc partons, while GUX to processes initiated 2 by gū. The S2 stands explicitly for the production of spin-2 X state. The input variables are: real matrix P(0:3,6) for four-momenta of incoming and outgoing particles, integers I3,I4 for the Particle Data Group (PDG) identifiers for final state parton flavours and integers H1,H2 for the outgoing τ helicity states. 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, i.e. indices/helicities H1 and H2. b) Moreover, since the subroutines and internal functions generated by MadGraph5 have the same names for all sub-processes, namely SMATRIX(P,ANS), the names had to be changed to be unique. As an example, the subroutine name for the subprocess ud → cd X, X → τ + τ − was changed to UDX_CDX_S2(P,H1,H2,ANS). These subroutines will be called by subroutine UDX_S2(P,4,-1,H1,H2,ANS). c) 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 configuration of outgoing 2 X in this context should not be confused with the spin-2 field X.
partons is unknown and it takes into account both possibilities: thus a compensating factor 1+δ l,k 2 has to be introduced due to the way of organizing the sum in Eq. (2) and in Ref. [7]. d) For calculation of matrix elements MadGraph5 is using ALOHA functions [26] stored in FORTRAN subroutines. Since some of these functions have originally names identical to functions in the TauSpinner source code for the implementation of the Standard Model VBF/Higgs-strahlung production, names of those functions have to be modified also to avoid any name conflicts. Therefore ALOHA functions stored in TauSpinner/examples/example-VBF/ SPIN2/ME/Spin2_functions.f are changed by adding "_S" suffix to the original names of subroutines, for example FFV4_0 is changed to FFV4_0_S. Table 1 summarizes the naming convention for the files. At the parton level each of the incoming or outgoing partons can be one of the flavours:csūd g d u s c, with Particle Data Group (PDG) identifiers: -4, -3, -2, -1, 21, 1, 2, 3, 4 respectively. For processes with two incoming partons, two outgoing τ -leptons and two outgoing partons, there are 9 4 possibilities, most of them evaluating to 0 or obtainable one from another, by relations following from CP symmetries and/or permutations of incoming and/or outgoing partons.
For each point in the parton level phase-space, consisting of all incoming and outgoing four momenta as well as their flavours, depending on the user choice, one of two variants of processes (i.e. pairs of matrix elements) may be used by TauSpinner executable. That is: the Drell-Yan variant (standard, and user provided New Physics matrix elements) or Higgs-like variant (again standard, and user provided one 3 ).
Certain limitations need to be kept in mind. In practice, it is simply impossible to obtain statistically significant distribution of weighted events for the particular model under study in the region of phase-space where original sample is sparse or possibly no events are present at all. In particular, the mass and width of the Higgs-like resonance need to coincide (be close) to those of the Higgs. Also, the algorithm is expected to be used in regions of the phase-space where kinematic distributions of the original and New Physics models are not massively different.

Tests of implementation of external matrix elements
Once the user-provided external matrix elements are prepared, numerical tests are necessary if it indeed has been implemented properly into the TauSpinner environment. In the following we discuss such tests, using spin-2 matrix elements of Xjj production as an example. We start from the technical one and continue with more physics oriented ones. Finally we will demonstrate limitations of the method.

Test of matrix elements using fixed kinematic configuration
For checking the consistency of the implemented codes generated with MadGraph5 and modified as explained in section 3.2, we have chosen a single event with fixed kinematic configuration at the parton level. We have calculated the matrix elements squared for that event and for all possible helicity and parton flavour configurations, using the code implemented as user example. We compared results with 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 external matrix element implementation, we have explored the reweighting procedure by comparing a number of kinematic distributions obtained directly or reweighted with wt A→B prod from series of 10M events generated by MadGraph5 for X particle and Higgs boson. Samples were generated for pp collisions at 13 TeV using CTEQ6L1 PDFs. The mass of both X particle and Higgs boson was set to 125 GeV and the width to 5.75 MeV. The details of cuts and MadGraph5 initialization used for the sample generation are given in Ref. [16]. On the generated events the following further selections were applied: m jjτ τ < 1500 GeV, p τ τ T < 600 GeV 3 The prototype is implemented in the example, see Appendix A.1.  and m jj < 800 GeV (loose selection) to eliminate excessive weight regions of the phase space, or eliminating also Z → jj or W → jj resonance peaks 100 < m jj < 800 GeV ( tight selection).
Before commenting on the actual results let us point to the size of statistical errors 4 which reflect comparability of the H ( process A) and X (process B) samples. Errors are always larger than what could be expected from weight-one samples of the similar size. This effect can be understood better with the weight distributions shown in Fig. 2. In both cases of reweighting: from H to X (top panel) and X to H (bottom panel), one can observe a constant slope on this double logarithmic plot with clear sharp upper end. With such spectrum of weights statistically sensible calculation of the cross sections and distributions may be still possible. If a tail of events with Cross-section (pb)  Table 2. Cross sections for the generated H production process and after its reweighting to the X production (H → τ τ block), and for the generated X production and after its reweighting to H production (X → τ τ block); acceptances with no, loose or tight selections applied for generated and reweighted event samples are also shown.  Selection cuts: invariant mass of outgoing particles m τ τ jj < 1500 GeV, invariant mass of jets system 100 < m jj < 800 GeV and p τ τ T < 600 GeV. Note that statistical errors for the distributions obtained with reweighting (red points) are much larger than for the case of Fig. 3. This is predominantly due to small acceptance of X sample: 1.7% only. But the agreement with the reference distribution (black histogram) remains within statistical fluctuation (dominated by large weight events). Variables on the x-axes as explained in Section 4.2.
ever higher weights would continue to form when increasing size of samples, statistical errors would never decrease. This happens, for example, if in some sub-dimensionalmanifold of the phase space the matrix element has a zero. Then with increasing statistics, events closer and closer to this zero are generated, and feature larger and larger weights. Even though contribution of such events to the weighted distribution is formally finite and integrable, the error estimate of the Monte Carlo generated distribution will not get reduced with the increasing statistical sample.
The tests were performed on a set of kinematic distributions: pseudorapidity of outgoing parton j, rapidity of τ τ and jj systems, invariant mass of τ τ system, pseudorapidity of τ τ system, opening angle between jets, opening angle between τ 's, angle between incoming parton and outgoing parton in the rest frame of jets and angle between resonance and outgoing parton in the rest frame of jets.
Plots for all these variables can be found on the web page [16]. Here, in Fig. 3 and Fig. 4, we present only plots for: the difference of jet's rapidities ∆η jj , the invariant mass of the jet pair m jj , the transverse momentum of τ pair p τ τ T and finally the invariant mass of τ -pair and jetpair combined m τ τ jj . In each plot the distribution Ref, for the reference process, is shown as a black histogram while the red histogram is the original distribution of generated events which are reweighted using TauSpinner wt A→B prod weight to obtain the distribution 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  Fig. 4), the distributions feature larger statistical errors than in the case of H to X reweighting (Fig. 3). This is simply because tight selection cuts leave only 1.7% of X events due to eliminating configurations with small m jj . For some bins the reweighted distribution lies below target (black) distribution, whereas the ones with big errors tend to lie above. If similar feature appears when sample size is increased it points to the possibility that original distribution had a zero along some hyperspace. Nevertheless, if in distributions normalized to cross section the neighbouring bins have no deficit of content, then the reweighting algorithm can be still used.
The tests validating reweighting algorithm are completed with the ones monitoring overall normalizations (integrated cross sections). For our samples and initializations, the resulting cross sections are shown in Table  2. Reasonable agreement between cross sections obtained from the MadGraph5 calculation and with reweighting was obtained, see Table 2 where the first line in the H → τ τ block should be compared to the second line in the X → τ τ block and vice-versa. Such a study has to be repeated for each new matrix element implemented and whenever selection cuts are changed sizably.

On reliability of the TauSpinner reweighting approach
The TauSpinner reweighting method is atypical compared to methods used in other tools, like REPOLO [28] PHYTIA [29], SHERPA. [30] or MadGraph [31]. Let us explain what are the advantages and disadvantages behind such a choice.
The advantage of our method is that it does not assume any knowledge of the initial and outgoing partons and tau leptons beyond their four momenta. Therefore it can be applied directly to the experimental data, e.g. of the embedded τ samples. We have demonstrated that our reweighting method is reliable for the hard process matrix elements convoluted with PDF's. The disadvantage is that it does not address the issue that both the parton shower and hadronization do depend on color configuration as well as on flavours of partons. Once event is reweighted, the reshuffle between categories of different color and/or flavour content takes place, inevitably leading to biases.
In experiment simulation production files [32] color information for the so called truth entries is not stored. Even in the data formats prepared and agreed on by the community [33], such information, at best, consists of a connected tree, navigation inside of which retains information on the event history including the parents of unstable particles. There is an important caveat here: the event generators are modeling quantum processes, and the event record has the structure of a classical decay chain. It is inevitable that compromises must be made and difficulties can arise from an over-literal interpretation of the tree structure. For the color it means that at best the so called flow approximation is pre-imposed. Even for such partial information, there are no detailed commonly accepted rules how it should be stored, see for example Section 2.3 in [34] or Section 4.4.1 in the HepMC manual [35] or [36].
In practice, in experiment production files of detector response simulated events, information on intermediate quantum states is generally not available. Usually only the 4-momenta of partons and their flavours are stored. We are not in a position to affect these experiments choices. Multitude of arguments have been raised for such choices, including the fact that distinct generators prepared by theorists, provide such information in different manner, or that it makes data files unnecessary large. We can only address the question if any useful solution for reweighting may be designed 5 and what kind of restrictions it implies have to be kept in mind.
There are two simulation steps which depend on hard process configurations of flavours and colors: parton shower and later hadronization. It is well known that even mainstream Monte Carlo programs do not match in this respect sufficiently well the experimental data for all required phase space regions [37]. This is a complex issue which we can not exhaust 6 .
The discussion of the resulting systematic errors of our method is out of scope of the present paper. It would require evaluation of how mismatches of the color and flavour input for parton shower and hadronization translates into reconstructed jets from simulated detector responses. This in turn would require the use of experimental detector response codes, not available publicly. In general, the TauSpinner application domain is restricted to observables where details of jets, resulting from partons flavours or color are not of importance. This has to be kept in mind.
Let us point that our study examples of the previous sub-sections are for the cases, where starting and target distributions are massively different. In practical applications we expect TauSpinner to be used in configurations where new contributions to matrix elements are at the edge of observability.
If required, it is possible to apply TauSpinner in the flavour savvy manner. Possible solution may follow the method described in Appendix A. Contributions from distinct flavour configurations can be treated separately only for cases when in experiment production files the flavour configurations are stored, or can be unwinded.

Spin dependent characteristics
So far we were discussing observables relying on the kinematics of final states consisting of four momenta of τ leptons and accompanying two jets. Inclusion of τ decay products increase the phase space dimensionality substantially, making the analysis much more difficult, especially when dependence on selection cuts is taken into account (as observed in previous sub-sections).
In the following, we will present a few spin dependent results obtained for the H and X samples within the tight selection cuts. Using TAUOLA ++ [40] we supplement these samples with τ decays in the simplest possible mode τ ± → π ± ν with no spin effects included. Spin effects are introduced with the help of TauSpinner weights, which are calculated according to the production and decay kinematics (see Refs. [1,41] for the spin weight definition). Figure 5 shows the spin weight histograms for the H and X samples. In both cases the spin weights are calculated first using the matrix element for X productions as described in Ref. [3], that is featuring effective Born 2 → 2 kinematic (open red circles), and compared with the new calculation in which amplitudes featuring two jet kinematics are taken into account (blue full circle points) 7 . In both cases the same X − τ τ couplings were used. As expected (see Eq. 8 from Ref. [1]), for the 2 → 2 case the range of spin weights is limited to [0, 2] since in this process there are no couplings which could lead to individual τ polarization. In the 2 → 4 case the spin weight distribution exhibits a tail which extends beyond 2 and covers most of the allowed [0, 4] range. This is due to, e.g., the presence of the subprocess W + W − → X → τ + τ − in which W 's radiated off quarks are polarized which has impact on τ polarization. The tail above 2, although not so much pronounced, will manifest itself in the distribution of τ decay products. Let us now turn to the standard spin sensitive distribution of the ratio E π /E τ (a fraction of τ energy carried by the decay pion) used in [41] for benchmarking τ polarization. In every case discussed below we will use again X production amplitudes to calculate τ pair density mat-rix. We will do that also for the sample generated with H production amplitudes 8 .
The τ polarization can originate from the X production via VBF process, which is asymmetric over the phasespace regions due to the asymmetry of valence u and d quark distributions in the proton. To exhibit the polarization effects we have to sort out events according to the τ polarization; otherwise the effects will average out. Since in the proton there are more u-type quarks than d-type, the X particle produced in the VBF preferentially will follow the direction of W + which are right-handed and impart their polarization on X bosons. One can then expect that τ lepton from X decay will have polarization dependent on τ direction with respect to the X flight direction correlated with its spin polarization. Thus it is suggestive to sort events according to positive and negative value of where Y X denotes the τ lepton pair rapidity and p τ − z , p τ + z are the z components of τ ± four-momenta. In Fig. 6 events with positive and negative C are plotted separately (the first bin for C > 0 is lower exhibiting the pion mass m π /m τ effect). We observe that spin weights, calculated with the X production amplitude, when applied to the H sample lead to a larger spin effect, than when applied to the X sample. In the second case the spin effect is barely visible 9 .
Our results illustrate the complexity of multidimensional distributions. Even within tight selection there is a sizable difference between events of X and H production, which is reflected in τ polarization effects greater for the H sample than for X sample, even though the same pp → τ τ jj matrix elements featuring intermediate X are used in both cases.
One could argue, that such small spin effect present in Fig. 6 for the X case is a consequence of substantial contribution from other than VBF channel in our samples, thus pointing that our cuts may need to be refined. However, because of the weight distribution, as seen in the lower lower plot of Fig. 5, such a refinement is unlikely to be found within our tight selection, since the tail of events with spin weight exceeding 2 is very small. It seems that a better discriminating power between the H → τ τ and X → τ τ hypotheses can be provided by longitudinal τ -τ spin correlation, the same as discussed already in Refs. [3,41]. Nonetheless τ polarization may offer (minor) help in exclusion of X hypothesis, even in the case when Xτ τ couplings are insensitive to parity.

Summary and outlook
The main purpose of the paper was to demonstrate how the new (with respect to the ones used for sample generation) matrix elements for the production of τ -lepton pair accompanied with two jets in pp collisions can be used 8 Note that the separate treatment of τ production from the distributions of τ decay products enables evaluation of how important spin effects can be in experimental analyses. 9 It can be used nevertheless to improve exclusion of spin-2 hypothesis. in TauSpinner environment to reweight events. For that purpose, the New Physics matrix element for spin-2 X particle was implemented as a user example.
We have provided numerical tests of the algorithm, demonstrating that starting from the H sample (or X sample), the other one can be obtained by applying eventby-event weight calculated from the implemented matrix elements. We have also addressed possible technical difficulties and limitations in implementing the user code for matrix elements. Even though the TauSpinner algorithm in the case of its native and external matrix elements works similarly, technical aspects due to e.g. rounding errors and other numerical complications may differ; thus require individual attention. The density of events to be reweighted may differ from the target one significantly, resulting in a few events with weights massively larger than the ones from other regions of the phase-space.
Limitations of the algorithm, discussed in sub-section 4.3, may be observed if for example colour or spin config-urations for original and new process play important role for the parton shower. Then, reweighting with matrix element of hard process only Eq. (2), may be too simplistic and factorization properties may need to be addressed. An effort into that direction can be found in Refs. [42,43].
Let us stress, that the TauSpinner reweighting can be repeated several times on the same event to obtain multiple variants of weights, e.g. due to several variants of coupling constants, or even completely distinct X interaction forms. In our examples we have used rather small sets of non-zero couplings, see Appendix B, in part to simplify differences in distributions of X and H mediated processes. The reweighting algorithm performed better when reduced region of the phase-space was used for comparisons.
To demonstrate effects sensitive to the τ lepton polarization we have chosen τ ± → π ± ν decay mode as a spin analyser. Spin effects originate from the X production vertex and are embedded in complexity of the multibody phase-space. They turn out to be rather small for our choice of the Xτ τ couplings. Nevertheless, they may turn useful in falsifying physics hypotheses alternative to Higgs production and decay processes.
This paper completes the description of TauSpinner functionality, initiated in [3] for the 2 → 2 matrix elements of New Physics, now also with the vector-bosonfusion 2 → 4 matrix elements. It supplements examples of TauSpinner applications for events with two jets accompanying τ -lepton pair production in pp collisions, discussed in [7].

Acknowledgments
We thank Tomasz Przedziński for fruitful discussions and help with programming aspects of the TauSpinner  A Example with X mediated processes: user installation prototype The purpose of this Appendix is to present how the reweighting with matrix elements of X mediated process, as available in the program distribution tar-ball, can be used. It is equally important to demonstrate how any other external matrix elements prepared by the user can be installed: the X case can serve as a prototype. The detailed instructions how the reweighting algorithm works and how it can be used for final states of τ lepton pair and additional two jets is given already in Appendix A2 of Ref. [7]. In that paper, the non Standard Model matrix elements were not discussed.
In the case of our example of spin-2 X mediated matrix elements, the command TauSpinner::set_-vbfdistrModif(SPIN2::spin2distr); is used to set the pointer to SPIN2::spin2distr(...,KEY). This can be done for the user own matrix element routines. These routines should over-load the prototype ones for the non-standard calculation.
At initialization, when command spin2init_() is executed the masses and coupling constants for SPIN2::spin2distr(...,KEY) calculation are set. Later, for every event the algorithm makes the choice for the actual matrix elements used in the weight calculation: it evaluates and pass to the user provided function its internal parameter KEY. The KEY = 0, 2 corresponds to Drell-Yan like processes of the Standard Model and anomalous (user provided) matrix element. Analogously KEY = 1, 3 is for the Higgs of the Standard Model and (user provided) matrix element for anomalous narrow resonance. The code will choose between Higgs and Drell-Yan background amplitudes on the basis of PDG identifier of the intermediate resonance found in the event record 10 . For X.pdgID = 25 it will set KEY=0 for Standard Model (A) -that is for denominator of wt A→B prod given in Eq.(2) and KEY=2 for (B), weight numerator. For X.pdgID = 25 it will set KEY=1 for (A) and KEY=3 for (B). This is why user-provided function for the matrix element calculation must have a KEY parameter among its arguments.
The interface assumes that the event sample is for the SM i.e. as of type (A) used in weight denominator while the user provided function is of type (B), and accordingly the weight wt A→B prod is calculated. If the analyzed event sample is (as for Fig. 4) of type (B), then inverse of the weight given by Eq. (2) and calculated by TauSpinner should be used for reweighting.
So far, we have not discussed spin polarization and spin correlations between outgoing τ leptons. If a generated sample features spin correlations and polarizations of τleptons, then it has to be taken into account. The event weight needs to be supplemented with spin weights of the Standard Model and New Physics calculation where WT and WT0 denote usual spin weight as of Eq. (5) from Ref. [1] calculated first for (B) and then for (A). They can be obtained as shown in the extract from the code of the demonstration program.
The settings and parameters for the MadGraph5 [15] generation of H and spin-2 X event samples used in Section 4 are collected in the web page [16]. We remind the reader that these two generation parameter sets must be carefully matched, with the corresponding ones in the initialization of native TauSpinner (file ../TauSpinner/src/VBF/ VBF_init.f) and external matrix elements. In our example, it is the file ../TauSpinner/examples/ example-VBF/SPIN2/ME/SPIN2_init.f . A similar input consistency check will have to be performed whenever user own New Physics matrix elements are installed. Let us recall some details how the code for matrix elements of X mediated processes were prepared. After the default (automatic) initialization with MadGraph5 settings set group_subprocesses Auto set ignore_six_quark_processes False set loop_optimized_output True set low_mem_multicore_nlo_generation False set loop_color_flows False set gauge unitary set complex_mass_scheme False set max_npoint_for_channel 0 import model sm define p = g u c d s u~c~d~sd efine j = g u c d s u~c~d~sd efine l+ = e+ mu+ define l-= e-mudefine vl = ve vm vt define vl~= ve~vm~vtw e produced the event generation code with the sets of commands given in Table 3. The corresponding spin-2 param_card.dat file contains an additional b lock called spin2 for model specific couplings 11 . Block spin2 4 1.000000e+00 # gXtautauM 5 1.000000e+00 # gXtautauP 6 0.000000e+00 # gXqqM 7 0.000000e+00 # gXqqP 8 0.000000e+00 # gXggEven 9 0.000000e+00 # gXggOdd 10 1.000000e+00 # gXWWEven 11 0.000000e+00 # gXWWOdd 12 0.000000e+00 # gXBBEven 13 0.000000e+00 # gXBBOdd as well as a block containing information about the quantum numbers [45] of the spin 2 object.
Block QNUMBERS 5000002 # x 1 0 # 3 times electric charge 2 5 # number of spin states (2S+1) 3 1 # colour rep (1: singlet, 3: triplet, 8: octet) 4 0 # Particle/Antiparticle distinction (0=own anti) 11 Of the couplings of X to gauge bosons only the CP-even ones were implemented in the SPIN2 library. They are denoted in the code with a with suffix Even and correspond to couplings gXii, i = B, W, g, in Eq. (3).  In addition, the block containing Wolfenstein parametrization of the CKM matrix is called ckmblock in the Spin-2 param_card.dat (as opposed to name wolfenstein in the SM one).
The run_card.dat files used subsequently for generation of both H and X events are available from [16].

B.1 Setting precision parameter of MadGraph initialization.
The weights calculated with external new matrix elements may produce a tail of rare, high weight events. This may lead to excessively large statistical errors, sometimes poorly monitored in the histograms. Their origin is from sparsely populated phase space regions in which some events receive high weight due to some resonance/collinear configuration of the new process, see. e.g. discussion of Figs. 3 and 4 in Section 4.2 of Ref. [7].
Let us now comment on seemingly similar observations but of a fundamentally different origin. In the early steps of tests, for reweighting of X to H samples, we have encountered a mismatch between reweighted and Higgs X import model sm-ckm import model spin2_w_CKM_UFO generate p p > j j h QED=99, generate p p > j j x QCD=0 QED=2 NPVV=1 h > ta+ ta-NPll=1, x > ta+ taoutput h_dir output spin2_dir Table 3. Commands to generate with the help of MadGraph5 the source code for H and X amplitudes calculation.
reference distributions generated in the gridpack mode. It turned out that the problem was related to the accuracy of samples generated by Madgraph5 based on MadEvent [46] algorithm. At MadEvent core lies the diagram enhancement method which separates the integration into a sum of integrals whose singularity structure is dictated by a single Feynman diagram topology. Each of these individual integrals, referred to as integration channel, is then integrated using an appropriate phase-space parametrisation for its underlying structure. A MadEvent run involves two steps. The first one is referred to as the survey and consists of computing crosssections for each integration channel down to a given accuracy. Calling the generate − events script invokes a series of four commands: survey, combine − events, store − events, create − gridpack, in which by default survey should reach precision of 0.01 within 8 iterations, where the first iteration consists of 2000 points. However, with this default iteration number the required accuracy of 0.01 was not achieved and the kinematical distributions for X samples did not show a good behavior. After invoking the four commands one by one with the number of iterations in survey 12 increased to 99, the reweighted distributions were matching the reference ones (as is the case of Figs. 3 and 4). One can see in Fig. 7, that the requirement for precision for chosen distributions η jet2 and θ pj turned out to be particularly demanding.