Study of variants for Monte Carlo generators of τ→3πν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau \rightarrow 3\pi \nu $$\end{document} decays

Low energy QCD (below 2 GeV) is a region of resonance dynamics, sometimes lacking a satisfactory description as compared to the precision of available experimental data. Hadronic τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} decays offer a probe for such an energy regime. In general, the predictions for decays are model dependent, with parameters fitted to experimental results. The parameterizations differ by the amount of assumptions and theoretical requirements taken into account. Both model distributions and acquired data samples used for the fits are the results of a complex effort. In this paper, we investigate the main parameterizations of τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} decay matrix elements for the one- and three-prong channels of three-pion τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} decays. The differences in analytical forms of the currents and resulting distributions used for comparison with the experimental data are studied. We use invariant mass spectra of all possible pion pairs and the whole three-pion system. Also three-dimensional histograms spanned over all distinct squared invariant masses are used to represent the results of models and experimental data. We present distributions from TAUOLA Monte Carlo generation and a semi-analytical calculation. These are necessary steps in the development for fitting in an as model-independent way as possible, and to explore multi-million event experimental data samples. This includes the response of distributions to model variants, and/or numerical values of the parameters. The interference effects of the currents’ parts are also studied. For technical purposes, weighted events are introduced. Even though we focus on 3πντ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3\pi \nu _\tau $$\end{document} modes, technical aspects of our study are relevant for all τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} decay modes into three hadrons.


Introduction
One aspect of phenomenological work is to represent the theory and experimental efforts in the form of Monte Carlo generators. Then the user should be able to obtain any disa e-mail: jakub.zaremba@ifj.edu.pl tribution of interest. All details of the detector response can be taken into account if needed. An introduction to the tasks of modeling with the help of Monte Carlo of low energy hadronic interactions can be found, e.g., in [1] and will not be repeated here. We will concentrate on τ lepton decays.
Let us point out that there are various aspects of the work for construction of Monte Carlo, which need to be combined. On the theoretical side, assumptions have to be translated into distributions which can be confronted with the data. This task is easy if the underlying theory is well established and its calculation methods, based, e.g., on perturbation expansions, are converging sufficiently fast. For theorists, it would be ideal if the experimental results were represented in the form of background-subtracted and detector corrected distributions, sufficient to constrain all constants (masses, widths, couplings) introduced by the theory (or model) used as a basis for the calculation.
In practice, such a task is often far from being straightforward. The assumptions behind theoretical calculations are often not well established and the number of phenomenological constants introduced may be too large to be constrained by experimental data.
From the point of view of Monte Carlo techniques, generating processes such as τ decays into 3πν is rather simple. However, it being necessary for that purpose, hadronic currents span distributions over eight-dimensional space and are the result of a massive effort. Taking into account the Lorentz invariance and properties of weak couplings of the τ -lepton to an intermediate virtual W boson, the predictions still require four complex scalar functions F i of three variables each (see the definition later in the text). Such a parameterization of matrix elements is commonly used by τ decay Monte Carlo generators, e.g., by TAUOLA; and it exists since its beginning [2,3].
In principle, all properties of the matrix element can be constrained from the experimental data. All F i can be fitted in a model-independent way as proposed in [4]. In practice, this is highly non-trivial. Even if necessary for that purpose, the direction of theτ neutrino can be reconstructed; one has to measure from the data at least seven distributions over three dimensions with sufficient detail. So far, this was never achieved without partial integration over some of the phasespace dimensions. This practical aspect of the phenomenology work should be separated from the previous two (theoretical and experimental ones), and at the same time offer comfortable and flexible usage.
In this paper, we investigate the differences and physics motivations for the presently used in TAUOLA hadronic currents of the τ → 3πν τ decay channels. Let us recall from Ref. [5] the comparisons of different parameterizations for the case of the decay mode τ → 2πν τ . In that case, all results are encapsulated in a single one-dimensional distribution. Our analysis, as we will see, must rely at least on a three-dimensional distribution. We compare its efficiency to constrain the model parameters first with the method when one-dimensional histograms are used only. To evaluate the full size of the model differences, an integral of |wt − 1| = |M model | 2 |M model | 2 − 1 over the event samples will be used. M model , M model denote the matrix elements for the two compared model variants. For completeness, we study the case when only the total width of the channel is used to constrain the model.
Theoretical principles embedded in the models and the remaining freedom for introducing new couplings and/or intermediate resonances has to be studied at the same time. We have to keep in mind that additional contributions to the currents can modify the shapes of the distributions in a nonphysical/uncontrolled way and in unexpected regions of the phase-space. This is especially important when the dimensionality of the experimental distributions is smaller than that of the decay phase-space.
The paper is organized as follows. In Sect. 2, we present the list of currents used for evaluation of the distributions. A detailed description is delegated to Sects. 3-6; the analytic forms of the currents are given. In Sect. 7, we present basic features and motivations behind the models used in data fitting. Section 8 consists of an introduction to numerical comparisons and commonly used representations of experimental data employed in the fits. Section 9 is devoted to models used by the CLEO collaboration. In the subsections, numerical results and comparisons are presented. Section 10 describes further variants of the currents that were never used outside the experimental collaboration, but their existence is mentioned in Ref. [6]. Two main CLEO models are then compared with the BaBar model in Sect. 11. Section 12 provides comparisons between the default TAUOLA model and two models based on the resonance chiral Lagrangian (RChL) [7]. For this task, we use mainly one-dimensional histograms; results of measurements in such a form are available not only from the CLEO but from the BaBar collabo-ration as well. Such distributions were used to constrain the model parameters in many of our studies. Only CLEO used the three-dimensional distributions as an input for fits.
In Sect. 13, contributions from different parts of the currents and their interferences are investigated. We discuss possible consequences of using model assumptions and limited data samples in Sect. 14. A summary, Sect. 15, closes the paper. In "Appendix A" we present arrangements for semi-analytical calculations which can be used for fits into three-dimensional histograms. In "Appendix B" numerical comparison of the results from different models which are of lesser importance are given. For completeness some histograms of three-dimensional form are given in "Appendix C".

List of currents
In the TAUOLA Monte Carlo generator [8,9] multiple versions of hadronic currents for the two channels of the τ → 3πν τ decay are available. They are the result of theorists' efforts but also of extensive work of the experimental collaborations. Over time, currents became available for use outside of the experiments. Unfortunately, sometimes these currents are poorly documented. We plan to explicitly present differences in formulas and values of numerical constants for the main currents available for TAUOLA. For others, we will skip some details as they are of lesser importance or are documented in the quoted publications.
Code of this current was obtained from the CLEO collaboration and distributed with TAUOLA (all its versions since Ref. [8]). See Sects. 4 and 9 for details. 3. TAUOLA CLEO isospin intricate-this is a variant where isospin rotation from the π 0 π 0 π − to the π − π − π + channel was performed for intermediate resonances.
Numerical constants remained unmodified. The code was obtained from the CLEO collaboration and distributed with TAUOLA (all versions since [8]) but it was never active. 1 See Sects. 4 and 9 for details. 4. TAUOLA BaBar-a parameterization numerically equivalent to the one used in the BaBar collaboration for basic simulations. Details are given in Sects. 5 and 11. 5. TAUOLA RChL 2012-a version of the RChL model introduced in TAUOLA in 2012. This model was based on the theoretical consideration of that time and fits to invariant mass distributions of π − π − π + and π − π + systems [13]. For details see Sects. 6 and 12. 6. TAUOLA RChL-this is the model motivated by further theoretical consideration and by comparison with experimental data also for the π − π − invariant mass; see Ref. [14]. 7. TAUOLA CPC-outdated parameterization of the Monte Carlo generator. Foundation of many technical benchmarks for phase-space and other verifications which are documented in [3]. Also, test distributions documented in Ref. [4] are available. 8. Pythia CLEO-the parameterization of Ref. [6] as implemented in Pythia [15]. A current nearly identical to the TAUOLA CLEO isospin intricate is used. Currents of Pythia differ only by numerical constants. In the present paper, we use numerical results of Pythia version 8.201. This current will be discussed in "Appendix B".
At present, TAUOLA CLEO is the default choice which should give the same result regardless of which version of TAUOLA is used. It remained unmodified for at least 14 years. It is identical for the FORTRAN and C++ implementations. The differences in distributions between CLEO publ., TAUOLA CLEO, TAUOLA CLEO isospin intricate, and Pythia CLEO are presented with all technical and numerical aspects. For completeness, we present numerical results from TAUOLA RChL. This last current is aimed to represent data with good precision and at the same time improve the relation with present day theoretical calculations. Numerical results for the other variants of the currents are given in less detail or are not presented at all.
where T μν = g μν − Q μ Q ν /Q 2 denotes the transverse projector, and Q μ = ( p 1 + p 2 + p 3 ) μ is the momentum of the hadronic system. 3 The decay products are ordered π − π − π + for the three-prong channel and π 0 π 0 π − for the one-prong channel and the pions' four-momenta 4 are denoted p 1 , p 2 , and p 3 , respectively. The same ordering is used in the further equations also for masses (m i ). The μ . νρσ is the Levi-Civita symbol. 5 In the equations of this and the following sections we use the notation: The constants c 1 , c 2 , etc. are Clebsch-Gordan coefficients, defined specifically for the particular hadronic current used.
Equations (2)-(4) describe the Breit-Wigner functions that are later used in the definition of the form factors: These are typical building blocks useful for hadronic currents parameterizations, such as of the Gounaris-Sakurai parametrization [16] for ρ → ππ.

Hadronic current in models of the CLEO category
Let us describe the default, hadronic currents of the CLEO modeling-TAUOLA CLEO. In Eqs. (5) and (6) the analytic form of the hadronic current for the τ → 2π 0 π − ν τ decay channel is presented. 6 Its functional form is extracted from TAUOLA code 7 but has a form exactly as described in "Appendix A" of Ref. [6]. The notations differ slightly. This form of the hadronic current is also used for the τ → 2π − π + ν τ decay channel in the default TAUOLA setup. It will be discussed later why such a choice is justified. In the PDG tables, Ref. [17], some mesons have the same name but different masses. In our paper a 1 (1260), σ [or f 0 (500)], f 2 (1270), f 0 (1370), ρ(770), ρ (1450), and K * (892) are used. Finally, 3 In TAUOLA code defined as PAA(4). 4 In TAUOLA code represented as PIM1(4), PIM2(4), PIM3(4). 5 The μ . νρσ p ν 1 p ρ 2 p σ 3 is coded in the subroutine PROD5. 6 The constants β, β 1 , β 2 etc. are coded in TAUOLA as BET, BT1, BT2, etc. 7 Such an approach checks if over the years any modifications were (intentionally or not) introduced.
we have the meson a 1 (1700) as in Ref. [6]. The masses in the names will be skipped from this point on. We have F 2 has the same functional form as F 1 . The only difference is an interchange for its arguments indices 1 and 2 in Eq. (1), and that the constant c 2 has a sign opposite to c 1 . Note that β is always set to 0. We suspect that it was introduced by the CLEO collaboration for studies of the a 1 influence. We have The form factors F 4 and F 5 are set to zero. 8 Note that the normalization factor used in the a 1 and a 1 widths is the same. The Q 2 dependence of the a 1 width in its propagator is given by Eq. (7). 9 It was obtained from an analysis performed by the CLEO collaboration [6]. The contributions from the three main a 1 decay channels are taken into account. A complicated form is determined by the decay channels: as the a 1 virtuality gets larger, thresholds are crossed, allowing more decay channels to open, therefore changing the Q 2 dependence of the effective width. Q 2 is given in GeV 2 units. In the equation below, C 3π = 0.2384 2 and C K * = 4.7621 2 C 3π . We have Let us now turn to the case of TAUOLA CLEO isospin intricate. It differs from the previous for the π − π − π + only. Equations (5) and (6) are still used for π 0 π 0 π − , but they are replaced by (8) and (9) in the three-prong channel. Those are also extracted from the TAUOLA code, but from a part which was not available for general use. We have Again, F 2 has the same functional form as F 1 . The only difference is an interchange of indices 1 and 2 for its arguments in Eq. (1), and that constant c 2 has a sign opposite to c 1 . We have One can ask why the isospin symmetry was imposed taking into account the final states (which means distributions only) in the version of the code used as default in TAUOLA for so many years, even though it is not the property of the model used; see e.g. Ref. [18]. We will return to this point later, when discussing the numerical results. The current of Eqs. (5) and (6) is used in TAUOLA for one-prong and three-prong of the three-pion channels in most of the applications. 10 Note that neither Ref. [11] nor [10] gave the functional form of the τ → 2π − π + ν τ current; also Ref. [6] presents the functional form of the hadronic current for τ → 2π 0 π − ν τ only. The parameters used for the currents will be presented later.

Analytic form of the hadronic current of TAUOLA BaBar
This model became public recently, it was developed by the experimental collaboration. As mentioned in Sect. 1, isospin symmetry and limitation of the quality of the data was probably behind the choice of form factors used by BaBar as well. They are the same for τ → 2π 0 π − ν τ and τ → 2π − π + ν τ . The BaBar model does not include resonances f 0 , f 2 , and σ . Note that the BaBar choice was introduced after the CLEO currents (both variants) had become known and available for them. The equations in this section are extracted from TAUOLA code, numerically equivalent to the one used by BaBar. 11 We have It has an opposite sign as well; F 3 = 0; F 4 = 0; F 5 = 0. We have Note that Eqs. (7) and (11) both aim at the a 1 phase-space dependence, but Eq. (11) does so in a simpler manner. The CLEO model uses a more sophisticated polynomial interpolation with inclusion of a K * K production threshold; above a 1 it can decay relatively abundantly into a pair of kaons and a pion.

Resonance chiral Lagrangian currents
Our first hadronic current based on the resonance chiral Lagrangian scheme, TAUOLA RChL 2012, was prepared for Ref. [20]. Similarly to the parameterization of BaBar, the same form of current was used for one-and three-prong τ decay modes into three pions. As an input for the fits, only two one-dimensional histograms of BaBar data were used. Later on, the invariant mass of π − π − pair became available. It turned out that contributions which, under isospin symmetry, provide different forms of the current for oneand three-prong decay modes were useful. We refer to this later parameterization as TAUOLA RChL; it is documented in Ref. [14].
We will return to the isospin symmetry aspect, when we will discuss the sensitivity of different distributions to specific parts of the current contributions. This may be a hint 11 The authors would like to acknowledge Swagato Banerjee and Tomasz Przedzinski for cross-validation of our modifications of the τ → 3πν hadronic current called TAUOLA BaBar, with the TAUOLA installation in the BaBar framework at a statistical level of 100 M events. For the comparisons, standard distributions prepared with MC-TESTER [19] were used. why isoscalar terms, introduced by CLEO many years ago, were actually dropped out in the currents, such as TAUOLA BaBar and TAUOLA RChL 2012.

Practical aspects of model construction
There are two main paths to be taken in model building. The empirical approach is to prepare something describing the data in the best possible way, without considering how physical processes work. The other choice is to start from theoretical principles, obtain distributions, and introduce only necessary adaptations on comparison with data. Of course, merging those two approaches is desirable and even essential. There are unquestioned rules, like Lorentz invariance, which have to be obeyed. All models used in TAUOLA and in our paper respect such rules.
In particle physics we expect Breit-Wigner distributions to describe the shapes of the experimental distributions. One should bear in mind that the precision of the experimental data of today can be very high. Multiple millions of events were collected for individual τ decay channels. On the other hand, once inspection of the assumptions behind theoretical models is performed, rather modest estimations for the systematic uncertainties of the predictions are obtained [21]. For example, for RChL models, in Ref. [20], Section 7.4, it was argued that on the basis of theoretical arguments alone, the predicted uncertainty can in principle be as high as 30 %. Even behind such a basic principle like isospin symmetry uncertainty at the level of 5-10 % can be expected at least in some cases; see, e.g., Ref. [1], Section 5.11. That is again more than one order of magnitude worse than the precision of the experimental data. Comparing to differences in the results of Sect. 9.2 it is of about the same order of magnitude.
That is why a comparison of the model's results with the data can be an exciting source of inspiration for further model developments, even if the attained agreement with the data is not perfect at the beginning. Such a development can be particularly promising for the models precisely representing distributions of several τ decay channels and/or final states in e + e − → hadrons.
In the case of TAUOLA CLEO, at the start, the model of Ref. [22] was used. With time, as a consequence of comparisons with the data, it evolved into a form where Breit-Wigner resonances seen in the experimental data were added. Some of the model assumptions, in particular numerical values of the chiral couplings in the soft limit, were compromised. One should keep in mind that only at very low energies the chiral dynamics is expected to determine unambiguously the hadron form factors probed by semileptonic τ decays [23][24][25].
At present, a more elegant approach is taken; the RChL model. It was started from theoretical calculations accounting the expansion in the number of colors behaving as N c → ∞ [26,27]. Later, this model required an empirical correction by the σ meson introduced in Ref. [14] as an additional Breit-Wigner contribution. The authors are still working out a way to include σ and possibly other isoscalars within the RChL scheme using unitarity and analyticity constraints to include this re-scattering effect.
Another interesting approach in modeling was taken, e.g., in Ref. [28]. It describes a hadronic current for τ → 4πν τ . This modeling is in part data driven. As measurements of e + e − annihilation into four pions arrived [29], they could be used to model the distribution of the invariant mass of the whole 4π system. 12 The authors introduced a correction into the matrix element to reproduce the experimental distribution. Then a 1 π dominance as the intermediate state was assumed. From that point on, any desired modeling of a 1 could be used. The one of Ref. [30] was chosen, as it was applied for the e + e − data as well.
Those are only a few of many possible approaches. There is no straightforward way of telling that one is better than the other. The actual choice has to be motivated by the quality and precision of the experimental data. Let us address this point in the next section.

Data representations and numerical methods
Experimentalists and theoreticians may prefer distinct forms for data representation in model building, comparisons, and fits. In this section we will concentrate on methods we use and those that were used for constraining investigated models in the past.
The simplest way to compare models and also to fit the experimental data is to look only at the total width of the channel 13 (we will call it option 1). Such an approach is used when there is no sufficient data available for a more advanced fit. Even though it is only one number, it can be used as a first test of any particular model. The total width is usually well documented by the PDG tables; see, e.g., Ref. [17], including a careful analysis of experimental statistical and systematical errors.
The next approach is to look at one-dimensional distributions. In this case, one should seek the possibility of having experimental data for all possible invariant mass distributions. This is not always possible. Using as an example the 12 It was done by relating the 2π + 2π − and π + π − 2π 0 systems produced in e + e − annihilation with π + 2π − π 0 and π − 3π 0 produced in τ decays. To this aim, isospin symmetry was applied to distributions measured from e + e − annihilation giving a prediction of the shapes of four-pion system invariant mass spectra in τ decays. 13 Option 1 was used, e.g., in Ref. [31] for 5π currents. π − π − π + channel, one must be aware that it is easy to have many models that fit perfectly, let us say the π − π − π + spectra, while having distributions of π − π − and π − π + off. We can articulate three options of one-dimensional distributions in use for the case of the τ → 3πν τ current input: π − π − π + system invariant mass only 14 (option 2).
In a more advanced approach of Dalitz distributions (option 5), slices in Q 2 were used already in Ref. [6]. This way of data representation is in nature three-dimensional. 17 Models originating from the CLEO collaboration were fitted in that way. To explore the whole structure of the decay channel, the methods of Ref. [4] have to be introduced (option 6). Up to now, no complete analysis was performed with their help. Only Ref. [6] mentions a cross-check with such methods.
Let us turn to practical aspects of our paper's comparisons for distributions with the options listed above. For that purpose we describe a re-weighting method. We use it often due to its numerical efficiency and simplicity. When generating events from a Monte Carlo generator for different models, statistically independent event samples are produced. But one may ask what the probability is of getting a particular event with a certain set of momenta depending on the model we use. Having all momenta of particles in an event, one can calculate the amplitude for such a specific configuration. Therefore, one can calculate the amplitudes for the event using each model and compare. This leads to the possibility of re-weighting events by using one model as a reference and instead of generating a new sample with a different model, calculating amplitudes from that other model for existing events. The ratio of the new amplitude squared to the old one can then be attributed to the event as its weight. Thanks to the use of weights, the samples are correlated. The statistical errors affect the difference for the compared models only. Using this strategy for a different model we may not only get the new distributions with a different shape, but also its integral, that is, the total width of the channel. The method can be applied to minuscule changes of the model parameters, which are needed, e.g., for the calculation of the derivatives. 14 Similar to option 2 a fit was performed for the 4π current in Ref. [28] with only the invariant mass of four pions used. 15 Option 3 was used in fits of TAUOLA RChL 2012. Probably also the BaBar model was fitted in this way. 16 Option 4 was used for fits of TAUOLA RChL.
Re-weighting events can be used to estimate the contribution from parts of the model amplitude as well. The ratio of the investigated part of the amplitude squared to the whole amplitude squared provides then the weight. Evaluation of the interferences between parts of the currents is also possible on the event by event basis. One has to evaluate how far from 1 is the sum of weights coming from all parts.
Let us look into that in more detail. For any integrated distributions, a bin always represents the average weight of events filling it. Once a histogram is filled, information from separate events is lost. Because of that, interference observed in any histogram is reduced with respect to its full size. The effect in the integrated distribution is dependent on the selection of bins. By comparing interference coming from the event by event analysis and integrated distributions one may evaluate how viable the used distributions are to constrain interference and if the model can be explored with the available data sufficiently well. The integrated contributions from parts of the current are identical for any distribution, but interferences may differ vastly. An analysis of option 1 is expected to give a very poor insight on interferences, and options 2-4 can give only limited information. Option 5 is expected to give the best way to explore models in that regard. Yet, it still does not explore the whole structure of multidimensional phase-space, contrary to the methods of Ref. [4].

TAUOLA CLEO, TAUOLA CLEO isospin intricate, CLEO publ.
Beside the functional form of the currents, the numerical parameters affect the shapes of the distributions. Among them, mainly masses, widths, and couplings of the intermediate resonances differ. That is why we collect all constants to the last detail, taking into account variants in the implementations. The set of parameters 18 is extracted from the TAUOLA code, to make sure that they were used in the TAUOLA CLEO and TAUOLA CLEO isospin intricate currents. These currents rely on an analysis of the data using option 5, but further cross-checks using option 6 were also mentioned in the CLEO publications. The published experimental data for τ → π 0 π 0 π − ν τ were used by the collaboration. The predictions for τ → π − π − π + ν τ were documented only in conference contributions [10,11]. It is not the theoreticians' or the MC authors' right to decide whether predictions for the π − π − π + channel should be obtained from the π 0 π 0 π − case at the level of distribution or taking into account isospin rotation of intermediate states. An actual choice must rely on an analysis of the systematic errors. Our default cur- Table 1 Parameters differing in TAUOLA and Ref. [6]. These differences in parameters are the reason of the differences visible in Fig to improve the distributions. We then use the distribution level isospin symmetry to get predictions for π − π − π + from π 0 π 0 π − ; we use the same current for the two cases. The second variant, TAUOLA isospin intricate, was also prepared by the collaboration but was not published. We have nonetheless distributed it together with all versions of TAUOLA but not made available for direct TAUOLA use. 19 To activate it, a modification of the code was always necessary. Let us stress that the numerical difference between the two options is comparable to/smaller than the differences to other parameterizations presented in later sections. We believe also that the choice lies within the systematic uncertainties. This would of course change with CLEO publication or with the future publications of other collaborations. For further discussion and a numerical comparison of different current options, see the following subsections and Sects. 11 and 12. The model described in Ref. [6] and the implementation used in TAUOLA differ only by the parameters stored in Table 1. One should stress that Ref. [6] is for the π 0 π 0 π − mode, whereas TAUOLA CLEO uses the same current for the π − π − π + mode as well.
The analytic form of the TAUOLA CLEO isospin intricate current for τ → π − π − π + ν τ can be understood as a simple consequence of isospin symmetry applied to the τ → π 0 π 0 π − ν τ case. But we must make sure that the hadronic current to be transformed by isospin rotation is sufficiently well established; see the discussion in Sect. 14.1.
One could wonder if it is not better to drop the terms of σ , f 0 , and f 2 from the current parameterization, so currents for both channels would have the same functional form. Indeed, that was often the case. For example, the parametrizations TAUOLA BaBar and TAUOLA RChL 2012, introduced later, used identical currents for π − π − π + and π 0 π 0 π − . In TAUOLA CLEO the same current for both channels was used disregarding the fact that it was seemingly introducing σ of a charge ±2. In fact, it was simply assumed that σ was just a feature of the current parametrization and not a real physical state. This is supported by the observation that on the Dalitz plot from Fig. 3 of Ref. [10] there seems to be a resonancelike structure in s 3 , the invariant mass of the π − π − pair. In the corresponding plot, Fig. 5 of Ref. [6] for the π 0 π 0 π − case, the same structure can be seen in π 0 π 0 . It is most likely originating in part from interferences, but it can be parameterized as a resonance. 20 Let us now look into the numerical differences between the basic variants.
9.1 Numerical differences of TAUOLA CLEO and CLEO publ.
The MC sample presented in [6] is not available now, therefore we can only try to recreate the results. To do so, we need to identify all differences in the initializations. It appears that only a few parameters differ; see Table 1. The resulting difference is small but clear; see Fig. 1. Note that these two variants used fixed values for the masses and widths of the resonances (while the RChL and BaBar models fitted them). We can only speculate if it is either due to stability 21 of the fit procedure or simply the belief that use of the results from other (non-τ ) measurements is better. Apart from extreme kinematic regions, the largest differences can be observed in the π − π − π + distribution. We can also investigate the differences in terms of threedimensional plots. Figure 2 presents the Dalitz plots obtained from TAUOLA CLEO. As we do not have data nor an MC generated sample from the CLEO collaboration, only a human eye comparison is possible with their published results. As mentioned Fig. 2 of our paper should coincide 22 with Fig. 5 of Ref. [6]. No convincing differences can be seen.

Numerical differences of TAUOLA CLEO and TAUOLA CLEO isospin intricate
These two currents have the same values of parameters but differ by the functional form of the hadronic current for τ → π − π − π + ν τ . TAUOLA CLEO is represented by Eqs. (5) and (6), and TAUOLA CLEO isospin intricate by Eqs. 20 See also Sect. 14.2 where we discuss some possible origins of unexpected structures. 21 Overpopulation of the parameters may affect the stability. If that was not the case, we may suspect lack of computing power as the additional parameters increase the CPU time needed for the fit to converge substantially. 22 Note that the results of Ref. [6] measurements are for the π 0 π 0 π − mode, but the same functional form of the hadronic current is used for the π − π − π + in TAUOLA CLEO, which is why such a comparison is valid.  Fig. 1 Invariant mass, GeV units, of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO with numerical constants shifted to values of Ref. [6] (red) and TAUOLA CLEO with default parameterization (green). For the black line, which representsthe ratio of the red plot to green one, the left hand side scale should be used. The right hand side scale represents the number of events per bin of 10 MeV Each Dalitz plot is a distribution for TAUOLA CLEO in the s 1 , s 2 variables (GeV 2 units). s 1 is taken to be the highest of the two possible values of M 2 π − π + in each event are at the 10 % level, that is, at the level of the theoretical uncertainty for use of the isospin symmetry [1].

Further variants of currents used by CLEO
The currents described in this paper are not the only ones developed by the experimental collaborations. We would like to hint why certain solutions were used and documented in publications, while others seem to be forgotten. Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO (red) and TAUOLA CLEO isospin intricate (green). For details of the histograms' definition see the caption of Fig. 1 We should start the discussion with CLEO publ. This current is fully described in Ref. [6]. But the very same publication mentions a couple of other different variants, stating for some of them even better agreement with the data. So why is it that none of those did become the main model of the collaboration?  Fig. 4 Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA RChL (red) and TAUOLA RChL 2012 (green). For details of the histograms' definition, see the caption of Fig. 1 First, we should mention the variant of the model where the σ mass and width were fit parameters. The best agreement 23 with the data was obtained with M σ = 555 MeV 23 Negative doubled logarithmic likelihood was 43 units below that of the nominal fit. and σ = 540 MeV (p. 22 of Ref. [6]). Such a fit caused a shift in all couplings not exceeding 20 %. We may again only speculate that it was believed that it is better to use theoretically predicted or experimentally well-established values. Fitting of sigma mass and width was probably performed in order to check how seriously this resonance should be treated. The result, different values than used for the nominal model (M σ = 860 MeV and σ = 880 MeV), was most likely not very encouraging.
Second, there is mentioned in Ref. [6] a study of the pseudo-scalar contribution from π (1300). It also improved the fit result, but with the available data its existence could not be established with unquestioned significance. The very same reasoning was present in a study of the a 1 contribution, which also was not included in the nominal model.
The Ph.D. thesis [12] presents a similar study of the pseudo-scalar π (1300) in τ → π − π − π + ν τ decay. Different ways of inclusion of the scalar current were tested there. Similarly agreement with the data got better, but it was not considered significant enough. In this study the hadronic current used for the non-scalar part was the same as TAUOLA CLEO isospin intricate. It suggests again that at some point the current without pseudo-scalar contribution was used by the collaboration.
Another thing under consideration by CLEO was the resonance radius. Their study showed little importance. Therefore, it was set to 0 in order to simplify the model. Keeping things as simple as possible seems to be important for many choices of the collaborations.
The BaBar collaboration had all CLEO currents available, yet decided to use simpler currents. Apart from simplicity, one can suspect another reason behind such a choice. The BaBar collaboration uses commonly one-dimensional histograms in the analysis even in our times [34]. The inclusion of any resonances other than ρ and ρ is not so easy to exploit without three-dimensional distributions. On Fig. 5 one can see that differences are not very compelling. They are mostly below 10 %. A comparison of the data and the theoretical model on scatter-grams like Fig. 6 can show more clearly the importance of other resonances, but it was never performed. On the latter figure, the ratio of Dalitz plots coming from BaBar model and CLEO model is shown. While the ratio of the invariant mass distributions was mostly between 0.9 and 1.1, the ratio in some areas of the Dalitz plots escapes the 0.5-2 range. Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO (red) and TAUOLA BaBar (green). For details of the histograms definition, see the caption of Fig. 1 code and how they were used. Table 2 presents those defined globally, outside the form factors. Table 3 presents those introduced inside the hadronic current only. The parameters in this table were kept constant during the fitting procedure of CLEO, but were fitted for BaBar modeling. Table 4 stores the parameters fitted to the data. They are also defined inside the hadronic current. Note that some parameters can be different for diverse parts of the code, for example the mass of the ρ resonance. This may be misleadingly worrisome, but it should not be considered as an error. Parameters used outside the hadronic form factors affect the optimization of the efficiency of phase-space generators, that is, essentially, the speed of generation only, or they are well-established constants like m π . Their variation does not affect the output distributions, or the variations are simply unphysical.
The TAUOLA CLEO code is essentially as described in Ref. [6] and named 'nominal fit'. It differs by some numerical constants only; see Table 1. The CLEO collaboration did not fit the masses and widths of the resonances but only     Table 4 Parameters used in the definition of the CLEO hadronic current in TAUOLA obtained by the CLEO collaboration from a fit to the data the coupling constants β i . Such a choice was driven by theoretical assumptions and measurements of the masses elsewhere. This may be the best solution when we lack good enough measurements of state used in the model. Let us recall, Sect. 10, how insensitive the model can be to the σ meson properties.

Numerical comparison of TAUOLA CLEO and TAUOLA BaBar
Here, we present a numerical comparison of TAUOLA CLEO and TAUOLA BaBar. On Fig. 5, the distributions of the invariant masses of π − π + , π − π − , π − π − π + are shown. Looking at the plots one may notice that, except for the lack of some resonances mentioned earlier, also the lack of the K * K threshold in the BaBar model is visible in the distribution of the π − π − π + mass (and even more visible on the ratio plot) in the 1.4 GeV region. Moreover, a different a 1 mass plays a role. In addition, on Fig. 6, the ratio of the Dalitz plots is shown. Lack of f 0 , f 2 , σ resonances in the model used by BaBar can also be seen there. The differences are greater than for comparison of TAUOLA CLEO isospin intricate and TAUOLA CLEO; see Fig. 1. The invariant mass of the π − π + system exhibits the largest difference in terms of the ratio.

Numerical comparison of TAUOLA CLEO isospin intricate and TAUOLA BaBar
For completeness we provide a comparison of TAUOLA CLEO isospin intricate and TAUOLA BaBar. The results are very similar to the previous comparison and can be seen on Fig. 7. We do not discuss those results as no new aspects appear. The largest difference is present in the π − π + invariant mass spectrum.

Comparison of TAUOLA CLEO and TAUOLA RChL
In this section, a numerical comparison of the CLEO and RChL models used in TAUOLA is shown. We will not recall the analytic form of the RChL current. A full description of the TAUOLA RChL 2012 current is given in Ref. [20] (see also [13]), and for TAUOLA RChL in Ref. [14]. The CLEO and RChL models differ a lot in the way of construction and also in the parameters, but one may easily guess what affects the spectra in the way as seen on Fig. 8. Most importantly, the inclusion and parameters of the σ resonance are different. The σ of lower mass and width than in the CLEO model affects mainly the low mass region of M(π − π + ) in RChL, while in the CLEO model it affects the whole spectra. Second, a clearly visible difference appears in the π − π − π + invariant mass distribution. It comes from different modeling of a 1 . CLEO translated theoretical modeling into a polynomial of Eq. (7). RChL keeps the direct model coded, but the parameters of K * meson were not fitted, as the π − π − π + channel should not be very sensitive to them. A reasonable  Fig. 7 Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO isospin intricate (red) and TAUOLA BaBar (green). For details of the histograms' definition see the caption of Fig. 1 modification of those parameters would require a simultaneous fit of the K K π channel. Such a fit could possibly reduce the difference in area of the K * K threshold mentioned in the previous section.  Fig. 8 Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO (red) and TAUOLA RChL (green). For details of the histograms' definition, see the caption of Fig. 1 In Fig. 9, we can find three-dimensional ratio of RChL to CLEO model results. On that plot, we can see areas where difference between the models exceeds factor of 2. This is much bigger difference than seen on one-dimensional distributions. As CLEO model was fitted to three-dimensional distribution it is a very worrisome observation from the RChL perspective. It may hint that even though one-dimensional BaBar data distributions are described well [14], the underlying physics model may be imperfect. An actual judgment should not be completed without comparison with experimental distributions of higher dimensionality.

Comparison of TAUOLA RChL 2012 and TAUOLA RChL
For completeness, we present a comparison of the first version of the RChL model and a new one. Because of the comparison with the data, especially inability to obtain satisfactory agreement in π − π + invariant mass, when all invariant mass distributions were fitted simultaneously TAUOLA RChL 2012 was modified. Those two models differ mainly by inclusion of the σ resonance in the TAUOLA RChL version. Other differences are the model parameters and a minor correction for the Coulomb interaction. The resulting differences in the distributions can be seen on Fig. 4. At this point we must emphasize that those two versions of RChL model were fitted to different data-sets. Invariant mass of π − π − system became available for fits only after development of TAUOLA RChL 2012. This model was fitted to invariant masses of π − π + , π − π − π + systems only, while TAUOLA RChL to invariant masses of π − π + , π − π − , π − π − π + . It was not the only difference. Also binning of histograms changed. New version was fitted to histograms of 10MeV bin width. 24 Distributions used in fit of TAUOLA RChL 2012 had bin width of 20MeV.

Numerical effects of interferences within currents
So far, we have concentrated on presentation of different currents used for Monte Carlo simulations of τ decays. We have looked into their analytical form as well as the numerical values of the parameters fixed by confrontations with experimental data. Comparisons obtained from different parameterizations were reported as well. In the present section let us look into numerical consequences of some features of the currents. We will concentrate on interferences between distinct parts of the current, which can be separated, with the help of isospin transformation.
Final states of long-living scalars with given multiplicity and charge combination, can be formed through intermediate resonances used in description of the hadronic current. One cannot think of them as sum of non-coherent contributions. Interferences represent important aspect of the hadronic current which should be understood and controlled. Even if used for fits experimental distributions are correctly reproduced, unnatural response to modifications of individual contributions to the current and interference can hint for modeling flaws.
In principle, addition of a new resonance to the hadronic current should result in negative interference below the peak and positive above. 25 There is no clear prediction on how strong such resonance should interfere. In case of a narrow intermediate resonance we expect such interference to cancel out in contributions to total rates or relatively inclusive distributions. Especially in the case of wide resonances, we may model a process with something that is not of a proper Table 5 Comparison of contributions from the isoscalars from the remaining part of the current and from interference, depending on the way of investigation. In each entry, the first number is from the spin one intermediate states, the second is from isoscalars, and the third one is from interference. For one-dimensional distributions, the interference effect of the histogram with the largest effect is taken. In the case of the "Absolute width differential", the average of (a module of) interference over all events is shown

Current
Absolute width 1D distributions 3D distributions Absolute width differential π 0 π 0 π − CLEO 75. In this section, we will investigate RChL and CLEO modeling of the three-pion channels. We concentrate on interference of the part depending on isoscalars 26 with the remaining part of the current. These two parts of the hadronic current transform differently under isospin rotation. In principle it is also a test on how justified (or not) is isospin rotation from π 0 π 0 π − to π − π − π + currents.
Let us look at the contribution to the rate. For π 0 π 0 π − , the ρ part, the σ, f 2 , f 0 part and their interference contribution are, respectively, 27 75.1, 4.4 and 20.5 %. Usage of TAUOLA CLEO isospin intricate current for π − π − π + gives 74, 5.8 and 20.2 %. Interference looks slightly different for fully differential distributions. Then, it is 22.6 % in TAUOLA CLEO, and 20.7 % for the TAUOLA CLEO isospin intricate case, when the method as explained in Sect. 8 is used. For the three-dimensional scatter-grams of the CLEO collaboration style, the interference is, respectively, 22.3 and 20.6 %. This demonstrates the nearly perfect sensitivity of such a method.
In Table 5 we collect such results also for other versions of the currents. One must remember that the interference visible in the histograms is somewhat dependent on binning. We do not have the same binning as the CLEO collaboration; they had less bins per energy unit in every distributions. Therefore, while they probably could have seen the same behavior, the interference observed was most likely smaller. Interference is an important numerical property of the distributions, which cannot be ignored while making the choices for the phenomenological approach. 28 26 The contribution only from σ in RChL, while for CLEO from σ, f 2 , f 0 . 27 When the same current is used for π − π − π + corresponding contributions are 75.3, 4.3 and 20.4 %. 28 Later, in particular in Sect. 14.2 we discuss how resonances are deformed and even can produce artificial peaks. Therefore, resonances are not always needed to emulate peaks present in the experimental distribution. Also, because the shapes of the resulting modifications to distributions are distinct with respect to, e.g., Gaussian peaks of the  Fig. 10 Invariant mass of (from the top) π − π + , π − π − , π − π − π + using TAUOLA CLEO isospin intricate (left) and TAUOLA CLEO (right). A green plot represents the contribution from the isoscalars (σ, f 2 , f 0 ) and blue from the remaining part. The red plot represents the whole current. Blue and green plots would sum up to the red one, in the absence of interference. The percentages of the three contributions are given in Table 5 Having discussed CLEO modeling, let us investigate the same problem in the RChL model. On Fig. 11, one can compare the contributions from different parts of the current on one-dimensional distributions. We can see that the plot follows the theoretical expectation of how interference should Footnote 28 continued resonance parameterizations, we may be concerned about the choice between TAUOLA CLEO and TAUOLA CLEO isospin intricate. This is especially so, as the dominant contribution of the isoscalars is the result of interference. A confrontation with the experimental data should be used to resolve this matter.  Fig. 11 Invariant mass of (from the top) π − π + , π − π − , π − π − π + using TAUOLA RChL current. The green plot represents the contribution from the isoscalars and the blue one from the remaining part. The red plot represents the whole current. Blue and green plots should sum up to a red one in the absence of interference. One can observe negative interference for the RChL model whereas it was almost everywhere positive for TAUOLA CLEO isospin intricate; see Fig. 10 and Table 5 look. Also, there is much less interference than for CLEO models; see Table 5. Moreover, the size of interference for the one-prong and three-prong channels is not as similar as in CLEO modeling. While both CLEO models have almost the same level of interference in both channels, RChL has twice more in the three-prong channel. Also, the size of the σ contribution differs a lot between the one-and the three-prong channels, while it is close to the same in both CLEO modelings. Another property worth noting is the limited sensitivity of the RChL model to interference when only one-dimensional distributions are used. It means that the most recent fit of this model does not fully control 29 the influence of σ and the fit on the three-dimensional distributions is strongly desired.
As both models aim at describing the same decay channels we should be concerned that the interplay of the internal parts of the current is quite different. It proves that we do not have a good understanding of the physics behind it, especially what exactly the σ state is, in at least one of the two approaches. Both models find it important [6,14] to include σ , but its description is far from consistent. 30

Interplay of experimental and theoretical input
In the process of creating a model for the hadronic current to be used in MC there are two aspects. Theoreticians provide a model for the hadronic currents and experimentalists use it and fit parameters. Usually, theoreticians concentrate on the consistency of the formulation and prefer to avoid ad hoc phenomenological corrections. Very often, model building is disconnected from data analysis. On the other hand, experimentalists would like to have a model perfectly describing the data. At the same time, they must have control over any possible sources of experimental errors. Therefore, they may prefer to sacrifice some model assumptions for the sake of simplicity and more straightforward systematics analysis. 31 As experiment controls the data sample it is experimentalists to give their "stamp of approval" for the hadronic current. 29 The reason why the interference effect is better visible in one-dimensional distributions in the case of CLEO model than in RChL is technically quite simple. In the first case the mass of the σ is larger, so the region of the phase-space above the peak is restricted, and no cancellation occurs. Nonetheless it is important to realize such features while evaluating the suitability of experimental input for the model fit. It should not be assumed that the conclusion will remain the same for the whole range of the fit parameters. 30 We may speculate that, in part, this difficulty in the control contribution from the isoscalars was the reason behind there being introduced later a much simpler current in the BaBar Monte Carlo simulations. Note that in Ref. [14] the parameters of the σ contributions were featuring errors from the fit larger than for the other resonances. We will return to this point in Sect. 14. 31 In particular that there are no strongly correlated parameters in the fits. Such a thing was done for the π 0 π 0 π − current of CLEO modeling but is lacking for the π − π − π + case. Only conference contributions describe some modelings of this channel. In the case of the CLEO modeling we may only speculate if there were some unresolved problems with this channel or it was simply lack of the manpower that stopped the collaboration publication from appearing. In the following subsections we will discuss some aspects of model building that can play a role.

Ambiguities of experimental inputs
The isospin symmetry between the τ decay channels can be restricted to a relation of the final states only. As a consequence, for intermediate states, such a choice could be understood as including a contribution of, e.g., σ → π − π − for the channel τ → ν τ π − π − π + . This, of course, would mean an error. However, as σ was not a well-established (and very broad) resonance at the time of the CLEO work, this may not be unreasonable. Back then, f 0 (400-1200) [35] with a width ranging from 400-1000 MeV was expected. Later on, in the years 2002-2012 f 0 (600) was expected [36]. Recently, in 2012 it was shifted to f 0 (500) [37], thanks to numerous calculations, e.g., [38,39], claiming both a mass and a width close to 500 MeV. In current measurements this resonance remains extremely wide, while the CLEO collaboration used an even wider (880 MeV) and more massive (860 MeV) σ , as of Ref. [40].
That is why, from the experimentalist point of view, one can think of this amplitude contribution as a heuristic modification for the current, of not much quantum number meaning. 32 Such an interpretation may hold, unless it is excluded by comparisons with the experimental data. In the meantime, such an option in the Monte Carlo can be helpful to study if the contribution is indeed the σ . Alternatively, it is possible to think that a broad σ is just a pretext to a correct current of a 1 , ρ intermediate states; thus it may be identical for one-and three-prong τ decay modes. Similarly one can think of f 0 . The f 2 is not as broad, but even in this case contributions to s 1 , s 2 spectra are not clearly localized in the phase-space. We should also mention that both f 2 and f 0 have limited phasespace for decay; therefore their Breit-Wigner distributions are strongly deformed. On top of that, in Sect. 13 we have shown that the isoscalars in the models we have presented do not contribute as separate resonances to the decay. They contribute positive interference over almost the whole spectra. That is why such an imposed isospin symmetry, understood as used for the same current for the π − π − π + and the 32 In Ref. [6] we read The form factors as defined in Eqn. A2 (Equation A2 of Ref. [6] correspond to our Eq. (1)) do not have simple correspondence with those that can be associated with specific resonant contributions to the hadronic current, which seems to support our statement.  Fig. 12 The invariant mass of the π − π + system in τ → π − π − π + decay using TAUOLA CLEO isospin intricate (green), current with ten times multiplied amplitude from isoscalars (red) and also with width of isoscalars reduced ten times (blue) π 0 π 0 π − decay channels of τ , can be justified. Let us recall that isospin symmetry may be broken and be of limited predictive power starting from 5 to 10 % precision level (Ref. [1], Section 5.11). Also, there is very small contribution to the currents from isospin violating decays. 33 The CLEO II detector [43] had only limited K ± /π ± separation at the time of work on Ref. [6]. The three-prong τ decays were challenging in terms of systematical errors. The TAUOLA CLEO isospin intricate model is based on an analysis of the τ → π 0 π 0 π − ν τ data. The contribution from isoscalars is expected to be bigger for the three-prong channel (see Sect. 4). Thus, one needs to extrapolate from a small effect to a larger one. This magnifies the errors on the isoscalars. We expect that such issues could have been behind the lack of a CLEO publication for the τ → π − π − π + ν τ decay channel.

Further concerns for model-data confrontation
Let us look into details of the TAUOLA CLEO isospin intricate current, and the decay channel τ → π − π − π + ν τ . The contribution from the isoscalars may be missed from the distribution of the π − π + invariant mass spectrum: no clean features are seen; Fig. 12, green line. If we artificially increase by a factor of 10 the isoscalars' amplitude, they dominate the shape of the π − π + invariant mass spectrum, see the red line, but clear peaks are still absent. The position of the peaks is 33 In the τ → 3πν τ decay mode isospin violating signals are expected [41]. There are two decays that can feed the 3π spectrum. For the threeprong channel it is τ → ωπ − ν τ with subsequent ω → π − π + decay. The one-prong channel can be fed by the decay τ → ηπ 0 π − ν τ (due to η-π 0 mixing). Those decays are estimated to give small contribution to the total rate at the level of 0.4% and 10 −3 %, respectively. Studies of τ → ηπ 0 π − ν τ case were also done within RChL scheme [42].
shifted with respect to the masses used in the propagator 34 of the hadronic current. One can see, blue line, that two peaks emerge clearly if the widths of the isoscalars are reduced by a factor of 10 as well. Looking at all three lines, one realizes that the appearing low energy peak/deformation present in TAUOLA CLEO isospin intricate (in other modelings as well) is for a large part indirect. There are two possible combinations of π − π + and both are put into the histogram. The lower peak-like bump comes from the small phase-space available for π − when the system of the other π − and π + is produced with a high invariant mass. That forces a low invariant mass bump for the π − π + system, even though those two pions do not originate from a common intermediate resonance.
The resonances at the end of the phase-space represent a challenging part of the description. As we could see from Fig. 12, their shapes are substantially deformed from the Breit-Wigner ones. Another effect is also possible. In Ref. [44] it was demonstrated how deformations due to resonances, acting as phase-space constraints, may result even in artificial resonance-like peaks. This observation is another reason why one may be hesitant to exploit isospin symmetries of currents in an uncritical way.
Another concern is the size of contribution due to interference. It may be much larger than the effect of extra resonances alone. At the same time, it is not fully explored by inclusive distributions. See Table 5; for RChL, the interference effect is much larger for the three-dimensional histograms than for the one-dimensional ones. For the CLEO case, the difference was rather minor. One has to remember such effects while introducing new resonances into the currents. It must be carefully studied if they affect shapes of the distributions in a controlled way.
Finally, if numerical effects are small and weakly depend on the model parameters, Monte Carlo methods may not be best suited for the studies-especially if systematic errors are to be taken into account [45]. If available, semi-analytical methods are better. Therefore, we have prepared a special arrangement, which is briefly presented in "Appendix A".

Summary
In this paper, we have investigated several variants of parameterizations of the hadronic currents used in Monte Carlo programs for simulations of the τ lepton decays. As a reference for comparisons we have used TAUOLA CLEO parameterization, because of its well-established status; it is supported by the collaboration publication. The first investigated variant was TAUOLA CLEO isospin intricate, as prepared for TAUOLA but not commonly used. This parameterization is documented by conference contributions only. We have evaluated the numerical size of the differences, which were found to be of the order of the systematic errors expected from the use of isospin symmetry alone. We have looked for the origin of numerical differences. We studied both analytic forms of the distributions as well as the dependence on numerical constants.
Later, we have presented comparisons with further parameterizations, in particular the variants of TAUOLA RChL 2012, TAUOLA RChL, and TAUOLA BaBar. All distributions were obtained with the help of Monte Carlo generations. We have prepared semi-analytical distributions to calculate such results as well. This is a necessary step toward fitting with the help of multidimensional distributions, using variants of the model-independent methods presented in Ref. [4], but adopted to the case when τ leptons are relativistic (that is, when the ν τ momentum cannot be reconstructed easily).
We have discussed the limitations and constraints on the use of the models and symmetries like isospin symmetry resulting from insufficiently differential distributions of experimental data. We have pointed to large numerical consequences of interferences. We have also pointed out problems related to phase-space constraints for the decay products of broad resonances.
We should accept that the predictions of the RChL models may agree with the data with a precision not worse than 1 N C , that is, about 30 % (with respect to the total rate) only. Indeed, all of the presented parameterizations pass such a condition and agree with each other at that level as well. The precision of the experimental data is at least one order of magnitude better. That is why one has to concentrate the effort on arrangements for convenient comparisons of the data with model predictions. It is important that models are prepared in a sufficiently flexible way (the necessary future adaptations can be introduced fast and the correlations of parameters are understood) before the fits to the experimental data are performed. Also the experimental distributions need to be investigated, if they are detailed enough to constrain the underlying physics.
We have discussed the consequences of the limited dimensionality of the distributions used in data representations for the results of the fits. We could see in Ref. [14] that the model of the τ → 3πν τ decay modes, prepared only a year earlier [20] and fitted to two one-dimensional, mass histograms, turned out to predict results for the invariant mas of π − π − pair rather poorly. It had to be updated. We are expecting to repeat this once three-dimensional experimental histograms are addressed. There is a factor of 2 discrepancy between RChL models and CLEO models results for some bins of the three-dimensional distributions (which were used by CLEO also as input to their fits). One-dimensional histograms may thus not be enough as experimental input for a precise parameterization of the hadronic currents. The three-dimensional ones of the type as used by CLEO [6] may be sufficient for τ → 3πν τ . However, this conjecture may need to be checked before serious work on the τ decay channels into three scalars with kaons, like that of Refs. [46,47], is completed. Then the full complexity of the hadronic currents is expected, and an analysis similar to the one of Ref. [4] may be indispensable. Nonetheless, attempts to construct currents on the basis of fits to one-dimensional (three-dimensional CLEO style) distributions are of interest. It can demonstrate the predictive power of models such as RChL beyond 1 N C precision level. On a practical side, once differences and options are understood in the context of one-dimensional distributions, we are better prepared to fit the high precision data of today.
If operators like those in Ref. [4] are needed, then they have to be adopted for conditions of τ leptons produced relativistically (making a reconstruction of the ν τ momentum non-trivial). Our paper documents also a step in such a direction; see also the web page [48].
Even though the purpose of our paper is to review the available options for hadronic currents used in a simulation of the τ decays into 3π , some speculations on why such choices were introduced in the past were unavoidable. This is a by far incomplete aspect of this work and may even be misleading at some points. There is little documentation on internal discussions of the experimental groups behind the actual choice of the analytical form of current parameterizations used in data analysis and publications. Often complicated and sometimes unfinished activities, where questions of great importance of the systematic errors for multidimensional distributions would rise, could not be referred to.
In general, we think experiments should be responsible for fits approving the models' quality. We hope our paper will contribute to the discussion how the relations between model developers, experimental physicists, and people working on Monte Carlo and fitting programs should look.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

A Analytical distributions
From Ref. [14], documenting the RChL current and its fit to the experimental data, it became obvious that having semi- analytical distributions available simplifies and speeds up the fitting procedure. As a part of this work we have prepared a program to calculate the distributions from a given hadronic current by integrating the structure function W A (as defined in Ref. [4]) within the limits of a chosen bin. The program uses the IntegralMultiple function of the ROOT [49] libraries, allowing for the acquisition of both one-dimensional and Dalitz distributions. We have checked the code to give proper results for TAUOLA CLEO, TAUOLA CLEO isospin intricate, TAUOLA BaBar, and TAUOLA RChL. The program offers an easy way for a re-fit of these models to the data of the form of three-dimensional distributions like in the CLEO publication. Also, it is much easier than if the MC method is used to perform an analysis of the statistical and systematical errors with the semi-analytical results free from statistical errors themselves. See Ref. [14] for details. The IntegralMultiple routine of the ROOT library has several input parameters. In particular, the minimum and maximum number of function evaluations requested for each bin. Those parameters affect the precision and the time of integration; therefore one has to find a balance between the two. We have found that integration has some stability problems if too low a number of requested points is chosen. Taking into account time, precision, and stability of the integration we choose a minimum 5000 points and a maximum 50,000 points per bin. With such a setup it takes 6 s (of CPU time on a typical PC processor) for TAUOLA CLEO and 16 s for TAUOLA RChL to get one-dimensional distributions and, respectively, 47 and 120 s to get Dalitz plots with the same binning as used in this paper. One has to add 2 s for the CLEO and 30 s for the RChL models for reinitialization whenever the model parameters are changed. If higher precision is needed, doubling the number of points per bin requires 50 % more time for three-dimensional plots, while it does not noticeably affect the calculation of the onedimensional distributions. We expect that the CPU time for a Dalitz plot calculation can still be optimized. A fit with Dalitz plots is more desirable; it gives a better test of the Each Dalitz plot is a distribution for Pythia CLEO in the s 1 , s 2 variables (GeV 2 units). s 1 is taken to be the highest of the two possible values of M 2 π − π + in each event investigated model and has a simpler calculation of the systematics; each individual event enters only one bin, contrary to one-dimensional histograms, where it enters each 35 of four histograms of the invariant masses. In Ref. [45] an optimization of the fitting procedure for TAUOLA RChL in the case of one-dimensional histograms is presented. We 35 Each event enters histograms for the invariant masses of π − π − , π − π − π + and twice for π − π + systems (as there are two possible pairs of π − π + ). hope that with a similar approach we could prepare a framework allowing for fast fitting to three-dimensional distributions for any theoretical model as well. Note that the current initialization of RChL is prepared for multiple channels of τ decays, including not yet published ones with kaons [46,47]. These have not been confronted with the data yet. For fitting to π − π − π + data, we can skip the running routines needed for these other channels, to speed it up to the level which was used in Ref. [45]. We have not checked if IntegralMultiple will provide smooth distributions with variation of the model parameters. It is important for fitting algorithms that the derivatives with respect to the model parameters are continuous, as we have found while working for Ref. [14].  Fig. 15 Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO (red) and Pythia CLEO (green). For details of the histograms' definition, see the caption of Fig. 1

B Comparison of TAUOLA CLEO and Pythia CLEO
In this "Appendix", we present a comparison of an alternative implementation of the CLEO model for τ → π − π − π + ν τ as done in PYTHIA 8.201, with the TAUOLA implementation.  16 Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO (red) and Pythia CLEO (green), when the same parameters are enforced for both generators. For details of the histograms' definition, see the caption of Fig. 1 Thanks to a comparison we can identify the source of the differences and estimate its impact on the generated event samples.
By investigation of the Pythia code, several differences in parameterization were found and stored in Table 6. The  Fig. 17 Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO isospin intricate (red) and Pythia CLEO (green), when the same parameters are enforced for both generators. For details of the histograms' definition, see the caption of Fig. 1 parameters not mentioned in that table are the same as for TAUOLA. The hadronic current of Pythia CLEO is of the same form as described in Sect. 4 for TAUOLA CLEO isospin intricate. The internal arrangements of the code are slightly  13,21) different. Due to this, in the structure of the code some constants are folded with others. The parameter a 1 is one of such. It is not coded explicitly, but we can extract its value: 0.784468371 GeV. Figures 2 and 13 represent Dalitz plots from TAUOLA and Pythia, respectively. Figure 14 shows the ratio of Figs. 2, 13. One-dimensional distributions of the invariant mass of π − π + , π − π − , π − π − π + for both generators can be seen on Fig. 15. There is a systematic difference between the generators. Disagreement may have appeared due to the different form of the hadronic current, or different numerical parameters. To test this, for both generators, constants of Ref. [6] Fig. 19 Invariant mass of (from the top) π − π + , π − π − , π − π − π + in TAUOLA CLEO isospin intricate (red) and Pythia CLEO (green). For details of the histograms' definition, see the caption of Fig. 1 were taken. Agreement was not improved substantially, with the exception of the invariant mass π − π − π + , where the influence of the a 1 mass and width turned out to be the main source of difference, as can be seen from Fig. 16. In general, agreement only to about 10 % can be concluded except for some less populated bins at the ends of the π − π + mass spectra.
It suggests that the functional form of the hadronic current is of higher numerical importance than the model constants.
To test this, we have used the TAUOLA CLEO isospin intricate current, while having the same parameters in both generators. As in the previous test, the constants were set that of [6]. The resulting distributions of the invariant masses as well as their ratios can be seen in Fig. 17: nearly perfect agreement can be concluded. 36 To complete our investigation, a comparison of TAUOLA CLEO isospin intricate and Pythia CLEO, with default initializations was performed. The generators use the hadronic current of Eqs. (8) and (9), and only the numerical constants differ; see those collected in Table 6. Figure 19 shows the resulting distributions of the invariant masses and their ratios. One can see the biggest difference in the invariant mass of π − π − π + due to the distinct a 1 mass and width, which are of the highest importance among all numerical constants.
Further, technical checks were performed: the same numerical constants and the same but simplified hadronic current were installed in both generators. Those included: -reducing the hadronic current to F 1 of Eq. (5) only (F 2 , F 3 = 0); -setting the hadronic current to constant; -reducing the hadronic current to the Breit-Wigner case of a 1 only.
In all these cases no discrepancies between results of TAUOLA and Pythia generators beyond statistical fluctuations were found. Samples of 3 million events were used. All considered, the differences in the distributions produced by TAUOLA and Pythia for the τ → π − π − π + ν τ decay channel result only from the parameters if the same variant of current is used in TAUOLA, otherwise the generators perform identically. All differences are insignificant in comparison to the uncertainties of the parameterization.

C CLEO style results of three-dimensional distribution
All over the paper we present three-dimensional distributions for three scalar final states, following the definition used by Such a representation does not constrain the full differential nature of the hadronic currents, but it was used in publications of experimental data of the most differential nature until now. It is an important reference point (Figs. 1-19). Each Dalitz plot is a distribution for TAUOLA CLEO isospin intricate in the s 1 , s 2 variables (GeV 2 units). s 1 is taken to be the highest of the two possible values of M 2 π − π + in each event