Merge pull request #81 from andreas01060/master

Adding Double Crystal Ball function

doc/api.rst

@@ -173,6 +173,13 @@ crystalball
 .. plot:: pyplots/pdf/crystalball.py
     :class: lightbox
 
+doublecrystalball
+^^^^^^^^^^^
+.. autofunction:: doublecrystalball
+
+.. plot:: pyplots/pdf/doublecrystalball.py
+    :class: lightbox
+
 cruijff
 ^^^^^^^
 .. autofunction:: cruijff

doc/index.rst

@@ -100,6 +100,7 @@ and send us pull request.
 .. autosummary::
     gaussian
     crystalball
+    doublecrystalball
     cruijff
     cauchy
     rtv_breitwigner

doc/pyplots/pdf/doublecrystalball.py

@@ -0,0 +1,14 @@
+from probfit.pdf import doublecrystalball
+from probfit.plotting import draw_normed_pdf
+import matplotlib.pyplot as plt
+
+bound = (-2, 12)
+
+arg = dict(alpha=1.,alpha2=2., n=2.,n2=4, mean=5, sigma=1)
+draw_normed_pdf(crystalball, arg=arg, bound=bound, label=str(arg), density=True)
+
+arg = dict(alpha=2, alpha2=1, n=7., n2=10., mean=5, sigma=1)
+draw_normed_pdf(crystalball, arg=arg, bound=bound, label=str(arg), density=True)
+
+plt.grid(True)
+plt.legend().get_frame().set_alpha(0.5)

probfit/__init__.py

@@ -5,9 +5,10 @@
 * Docs: http://probfit.readthedocs.io
 """
 
+
 from .costfunc import UnbinnedLH, BinnedLH, Chi2Regression, BinnedChi2,\
  SimultaneousFit
-from .pdf import doublegaussian, ugaussian, gaussian, crystalball, \
+from .pdf import doublegaussian, doublecrystalball, ugaussian, gaussian, crystalball, \
     argus, cruijff, linear, poly2, poly3, novosibirsk, \
     Polynomial, HistogramPdf, cauchy, rtv_breitwigner, johnsonSU
 from .toy import gen_toy, gen_toyn
@@ -41,6 +42,7 @@
     'cauchy',
     'rtv_breitwigner',
     'crystalball',
+    'doublecrystalball',
     'describe',
     'doublegaussian',
     'johnsonSU',

probfit/pdf.pxd

@@ -21,3 +21,5 @@ cpdef double novosibirsk(double x, double width, double peak, double tail)
 cpdef double rtv_breitwigner(double x, double m, double gamma)
 
 cpdef double cauchy(double x, double m, double gamma)
+
+cpdef double doublecrystalball(double x, double alpha, double alpha2, double n, double n2, double mean, double sigma)

probfit/pdf.pyx

@@ -312,6 +312,43 @@ cpdef double crystalball(double x, double alpha, double n, double mean, double s
             ret = A*pow(B-d,-n)
     return ret
 
+cpdef double doublecrystalball(double x, double alpha, double alpha2, double n, double n2, double mean, double sigma):
+    """
+    Unnormalized double crystal ball function
+    A gaussian core with two power tails
+
+    
+
+    """
+    cdef double d = 0.
+    cdef double ret = 0
+    cdef double A = 0
+    cdef double B = 0
+    cdef double A2 = 0
+    cdef double B2 = 0
+    if sigma < smallestdiv:
+        ret = badvalue
+    elif fabs(alpha) < smallestdiv:
+        ret = badvalue
+    elif n<1.:
+        ret = badvalue
+    else:
+        d = (x-mean)/sigma
+        if d<-alpha:
+            al = fabs(alpha)
+            A=pow(n/al,n)*exp(-al**2/2.)
+            B=n/al-al
+            ret = A*pow(B-d,-n)
+
+        elif d < alpha2 :
+            ret = exp(-0.5*d**2)
+        else:
+            al2 = fabs(alpha2)
+            A2=pow(n2/al2,n2)*exp(-al2**2/2.)
+            B2=n2/al2-al2
+            ret = A2*pow(B2+d,-n2)
+    return ret
+
 
 #Background stuff
 cpdef double argus(double x, double c, double chi, double p):

tests/test_func.py

@@ -66,7 +66,22 @@ def test_crystalball():
     assert_allclose(pdf.crystalball(14, 1, 2, 10, 2), 0.1353352832366127)
     assert_allclose(pdf.crystalball(6, 1, 2, 10, 2), 0.26956918209450376)
 
-
+# cpdef double doubecrystalball(double x,double alpha,double alpha2, double n,double n2, double mean,double sigma)
+def test_doublecrystalball():
+    assert describe(pdf.doublecrystalball) == ['x', 'alpha', 'alpha2', 'n', 'n2', 'mean', 'sigma']
+    assert_allclose(pdf.doublecrystalball(10, 1, 1, 2, 2, 10, 2), 1.)
+    assert_allclose(pdf.doublecrystalball(11, 1, 1, 2, 2, 10, 2), 0.8824969025845955)
+    assert_allclose(pdf.doublecrystalball(12, 1, 1, 2, 2, 10, 2), 0.6065306597126334)
+    assert_allclose(pdf.doublecrystalball(14, 1, 1, 2, 2, 10, 2), 0.26956918209450376)
+    assert_allclose(pdf.doublecrystalball(6, 1, 1, 2, 2, 10, 2), 0.26956918209450376)
+    assert_allclose(pdf.doublecrystalball(-10, 1, 5, 3, 4, 10, 2), 0.00947704155801)
+    assert_allclose(pdf.doublecrystalball(0, 1, 5, 3, 4, 10, 2), 0.047744395954055)
+    assert_allclose(pdf.doublecrystalball(11, 1, 5, 3, 4, 10, 2), 0.8824969025846)
+    assert_allclose(pdf.doublecrystalball(20, 1, 5, 3, 4, 10, 2), 0.0000037266531720786)
+    assert_allclose(pdf.doublecrystalball(25, 1, 5, 3, 4, 10, 2), 0.00000001287132228271)
+
+
+
 # cpdef double argus(double x, double c, double chi, double p)
 def test_argus():
     assert describe(pdf.argus) == ['x', 'c', 'chi', 'p']

tests/test_plotting.py

@@ -58,6 +58,7 @@ def test_draw_pdf_linear():
     f = linear
     draw_pdf(f, {'m': 1., 'c': 2.}, bound=(-10, 10))
 
+
 # There is a slight difference in the x-axis tick label positioning for this
 # plot between Python 2 and 3, it's not important here so increase the RMS
 # slightly such that it's ignored

tests/test_toy.py

@@ -16,7 +16,6 @@ def test_gen_toy():
     ntoy = 100000
     toy = gen_toy(crystalball, ntoy, bound=bound,
                   alpha=1., n=2., mean=1., sigma=0.3, quiet=False)
-
     assert len(toy) == ntoy
 
     htoy, bins = np.histogram(toy, bins=1000, range=bound)