import ROOT,math,os,sys
from array import array
import copy
pandapath=os.environ.get('PANDAPATH')
sys.path.append(pandapath+'/macro/tpc/FOPI/python/argparse-1.2.1')
sys.path.append(pandapath+'/macro/tpc/FOPI/mberger')
from PyChainMaker import PyChainMaker
#from ROOT import std,TMath

colors = [ROOT.kBlack,ROOT.kRed+1,ROOT.kGreen+3,ROOT.kBlue+1,ROOT.kMagenta+1,ROOT.kOrange+1,ROOT.kViolet-7]

def calcSigErr2(s1,ds1,c1,dc1,s2,ds2,c2,dc2) :

    divider=c1*s1+c2*s2
    subst=(c1*(s1**2) + c2*(s2**2)) / (divider**2)
    df_ds1 = ( (2*c1*s1)/divider - (c1*subst) )**2
    df_ds2 = ( (2*c2*s2)/divider - (c2*subst) )**2
    df_dc1 = ( (s1**2)/divider - (s1*subst) )**2
    df_dc2 = ( (s2**2)/divider - (s2*subst) )**2 
    
    sigerr=math.sqrt( df_ds1*(ds1**2) + df_ds2*(ds2**2) + df_dc1*(dc1**2) + df_dc2*(dc2**2) )
    
    return sigerr
    ''' old calculation (changed h to c)
    divider=(h1*s1+h2*s2)
    divider *= divider

    ddh1  = math.fabs(h2*s1*s1*s2-h2*s2*s2*s1)
    ddh1 /= divider
    ddh1 *= dh1

    dds1  = math.fabs(h1*h1*s1*s1+2*h1*h2*s1*s2-h2*h1*s2*s2)
    dds1 /=divider
    dds1 *= ds1

    ddh2  = math.fabs(h1*s1*s2*s2-h1*s1*s1*s2*s2)
    ddh2 /= divider
    ddh2 *= dh2

    dds2  = math.fabs(h2*h2*s2*s2+2*h2*h1*s1*s2 - h1*h2*s1*s1)
    dds2 /= divider
    dds2 *= ds2

    meanerr = math.sqrt(ddh1+ddh2+dds1+dds2)

    return meanerr
    '''
def calcMeanErr(s1,ds1,c1,dc1,m1,dm1,s2,ds2,c2,dc2,m2,dm2):
    
    divider=c1*s1+c2*s2
    subst= (c1*m1*s1 + c2*m2*s2)/(divider**2)
    
    df_dm1 = (c1**2 + s1**2)/(divider**2)
    df_dm2 = (c2**2 + s2**2)/(divider**2)
    df_ds1 = ((c1*m1)/divider - c1*subst )**2
    df_ds2 = ((c2*m2)/divider - c2*subst )**2
    df_dc1 = ((m1*s1)/divider - s1*subst )**2
    df_dc2 = ((m2*s2)/divider - s2*subst )**2
    
    meanerr = math.sqrt( df_dm1*(dm1**2) + df_dm2*(dm2**2) + df_ds1*(ds1**2) + df_ds2*(ds2**2) + df_dc1*(dc1**2) + df_dc2*(dc2**2) )
    return meanerr
    

def fitres(h,_debug=0,_df=0,_startval=[],_frange=[-10000.,10000.],_factors=[1.,1.,1.,5.]) :
    debug=_debug
    df=_df
    startval=[]
    for val in _startval:
        startval.append(val)
    frange=[]
    for val in _frange:
        frange.append(val)
    factors=[]
    for val in _factors:
        factors.append(val)
            
    data={}
    #if h.GetName().find("hpull_cl_a1VSz")!=-1:
    #    debug=1
    
    if debug:
        print h.GetTitle()
        
    histo=h.Clone()
    histo.SetFillColor(ROOT.kAzure-8)

    if frange[0]==-10000:
        #print 'change range: ',frange[0],'->',histo.GetMean()-5.*histo.GetRMS(),frange[1],'->',histo.GetMean()+5.*histo.GetRMS()
        frange[0]=histo.GetMean()-5.*histo.GetRMS()
        frange[1]=histo.GetMean()+5.*histo.GetRMS()
     
    testfit = ROOT.TF1("testfitY"+str(1),"gaus",frange[0],frange[1])
    if len(startval)==3:
        testfit.SetParameter(0,startval[0])
        testfit.SetParLimits(0,0,startval[0]*1.2)
        testfit.SetParameter(1,startval[1])
        testfit.SetParLimits(1,startval[1]-3*startval[2],startval[1]+3*startval[2])
        testfit.SetParameter(2,startval[2])
        testfit.SetParLimits(2,-2*startval[2],2*startval[2])
    else:
        testfit.SetParameter(1,histo.GetMean())
        testfit.SetParameter(2,histo.GetRMS())
        startval.append(-1)
        startval.append(histo.GetMean())
        startval.append(histo.GetRMS())
        
    if debug:
        print "startval:",startval[0],startval[1],startval[2],frange[0],frange[1]
        print 'bounds 0:',0,startval[0]*1.2
        print 'bounds 1:',startval[1]-3*startval[2],startval[1]+3*startval[2]
        print 'bounds 2:',-2*startval[2],2*startval[2]
        
    testfit.SetLineColor(2)
    
    if debug:
        print 'testfit'
    
    fitopt="QRN+"
    
    if debug:
        fitopt="R+"
    
    if debug:
        print 'ranges testfit:',frange[0],frange[1],histo.GetRMS(),histo.GetMean()

    histo.Fit(testfit, fitopt, "", frange[0],frange[1])
    
    data['tconst']    =testfit.GetParameter(0)    
    data['tmean']     =testfit.GetParameter(1)
    data['tsig']      =testfit.GetParameter(2)
    data['tconsterr'] =testfit.GetParError(0)
    data['tmeanerr']  =testfit.GetParError(1)
    data['tsigerr']   =testfit.GetParError(2)
    
    
    if debug:
        print 'done testfit'

    rfit = ROOT.TF1("fitfuncY"+str(1),"gaus + gaus(3)",frange[0],frange[1])
    rfit1 = ROOT.TF1("gaus1","gaus",frange[0],frange[1])
    rfit2 = ROOT.TF1("gaus2","gaus",frange[0],frange[1])

    rfit.SetNpx(1000)
    rfit.SetParName(0,"const1")
    rfit.SetParName(1,"mean1")
    rfit.SetParName(2,"sigma1")
    rfit.SetParName(3,"const2")
    rfit.SetParName(4,"mean2")
    rfit.SetParName(5,"sigma2")
     
    tconst = [None]*2
    tmean  = [None]*2
    tsigma = [None]*2

    tconst[0] = tconst[1] = testfit.GetParameter(0)
    tmean[0]  = tmean[1]  = testfit.GetParameter(1)
    tsigma[0] = tsigma[1] = testfit.GetParameter(2)

    #broad gaus
    rfit.SetParameter(0,tconst[0]/1.5)
    rfit.SetParameter(1,tmean[0])
    rfit.SetParameter(2,tsigma[0]*2.5)
    rfit.SetParLimits(2,tsigma[0]*factors[2], tsigma[0]*factors[3])
    
    #narrow gaus
    rfit.SetParameter(3, tconst[1])
    rfit.SetParameter(4, tmean[0])
    rfit.SetParameter(5, tsigma[1]*0.5)
    rfit.SetParLimits(5, 0.0*factors[0], tsigma[1]*factors[1])

    rfit.SetLineWidth(1)
    histo.Fit(rfit, fitopt,"", -3*histo.GetRMS(),3*histo.GetRMS())
    
    '''
    if(rfit.GetParError(2)/rfit.GetParameter(2) > 0.1):
        for i in range(1):
            histo.Fit(rfit, fitopt,"", histo.GetMean()-3*histo.GetRMS(),histo.GetMean()+3*histo.GetRMS())
        
    if(rfit.GetParError(2)/rfit.GetParameter(2) > 0.1):
        histo.Fit(rfit, fitopt,"", histo.GetMean()-3*histo.GetRMS(),histo.GetMean()+3*histo.GetRMS())   
    '''    
    if debug:
        print 'ranges fit:',histo.GetMean()-3*histo.GetRMS(),histo.GetMean()+3*histo.GetRMS()

    if rfit.GetNDF()>0:
        scaleu=rfit.GetChisquare()/rfit.GetNDF()
        scaled=1/scaleu
    else:
        scaleu=2.
        scaled=0.5

    if df==1:
        for j in range(2) :
            tconst[j]=rfit.GetParameter(j*2)
            tmean[j]=rfit.GetParameter(j*2+1)
            tsigma[j]=rfit.GetParameter(j*2+2)

        #broad gaus
        rfit.SetParameter(0,tconst[0])
        rfit.SetParameter(2,tsigma[0])
        rfit.SetParLimits(2,tsigma[0]*scaled, tsigma[0]*scaleu)
    
        #narrow gaus
        rfit.SetParameter(3, tconst[1])
        rfit.SetParameter(5, tsigma[1])
        rfit.SetParLimits(5, tsigma[1]*scaled, tsigma[1]*scaleu)
    
        rfit.SetLineWidth(2)
        histo.Fit(rfit, "+", "", frange[0],frange[1])

    nconst=rfit.GetParameter(3)
    dnconst=rfit.GetParError(3)
    bconst=rfit.GetParameter(0)
    dbconst=rfit.GetParError(0)
    
    nmean=rfit.GetParameter(4)
    dnmean=rfit.GetParError(4)
    bmean=rfit.GetParameter(1)
    dbmean=rfit.GetParError(1)

    nsig=rfit.GetParameter(5)
    dnsig=rfit.GetParError(5)
    bsig=rfit.GetParameter(2)
    dbsig=rfit.GetParError(2)

    nint=abs(nsig*nconst)
    bint=abs(bsig*bconst)

    meanerr = 0
    if nint!=0 or bint!=0 :
        data['csig'] = (nint*nsig+bint*bsig)/(nint+bint)
        data['csigerr'] = calcSigErr2(abs(nsig), abs(dnsig), nconst, dnconst, abs(bsig), abs(dbsig), bconst, dbconst)
        data['cmean']= (nint*nmean+bint*bmean)/(nint+bint)
        data['cmeanerr'] = calcMeanErr(abs(nsig), abs(dnsig), nconst, dnconst, nmean, dnmean, abs(bsig), abs(dbsig), bconst, dbconst, bmean, dbmean)
        data['cconst']=nconst+bconst
        data['cconsterr']=math.sqrt(dnconst**2+dbconst**2)
    else:
        data['csig'] = -9999
        data['csigerr']=0.
        data['cmean']= -9999
        data['cmeanerr']=0
        data['cconst']=0
        data['cconsterr']=0


    data['nsig']    = nsig
    data['nsigerr'] = dnsig
    data['bsig']    = bsig
    data['bsigerr'] = dbsig

    data['nconst']   = nconst
    data['nconsterr']= dnconst
    data['bconst']   = bconst
    data['bconsterr']= dbconst

    data['nmean']   = nmean 
    data['nmeanerr']= dnmean
    data['bmean']   = bmean
    data['bmeanerr']= dbmean
        
    if debug:
        print abs(nsig), abs(dnsig), nconst, dnconst, nmean, dnmean, abs(bsig), abs(dbsig), bconst, dbconst, bmean, dbmean,data['csigerr']
        print 
    
    if rfit.GetNDF()>0:
        data['chindf']=rfit.GetChisquare()/rfit.GetNDF()
    else:
        data['chindf']=-1

    if debug:
        canv=ROOT.gROOT.FindObject("debug_canv")
        if type(canv)!=ROOT.TCanvas:
            canv=ROOT.TCanvas("debug_canv",h.GetTitle())
        leg=ROOT.TLegend(0.7,.7,0.9,.9,"","NDC")
        canv.cd()
        histo.Draw()
        rfit1.SetParameter(0,rfit.GetParameter(0))
        rfit1.SetParameter(1,rfit.GetParameter(1))
        rfit1.SetParameter(2,rfit.GetParameter(2))
        rfit1.SetLineColor(ROOT.kRed)
        leg.AddEntry(rfit1,"Broad")
        rfit1.DrawCopy("same")
        rfit2.SetParameter(0,rfit.GetParameter(3))
        rfit2.SetParameter(1,rfit.GetParameter(4))
        rfit2.SetParameter(2,rfit.GetParameter(5))
        rfit2.SetLineColor(ROOT.kBlue)
        leg.AddEntry(rfit2,"Narrow")
        rfit2.DrawCopy("same")
        testfit.SetLineColor(ROOT.kGreen)
        testfit.Draw("same")
        leg.AddEntry(testfit,"testfit")
        leg.Draw()
        print "*******************************************"
        print tsigma[0]*10.2, tsigma[1]*0.05, tsigma[1]*1.5
        print "*******************************************"
        #drawtext(["broad"+str(rfit.GetParameter(2)),"small"+str(rfit.GetParameter(5))])
        drawresdata(data,"sigma chi/ndf S/B ngaus bgaus cgaus tgaus")
        canv.Update()
        nix=raw_input("blahorst")

    return data

def drawresdata(data,opt="sigma chi/ndf S/B ngaus bgaus cgaus",unit="#mum",roundto=1,yoff=0.8,ystep=0.05):
    text=ROOT.TLatex()
    text.SetTextColor(ROOT.kBlack)
    text.SetTextSize(0.05)
    text.SetNDC()
    opts=opt.split()
    if opt.find('times')!=-1:
        text.SetTextFont(132)
    if opt.find('tsize')!=-1:
        for o in opts:
            if o.find('tsize')!=-1:
                size=float(o.split('=')[1])
                text.SetTextSize(size)
        
    y=yoff
    x=0.24
    #ystep=0.05
    y+=ystep
    if opt.find("title"):
        for o in opts:
            if o.find("title")!=-1:
                title=o.split('=')[1]
                y-=ystep
                text.DrawLatex(x,y,title)
    
    if opt.find("sigma")!=-1:
        if opt.find("ngaus")!=-1:
            y-=ystep
            text.DrawLatex(x,y,"#color[{3}]{{#sigma_{{narrow}}={0:8.3f} #pm{1:6.3f} {2} }}".format(data['nsig'],data['nsigerr'],unit,ROOT.kGreen+2))
        if opt.find("bgaus")!=-1:
            y-=ystep
            text.DrawLatex(x,y,"#color[2]{{#sigma_{{broad}} ={0:8.3f} #pm{1:6.3f} {2} }}".format(data['bsig'],data['bsigerr'],unit))
        if opt.find("cgaus")!=-1:
            y-=ystep
            text.DrawLatex(x,y,"#color[{3}]{{#sigma_{{mean}} ={0:8.3f} #pm{1:6.3f} {2} }}".format(data['csig'],data['csigerr'],unit,ROOT.kAzure+7))
        if opt.find("tgaus")!=-1:
            y-=ystep
            text.DrawLatex(x,y,"#color[{3}]{{#sigma ={0:8.3f} #pm{1:6.3f} {2} }}".format(data['tsig'],data['tsigerr'],unit,ROOT.kBlack))
            
    if opt.find("mean")!=-1:
        y-=ystep
        text.DrawLatex(x,y,"Mean={0:3.5f}".format(data['nmean']))
    if opt.find('chi/ndf')!=-1:
        y-=ystep
        text.DrawLatex(x,y,"#chi^{{2}}/ndf={0:8.5f}".format(data['chindf']))
    if opt.find("S/B")!=-1:
        y-=ystep
        if (data['bsig']*data['bconst'])>0:
            text.DrawLatex(x,y,"S/B="+str(round((data['nsig']*data['nconst'])/(data['bsig']*data['bconst']),2*roundto)))
            
    if opt.find("const")!=-1:
        y-=ystep
        text.DrawLatex(x,y,"c_{{narrow}}={0:8.3f}#pm{1:8.3f}".format(data['nconst'],data['nconsterr']))
        y-=ystep
        text.DrawLatex(x,y,"c_{{broad}}={0:8.3f}#pm{1:8.3f}".format(data['bconst'],data['bconsterr']))
        
    
    return

def drawGaussians(data,opt="ngaus bgaus cgaus"):
    rfit1 = ROOT.TF1("gaus1","gaus",-10000.,10000.)
    if data.get('bconst',-1)==-1 or data.get('bmean',-1)==-1 or data.get('bsig',-1)==-1:
        return
    if opt.find("bgaus")!=-1:
        rfit1.SetParameter(0,data['bconst'])
        rfit1.SetParameter(1,data['bmean'])
        rfit1.SetParameter(2,data['bsig'])
        rfit1.SetLineColor(ROOT.kRed)
        rfit1.DrawCopy("same")
   
    if opt.find("ngaus")!=-1:
        rfit2 = ROOT.TF1("gaus2","gaus",-10000.,10000.) 
        rfit2.SetParameter(0,data['nconst'])
        rfit2.SetParameter(1,data['nmean'])
        rfit2.SetParameter(2,data['nsig'])
        #rfit2.SetLineColor(ROOT.kBlue)
        rfit2.SetLineColor(ROOT.kGreen+2)
        rfit2.DrawCopy("same")

    if opt.find("cgaus")!=-1:
        rfit3 = ROOT.TF1("gaus3","gaus+gaus(3)",-10000.,10000.)
        rfit3.SetParameter(0,data['bconst'])
        rfit3.SetParameter(1,data['bmean'])
        rfit3.SetParameter(2,data['bsig'])
        rfit3.SetParameter(3,data['nconst'])
        rfit3.SetParameter(4,data['nmean'])
        rfit3.SetParameter(5,data['nsig'])
        #rfit3.SetLineColor(ROOT.kGreen+2)
        rfit3.SetLineColor(1)
        rfit3.SetLineStyle(2)
        rfit3.SetLineWidth(3)
        rfit3.DrawCopy("same")
        
    if opt.find("tgaus")!=-1:
        rfit4 = ROOT.TF1("gaus4","gaus",-10000.,10000.)
        rfit4.SetParameter(0,data['tconst'])
        rfit4.SetParameter(1,data['tmean'])
        rfit4.SetParameter(2,data['tsig'])
        rfit4.SetLineColor(ROOT.kBlack)
        rfit4.SetLineColor(1)
        rfit4.DrawCopy("same")
    return

def drawPrelim() :
    prelim=ROOT.TLatex()
    prelim.SetTextColor(ROOT.kGray)
    prelim.SetTextSize(0.07)
    prelim.SetTextAngle(20)
    prelim.SetNDC()
    prelim.DrawLatex(0.38,0.45,"preliminary")
    return

def drawDiff(noB,ar,ne):
    if noB==1:
        if ar==1:
            DIFFX = 0.0227*10000   #(mu m/ mu s) without B-Field
    else:
        if ar==1:
            DIFFX = 0.020534*10000 #(mu m/ mu s)Argon with B-Field
        if ne==1:
            DIFFX = 0.019833*10000 #(mu m/ mu s)Neon with B-Field

    diffT = ROOT.TF1("diffT","[0]*TMath::Sqrt(x)",0,75)
    diffT.SetNpx(1000)
    diffT.SetLineWidth(2)
    diffT.SetLineStyle(4)
    diffT.SetParameter(0,DIFFX)
    diffT.SetTitle("")

    return diffT

def sort_inner(inner):
#inner is each inner list in the list of lists to be sorted
#(here item at index 1 of each inner list is to be sorted)
    return inner[0]

def sort_inner2(inner):
    return inner[2]

def Drawh(hist,hists,legar,canv):
    
    col=0
    colors = [ROOT.kBlack,ROOT.kRed+1,ROOT.kRed-7,ROOT.kGreen+3,ROOT.kBlue+1,ROOT.kMagenta+1,ROOT.kOrange+1,ROOT.kViolet-7]
    leg=ROOT.TLegend(0.7,.7,0.9,.9,"","NDC")
    canv.cd()
    maxbin=0
    for h in hist:
        col+=1
        if col < len(colors):
            color=colors[col]
        else:
            print "no color available"
            color=1
#        hist[h].SetTitle(hists)
        if hist[h].GetMaximum() > maxbin:
            maxbin=hist[h].GetMaximum()
        hist[h].SetLineColor(color)
        if col == 1:
            hist[h].GetXaxis().SetTitle(axlabel[0].get(hists,""))
            hist[h].GetYaxis().SetTitle(axlabel[1].get(hists,""))
            hist[h].Draw()
            first=hist[h]
        else:
            hist[h].Draw("same")
        leg.AddEntry(hist[h],h)
    first.SetMaximum(maxbin*1.02)
    leg.SetFillColor(0)
    leg.Draw()
    legar.append(leg)
#    else:
#        hist.Draw()

    return 

def Divideh(h1,h2,name,title):
    temphist=h1.Clone(name)
    temphist.SetTitle(title)
    temphist.Divide(h2)
    return temphist

def Addh(h1,h2,factor,name,title=""):
    temphist=h1.Clone(name)
    temphist.SetTitle(title)
    temphist.Add(h2,factor)
    return temphist

def fillcutred(cut,event,hist):
#    hist=histsdict['cut reduction']
    hist.Fill(1,cut)

    if event['clmeanamp']<400:
        hist.Fill(2,cut)
    if event['theta']<10 or event['theta']>170:
        hist.Fill(3,cut)
    if abs(event['phi'])<2:
        hist.Fill(4,cut)
    if event['clmeanamp']>400 and (event['theta']<10 or event['theta']>170):
        hist.Fill(5,cut)
    return

def plotfullevent(thevent):
    trackcounter=0
    digamp=[]
    digpos=[]
    canvamp=[]
    canvpos=[]
    lines=[]
    lines2d=[]
    lines2d2=[]
    for tr in range(len(thevent)):
        trackcounter+=1
        vec=ROOT.TVector3(1.,1.,1.)
        vec.SetPhi(thevent[tr]['phi'])
        vec.SetTheta(thevent[tr]['theta'])
        vec.SetMag(20.)
        line=ROOT.TPolyLine3D(2,"")
        x2=thevent[tr]['hits'][len(thevent[tr]['hits'])-1]['pos'][0]
        y2=thevent[tr]['hits'][len(thevent[tr]['hits'])-1]['pos'][1]
        z2=thevent[tr]['hits'][len(thevent[tr]['hits'])-1]['pos'][2]
        x1=thevent[tr]['hits'][0]['pos'][0]
        y1=thevent[tr]['hits'][0]['pos'][1]
        z1=thevent[tr]['hits'][0]['pos'][2]
        line.SetPoint(0,x1,y1,z1)
        line.SetPoint(1,x2,y2,z2)
        lines2d.append(ROOT.TLine(x1,y1,x2,y2))
        lines2d2.append(ROOT.TLine(z1,y1,z2,y2))
        lines2d[trackcounter-1].SetLineColor(trackcounter)
        lines2d2[trackcounter-1].SetLineColor(trackcounter)

        print "appending line, size:", len(lines2d), len(lines2d2)                    
                    #tcn2.cd()
                    #hits3d.Draw()
                    #vec.Draw("same")
                    #line.Draw("same");
                    #tcn2.Update()

        canvamp.append(ROOT.TCanvas("samples"+str(trackcounter),"samples"+str(trackcounter),720,10,700,500))
        canvpos.append(ROOT.TCanvas("samples2D"+str(trackcounter),"samples2D"+str(trackcounter),720,520,700,500))
        print int(math.ceil(math.sqrt(len(thevent[tr]['hits']))))
        canvamp[trackcounter-1].Divide(int(math.ceil(math.sqrt(len(thevent[tr]['hits'])))),int(math.ceil(math.sqrt(len(thevent[tr]['hits'])))))
        canvpos[trackcounter-1].Divide(int(math.ceil(math.sqrt(len(thevent[tr]['hits'])))),int(math.ceil(math.sqrt(len(thevent[tr]['hits'])))))
        clcount=1
        hcount=0
#        his=[]
#        his2d=[]
        digpos.append([])
        digamp.append([])
        for clust in thevent[tr]['hits']:
            histsdict['hxy'].Fill(clust['pos'][0],clust['pos'][1])#,clust['res'][0])
            histsdict['hyz'].Fill(clust['pos'][2],clust['pos'][1])#,clust['res'][0])
            canvamp[trackcounter-1].cd(clcount)
            clcount+=1
            same=""
            i=1
            digpos[trackcounter-1].append(ROOT.TH2D("cl2d"+str(clcount-1)+"_"+str(hcount),"cluster "+str(clcount-1),200,-16.,16.,200,-16.,16.))

            for dig in clust['digis']:
                digamp[trackcounter-1].append(ROOT.TH1D("cl"+str(clcount-1)+"_"+str(hcount),"cluster "+str(clcount-1),511,0.,511.))

                digamp[trackcounter-1][hcount].Fill(dig['time'],dig['amp'])
                digpos[trackcounter-1][clcount-2].Fill(mapping[int(dig['pad'])][0],mapping[int(dig['pad'])][1])

                digamp[trackcounter-1][hcount].SetLineColor(i)
                digamp[trackcounter-1][hcount].Draw(same)
                same="same"
                hcount+=1
                i+=1
                print "pads:",int(dig['pad']), mapping[int(dig['pad'])][0],mapping[int(dig['pad'])][1]
            canvpos[trackcounter-1].cd(clcount-1)
            digpos[trackcounter-1][clcount-2].Draw("colz")
            print ""

        print "phi angle:",thevent[tr]['phi']
        print "theta angle:",thevent[tr]['theta']
        print "momentum:",thevent[tr]['mom']
        print "chi2",thevent[tr]['chi2']," ndf",thevent[tr]['ndf']
        print "clmeanamp:",thevent[tr]['clmeanamp']
        print "event number:", evt
        print "vector",vec.X(),vec.Y(),vec.Z()
        print "cluster1",thevent[tr]['hits'][0]['pos'][0],thevent[tr]['hits'][0]['pos'][1],thevent[tr]['hits'][0]['pos'][2]
        print "cluster2",clust['pos'][0],clust['pos'][1],clust['pos'][2]
        print "trackcounter: ",trackcounter

        tcn.cd(1)
        histsdict['hxy'].Draw("colz")
        tcn.cd(2)
        histsdict['hyz'].Draw("colz")
        tcn.Update()

        canvamp[trackcounter-1].Update()
        canvpos[trackcounter-1].Update()
#        if trackcounter==len(event):
    for l in range(len(lines2d)):
        tcn.cd(1)
        lines2d[l].Draw("same")
        tcn.cd(2)
        lines2d2[l].Draw("same")
        tcn.Update()                            
    raw_input("next event")

#    tcn3.Delete()
#    tcn4.Delete()
#    for i in range(len(his)):
#        his[i].Delete()
#    del his
    histsdict['hxy'].Reset()
    histsdict['hyz'].Reset()

    return

def parseRuns(theruns):
    runList=[]
    colindex = theruns.find(",")
    if colindex > 0:
        runs = theruns.split(",")
        for i in range(len(runs)) :
            dashIndex = runs[i].find("-")
            if dashIndex > 0 :
                runs[i] = runs[i].split("-")
                for j in range(int(runs[i][0]), int(runs[i][1])+1) :
                    runList.append(j)
            if dashIndex < 0 and len(runs) > 1 :
                runList.append(int(runs[i]))
    else :
        dashIndex = theruns.find("-")
        if dashIndex > 0 :
            runs = theruns.split("-")
            runList = range(int(runs[0]), int(runs[1])+1)
        else:
            runList = [theruns]
    return runList 

def set_palette(name="palette", ncontours=99):
    """Set a color palette from a given RGB list
    stops, red, green and blue should all be lists of the same length
    see set_decent_colors for an example"""

    if name == "gray" or name == "grayscale":
        stops = [0.00, 0.34, 0.61, 0.84, 1.00]
        red   = [1.00, 0.84, 0.61, 0.34, 0.00]
        green = [1.00, 0.84, 0.61, 0.34, 0.00]
        blue  = [1.00, 0.84, 0.61, 0.34, 0.00]
    # elif name == "whatever":
        # (define more palettes)
    else:
        # default palette, looks cool
        print 'defining new palette'
        stops = [0.00, 0.34, 0.61, 0.84, 1.00]
        red   = [0.00, 0.00, 0.87, 1.00, 0.51]
        green = [0.00, 0.81, 1.00, 0.20, 0.00]
        blue  = [0.51, 1.00, 0.12, 0.00, 0.00]

    s = array('d', stops)
    r = array('d', red)
    g = array('d', green)
    b = array('d', blue)

    npoints = len(s)
    ROOT.TColor.CreateGradientColorTable(npoints, s, r, g, b, ncontours)
    ROOT.gStyle.SetNumberContours(ncontours)
    ROOT.gROOT.ForceStyle()
    print "defined new palette",name,ncontours

def setRange(h,mini=0,maxi=0):
    for i in range(h.GetNbinsX()+1):
        for j in range(h.GetNbinsY()+1):
            if (maxi!=mini):
                cont=h.GetBinContent(i,j)
                if (cont>maxi):
                    h.SetBinContent(i,j,maxi)
                if (cont<mini):
                    h.SetBinContent(i,j,mini)
    if(mini!=maxi):
        h.GetZaxis().SetRangeUser(mini,maxi)
    return

def createGraph(hists,points,toplot="mean",color=1,style=20):
    graph=ROOT.TGraph()
    counter=-1
    mini=100000
    maxi=0
    for h in hists:
        counter+=1
        val=5
        if toplot=="mean":
            val=h.GetMean()
        if toplot=="int2":
            nb=h.GetNbinsX()
            val=0
            for i in range(1,nb):
                val+=h.GetBinContent(i)*h.GetXaxis().GetBinCenter(i)
        maxi=max(val,maxi)
        mini=min(val,mini)
        graph.SetPoint(counter,points[counter],val)

    graph.SetMarkerStyle(style)
    graph.SetMarkerColor(color)
    graph.SetMaximum(maxi)
    graph.SetMinimum(mini)
    return graph

def setGraphLayout(g,xtit,ytit,marker=20,size=1,col=1):

    return g

def getColors():
    col=[]
    start=-9
    stop=10
    step=3
    for i in range(start,stop,step):
        col.append(ROOT.kViolet+i)
    for i in range(start,stop,step):
        col.append(ROOT.kOrange+i)
    for i in range(start,stop,step):
        col.append(ROOT.kSpring+i)
    for i in range(start,stop,step):
        col.append(ROOT.kAzure+i)
    for i in range(start,stop,step):
        col.append(ROOT.kPink+i)
    for i in range(start,stop,step):
        col.append(ROOT.kTeal+i)
    return col

def getColors2():
    col=[]
    start=-9
    stop=10
    step=3
    for i in range(start,stop,step):
        col.append(ROOT.kViolet+i)
        col.append(ROOT.kMagenta+i)
        col.append(ROOT.kOrange+i)
        col.append(ROOT.kAzure+i)
        col.append(ROOT.kSpring+i)
        col.append(ROOT.kTeal+i)
    return col
    
def getSlices(h,projs,opt=0,func=0,ropt='c',rebin=-1):
    graph={}

    if opt==0:
        tmp=fitSlices(h,graph,func,0,{},{},3,ropt,rebin)

    elif opt==1:
        tmp=getSlicesHist(h,graph)

    if projs!=None:
        for p in tmp:
            projs.append(copy.deepcopy(p))

    return graph

def fitSlices(hist,graph,func=0,fnum=0,mini={},maxi={},binwidth=3,ropt='c',rebin=-1):
    if func==1:
        func=ROOT.TF1("landau","landau",0,hist.GetRMS(2)*5)
        func.SetParLimits(1,0,10000)

    nbins=hist.GetNbinsX()
    projs=[]
    #graph.append({})
    graph['RMS']=ROOT.TGraphErrors()
    setGraph(graph['RMS'],'Sigma '+hist.GetTitle(),fnum+1)
    graph['Mean']=ROOT.TGraphErrors()
    setGraph(graph['Mean'],'Mean '+hist.GetTitle(),fnum+1)
    graph['Const']=ROOT.TGraphErrors()
    setGraph(graph['Const'],'Const '+hist.GetTitle(),fnum+1)

    pcounter=0
    i=binwidth
    while i<nbins:
        projs.append(hist.ProjectionY(hist.GetName()+str(i),i-binwidth,i))
        projs[-1].SetTitle(hist.GetTitle()+'_'+str(i))
        if rebin!=-1:
            #print 'rebinning',rebin
            projs[-1].Rebin(rebin)
        if projs[-1].GetEntries()==0:
            i+=1
            continue
        if func==1:
            projs[-1].Fit(func,"QRN+")
            data={}
            data[ropt+'sig']=func.GetParameter(2)
            data[ropt+'sigerr']=func.GetParError(2)
            data[ropt+'mean']=func.GetParameter(1)
            data[ropt+'meanerr']=func.GetParError(1)
            data[ropt+'const']=func.GetParameter(0)
            data[ropt+'consterr']=func.GetParError(0)
        else:
            data=fitres(projs[-1])

        graph['RMS'].SetPoint(pcounter,i-int(binwidth/2.),data[ropt+'sig'])
        graph['RMS'].SetPointError(pcounter,float(1./2.),data[ropt+'sigerr'])
        graph['Mean'].SetPoint(pcounter,i-int(binwidth/2.),data[ropt+'mean'])
        graph['Mean'].SetPointError(pcounter,float(1/2.),data[ropt+'meanerr'])
        graph['Const'].SetPoint(pcounter,i-int(binwidth/2.),data[ropt+'const'])
        graph['Const'].SetPointError(pcounter,float(1/2.),data[ropt+'consterr'])

        mini['RMS']=min(mini.get('RMS',1000),data[ropt+'sig'])
        maxi['RMS']=max(maxi.get('RMS',0),data[ropt+'sig'])
        mini['Mean']=min(mini.get('Mean',1000),data[ropt+'mean'])
        maxi['Mean']=max(maxi.get('Mean',0),data[ropt+'mean'])
        pcounter+=1
        i+=1
        #if (i>=nbins):
        #    break
    graph['RMS'].SetName(hist.GetName()+'_Sigma')
    graph['Mean'].SetName(hist.GetName()+'_Mean')
    graph['Const'].SetName(hist.GetName()+'_Const')

    return projs

def getSlicesHist(hist,graph,fnum=0,mini={},maxi={},binwidth=2):
    nbins=hist.GetNbinsX()
    projs=[]
    #graph.append({})
    graph['RMS']=ROOT.TGraphErrors()
    setGraph(graph['RMS'],'RMS '+hist.GetTitle(),1,getColors()[fnum+1])
    graph['Mean']=ROOT.TGraphErrors()
    setGraph(graph['Mean'],'Mean '+hist.GetTitle(),1,getColors()[fnum+1])
    pcounter=0
    binwidth
    i=binwidth
    while i<nbins:
        projs.append(hist.ProjectionY(hist.GetName()+str(i),i-binwidth,i))
        if (projs[-1].GetEntries()==0):
            i+=1
            continue
        graph['RMS'].SetPoint(pcounter,i-binwidth+1,projs[-1].GetRMS())
        graph['RMS'].SetPointError(pcounter,float(1./2.),projs[-1].GetRMSError())
        graph['Mean'].SetPoint(pcounter,i-binwidth+1,projs[-1].GetMean())
        graph['Mean'].SetPointError(pcounter,float(1/2.),projs[-1].GetMeanError())
        graph['RMS'].SetName(hist.GetName()+'_RMS')
        graph['Mean'].SetName(hist.GetName()+'_Mean')
        mini['RMS']=min(mini.get('RMS',1000),projs[-1].GetRMS())
        maxi['RMS']=max(maxi.get('RMS',0),projs[-1].GetRMS())
        mini['Mean']=min(mini.get('Mean',1000),projs[-1].GetMean())
        maxi['Mean']=max(maxi.get('Mean',0),projs[-1].GetMean())
        pcounter+=1
        i+=1
        #if (i>=nbins):
        #    break
        

    return projs

def setGraph(g,title,style,size=1,col=1,xtit='',ytit='',xtoff=1,ytoff=1):
    g.SetTitle(title)
    g.SetMarkerStyle(style)
    g.SetMarkerSize(size)
    g.SetMarkerColor(col)
    g.GetXaxis().SetTitle(xtit)
    g.GetYaxis().SetTitle(ytit)
    g.GetXaxis().SetTitleOffset(xtoff)
    g.GetXaxis().SetTitleOffset(ytoff)

    return

def getGraphYMaxMin(glist):
    px,py=ROOT.Double(0),ROOT.Double(0)
    mini=9999
    maxi=0
    for g in glist:
        for i in range(g.GetN()):
            g.GetPoint(i,px,py)
            mini=min(mini,float(py))
            maxi=max(maxi,float(py))

    if mini<0:
        mini*=1.1
    else:
        mini*=0.9
    if maxi<0:
        maxi*=0.9
    else:
        maxi*=1.1
    return maxi,mini

def thisIsTheEnd():
    u='w'
    while u!='q':
        u=raw_input('Done?')
    
def calcCovScale(cov,res,cut=True):
    cov2=ROOT.TMatrixDSym(3,cov.GetMatrixArray())
    cov_i=ROOT.TMatrixDSym(cov2)
    cov_i.InvertFast()

    residual=ROOT.TMatrixD(3,1)

    if res.Mag()==0:
        return -1,-1

    for i in range(3):
        residual[i][0]=res[i]/res.Mag()

    residual_t=ROOT.TMatrixD(ROOT.TMatrixD.kTransposed,residual)

    if cut:
        temp2=ROOT.TMatrixD(cov_i,ROOT.TMatrixD.kMult,residual)
    else:
        temp2=ROOT.TMatrixD(cov,ROOT.TMatrixD.kMult,residual)

    lam=ROOT.TMatrixD(residual_t,ROOT.TMatrixD.kMult,temp2)
    
    if cut:
        lam=1./math.sqrt(abs(lam[0][0]))
    else:
        lam=math.sqrt(abs(lam[0][0]))

    num=res.Mag()/abs(lam)

    return num,abs(lam)

def openTree(recofile,runs='',pattern=''):
    if recofile=='':
        print 'no recofile given'
        return False
    
    pyChain=PyChainMaker()     
    #pyChain.setDebug()
    if os.path.isdir(recofile):
        pyChain.setPath(recofile)
        if runs!='' and pattern!='':
            pyChain.setRuns(runs)
            pyChain.setPattern(pattern)
        else:
            pyChain.setMC()
    else:
        pyChain.setMC()
        pyChain.setFiles([recofile])
        
    tree=pyChain.getChain()
             
    if tree==None:
        print 'no tree found'
        return False
    
    return tree