diff --git a/tf_pwa/amp/ampgen_FOCUS.py b/tf_pwa/amp/ampgen_FOCUS.py
new file mode 100644
index 0000000..00bc313
--- /dev/null
+++ b/tf_pwa/amp/ampgen_FOCUS.py
@@ -0,0 +1,95 @@
+import numpy as np
+import tensorflow as tf
+
+from tf_pwa.amp import Particle, register_particle
+
+from .ampgen_pipi_swave import constructKMatrix, phsp_FOCUS, pol, poleConfig
+
+
+@register_particle("Kpi_Swave")
+class KPiSwaveKmatrix(Particle):
+    """
+
+    Kpi S wave model from AmpGen (https://github.com/GooFit/AmpGen/blob/master/src/Lineshapes/FOCUS.cpp).
+
+    .. plot::
+
+        >>> import matplotlib.pyplot as plt
+        >>> plt.clf()
+        >>> from tf_pwa.utils import plot_particle_model
+        >>> ax = plot_particle_model("Kpi_Swave")
+        >>> ax = plot_particle_model("Kpi_Swave", params={"lineshape_modifier": "KEta"}, axis=ax)
+        >>> ax = plot_particle_model("Kpi_Swave", params={"lineshape_modifier": "I32"}, axis=ax)
+        >>> _ = ax[1].legend(["Kpi", "KEta", "I32"])
+
+    """
+
+    def __init__(self, name, **kwargs):
+        self.lineshape_modifier = "Kpi"
+        super().__init__(name, **kwargs)
+
+    def __call__(self, m, **kwargs):
+        s = m * m
+        ret = FOCUS_fun(s, self.lineshape_modifier)
+        return ret
+
+    def get_amp(self, data, data_c, **kwargs):
+        m = data["m"]
+        return self(m)
+
+
+def FOCUS_fun(s, lineshapeModifier="Kpi"):
+    I = 1.0j
+    sInGeV = s
+    mK = 0.493677  # ParticlePropertiesList::get( "K+" )->mass() ;
+    mPi = 0.13957018  # ParticlePropertiesList::get( "pi+" )->mass() ;
+    mEtap = 0.95766  # # ParticlePropertiesList::get( "eta'(958)0" )->mass();
+    sNorm = mK * mK + mPi * mPi
+    s12 = 0.23
+    s32 = 0.27
+    I12_adler = (sInGeV - s12) / sNorm
+    I32_adler = (sInGeV - s32) / sNorm
+
+    poleConfigs = [poleConfig(np.array(1.7919 + 0j), [0.31072, -0.02323])]
+
+    rho1 = phsp_FOCUS(sInGeV, mK, mPi)
+    rho2 = phsp_FOCUS(sInGeV, mK, mEtap)
+    X = (sInGeV / sNorm) - 1
+
+    kMatrix = constructKMatrix(sInGeV, 2, poleConfigs)
+    scattPart = [
+        pol(X, [0.79299, -0.15099, 0.00811]),
+        pol(X, [0.15040, -0.038266, 0.0022596]),
+        pol(X, [0.15040, -0.038266, 0.0022596]),
+        pol(X, [0.17054, -0.0219, 0.00085655]),
+    ]
+    scattPart = tf.reshape(tf.stack(scattPart, axis=-1), kMatrix.shape)
+
+    kMatrix = kMatrix + scattPart
+
+    I12_adler = tf.cast(I12_adler, kMatrix.dtype)
+    I32_adler = tf.cast(I32_adler, kMatrix.dtype)
+
+    K11 = I12_adler * kMatrix[..., 0, 0]
+    K12 = I12_adler * kMatrix[..., 0, 1]
+    K22 = I12_adler * kMatrix[..., 1, 1]
+    K32 = I32_adler * pol(X, [-0.22147, 0.026637, -0.00092057])
+
+    detK = K11 * K22 - K12 * K12
+
+    del_ = 1 - rho1 * rho2 * detK - I * (rho1 * K11 + rho2 * K22)
+
+    T11 = 1.0 - I * rho2 * K22
+    T22 = 1.0 - I * rho1 * K11
+    T12 = I * rho2 * K12
+
+    T32 = 1 / (1 - I * K32 * rho1)
+    if lineshapeModifier == "Kpi":
+        return (K11 - I * rho2 * detK) / del_
+    elif lineshapeModifier == "KEta":
+        return K12 / del_
+    elif lineshapeModifier == "I32":
+        return T32
+    else:
+        print("P-vector component : ", lineshapeModifier, " is not recognised")
+        return tf.complex(tf.ones_like(s), tf.zeros_like(s))
diff --git a/tf_pwa/amp/ampgen_pipi_swave.py b/tf_pwa/amp/ampgen_pipi_swave.py
new file mode 100644
index 0000000..3cfabec
--- /dev/null
+++ b/tf_pwa/amp/ampgen_pipi_swave.py
@@ -0,0 +1,308 @@
+import numpy as np
+import tensorflow as tf
+
+from tf_pwa.amp import Particle, register_particle
+
+
+@register_particle("pipi_Swave")
+class PiPiSwaveKmatrix(Particle):
+    """
+
+    pipi S wave model from AmpGen (https://github.com/GooFit/AmpGen/blob/master/src/Lineshapes/kMatrix.cpp).
+    using the parameters from `DtoKKpipi_v2.opt` (https://github.com/GooFit/AmpGen/blob/master/options/DtoKKpipi_v2.opt)
+
+    .. plot::
+
+        >>> import matplotlib.pyplot as plt
+        >>> plt.clf()
+        >>> from tf_pwa.utils import plot_particle_model
+        >>> ax = plot_particle_model("pipi_Swave", params={"all_tokens": ["pole.0"]})
+        >>> ax = plot_particle_model("pipi_Swave", params={"all_tokens": ["poleKK.1"]}, axis=ax)
+        >>> ax = plot_particle_model("pipi_Swave", params={"all_tokens": ["prodKK.0"]}, axis=ax)
+        >>> ax = plot_particle_model("pipi_Swave", params={"all_tokens": ["prod.1"]}, axis=ax)
+        >>> _ = ax[1].legend(["pole.0", "poleKK.1", "prodKK.0", "prod.1"])
+
+    """
+
+    def __init__(self, name, **kwargs):
+        self.all_tokens = []
+        for i in range(5):
+            self.all_tokens.append(f"pole.{i}")
+        for i in range(5):
+            self.all_tokens.append(f"prod.{i}")
+        self.kmatrix_params = {}
+        self.particleName = "PiPi00"
+        super().__init__(name, **kwargs)
+
+    def init_params(self):
+        super().init_params()
+        self.all_var = []
+        for tokens in self.all_tokens:
+            self.all_var.append(self.add_var(tokens, is_complex=True))
+        self.all_var[0].fixed(1.0)
+
+    def get_params_vecter(self):
+        return tf.stack([i() for i in self.all_var], axis=-1)
+
+    def __call__(self, m, **kwargs):
+        s = m * m
+        params_vector = self.get_params_vecter()
+        ret = kMatrix_fun(
+            s, self.all_tokens, self.kmatrix_params, self.particleName
+        )
+        return tf.reduce_sum(ret * params_vector, axis=-1)
+
+    def get_amp(self, data, data_c, **kwargs):
+        m = data["m"]
+        return self(m)
+
+
+def to_matrix(x):
+    return tf.expand_dims(tf.expand_dims(x, axis=-1), axis=-1)
+
+
+def d_prod(a, b):
+    return tf.expand_dims(a, axis=-1) * tf.expand_dims(b, axis=-2)
+
+
+def pol(x, p):
+    ret = tf.math.polyval([tf.cast(i, x.dtype) for i in p], x)
+    return tf.complex(ret, tf.zeros_like(ret))
+    # F = tf.zeros_like(x)
+    # L = 1
+    # for ip in p:
+    # F = F + ip * L;
+    # L = L * x;
+    # return tf.complex(F, tf.zeros_like(F))
+
+
+def complex_sqrt(x):
+    return tf.sqrt(tf.complex(x, tf.zeros_like(x)))
+
+
+def phsp_twoBody(s, m0, m1):
+    return complex_sqrt(1.0 - (m0 + m1) * (m0 + m1) / s)
+
+
+def phsp_fourPi(s):
+    """
+    Parameterisation of the 4pi phase space taken from `Laura` (https://laura.hepforge.org/  or Ref. https://arxiv.org/abs/1711.09854)
+    """
+    mPiPlus = 0.139570
+
+    rho_4pi = pol(
+        s, [0.00051, -0.01933, 0.13851, -0.20840, -0.29744, 0.13655, 1.07885]
+    )
+
+    return tf.where(s > 1, phsp_twoBody(s, 2 * mPiPlus, 2 * mPiPlus), rho_4pi)
+
+
+def phsp_FOCUS(s, m0, m1):
+    mp = m0 + m1
+    mm = m0 - m1
+    return complex_sqrt((1.0 - mp * mp / s) * (1.0 - mm * mm / s))
+
+
+def gFromGamma(m, gamma, rho):
+    return tf.sqrt(m * gamma / rho)
+
+
+def getPropagator(kMatrix, phaseSpace):
+    nChannels = kMatrix.shape[-1]
+    Id = tf.eye(nChannels, dtype=kMatrix.dtype)
+    T = Id - 1j * kMatrix * tf.expand_dims(phaseSpace, axis=-2)
+    return tf.linalg.inv(T)
+
+
+def constructKMatrix(this_s, nChannels, poleConfigs):
+    kMatrix = tf.zeros((nChannels, nChannels), dtype=tf.complex128)
+    for pole in poleConfigs:
+        num = d_prod(pole.couplings, pole.couplings)
+        dom = pole.s - tf.cast(this_s, pole.s.dtype)
+        term = tf.cast(num, dom.dtype) / to_matrix(dom)
+        kMatrix = kMatrix + term
+    return kMatrix
+
+
+class poleConfig:
+    def __init__(self, s, couplings=None):
+        self.s = s
+        if couplings is None:
+            self._couplings = []
+        else:
+            self._couplings = couplings
+
+    def add(self, x):
+        self._couplings.append(x)
+
+    @property
+    def couplings(self):
+        return tf.stack(self._couplings, axis=-1)
+
+
+def kMatrix_fun(
+    s, all_tokens=["pole.0"], new_params={}, particleName="PiPi00"
+):
+    sInGeV = s
+    nPoles = 5
+    nChannels = 5
+    channels = ["pipi", "KK", "4pi", "EtaEta", "EtapEta"]
+
+    mPiPlus = 0.139570
+    mKPlus = 0.493677
+    mEta = 0.547862
+    mEtap = 0.967780
+
+    all_tokens = [i.split(".") for i in all_tokens]
+
+    params = default_params.copy()
+    params.update(new_params)
+
+    def Parameter(name, value=0.0):
+        return params.get(name, value)
+
+    def paramVector(name, n):
+        return tf.stack([Parameter(name + str(i)) for i in range(n)])
+
+    sA0 = Parameter("sA0", -0.15)
+    sA = Parameter("sA", 1.0)
+    s0_prod = Parameter(particleName + "_s0_prod", -0.07)
+    s0_scatt = Parameter("s0_scatt", -3.92637)
+    fScatt = paramVector("f_scatt", nChannels)
+    poleConfigs = []
+    addImaginaryMass = Parameter("kMatrix::fp", True)
+
+    for pole in range(1, nPoles + 1):
+        stub = "IS_p" + str(pole) + "_"
+        mass = Parameter(stub + "mass")
+        # add a tiny imaginary part to the mass to avoid floating point errors //
+        if addImaginaryMass:
+            p = poleConfig(
+                tf.cast(mass * mass, dtype=tf.complex128) + (1.0e-6j)
+            )
+        else:
+            p = poleConfig(tf.cast(mass * mass, dtype=tf.complex128))
+        for ch in range(nChannels):
+            p.add(Parameter(stub + channels[ch]))
+        poleConfigs.append(p)
+
+    phaseSpace = tf.stack(
+        [
+            phsp_twoBody(sInGeV, mPiPlus, mPiPlus),
+            phsp_twoBody(sInGeV, mKPlus, mKPlus),
+            phsp_fourPi(sInGeV),
+            phsp_twoBody(sInGeV, mEta, mEta),
+            phsp_twoBody(sInGeV, mEta, mEtap),
+        ],
+        axis=-1,
+    )
+
+    kMatrix = constructKMatrix(sInGeV, nChannels, poleConfigs)
+
+    scattPart = tf.zeros((nChannels, nChannels))
+    slow_term = (1 - s0_scatt) / (s - s0_scatt)
+    fScatt_scale = fScatt * np.array([0.5] + [1.0] * (nChannels - 1))
+    scattPart = tf.expand_dims(fScatt_scale, axis=-1) + tf.expand_dims(
+        fScatt_scale, axis=-2
+    )
+    scattPart = tf.cast(scattPart, slow_term.dtype) * to_matrix(slow_term)
+
+    kMatrix = kMatrix + tf.cast(scattPart, kMatrix.dtype)
+    # kMatrix.imposeSymmetry(0,1);
+
+    adlerTerm = (
+        (1.0 - sA0) * (sInGeV - sA * mPiPlus * mPiPlus / 2) / (sInGeV - sA0)
+    )
+
+    kMatrix = tf.cast(to_matrix(adlerTerm), kMatrix.dtype) * kMatrix
+    F = getPropagator(kMatrix, phaseSpace)
+
+    all_ret = []
+    for tokens in all_tokens:
+        if tokens[0] == "scatt":
+            i = int(tokens[1])
+            j = int(tokens[2])
+            M = tf.mathmul(F, kMatrix)
+            ret_i = M[..., j, i]
+        elif tokens[0].startswith("pole"):
+            if len(tokens[0]) == 4:
+                idx = 0
+            else:
+                idx = channels.index(tokens[0][4:])
+            pTerm = int(tokens[1])
+            pole = poleConfigs[pTerm]
+            M = tf.reduce_sum(
+                F[..., idx, :] * tf.cast(pole.couplings, F.dtype), axis=-1
+            )
+            ret_i = M / (pole.s - tf.cast(sInGeV, pole.s.dtype))
+        elif tokens[0].startswith("prod"):
+            if len(tokens[0]) == 4:
+                idx = 0
+            else:
+                idx = channels.index(tokens[0][4:])
+            pTerm = int(tokens[1])
+            pd = (1 - s0_prod) / (sInGeV - s0_prod)
+            ret_i = F[..., idx, pTerm] * tf.cast(pd, F.dtype)
+        else:
+            print(
+                "Modifier not found: , expecting one of {scatt, pole, poleKK, prod, prodKK}"
+            )
+            ret_i = tf.complex(tf.zeros_like(s), tf.zeros_like(s))
+        all_ret.append(ret_i)
+    return tf.stack(all_ret, axis=-1)
+
+
+# from DtoKKpipi_v2.opt
+params_str = """
+KK00_s0_prod	2	-0.165753	0
+KK10_s0_prod	2	-0.165753	0
+PiPi00_s0_prod	2	-0.165753	0
+PiPi20_s0_prod	2	-0.165753	0
+
+IS_p1_4pi	2	0	0
+IS_p1_EtaEta	2	-0.39899	0
+IS_p1_EtapEta	2	-0.34639	0
+IS_p1_KK	2	-0.55377	0
+IS_p1_mass	2	0.651	0
+IS_p1_pipi	2	0.22889	0
+IS_p2_4pi	2	0	0
+IS_p2_EtaEta	2	0.39065	0
+IS_p2_EtapEta	2	0.31503	0
+IS_p2_KK	2	0.55095	0
+IS_p2_mass	2	1.2036	0
+IS_p2_pipi	2	0.94128	0
+IS_p3_4pi	2	0.55639	0
+IS_p3_EtaEta	2	0.1834	0
+IS_p3_EtapEta	2	0.18681	0
+IS_p3_KK	2	0.23888	0
+IS_p3_mass	2	1.55817	0
+IS_p3_pipi	2	0.36856	0
+IS_p4_4pi	2	0.85679	0
+IS_p4_EtaEta	2	0.19906	0
+IS_p4_EtapEta	2	-0.00984	0
+IS_p4_KK	2	0.40907	0
+IS_p4_mass	2	1.21	0
+IS_p4_pipi	2	0.3365	0
+IS_p5_4pi	2	-0.79658	0
+IS_p5_EtaEta	2	-0.00355	0
+IS_p5_EtapEta	2	0.22358	0
+IS_p5_KK	2	-0.17558	0
+IS_p5_mass	2	1.82206	0
+IS_p5_pipi	2	0.18171	0
+f_scatt0	2	0.23399	0
+f_scatt1	2	0.15044	0
+f_scatt2	2	-0.20545	0
+f_scatt3	2	0.32825	0
+f_scatt4	2	0.35412	0
+s0_prod	2	-1	0
+s0_scatt	2	-3.92637	0
+sA	2	1	0
+sA0	2	-0.15	0
+"""
+
+default_params = {}
+for i in params_str.splitlines():
+    if len(i.strip()) == 0:
+        continue
+    name = i.strip().split("\t")
+    default_params[name[0]] = float(name[2])
diff --git a/tf_pwa/amp/core.py b/tf_pwa/amp/core.py
index 0c78b63..8cfd14b 100644
--- a/tf_pwa/amp/core.py
+++ b/tf_pwa/amp/core.py
@@ -1242,10 +1242,8 @@ def ls_selector_qr(decay, ls_list):
     for l1 in p1.spins:
         for l2 in p2.spins:
             if abs(l1 - l2) <= p0.J:
-                if decay.p_break:
-                    if p1.J > p2.J and l1 > 0:
-                        continue
-                    if p1.J <= p2.J and l2 > 0:
+                if not decay.p_break:
+                    if (-l1, -l2) in hel_list:
                         continue
                 hel_list.append((l1, l2))
     from sympy import Matrix