In [5]:
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.stats import norm
In [6]:
# Cruijff function
# http://rarfit.sourceforge.net/html/RooCruijff_8cc-source.html
# 
# It is a bifurcated Gaussian (different and high/low sigmas)
# with sigma^2 replaced by sigma_R/L^2 + alpha_R/L*(x-mu)^2
# 
# Note that the integral over the full range is unbounded, 
# so cannot be normalized.  We normalize it here arbitrarily
# so that the peak is the same as the peak of a Gaussian
# of sigma=(sigmaL+sigmaR)/2.
#
def cruijff(x, mu, sigmaL, sigmaR, alphaL, alphaR):
    sigmaAv = (sigmaL+sigmaR)/2.
    N       = 1./(math.sqrt(2*math.pi)*sigmaAv)
    dx      = x - mu
    if hasattr(x, "__len__"):  # x is a vector, assumed numpy
        fL   = 2*sigmaL*sigmaL + alphaL*dx*dx
        fR   = 2*sigmaR*sigmaR + alphaR*dx*dx
        LorR = x < mu
        f    = LorR * fL + (1-LorR) * fR
        return N*np.exp(-dx*dx/f)
    
    else:   # x is a number
        if x<mu:
            f = 2*sigmaL*sigmaL + alphaL*dx*dx
        else:
            f = 2*sigmaR*sigmaR + alphaR*dx*dx
        return N*math.exp(-dx*dx)
In [9]:
mean   = 100
sigmaL = 7
sigmaR = 7
alphaL = 0.2
alphaR = 0.3
In [10]:
x  = np.linspace(max(0,mean-4*sigmaL), mean+4*sigmaR, 1000) # plot range
ax = plt.subplot(111)
ax.plot(x, cruijff(x, mean, sigmaL, sigmaR, alphaL, alphaR),  label="Cruijff")
ax.plot(x, norm.pdf(x, scale=0.5*(sigmaL+sigmaR), loc=mean),  label="Gaussian")
ax.grid()
ax.legend()
_ = ax.set_ylim(bottom=0)
In [ ]: