import  operator, numpy, os, sys
pandapath=os.environ.get('PANDAPATH')
sys.path.append(pandapath+'/macro/tpc/FOPI/python/argparse-1.2.1')
import argparse

parser=argparse.ArgumentParser(description='Plot stuff along a track')
parser.add_argument('filename',help='filename of the reco file',type=str)
parser.add_argument('--ev',help='which event to go to',type=int)
parser.add_argument('--events',help='number of events to process',type=int,default=-1)
parser.add_argument('--zcut',help='zcut',type=int,default=-80)
parser.add_argument('--data',help='analyse real data',action='store_true')
parser.add_argument('--pfile',help='the picture file',type=str,default="not")
parser.add_argument('--batch',help='batch mode',action='store_true')
args=parser.parse_args()

import ROOT
if args.batch:
    ROOT.gROOT.SetBatch(True)
ROOT.gROOT.ProcessLine(".x rootlogon.C") 
ROOT.gROOT.ProcessLine('gStyle->SetPalette(1)')
ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")')
ROOT.gROOT.LoadMacro("stlPYROOT.h+")


rfile=ROOT.TFile.Open(args.filename,"read")
if rfile.IsOpen()==ROOT.kTRUE:
    tree = rfile.Get("cbmsim")
if tree==None:
    print "No cbmsim found"
    if type(Rfile)==ROOT.TFile:
        Rfile.Close()
    exit()

tree.SetBranchStatus('*',0)
tree.SetBranchStatus('CutCosmics.*',1)

zSlices = (  12.0, 24.0, 36.0,  48.0, 60.0, 73.0 )

#Create Hisogramms---------------------------------------------------------
hcpt = ROOT.TH1D("hcpt", "Cluster per track",45, 0, 45 )
hcpcAll = ROOT.TH1D("hcpcAll", "Cluster per cm", 200, 0,2)
#hcpcVsZ = ROOT.TH2D("hcpcVsZ",'Cluster per cm Vs z',75,0,75,200,0,2)
hldt = ROOT.TH1D("hldt", "Length of the tracks", 40, 0, 40)
hnof = ROOT.TH1D("hnof", "number of tracks", 2,0,2)
hdbc = ROOT.TH2D("hdbc", "distance between cluster", 75,0,75,200, 0, 4)

hdpcz = ROOT.TH2D("hdpcz", "Digis per Cluster Vs z", 75, 0, 75, 45, 0, 45)
hppcz = ROOT.TH2D("hppcz", "Pads per Cluster Vs z", 75, 0, 75, 45, 0, 45)
hresz = ROOT.TH2D("hresz", "Residuals Vs z", 75, 0, 75, 2000, -0.1, 0.1)

hdpc=[]
hcpc=[]
hdpcm=[]
for i in range(len(zSlices)):
    hdpc.append(ROOT.TH1D("hdpc{0}".format(i), "Digis per Cluster, {0}st slice".format(i), 45, 0, 45))
    hcpc.append(ROOT.TH1D("hcpc{0}".format(i), "Cluster per cm, {0}st slice".format(i), 200, 0, 2))
    hdpcm.append(ROOT.TH1D("hdpcm{0}".format(i), "Digis per cm, {0}st slice".format(i), 40, 0, 40))


c1 = ROOT.TCanvas("c1", "Cluster 1", 700,1000) 
c2 = ROOT.TCanvas("c2", "Cluster 2", 700,1000)
c3 = ROOT.TCanvas("c3", "Cluster 3", 1000,700)
c4 = ROOT.TCanvas("c4", "Cluster 4", 1000,700)
c5 = ROOT.TCanvas("c5", "Cluster 5", 1000,700)
c6 = ROOT.TCanvas("c6", "Cluster 6", 700,1000)
c7 = ROOT.TCanvas("c7", "Cluster 7", 700,1000)
c8 = ROOT.TCanvas("c8", "Cluster 8", 700,1000)
c9 = ROOT.TCanvas("c9", "Cluster 9", 700,1000)

c1.Divide(1,2)
c2.Divide(1,2)
c3.Divide(3,2)
c4.Divide(3,2)
c5.Divide(3,2)
c6.Divide(1,2)
c7.Divide(1,2)
c8.Divide(1,2)
c9.Divide(1,2)

#-------------------------------------------------------------------------

#-------------------------------------------------------------------------
oldpos = ROOT.TVector3(0,0,0)
evnum = -1
dist = 0
summe_z = 0
digisize = 0
entries = 0
print "events:", tree.GetEntries()
for e in tree:
    evnum += 1
    trackindex = -1
    if evnum%100==0:
        print 'Events: {0}'.format(evnum)

    if e.CutCosmics.GetEntries() == 1:
        entries = 1
        hnof.Fill(entries)
    if evnum < args.ev:
        continue
    if args.events!=-1 and evnum>args.events:
        break

    for tr in e.CutCosmics:
        trackindex += 1
        resX=tr.GetResX()
        resY=tr.GetResY()
        resZ=tr.GetResZ()
        resU=tr.GetResU()
        resV=tr.GetResV()

        hitIds = tr.GetHitIDs()

        hitindex=0
        if resX.size()==0:
            continue

        dist = 0
        laenge=0
        cluster = 0
        dist1 = 0
        oldpos = ROOT.TVector3(0,0,0)
        for ihit in hitIds:
            cluster +=1
            hitindex +=1
            if hitindex>=resX.size():
                continue            
            #cl=e.TpcCluster.At(clid)
            pos=tr.GetHitPosition(hitindex)
            size = tr.GetClusterSize(hitindex)
            size2D=tr.Get2DClusterSize(hitindex)
        #    print "Y-Pos:", pos.Y()
     
            if pos.Z()<args.zcut:
                continue
            if hitindex>0:
                distance = tr.GetClusterDist(hitindex-1)
                hdbc.Fill(pos.Z(),distance)

            digisize += size
            res=tr.GetRes(hitindex)
            err=tr.GetSigXYZ(hitindex)
            hresz.Fill(pos.Z(),res.X())
           
            #Digis per Cluster - Fill Histogramm
            hdpcz.Fill(pos.Z(),size)
            hppcz.Fill(pos.Z(),size2D)
            for i in range(len(zSlices)):
                if pos.Z()<zSlices[i]:
                    hdpc[i].Fill(size)
                    break
            

            
        laenge = tr.GetLength()
        

        gesdigi = digisize
        digisize = 0
        
#print "******Event:", evnum
#print "Laenge:", laenge
#print "Anzahl der Cluster:", anzahl
#if laenge > 50:
#print  "*********zu lang:", evnum
        
        mw_pos=tr.GetMeanZ()

        if laenge != 0:
            cpc = len(hitIds)/laenge
        else:
            cpc = -1

        if laenge != 0:
            dpc = gesdigi/laenge
        else:
            dpc = 0
        
#Fill Histogramms----------------------------------------------
        hcpt.Fill(len(hitIds))
        hcpcAll.Fill(cpc)
        hldt.Fill(laenge)
#        hnof.Fill(entries)
        #hdpcz.Fill(pos.Z(), dpc)

        for i in range(len(zSlices)):
            if mw_pos<zSlices[i]:
                hcpc[i].Fill(cpc)
                hdpcm[i].Fill(dpc)
                break

#lable Histogramms-------------------------------------------------

hcpt.GetXaxis().SetTitle("count of cluster per track")
hcpt.GetYaxis().SetTitle("count")

hcpcAll.GetXaxis().SetTitle("cluster per cm [1/cm]")
hcpcAll.GetYaxis().SetTitle("count")

hldt.GetXaxis().SetTitle("length of the track [cm]")
hldt.GetYaxis().SetTitle("count")

hnof.GetYaxis().SetTitle("count")

for h in hcpc:
    h.GetXaxis().SetTitle("cluster per cm [1/cm]")
    h.GetYaxis().SetTitle("count")

for h in hdpc:
    h.GetXaxis().SetTitle("digis per cluster")
    h.GetYaxis().SetTitle("count")

for h in hdpcm:
    h.GetXaxis().SetTitle("digis per cm [1/cm]")
    h.GetYaxis().SetTitle("count")

hdbc.GetXaxis().SetTitle("Z (cm)")
hdbc.GetYaxis().SetTitle("Distance (cm)")

hdpcz.GetXaxis().SetTitle("Z (cm)")
hdpcz.GetYaxis().SetTitle("Clustersize")

hppcz.GetXaxis().SetTitle("Z (cm)")
hppcz.GetYaxis().SetTitle("Clustersize 2D")

hresz.GetXaxis().SetTitle("cm")
hresz.GetYaxis().SetTitle("cm")

#Output--------------------------------------------------------
print "Events:", evnum        
print 'done'
#Fit with gauss-------------------------------------------------
gaussfit=[]
#for h in hdpc:
    #gaussfit.append(ROOT.TF1("gausfit1", "landau", 0., 45.))
    #h.Fit(gaussfit[-1], "N+")

#"[0]*exp(-0.5*((x-[1])/[2])^2)"
#Draw Histogramms-----------------------------------------------

c1.cd(1)
hcpt.Draw()
c1.cd(2)
hcpcAll.Draw()
c1.Update()

c2.cd(1)
hldt.Draw()
c2.cd(2)
hnof.Draw()
c2.Update()

i=0
for h in hcpc:
    i+=1
    c3.cd(i)
    h.Draw()
c3.Update()

i=0
for h in hdpc:
    i+=1
    c4.cd(i)
    h.Draw()
    #gaussfit[i-1].Draw("SAME")
c4.Update()

i=0
for h in hdpcm: 
    i+=1
    c5.cd(i)
    h.Draw()
c5.Update()

c7.cd(1)
(c7.cd(1)).SetLogz()
hdpcz.Draw("colz")
c7.cd(2)
hdpcz.FitSlicesY(0,0,-1,0,"QNRL")
ROOT.gDirectory.Get("hdpcz_1").SetStats(0)
ROOT.gDirectory.Get("hdpcz_1").Draw()
c7.Update()

c9.cd(1)
(c9.cd(1)).SetLogz()
hppcz.Draw("colz")
c9.cd(2)
hppcz.FitSlicesY(0,0,-1,0,"QNRL")
ROOT.gDirectory.Get("hppcz_1").SetStats(0)
ROOT.gDirectory.Get("hppcz_1").Draw()
c9.Update()

c8.cd(1)
c8.cd(1).SetLogz()
hdbc.Draw('colz')
c8.cd(2)
c8.cd(2).SetLogz()
hresz.Draw("colz")
c8.Update()


#Graph mean----------------------------------------------------


#Get means
clmean=[]
for h in hcpc:
    clmean.append(h.GetMean())

digimean=[]
for h in hdpcm:
    digimean.append(h.GetMean())


#Create Graphs
clmeangraph = ROOT.TGraph()
digimeangraph = ROOT.TGraph()

#Fill Graphs
for i in range(len(zSlices)):
    clmeangraph.SetPoint(i,zSlices[i], clmean[i])
    digimeangraph.SetPoint(i,zSlices[i], digimean[i])

#Style of the Graphs
clmeangraph.SetTitle("cluster per cm - mean")
digimeangraph.SetTitle("digi per cm - mean")

clmeangraph.GetXaxis().SetTitle("z-pos [cm]")
clmeangraph.GetYaxis().SetTitle("mean [1/cm]")

clmeangraph.SetMarkerStyle(20)
clmeangraph.SetMarkerSize(1)
clmeangraph.SetMarkerColor(2)

digimeangraph.GetXaxis().SetTitle("z-pos [cm]")
digimeangraph.GetYaxis().SetTitle("mean [1/cm]")

digimeangraph.SetMarkerStyle(20)
digimeangraph.SetMarkerSize(1)
digimeangraph.SetMarkerColor(2)

Gerade = ROOT.TF1("Gerade", "[0]*x + [1]", -60., 10.);
digimeangraph.Fit(Gerade,"N")
steigung=Gerade.GetParameter(0)
yAbsch=Gerade.GetParameter(1)

text=ROOT.TLatex()
text.SetTextColor(ROOT.kBlack)
text.SetTextSize(0.03)
text.SetNDC()
y=0.8
ystep=0.05
#text.DrawLatex(0.15,y,str(steigung))


#Draw Graphs
c6.cd(1)
clmeangraph.Draw("AP")
c6.cd(2)
digimeangraph.Draw("AP")
Gerade.Draw("SAME")
text.DrawLatex(0.15,y,"y="+str(steigung)+"x"+"+"+str(yAbsch))
c6.Update()

#Save Pictures-----------------------------------------------
if args.pfile!="not":
    c1.SaveAs(args.pfile+'(')
    c2.SaveAs(args.pfile)
    c3.SaveAs(args.pfile)
    c4.SaveAs(args.pfile)
    c5.SaveAs(args.pfile)
    c6.SaveAs(args.pfile)
    c7.SaveAs(args.pfile)
    c9.SaveAs(args.pfile)
    c8.SaveAs(args.pfile+')')

#------------------------------------------------------------
u='w'
while u!='q' and not args.batch:
    u = raw_input("Exit with q")