From 7ef8018c6e0bf1b648822a3814f30c158ec8c6c8 Mon Sep 17 00:00:00 2001
From: Matej Hudec <77328256+hoodyn@users.noreply.github.com>
Date: Fri, 29 Jan 2021 01:17:58 +0100
Subject: [PATCH 1/7] Fixing the LFV decays KL,KS -> l1+ l2-

The original `amplitudes` function renamed to `amplitudes_weak_eigst` as it actually only yields amplitudes relevant for decays of weak eigenstates K0, K0bar.

A new `amplitudes` function, which newly takes also `K` (to be `'KL'` or `'KS'`) as an argument, has been implemented. In merely makes proper linear combinations of amplitudes for K0 and K0bar.

For l1=l2, `amplitudes` yields the imag/real parts of the S, P weak_eigs coefficients, as it should. Thus, taking real/imag parts is no longer tackled in `amplitudes_eff`, which now serves merely to add up short- and long-distance contributions, as it should, I believe.

The bug which lead to error in `amplitudes_eff` for 'KS' and l1!=l2, has been fixed.
---
 flavio/physics/kdecays/kll.py | 54 ++++++++++++++++++++++++++++-------
 1 file changed, 43 insertions(+), 11 deletions(-)

diff --git a/flavio/physics/kdecays/kll.py b/flavio/physics/kdecays/kll.py
index 2d873b31..bda9f2d6 100644
--- a/flavio/physics/kdecays/kll.py
+++ b/flavio/physics/kdecays/kll.py
@@ -9,15 +9,13 @@
 from math import pi, sqrt
 
 
-def amplitudes(par, wc, l1, l2):
-    r"""Amplitudes P and S entering the $K\to\ell_1^+\ell_2^-$ observables.
+def amplitudes_weak_eigst(par, wc, l1, l2):
+    r"""Amplitudes P and S for the decay $\bar K^0\to\ell_1^+\ell_2^-$.
 
     Parameters
     ----------
-
     - `par`: parameter dictionary
     - `wc`: Wilson coefficient dictionary
-    - `K`: should be `'KL'` or `'KS'`
     - `l1` and `l2`: should be `'e'` or `'mu'`
     """
     # masses
@@ -37,6 +35,40 @@ def amplitudes(par, wc, l1, l2):
     return xi_t * P, xi_t * S
 
 
+def amplitudes(par, wc, K, l1, l2):
+    r"""Amplitudes P and S entering the $K_{L,S}\to\ell_1^+\ell_2^-$ observables.
+    
+    Parameters
+    ----------
+    
+    - `par`: parameter dictionary
+    - `wc`: Wilson coefficient dictionary
+    - `K`: should be `'KL'` or `'KS'`
+    - `l1` and `l2`: should be `'e'` or `'mu'`
+    """
+    # KL, KS are linear combinations of K0, K0bar. So are the amplitudes.
+    S_K0bar, P_K0bar = amplitudes_weak_eigs(par, wc, l1, l2)
+    if l1 != l2:
+        S_aux, P_aux = amplitudes_weak_eigs(par, wc, l2, l1)
+        S_K0 = -S_aux.conjugate()
+        P_K0 = P_aux.conjugate()
+        if 'K' == 'KL':
+            sig = +1
+        elif 'K' == 'KS':
+            sig = -1
+        S = (S_K0 + sig * S_K0bar) / 2
+        P = (P_K0 + sig * P_K0bar) / 2
+    # Save computing time for special cases. See also arXiv:1711.11030.
+    elif l1 == l2:
+        if 'K' == 'KL':
+            S = -1j * S_K0bar.imag 
+            P = P_K0bar.real
+        elif 'K' == 'KS':
+            S = -S_K0bar.real
+            P = -1j * P_K0bar.imag
+    return P, S
+
+
 def amplitudes_LD(par, K, l):
     r"""Long-distance amplitudes entering the $K\to\ell^+\ell^-$ observables."""
     ml = par['m_' + l]
@@ -53,18 +85,18 @@ def amplitudes_LD(par, K, l):
 
 def amplitudes_eff(par, wc, K, l1, l2, ld=True):
     r"""Effective amplitudes entering the $K\to\ell_1^+\ell_2^-$ observables."""
-    P, S = amplitudes(par, wc, l1, l2)
+    P, S = amplitudes(par, wc, K, l1, l2)
     if l1 != l2 or not ld:
         SLD = 0
         PLD = 0
     else:
         SLD, PLD = amplitudes_LD(par, K, l1)
-    if K == 'KS' and l1 == l2:
-        Peff = P.imag
-        Seff = S.real + SLD
-    if K == 'KL':
-        Peff = P.real + PLD
-        Seff = S.imag
+    if K == 'KS':
+        Peff = P
+        Seff = S + SLD
+    elif K == 'KL':
+        Peff = P + PLD
+        Seff = S
     return Peff, Seff
 
 

From f202571b40b67d7c5671eaf928ca03477cf7aaba Mon Sep 17 00:00:00 2001
From: Matej Hudec <77328256+hoodyn@users.noreply.github.com>
Date: Fri, 29 Jan 2021 01:24:28 +0100
Subject: [PATCH 2/7] Update comments in kll.py

There was an odd comment above a normalization coefficient, so I've updated it.
---
 flavio/physics/kdecays/kll.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/flavio/physics/kdecays/kll.py b/flavio/physics/kdecays/kll.py
index bda9f2d6..dbe84f2d 100644
--- a/flavio/physics/kdecays/kll.py
+++ b/flavio/physics/kdecays/kll.py
@@ -31,6 +31,7 @@ def amplitudes_weak_eigst(par, wc, l1, l2):
     CSm = wc['CS_'+qqll] - wc['CSp_'+qqll]
     P = (ml2 + ml1)/mK * C10m + mK * CPm  # neglecting mu, md
     S = (ml2 - ml1)/mK * C9m  + mK * CSm  # neglecting mu, md
+    # Include complex part of the eff. operator prefactor. Phases matter.
     xi_t = ckm.xi('t', 'sd')(par)
     return xi_t * P, xi_t * S
 
@@ -127,7 +128,7 @@ def br_kll(par, wc_obj, K, l1, l2, ld=True):
     mK = par['m_K0']
     tauK = par['tau_'+K]
     fK = par['f_K0']
-    # appropriate CKM elements
+    # CKM part of the eff. operator prefactor included in Peff and Seff
     N = 4 * GF / sqrt(2) * alphaem / (4 * pi)
     beta = sqrt(lambda_K(mK**2, ml1**2, ml2**2)) / mK**2
     beta_p = sqrt(1 - (ml1 + ml2)**2 / mK**2)

From f04c0da793d754ab844ff68e746f9af10ffd92af Mon Sep 17 00:00:00 2001
From: Matej Hudec <77328256+hoodyn@users.noreply.github.com>
Date: Fri, 29 Jan 2021 01:36:19 +0100
Subject: [PATCH 3/7] Update sending long-distance contibutions to KL,KS->l+l-
 in kll.py

Originally, the function `amplitudes_LD(par, K, l)` was always sending both long-distance coefficients - for KL and KS, despite the fact that it takes `K` as an argument; `amplitudes_eff` then always thrown one of them away.
Now, it should be clearer.
---
 flavio/physics/kdecays/kll.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/flavio/physics/kdecays/kll.py b/flavio/physics/kdecays/kll.py
index dbe84f2d..a4b74424 100644
--- a/flavio/physics/kdecays/kll.py
+++ b/flavio/physics/kdecays/kll.py
@@ -77,11 +77,15 @@ def amplitudes_LD(par, K, l):
     s2w = par['s2w']
     pre = 2 * ml / mK / s2w
     # numbers extracted from arXiv:1711.11030
-    ASgaga = 2.49e-4 * (-2.821 + 1.216j)
-    ALgaga = 2.02e-4 * (par['chi_disp(KL->gammagamma)'] - 5.21j)
-    S = pre * ASgaga
-    P = -pre * ALgaga
-    return S, P
+    if K == 'KS':
+        ASgaga = 2.49e-4 * (-2.821 + 1.216j)
+        SLD = +pre * ASgaga
+        PLD = 0
+    elif K == 'KL':
+        ALgaga = 2.02e-4 * (par['chi_disp(KL->gammagamma)'] - 5.21j)
+        SLD = 0
+        PLD = -pre * ALgaga
+    return SLD, PLD
 
 
 def amplitudes_eff(par, wc, K, l1, l2, ld=True):
@@ -92,12 +96,8 @@ def amplitudes_eff(par, wc, K, l1, l2, ld=True):
         PLD = 0
     else:
         SLD, PLD = amplitudes_LD(par, K, l1)
-    if K == 'KS':
-        Peff = P
-        Seff = S + SLD
-    elif K == 'KL':
-        Peff = P + PLD
-        Seff = S
+    Peff = P + PLD
+    Seff = S + SLD
     return Peff, Seff
 
 

From 5425de8eccb54cb3180ada56bc3d70cc3a4d97c9 Mon Sep 17 00:00:00 2001
From: Matej Hudec <77328256+hoodyn@users.noreply.github.com>
Date: Fri, 29 Jan 2021 01:48:29 +0100
Subject: [PATCH 4/7] Standardize normalization of KL,KS->l1l2 amplitude

Taking an extra factor of 2 from the prefactor in `br_kll` and putting a factor of sqrt(2) to both short distance `amplitudes` and `amplitudes_LD`. This does not affect the result, of course. Now, the formula for the BR is the more standard one (the same that is used for B->l1 l2 in flavio), and also combining weak-K0-eigenstate amplitudes to those of KL or KS is more intuitive now.
---
 flavio/physics/kdecays/kll.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/flavio/physics/kdecays/kll.py b/flavio/physics/kdecays/kll.py
index a4b74424..5f752471 100644
--- a/flavio/physics/kdecays/kll.py
+++ b/flavio/physics/kdecays/kll.py
@@ -57,16 +57,16 @@ def amplitudes(par, wc, K, l1, l2):
             sig = +1
         elif 'K' == 'KS':
             sig = -1
-        S = (S_K0 + sig * S_K0bar) / 2
-        P = (P_K0 + sig * P_K0bar) / 2
+        S = (S_K0 + sig * S_K0bar) / sqrt(2)
+        P = (P_K0 + sig * P_K0bar) / sqrt(2)
     # Save computing time for special cases. See also arXiv:1711.11030.
     elif l1 == l2:
         if 'K' == 'KL':
-            S = -1j * S_K0bar.imag 
-            P = P_K0bar.real
+            S = -1j * S_K0bar.imag * sqrt(2)
+            P = P_K0bar.real * sqrt(2)
         elif 'K' == 'KS':
-            S = -S_K0bar.real
-            P = -1j * P_K0bar.imag
+            S = -S_K0bar.real * sqrt(2)
+            P = -1j * P_K0bar.imag * sqrt(2)
     return P, S
 
 
@@ -75,7 +75,7 @@ def amplitudes_LD(par, K, l):
     ml = par['m_' + l]
     mK = par['m_' + K]
     s2w = par['s2w']
-    pre = 2 * ml / mK / s2w
+    pre = sqrt(2) * 2 * ml / mK / s2w
     # numbers extracted from arXiv:1711.11030
     if K == 'KS':
         ASgaga = 2.49e-4 * (-2.821 + 1.216j)
@@ -133,7 +133,7 @@ def br_kll(par, wc_obj, K, l1, l2, ld=True):
     beta = sqrt(lambda_K(mK**2, ml1**2, ml2**2)) / mK**2
     beta_p = sqrt(1 - (ml1 + ml2)**2 / mK**2)
     beta_m = sqrt(1 - (ml1 - ml2)**2 / mK**2)
-    prefactor = 2 * abs(N)**2 / 32. / pi * mK**3 * tauK * beta * fK**2
+    prefactor = abs(N)**2 / 32. / pi * mK**3 * tauK * beta * fK**2
     Peff, Seff = amplitudes_eff(par, wc, K, l1, l2, ld=ld)
     return prefactor * (beta_m**2 * abs(Peff)**2 + beta_p**2 * abs(Seff)**2)
 

From 77509df0335ecdbcb2adc39f56b9ca9a231a413d Mon Sep 17 00:00:00 2001
From: Matej Hudec <77328256+hoodyn@users.noreply.github.com>
Date: Fri, 29 Jan 2021 01:54:45 +0100
Subject: [PATCH 5/7] Unify m_K0 in amplitudes_LD in bll.py

Like almost everywhere else, there is no need to distinguish between mass of KL and KS in the prefactors in amplitudes_LD.
---
 flavio/physics/kdecays/kll.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/flavio/physics/kdecays/kll.py b/flavio/physics/kdecays/kll.py
index 5f752471..f982822b 100644
--- a/flavio/physics/kdecays/kll.py
+++ b/flavio/physics/kdecays/kll.py
@@ -73,7 +73,7 @@ def amplitudes(par, wc, K, l1, l2):
 def amplitudes_LD(par, K, l):
     r"""Long-distance amplitudes entering the $K\to\ell^+\ell^-$ observables."""
     ml = par['m_' + l]
-    mK = par['m_' + K]
+    mK = par['m_K0']
     s2w = par['s2w']
     pre = sqrt(2) * 2 * ml / mK / s2w
     # numbers extracted from arXiv:1711.11030

From 676d8d22dce6c186cb0ecb6f37a6b7adc6a5238a Mon Sep 17 00:00:00 2001
From: Matej Hudec <77328256+hoodyn@users.noreply.github.com>
Date: Fri, 29 Jan 2021 02:04:11 +0100
Subject: [PATCH 6/7] Fix my recent typo in kll.py

---
 flavio/physics/kdecays/kll.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/flavio/physics/kdecays/kll.py b/flavio/physics/kdecays/kll.py
index f982822b..d5fdf5b7 100644
--- a/flavio/physics/kdecays/kll.py
+++ b/flavio/physics/kdecays/kll.py
@@ -48,9 +48,9 @@ def amplitudes(par, wc, K, l1, l2):
     - `l1` and `l2`: should be `'e'` or `'mu'`
     """
     # KL, KS are linear combinations of K0, K0bar. So are the amplitudes.
-    S_K0bar, P_K0bar = amplitudes_weak_eigs(par, wc, l1, l2)
+    S_K0bar, P_K0bar = amplitudes_weak_eigst(par, wc, l1, l2)
     if l1 != l2:
-        S_aux, P_aux = amplitudes_weak_eigs(par, wc, l2, l1)
+        S_aux, P_aux = amplitudes_weak_eigst(par, wc, l2, l1)
         S_K0 = -S_aux.conjugate()
         P_K0 = P_aux.conjugate()
         if 'K' == 'KL':

From 27b5606db2a0bc62eadceef0aa1848ae2cf3f20b Mon Sep 17 00:00:00 2001
From: Matej Hudec <77328256+hoodyn@users.noreply.github.com>
Date: Fri, 29 Jan 2021 15:31:33 +0100
Subject: [PATCH 7/7] Debugging the typos in comparisons

As @peterstangl noted, comparisons like
```
'K' == 'KL'
```
were present, while the correct version is
```
K == 'KL'
```
---
 flavio/physics/kdecays/kll.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/flavio/physics/kdecays/kll.py b/flavio/physics/kdecays/kll.py
index d5fdf5b7..db19a36d 100644
--- a/flavio/physics/kdecays/kll.py
+++ b/flavio/physics/kdecays/kll.py
@@ -53,18 +53,18 @@ def amplitudes(par, wc, K, l1, l2):
         S_aux, P_aux = amplitudes_weak_eigst(par, wc, l2, l1)
         S_K0 = -S_aux.conjugate()
         P_K0 = P_aux.conjugate()
-        if 'K' == 'KL':
+        if K == 'KL':
             sig = +1
-        elif 'K' == 'KS':
+        elif K == 'KS':
             sig = -1
         S = (S_K0 + sig * S_K0bar) / sqrt(2)
         P = (P_K0 + sig * P_K0bar) / sqrt(2)
     # Save computing time for special cases. See also arXiv:1711.11030.
     elif l1 == l2:
-        if 'K' == 'KL':
+        if K == 'KL':
             S = -1j * S_K0bar.imag * sqrt(2)
             P = P_K0bar.real * sqrt(2)
-        elif 'K' == 'KS':
+        elif K == 'KS':
             S = -S_K0bar.real * sqrt(2)
             P = -1j * P_K0bar.imag * sqrt(2)
     return P, S