import ROOT

ROOT.gROOT.ProcessLine(".x rootlogon.C") 
#rfile=ROOT.TFile('outfiles_e12/MC/sphere/drift_86.0/new_error/Cosmics_FopiMC_86ar85_1_scaled_smallCorr2_recl.reco.root','read')
#rfile=ROOT.TFile('outfiles/MC/sphere/drift_86.0/new_error/Cosmics_FopiMC_86ar85_1_test2_recl.reco.root','read')
rfile=ROOT.TFile('outfiles/Data/Cosmics/run_3888_planeErr_scaled.reco.root','read')

tree=rfile.Get('cbmsim')
tree.SetBranchStatus('*',0)
tree.SetBranchStatus('CutCosmics.*',1)
tree.SetBranchStatus('TpcSPHit.*',1)
tree.SetBranchStatus('TpcCluster.*',1)
tree.SetBranchStatus('TpcTrackPreFit.*',1)

ev=-1
numev=500

c1=ROOT.TCanvas("c1","c1",2000,2000)
c1.Divide(3,3)
dtheta=ROOT.TH1D("dtheta","ThetaAxis-ThetaPlane",360,-180,180)
dtheta2=ROOT.TH1D("dtheta2","ThetaMom-ThetaPlane",360,-180,180)
dphi=ROOT.TH1D("dphi","PhiAxis-PhiPlane",720,-360,360)
dphi2=ROOT.TH1D("dphi2","PhiMom-PhiPlane",720,-360,360)
dangle=ROOT.TH1D("dangle","AngleAxis-AnglePlane",720,-360,360)
dangle2=ROOT.TH1D("dangle2","AngleMom-AnglePlane",720,-360,360)
scalprod=ROOT.TH1D('scalprod','scalar product Axis*Plane',100,-1,1)
scalprod2=ROOT.TH1D('scalprod2','scalar product Mom*Plane',100,-1,1)
for e in tree:
    ev+=1
    if ev>numev:
        break

    if ev%10==0:
        print ev
        
    for tr in e.CutCosmics:
        if e.CutCosmics.GetEntries()>1:
            #continue
            pass
        hitids=tr.GetHitIDs()

        for icl in range(len(hitids)):
            sphit=e.TpcSPHit.At(hitids[icl])
            
            cl=e.TpcCluster.At(sphit.getClusterID())
            if cl.index()!=icl:
                print "clindex:", cl.index(),icl
            midplane=cl.midPlane()
            midplaneNorm=midplane.getNormal()
            #midplaneNorm=tr.GetMidPlaneW(icl)
            axis=tr.GetAxis1(icl)
            #axis.SetMag(1)
            #axis=cl.axis(0)
            refposMom=tr.GetProjectionPointsMom(icl)

            #if midplane.getO().Z==-10:
            #    print 'midplane bad'
            #    continue
            if tr.GetIsGeomErr(icl).Mag()!=1:
                #print "geom err",(axis.Theta()-midplane.getNormal().Theta())*ROOT.TMath.RadToDeg(),(refposMom.Theta()-midplane.getNormal().Theta())*ROOT.TMath.RadToDeg()
                continue
                #pass
            if not tr.GetClusterFitGood(icl):
                #print "bad fit",(axis.Theta()-midplane.getNormal().Theta())*ROOT.TMath.RadToDeg(),(refposMom.Theta()-midplane.getNormal().Theta())*ROOT.TMath.RadToDeg()
                continue
                #pass
            if(abs(axis.Theta()-midplaneNorm.Theta())*ROOT.TMath.RadToDeg()-90)>100:
                   print ev,icl

            if(axis.Angle(midplaneNorm)*ROOT.TMath.RadToDeg()>170):
                #continue
                pass

            dtheta.Fill( (axis.Theta()-midplaneNorm.Theta())*ROOT.TMath.RadToDeg())
            dtheta2.Fill( (refposMom.Theta()-midplaneNorm.Theta())*ROOT.TMath.RadToDeg())

            dphi.Fill( (axis.Phi()-midplaneNorm.Phi())*ROOT.TMath.RadToDeg())
            dphi2.Fill( (refposMom.Phi()-midplaneNorm.Phi())*ROOT.TMath.RadToDeg())

            dangle.Fill(axis.Angle(midplaneNorm)*ROOT.TMath.RadToDeg())
            dangle2.Fill(refposMom.Angle(midplaneNorm)*ROOT.TMath.RadToDeg())

            axis.SetMag(1)
            normal=midplaneNorm
            normal.SetMag(1)
            refposMom.SetMag(1)
            scalprod.Fill(axis*normal)
            scalprod2.Fill(refposMom*normal)


c1.cd(1)
c1.cd(1).SetLogy()
dtheta.Draw()
c1.cd(2)
c1.cd(2).SetLogy()
dtheta2.Draw()
c1.cd(3)
c1.cd(3).SetLogy()
dphi.Draw()
c1.cd(4)
c1.cd(4).SetLogy()
dphi2.Draw()
c1.cd(5)
c1.cd(5).SetLogy()
dangle.Draw()
c1.cd(6)
c1.cd(6).SetLogy()
dangle2.Draw()
c1.cd(7)
c1.cd(7).SetLogy()
scalprod.Draw()
c1.cd(8)
c1.cd(8).SetLogy()
scalprod2.Draw()
c1.Update()
u=raw_input("done?")