import numpy as np
import pylab as pl

pl.ion()
def gauss(x,x0=0,sig=10,A=1.):
    return A/np.sqrt(2.*np.pi*sig**2)*np.exp(-(x-x0)**2/(2.*sig**2))

x0 = 0
sig = 10
A = 3

x = np.arange(-50.,51,1.)

def plot_gauss(x,x0,sig,A,plot=True):
    y1 = gauss(x,x0 = x0,sig = sig, A = A)
    y2 = np.random.poisson(gauss(x,x0 = x0,sig =sig, A = A))

    A1 = sum(y1)
    A2 = sum(y2)
    if A1!=0:
        x1 = np.average(x,weights=y1)
        sig1 = np.sqrt(np.average((x-x1)**2,weights=y1))
    else:
        x1=0
        sig1=0
    if A2!=0:
        x2 = np.average(x,weights=y2)
        sig2 = np.sqrt(np.average((x-x2)**2,weights=y2))
    else:
        x2=0
        sig2=0
    if plot:
        pl.figure()
        pl.plot(x,y1)
        pl.plot(x,y2)
        print "g :",x0,sig,A
        print "g1 :",x1,sig1,A1
        print "g2 :",x2,sig2,A2
    return (x1,sig1,A1),(x2,sig2,A2)

x1a = []
sig1a = []
A1a = []
x2a = []
sig2a = []
A2a = []

x0 = 0
sig = 10
A = 10

for i in range(10000):
    g1,g2 = plot_gauss(x,x0,sig,A,plot=False)
    x1a.append(g1[0])
    x2a.append(g2[0])
    sig1a.append(g1[1])
    sig2a.append(g2[1])
    A1a.append(g1[2])
    A2a.append(g2[2])


pl.figure()
pl.title("x0")
pl.hist(x2a,bins=np.arange(-30,30.1,0.1))
pl.figure()
pl.title("sig")
pl.hist(sig2a,bins=np.arange(0.,30.1,1.))
pl.hist(np.random.poisson(sig,10000),bins=np.arange(0.,4*A,1.),alpha=0.5)    
pl.hist(np.random.normal(np.mean(np.array(sig2a)),np.std(np.array(sig2a)),10000),bins=np.arange(0.,4*A,1.),alpha=0.5)    
pl.figure()
pl.title("A")
pl.hist(A2a,bins=np.arange(0.,3*A,1.))
pl.hist(np.random.poisson(A,10000),bins=np.arange(0.,4*A,1.),alpha=0.5)