From 6e4981b9574a1a5377c4da2184703f6cc7aa911f Mon Sep 17 00:00:00 2001
From: iasonkrom <iason.krom@gmail.com>
Date: Fri, 22 Mar 2024 16:58:57 +0100
Subject: [PATCH 01/10] add erfexp pdf

---
 tests/test_pdf_cruijff.py          |  2 +-
 tests/test_pdf_erfexp.py           | 89 ++++++++++++++++++++++++++++++
 zfit_physics/models/pdf_cruijff.py |  4 +-
 zfit_physics/models/pdf_erfexp.py  | 89 ++++++++++++++++++++++++++++++
 zfit_physics/pdf.py                |  3 +-
 5 files changed, 183 insertions(+), 4 deletions(-)
 create mode 100644 tests/test_pdf_erfexp.py
 create mode 100644 zfit_physics/models/pdf_erfexp.py

diff --git a/tests/test_pdf_cruijff.py b/tests/test_pdf_cruijff.py
index e81eca2..011fbc1 100644
--- a/tests/test_pdf_cruijff.py
+++ b/tests/test_pdf_cruijff.py
@@ -54,7 +54,7 @@ def test_cruijff_pdf():
     assert cruijff.pdf(tf.range(50.0, 130, 10_000)) <= cruijff.pdf(90.0)
 
 
-def test_cruihff_integral():
+def test_cruijff_integral():
     # Test CDF and integral here
     cruijff, obs = create_cruijff(
         mu=mu_true, sigmal=sigmal_true, alphal=alphal_true, sigmar=sigmar_true, alphar=alphar_true, limits=(50, 130)
diff --git a/tests/test_pdf_erfexp.py b/tests/test_pdf_erfexp.py
new file mode 100644
index 0000000..7a4eab9
--- /dev/null
+++ b/tests/test_pdf_erfexp.py
@@ -0,0 +1,89 @@
+"""Tests for ErfExp PDF."""
+
+import numpy as np
+import pytest
+import tensorflow as tf
+import zfit
+from scipy import integrate, special
+from zfit.core.testing import tester
+
+import zfit_physics as zphys
+
+alpha_true = 90.0
+beta_true = 0.08
+gamma_true = -1.0
+n_true = 0.2
+
+alpha_range = (65.0, 90.0)
+gamma_range = (-10, 10)
+beta_range = (0.01, 10)
+n_range = (0.1, 0.5)
+
+
+def _erfexp_numpy(x, alpha, beta, gamma, n):
+    return special.erfc((x - alpha) * beta) * np.exp(-gamma * (np.power(x, n) - np.power(alpha, n)))
+
+
+erfexp_numpy = np.vectorize(_erfexp_numpy, excluded=["alpha", "beta", "gamma", "n"])
+
+
+def create_erfexp(alpha, beta, gamma, n, limits):
+    obs = zfit.Space("obs1", limits=limits)
+    erfexp = zphys.pdf.ErfExp(alpha=alpha, beta=beta, gamma=gamma, n=n, obs=obs, norm=False)
+    return erfexp, obs
+
+
+def test_erfexp_pdf():
+    # Test PDF here
+    erfexp, _ = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
+    assert erfexp.pdf(90.0).numpy().item() == pytest.approx(
+        erfexp_numpy(90.0, alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true), rel=1e-8
+    )
+    np.testing.assert_allclose(
+        erfexp.pdf(tf.range(50.0, 130, 10_000)),
+        erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),
+        rtol=1e-8,
+    )
+
+
+def test_erfexp_pdf_random_params():
+    # Test PDF here in a loop with random parameters
+    for _ in range(1000):
+        alpha_true = np.random.uniform(*alpha_range)
+        beta_true = np.random.uniform(*beta_range)
+        gamma_true = np.random.uniform(*gamma_range)
+        n_true = np.random.uniform(*n_range)
+
+        erfexp, __ = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
+        np.testing.assert_allclose(
+            erfexp.pdf(tf.range(50.0, 130, 10_000)),
+            erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),
+            rtol=1e-5,
+        )
+
+
+def test_erfexp_integral():
+    # Test CDF and integral here
+    erfexp, obs = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
+    full_interval_numeric = erfexp.numeric_integrate(obs, norm=False).numpy()
+    true_integral = 71.18838
+    numpy_full_integral = integrate.quad(erfexp_numpy, 50, 130, args=(alpha_true, beta_true, gamma_true, n_true))[0]
+    assert full_interval_numeric == pytest.approx(true_integral, 1e-7)
+    assert full_interval_numeric == pytest.approx(numpy_full_integral, 1e-7)
+
+    numeric_integral = erfexp.numeric_integrate(limits=(80, 100), norm=False).numpy()
+    numpy_integral = integrate.quad(erfexp_numpy, 80, 100, args=(alpha_true, beta_true, gamma_true, n_true))[0]
+    assert numeric_integral == pytest.approx(numpy_integral, 1e-7)
+
+
+# register the pdf here and provide sets of working parameter configurations
+def erfexp_params_factory():
+    alpha = zfit.Parameter("alpha", alpha_true)
+    beta = zfit.Parameter("beta", beta_true)
+    gamma = zfit.Parameter("gamma", gamma_true)
+    n = zfit.Parameter("n", n_true)
+
+    return {"alpha": alpha, "beta": beta, "gamma": gamma, "n": n}
+
+
+tester.register_pdf(pdf_class=zphys.pdf.ErfExp, params_factories=erfexp_params_factory)
diff --git a/zfit_physics/models/pdf_cruijff.py b/zfit_physics/models/pdf_cruijff.py
index b0ce2de..4a85977 100644
--- a/zfit_physics/models/pdf_cruijff.py
+++ b/zfit_physics/models/pdf_cruijff.py
@@ -62,8 +62,8 @@ def __init__(
         .. math:
 
             f(x; \\mu, \\sigma_{L}, \\alpha_{L}, \\sigma_{R}, \\alpha_{R}) = \\begin{cases}
-            \\exp(- \\frac{(x-\\mu)^2}{2 \\sigma_{L}^2 + \\alpha_{L} (x-\\mu)^2)}, \\mbox{for} x \\leqslant mu \\newline
-            \\exp(- \\frac{(x-\\mu)^2}{2 \\sigma_{R}^2 + \\alpha_{R} (x-\\mu)^2)}, \\mbox{for} x > mu \\end{cases}
+            \\exp{\\left(- \\frac{(x-\\mu)^2}{2 \\sigma_{L}^2 + \\alpha_{L} (x-\\mu)^2}\\right)}, \\mbox{for} x \\leqslant mu \\newline
+            \\exp{\\left(- \\frac{(x-\\mu)^2}{2 \\sigma_{R}^2 + \\alpha_{R} (x-\\mu)^2}\\right)}, \\mbox{for} x > mu
 
         Args:
             mu: Mean value
diff --git a/zfit_physics/models/pdf_erfexp.py b/zfit_physics/models/pdf_erfexp.py
new file mode 100644
index 0000000..6917a6c
--- /dev/null
+++ b/zfit_physics/models/pdf_erfexp.py
@@ -0,0 +1,89 @@
+from typing import Optional
+
+import tensorflow as tf
+import zfit
+from zfit import z
+from zfit.util import ztyping
+from zfit.z import numpy as znp
+
+
+@z.function(wraps="tensor")
+def erfexp_pdf_func(x, alpha, beta, gamma, n):
+    """Calculate the ErfExp PDF.
+
+    Args:
+        x: value(s) for which the PDF will be calculated.
+        alpha: Location parameter.
+        beta: Scale parameter.
+        gamma: Shape parameter.
+        n: Shape parameter.
+
+    Returns:
+        `tf.Tensor`: The calculated PDF values.
+
+    Notes:
+        Implementation from https://gitlab.cern.ch/cms-muonPOG/spark_tnp/-/blob/Spark3/RooErfExp.cc
+        The parameters beta and gamma are given in reverse order in this c++ implementation.
+    """
+    return tf.math.erfc((x - alpha) * beta) * znp.exp(-gamma * (znp.power(x, n) - znp.power(alpha, n)))
+
+
+class ErfExp(zfit.pdf.BasePDF):
+    _N_OBS = 1
+
+    def __init__(
+        self,
+        alpha: ztyping.ParamTypeInput,
+        beta: ztyping.ParamTypeInput,
+        gamma: ztyping.ParamTypeInput,
+        n: ztyping.ParamTypeInput,
+        obs: ztyping.ObsTypeInput,
+        *,
+        extended: Optional[ztyping.ExtendedInputType] = False,
+        norm: Optional[ztyping.NormInputType] = None,
+        name: str = "ErfExp",
+    ):
+        """ErfExp PDF, the product of a complementary error function and an exponential function.
+
+        Implementation from https://gitlab.cern.ch/cms-muonPOG/spark_tnp/-/blob/Spark3/RooErfExp.cc
+
+        .. math:
+
+            f(x; \\alpha, \\beta, \\gamma, n) = \\text{erfc}(\\beta (x - \\alpha)) \\exp{(-\\gamma (x^n - \\alpha^n))}
+
+        Args:
+            alpha: Location parameter.
+            beta: Scale parameter.
+            gamma: Shape parameter.
+            n: Shape parameter.
+            obs: |@doc:pdf.init.obs| Observables of the
+               model. This will be used as the default space of the PDF and,
+               if not given explicitly, as the normalization range.
+
+               The default space is used for example in the sample method: if no
+               sampling limits are given, the default space is used.
+
+               The observables are not equal to the domain as it does not restrict or
+               truncate the model outside this range. |@docend:pdf.init.obs|
+            extended: |@doc:pdf.init.extended| The overall yield of the PDF.
+               If this is parameter-like, it will be used as the yield,
+               the expected number of events, and the PDF will be extended.
+               An extended PDF has additional functionality, such as the
+               ``ext_*`` methods and the ``counts`` (for binned PDFs). |@docend:pdf.init.extended|
+            norm: |@doc:pdf.init.norm| Normalization of the PDF.
+               By default, this is the same as the default space of the PDF. |@docend:pdf.init.norm|
+            name: |@doc:pdf.init.name| Human-readable name
+               or label of
+               the PDF for better identification.
+               Has no programmatical functional purpose as identification. |@docend:pdf.init.name|
+        """
+        params = {"alpha": alpha, "beta": beta, "gamma": gamma, "n": n}
+        super().__init__(obs=obs, params=params, extended=extended, norm=norm)
+
+    def _unnormalized_pdf(self, x):
+        alpha = self.params["alpha"]
+        beta = self.params["beta"]
+        gamma = self.params["gamma"]
+        n = self.params["n"]
+        x = z.unstack_x(x)
+        return erfexp_pdf_func(x=x, alpha=alpha, beta=beta, gamma=gamma, n=n)
diff --git a/zfit_physics/pdf.py b/zfit_physics/pdf.py
index 3791657..c88e44a 100644
--- a/zfit_physics/pdf.py
+++ b/zfit_physics/pdf.py
@@ -1,6 +1,7 @@
 from .models.pdf_argus import Argus
 from .models.pdf_cmsshape import CMSShape
 from .models.pdf_cruijff import Cruijff
+from .models.pdf_erfexp import ErfExp
 from .models.pdf_relbw import RelativisticBreitWigner
 
-__all__ = ["Argus", "RelativisticBreitWigner", "CMSShape", "Cruijff"]
+__all__ = ["Argus", "RelativisticBreitWigner", "CMSShape", "Cruijff", "ErfExp"]

From 3dc9b940cef4bda9f9bba6c9faf488295f66c620 Mon Sep 17 00:00:00 2001
From: iasonkrom <iason.krom@gmail.com>
Date: Fri, 22 Mar 2024 17:02:01 +0100
Subject: [PATCH 02/10] update changelog

---
 CHANGELOG.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 0597d43..f99edac 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -9,6 +9,7 @@ Major Features and Improvements
 -------------------------------
 - added CMSShape PDF
 - added Cruijff PDF
+- added ErExp PDF
 
 Breaking changes
 ------------------

From 909fa4d05628e7fe7153e3b5d15d9aafcb0b732d Mon Sep 17 00:00:00 2001
From: iasonkrom <iason.krom@gmail.com>
Date: Fri, 22 Mar 2024 17:03:48 +0100
Subject: [PATCH 03/10] accidentally removed end cases

---
 zfit_physics/models/pdf_cruijff.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/zfit_physics/models/pdf_cruijff.py b/zfit_physics/models/pdf_cruijff.py
index 4a85977..4619a47 100644
--- a/zfit_physics/models/pdf_cruijff.py
+++ b/zfit_physics/models/pdf_cruijff.py
@@ -64,6 +64,7 @@ def __init__(
             f(x; \\mu, \\sigma_{L}, \\alpha_{L}, \\sigma_{R}, \\alpha_{R}) = \\begin{cases}
             \\exp{\\left(- \\frac{(x-\\mu)^2}{2 \\sigma_{L}^2 + \\alpha_{L} (x-\\mu)^2}\\right)}, \\mbox{for} x \\leqslant mu \\newline
             \\exp{\\left(- \\frac{(x-\\mu)^2}{2 \\sigma_{R}^2 + \\alpha_{R} (x-\\mu)^2}\\right)}, \\mbox{for} x > mu
+            \\end{cases}
 
         Args:
             mu: Mean value

From f95b2b1dbe965699b4ec704c1adfc1adbf6db848 Mon Sep 17 00:00:00 2001
From: iasonkrom <iason.krom@gmail.com>
Date: Fri, 22 Mar 2024 18:40:17 +0100
Subject: [PATCH 04/10] add integral comment

---
 zfit_physics/models/pdf_erfexp.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/zfit_physics/models/pdf_erfexp.py b/zfit_physics/models/pdf_erfexp.py
index 6917a6c..3dc1a80 100644
--- a/zfit_physics/models/pdf_erfexp.py
+++ b/zfit_physics/models/pdf_erfexp.py
@@ -28,6 +28,8 @@ def erfexp_pdf_func(x, alpha, beta, gamma, n):
     return tf.math.erfc((x - alpha) * beta) * znp.exp(-gamma * (znp.power(x, n) - znp.power(alpha, n)))
 
 
+# Note: There is no analytic integral for the ErfExp PDF
+# We tried with sympy, Mathematica, Wolfram Alpha and https://www.integral-calculator.com/
 class ErfExp(zfit.pdf.BasePDF):
     _N_OBS = 1
 

From 2e40ef81379a383a11565d90cbd67e02b1566dbc Mon Sep 17 00:00:00 2001
From: Iason Krommydas <iason.krom@gmail.com>
Date: Fri, 22 Mar 2024 23:20:36 +0100
Subject: [PATCH 05/10] typo in changelog

Co-authored-by: Jonas Eschle <jonas.eschle@cern.ch>
---
 CHANGELOG.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index f99edac..273778c 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -9,7 +9,7 @@ Major Features and Improvements
 -------------------------------
 - added CMSShape PDF
 - added Cruijff PDF
-- added ErExp PDF
+- added ErfExp PDF
 
 Breaking changes
 ------------------

From e30818c9479f7f22c4776a0f48797874938563b0 Mon Sep 17 00:00:00 2001
From: iasonkrom <iason.krom@gmail.com>
Date: Fri, 22 Mar 2024 23:22:17 +0100
Subject: [PATCH 06/10] add commented sympy code

---
 zfit_physics/models/pdf_erfexp.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/zfit_physics/models/pdf_erfexp.py b/zfit_physics/models/pdf_erfexp.py
index 3dc1a80..95067bf 100644
--- a/zfit_physics/models/pdf_erfexp.py
+++ b/zfit_physics/models/pdf_erfexp.py
@@ -30,6 +30,14 @@ def erfexp_pdf_func(x, alpha, beta, gamma, n):
 
 # Note: There is no analytic integral for the ErfExp PDF
 # We tried with sympy, Mathematica, Wolfram Alpha and https://www.integral-calculator.com/
+# import sympy as sp
+#
+# # Define symbols
+# x, alpha, beta, gamma, n = sp.symbols('x alpha beta gamma n', real=True)
+#
+# # Define the function
+# func = sp.erfc((x - alpha) * beta) * sp.exp(-gamma * (x**n - alpha**n))
+# sp.integrate(func, x)
 class ErfExp(zfit.pdf.BasePDF):
     _N_OBS = 1
 

From 0bf4ec7ebbdec1e046b7ae68c6de8590080e573a Mon Sep 17 00:00:00 2001
From: iasonkrom <iason.krom@gmail.com>
Date: Fri, 22 Mar 2024 23:28:04 +0100
Subject: [PATCH 07/10] update tests

---
 tests/test_pdf_erfexp.py | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/tests/test_pdf_erfexp.py b/tests/test_pdf_erfexp.py
index 7a4eab9..307f0ea 100644
--- a/tests/test_pdf_erfexp.py
+++ b/tests/test_pdf_erfexp.py
@@ -29,21 +29,28 @@ def _erfexp_numpy(x, alpha, beta, gamma, n):
 
 def create_erfexp(alpha, beta, gamma, n, limits):
     obs = zfit.Space("obs1", limits=limits)
-    erfexp = zphys.pdf.ErfExp(alpha=alpha, beta=beta, gamma=gamma, n=n, obs=obs, norm=False)
+    erfexp = zphys.pdf.ErfExp(alpha=alpha, beta=beta, gamma=gamma, n=n, obs=obs)
     return erfexp, obs
 
 
 def test_erfexp_pdf():
     # Test PDF here
     erfexp, _ = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
-    assert erfexp.pdf(90.0).numpy().item() == pytest.approx(
+    assert erfexp.pdf(90.0, norm=False).numpy().item() == pytest.approx(
         erfexp_numpy(90.0, alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true), rel=1e-8
     )
     np.testing.assert_allclose(
-        erfexp.pdf(tf.range(50.0, 130, 10_000)),
+        erfexp.pdf(tf.range(50.0, 130, 10_000), norm=False),
         erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),
         rtol=1e-8,
     )
+    np.testing.assert_allclose(
+        erfexp.pdf(tf.range(50.0, 130, 10_000)),
+        erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true)
+        / 71.18838,
+        rtol=1e-8,
+        atol=1e-8,
+    )
 
 
 def test_erfexp_pdf_random_params():
@@ -56,7 +63,7 @@ def test_erfexp_pdf_random_params():
 
         erfexp, __ = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
         np.testing.assert_allclose(
-            erfexp.pdf(tf.range(50.0, 130, 10_000)),
+            erfexp.pdf(tf.range(50.0, 130, 10_000), norm=False),
             erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),
             rtol=1e-5,
         )

From ca8a1e14f46aa21aa66808f755bb68593f7b187e Mon Sep 17 00:00:00 2001
From: iasonkrom <iason.krom@gmail.com>
Date: Fri, 22 Mar 2024 23:35:05 +0100
Subject: [PATCH 08/10] also update cruijff tests in the same manner

---
 tests/test_pdf_cruijff.py | 31 ++++++++++++++++++++++++++++---
 tests/test_pdf_erfexp.py  |  3 +++
 2 files changed, 31 insertions(+), 3 deletions(-)

diff --git a/tests/test_pdf_cruijff.py b/tests/test_pdf_cruijff.py
index 011fbc1..fdcce7a 100644
--- a/tests/test_pdf_cruijff.py
+++ b/tests/test_pdf_cruijff.py
@@ -19,7 +19,7 @@
 
 def create_cruijff(mu, sigmal, alphal, sigmar, alphar, limits):
     obs = zfit.Space("obs1", limits=limits)
-    cruijff = zphys.pdf.Cruijff(mu=mu, sigmal=sigmal, alphal=alphal, sigmar=sigmar, alphar=alphar, obs=obs, norm=False)
+    cruijff = zphys.pdf.Cruijff(mu=mu, sigmal=sigmal, alphal=alphal, sigmar=sigmar, alphar=alphar, obs=obs)
     return cruijff, obs
 
 
@@ -28,7 +28,7 @@ def test_cruijff_pdf():
     cruijff, _ = create_cruijff(
         mu=mu_true, sigmal=sigmal_true, alphal=alphal_true, sigmar=sigmar_true, alphar=alphar_true, limits=(50, 130)
     )
-    assert cruijff.pdf(90.0).numpy() == pytest.approx(
+    assert cruijff.pdf(90.0, norm=False).numpy() == pytest.approx(
         cruijff_numba.density(
             90.0,
             beta_left=alphal_true,
@@ -39,8 +39,20 @@ def test_cruijff_pdf():
         ).item(),
         rel=1e-8,
     )
+    assert cruijff.pdf(90.0).numpy() == pytest.approx(
+        cruijff_numba.density(
+            90.0,
+            beta_left=alphal_true,
+            beta_right=alphar_true,
+            loc=mu_true,
+            scale_left=sigmal_true,
+            scale_right=sigmar_true,
+        ).item()
+        / 67.71494,
+        rel=1e-7,
+    )
     np.testing.assert_allclose(
-        cruijff.pdf(tf.range(50.0, 130, 10_000)),
+        cruijff.pdf(tf.range(50.0, 130, 10_000), norm=False),
         cruijff_numba.density(
             tf.range(50.0, 130, 10_000).numpy(),
             beta_left=alphal_true,
@@ -51,6 +63,19 @@ def test_cruijff_pdf():
         ),
         rtol=1e-8,
     )
+    np.testing.assert_allclose(
+        cruijff.pdf(tf.range(50.0, 130, 10_000)),
+        cruijff_numba.density(
+            tf.range(50.0, 130, 10_000).numpy(),
+            beta_left=alphal_true,
+            beta_right=alphar_true,
+            loc=mu_true,
+            scale_left=sigmal_true,
+            scale_right=sigmar_true,
+        )
+        / 67.71494,
+        rtol=1e-8,
+    )
     assert cruijff.pdf(tf.range(50.0, 130, 10_000)) <= cruijff.pdf(90.0)
 
 
diff --git a/tests/test_pdf_erfexp.py b/tests/test_pdf_erfexp.py
index 307f0ea..5eeada8 100644
--- a/tests/test_pdf_erfexp.py
+++ b/tests/test_pdf_erfexp.py
@@ -39,6 +39,9 @@ def test_erfexp_pdf():
     assert erfexp.pdf(90.0, norm=False).numpy().item() == pytest.approx(
         erfexp_numpy(90.0, alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true), rel=1e-8
     )
+    assert erfexp.pdf(90.0).numpy().item() == pytest.approx(
+        erfexp_numpy(90.0, alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true) / 71.18838, rel=1e-8
+    )
     np.testing.assert_allclose(
         erfexp.pdf(tf.range(50.0, 130, 10_000), norm=False),
         erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),

From ad1be34680a1b302ad49c5501636dfb1d3c10ae3 Mon Sep 17 00:00:00 2001
From: iasonkrom <iason.krom@gmail.com>
Date: Sat, 23 Mar 2024 17:37:21 +0100
Subject: [PATCH 09/10] rename alpha to mu

---
 tests/test_pdf_erfexp.py          | 41 +++++++++++++++----------------
 zfit_physics/models/pdf_erfexp.py | 22 ++++++++---------
 2 files changed, 31 insertions(+), 32 deletions(-)

diff --git a/tests/test_pdf_erfexp.py b/tests/test_pdf_erfexp.py
index 5eeada8..c3b8ef3 100644
--- a/tests/test_pdf_erfexp.py
+++ b/tests/test_pdf_erfexp.py
@@ -9,48 +9,47 @@
 
 import zfit_physics as zphys
 
-alpha_true = 90.0
+mu_true = 90.0
 beta_true = 0.08
 gamma_true = -1.0
 n_true = 0.2
 
-alpha_range = (65.0, 90.0)
+mu_range = (65.0, 90.0)
 gamma_range = (-10, 10)
 beta_range = (0.01, 10)
 n_range = (0.1, 0.5)
 
 
-def _erfexp_numpy(x, alpha, beta, gamma, n):
-    return special.erfc((x - alpha) * beta) * np.exp(-gamma * (np.power(x, n) - np.power(alpha, n)))
+def _erfexp_numpy(x, mu, beta, gamma, n):
+    return special.erfc((x - mu) * beta) * np.exp(-gamma * (np.power(x, n) - np.power(mu, n)))
 
 
-erfexp_numpy = np.vectorize(_erfexp_numpy, excluded=["alpha", "beta", "gamma", "n"])
+erfexp_numpy = np.vectorize(_erfexp_numpy, excluded=["mu", "beta", "gamma", "n"])
 
 
-def create_erfexp(alpha, beta, gamma, n, limits):
+def create_erfexp(mu, beta, gamma, n, limits):
     obs = zfit.Space("obs1", limits=limits)
-    erfexp = zphys.pdf.ErfExp(alpha=alpha, beta=beta, gamma=gamma, n=n, obs=obs)
+    erfexp = zphys.pdf.ErfExp(mu=mu, beta=beta, gamma=gamma, n=n, obs=obs)
     return erfexp, obs
 
 
 def test_erfexp_pdf():
     # Test PDF here
-    erfexp, _ = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
+    erfexp, _ = create_erfexp(mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
     assert erfexp.pdf(90.0, norm=False).numpy().item() == pytest.approx(
-        erfexp_numpy(90.0, alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true), rel=1e-8
+        erfexp_numpy(90.0, mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true), rel=1e-8
     )
     assert erfexp.pdf(90.0).numpy().item() == pytest.approx(
-        erfexp_numpy(90.0, alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true) / 71.18838, rel=1e-8
+        erfexp_numpy(90.0, mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true) / 71.18838, rel=1e-8
     )
     np.testing.assert_allclose(
         erfexp.pdf(tf.range(50.0, 130, 10_000), norm=False),
-        erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),
+        erfexp_numpy(tf.range(50.0, 130, 10_000), mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true),
         rtol=1e-8,
     )
     np.testing.assert_allclose(
         erfexp.pdf(tf.range(50.0, 130, 10_000)),
-        erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true)
-        / 71.18838,
+        erfexp_numpy(tf.range(50.0, 130, 10_000), mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true) / 71.18838,
         rtol=1e-8,
         atol=1e-8,
     )
@@ -59,41 +58,41 @@ def test_erfexp_pdf():
 def test_erfexp_pdf_random_params():
     # Test PDF here in a loop with random parameters
     for _ in range(1000):
-        alpha_true = np.random.uniform(*alpha_range)
+        mu_true = np.random.uniform(*mu_range)
         beta_true = np.random.uniform(*beta_range)
         gamma_true = np.random.uniform(*gamma_range)
         n_true = np.random.uniform(*n_range)
 
-        erfexp, __ = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
+        erfexp, __ = create_erfexp(mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
         np.testing.assert_allclose(
             erfexp.pdf(tf.range(50.0, 130, 10_000), norm=False),
-            erfexp_numpy(tf.range(50.0, 130, 10_000), alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true),
+            erfexp_numpy(tf.range(50.0, 130, 10_000), mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true),
             rtol=1e-5,
         )
 
 
 def test_erfexp_integral():
     # Test CDF and integral here
-    erfexp, obs = create_erfexp(alpha=alpha_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
+    erfexp, obs = create_erfexp(mu=mu_true, beta=beta_true, gamma=gamma_true, n=n_true, limits=(50, 130))
     full_interval_numeric = erfexp.numeric_integrate(obs, norm=False).numpy()
     true_integral = 71.18838
-    numpy_full_integral = integrate.quad(erfexp_numpy, 50, 130, args=(alpha_true, beta_true, gamma_true, n_true))[0]
+    numpy_full_integral = integrate.quad(erfexp_numpy, 50, 130, args=(mu_true, beta_true, gamma_true, n_true))[0]
     assert full_interval_numeric == pytest.approx(true_integral, 1e-7)
     assert full_interval_numeric == pytest.approx(numpy_full_integral, 1e-7)
 
     numeric_integral = erfexp.numeric_integrate(limits=(80, 100), norm=False).numpy()
-    numpy_integral = integrate.quad(erfexp_numpy, 80, 100, args=(alpha_true, beta_true, gamma_true, n_true))[0]
+    numpy_integral = integrate.quad(erfexp_numpy, 80, 100, args=(mu_true, beta_true, gamma_true, n_true))[0]
     assert numeric_integral == pytest.approx(numpy_integral, 1e-7)
 
 
 # register the pdf here and provide sets of working parameter configurations
 def erfexp_params_factory():
-    alpha = zfit.Parameter("alpha", alpha_true)
+    mu = zfit.Parameter("mu", mu_true)
     beta = zfit.Parameter("beta", beta_true)
     gamma = zfit.Parameter("gamma", gamma_true)
     n = zfit.Parameter("n", n_true)
 
-    return {"alpha": alpha, "beta": beta, "gamma": gamma, "n": n}
+    return {"mu": mu, "beta": beta, "gamma": gamma, "n": n}
 
 
 tester.register_pdf(pdf_class=zphys.pdf.ErfExp, params_factories=erfexp_params_factory)
diff --git a/zfit_physics/models/pdf_erfexp.py b/zfit_physics/models/pdf_erfexp.py
index 95067bf..d21b859 100644
--- a/zfit_physics/models/pdf_erfexp.py
+++ b/zfit_physics/models/pdf_erfexp.py
@@ -8,12 +8,12 @@
 
 
 @z.function(wraps="tensor")
-def erfexp_pdf_func(x, alpha, beta, gamma, n):
+def erfexp_pdf_func(x, mu, beta, gamma, n):
     """Calculate the ErfExp PDF.
 
     Args:
         x: value(s) for which the PDF will be calculated.
-        alpha: Location parameter.
+        mu: Location parameter.
         beta: Scale parameter.
         gamma: Shape parameter.
         n: Shape parameter.
@@ -25,7 +25,7 @@ def erfexp_pdf_func(x, alpha, beta, gamma, n):
         Implementation from https://gitlab.cern.ch/cms-muonPOG/spark_tnp/-/blob/Spark3/RooErfExp.cc
         The parameters beta and gamma are given in reverse order in this c++ implementation.
     """
-    return tf.math.erfc((x - alpha) * beta) * znp.exp(-gamma * (znp.power(x, n) - znp.power(alpha, n)))
+    return tf.math.erfc((x - mu) * beta) * znp.exp(-gamma * (znp.power(x, n) - znp.power(mu, n)))
 
 
 # Note: There is no analytic integral for the ErfExp PDF
@@ -33,17 +33,17 @@ def erfexp_pdf_func(x, alpha, beta, gamma, n):
 # import sympy as sp
 #
 # # Define symbols
-# x, alpha, beta, gamma, n = sp.symbols('x alpha beta gamma n', real=True)
+# x, mu, beta, gamma, n = sp.symbols('x mu beta gamma n', real=True)
 #
 # # Define the function
-# func = sp.erfc((x - alpha) * beta) * sp.exp(-gamma * (x**n - alpha**n))
+# func = sp.erfc((x - mu) * beta) * sp.exp(-gamma * (x**n - mu**n))
 # sp.integrate(func, x)
 class ErfExp(zfit.pdf.BasePDF):
     _N_OBS = 1
 
     def __init__(
         self,
-        alpha: ztyping.ParamTypeInput,
+        mu: ztyping.ParamTypeInput,
         beta: ztyping.ParamTypeInput,
         gamma: ztyping.ParamTypeInput,
         n: ztyping.ParamTypeInput,
@@ -59,10 +59,10 @@ def __init__(
 
         .. math:
 
-            f(x; \\alpha, \\beta, \\gamma, n) = \\text{erfc}(\\beta (x - \\alpha)) \\exp{(-\\gamma (x^n - \\alpha^n))}
+            f(x; \\mu, \\beta, \\gamma, n) = \\text{erfc}(\\beta (x - \\mu)) \\exp{(-\\gamma (x^n - \\mu^n))}
 
         Args:
-            alpha: Location parameter.
+            mu: Location parameter.
             beta: Scale parameter.
             gamma: Shape parameter.
             n: Shape parameter.
@@ -87,13 +87,13 @@ def __init__(
                the PDF for better identification.
                Has no programmatical functional purpose as identification. |@docend:pdf.init.name|
         """
-        params = {"alpha": alpha, "beta": beta, "gamma": gamma, "n": n}
+        params = {"mu": mu, "beta": beta, "gamma": gamma, "n": n}
         super().__init__(obs=obs, params=params, extended=extended, norm=norm)
 
     def _unnormalized_pdf(self, x):
-        alpha = self.params["alpha"]
+        mu = self.params["mu"]
         beta = self.params["beta"]
         gamma = self.params["gamma"]
         n = self.params["n"]
         x = z.unstack_x(x)
-        return erfexp_pdf_func(x=x, alpha=alpha, beta=beta, gamma=gamma, n=n)
+        return erfexp_pdf_func(x=x, mu=mu, beta=beta, gamma=gamma, n=n)

From ef77b134089ae20267201c68f36de568f514c9fc Mon Sep 17 00:00:00 2001
From: Jonas Eschle <jonas.eschle@cern.ch>
Date: Sat, 23 Mar 2024 16:09:24 -0400
Subject: [PATCH 10/10] Update pdf_erfexp.py

---
 zfit_physics/models/pdf_erfexp.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/zfit_physics/models/pdf_erfexp.py b/zfit_physics/models/pdf_erfexp.py
index d21b859..c948dbd 100644
--- a/zfit_physics/models/pdf_erfexp.py
+++ b/zfit_physics/models/pdf_erfexp.py
@@ -55,7 +55,7 @@ def __init__(
     ):
         """ErfExp PDF, the product of a complementary error function and an exponential function.
 
-        Implementation from https://gitlab.cern.ch/cms-muonPOG/spark_tnp/-/blob/Spark3/RooErfExp.cc
+        Implementation following closely `C++ version of custom RooErfExp.cc <https://gitlab.cern.ch/cms-muonPOG/spark_tnp/-/blob/Spark3/RooErfExp.cc>`_
 
         .. math:
 
@@ -64,8 +64,8 @@ def __init__(
         Args:
             mu: Location parameter.
             beta: Scale parameter.
-            gamma: Shape parameter.
-            n: Shape parameter.
+            gamma: Shape parameter, scale of exponential term.
+            n: Shape parameter, power in exponential term.
             obs: |@doc:pdf.init.obs| Observables of the
                model. This will be used as the default space of the PDF and,
                if not given explicitly, as the normalization range.