import ROOT,math


graph=ROOT.TGraph()
graph2=ROOT.TGraph()
hist_s=ROOT.TH2D("hist_s","sigma",100,0,10,100,0,100)
hist_n=ROOT.TH2D('hist_n','N',100,0,10,40,0,40)
hist_c=ROOT.TH2D('hist_c','C',100,0,10,100,0,1000)
func=ROOT.TF1("f","gaus")
func2=ROOT.TF1("f2","gaus")
c1=ROOT.TCanvas('c1','c1',1000,700)
c2=ROOT.TCanvas('c2','c2',1000,1000)
c3=ROOT.TCanvas('c3','c3',1000,1000)

do_draw=False
#do_draw=True
c1.Divide(1,2)
c1.cd()
for i in range(5000):
    if i>0:
        del histo
        del histo2
        del testgraph
    #histo.Reset()    
    if i%100==0:
        print i

    sigma=ROOT.gRandom.Uniform(1,100)
    mean=ROOT.gRandom.Uniform(-0.5*sigma,0.5*sigma)
    const=ROOT.gRandom.Uniform(1,1000)

    histo=ROOT.TH1D("h","h",50,-5*sigma,5*sigma)
    histo2=ROOT.TH1D("h2","h2",50,-5*sigma,5*sigma)
    testgraph=ROOT.TGraphErrors()

    #N=int(ROOT.gRandom.Uniform(1,1000))
    #for j in range(N):
    #    histo.Fill(ROOT.gRandom.Gaus(0,sigma))

    N2=int(ROOT.gRandom.Uniform(3,40))
    func2.SetParameter(0,const)
    func2.SetParameter(1,mean)
    func2.SetParameter(2,sigma)

    posval=[]
    for j in range(N2):
        pos=(ROOT.gRandom.Uniform(-3*sigma,3*sigma))
        val=func2.Eval(pos)#*ROOT.gRandom.Uniform(0.9,1.1)
        posval.append([pos,val])
        histo2.Fill(pos,val)

    posval.sort(key=lambda x: x[0])
    for j in range(len(posval)):
        pos=posval[j][0]
        val=posval[j][1]
        if val>0:
            testgraph.SetPoint(j,pos,val)
            testgraph.SetPointError(j,20,0.05*val)

    #histo.Scale(1./(math.sqrt(2*ROOT.TMath.Pi())*sigma))
    #histo.Draw()

    N=N2
    if do_draw:
        c1.cd(1)
        testgraph.Draw("AP")
        testgraph.SetMarkerStyle(20)

        c1.cd(2)
        histo2.Draw()
        c1.Update()
        #u=raw_input()
    
    func.SetParameter(0,histo.GetBinContent(histo.GetMaximumBin()))
    func.SetParameter(1,histo.GetMean())
    #func.SetParameter(2,histo.GetRMS())
    func.SetParameter(2,sigma)
    #histo.Fit(func,"Q")
    testgraph.Fit(func,"Q")
    #histo2.Fit(func,"QL")
    if do_draw:
        c1.Update()
        c1.cd(1).Update()
        print func.GetParameter(0),func.GetParameter(1),func.GetParameter(2),
        u=raw_input()

    const_f=func.GetParameter(0)
    sigma_f=func.GetParameter(2)
    mean_err=func.GetParError(1)
    
    #graph.SetPoint(i,sigma,mean_err)
    graph.SetPoint(i,mean_err,sigma_f*const_f)
    hist_s.Fill(mean_err,sigma_f)
    graph2.SetPoint(i,mean_err,N)
    hist_n.Fill(mean_err,N)
    hist_c.Fill(mean_err,const_f)


c2.Divide(1,2)
c2.cd(1)
graph.SetMarkerStyle(20)
graph.SetMarkerSize(0.5)
graph.GetXaxis().SetTitle("#Delta Mean")
graph.GetYaxis().SetTitle("#sigma")
graph.Draw("AP")
c2.cd(2)
graph2.SetMarkerStyle(20)
graph2.SetMarkerSize(0.5)
graph2.GetXaxis().SetTitle("#Delta Mean")
graph2.GetYaxis().SetTitle("N")
graph2.Draw("AP")
c2.Update()



c3.Divide(1,3)
c3.cd(1)
hist_s.GetXaxis().SetTitle("#Delta Mean")
hist_s.GetYaxis().SetTitle("#sigma")
hist_s.Draw("colz")
c3.cd(2)
hist_n.GetXaxis().SetTitle("#Delta Mean")
hist_n.GetYaxis().SetTitle("N")
hist_n.Draw('colz')
c3.cd(3)
hist_c.GetXaxis().SetTitle("#Delta Mean")
hist_c.GetYaxis().SetTitle("Const")
hist_c.Draw("colz")
c3.Update()

u='w'
while u!='q':
    u=raw_input("Done?")