import ROOT,sys,os
import random

#----------------------------------------------------------------------------------------------------------------------------------
#
#ist=ROOT.TH1D("hist","hist",20,-10,10)
#ist2=ROOT.TH1D("hist2","hist treshold",20,-10,10)
#
#cog1=ROOT.TH1D("hcog1","cog1",   1000,-2,2)
#cog2=ROOT.TH1D("hcog2","cog2",   1000,-2,2)
#mean1=ROOT.TH1D("hmean1","mean1",1000,-2,2)
#mean2=ROOT.TH1D("hmean2","mean2",1000,-2,2)
#=ROOT.TCanvas('c1','c',1000,1000)
#.Divide(2,2)
#
#gaus=ROOT.TF1("fgaus","gaus")
#gaus.FixParameter(2,2)
#verbose=False
#erbose=True
#or j in range(10000):
#   cog1=0
#   amp=0
#   hist.Reset()
#   hist2.Reset()
#   fgaus.SetParameter(0,200)
#   fgaus.SetParameter(1,0)
#   fgaus.SetParameter(2,2)
#
#   vararr=[]
#
#   for i in range(10):
#       var=random.uniform(-100,100)
#       var=int(var/10)
#       while var in vararr:
#           var=random.uniform(-100,100)
#           var=int(var/10)
#       vararr.append(var)
#       amp+=fgaus.Eval(var)
#       hist.Fill(var,fgaus.Eval(var))
#       cog1+=var*fgaus.Eval(var)
#       
#   cog1/=amp
#
#
#   if verbose:
#       c.cd(1)
#       hist.Fit(fgaus,"B")
#   else:
#       hist.Fit(fgaus,"BQN")
#   mean1=fgaus.GetParameter(1)
#   
#   cog2=0
#   amp2=0
#   fgaus.SetParameter(0,200)
#   fgaus.SetParameter(1,0)
#   fgaus.SetParameter(2,2)
#   for i in range(hist.GetNbinsX()):
#       if hist.GetBinContent(i)>20:
#           hist2.SetBinContent(i,hist.GetBinContent(i))
#           cog2+=hist.GetBinCenter(i)*hist.GetBinContent(i)
#           amp2+=hist.GetBinContent(i)
#   if amp2!=0:
#       cog2/=amp2
#       if verbose:
#           c.cd(2)
#           hist2.Fit(fgaus,'B')
#       else:
#           hist2.Fit(fgaus,'BQN')
#       mean2=fgaus.GetParameter(1)
#   else:
#       cog2=-10
#       mean2=-10
#   
#   if j%1000==0:
#       print 'at',j
#
#
#   hcog1.Fill(cog1)
#   hcog2.Fill(cog2)
#   hmean1.Fill(mean1)
#   hmean2.Fill(mean2)
#   
#   if verbose:
#       c.cd(1)
#       hist.Draw()
#       c.cd(2)
#       hist2.Draw()
#       print 'j=',j,'cog1=',cog1, 'mean1=',mean1,'cog2=',cog2,'mean2=',mean2
#       c.Update()
#       u=raw_input('next?')    
#   
#.cd(1)
#cog1.Draw()
#mean1.SetLineColor(2)
#mean1.Draw("same")
#
#.cd(2)
#cog2.SetLineWidth(2)
#cog2.Draw()
#mean2.SetLineColor(2)
#mean2.Draw("same")
#
#diff1=hcog1.Clone("hdiff1")
#diff1.Add(hmean1,-1)
#.cd(3)
#diff1.Draw()
#
#diff2=hcog2.Clone("hdiff2")
#diff2.Add(hmean2,-1)
#.cd(4)
#diff2.Draw()
#
#.Update()
#----------------------------------------------------------------------------------------------------------------------------------

#S=ROOT.TMatrixD(3,3)
#S[0][0]=1
#S[1][1]=2
#S[2][2]=3
#
#D=ROOT.TMatrixD(3,3)
#D[0][0]=1
#D[1][1]=1
#D[2][2]=1
#D_i=ROOT.TMatrixD(3,3,D.GetMatrixArray())
#D.Invert()
#
#
#M1=ROOT.TMatrixD(3,3, (D*S*D_i).GetMatrixArray())
#M1.Print()
#
#eigenprop=ROOT.TMatrixDEigen(M1)
#S2=eigenprop.GetEigenValues()
#D2=eigenprop.GetEigenVectors()
#D2_i=ROOT.TMatrixD(3,3,D2.GetMatrixArray())
#D2_i.Invert()
#
#M2=ROOT.TMatrixD(3,3, (D2_i*S2*D2).GetMatrixArray())
#M2.Print()

h1=ROOT.TH2D('h1','h1',100,-5,5,100,-5,5)
h2=ROOT.TH2D('h2','h2',100,-5,5,100,-5,5)
h3=ROOT.TH2D('h3','h3',100,-5,5,100,-5,5)

for i in range(h1.GetNbinsX()+2):
    for j in range(h1.GetNbinsY()+2):
        if h1.GetXaxis().GetBinCenter(i)<1:
            h1.Fill(h1.GetXaxis().GetBinCenter(i),h1.GetYaxis().GetBinCenter(j),1)
            h1.Fill(h1.GetXaxis().GetBinCenter(i),h1.GetYaxis().GetBinCenter(j),1)
        if h1.GetXaxis().GetBinCenter(i)>-1:
            h2.Fill(h2.GetXaxis().GetBinCenter(i),h2.GetYaxis().GetBinCenter(j),1)
            h2.Fill(h2.GetXaxis().GetBinCenter(i),h2.GetYaxis().GetBinCenter(j),1)
        h3.Fill(h3.GetXaxis().GetBinCenter(i),h3.GetYaxis().GetBinCenter(j),5)
        h3.Fill(h3.GetXaxis().GetBinCenter(i),h3.GetYaxis().GetBinCenter(j),5)

c1=ROOT.TCanvas('c1','c1',1500,1500)
c1.Divide(3,3)
c1.cd(1)
h1.Draw('colz')
c1.cd(2)
h2.Draw('colz')
c1.cd(3)
h3.Draw('colz')

h4=h1.Clone("h4")
h4.Add(h2)
c1.cd(4)
h4.Draw('colz')

h5=h3.Clone('h5')
h5.Divide(h4)
c1.cd(5)
h5.Draw('colz')

file1=ROOT.TFile("testfile1.root",'recreate')
h1.Write()
file1.Close()

file2=ROOT.TFile("testfile2.root",'recreate')
h2.SetName("h1")
h2.SetTitle("h1")
h2.Write()
file2.Close()

os.system("hadd -f testfile.root testfile1.root testfile2.root")

file3=ROOT.TFile("testfile.root","read")
h6=file3.Get("h1")
h6.SetTitle("h6")

c1.cd(6)
h6.Draw('colz')

h7=h3.Clone('h7')
h7.Divide(h6)
c1.cd(7)
h7.Draw('colz')


u=raw_input("done?")