HEP-KBFI / tautau-Entanglement

Code for studies of quantum entanglement in Z->tautau and H->tautau events
0 stars 1 forks source link

Discrepancies due to the choice of tau decay software #2

Closed ktht closed 1 year ago

ktht commented 1 year ago

It turns out that depending on the software one uses to decay the tau pair -- Pythia8, Tauola or TauDecay --, the spin correlation matrix can vary quite a bit even if no event selection cuts are applied. Here are the matrices for the simplest tau decay mode (tau->pi+nu) computed by averaging over the scalar products of polarimetric vectors, but the results change very little when other methods for computing the spin correlation matrix are used:

Pythia8 (115k events):

[    -0.249686  0.00508408  0.00391348
    0.00653464    0.138683   -0.276996
   -0.00186421    0.243658    0.885224 ]
+/-
[     0.134348    0.133983    0.135247
      0.144723    0.132779    0.130972
      0.133501    0.135441    0.120462 ]

Tauola (126k events):

[     0.015075  -0.0591435   0.0364924
    -0.0281235  -0.0268309   0.0240461
     0.0477844-0.000387594    0.848719 ]
+/-
[     0.143917    0.131293    0.134135
       0.13149    0.136686      0.1367
       0.14108    0.139452    0.128245 ]

TauDecay (10M events):

[    -0.412114 -0.00075686 0.000226434
   -0.00180427    0.526222  0.00374591
  -0.000500454  0.00162725    0.881012 ]
+/-
[     0.134427    0.137754    0.130989
      0.133747    0.122811    0.124256
      0.130581    0.126311    0.130846 ]
ktht commented 1 year ago

TLDR is given at the very end.

I reran the analysis for the following three cases, where the tau decays are delegated to

  1. Pythia 8.306 (in Belle II conditions);
  2. Tauola++ 1.1.8 out-of-the-box;
  3. Tauola++ 1.1.8 with the SANC library.

All results are extracted from generator-level $\tau\to\pi+\nu$ decays, so to rule out any diluting effects caused by kinematic fit. In all cases, the tau decays are inclusive, meaning that about 110k events survives the requirement that both taus (which are produced in DY events) decay to single pions. Here are the results presented in beam coordinate axis, where the 3rd axis ($\mathbf{k}$) is along tau momentum, the 2nd axis ($\mathbf{r}$) lies on the plane spanned by tau and electron momenta, and the 1st axis ($\mathbf{n}$) orthogonal to the other two, all in c.o.m. frame:

PYTHIA (114631 events)
==================================

Asymmetry:
[    -0.245623  0.00659507  0.00603676
    0.00401288    0.146801   -0.285333
    0.00464098    0.237178    0.871981 ]

Summation:
[    -0.249686  0.00508408  0.00391348
    0.00653464    0.138683   -0.276996
   -0.00186421    0.243658    0.885224 ]

ML fit:
[    -0.249671  0.00921341 0.000212506
    0.00753089    0.139364   -0.276324
   -0.00312754     0.23755    0.886287 ]

1D diff XS:
[    -0.249255  0.00454215  0.00443124
    0.00593798    0.138174   -0.276008
   -0.00217371    0.244018    0.887989 ]

2D diff XS:
[     -0.24979  0.00523523  0.00417337
    0.00672311     0.13868   -0.277246
   -0.00190086     0.24378     0.88803 ]

TAUOLA++ (125591 events)
==================================

Asymmetry:
[    0.0111791  -0.0478697   0.0527108
    -0.0635396  -0.0205429   0.0736677
     0.0388244  -0.0267217    0.842576 ]

Summation:
Matrix C:
[     0.015075  -0.0591435   0.0364924
    -0.0281235  -0.0268309   0.0240461
     0.0477844-0.000387594    0.848719 ]

ML fit:
[    0.0155621  -0.0607487   0.0326157
    -0.0233263  -0.0256499    0.026586
     0.0424554  0.00254594    0.852014 ]

1D diff XS:
[    0.0152708  -0.0572985   0.0382227
    -0.0298607  -0.0268092   0.0243901
     0.0476596 0.000498932    0.849958 ]

2D diff XS:
[    0.0145587  -0.0598416   0.0367065
    -0.0284408  -0.0269722   0.0240464
     0.0479108 0.000526788    0.852377 ]

TAUOLA++ with SANC (125285 events)
==================================

Asymmetry:
Matrix C:
[    -0.353658  -0.0329808   0.0666959
    -0.0337471    0.446566   0.0746139
     0.0401325 -0.00239454    0.736337 ]

Summation:
[    -0.340833  -0.0454045   0.0333969
   -0.00858326    0.444391   0.0194261
     0.0354548  -0.0014837    0.744016 ]

ML fit:
[    -0.343539  -0.0384257   0.0292925
    -0.0120163    0.449523   0.0277464
     0.0359525 -0.00182352     0.75129 ]

1D diff XS:
[    -0.344145  -0.0433198   0.0348278
    -0.0104284    0.445861   0.0210741
     0.0351133-0.000915811    0.744652 ]

2D diff XS:
[     -0.34602  -0.0457655   0.0330568
   -0.00811796    0.446457   0.0196809
     0.0361367 -0.00180716     0.74659 ]

For context, the SM expectation for the spin correlation matrix in Belle II conditions is $\mathrm{diag}(-0.42, 0.53, 0.89)$: tau_spin_correlation_matrix

The elements of matrix R are taken from Eq. 16 of this paper.

Now there's also the SANC library, which lets users generate spin correlation matrices for a number of c.o.m. energies and scattering angles of taus that are produced in DY events. Each spin correlation matrix is given two event weights, one with EW corrections applied and the other without the said corrections. Additional NLO QED corrections can be applied as well. The spin correlation matrices are saved to a text file with some metadata like the SM parameters, which NLO corrections are applied and the author credits. AFAICT, the library generates the spin correlation matrices only for three initial states: $d\bar{d}$, $u\bar{u}$ and $e^+e^-$. If one wants to pass those tables to Tauola++, then the following conventions must be followed:

Tauola++ interpolates between the energies and scattering angles if a given event is generated at slightly different c.o.m. energies or scattering angles than those for which the spin correlation matrix is known. Official documentation of Tauola++ and SANC can be found here (*) (also published here). I used f95 compiler with -std=legacy -ffixed-line-length-132 options to build the auxiliary LoopTools libary. The C object files had to be compiled with -DHAVE_UNDERSCORE option because C functions called from Fortran need to end with an underscore (a detail that I had long forgotten).

As a cross-check, I generated the spin correlation matrices with SANC at c.o.m. energy of 10.58 GeV for 100 different scattering angles and computed their weighted average:

LO:
[[-0.41986834 -0.00002869  0.00000000]
 [-0.00002869  0.52668345 -0.00264794]
 [ 0.00000000 -0.00264794  0.89318489]]

With EW corrections:
[[-0.41986688 -0.00002869  0.00000001]
 [-0.00002869  0.52668162 -0.00271822]
 [ 0.00000001 -0.00271822  0.89318526]]

With QED corrections:
[[-0.41957273  0.00059840  0.00030224]
 [ 0.00059840  0.52587721 -0.01316196]
 [ 0.00030224 -0.01316196  0.89369552]]

Here's my summary:

What else can we do?

Besides the last point I won't make any promises on whether or not I actually deliver them.

ktht commented 1 year ago

So I've been debugging these issues for a solid week. I'll try to summarize my findings here. TLDR is provided at the very end.

Pythia8 In order to understand the issue better, I built Pythia from scratch and linked against CMSSW, which allowed me to add logging to the code. The logs revelead that Pythia tries to first guess the production mode based on the final state: if it sees a pair of OS taus, it assumes that the pair originates from a Z boson. The code then constructs its 4-momentum from the 4-momentum sum of the two taus and sets its PDG ID to 23. By default the TauDecays:mode option is set to 1, which means that this function is triggered first. The function deduces that the mediator particle is a Z boson and chooses the ME of $Z\to\tau^+\tau^-$ to represent the ME of the hard process. In other words, when passing DY events to Pythia, it assumes that the final state leptons originate from an on-shell Z boson. This is the default behavior.

Fortunately, there's an easy fix that skips the erroneous deduction method in favor of another function if the TauDecays:mode option is set to 4. The alternative method checks the input states of the hard process, i.e., the electron-positron pair, and correctly deduces that the hard ME corresponds to $Z/\gamma^*\to\tau^+\tau^-$.

However, it still confuses me why the resulting spin correlation matrix features sizable elements in its off-diagonal. According to Eq. 21 (which keeps only terms linear in $Z\to\tau^+\tau^-$ coupling, $v_\tau\approx=-0.038$) of arXiv:2307.03526 (1), the transverse correlations should have the same magnitude and opposite sign while the off-diagonal elements develop only if the anomalous couplings are substantial. None of these features are present in the spin correlation matrix estimated from the original Pythia sample. That said, I don't think it's worth my time to debug a scenario that we don't need.

Here are the spin correlation matrices extracted from the Pythia sample after the fix:

PYTHIA (116073 events)
==================================

Asymmetry:
[    -0.412947 -0.00844296  -0.0195394
     0.0189881    0.550722   0.0203665
     0.0272587   0.0438689    0.882169 ]

Summation:
[    -0.418821  -0.0104138 -0.00688308
    0.00378321    0.540804    0.009921
     0.0185955    0.016996    0.886184 ]

ML fit:
[    -0.415289  -0.0112852 -0.00506105
    0.00446136    0.541501  0.00737583
     0.0147897    0.018047      0.8845 ]

1D diff XS:
[    -0.415244  -0.0102197 -0.00654984
    0.00392938    0.540056  0.00981594
     0.0186572   0.0166514    0.884126 ]

2D diff XS:
[    -0.416881  -0.0113028 -0.00610333
    0.00355741    0.541128  0.00970768
      0.018666   0.0169275    0.885356 ]

TAUOLA++ with SANC Similarly to Pythia, I rebuilt Tauola++ from scratch. This allowed me to increase the number of bins in the scattering angle from 20 to 200. The results did not change by much:

TAUOLA++ with SANC (125412 events)
==================================

Asymmetry:
[     -0.36947  -0.0468855    0.064938
    -0.0317035    0.441042   0.0717635
     0.0503939 -0.00338086    0.736516 ]

Summation:
[    -0.351581  -0.0616047   0.0289145
   -0.00626352    0.428284   0.0144325
     0.0497008  0.00293103    0.747394 ]

ML fit:
[    -0.355704  -0.0545803   0.0306469
   -0.00901486    0.432133   0.0208993
     0.0506257 -0.00297756    0.752322 ]

1D diff XS:
[    -0.354998   -0.058655    0.029896
   -0.00749712    0.431129   0.0161861
     0.0491126  0.00296605    0.747999 ]

2D diff XS:
[    -0.356939  -0.0617213   0.0287699
   -0.00630209    0.431743   0.0144143
     0.0500239  0.00280499    0.749387 ]

I could not find any errors in the code that does the interpolation between the spin correlation matrices. It still remains a mystery as to why the elements appear off by a constant factor of $1.2$ from theory expectation.

TAUOLA++ with TauSpinner After having a closer look at this paper (2) and the code (3), my conclusion is that like Tauola++ ignores transverse spin correlations, TauSpinner similarly does not calculate weights that includes the effects of transverse spin correlations by default. The main justification presented in the paper is that boosted taus lose the transverse correlation anyways. Obviously for us that's not the case.

Correcting for transverse spin correlations is possible, but analogously to Tauola++ it would still rely on SANC tables for getting the correct spin correlation matrix. Alternatively, one could also implement the analytic formulas directly in the code like it was done in (1). Not only does the latest release of Tauola++ exclude transverse spin effects out-of-the-box, cross-correlations between longitudinal and transverse spins are ignored according to Eq. 28 of (2). The former was fixed in the upstream (3, L606) but there's no release that includes the latest changes (see the changelogs of the development version and the latest release). The code repository is not publicly accessible (but one could just scrape the doxygen documentation for code). It all begs the question why would we go the extent of producing a sample of unpolarized taus to later reweight them if we can get an unweighted sample with the intended spin correlations and spend less effort.

KKMC I tried to build KKMC generator as part of Belle II Analysis Software Framework (basf2) by following these instructions. Turns out that it's a monumental task, because it requires building the complete toolchain (including gcc, python, Geant4, ROOT, Pythia etc) from scratch. This was complicated by the fact that some libraries were missing from the host machine (libxml2, ncurses, redline, texinfo) and had to be rebuilt manually, or that some software packages were fetched form dead URLs (like zlib or MadGraph v3.4.0, which is basically nowhere to be found) or from Belle II mirrors, which I don't have access to. Some packages (like BELLE_FLC and belle_legacy) are not publicly available at all. While I managed to build the toolchain by cutting corners here and there (by e.g. disabling X11 when building ROOT) after spending more than a day on this, I failed to set up the environment correctly, because I cannot seem to instruct gcc to pick up the correct headers when building basf2. I've tried pretty much everything at this point.

The last option would be to build the software in docker, but it's going to be very slow. Unfortunately, I wasn't able to find any recent docker images that already have everything built in. It's worth pointing out that one of the authors of (1) is a full member of Belle II collaboration, which probably explains why they were able to run KKMC with ease.

Comparison of different generators Comparison of Pythia samples before (blue) and after (orange) the fix: compare_pythia

Comparison of Taula++ samples with the default settings (blue), with SANC (orange) and with fine binning in $\cos\theta^*$ (green): compare_tauola

Comparison of Pythia, Tauola++and TauDecay samples: compare_all

TLDR

While TauDecay is extremely fast in generating MC events, if we had to pick a single generator, I'd nevertheless go with Pythia simply because it's the industry standard.

veelken commented 1 year ago

Hi Karl,

thank you very much for the detailed investigation. In my opinion, the case is settled now. I am also in favor of using Pythia8 for the paper. We can mentioned that we did a cross check with TauDecay and Tauola.

I find the distributions of C_nn, C_rr, and C_kk for Tauola a bit worrisome: The factor 1.2 difference in the elements of the C matrix wrt Pythia8 and TauDecay does not seem to be a simple normalization issue, but seems to be caused by subtle differences in the tails of the distribution. I think the difference could be either caused by a bug in the Tauola code, by some approximation of the theory calculations on which Tauola is based, or by a physics effect. Would be great if you could contact the Tauola author Zbigniew Was (z.was@cern.ch) about this difference and ask him for his opinion concerning the cause of the difference.

ktht commented 1 year ago

I finally managed to get basf2 working in docker and generate a bunch of events with KKMC. Skip to the very end if you just want to read the summary.

Setup

Even though I mentioned in my earlier post that there are no docker images available with everything built-in, I did find this image, which has the complete toolchain (i.e., this repository) built and ready to use. This was great because all I had to do was to build basf2. However, it was a bit more complicated than that because of versioning issues: the github repository has no releases that support the version of toolchain that I'm using. In fact, it seems that the repositories hosted on github are mere mirrors of the private repositories hosted on DESY servers where the actual development takes place. The mirrors have not received any updates since May 5th so the developers likely stopped maintaining them, which explains why there are no releases of basf2 that match the externals package. Thus, I had to come up with my own release. Fortunately, the basf2 repository has this file, which encodes the required version of the toolchain. I created a pseudo-release of the framework from the commit that sets the version of the toolchain in that .externals file to what I'm using. After massaging the installation scripts for half a day I came up with these instructions that anyone can follow to build their own basf2 framework. To save an hour of build time I created a new image of my container. Anyone can run it in docker or singularity/apptainer with:

docker run -it belle2/externals-centos7:v01-12-01 /bin/bash
# or
singularity exec docker://ktht/basf2-centos7:v01-12-01-1df2691 /bin/bash

and set up the environment with source /belle2/tools/b2setup $(ls /belle2/releases) inside the container. Here's a screenshot of how it looks like after it's all said and done: belle2_greeting

It's interesting to note that Belle II software stack shares striking similarities to CMSSW. In some cases they even use the same concepts (like GTs and IOVs) and software solutions (CVMFS).

Simulation

I should emphasize here that the basf2 framework is still quite opaque to me; the following should be taken with a slight grain of salt.

I found these three files from the basf2 repository that run the simulation of $e^+e^-\to\tau^+\tau^-$ events with KKMC:

  1. generators/examples/KKGenGenerationOnly.py;
  2. generators/examples/KKGenGenerationOnly_tauola_orig.py;
  3. generators/examples/KKGenGenerationOnly_tauola_bbb.py.

The files are basically the same, except that

As much as I can gather from the context clues, tauinputFile should specify the conditions for generating and decaying tau pairs, while taudecaytableFile specifies only the tau BRs. In short:

Given this information, I'd say that the only difference between the 1st and the 2nd MC sample are the tau BRs (which I expect to be very subtle), while the 3rd sample that decays taus with Tauola-bbb is specifically fine-tuned for Belle II conditions and employs state-of-the-art approximations (which again I don't expect to play a huge role in the final results).

A few comments regarding technicalities:

Results

For every Python configuration file I generated 10M events inclusive in tau decays.

1. KKGenGenerationOnly.py (108623 events)
==================================================

Asymmetry:
[    -0.402309  0.00371929  -0.0198485
     -0.014693    0.543936 -0.00710715
    0.00173076 -0.00364564     0.85643 ]

Summation:
[    -0.397741  0.00844638  0.00261266
    -0.0168413    0.543477  0.00283693
  -0.000526955  0.00330999    0.867469 ]

ML fit:
[    -0.397699  0.00917687   0.0101992
    -0.0157196    0.545113-0.000700433
   -0.00206537  0.00161711    0.867004 ]

1D diff XS:
[    -0.398335  0.00762225  0.00254002
    -0.0164379    0.546706  0.00306423
  -0.000704576  0.00311764    0.865229 ]

2D diff XS:
[    -0.399248  0.00790543  0.00327732
    -0.0169679    0.547768  0.00253293
  -0.000297986   0.0031202    0.867204 ]

2. KKGenGenerationOnly_tauola_orig.py (110391 events)
==================================================

Asymmetry:
[    -0.387387  -0.0132982  0.00590628
    0.00257267    0.529717 -0.00503664
   -0.00307996 -0.00684838    0.872426 ]

Summation:
[    -0.384804 -0.00377286   0.0108233
  -0.000193635    0.520684  -0.0120182
   -0.00904017 -0.00446079    0.876573 ]

ML fit:
[    -0.390271  -0.0027375  0.00705376
  -0.000290429    0.522723   -0.007865
   -0.00727604 -0.00155825    0.875298 ]

1D diff XS:
[    -0.385552 -0.00488933   0.0108551
   -0.00044397    0.518597  -0.0109927
   -0.00905174 -0.00430057    0.870575 ]

2D diff XS:
[    -0.385948 -0.00352788   0.0106701
  -0.000373482    0.520494  -0.0118629
   -0.00942096 -0.00426226    0.872737 ]

3. KKGenGenerationOnly_tauola_bbb.py (103874 events)
==================================================

Asymmetry:
[    -0.399715 -0.00446695   0.0109363
    0.00107823    0.495215   0.0164815
    -0.0119375  0.00808672    0.857269 ]

Summation:
[    -0.388131  0.00147889  0.00192074
   -0.00428194    0.514161  0.00626645
    0.00289267  0.00579045    0.868392 ]

ML fit:
[     -0.38984  0.00617673  0.00445032
  -9.83128e-05    0.513975  0.00198999
    0.00207316 5.48579e-05    0.864278 ]

1D diff XS:
[    -0.388388  0.00101536  0.00270942
   -0.00535866    0.513965  0.00616118
    0.00272844  0.00568301     0.86164 ]

2D diff XS:
[    -0.388949  0.00160133  0.00193946
   -0.00459997    0.515191   0.0067071
    0.00300906   0.0061322    0.863625 ]

Here are the distributions of diagonal elements: compare_kkmc compare_all

The full spin correlations are there, but in most extreme cases can deviate from theory expecation quite a bit: up to 8% in $C{nn}$, up to 6% in $C{rr}$ and up to 4% in $C_{kk}$. There are no striking differences between the samples themselves.


Summary

ktht commented 1 year ago

A small comment regarding the discrepancies seen in the KKMC sample (before I forget the details): one likely reason why those numbers don't match our theory prediction exactly is because KKMC considers bremsstrahlung in spin correlations. Bremsstrahlung (ISR, to be specific) has the effect of reducing c.o.m. energy by 6% down to 9.94 GeV on average and giving a small ~0.5 GeV transverse boost to the ditau system (see the plots below). By making the same assumptions in our theory prediction but reducing the c.o.m. energy by 6% the expectation changes from $\mathrm{diag}(-0.42, 0.53, 0.89)$ to $\mathrm{diag}(-0.41, 0.53, 0.88)$. I guess the residual differences are due to other effects caused by ISR. It would be interesting to see if MG5 v3 simulates those effects (and more) as claimed in this paper. kkmc

ktht commented 1 year ago

One of Tauola++ authors confirmed privately that the C++ interface to Tauola was tested only in HEP conditions that are similar to LHC and FCC, and it's not guaranteed that the interface produces correct results at low energy scales where the effects due to tau lepton mass become more important. However, finding out the actual reason turned out to be a lot more time-consuming than I initially anticipated. So far I've ruled out errors caused by poor discretization in cosTheta (which go away with higher number of bins in cosTheta) and bremsstrahlung effects (which are not present in initial state and identical to Pythia in final state, see the plots below).

cosTheta_beforeAfterFSR mTauTau_beforeAfterFSR

I'm closing the issue now because any further investigation goes well beyond the scope of this project. I'm perfectly fine with the explanation that Tauola++ isn't simply meant for low-energy applications. If we ever want to use Tauola again, we should use it in conjuction with KKMC that is available in the Belle II software framework. It not only models the transverse spin correlations correctly, but also accounts for the effects caused by ISR.