import sys,os
sys.path.append('/home/mberger/fopiroot/fopiroot_dev/macro/tpc/FOPI/mberger')
pandapath=os.environ.get('PANDAPATH')
sys.path.append(pandapath+'/macro/tpc/FOPI/python/argparse-1.2.1')
import glob,math,argparse
from functions import Divideh,set_palette,setRange,getColors,getColors2
from anaFile import anaFile

def writeHist(histos):
    for hl in histos:
        for h in hl:
            h.Write()

def drawAngles(c,h,i,col,opt):
    for j in range(len(h)):
        c.cd(j+1)
        h[j][i].SetLineColor(col[i])
        h[j][i].Draw(opt)

def draw2DAngles(c,i,h):
    for k in h:
        for j in range(3):
            c[k].cd(j+1)
            c[k].cd(j+1).SetLogz()
            h[k][j][i].Draw('colz')

def drawClSizes(c,graphs,opt):
    gcounter=-1
    for g in graphs:
        gcounter+=1
        c.cd(gcounter*2+1)
        g['RMS'].Draw(opt)
        c.cd(gcounter*2+2)
        g['Mean'].Draw(opt)
        
def sign(x):
    return(x/abs(x))

def setCanvas(c,d1,d2,leg=None):
    c.Divide(2,1)
    c.cd(1).SetPad(0,0,0.8,1)
    c.cd(1).Divide(d1,d2)
    c.cd(2).SetPad(0.8,0,1,1)
    if leg!=None:
        c.cd(2)
        leg.Draw()

def setGraph(g,title,col=1):
    #colors=getColors()

    g.SetTitle(title)
    g.SetMarkerStyle(20)
    g.SetMarkerSize(0.5)
    g.SetMarkerColor(colors[col])
    return

def fillGraph(hist,graph,fnum=0,mini=[],maxi=[]):
    nbins=hist.GetNbinsX()
    projs=[]
    graph.append({})
    graph[-1]['RMS']=ROOT.TGraph()
    setGraph(graph[-1]['RMS'],'RMS '+hist.GetTitle(),fnum)
    graph[-1]['Mean']=ROOT.TGraph()
    setGraph(graph[-1]['Mean'],'Mean '+hist.GetTitle(),fnum)

    pcounter=0
    i=args.binwidth
    while 1:
        projs.append(hist.ProjectionY(hist.GetName()+titles[fnum]+'_'+str(i),i-args.binwidth,i))
        graph[-1]['RMS'].SetPoint(pcounter,i-args.binwidth+1,projs[-1].GetRMS())
        graph[-1]['Mean'].SetPoint(pcounter,i-args.binwidth+1,projs[-1].GetMean())
        mini['RMS']=min(mini['RMS'],projs[-1].GetRMS())
        maxi['RMS']=max(maxi['RMS'],projs[-1].GetRMS())
        mini['Mean']=min(mini['Mean'],projs[-1].GetMean())
        maxi['Mean']=max(maxi['Mean'],projs[-1].GetMean())
        pcounter+=1
        i+=1
        if (i>=nbins):
            break

    '''
    for i in range(args.binwidth,nbins,args.binwidth):
        projs.append(hist.ProjectionY(hist.GetName()+'_'+str(i),i-args.binwidth,i))
        graph[-1]['RMS'].SetPoint(pcounter,i,projs[-1].GetRMS())
        graph[-1]['Mean'].SetPoint(pcounter,i,projs[-1].GetMean())
        mini['RMS']=min(mini['RMS'],projs[-1].GetRMS())
        maxi['RMS']=max(maxi['RMS'],projs[-1].GetRMS())
        mini['Mean']=min(mini['Mean'],projs[-1].GetMean())
        maxi['Mean']=max(maxi['Mean'],projs[-1].GetMean())

        pcounter+=1
        '''
    return projs

def fillClusterSizes(cl,e,phi=-1):

    #xsize=[cl.pos().X(),cl.pos().X()]
    #ysize=[cl.pos().Y(),cl.pos().Y()]
    xsize=[0,99]
    ysize=[0,99]
    #zsize=[cl.pos().Z(),cl.pos().Z()]
    zsize=[0,511]
    lsize=[0,99]#[cl.pos().Z(),cl.pos().Z()]
    p1size=[0,99]#[cl.pos().Z(),cl.pos().Z()]
    p2size=[0,99]#[cl.pos().Z(),cl.pos().Z()]
    
    a1size=[0,99]
    a2size=[0,99]
    a3size=[0,99]

    if args.recalcClSize:
        digimap=cl.getDigiMap()
        for dig in digimap:
            digi=e.TpcDigi.At(dig.first)
            if type(digi)!=ROOT.TpcDigi:
                print 'no digi',type(digi)
                continue
            dx=float(DigiMapper[digi.padId()][0])
            dy=float(DigiMapper[digi.padId()][1])
            #dz=(digi.t()*args.ftbin+args.ft0)*args.vd+args.fzGem
            dz=digi.t()
            xsize[0]=max(xsize[0],dx)
            ysize[0]=max(ysize[0],dy)
            zsize[0]=max(zsize[0],dz)
            xsize[1]=min(xsize[1],dx)
            ysize[1]=min(ysize[1],dy)
            zsize[1]=min(zsize[1],dz)

            dPos=ROOT.TVector3(dx,dy,dz)
            dist=dPos-cl.pos()
            if phi!=-1:
                alpha=dist.Phi()-phi
    #            print 'alpha=',alpha*ROOT.TMath.RadToDeg(),'dphi=',dist.Phi()*ROOT.TMath.RadToDeg(),'phi=',phi*ROOT.TMath.RadToDeg()
    #            dPos.Print()
    #            cl.pos().Print()
                longi=math.cos(alpha)*dist.Perp()
                perp=math.sin(alpha)*dist.Perp()
                lsize[0]=max(lsize[0],longi)
                p1size[0]=max(p1size[0],perp)
                lsize[1]=min(lsize[1],longi)
                p1size[1]=min(p1size[1],perp)


    #        cl.axis(1).Print()
    #        cl.cov().Print()
    #        if cl.cov().Determinant()==0:
    #            cl.cov().Print()
    #            cl.pos().Print()
    #            dPos.Print()
            a1=cl.axis(0)*dist
            a2=cl.axis(1)*dist
            a3=cl.axis(2)*dist

            a1size[0]=max(a1size[0],a1)
            a1size[1]=min(a1size[1],a1)
            a2size[0]=max(a2size[0],a2)
            a2size[1]=min(a2size[1],a2)
            a3size[0]=max(a3size[0],a3)
            a3size[1]=min(a3size[1],a3)

        xlen=0
        ylen=0
        zlen=0
        llen=0
        p1len=0
        p2len=0
        a1len=0
        a2len=0
        a3len=0

        xlen=xsize[0]-xsize[1]
        ylen=ysize[0]-ysize[1]
        zlen=zsize[0]-zsize[1]
    #    if zsize[0]==zsize[1]:
    #        zlen=-0.5
    #    print zsize[0],zsize[1],zlen

        if lsize[0]*lsize[1]<0:
            llen=lsize[0]-lsize[1]
        else:
            llen=max(lsize)
        if p1size[0]*p1size[1]<0:
            p1len=p1size[0]-p1size[1]
        else:
            p1len=max(p1size)
        p2len=zlen

        a1len=a1size[0]-a1size[1]
        a2len=a2size[0]-a2size[1]
        a3len=a3size[0]-a3size[1]
    else:
        xlen=0
        ylen=0
        zlen=0
        llen=0
        p1len=0
        p2len=0
        a1len=0
        a2len=0
        a3len=0

        xyzlen=tr.GetFullClusterSize()[clcounter]
        xlen=xyzlen.X()
        ylen=xyzlen.Y()
        zlen=xyzlen.Z()
        lpplen=tr.GetFullClusterSizeTrack()[clcounter]
        llen=lpplen.X()
        p1len=lpplen.Y()
        p2len=xyzlen.Z()
        
        axlen=tr.GetFullClusterSizeAxis()[clcounter]
        a1len=axlen.X()
        a2len=axlen.Y()
        a3len=axlen.Z()
        

    hClustersizes[0][fcounter].Fill(cl.pos().Z(),xlen)
    hClustersizes[1][fcounter].Fill(cl.pos().Z(),ylen)
    hClustersizes[2][fcounter].Fill(cl.pos().Z(),zlen)

    if lsize[1]==99 and args.recalcClSize:
        pass
    else:
        hClustersizes[3][fcounter].Fill(cl.pos().Z(),llen)
        hClustersizes[4][fcounter].Fill(cl.pos().Z(),p1len)
        hClustersizes[5][fcounter].Fill(cl.pos().Z(),p2len)

    hClustersizes[6][fcounter].Fill(cl.pos().Z(),a1len)
    hClustersizes[7][fcounter].Fill(cl.pos().Z(),a2len)
    hClustersizes[8][fcounter].Fill(cl.pos().Z(),a3len)
    return

def fill2dAngles(key,fcounter,Z,vals):
    for i in range(len(vals)):
        all2DAngles[key][i][fcounter].Fill(Z,vals[i])

parser=argparse.ArgumentParser(description="plots track related infos like cluster per cm, cluster per track ...")
parser.add_argument('files',help='the histo file',type=str,nargs="+")
parser.add_argument('--digifiles',help='the digifiles',type=str,nargs='+')
parser.add_argument('--titles',help='the titles to use',type=str,default=['None'],nargs='+')
parser.add_argument('--sqrt',help='plot a sqrt dependence',action='store_true')
parser.add_argument('--nevts',help='number of events to analyze',type=int,default=-1)
parser.add_argument('--list',help='do for all files in list',action='store_true')
parser.add_argument('--raw',help='use raw cluster information, not TpcCosmics',action='store_true')
parser.add_argument('--new',help='create new analysis',action='store_true')
parser.add_argument('--histfile',help='name of the histfile',type=str,default='clusdist_hists.root')
parser.add_argument('--basepath',help='the path to save all hists',type=str,default='./')
parser.add_argument('--pfile',help='the file to save all canvas',type=str,default='None')
parser.add_argument('--padf',help='the padplane file to use',type=str,default='tpc/parfiles/padPlane_FOPI.dat')
parser.add_argument('--ft0',help='ft for digi z calc',type=float,default=0)
parser.add_argument('--ftbin',help='ftbin for digi z calc',type=float,default=64.3087)
parser.add_argument('--vd',help='drift velocity for digi z calc',type=float,default=0.00237708)
parser.add_argument('--fzGem',help='z gem for digi z calc',type=float,default=0)
parser.add_argument('--binwidth',help='width of bins',type=int,default=5)
parser.add_argument('--anima',help='create phi and theta animation',action='store_true')
parser.add_argument('--mc',help='user MC angles',action='store_true')
parser.add_argument('--phicut',help='window for cut. only inside is taken',type=float,nargs=2)
parser.add_argument('--clsizes',help='plot all clustersizes',action='store_true')
parser.add_argument('--recalcClSize',help='recalculate cluster sizes',action='store_true')
parser.add_argument('--correct',help='do corrections',action='store_true')
parser.add_argument('--batch',help='no plotting',action='store_true')
parser.add_argument('--hlp',help='display help',action='store_true')
args=parser.parse_args()

if args.hlp:
    parser.print_help()
    exit(0)

if args.batch:
    sys.argv.append( '-b-' )
import ROOT
if args.batch:
    ROOT.gROOT.SetBatch(True)

ROOT.gROOT.ProcessLine(".x rootlogon.C") 
ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")')
ROOT.gROOT.LoadMacro("stlPYROOT.h+")
set_palette()

if args.mc and args.histfile=='clusdist_hists.root':
    args.histfile.replace('.root','_mc.root')

DigiMapper = []
file=open(args.padf,'r')
lines=file.readlines()
file.close()
lines.reverse()
for line in lines:
    koord = line.split(" ") 
    DigiMapper.append([koord[3], koord[4]])

hClustersizes=[]
for i in range(9):
    hClustersizes.append([])

hdphi=[]
hdphi2=[]
hdtheta=[]
hdangle=[]
hdangleAZ=[]
hdsorted=[]
hphiA=[]
hthetaA=[]
hphiT=[]
hthetaT=[]
#hthetaAZ=[]
#hphiAZ=[]

all2DAngles={'phiVsZ':[],'thetaVsZ':[],'dphiVsZ':[],'dthetaVsZ':[],'dangleVsZ':[],'dangleAZVsZ':[]}
for a in all2DAngles.values():
    for i in range(3):
        a.append([])

for i in range(3):
    hdphi.append([])
    hdtheta.append([])
    hdsorted.append([])
    hphiA.append([])
    hthetaA.append([])
    hphiT.append([])
    hthetaT.append([])
    hdphi2.append([])
    hdangle.append([])
    hdangleAZ.append([])


    #hthetaAZ.append([])
    #hphiAZ.append([])

hphiUV=[]
hthetaUV=[]
hdphiUV=[]
hdthetaUV=[]
hdangleUV=[]
for i in range(2):
    hphiUV.append([])
    hthetaUV.append([])
    hdphiUV.append([])
    hdthetaUV.append([])
    hdangleUV.append([])


if args.list:
    listfile=open(args.files[0],'r')
    files=[]
    digifiles=[]
    titles=[]
    for line in listfile:
        if line[0]=='#' or line[0]=='\n':
            continue
        words=line.split(';')
        files.append(str(words[0]))
        titles.append(str(words[1]).replace('\n',''))
        digifiles.append(str(words[2]).replace('\n',''))
    listfile.close()
else:
    files=args.files
    digifiles=args.digifiles
    titles=[]    
    if args.titles[0]=='None':
        for i in range(len(files)):
            titles=str(i)
    elif len(titles)<len(files):
        titles=args.titles

Xaxis=ROOT.TVector3(1,0,0)
Yaxis=ROOT.TVector3(0,1,0)
Zaxis=ROOT.TVector3(0,0,1)

if args.new:
    fcounter=-1
    for f in files:
        print f
        fcounter+=1
        #RFile=ROOT.TFile()
        Rfile = ROOT.TFile.Open(f,"read")
        tree = Rfile.Get("cbmsim")
        if tree == None :
            print "No cbmsim in file", file
            if type(Rfile)==ROOT.TFile:
                Rfile.Close()
            exit()

        tree.AddFriend('cbmsim',digifiles[fcounter])
        tree.SetBranchStatus("*", 0)
        tree.SetBranchStatus("CutCosmics.*",1)
        tree.SetBranchStatus("TpcCluster.*",1)
        tree.SetBranchStatus("TpcDigi.*",1)
        tree.SetBranchStatus("TrackFitStat_0.*",1)
        tree.SetBranchStatus("TpcSPHit.*",1)
        tree.SetBranchStatus("TrackPostFit.*",1)

        hClustersizes[0].append(ROOT.TH2D('hXsize_'+titles[fcounter],'X'+titles[fcounter],73,0,73,25,0,2.5))
        hClustersizes[1].append(ROOT.TH2D('hYsize_'+titles[fcounter],'Y'+titles[fcounter],73,0,73,25,0,2.5))
        hClustersizes[2].append(ROOT.TH2D('hZsize_'+titles[fcounter],'Z'+titles[fcounter],73,0,73,31,0,2.))

        hClustersizes[3].append(ROOT.TH2D('hLsize_'+titles[fcounter],'L'+titles[fcounter],73,0,73,25,0,2.5))
        hClustersizes[4].append(ROOT.TH2D('hP1size_'+titles[fcounter],'P1'+titles[fcounter],73,0,73,25,0,2.5))
        hClustersizes[5].append(ROOT.TH2D('hP2size_'+titles[fcounter],'P2'+titles[fcounter],73,0,73,25,0,2.5))

        hClustersizes[6].append(ROOT.TH2D('hA1size_'+titles[fcounter],'A1'+titles[fcounter],73,0,73,25,0,2.5))
        hClustersizes[7].append(ROOT.TH2D('hA2size_'+titles[fcounter],'A2'+titles[fcounter],73,0,73,25,0,2.5))
        hClustersizes[8].append(ROOT.TH2D('hA3size_'+titles[fcounter],'A3'+titles[fcounter],73,0,73,25,0,2.5))

        hdphi[0].append(ROOT.TH1D('hdphi0_'+titles[fcounter],'Dphi0'+titles[fcounter],720,-360,360))
        hdphi[1].append(ROOT.TH1D('hdphi1_'+titles[fcounter],'Dphi1'+titles[fcounter],720,-360,360))
        hdphi[2].append(ROOT.TH1D('hdphi2_'+titles[fcounter],'Dphi2'+titles[fcounter],720,-360,360))

        for i in range(3):
            hphiA[i].append(ROOT.TH1D('hphiA'+str(i)+'_'+titles[fcounter],'Phi A '+str(i)+titles[fcounter],362,-181,181))
            hphiT[i].append(ROOT.TH1D('hphiT'+str(i)+'_'+titles[fcounter],'Phi T '+str(i)+titles[fcounter],362,-181,181))
            hthetaA[i].append(ROOT.TH1D('hthetaA'+str(i)+'_'+titles[fcounter],'Theta A '+str(i)+titles[fcounter],362,-181,181))
            hthetaT[i].append(ROOT.TH1D('hthetaT'+str(i)+'_'+titles[fcounter],'Theta T '+str(i)+titles[fcounter],362,-181,181))                
            hdangle[i].append(ROOT.TH1D('hdangle'+str(i)+'_'+titles[fcounter],'DAngle'+str(i)+titles[fcounter],362,-1,181))
            hdangleAZ[i].append(ROOT.TH1D('hdangleAZ'+str(i)+'_'+titles[fcounter],'DAngleAZ'+str(i)+titles[fcounter],362,-1,181))
#            hphiAZ[i].append(ROOT.TH1D('hphiAZ'+str(i)+'_'+titles[fcounter],'Phi AZ '+str(i)+titles[fcounter],722,-361,361))
#            hthetaAZ[i].append(ROOT.TH1D('hThetaAZ'+str(i)+'_'+titles[fcounter],'Theta AZ '+str(i)+titles[fcounter],722,-361,361))
            
        for i in range(2):
            hphiUV[i].append(ROOT.TH1D('hphiUV'+str(i)+'_'+titles[fcounter],'Phi UV '+str(i)+titles[fcounter],720,-360,360))
            hdphiUV[i].append(ROOT.TH1D('hdphiUV'+str(i)+'_'+titles[fcounter],'Dphi UV'+str(i)+titles[fcounter],720,-360,360))
            hthetaUV[i].append(ROOT.TH1D('hthetaUV'+str(i)+'_'+titles[fcounter],'Thteta UV'+str(i)+titles[fcounter],720,-360,360))
            hdthetaUV[i].append(ROOT.TH1D('hdthetaUV'+str(i)+'_'+titles[fcounter],'Dtheta UV '+str(i)+titles[fcounter],720,-360,360))
            hdangleUV[i].append(ROOT.TH1D('hdangleUV'+str(i)+'_'+titles[fcounter],'Dangle UV'+str(i)+titles[fcounter],360,0,360))
            


        hdphi2[0].append(ROOT.TH1D('hdphi2_0_'+titles[fcounter],'Dphi v2 0 '+titles[fcounter],360,-90,90))
        hdphi2[1].append(ROOT.TH1D('hdphi2_1_'+titles[fcounter],'Dphi v2 1 '+titles[fcounter],360,0,180))
        hdphi2[2].append(ROOT.TH1D('hdphi2_2_'+titles[fcounter],'Dphi v2 2 '+titles[fcounter],360,-90,90))

        hdtheta[0].append(ROOT.TH1D('hdtheta0_'+titles[fcounter],'Dtheta1'+titles[fcounter],360,-180,180))
        hdtheta[1].append(ROOT.TH1D('hdtheta1_'+titles[fcounter],'Dtheta2'+titles[fcounter],360,-180,180))
        hdtheta[2].append(ROOT.TH1D('hdtheta2_'+titles[fcounter],'Dtheta3'+titles[fcounter],360,-180,180))

#        hdsorted[0].append(ROOT.TH1D('hdsorted1_'titles[fcounter],'

        for k,v in all2DAngles.iteritems():
            for i in range(3):
                v[i].append(ROOT.TH2D('h'+k+str(i)+titles[fcounter],k+str(i)+titles[fcounter],73,0,73,360,-180,180))


        evt=-1
        done=0
        olddone=0
        for e in tree:
            evt+=1
            if args.nevts!=-1 and evt>args.nevts:
                break
            if evt%100==0:
                print evt

            if args.raw:
                for cl in e.TpcCluster:
                    fillClusterSizes(cl,e)

            else:
                if e.TrackFitStat_0.GetEntries()>1:
                    continue
                trackCounter=-1
                for tr in e.TrackFitStat_0:
                    trackCounter+=1
#                    if tr.GetStatFlag()!=0:
#                        continue
                    #track=e.TrackPostFit.At(trackCounter)
                    #rep = track.getTrackRep(0)
                    posX=tr.GetHitPositionsX()
                    posY=tr.GetHitPositionsY()
                    planeU=tr.GetPlaneU()
                    planeV=tr.GetPlaneV()
                    if args.mc:
                        PosMoms=tr.GetMcPosMom()
                    else:
                        PosMoms=tr.GetProjectionPointsMom()
                    
                    clcounter=-1
                    if len(tr.GetHitIDs())!=len(PosMoms):
                        print len(tr.GetHitIDs()),len(PosMoms)
                        continue
                    '''
                    if not args.recalcClSize:
                        for clsize in range(len(tr.GetFullClusterSize())):
                            fillClusterSizes(clsize,e,PosMoms[clcounter].Phi())
                            '''

                    for clid in tr.GetHitIDs():
                        clcounter+=1
                        SpHit=e.TpcSPHit.At(clid)
                        cl=e.TpcCluster.At(SpHit.getClusterID())
                        if args.recalcClSize:
                            fillClusterSizes(cl,e,PosMoms[clcounter].Phi())
                        else:
                            fillClusterSizes(cl,e,PosMoms[clcounter].Phi())
                       
                        phidiff1=(cl.axis(0).Phi()-PosMoms[clcounter].Phi())*ROOT.TMath.RadToDeg()
                        phidiff2=(cl.axis(1).Phi()-PosMoms[clcounter].Phi())*ROOT.TMath.RadToDeg()
                        phidiff3=(cl.axis(2).Phi()-PosMoms[clcounter].Phi())*ROOT.TMath.RadToDeg()
                        hdphi[0][-1].Fill(phidiff1)
                        hdphi[1][-1].Fill(phidiff2)
                        hdphi[2][-1].Fill(phidiff3)

                        if cl.axis(0)*PosMoms[clcounter]<0 and args.correct:
                            thetadiff1=((ROOT.TMath.Pi()-cl.axis(0).Theta())-PosMoms[clcounter].Theta())*ROOT.TMath.RadToDeg()
                            thetadiff1*=-1
                        else:
                            thetadiff1=(cl.axis(0).Theta()-PosMoms[clcounter].Theta())*ROOT.TMath.RadToDeg()


                        thetadiff2=(cl.axis(1).Theta()-PosMoms[clcounter].Theta())*ROOT.TMath.RadToDeg()
                        if cl.axis(2)*Yaxis>0:
                            thetadiff3=((ROOT.TMath.Pi()-cl.axis(2).Theta())-PosMoms[clcounter].Theta())*ROOT.TMath.RadToDeg()
                        else:
                            thetadiff3=(cl.axis(2).Theta()-PosMoms[clcounter].Theta())*ROOT.TMath.RadToDeg()

                        hdtheta[0][-1].Fill(thetadiff1)
                        hdtheta[1][-1].Fill(thetadiff2)
                        hdtheta[2][-1].Fill(thetadiff3)
                        fill2dAngles('dthetaVsZ',fcounter,cl.pos().Z(),[thetadiff1,thetadiff2,thetadiff3])
                        

                        if abs(cl.axis(0).DeltaPhi(PosMoms[clcounter]))<ROOT.TMath.Pi()/2:
                            phidiff1_v2=cl.axis(0).DeltaPhi(PosMoms[clcounter])*ROOT.TMath.RadToDeg()
                        else:
                            phidiff1_v2=(sign(cl.axis(0).DeltaPhi(PosMoms[clcounter]))*180
                                         -cl.axis(0).DeltaPhi(PosMoms[clcounter])*ROOT.TMath.RadToDeg())

                        phidiff2_v2=abs(cl.axis(1).DeltaPhi(PosMoms[clcounter])*ROOT.TMath.RadToDeg())
                        
                        if abs(cl.axis(2).DeltaPhi(PosMoms[clcounter]))<ROOT.TMath.Pi()/2:
                            phidiff3_v2=cl.axis(2).DeltaPhi(PosMoms[clcounter])*ROOT.TMath.RadToDeg()
                        else:
                            phidiff3_v2=(sign(cl.axis(2).DeltaPhi(PosMoms[clcounter]))*180
                                         -cl.axis(2).DeltaPhi(PosMoms[clcounter])*ROOT.TMath.RadToDeg())


                        hdphi2[0][-1].Fill( phidiff1_v2 )
                        hdphi2[1][-1].Fill( phidiff2_v2 )
                        hdphi2[2][-1].Fill( phidiff3_v2 )
                        fill2dAngles('dphiVsZ',fcounter,cl.pos().Z(),[phidiff1_v2,phidiff2_v2,phidiff3_v2])
                        
                        dangle=[]
                        dangleAZ=[]
                        for i in range(3):
                            dangle.append(cl.axis(i).Angle(PosMoms[clcounter])*ROOT.TMath.RadToDeg())
                            dangleAZ.append(cl.axis(i).Angle(Zaxis)*ROOT.TMath.RadToDeg())

                        hdangle[0][-1].Fill(dangle[0])
                        hdangle[1][-1].Fill(dangle[1])
                        hdangle[2][-1].Fill(dangle[2])
                        fill2dAngles('dangleVsZ',fcounter,cl.pos().Z(),dangle)
                        
                        hdangleAZ[0][-1].Fill(dangleAZ[0])
                        hdangleAZ[1][-1].Fill(dangleAZ[1])
                        hdangleAZ[2][-1].Fill(dangleAZ[2])
                        fill2dAngles('dangleAZVsZ',fcounter,cl.pos().Z(),dangleAZ)


                        hphiUV[0][-1].Fill(planeU[clcounter].Phi()*ROOT.TMath.RadToDeg())
                        hphiUV[1][-1].Fill(planeV[clcounter].Phi()*ROOT.TMath.RadToDeg())

                        hthetaUV[0][-1].Fill(planeU[clcounter].Theta()*ROOT.TMath.RadToDeg())
                        hthetaUV[1][-1].Fill(planeV[clcounter].Theta()*ROOT.TMath.RadToDeg())

                        hdphiUV[0][-1].Fill(planeU[clcounter].DeltaPhi(PosMoms[clcounter])*ROOT.TMath.RadToDeg())
                        hdphiUV[1][-1].Fill(planeV[clcounter].DeltaPhi(PosMoms[clcounter])*ROOT.TMath.RadToDeg())

                        thetadiffU=planeU[clcounter].Theta()-PosMoms[clcounter].Theta()
                        thetadiffV=planeV[clcounter].Theta()-PosMoms[clcounter].Theta()
                        hdthetaUV[0][-1].Fill(thetadiffU*ROOT.TMath.RadToDeg())
                        hdthetaUV[1][-1].Fill(thetadiffV*ROOT.TMath.RadToDeg())

                        hdangleUV[0][-1].Fill(planeU[clcounter].Angle(PosMoms[clcounter])*ROOT.TMath.RadToDeg())
                        hdangleUV[1][-1].Fill(planeV[clcounter].Angle(PosMoms[clcounter])*ROOT.TMath.RadToDeg())

                        phiA=[]
                        thetaA=[]
                        for i in range(3):
                            phiA.append(cl.axis(i).Phi()*ROOT.TMath.RadToDeg())
                            thetaA.append(cl.axis(i).Theta()*ROOT.TMath.RadToDeg())

                        for i in range(3):
                            hphiA[i][-1].Fill(phiA[i])
                            hphiT[i][-1].Fill(PosMoms[clcounter].Phi()*ROOT.TMath.RadToDeg())
                            hthetaA[i][-1].Fill(thetaA[i])
                            hthetaT[i][-1].Fill(PosMoms[clcounter].Theta()*ROOT.TMath.RadToDeg())

                        fill2dAngles('phiVsZ',fcounter,cl.pos().Z(),phiA)
                        fill2dAngles('thetaVsZ',fcounter,cl.pos().Z(),thetaA)

    histfile=ROOT.TFile(args.basepath+args.histfile,'recreate')
    writeHist(hdphi)
    writeHist(hdtheta)
    writeHist(hphiA)
    writeHist(hphiT)
    writeHist(hthetaA)
    writeHist(hthetaT)
    writeHist(hdphi2)
    writeHist(hdangle)            
    writeHist(hdangleAZ)
    writeHist(hphiUV)
    writeHist(hthetaUV)
    writeHist(hdphiUV)
    writeHist(hdthetaUV)
    writeHist(hdangleUV)
    writeHist(hClustersizes)
    for v in all2DAngles.values():
        writeHist(v)
    histfile.Close()
else:
    histfile=ROOT.TFile(args.basepath+args.histfile,'update')
    histfile.GetListOfKeys().Print()
    print 'using histfile:',args.basepath+args.histfile
    for tit in titles:
        print 'getting histos for:',tit
        hClustersizes[0].append(histfile.Get('hXsize_'+tit))
        hClustersizes[1].append(histfile.Get('hYsize_'+tit))
        hClustersizes[2].append(histfile.Get('hZsize_'+tit))
        hClustersizes[3].append(histfile.Get('hLsize_'+tit))
        hClustersizes[4].append(histfile.Get('hP1size_'+tit))
        hClustersizes[5].append(histfile.Get('hP2size_'+tit))
        hClustersizes[6].append(histfile.Get('hA1size_'+tit))
        hClustersizes[7].append(histfile.Get('hA2size_'+tit))
        hClustersizes[8].append(histfile.Get('hA3size_'+tit))
        for i in range(3):
            hdphi[i].append(histfile.Get('hdphi'+str(i)+'_'+tit))
            hdtheta[i].append(histfile.Get('hdtheta'+str(i)+'_'+tit))
            hphiA[i].append(histfile.Get('hphiA'+str(i)+'_'+tit))
            hphiT[i].append(histfile.Get('hphiT'+str(i)+'_'+tit))
            hthetaA[i].append(histfile.Get('hthetaA'+str(i)+'_'+tit))
            hthetaT[i].append(histfile.Get('hthetaT'+str(i)+'_'+tit))
            hdphi2[i].append(histfile.Get('hdphi2_'+str(i)+'_'+tit))
            hdangle[i].append(histfile.Get('hdangle'+str(i)+'_'+tit))
            hdangleAZ[i].append(histfile.Get('hdangleAZ'+str(i)+'_'+tit))
        for i in range(2):
            hphiUV[i].append(histfile.Get('hphiUV'+str(i)+'_'+tit))
            hthetaUV[i].append(histfile.Get('hthetaUV'+str(i)+'_'+tit))
            hdphiUV[i].append(histfile.Get('hdphiUV'+str(i)+'_'+tit))
            hdthetaUV[i].append(histfile.Get('hdthetaUV'+str(i)+'_'+tit))
            hdangleUV[i].append(histfile.Get('hdangleUV'+str(i)+'_'+tit))
        for k,v in all2DAngles.iteritems():
            for i in range(3):
                v[i].append(histfile.Get('h'+k+str(i)+tit))



chClustersizes=[]
gClustersizes=[]
gClustersizes2=[]
minClustersizes=[]
maxClustersizes=[]

for i in range(9):
    gClustersizes.append([])
    minClustersizes.append({})
    minClustersizes[-1]['RMS']=1000
    minClustersizes[-1]['Mean']=1000
    maxClustersizes.append({})
    maxClustersizes[-1]['RMS']=0
    maxClustersizes[-1]['Mean']=0

leg=ROOT.TLegend(0,0,1,1)
leg.SetFillColor(0)
leg.SetTextSize(0.1)
legp=[]
counter=-1
if len(files)<10:
    colors=getColors2()
else:
    colors=getColors()

for tit in titles: 
    print 'plotting histos for:',tit
    counter+=1
    if args.clsizes:
        chClustersizes.append(ROOT.TCanvas('chClustersizes'+tit,'Clustersizes H '+tit,500,1000))
        chClustersizes[-1].Divide(1,3)
        chClustersizes[-1].cd(1)
        hClustersizes[0][counter].Draw('colz')
        chClustersizes[-1].cd(2)
        hClustersizes[1][counter].Draw('colz')
        chClustersizes[-1].cd(3)
        hClustersizes[2][counter].Draw('colz')
        chClustersizes[-1].Update()

        chClustersizes.append(ROOT.TCanvas('chClustersizes2'+tit,'Clustersizes2 H '+tit,500,1000))
        chClustersizes[-1].Divide(1,3)
        chClustersizes[-1].cd(1)
        hClustersizes[3][counter].Draw('colz')
        chClustersizes[-1].cd(2)
        hClustersizes[4][counter].Draw('colz')
        chClustersizes[-1].cd(3)
        hClustersizes[5][counter].Draw('colz')
        chClustersizes[-1].Update()

        chClustersizes.append(ROOT.TCanvas('chClustersizesA'+tit,'Clustersizes A '+tit,500,1000))
        chClustersizes[-1].Divide(1,3)
        chClustersizes[-1].cd(1)
        hClustersizes[6][counter].Draw('colz')
        chClustersizes[-1].cd(2)
        hClustersizes[7][counter].Draw('colz')
        chClustersizes[-1].cd(3)
        hClustersizes[8][counter].Draw('colz')
        chClustersizes[-1].Update()

    fillGraph(hClustersizes[0][counter],gClustersizes[0],counter,minClustersizes[0],maxClustersizes[0])
    fillGraph(hClustersizes[1][counter],gClustersizes[1],counter,minClustersizes[1],maxClustersizes[1])
    fillGraph(hClustersizes[2][counter],gClustersizes[2],counter,minClustersizes[2],maxClustersizes[2])

    fillGraph(hClustersizes[3][counter],gClustersizes[3],counter,minClustersizes[3],maxClustersizes[3])
    fillGraph(hClustersizes[4][counter],gClustersizes[4],counter,minClustersizes[4],maxClustersizes[4])
    fillGraph(hClustersizes[5][counter],gClustersizes[5],counter,minClustersizes[5],maxClustersizes[5])

    fillGraph(hClustersizes[6][counter],gClustersizes[6],counter,minClustersizes[6],maxClustersizes[6])
    fillGraph(hClustersizes[7][counter],gClustersizes[7],counter,minClustersizes[7],maxClustersizes[7])
    fillGraph(hClustersizes[8][counter],gClustersizes[8],counter,minClustersizes[8],maxClustersizes[8])


    opt='P'
    if counter==0:
        opt='AP'
        cgClustersizes=ROOT.TCanvas('cgClustersizes'+tit,'Clustersizes G '+tit,1000,1000)
        setCanvas(cgClustersizes,2,3,leg)
        #cgClustersizes.Divide(2,3)
        cgClustersizes2=ROOT.TCanvas('cgClustersizes2'+tit,'Clustersizes2 G '+tit,1000,1000)
        setCanvas(cgClustersizes2,2,3,leg)
        cgClustersizesA=ROOT.TCanvas('cgClustersizesA'+tit,'Clustersizes A '+tit,1000,1000)
        setCanvas(cgClustersizesA,2,3,leg)
        


    drawClSizes(cgClustersizes.cd(1),[ gClustersizes[0][-1], gClustersizes[1][-1], gClustersizes[2][-1] ],opt)
    #draw the graphs for along track cluster
    drawClSizes( cgClustersizes2.cd(1), [ gClustersizes[3][-1], gClustersizes[4][-1], gClustersizes[5][-1] ], opt)
    #draw the graphs for along cl.cov axis
    drawClSizes(cgClustersizesA.cd(1),[ gClustersizes[6][-1],gClustersizes[7][-1],gClustersizes[8][-1] ],opt)

    legp.append(ROOT.TPolyMarker())
    legp[-1].SetMarkerStyle(gClustersizes[2][-1]['Mean'].GetMarkerStyle())
    legp[-1].SetMarkerSize(3)
    legp[-1].SetMarkerColor(gClustersizes[2][-1]['Mean'].GetMarkerColor())
    leg.AddEntry(legp[-1],tit,"P")

    opt2='same'
    if counter==0:
        opt2=''
        cphi=ROOT.TCanvas('cphi'+tit,'Phi_'+tit,700,1000)
        setCanvas(cphi,1,3,leg)
        ctheta=ROOT.TCanvas('ctheta'+tit,'Theta_'+tit,700,1000)
        setCanvas(ctheta,1,3,leg)
        cphiA=ROOT.TCanvas('cphiA'+tit,'Phi_A_'+tit,700,1000)
        setCanvas(cphiA,1,3,leg)
        cphiT=ROOT.TCanvas('cphiT'+tit,'Phi_T_'+tit,700,1000)
        setCanvas(cphiT,1,3,leg)
        cthetaA=ROOT.TCanvas('cthetaA'+tit,'Theta_A_'+tit,700,1000)
        setCanvas(cthetaA,1,3,leg)
        cthetaT=ROOT.TCanvas('cthetaT'+tit,'Theta_T_'+tit,700,1000)
        setCanvas(cthetaT,1,3,leg)
        cphi2=ROOT.TCanvas('cphi2'+tit,'Phi2_'+tit,700,1000)
        setCanvas(cphi2,1,3,leg)
        cangle=ROOT.TCanvas('cangle'+tit,'Angle_'+tit,700,1000)
        setCanvas(cangle,1,3,leg)
        cphiUV=ROOT.TCanvas('cphiUV'+tit,'Phi_UV_'+tit,700,1000)
        setCanvas(cphiUV,1,2,leg)
        cthetaUV=ROOT.TCanvas('cthetaUV'+tit,'Theta_UV_'+tit,700,1000)
        setCanvas(cthetaUV,1,2,leg)
        cdphiUV=ROOT.TCanvas('cdphiUV'+tit,'D_phi_UV_'+tit,700,1000)
        setCanvas(cdphiUV,1,2,leg)
        cdthetaUV=ROOT.TCanvas('cdthetaUV'+tit,'D_theta_UV_'+tit,700,1000)
        setCanvas(cdthetaUV,1,2,leg)
        cdangleUV=ROOT.TCanvas('cdangleUV'+tit,'D_angle_UV_'+tit,700,1000)
        setCanvas(cdangleUV,1,2,leg)
        cdangleAZ=ROOT.TCanvas('cdangleAZ'+tit,'D_angle_AZ_'+tit,700,1000)
        setCanvas(cdangleAZ,1,3,leg)
        
        c2dAngles={}

        if args.anima:
            cphi2=ROOT.TCanvas('cphi2'+tit,'Phi_'+tit,700,1000)
            setCanvas(cphi2,1,3,leg)
            ctheta2=ROOT.TCanvas('ctheta2'+tit,'Theta_'+tit,700,1000)
            setCanvas(ctheta2,1,3,leg)

    c2dAngles[tit]={}
    for k in all2DAngles:
        print 'k=',k
        c2dAngles[tit][k]=ROOT.TCanvas('c'+k+tit,k+tit,500,1000)
        c2dAngles[tit][k].Divide(1,3)

    draw2DAngles(c2dAngles[tit],counter,all2DAngles)

    drawAngles(cphiUV.cd(1),hphiUV,counter,colors,opt2)
    drawAngles(cthetaUV.cd(1),hthetaUV,counter,colors,opt2)
    drawAngles(cdphiUV.cd(1),hdphiUV,counter,colors,opt2)
    drawAngles(cdthetaUV.cd(1),hdthetaUV,counter,colors,opt2)
    drawAngles(cdangleUV.cd(1),hdangleUV,counter,colors,opt2)
    drawAngles(cphi.cd(1),hdphi,counter,colors,opt2)
    drawAngles(ctheta.cd(1),hdtheta,counter,colors,opt2)
    drawAngles(cphiA.cd(1),hphiA,counter,colors,opt2)
    drawAngles(cthetaA.cd(1),hthetaA,counter,colors,opt2)
    drawAngles(cphiT.cd(1),hphiT,counter,colors,opt2)
    drawAngles(cthetaT.cd(1),hthetaT,counter,colors,opt2)
    drawAngles(cphi2.cd(1),hdphi2,counter,colors,opt2)
    drawAngles(cangle.cd(1),hdangle,counter,colors,opt2)
    drawAngles(cdangleAZ.cd(1),hdangleAZ,counter,colors,opt2)

    if args.anima:
        cphi2.cd(1)
        hdphi[0][counter].SetLineColor(colors[counter])
        hdphi[0][counter].Draw()
        cphi2.cd(2)
        hdphi[1][counter].SetLineColor(colors[counter])
        hdphi[1][counter].Draw()
        cphi2.cd(3)
        hdphi[2][counter].SetLineColor(colors[counter])
        hdphi[2][counter].Draw()
        cphi2.Update()
        cphi2.Print(args.basepath+'/phi.gif+100')

        ctheta2.cd(1)
        hdtheta[0][counter].SetLineColor(colors[counter])
        hdtheta[0][counter].Draw()
        ctheta2.cd(2)
        hdtheta[1][counter].SetLineColor(colors[counter])
        hdtheta[1][counter].Draw()
        ctheta2.cd(3)
        hdtheta[2][counter].SetLineColor(colors[counter])
        hdtheta[2][counter].Draw()
        ctheta2.Update()
        ctheta2.Print(args.basepath+'/theta.gif+100')

gClustersizes[0][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[0]['RMS']*0.9,maxClustersizes[0]['RMS']*1.1)
gClustersizes[1][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[1]['RMS']*0.9,maxClustersizes[1]['RMS']*1.1)
gClustersizes[2][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[2]['RMS']*0.9,maxClustersizes[2]['RMS']*1.1)
gClustersizes[0][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[0]['Mean']*0.9,maxClustersizes[0]['Mean']*1.1)
gClustersizes[1][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[1]['Mean']*0.9,maxClustersizes[1]['Mean']*1.1)
gClustersizes[2][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[2]['Mean']*0.9,maxClustersizes[2]['Mean']*1.1)
cgClustersizes.Update()

gClustersizes[3][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[3]['RMS']*0.9,maxClustersizes[3]['RMS']*1.1)
gClustersizes[4][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[4]['RMS']*0.9,maxClustersizes[4]['RMS']*1.1)
gClustersizes[5][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[5]['RMS']*0.9,maxClustersizes[5]['RMS']*1.1)
gClustersizes[3][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[3]['Mean']*0.9,maxClustersizes[3]['Mean']*1.1)
gClustersizes[4][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[4]['Mean']*0.9,maxClustersizes[4]['Mean']*1.1)
gClustersizes[5][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[5]['Mean']*0.9,maxClustersizes[5]['Mean']*1.1)
cgClustersizes2.Update()

gClustersizes[6][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[6]['RMS']*0.9,maxClustersizes[6]['RMS']*1.1)
gClustersizes[7][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[7]['RMS']*0.9,maxClustersizes[7]['RMS']*1.1)
gClustersizes[8][0]['RMS'].GetYaxis().SetRangeUser(minClustersizes[8]['RMS']*0.9,maxClustersizes[8]['RMS']*1.1)
gClustersizes[6][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[6]['Mean']*0.9,maxClustersizes[6]['Mean']*1.1)
gClustersizes[7][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[7]['Mean']*0.9,maxClustersizes[7]['Mean']*1.1)
gClustersizes[8][0]['Mean'].GetYaxis().SetRangeUser(minClustersizes[8]['Mean']*0.9,maxClustersizes[8]['Mean']*1.1)
cgClustersizesA.Update()

allSizeCanv=[cgClustersizes,cgClustersizes2,cgClustersizesA]
if args.clsizes:
    for c in range(len(chClustersizes)):
        allSizeCanv.append(chClustersizes[c])

cphi.Update()
ctheta.Update()
cphiA.Update()
cthetaA.Update()
cphiT.Update()
cthetaT.Update()
cangle.Update()
cphi2.Update()
for tit in c2dAngles:
    for k in c2dAngles[tit]:
        c2dAngles[tit][k].Update()

cLegend=ROOT.TCanvas('cleg','Legend',500,1000)
leg.Draw()

allAngleCanv=[cLegend,cphi,ctheta,cphiA,cthetaA,cphiT,cthetaT,cangle,cphi2,cphiUV,cthetaUV,cdphiUV,cdthetaUV,cdangleUV,cdangleAZ]

histfile.cd()
#histfile=ROOT.TFile(args.basepath+args.histfile,'update')
for c in allAngleCanv:
    c.Write()
for c in allSizeCanv:
    c.Write()
for tit in c2dAngles:
    for k in c2dAngles[tit]:
        c2dAngles[tit][k].Update()
        c2dAngles[tit][k].Write()



if args.pfile!='None':
    allAngleCanv[0].SaveAs(args.pfile+'[')
    for c in allAngleCanv:
        c.SaveAs(args.pfile)
    for c in allSizeCanv:
        c.SaveAs(args.pfile)
    for tit in c2dAngles:
        for k in c2dAngles[tit]:
            c2dAngles[tit][k].SaveAs(args.pfile)

    allAngleCanv[0].SaveAs(args.pfile+']')



if not args.batch:
    u=raw_input("done?")
histfile.Close()