import ROOT,math

graph=ROOT.TGraph()
graph_err=ROOT.TGraphErrors()



c1=ROOT.TCanvas('c1','c1',1000,700)
c2=[]
c3=ROOT.TCanvas('c3','c3')

func=ROOT.TF1("f","gaus")
func.SetLineColor(1)
func2=ROOT.TF1("f2","gaus")
func2.SetLineColor(3)

draw_single=False
#draw_single=True
hmeanErr=[]
hmeanErrnorm=[]
hsigma=[]
meanErrProf=[]
fitfunc3=[]
for ii in range(20):
    sigma=2*(ii+1)
    c2.append(ROOT.TCanvas('c2_{0}'.format(ii),'c2 {0}'.format(sigma),1000,1000))
    c2[-1].Divide(1,3)

for ii in range(20):
    sigma=2*(ii+1)
    mean=0
    const=100

    hmeanErr.append(ROOT.TH2D("hmeanErr{0}".format(ii),
                              "error of mean as function of measurement error (sigma: {0})".format(sigma),
                              200,0,ii+1,300,0,3))
    hmeanErr[-1].GetXaxis().SetTitle("measurement error")
    hmeanErr[-1].GetYaxis().SetTitle("fit error on mean")

    hmeanErrnorm.append(ROOT.TH2D("hmeanErrnorm{0}".format(ii),
                                  "(error of mean)/sigma as function of measurement error (sigma:{0})".format(sigma),
                                  (ii+1)*10,0,ii+1,300,0,3./20.))
    hmeanErrnorm[-1].GetXaxis().SetTitle("measurement error")
    hmeanErrnorm[-1].GetYaxis().SetTitle("fit error on mean / sigma")
    
    hsigma.append(ROOT.TH2D("hsigma{0}".format(ii),"sigma as function of measurement error (sigma:{0})".format(sigma),(ii+1)*10,0,(ii+1),1000,sigma-5,sigma+5))
    hsigma[-1].GetXaxis().SetTitle("measurement error")
    hsigma[-1].GetYaxis().SetTitle("sigma from fit")
    
    for i in range(100000):
        testgraph=ROOT.TGraph()
        testgraph.SetMarkerStyle(20)

        if i%1000==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)


        func.SetParameter(0,const)
        func.SetParameter(1,mean)
        func.SetParameter(2,sigma)

        func2.SetParameter(0,const)
        func2.SetParameter(1,mean)
        func2.SetParameter(2,sigma)

        #histo=ROOT.TH1D("h","h",50,-5*sigma,5*sigma)

        #N=int(ROOT.gRandom.Uniform(3,40))
        N=20

        error=ROOT.gRandom.Uniform(0,sigma/2.)
        #print 'the error is:',error

        posval=[]
        for j in range(0,N):
            pos=(ROOT.gRandom.Uniform(-3*sigma,3*sigma))
            #print 'pos=',pos
            pos=int(pos/error)*error
            #print 'pos binned=',pos
            if pos in (x[0] for x in posval):
                continue
            val=func.Eval(pos)
            posval.append([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)

        if draw_single:
            c1.cd()
            testgraph.Draw('AP')
            testgraph.Fit(func)


        testgraph2=ROOT.TGraph()
        testgraph2.SetMarkerStyle(20)
        testgraph2.SetMarkerColor(3)

        #error=ROOT.gRandom.Uniform(0,10)

        for j in range(len(posval)):
            pos=posval[j][0]+ROOT.gRandom.Uniform(-error/2.,error/2.)
            val=posval[j][1]
            if val>0:
                testgraph2.SetPoint(j,pos,val)
                #testgraph2.SetPointError(j,

        opt='QN'
        if draw_single:
            opt=''
        testgraph2.Fit(func2,opt)

        if draw_single:
            testgraph2.Draw('P')
            c1.Update()
            u=raw_input()

        hmeanErr[ii].Fill(error,func2.GetParError(1))
        hmeanErrnorm[ii].Fill(error,func2.GetParError(1)/func2.GetParameter(2))
        hsigma[ii].Fill(error,func2.GetParameter(2))


    c2[ii].cd(1)
    hmeanErr[ii].Draw('colz')
    meanErrProf.append(hmeanErr[ii].ProfileX(hmeanErr[ii].GetTitle()+"prof",1,-1,"OP"))
    meanErrProf[ii].Draw('same')
    fitfunc3.append(ROOT.TF1('func3_{0}'.format(ii),"[0]+[1]*x",0,(ii+1)))
    fitfunc3[-1].SetParameter(0,0)
    fitfunc3[-1].FixParameter(0,0)
    fitfunc3[-1].SetParameter(1,0.075)
    meanErrProf[ii].Fit(fitfunc3[ii],"QNMW","",0,(ii+1))
    graph_err.SetPoint(ii,sigma,fitfunc3[ii].GetParameter(1))
    #graph_err.SetPointError(ii,0,fitfunc3[ii].GetParError(1))
    c2[ii].cd(2)
    meanErrProf[ii].Draw()
    fitfunc3[ii].Draw('same')
    #hmeanErrnorm[ii].Draw('colz')
    c2[ii].cd(3)
    hsigma[ii].Draw('colz')
    c2[ii].Update()


c3.cd()
graph_err.SetMarkerStyle(7)
graph_err.Draw("AP")
c3.Update()
u='w'
while u!='q':
    u=raw_input('done?')