import ROOT


eigenvalgraphs=[]
eigenvalgraphs.append(ROOT.TGraph())
eigenvalgraphs.append(ROOT.TGraph())

testgraph=ROOT.TGraph()
for thetaDeg in range(1,89):
    print "at theta:",thetaDeg
    theta=thetaDeg*ROOT.TMath.DegToRad()    

    thecov=ROOT.TMatrixDSym(3)
    thecov[0][0]=2
    thecov[1][1]=2+2*ROOT.TMath.Cos(theta)#*ROOT.TMath.Sin(theta)
    thecov[2][2]=4
    
    #testgraph.SetPoint(testgraph.GetN(),theta,ROOT.TMath.Sin(theta)*4+ROOT.TMath.Cos(theta)*2)
    
    normal=ROOT.TVector3(1,0,0)
    normal.SetTheta(theta)
    normal.SetPhi(90*ROOT.TMath.DegToRad())
    #normal.SetPhi(theta)
    #normal.SetTheta(45*ROOT.TMath.DegToRad())
    print 'normal vector:'
    normal.Print()
    planeU=normal.Orthogonal()
    planeV=normal.Cross(planeU)
    planeM=ROOT.TMatrixD(3,2)
    planeU.SetMag(1)
    planeV.SetMag(1)
    print 'planeU vector:'
    planeU.Print()
    print 'planeV vector:'
    planeV.Print()
    planeM[0][0]=planeU[0]
    planeM[1][0]=planeU[1]
    planeM[2][0]=planeU[2]
    planeM[0][1]=planeV[0]
    planeM[1][1]=planeV[1]
    planeM[2][1]=planeV[2]
    
    planeUD=ROOT.TVectorD(3)
    for i in range(3):
        planeUD[i]=planeU[i]
    
    planeM.Print()
    print 'thecov:'
    thecov.Print()
    projCov=ROOT.TMatrixDSym(3,thecov.GetMatrixArray())
    #planeM_T=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,planeM)
    #tmp=ROOT.TMatrixD(thecov,ROOT.TMatrixD.kMult,planeM)
    #projCov=ROOT.TMatrixDSym(2,  ROOT.TMatrixD(planeM_T,ROOT.TMatrixD.kMult,tmp).GetMatrixArray())
    print 'projcov before:'
    projCov.Print()
    #projCov.Invert()

    projCov.SimilarityT(planeM)
    #projCov.Invert()
    print 'projcov after:'
    projCov.Print()
    
    #thecov.Invert()
    projCovU=1./thecov.Similarity(planeUD)
    #thecov.Invert()
    testgraph.SetPoint(testgraph.GetN(),theta,projCovU)
                       
    eigenprob=ROOT.TMatrixDSymEigen(projCov)
    eigenval=eigenprob.GetEigenValues()
    eigenvect=eigenprob.GetEigenVectors()
    print 'projCov eigenvalues:'
    eigenval.Print()
    print 'projcov eigenvectors:'
    eigenvect.Print()
    tmp1=ROOT.TVector3(planeU)
    tmp1*=eigenvect[0][0]
    tmp2=ROOT.TVector3(planeV)
    tmp2*=eigenvect[1][0]
    ev1=ROOT.TVector3(tmp1+tmp2)
    tmp1=ROOT.TVector3(planeU)
    tmp1*=eigenvect[0][1]
    tmp2=ROOT.TVector3(planeV)
    tmp2*=eigenvect[1][1]
    ev2=ROOT.TVector3(tmp1+tmp2)
    #if ev1.Z()<ev2.Z():
    if(1):
        eigenvalgraphs[0].SetPoint(eigenvalgraphs[0].GetN(),theta,eigenval[0])
        eigenvalgraphs[1].SetPoint(eigenvalgraphs[1].GetN(),theta,eigenval[1])
    #else:
    #    eigenvalgraphs[0].SetPoint(eigenvalgraphs[0].GetN(),theta,eigenval[1])
    #    eigenvalgraphs[1].SetPoint(eigenvalgraphs[1].GetN(),theta,eigenval[0])
    
eigenvalgraphs[0].GetYaxis().SetRangeUser(-4.5,4.5)
eigenvalgraphs[0].Draw('APL')
eigenvalgraphs[1].SetMarkerColor(2)
eigenvalgraphs[1].Draw('PL')
testgraph.SetMarkerColor(4)
testgraph.Draw('PL')
u=raw_input('done?')