Monte Carlo, fitting and Machine Learning for Tau leptons

Status of tau lepton decay Monte Carlo generator TAUOLA, and its main recent applications are reviewed. It is underlined, that in recent efforts on development of new hadronic currents, the multi-dimensional nature of distributions of the experimental data must be taken with a great care. Studies for H to tau tau; tau to hadrons indeed demonstrate that multi-dimensional nature of distributions is important and available for evaluation of observables where tau leptons are used to constrain experimental data. For that part of the presentation, use of the TAUOLA program for phenomenology of H and Z decays at LHC is discussed, in particular in the context of the Higgs boson parity measurements with the use of Machine Learning techniques. Some additions, relevant for QED lepton pair emission and electroweak corrections are mentioned as well.


Introduction
It is thirty years since the first versions of TAUOLA package [1][2][3][4] for simulation of τ -lepton decays and PHOTOS [5][6][7] for simulation of QED radiative corrections in decays became public. The interfaces of these programs are designed in a way to fulfil needs of different groups of users. The bulk of the code remain written and maintained by the same (main) authors, but contributions originating from other researchers became integrated over the time too. These programs became parts of a wide range of applications, sometimes consisting of quite complicated simulation chains. Some versions of these codes at the moment require an independent maintenance. Some of those code modifications are not prepared and maintained by the program authors. This is often the case of TAUOLA. From the design point of view new versions differ little, usually in numerical values of some constants or analytical form of hadronic currents. Nonetheless, these changes are of a great value because they often depend on fits to experimental data, sometimes even not fully made public by the experiments. Variants may archive details of some other τ -lepton phenomenological projects, these topics were presented in previous τ conferences.
Our presentation is organized as follows: Section 2 is devoted to the discussion of recent changes in TAUOLA package. This point was already discussed in [8] and in [9], that is why, we will be brief. In Section 3 we summarize recent developments for PHOTOS Monte Carlo of radiative corrections in decays. Section 4 is devoted to applications of TAUOLA for hard processes with final state τ leptons. In particular, for construction of spin observables and evaluation of their sensitivity. The modern techniques like Machine Learning (ML), were found to be important in analyses of HEP data. In this context we discuss TauSpinner algorithm, which was found to be useful for evaluation of observable for Higgs boson parity measurement. We mention other applications or tests for this tool; in particular in the domain of algorithm of calculating spin states of τ pairs in events where high p T jets are present in pp collisions. The following two sections provide examples for applications in the domain of precision measurements and Higgs boson CP parity evaluation at LHC. Summary Section 7 closes the presentation.

Currents and structure of TAUOLA Monte Carlo
The TAUOLA and PHOTOS represent long term projects, that is why no substantial changes were introduced since the τ conferences of 2014 and 2016 [8,10]. Nowadays, the C++ became dominant in many segments of the code for programs and for the tests as well. The structure defining part of the code was modified, to prepare for more modular code organization. In fact only small changes in the text of the code were necessary, that is why, numerical tests if introduced changes were bug free were simple. Now, once decision to re-write the whole code into C++ is finally taken, it can be performed in quick and well controlled steps, each to be tested separately. See Ref. [9] for more explanations. Let us point, that language change for TAUOLA is more complex than it was for PHOTOS [11]. This is because of the three phenomenology aspects of the work which need to be taken into account: • Development of physic assumptions and later of the code, for new versions of hadronic currents of all or some τ decay channels.
• Preparation of experimental data for fits. Question of background control or dimensionality of the distributions is of a great importance. For the optimal performance participation of physicists involved in experiments data analysis is important.
• Preparations of algorithms and choice of distributions to be used in fits. We will return to this point later in the presentation.
Synchronization of work on these tasks and evolution of the code need to be assured. That is not easy, as different researchers are involved in each of these activities.
In the following sections, Optimal Variables and ML techniques are evaluated for purposes of Higgs and Z phenomenology. We expect it of importance for modeling of τ decays as well. Especially in this case, proper evaluation of multidimensional distributions used for fits and in particular systematic uncertainties may be of a great importance. We study if Monte Carlo matrix element dependent event weights may be useful.
Finally let us recall some limitations on physics precision of the TAUOLA resulting from the program design. Up to a precision level of about 0.2% choice of hadronic currents play the central role for systematic uncertainties. That is why, corresponding parts of the program represent well defined, easy to replace, building block. Confrontation of model's predictions with experimental data with the help of multi-dimensional distribution has to be central for the future project developments. This was pointed already in Ref. [12] and its importance is clear from our experience of work on hadronic currents for τ → 3πν decay modes as well [8,10].

PHOTOS Monte Carlo for bremsstrahlung: its systematic uncertainties
Over the last two years no major upgrades for functionalities were introduced into PHOTOS Monte Carlo, except introduction of emission of lepton pairs. Published documentation of the program [11], correspond to the present day head version [13]. Numerical tests for pair emission algorithm are published [14], and updates on the program necessary due to evolution how event records contents have to be understood are continuously updated, see changelog.txt of Ref. [13]. A possibility for the user modification of the Z → l + l − γ matrix element is prepared. The appropriate matrix element can be replaced by the user own. In this way, e.g. effects of non-standard-model matrix elements or leading contributions of the loop electroweak corrections, can be studied.

TAUOLA -hard process -TauSpinner algorithm
The packages such as TAUOLA or PHOTOS are rarely used alone. Also tests of the programs can not be performed fully independently of users' projects. Thus stand-alone tests are not sufficient. Typically, user applications rely on other libraries of programs as well, which in size may surpass largely TAUOLA or PHOTOS themselves. The complete simulation chains consist not only of segments for so called 'truth' i.e. physics processes based on theoretical predictions, but emulate detector response as well. Physics interest is usually not focused on properties of simulation chains, but on intermediate state properties such as mass, coupling constants or parity. If intermediate state properties can be modified by the experimental user and consequences of such changes for detector responses as well, work for defining sensitive observables can be simplified.
For the final states with τ leptons, methods to manipulate their spin state can be used to optimize the measurements. The τ is the only lepton of its spin accessible to the measurements. From the first paper [15] TauSpinner was oriented for such applications. The response due to changes in Z, W or H decays, represents a valuable information. Let us review the status. General structure of the tool was recently documented in [16], but let us nonetheless point to the general idea. The program is calculating weights corresponding to changes of the physics assumption. Ratios of matrix elements squared for the compared options is used for that purpose. As an input, events stored on the data file are used. No changes for kinematical configurations are introduced. In that respect nothing has changed since 2016 τ conference and remain as in Ref. [10], however new applications were documented in Refs. [17,18]. With these publications, possibility to use matrix elements, for parton level processes with the two outgoing jets accompanying τ pair, was introduced. Another new option was to introduce, with the help of TauSpinner, weights that account for the electroweak loop corrections to Drell-Yan processes [19].
With the help of TauSpinner we could evaluate observable to study Higgs boson parity, in its cascade decay with intermediate τ leptons [20]. More precisely; we have investigated the potential for measuring the CP state of the Higgs boson in H → τ τ decay, with consecutive τ -lepton decays in the channels: τ ± → ρ ± ν τ and τ ± → a ± 1 ν τ combined. Subsequent decays ρ ± → π ± π 0 , a ± 1 → ρ 0 π ± and ρ 0 → π + π − were taken into account. We have extended method of Ref. [21] first. Also in the case of the cascade decays of τ → a 1 ν, information on the CP state of Higgs boson can be extracted from the acoplanarity angles. However for a ± 1 → ρ 0 π ± , ρ 0 → π + π − cascade, up to four planes can be defined, thus 16 distinct acoplanarity angle distributions are available for H → τ τ → a + 1 a − 1 νν. Each of such distributions carry some information but it is laborious to evaluate overall sensitivity using previous approach.
We have investigated the sensitivity with the help of Machine Learning (ML) techniques. Since the 2016 τ conference, we have in part taken into account ambiguities resulting from the detector uncertainties and background contamination, see Ref. [22]. Sensitivity largely survived.
Usefulness of ML methods and τ → 3πν decays for Higgs boson parity measurement can be understood as an example. It required analysis of multi-dimensional distributions simultaneously for signatures and background processes. The case of Higgs boson parity remained physics-wise simple, the signatures to be distinguished, were described by analytically clear and simple weights. The examples which we describe below, of interest in themselves, may also lead to helpful solutions for physics of τ decays. There, of course theoretical description is challenging, but thanks to better control of interferences (thanks to multidimensional distributions) new insight may become available. Let us now return to the examples of the high energy domain. Here, the ML and more classical (but leading to intuition) solutions based on so called Optimal Variables techniques were used.

Towards Optimal Variables for Z/H → τ τ spin observables
The knowledge of the τ lepton kinematic is essentially important for the spin analysis. An angle between the τ lepton flight direction and a neutrino in its decay is a powerful spin analyzer. This angle can be reconstructed unambiguously only if both the total momentum and the flight direction of the τ lepton are reconstructed. Unlike in e + e − collisions in pp collisions at LHC there is no beam energy constraint, however it is possible to place a reasonably good estimate on the neutrino momenta in the decays Z/H → τ τ [23][24][25]. In this section two examples of the TauSpinner application to the analysis of the longitudinal and transverse spin of τ leptons in the decays Z/H → τ τ assuming that the whole kinematic of the decay is available are discussed.

Longitudinal τ polarization in the decay Z → τ τ
Let us start with the example of longitudinal τ polarization, the observable important for the precision tests of the Standard Model. The difference of the neutral weak couplings to the right-and left-handed fermions results in the polarization of the fermion-antifermion pairs produced in the decay of the Z bosons. A measurement of the τ leptons polarization at LHC provides an independent and complementary determination of the effective weak mixing angle sin 2 θ τ ef f as well as a test of the lepton universality of the weak neutral current. The measurement of the τ polarization requires a knowledge of the τ spin state, this can be concluded by analyzing the angular distributions of the τ decay products with respect to the τ flight direction. Following the notations of [3] the partial width of any τ decay is given by: where M is the τ mass, |M | 2 spin averaged matrix element, s -four-vector of the τ polarization and polarimetric vector h is a function of the τ momentum and the decay products. The standard abbreviation for Lorentz invariant phase-space integration element dLips is used. And angle cos θ h (further referred as ω h ) between the polarimetric vector h and the τ flight direction, as seen from the τ rest frame, carries the full information about the spin of the τ (assuming only longitudinal spin component of the τ leptons). An explicit form of the polarimetric vector for τ decays, can be well identified by CMS [26] or ATLAS [27] detectors (τ − → π − ν, π − π 0 ν and π − π − π + ν), it can be found in [3]. An analyzing power of the τ polarization can be further gained noting that the helicity states of both τ leptons in the decay Z → τ τ are almost 100% anti-correlated. Denoting ω 1 h and ω 2 h to be the observables for both τ leptons a combined observable is given by [28]: The distribution of a single ω h and Ω hh for the decays τ − → π − ν, π − π 0 ν and π − π − π + ν is shown in Fig. 1. It should be noted that only in the decay τ → π − π − π + ν there is a model dependence that comes from the imperfect knowledge of the a − 1 → π − π − π + decay structure. It has been shown that assuming the CLEO parametrization of the a − 1 → π − π − π + decay [29] the modeling systematic uncertainty on the the τ polarization measurement is negligibly small and will not limit the precision [30].

Transverse spin correlation in the decay H → τ τ and a quest for an Optimal Observable of the Higgs boson CP parity
There is no longitudinal net polarization of τ leptons in the decay H → τ τ , however the transverse spin correlation of τ leptons in the decay H → τ τ might reveal the information on the CP structure of this decay.
Similarly to the longitudinal spin analysis the analyzing power of the transverse spin correlation can be maximized by considering the full kinematic of the τ pair. Denoting h 1 and h 2 as the polarimetric vectors of both τ leptons in their rest frames one can introduce two observables (the so-called acollinearity and acoplanarity angles): where k 1,2 = h 1,2 × n 1,2 | h 1,2 × n 1,2 | , and n 1,2 are the unit vectors pointing along the direction of the first and second τ lepton in the H rest frame. The distribution of δ and φ for even and odd Higgs boson CP state are shown in Fig. 2. The distributions are uniform for all considered τ decays and carry the full analyzing power. However, the analyzing power in the experiment will be naturally diluted by the detector effects and limited resolution. The best performance can be expected when both τ leptons decay into three charged pions. Three charged tracks in the detector offer robust reconstruction of the point of each τ decay (no flight length of the quickly decaying a 1 resonance) and the high multiplicity of pp collisions allows to determine the point of interaction. This information, imposed as an additional constraint, can significantly improve the performance of algorithms for kinematic reconstruction of the τ leptons momenta.
In this approach all experimental complexity is hidden in the measurement of h 1 and h 2 . The δ or φ take the role of Optimal Variables. That is very helpful for physics intuition.

From Optimal Variables to Machine Learning approach
Let us quote an example of numerical result taken from Ref. [20] and Ref. [22]. For the analysis of the sensitivity of H → τ τ ; τ → ν τ 2(3)π we have used ML technique. This path was explored, because form the previous studies [21] we could expect that the signature will depend of the distribution features embedded over multidimensional space. As a part of our approach was to evaluate reliability of the new for us techniques, we performed study of prototypes for the possible elements of construction of Optimal Variables, defining hyperspace of this multidimensional space. We recall results, the Table 1 for the case when detector smearings were not taken into account, and Table 2 where they were taken into account. The encouraging results have been obtained.
Let us now recall a few details, that may be useful if an attempt to apply similar techniques will be directed, e.g. toward evaluation of models and fits of the τ decay matrix elements. The discussion of the systematic uncertainties of the results could be simplified, because for each event we could calculate matrix elements for the confronted assumption. Complexity of the resulting weights was the results of τ decays, but the part describing vertex H → τ τ was rather simple and of clear character.
In case of τ decays complexity and ambiguities in definition of matrix elements as a function of details of the intermediate energy strong interaction models will be inevitably larger. However, even in this case corresponding to the confronted assumption weights will be available. Already now, this is technically prepared and embedded in TauSpinner algorithm. Note that different assumption of τ decay models may lead to relatively similar predictions for one dimensional distribution, but largely differ for the multidimensional ones, see for example [31].

Features/variables
Decay mode: ρ ± − ρ ∓ Decay mode: a ± 1 − ρ ∓ Decay mode:  Table 1: Average probability p i (calculated as explained in Ref. [20]) that a model predicts correctly event x i to be of a type A (scalar), with training being performed for separation between type A and B (pseudo-scalar). Different sets of variables were used as an input. The ϕ * i,k , y i , y k can be understood as approaches to construct expert variables, but as we could observe, their use did not brought improvements with respect to the case when 4-vectors of the τ decay products were used alone, provided proper reference frame was used. Results taken from Ref. [20].
In our applications for the Higgs boson we have confronted two assumptions for the Higgs boson CP parity states, in most cases it was choice between scalar and pseudo-scalar variants only. However, the forthcoming solutions based on ML approach will have the possibility of evaluating best performing parameters of some models too. This will open the way for fits of τ decays themselves. We leave such solutions for the forthcoming works.

Summary and future possibilities
There was not much of the development for the implementation of τ lepton decays embedded in TAUOLA library. The new version of the program implementation was archived and finally published [9]. The status of associated projects: TAUOLA universal interface and TauSpinner was reviewed. Also new results for the high-precision version of PHOTOS for QED radiative corrections in decays, were presented, in particular algorithm for emission of additional light lepton pairs was supported with documented tests.
Some details of presentation of the TAUOLA general-purpose C++ interface was given and its applications were shown. Use of weighted events was demonstrated to be useful for studies of LHC phenomenology in domain of τ lepton polarization observables for precision tests of Standard Model as well as for Higgs boson CP parity. Complementary methods of Optimal Variables as well as of ML approaches, both useful for phenomenological applications were presented and example results were shown.

Acknowledgments
The work on TAUOLA would not be possible without continuous help an encouragements from experimental colleagues. The work of all co-authors of the papers devoted to TAUOLA development was of great importance. I hope, that it is clearly visible from my contribution.  Table 2: Area Under the Curve is shown for the Neural Networks (NN) studied to separate scalar and pseudo-scalar hypotheses. Inputs with a are used. Results in column "Ideal" -from NNs trained/used with particle-level simulation, in column "Smeared" -from NNs trained/used with smearing. NN trained on smeared samples when used on exact samples give similar results as "Ideal" what in this case mean no detector smearing. Result taken from Ref. [22]. Reference frame as for Table 1.