'''
Created on Apr 25, 2015

@author: mberger
'''
import random
import ROOT
from math import log,sqrt,exp
from functions import thisIsTheEnd

n=10000

data=[]
hist=ROOT.TH1D("hist","gaus",100,0,0)
hist2=ROOT.TH1D("hist2","scaled gaus",100,0,0)
hist2.SetLineColor(2)
b='w'
graph=ROOT.TGraph()
graph3=ROOT.TGraph()
graph4=ROOT.TGraph()
graph5=ROOT.TGraph()

c1=ROOT.TCanvas("sigma","sigma")
c3=ROOT.TCanvas("mean","mean")
c4=ROOT.TCanvas("sigma_mean","sigma mean")
c5=ROOT.TCanvas("chi2Test","chi2 Test")
j=-1
rms0=2
width=0.6

'''
while width<3:
    j+=1     
    hist.Reset('M')
    llsum=0
    data=[]
    mean=0
    mtest=random.uniform(-2,2)
    for i in range(n):
        x=random.gauss(mtest,width)
        data.append(x)
        mean+=x

    mean/=n
    rms=0
    for x in data:
        rms+=(x-mean)**2
        llsum+=(x-0)**2 
    rms/=n
    rms=sqrt(rms)


    ll=n/2.+n*(log(rms)-log(rms0))-n/2*(rms**2/rms0**2)
    ll*=-1

    graph.SetPoint(j,width,ll)
    
    print 'log likelihood:',ll, 'width=',width
    width+=0.005
c1.cd()
graph.Draw('AP')
c1.Update()
'''

'''
mean=-2
mean0=0
j=-1
while mean<2:
    j+=1
    data=[]
    mean_hat=0
    wtest=1#random.uniform(0,5)
    for i in range(n):
        x=random.gauss(mean,wtest)
        data.append(x)
        mean_hat+=x
    mean_hat/=n
    
    rms1=0
    rms2=0
    for x in data:
        rms1+=(x-mean_hat)**2
        rms2+=(x-mean0)**2
#    rms/=n
#    rms=sqrt(rms)

    ll=-n/2.*(log(rms1)-log(rms2))
    mean+=0.005
    print 'll ',ll,'mean',mean
    graph3.SetPoint(j,mean,ll)
    
c3.cd()
graph3.Draw('AP')
c3.Update()

data2=[]
n=10000
for i in range(n):
    data2.append(random.gauss(0,2))
    
graph2=ROOT.TGraph()
width=1.8
counter=-1
while width<2.2:
    counter+=1
    llsum=0
    for x in data2:
        llsum+=(x-0)**2
    ll=0.5*n*log(2*3.141592654) + 0.5*n*log(width**2) + 0.5/(width**2)*llsum
    graph2.SetPoint(counter,width,ll)
    width+=0.005
    
c2=ROOT.TCanvas()
graph2.Draw("AP")
'''


mean=0
mean0=0
rms0=1
wtest=0.1
j=-1
hist4=ROOT.TH1D("hist4","hist4",200,-10,10)
hist5=ROOT.TH1D("hist5","hist5",200,-10,10)
hist5.SetLineColor(2)
#while mean<=1:
while wtest<=2:
    j+=1
    data=[]
    mean_hat=0
    
    llsum=0
    hist4.Reset("M")
    for i in range(n):
        #if i%10==0:
        #    x=random.gauss(mean,wtest*4)
        #else:
        x=random.gauss(mean,wtest)

        data.append(x)
        mean_hat+=x
        hist4.Fill(x)
    mean_hat/=n
    
    rms1=0
    rms2=0
    neff=0
    
    for x in data:
        if abs(x)<3:
            rms1+=(x-mean_hat)**2
            rms2+=(x-mean0)**2
            neff+=1
        
    ll1=-neff/2.*(log(rms1)-log(rms2))
    rms1/=neff
    rms1=sqrt(rms1)
    ll2=neff/2.+neff*(log(rms1)-log(rms0))-neff/2*(rms1**2/rms0**2)
    ll2*=-1
    ll=ll1+ll2
   #llerr1=n/2.*abs(2*resSum/rms1 * rms1)
    #mean+=0.005
    wtest+=0.005
    #graph4.SetPoint(j,mean,ll)
    graph4.SetPoint(j,wtest,ll)
    
    meanBin=hist4.GetXaxis().FindBin(hist4.GetMean());
    height=max(hist4.GetBinContent(meanBin-1),hist4.GetBinContent(meanBin),hist4.GetBinContent(meanBin+1))
    chi2Test=0
    #c5.cd()
    #hist4.Draw()
    #c5.Update()
    for ibi in range(hist4.GetNbinsX()):
        x=hist4.GetBinCenter(ibi)
        expected=height*exp(- ((x-mean0)**2)/2. )
        measured=hist4.GetBinContent(ibi)
        #if expected>0:
        if abs(x)<3:
            chi2Test+= ((expected-measured)**2)/expected
            hist5.SetBinContent(ibi,expected)
     #       print 'chi2testaa',((expected-measured)**2)/expected
     #       print 'expected:',expected, 'measured:',measured,'ibi:',ibi, 'X=',x
            
    print 'll ',ll, 'chi2Test ',chi2Test,'mean',mean,'sig',wtest
    #graph5.SetPoint(j,mean,chi2Test)
    graph5.SetPoint(j,wtest,chi2Test)
    #hist5.Draw("same")
    #c5.Update()
    #u=raw_input()
c4.cd()
graph4.Draw('AP')
c4.Update()

c5.cd()
graph5.Draw('AP')
c5.Update()

'''
data2=[]
n=10000
for i in range(n):
    data2.append(random.gauss(0,2))
    
graph2=ROOT.TGraph()
width=1.8
counter=-1
while width<2.2:
    counter+=1
    llsum=0
    for x in data2:
        llsum+=(x-0)**2
    ll=0.5*n*log(2*3.141592654) + 0.5*n*log(width**2) + 0.5/(width**2)*llsum
    graph2.SetPoint(counter,width,ll)
    width+=0.005
    
c2=ROOT.TCanvas()
graph2.Draw("AP")
'''


thisIsTheEnd()