import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.stats import norm
# 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)
mean = 100
sigmaL = 7
sigmaR = 7
alphaL = 0.2
alphaR = 0.3
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)