import ROOT,sys,math,os,time,glob
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')
import argparse
from functions import *
from  cosmic_histos import *

parser=argparse.ArgumentParser(description='analyze cosmics from ana file')
#parser.add_argument("--branch",help="which branch to use",default="TrackFitStat_0")
parser.add_argument("filename",help="the file to analyze")
parser.add_argument("--slicef",help='a root file with the slices histo',default="")
parser.add_argument("--nEvts",help='number of events to analyze',default=-1,type=int)
parser.add_argument("--outfname",help="name of the output file",default="",type=str)
parser.add_argument('--runs', help='if data, the runs to use (%(default)s)', type=str, default='3888')
parser.add_argument('--pattern', help='if data, the file name patters ('')', type=str, default='')
parser.add_argument("--drawSingle",help="draw a single histo or dict",default="",type=str)
parser.add_argument("--covErr",help="use errors derived from shapematrix",action="store_true")
parser.add_argument("--axErr",help='use errors from cov matrix along axis',action="store_true")
parser.add_argument("--cuts",help='enable cosmics cuts (normally not used)',action='store_true')
parser.add_argument("--hlp",help="display help",action='store_true')
parser.add_argument('--batch',help="batchmode",action='store_true')
parser.add_argument('--rangefile',help='file with the ranges for the histograms',type=str,default='defaultRanges.txt')
parser.add_argument('--nopostprocess',help='do not post process histograms',action='store_true')
parser.add_argument('--parfile',help='the param file',type=str,default='')
parser.add_argument('--info',help='additional information. a string: name=val,name2=val2,...',type=str,default='')

args=parser.parse_args()

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

if args.batch:
    ROOT.gROOT.SetBatch(True)

print 'the pandapath=',pandapath

ROOT.gROOT.ProcessLine(".x rootlogon.C") 
ROOT.gROOT.ProcessLine('gStyle->SetPalette(1)')
ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")')

if args.outfname=="":
    args.outfname=args.filename.replace('.reco.root','.hist.root')

print 'histos will be written into:', args.outfname
histsdict=cosmic_histos(args.rangefile)
histsopt=cosmic_histopt()
canvasopt=cosmic_canvasopt()
axlabel=cosmic_axlabel()

zSlices = []
if args.slicef!="":
    print "slicename", args.slicef
    infiles=ROOT.TFile(args.slicef)
    zSlices = []
    histsdict['hslices'] = infiles.Get("hslices")
    for k in range(1,7):
        print histsdict['hslices'].GetBinContent(k)
        zSlices.append(histsdict['hslices']).GetBinContent(k)
else:
    zSlices = (  12.0, 24.0, 36.0,  48.0, 60.0, 73.0 ) # cosmic ne
    #zSlices =  (   20,  30,  40,  49,  57, 70 ) # beam
    for i in range(len(zSlices)):
        histsdict['hslices'].SetBinContent(i+1,zSlices[i])

phiSlices = (30,60,90,120,150,180)
thetaSlices=[i*20 for i in range(9)]
hslices=cosmic_hslices(zSlices,args.rangefile)


#RFile=ROOT.TFile()
#Rfile = ROOT.TFile.Open(args.filename,"read")
#tree = Rfile.Get("cbmsim")
tree=openTree(args.filename,args.runs,args.pattern)
if tree == None :
    print "No cbmsim in file", file
    if type(Rfile)==ROOT.TFile:
        Rfile.Close()
        exit()

haspar=False
if args.parfile!='':
    parfilename=args.parfile
else:
    parfilename=args.filename.replace('reco','param')

if os.path.isfile(parfilename):
    parfile=ROOT.TFile.Open(parfilename)
    digiPar=parfile.Get('TpcDigiPar')
    haspar=True
if not haspar:
    paths=parfilename.split('/')
    path_idx=[i for i,p in enumerate(paths) if 'parts' in p]
    for idx in path_idx:
        paths.pop(idx)
    paths.pop(-1)
    paths.pop(-1)
    parfilepath="/".join(paths)
    parfilename=glob.glob(parfilepath+'/*.param.root')
    if len(parfilename)!=0:
        parfilename=parfilename[0]
    else:
        parfilename=''
    print parfilename
    if os.path.isfile(parfilename):
        parfile=ROOT.TFile.Open(parfilename)
        digiPar=parfile.Get('TpcDigiPar')
        haspar=True

#print 'gas test:',digiPar.getGas().Dt()

tree.SetBranchStatus("*", 0)
#tree.SetBranchStatus(args.branch+".*",1)
tree.SetBranchStatus("CutCosmics.*",1)

trackbranch=ROOT.TBranch()
evt=-1
step_ev=-1
done=0
olddone=0
if args.nEvts==-1 or tree.GetEntries()<args.nEvts:
    todoev=tree.GetEntries()
else:
    todoev=args.nEvts
print "starting:", args.filename, "("+str(todoev)+")"

main_start=time.time()
ev_start=time.time()
digiPerpVariance=[]
digiNum=[]
for iz in range(75):
    digiPerpVariance.append(0)
    digiNum.append(0)

for e in tree:
    step_ev+=1
    evt+=1
    if args.nEvts!=-1 and args.nEvts<=evt:
        break
    trk=e.CutCosmics.GetEntries()

    if trk != 1:
        continue

    trackbranch=e.CutCosmics
            
    for tr in trackbranch:
        if args.cuts and tr.GetStatFlag()!=0:
            continue

        if args.cuts and (tr.GetTheta()<40 or tr.GetTheta()>140):
            continue
        
        if args.cuts and (tr.GetClMeanAmp()<400):
            continue

        if args.cuts and (tr.GetClMeanAmp()>5000 and tr.GetLength()<2):
            continue


        phi=tr.GetPhi()
        theta=tr.GetTheta()
        clMeanAmp=tr.GetClMeanAmp()
        histsdict['theta']['no cuts'].Fill(theta)
        histsdict['phi']['no cuts'].Fill(phi)
        histsdict['theta vs clmeanamp'].Fill(clMeanAmp,theta)
        histsdict['clmeanamp'].Fill(clMeanAmp)

        resX=tr.GetResX()
        resY=tr.GetResY()
        resZ=tr.GetResZ()
        resU=tr.GetResU()
        resV=tr.GetResV()

        resA=tr.GetResA()
        resAMC=tr.GetResAMC()
        resAMC_bad=False
        if len(resAMC)==0:
            resAMC_bad=True
            resAMC=[]
        
        axis1=tr.GetAxis1()
        axis2=tr.GetAxis2()
        axis3=tr.GetAxis3()

        resR=tr.GetResR()
        resP=tr.GetResP()

        posX=tr.GetHitPositionsX()
        posY=tr.GetHitPositionsY()
        posZ=tr.GetHitPositionsZ()

        errA=tr.GetSigA()
        errC=tr.GetSigC()
        errX=tr.GetSigX()
        errY=tr.GetSigY()
        errZ=tr.GetSigZ()
        errU=tr.GetSigU()
        errV=tr.GetSigV()


        refPos=tr.GetProjectionPoints()
        refPosErr=tr.GetProjectionPointsErr()

        trackPos=tr.GetPos()
        trackPosErr=tr.GetPosErr()
        trackMCpos=tr.GetMCpos()
        trackMom=tr.GetMom()
        trackMomErr=tr.GetMomErr()
        trackMCmom=tr.GetMCmom()

        size2D=tr.Get2DClSizes()
        size3D=tr.GetClSizes()
        size3DU=tr.GetClSizesU()
        
        length=tr.GetLength()
        distance=tr.GetClusterDists()
        
        digiSNR=tr.GetDigiSNR()
        clusterSNR=tr.GetClusterSNR()
        sampleSNR=tr.GetSampleSNR()
        
        chi2=tr.getChi2()
        NDF=tr.getNDF() 
        meanz=tr.GetMeanZ()

        Chi2X=tr.GetChi2X()

        resMcX=tr.GetResMcX()  
        resMcY=tr.GetResMcY()
        resMcZ=tr.GetResMcZ()
        resMcR=tr.GetResMcR()
        resMcP=tr.GetResMcP()
        resMcU=tr.GetResMCU()
        resMcV=tr.GetResMCV()

        refPosXMC=tr.GetMcRefPosX()
        refPosYMC=tr.GetMcRefPosY()
        refPosZMC=tr.GetMcRefPosZ()

        resMCTrack=tr.GetResMCTrack()
        resMCTrackUV=tr.GetResMCTrackUV()
        resMCTrackA=tr.GetResMCTrackA()
        resMCTrackA_bad=False
        if len(resMCTrackA)==0:
            resMCTrackA=[]
            resMCTrackA_bad=True
        digis_per_track=0
        resMcA0=[]
        resMcA1=[]
        resMcA2=[]
        
        #verena code, TODO clean this up
        for i in range(len(size3D)):
            '''
            histsdict['hdbc'].Fill(posZ[i],distance[i])
            histsdict['hcs2D'].Fill(posZ[i],size2D[i])
            histsdict['hcs3D'].Fill(posZ[i],size3D[i])
            diff = abs(size3D[i]-size2D[i])
            histsdict['hdsoc'].Fill(posZ[i],diff)
            if posZ[i]<-48.0:
                histsdict['hdpc1'].Fill(size3D[i])
            if -48<posZ[i]<-36.0:
                histsdict['hdpc2'].Fill(size3D[i])
            if -36.0<posZ[i]<-24.0:
                histsdict['hdpc3'].Fill(size3D[i]) 
            if -24.0<posZ[i]<-12.0:
                histsdict['hdpc4'].Fill(size3D[i])
            if -12.0<posZ[i]<0.0:
                histsdict['hdpc5'].Fill(size3D[i])
            if 0<posZ[i]:
                histsdict['hdpc6'].Fill(size3D[i])
            '''
            for i in range(len(resMcX)):
                if abs(phi)>70 and abs(phi)<110:
                    histsdict['hresZ'].Fill(posZ[i],resX[i])


        for icl in range(len(resX)):
 
            resX[icl]*=  10000
            resY[icl]*=  10000
            resZ[icl]*=  10000

            resMcX[icl]*=10000
            resMcY[icl]*=10000
            resMcZ[icl]*=10000
            
            resA[icl]*=  10000

            resU[icl]*=  10000
            resV[icl]*=  10000
            resMcU[icl]*=10000
            resMcV[icl]*=10000
            
            errX[icl]*=  10000
            errY[icl]*=  10000
            errZ[icl]*=  10000
            errA[icl]*=  10000
            errC[icl]*=  10000
            errU[icl]*=  10000
            errV[icl]*=  10000

            resP[icl]*=  10000
            resR[icl]*=  10000
            
            resMcR[icl]*=10000
            resMcP[icl]*=10000


            resMCTrack[icl]*=10000
            resMCTrackUV[icl]*=10000
            
            
            refPosErr[icl]*=10000
            
            #this if block is only due to a bug in cosmics task where resAMC was not set properly and resA was set to resAMC
            if resAMC_bad:
                #project mc position onto axis
                resMC=ROOT.TVector3(resMcX[icl],resMcY[icl],resMcZ[icl])
                res=ROOT.TVector3(resX[icl],resY[icl],resZ[icl])
                if axis1[icl].Mag()!=0:
                    resMcA0.append((resMC*axis1[icl])/axis1[icl].Mag())
                    resA[icl].SetX(res*axis1[icl]/axis1[icl].Mag())
                else:
                    resMcA0.append(-99999)
                    resA[icl].SetX(-99999)

                if axis2[icl].Mag()!=0:
                    resMcA1.append((resMC*axis2[icl])/axis2[icl].Mag())
                    resA[icl].SetY(res*axis2[icl]/axis2[icl].Mag())
                else:
                    resMcA1.append(-99999)
                    resA[icl].SetY(-99999)

                if axis3[icl].Mag()!=0:
                    resMcA2.append((resMC*axis3[icl])/axis3[icl].Mag())
                    resA[icl].SetZ(res*axis3[icl]/axis3[icl].Mag())
                else:
                    resMcA2.append(-99999)
                    resA[icl].SetZ(-99999)
                resAMC.append(ROOT.TVector3(resMcA0[-1],resMcA1[-1],resMcA2[-1]))
                
            else:
                resAMC[icl]*=10000
                resMcA0.append(resAMC[icl].X())
                resMcA1.append(resAMC[icl].Y())
                resMcA2.append(resAMC[icl].Z())

            #this block is for older files where resMCTrackA was not calculated
            if resMCTrackA_bad:
                #project residual onto axis
                res=ROOT.TVector3(resX[icl],resY[icl],resZ[icl])
                resmctrackA0=0
                resmctrackA1=0
                resmctrackA2=0
                if axis1[icl].Mag()!=0:
                    resmctrackA0=(resMCTrack[icl]*axis1[icl])/axis1[icl].Mag() 
                if axis2[icl].Mag()!=0:   
                    resmctrackA1=(resMCTrack[icl]*axis2[icl])/axis2[icl].Mag()
                if axis3[icl].Mag()!=0:   
                    resmctrackA2=(resMCTrack[icl]*axis3[icl])/axis3[icl].Mag() 
                resMCTrackA.append(ROOT.TVector3(resmctrackA0,resmctrackA1,resmctrackA2))
            else:
                resMCTrackA[icl]*=10000

            rad=math.sqrt(posX[icl]**2+posY[icl]**2)
            refRad=math.sqrt(refPos[icl].X()**2+refPos[icl].Y()**2)
            refRadMC=math.sqrt(refPosXMC[icl]**2+refPosYMC[icl]**2)

            digis_per_track+=size3DU[icl]
            refPosMC=ROOT.TVector3(refPosXMC[icl],refPosYMC[icl],refPosZMC[icl])

            #2D residual histograms (nt as function of Z)
            if 5. < rad and rad < 15.:
                histsdict['hxy'].Fill(posX[icl],posY[icl])
                histsdict['hxyz'].Fill(posX[icl],posY[icl],posZ[icl])
                histsdict['hxz'].Fill(posZ[icl],posX[icl])
                histsdict['hyz'].Fill(posZ[icl],posY[icl])

                histsdict['hxy MC'].Fill(refPosXMC[icl],refPosYMC[icl])
                histsdict['hxz MC'].Fill(refPosZMC[icl],refPosXMC[icl])
                histsdict['hyz MC'].Fill(refPosZMC[icl],refPosYMC[icl])
                
                #X-Residuals
                histsdict['hxy res'].Fill(posX[icl],posY[icl],resX[icl])
                histsdict['hxy MC res'].Fill(refPosXMC[icl],refPosYMC[icl],resMcX[icl])
                
                histsdict['hyz res'].Fill(posZ[icl],posY[icl],resX[icl])
                histsdict['hyz MC res'].Fill(refPosZMC[icl],refPosYMC[icl],resMcX[icl])

                #Y-Residuals
                histsdict['hxy resY'].Fill(posX[icl],posY[icl],resY[icl])
                histsdict['hxy MC resY'].Fill(refPosXMC[icl],refPosYMC[icl],resMcY[icl])
                
                histsdict['hyz resY'].Fill(posZ[icl],posY[icl],resY[icl])
                histsdict['hyz MC resY'].Fill(refPosZMC[icl],refPosYMC[icl],resMcY[icl])

                

                #Z-Residuals
                histsdict['hxy resZ'].Fill(posX[icl],posY[icl],resZ[icl])
                histsdict['hxy MC resZ'].Fill(refPosXMC[icl],refPosYMC[icl],resMcZ[icl])
                
                histsdict['hyz resZ'].Fill(posZ[icl],posY[icl],resZ[icl])
                histsdict['hyz MC resZ'].Fill(refPosZMC[icl],refPosYMC[icl],resMcZ[icl])

                
                
                histsdict['hxyz res'].Fill(posX[icl],posY[icl],posZ[icl],resX[icl])
                histsdict['hxyz resY'].Fill(posX[icl],posY[icl],posZ[icl],resY[icl])
                histsdict['hxyz resZ'].Fill(posX[icl],posY[icl],posZ[icl],resZ[icl])
                histsdict['hxyz resU']    .Fill(posX[icl],posY[icl],posZ[icl],resU[icl])    
                histsdict['hxyz resV']    .Fill(posX[icl],posY[icl],posZ[icl],resV[icl])
                histsdict['hxyz resA1']   .Fill(posX[icl],posY[icl],posZ[icl],resA[icl].X())
                histsdict['hxyz resA2']   .Fill(posX[icl],posY[icl],posZ[icl],resA[icl].Y())
                histsdict['hxyz resA3']   .Fill(posX[icl],posY[icl],posZ[icl],resA[icl].Z())
                                          
                histsdict['hxyz MC res'].Fill(posX[icl],posY[icl],posZ[icl] ,resMcX[icl])
                histsdict['hxyz MC resY'].Fill(posX[icl],posY[icl],posZ[icl],resMcY[icl])
                histsdict['hxyz MC resZ'].Fill(posX[icl],posY[icl],posZ[icl],resMcZ[icl])                                          
                histsdict['hxyz MC resU'] .Fill(posX[icl],posY[icl],posZ[icl],resMcU[icl])
                histsdict['hxyz MC resV'] .Fill(posX[icl],posY[icl],posZ[icl],resMcV[icl])
                histsdict['hxyz MC resA1'].Fill(posX[icl],posY[icl],posZ[icl],resAMC[icl].X())
                histsdict['hxyz MC resA2'].Fill(posX[icl],posY[icl],posZ[icl],resAMC[icl].Y())
                histsdict['hxyz MC resA3'].Fill(posX[icl],posY[icl],posZ[icl],resAMC[icl].Z())


                histsdict['hxy chi2X'].Fill(posX[icl],posY[icl],Chi2X[icl])
               
                histsdict['hxy resP'].Fill(refPos[icl].X(),refPos[icl].Y(),resP[icl])
                histsdict['hxy MC resP'].Fill(refPosXMC[icl],refPosYMC[icl],resMcP[icl])
                histsdict['hxy resR'].Fill(refPos[icl].X(),refPos[icl].Y(),resR[icl])
                histsdict['hxy MC resR'].Fill(refPosXMC[icl],refPosYMC[icl],resMcR[icl]) 
                histsdict['hxy resU'].Fill(refPos[icl].X(),refPos[icl].Y(),resU[icl])
                histsdict['hxy MC resU'].Fill(refPosXMC[icl],refPosYMC[icl],resMcU[icl])
                histsdict['hxy resV'].Fill(refPos[icl].X(),refPos[icl].Y(),resV[icl])
                histsdict['hxy MC resV'].Fill(refPosXMC[icl],refPosYMC[icl],resMcV[icl]) 


            
                for j in range(len(zSlices)):
                    if refPos[icl].Z()<zSlices[j]:
                        hslices['hxy slice'][j].Fill(posX[icl],posY[icl])
                        hslices['hxy res slice'][j].Fill(posX[icl],posY[icl],resX[icl])
                        hslices['residualsX'][j].Fill(resX[icl])
                        break

                for j in range(len(zSlices)):
                    if refPosZMC[icl]<zSlices[j]:
                        hslices['hxy MC res slice'][j].Fill(refPosXMC[icl],refPosYMC[icl],resMcX[icl])
                        hslices['MCresidualsX'][j].Fill(resMcX[icl])
                        break

                for j in range(len(zSlices)):
                    if refPos[icl].Z()<zSlices[j]:
                        hslices['hxy Y slice'][j].Fill(posX[icl],posY[icl])
                        hslices['hxy resY slice'][j].Fill(posX[icl],posY[icl],resY[icl])
                        hslices['residualsY'][j].Fill(resY[icl])
                        break

                for j in range(len(zSlices)):
                    if refPosZMC[icl]<zSlices[j]:
                        hslices['hxy MC resY slice'][j].Fill(refPosXMC[icl],refPosYMC[icl],resMcY[icl])
                        hslices['MCresidualsY'][j].Fill(resMcY[icl])
                        break

                for j in range(len(zSlices)):
                    if refPos[icl].Z()<zSlices[j]:
                        hslices['hxy Z slice'][j].Fill(posX[icl],posY[icl])
                        hslices['hxy resZ slice'][j].Fill(posX[icl],posY[icl],resZ[icl])
                        hslices['residualsZ'][j].Fill(resZ[icl])
                        break

                for j in range(len(zSlices)):
                    if refPosZMC[icl]<zSlices[j]:
                        hslices['hxy MC resZ slice'][j].Fill(refPosXMC[icl],refPosYMC[icl],resMcZ[icl])
                        hslices['MCresidualsZ'][j].Fill(resMcZ[icl])
                        break
                
                histsdict['res vs theta']['cluster MC X'].Fill(theta,resMcX[icl])                          
                histsdict['res vs Z']['cluster track X'].Fill(refPos[icl].Z(),resX[icl])
                histsdict['res vs Z']['cluster MC X'].Fill(refPosZMC[icl],resMcX[icl])
                histsdict['res vs Z']['track MC X'].Fill(refPosZMC[icl],resMCTrack[icl].X())
                
                histsdict['res vs theta']['cluster MC Y'].Fill(theta,resMcY[icl])
                histsdict['res vs Z']['cluster MC Y'].Fill(refPosZMC[icl],resMcY[icl])
                histsdict['res vs Z']['cluster track Y'].Fill(refPos[icl].Z(),resY[icl])
                histsdict['res vs Z']['track MC Y'].Fill(refPosZMC[icl],resMCTrack[icl].Y())

                histsdict['res vs theta']['cluster MC Z'].Fill(theta,resMcZ[icl])
                histsdict['res vs Z']['cluster MC Z'].Fill(refPosZMC[icl],resMcZ[icl])
                histsdict['res vs Z']['cluster track Z'].Fill(refPos[icl].Z(),resZ[icl])
                histsdict['res vs Z']['track MC Z'].Fill(refPosZMC[icl],resMCTrack[icl].Z())

                histsdict['res vs Z']['cluster track U'].Fill(refPos[icl].Z(),resU[icl])
                histsdict['res vs Z']['cluster track V'].Fill(refPos[icl].Z(),resV[icl])
                histsdict['res vs Z']['cluster MC U'].Fill(refPosZMC[icl],resMcU[icl])
                histsdict['res vs Z']['cluster MC V'].Fill(refPosZMC[icl],resMcV[icl])
                histsdict['res vs Z']['track MC U'].Fill(refPosZMC[icl],resMCTrackUV[icl].X())
                histsdict['res vs Z']['track MC V'].Fill(refPosZMC[icl],resMCTrackUV[icl].Y())

                histsdict['res vs Z']['cluster track A0'].Fill(refPos[icl].Z(),resA[icl].X())
                histsdict['res vs Z']['cluster track A1'].Fill(refPos[icl].Z(),resA[icl].Y())
                histsdict['res vs Z']['cluster track A2'].Fill(refPos[icl].Z(),resA[icl].Z())
                histsdict['res vs Z']['cluster MC A0'].Fill(refPosZMC[icl],resAMC[icl].X())
                histsdict['res vs Z']['cluster MC A1'].Fill(refPosZMC[icl],resAMC[icl].Y())
                histsdict['res vs Z']['cluster MC A2'].Fill(refPosZMC[icl],resAMC[icl].Z())
                histsdict['res vs Z']['track MC A0'].Fill(refPosZMC[icl],resMCTrackA[icl].X())
                histsdict['res vs Z']['track MC A1'].Fill(refPosZMC[icl],resMCTrackA[icl].Y())
                histsdict['res vs Z']['track MC A2'].Fill(refPosZMC[icl],resMCTrackA[icl].Z())

                if(refPos[icl].Perp()>=7 and refPos[icl].Perp()<=13):
                    if (abs(abs(phi)-90)<20):
                        histsdict['res vs Z']['cluster track X Cut20'].Fill(refPos[icl].Z(),resX[icl])
                    if (abs(phi)<20):
                        histsdict['res vs Z']['cluster track Y Cut20'].Fill(refPos[icl].Z(),resY[icl])
                    if (abs(abs(theta)-90)<20):
                        histsdict['res vs Z']['cluster track Z Cut20'].Fill(refPos[icl].Z(),resZ[icl])

                    if (abs(abs(phi)-90)<10):
                        histsdict['res vs Z']['cluster track X Cut10'].Fill(refPos[icl].Z(),resX[icl])
                    if (abs(phi)<10):
                        histsdict['res vs Z']['cluster track Y Cut10'].Fill(refPos[icl].Z(),resY[icl])
                    if (abs(abs(theta)-90)<10):
                        histsdict['res vs Z']['cluster track Z Cut10'].Fill(refPos[icl].Z(),resZ[icl])
                        
                    histsdict['res vs Z']['cluster track U Cut'].Fill(refPos[icl].Z(),resU[icl])
                    histsdict['res vs Z']['cluster track V Cut'].Fill(refPos[icl].Z(),resV[icl])
                
                fitgood=tr.GetClusterFitGood(icl)
                isGeomErr=tr.GetIsGeomErr(icl)

                if (fitgood and isGeomErr.Mag2()==1):
                    histsdict['res vs Z']['cluster track A0 fg'].Fill(refPos[icl].Z(),resA[icl].X())
                    histsdict['res vs Z']['cluster track A1 fg'].Fill(refPos[icl].Z(),resA[icl].Y())
                    histsdict['res vs Z']['cluster track A2 fg'].Fill(refPos[icl].Z(),resA[icl].Z())
                    histsdict['res vs Z']['cluster MC A0 fg'].Fill(refPosZMC[icl],resAMC[icl].X())
                    histsdict['res vs Z']['cluster MC A1 fg'].Fill(refPosZMC[icl],resAMC[icl].Y())
                    histsdict['res vs Z']['cluster MC A2 fg'].Fill(refPosZMC[icl],resAMC[icl].Z())


                histsdict['clusterdists'].Fill(distance[icl])
                histsdict['clusterdists vs z'].Fill(refPos[icl].Z(),distance[icl])
                histsdict['hclusterSize2DVsZ'].Fill(refPos[icl].Z(),size2D[icl])
                histsdict['hclusterSize3DVsZ'].Fill(refPos[icl].Z(),size3D[icl])
                histsdict['hclusterSize3DUVsZ'].Fill(refPos[icl].Z(),size3DU[icl])
                histsdict['res vs theta']['cluster MC rp'].Fill(theta,resMcP[icl]*refRadMC)
                histsdict['res vs Z']['cluster MC rp'].Fill(refPosZMC[icl],resMcP[icl]*refRadMC)
                
                histsdict['hsnr'].Fill(digiSNR[icl])
                histsdict['csnr'].Fill(clusterSNR[icl])
                histsdict['ssnr'].Fill(sampleSNR[icl])

                if errU[icl]!=0:
                    histsdict['pulls']['cluster MCU vs z'].Fill(refPos[icl].Z(),resMcU[icl]/errU[icl])
                    histsdict['pulls']['cluster U vs z'].Fill(refPos[icl].Z(),resU[icl]/errU[icl])
                if errV[icl]!=0:
                    histsdict['pulls']['cluster MCV vs z'].Fill(refPos[icl].Z(),resMcV[icl]/errV[icl])
                    histsdict['pulls']['cluster V vs z'].Fill(refPos[icl].Z(),resV[icl]/errV[icl])
                
                if errX[icl]!=0:
                    histsdict['pulls']['cluster x vs z'].Fill(refPos[icl].Z(),resX[icl]/errX[icl])
                    histsdict['pulls']['cluster MCx vs z'].Fill(refPosZMC[icl],resMcX[icl]/errX[icl])
                    histsdict['error']['cluster err X vs z'].Fill(refPos[icl].Z(),errX[icl])
                if errY[icl]!=0:
                    histsdict['pulls']['cluster y vs z'].Fill(refPos[icl].Z(),resY[icl]/errY[icl])
                    histsdict['pulls']['cluster MCy vs z'].Fill(refPosZMC[icl],resMcY[icl]/errY[icl])
                    histsdict['error']['cluster err Y vs z'].Fill(refPos[icl].Z(),errY[icl])
                if errZ[icl]!=0:
                    histsdict['pulls']['cluster z vs z'].Fill(refPos[icl].Z(),resZ[icl]/errZ[icl])
                    histsdict['pulls']['cluster MCz vs z'].Fill(refPosZMC[icl],resMcZ[icl]/errZ[icl])
                    histsdict['error']['cluster err Z vs z'].Fill(refPos[icl].Z(),errZ[icl])

                if errA[icl].X()!=0:
                    histsdict['pulls']['cluster a0 vs z'].Fill(refPos[icl].Z(),resA[icl].X()/errA[icl].X())
                    histsdict['pulls']['cluster MCa0 vs z'].Fill(refPosZMC[icl],resMcA0[icl]/errA[icl].X())
                    histsdict['error']['cluster err A0 vs z'].Fill(refPos[icl].Z(),errA[icl].X())
                    if (fitgood and isGeomErr.Mag2()==1):
                        histsdict['pulls']['cluster MCa0 vs z fg'].Fill(refPosZMC[icl],resMcA0[icl]/errA[icl].X())
                if errA[icl].Y()!=0:
                    histsdict['pulls']['cluster a1 vs z'].Fill(refPos[icl].Z(),resA[icl].Y()/errA[icl].Y())
                    histsdict['pulls']['cluster MCa1 vs z'].Fill(refPosZMC[icl],resMcA1[icl]/errA[icl].Y())
                    histsdict['error']['cluster err A1 vs z'].Fill(refPos[icl].Z(),errA[icl].Y())
                    if (fitgood and isGeomErr.Mag2()==1):
                        histsdict['pulls']['cluster MCa1 vs z fg'].Fill(refPosZMC[icl],resMcA1[icl]/errA[icl].Y())
                if errA[icl].Z()!=0:
                    histsdict['pulls']['cluster a2 vs z'].Fill(refPos[icl].Z(),resA[icl].Z()/errA[icl].Z())
                    histsdict['pulls']['cluster MCa2 vs z'].Fill(refPosZMC[icl],resMcA2[icl]/errA[icl].Z())
                    histsdict['error']['cluster err A2 vs z'].Fill(refPos[icl].Z(),errA[icl].Z())
                    if (fitgood and isGeomErr.Mag2()==1):
                        histsdict['pulls']['cluster MCa2 vs z fg'].Fill(refPosZMC[icl],resMcA2[icl]/errA[icl].Z())

                

                trackRes=resMCTrack[icl]
                if refPosErr[icl].X()!=0:
                    histsdict['pulls']['clustert x vs z'].Fill(refPos[icl].Z(),resX[icl]/(refPosErr[icl].X()))
                    histsdict['pulls']['trackp x vs z'].Fill(refPos[icl].Z(),trackRes.X()/refPosErr[icl].X())

                if refPosErr[icl].Y()!=0:
                    histsdict['pulls']['clustert y vs z'].Fill(refPos[icl].Z(),resY[icl]/(refPosErr[icl].Y()))
                    histsdict['pulls']['trackp y vs z'].Fill(refPos[icl].Z(),trackRes.Y()/refPosErr[icl].Y())
                
                if refPosErr[icl].Z()!=0:
                    histsdict['pulls']['clustert z vs z'].Fill(refPos[icl].Z(),resZ[icl]/(refPosErr[icl].Z()))
                    histsdict['pulls']['trackp z vs z'].Fill(refPos[icl].Z(),trackRes.Z()/refPosErr[icl].Z())

                #fill cov cut variables
                cov=tr.GetClusterCov(icl)
                cov*=10000**2
                covPullMC,covErrMC=calcCovScale(cov,resAMC[icl])
                histsdict['pulls']['cov cut MC'].Fill(refPos[icl].Z(),covPullMC)
                histsdict['error']['cov res cut MC'].Fill(refPos[icl].Z(),covErrMC)
                histsdict['res vs Z']['residualMagMC'].Fill(refPos[icl].Z(),resAMC[icl].Mag())
                covPull,covErr=calcCovScale(cov,resA[icl])
                histsdict['pulls']['cov cut'].Fill(refPos[icl].Z(),covPull)
                histsdict['error']['cov res cut'].Fill(refPos[icl].Z(),covErr)
                histsdict['res vs Z']['residualMag'].Fill(refPos[icl].Z(),resA[icl].Mag())
                    
                histsdict['error']['cluster err RP vs z'].Fill(refPos[icl].Z(),9999)
                histsdict['error']['track err X vs z'].Fill(refPos[icl].Z(),refPosErr[icl].X())
                histsdict['error']['track err Y vs z'].Fill(refPos[icl].Z(),refPosErr[icl].Y())
                histsdict['error']['track err Z vs z'].Fill(refPos[icl].Z(),refPosErr[icl].Z())

                histsdict['hclusteramp'].Fill(tr.GetAmp(icl))

                for j in range(len(zSlices)):
                    if refPos[icl].Z()<zSlices[j]:
                        hslices['clustersize'][j].Fill(size3D[icl])
                        hslices['clustersizeU'][j].Fill(size3DU[icl])
                        hslices['clustersize2D'][j].Fill(size2D[icl])
                        hslices['nsnr'][j].Fill(digiSNR[icl])
                        hslices['csnr'][j].Fill(clusterSNR[icl])
                        hslices['hxy resR slice'][j].Fill(refPos[icl].X(),refPos[icl].Y(),resR[icl])
                        hslices['hxy resP slice'][j].Fill(refPos[icl].X(),refPos[icl].Y(),resP[icl])
                        hslices['hxy resU slice'][j].Fill(refPos[icl].X(),refPos[icl].Y(),resU[icl])
                        hslices['hxy resV slice'][j].Fill(refPos[icl].X(),refPos[icl].Y(),resV[icl])
                        hslices['hxy all slice'][j].Fill(refPos[icl].X(),refPos[icl].Y())
                        break
                for j in range(len(zSlices)):
                    if refPosZMC[icl]<zSlices[j]:
                        hslices['hxy MC resR slice'][j].Fill(refPosXMC[icl],refPosYMC[icl],resMcR[icl])
                        hslices['hxy MC resP slice'][j].Fill(refPosXMC[icl],refPosYMC[icl],resMcP[icl])
                        break
                for j in range(len(phiSlices)):
                    if abs(phi)<phiSlices[j]:
                        hslices['hxy slice phi'][j].Fill(refPos[icl].X(),refPos[icl].Y())
                        hslices['hxy resX phi slice'][j].Fill(refPos[icl].X(),refPos[icl].Y(),resX[icl])
                        hslices['hxy MC resX phi slice'][j].Fill(refPosXMC[icl],refPosYMC[icl],resMcX[icl])
                        break
                
            
                #start filling DIGI based histograms
                digiTrackRes=tr.GetDigiTrackRes(icl)
                digiTrackMom=tr.GetDigiMom(icl)
                digiPos=tr.GetDigiPos(icl)
                digilen=tr.GetDigiLen(icl)
                  
                for idi in range(len(digiTrackRes)):
                    sign=1
                    if (digiTrackMom[idi].DeltaPhi(digiTrackRes[idi])<0):
                        sign=-1
                    histsdict['digi residuals']['digi-track perp vs z'].Fill(digiPos[idi].Z(),digiTrackRes[idi].Perp())
                    histsdict['digi residuals']['digi-track perp sign vs z'].Fill(digiPos[idi].Z(),digiTrackRes[idi].Perp()*sign)
                    histsdict['digi residuals']['digi-track x vs z'].Fill(digiPos[idi].Z(),digiTrackRes[idi].X())
                    histsdict['digi residuals']['digi-track y vs z'].Fill(digiPos[idi].Z(),digiTrackRes[idi].Y())
                    histsdict['digi residuals']['digi-track z vs z'].Fill(digiPos[idi].Z(),digiTrackRes[idi].Z())
                    
                    histsdict['digi residuals']['digilenVsZ'].Fill(digiPos[idi].Z(),digilen[idi])
                    if abs(theta-90)<20:
                        histsdict['digi residuals']['digilenVsZres'].Fill(digilen[idi],digiTrackRes[idi].Z())
                        
                    if int(digiPos[idi].Z())>=0 and int(digiPos[idi].Z())<75:
                        digiPerpVariance[int(digiPos[idi].Z())]+=digiTrackRes[idi].Perp()
                        digiNum[int(digiPos[idi].Z())]+=1
                    histsdict['hdigiamp'].Fill(abs(tr.GetDigisOnPlaneAmp(icl,idi)))
                    if(abs(tr.GetDigisOnPlaneAmp(icl,idi)))<500:
                        histsdict['hdigiamp_low'].Fill(abs(tr.GetDigisOnPlaneAmp(icl,idi)))
                        
                #done digi histos
        #closing the cluster loop -> back inside the loop over the tracks
        if tr.GetLength()>0:
            histsdict['hclusterpercm vs Z'].Fill(tr.GetMeanZ(),len(resX)/tr.GetLength())
            histsdict['hdigipercm vs z'].Fill(tr.GetMeanZ(),digis_per_track/tr.GetLength())
            histsdict['hclusterpercm'].Fill(len(resX)/tr.GetLength())
            histsdict['hdigipercm'].Fill(digis_per_track/tr.GetLength())
        histsdict['hclusterpertrack'].Fill(len(resX))
        histsdict['hdigipertrack'].Fill(digis_per_track)
        chi2cdf=ROOT.Math.chisquared_cdf(tr.getChi2(),tr.getNDF())
        histsdict['chi2cdf'].Fill(chi2cdf)

        phierr=math.sqrt( ( (trackMom.X() / (trackMom.X()**2+trackMom.Y()**2) ) *trackMomErr.Y())**2 +
                          ( (trackMom.Y() / (trackMom.X()**2+trackMom.Y()**2) ) *trackMomErr.X())**2)
        MCphi=trackMCmom.Phi()*ROOT.TMath.RadToDeg()
        histsdict['pulls']['track phi vs z'].Fill(tr.GetMeanZ(),(phi-MCphi)/phierr)
        histsdict['pulls']['phierr'].Fill(phierr)
        histsdict['pulls']['MCphi'].Fill(MCphi)
        histsdict['phi']['phi'].Fill(phi)
        Track_long_mom=trackMom.Z()
        Track_trans_mom=trackMom.Perp()
        Track_long_MCmom=trackMCmom.Z()
        Track_trans_MCmom=trackMCmom.Perp()

        Track_long_momErr=trackMomErr.Z()
        Track_trans_momErr=math.sqrt( ((trackMom.X()/Track_trans_mom)*trackMomErr.X())**2 +
                                      ((trackMom.Y()/Track_trans_mom)*trackMomErr.Y())**2 )

        pull_mom_long =Track_long_mom-Track_long_MCmom
        pull_mom_long/=Track_long_momErr
        pull_mom_trans =Track_trans_mom-Track_trans_MCmom
        pull_mom_trans/=Track_trans_momErr

        histsdict['pulls']['track moml vs z'].Fill(tr.GetMeanZ(), pull_mom_long)
        histsdict['pulls']['track momt vs z'].Fill(tr.GetMeanZ(), pull_mom_trans)
        for j in range(len(zSlices)):
            if tr.GetMeanZ()<zSlices[j] and tr.GetLength()>0:
                hslices['hcluster per cm'][j].Fill(len(resX)/tr.GetLength())
                hslices['hdigi per cm'][j].Fill(digis_per_track/tr.GetLength())
                hslices['htracklength'][j].Fill(tr.GetLength())
                hslices['hcluster per track'][j].Fill(len(resX))
                hslices['hdigi per track'][j].Fill(digis_per_track)
                hslices['hchi2ndf'][j].Fill(chi2cdf)
                break
        
    #closing the track loop -> back inside the loop over the events
    done=int(float(evt*100./todoev))
    if done!=olddone:
        if step_ev==0:
            continue
        t_per_ev=(time.time()-ev_start)/step_ev
        eta_tot=(todoev-evt)*t_per_ev
        eta_h=int(eta_tot/(60*60))
        eta_min=int((eta_tot-60*60*eta_h)/60)
        eta_s=eta_tot-60*60*eta_h-60*eta_min
        print "{0}% done, Good event: {1} (t_tot={2:.2}min; t_ev={3:.2}s; eta={4:02}:{5:02}:{6:02})".format(done,evt,(time.time()-main_start)/60,t_per_ev,eta_h,eta_min,int(eta_s))
        ev_start=time.time()
        step_ev=-1
        olddone=done

if step_ev!=0:
    print "{0}% done, Good event: {1} (t_tot={2:.2} min; t_ev={3:.2} s)".format(done,evt,(time.time()-main_start)/60,(time.time()-ev_start)/step_ev)

print 'all events processed'
#closing the event loop -> back in no loop at all

#X-Res-Norm
histsdict['hxy res norm'] =Divideh(histsdict['hxy res'] ,histsdict['hxy']  ,"hxyresnorm" ,"xy X-Residual normalized")
histsdict['hxy MCres norm'] =Divideh(histsdict['hxy MC res'] ,histsdict['hxy MC']  ,"hxyMCresnorm" ,"xy MC X-Residual normalized")

histsdict['hyz res norm'] =Divideh(histsdict['hyz res'] ,histsdict['hyz']  ,"hyzresnorm" ,"yz X-Residual normalized")
histsdict['hyz MCres norm'] =Divideh(histsdict['hyz MC res'] ,histsdict['hyz MC']  ,"hyzMCresnorm" ,"yz MC X-Residual normalized")

histsdict['hxyz res norm'] =Divideh(histsdict['hxyz res'] ,histsdict['hxyz']  ,"hxyzresnorm" ,"xyz X-Residual normalized")

#Y-Res-Norm
histsdict['hxy resY norm']=Divideh(histsdict['hxy resY'],histsdict['hxy'],"hxyresYnorm","xy Y-Residual normalized")
histsdict['hxy MCresY norm']=Divideh(histsdict['hxy MC resY'],histsdict['hxy MC'],"hxyMCresYnorm","xy MC Y-Residual normalized")

histsdict['hyz resY norm'] =Divideh(histsdict['hyz resY'] ,histsdict['hyz']  ,"hyzresYnorm" ,"yz Y-Residual normalized")
histsdict['hyz MCresY norm'] =Divideh(histsdict['hyz MC resY'] ,histsdict['hyz MC']  ,"hyzMCresYnorm" ,"yz MC Y-Residual normalized")

histsdict['hxyz resY norm'] =Divideh(histsdict['hxyz resY'] ,histsdict['hxyz']  ,"hxyzresnormY" ,"xyz Y-Residual normalized")

#Y-Res-Norm
histsdict['hxy resZ norm']=Divideh(histsdict['hxy resZ'],histsdict['hxy'],"hxyresZnorm","xy Z-Residual normalized")
histsdict['hxy MCresZ norm']=Divideh(histsdict['hxy MC resZ'],histsdict['hxy MC'],"hxyMCresZnorm","xy MC Z-Residual normalized")

histsdict['hyz resZ norm'] =Divideh(histsdict['hyz resZ'] ,histsdict['hyz']  ,"hyzresZnorm" ,"yz Z-Residual normalized")
histsdict['hyz MCresZ norm'] =Divideh(histsdict['hyz MC resZ'] ,histsdict['hyz MC']  ,"hyzMCresZnorm" ,"yz MC Z-Residual normalized")

histsdict['hxyz resZ norm'] =Divideh(histsdict['hxyz resZ'] ,histsdict['hxyz']  ,"hxyzresnormZ" ,"xyz Z-Residual normalized")


histsdict['hxy chi2X norm'] =Divideh(histsdict['hxy chi2X'] ,histsdict['hxy']  ,"hxychi2Xnorm" ,"xy chi2X normalized")
histsdict['hxy resP norm']=Divideh(histsdict['hxy resP'] ,histsdict['hxy']  ,"hxyresnormP" ,"xy P-Residual normalized")
histsdict['hxy resR norm']=Divideh(histsdict['hxy resR'] ,histsdict['hxy']  ,"hxyresnormR" ,"xy R-Residual normalized")


histsdict['hxy resU norm']=Divideh(histsdict['hxy resU'] ,histsdict['hxy']  ,"hxyresnormU" ,"xy U-Residual normalized")
histsdict['hxy resV norm']=Divideh(histsdict['hxy resV'] ,histsdict['hxy']  ,"hxyresnormV" ,"xy V-Residual normalized")


histsdict['hxy MCresP norm']=Divideh(histsdict['hxy MC resP'] ,histsdict['hxy']  ,"hxyMCresnormP" ,"xy MC P-Residual normalized")
histsdict['hxy MCresR norm']=Divideh(histsdict['hxy MC resR'] ,histsdict['hxy']  ,"hxyMCresnormR" ,"xy MC R-Residual normalized")


histsdict['hxy MCresU norm']=Divideh(histsdict['hxy MC resU'] ,histsdict['hxy']  ,"hxyMCresnormU" ,"xy MC U-Residual normalized")
histsdict['hxy MCresV norm']=Divideh(histsdict['hxy MC resV'] ,histsdict['hxy']  ,"hxyMCresnormV" ,"xy MC V-Residual normalized")


hslices['hxy res norm slice']={}
hslices['hxy resY norm slice']={}
hslices['hxy resZ norm slice']={}

hslices['hxy resR norm slice']={}
hslices['hxy resP norm slice']={}
hslices['hxy resU norm slice']={}
hslices['hxy resV norm slice']={}

hslices['hxy MC res norm slice']={}
hslices['hxy MC resY norm slice']={}
hslices['hxy MC resZ norm slice']={}

hslices['hxy MC resR norm slice']={}
hslices['hxy MC resP norm slice']={}
hslices['hxy MC resU norm slice']={}
hslices['hxy MC resV norm slice']={}
hslices['hxy resX norm phi slice']={}
hslices['hxy MC resX norm phi slice']={}
for i in range(len(zSlices)):
    hslices['hxy res norm slice'][i]=Divideh(hslices['hxy res slice'][i],hslices['hxy slice'][i],"hxyresnorm"+str(zSlices[i]),"xy X-residual normalized "+str(zSlices[i]))

    hslices['hxy resY norm slice'][i]=Divideh(hslices['hxy resY slice'][i],hslices['hxy Y slice'][i],"hxyresYnorm"+str(zSlices[i]),"xy Y-residual normalized")
    hslices['hxy resZ norm slice'][i]=Divideh(hslices['hxy resZ slice'][i],hslices['hxy Z slice'][i],"hxyresZnorm"+str(zSlices[i]),"xy Z-residual normalized")
    hslices['hxy resR norm slice'][i]=Divideh(hslices['hxy resR slice'][i],hslices['hxy all slice'][i],"hxyresRnorm"+str(zSlices[i]),"xy R-residual normalized")
    hslices['hxy resP norm slice'][i]=Divideh(hslices['hxy resP slice'][i],hslices['hxy all slice'][i],"hxyresPnorm"+str(zSlices[i]),"xy P-residual normalized")
    hslices['hxy resU norm slice'][i]=Divideh(hslices['hxy resU slice'][i],hslices['hxy all slice'][i],"hxyresUnorm"+str(zSlices[i]),"xy U-residual normalized")
    hslices['hxy resV norm slice'][i]=Divideh(hslices['hxy resV slice'][i],hslices['hxy all slice'][i],"hxyresVnorm"+str(zSlices[i]),"xy V-residual normalized")

    hslices['hxy MC res norm slice'][i]=Divideh(hslices['hxy MC res slice'][i],hslices['hxy slice'][i],"hxyMCresnorm"+str(zSlices[i]),"xy MC X-residual normalized "+str(zSlices[i]))
    hslices['hxy MC resY norm slice'][i]=Divideh(hslices['hxy MC resY slice'][i],hslices['hxy slice'][i],"hxyMCresYnorm"+str(zSlices[i]),"xy MC Y-residual normalized "+str(zSlices[i]))
    hslices['hxy MC resZ norm slice'][i]=Divideh(hslices['hxy MC resZ slice'][i],hslices['hxy slice'][i],"hxyMCresZnorm"+str(zSlices[i]),"xy MC Z-residual normalized "+str(zSlices[i]))
    hslices['hxy MC resR norm slice'][i]=Divideh(hslices['hxy MC resR slice'][i],hslices['hxy all slice'][i],"hxyMCresRnorm"+str(zSlices[i]),"xy MC R-residual normalized "+str(zSlices[i]))
    hslices['hxy MC resP norm slice'][i]=Divideh(hslices['hxy MC resP slice'][i],hslices['hxy all slice'][i],"hxyMCresPnorm"+str(zSlices[i]),"xy MC P-residual normalized "+str(zSlices[i]))
    hslices['hxy MC resU norm slice'][i]=Divideh(hslices['hxy MC resU slice'][i],hslices['hxy all slice'][i],"hxyMCresUnorm"+str(zSlices[i]),"xy MC U-residual normalized "+str(zSlices[i]))
    hslices['hxy MC resV norm slice'][i]=Divideh(hslices['hxy MC resV slice'][i],hslices['hxy all slice'][i],"hxyMCresVnorm"+str(zSlices[i]),"xy MC V-residual normalized "+str(zSlices[i]))

for i in range(len(phiSlices)):
    hslices['hxy resX norm phi slice'][i]=Divideh(hslices['hxy resX phi slice'][i],hslices['hxy slice phi'][i],"hxyresnorm_phi"+str(phiSlices[i]),"xy X-residual normalized phi"+str(phiSlices[i]))
    hslices['hxy MC resX norm phi slice'][i]=Divideh(hslices['hxy MC resX phi slice'][i],hslices['hxy slice phi'][i],"hxyMCresnorm_phi"+str(phiSlices[i]),"xy MC X-residual normalized phi"+str(phiSlices[i]))
 

for ivar in range(len(digiPerpVariance)):
    histsdict['digi residuals']['digi variance'].SetBinContent(ivar,digiPerpVariance[ivar])
    histsdict['digi residuals']['ndigis'].SetBinContent(ivar,digiNum[ivar])

outfile = ROOT.TFile(args.outfname,"recreate")
outfile.cd()

outfile.mkdir('info')
outfile.cd('info')

infos=args.info.split(',')
rinfos={}
rinfos['generator'] =ROOT.TNamed('generator','undefined')
rinfos['angle']     =ROOT.TNamed('angle','undefined')
rinfos['error']     =ROOT.TNamed('error','undefined')
rinfos['correction']=ROOT.TNamed('correction','undefined')
rinfos['field']     =ROOT.TNamed('field','undefined')
rinfos['gas']       =ROOT.TNamed('gas','undefined')
rinfos['gain']      =ROOT.TNamed('gain','undefined')
rinfos['name']      =ROOT.TNamed('name','undefined')
if haspar:
    rinfos['Dl']        =ROOT.TNamed('Dl',str(digiPar.getGas().Dl()))
    rinfos['Dt']        =ROOT.TNamed('Dt',str(digiPar.getGas().Dt()))
    rinfos['vd']        =ROOT.TNamed('vd',str(digiPar.getGas().VDrift()))
    rinfos['gain']      =ROOT.TNamed('gain',str(digiPar.getGain()))
    rinfos['field']     =ROOT.TNamed('field',str(digiPar.getEField()))

if len(infos)>2:
    for inf in infos:
        inf=inf.split('=')
        rinfos[inf[0]]=ROOT.TNamed(inf[0],inf[1])
        
for rinf in rinfos:
    rinfos[rinf].Write()
outfile.cd()

histcounter = 0
outfile.mkdir('single')
for hists in histsdict:
    histcounter+=1
    hist=histsdict[hists]
    outfile.cd()
    if type(hist)==dict:
        outfile.mkdir(hists)
        outfile.cd(hists)
        for h in hist:
            hist[h].Write()
    else:
        outfile.cd('single')
        hist.Write()


outfile.cd()
outfile.mkdir('Z-slices')
for d in hslices:
    sdict=hslices[d]
    outfile.cd()
    outfile.mkdir('Z-slices/'+d)
    outfile.cd('Z-slices/'+d)
    for i in range(len(zSlices)):
        h=sdict[i]
        h.Write()


outfile.cd()

'''
if not args.nopostprocess:
    #slice through pulls, fit with double gauss and save the results:
    pullprojs={}
    pullgraphs={}
    outfile.mkdir('pulls/')
    outfile.mkdir('pulls/graphs')
    outfile.mkdir('pulls/projs')
    for h in histsdict['pulls']:
        pullprojs[h]=[]
        sliceopt=0
        if h.find('err')!=-1:
            sliceopt=1    
        if h.find('vs z')==-1:
            continue
        pullgraphs[h]=getSlices(histsdict['pulls'][h],pullprojs[h],sliceopt)
        outfile.cd('pulls/graphs/')
        pullgraphs[h]['Mean'].Write()
        pullgraphs[h]['RMS'].Write()
        if( pullgraphs[h].get('Const',None)!=None):
            pullgraphs[h]['Const'].Write()
        outfile.cd()
        outfile.mkdir('pulls/projs/'+h)
        outfile.cd('pulls/projs/'+h)
        for proj in pullprojs[h]:
            proj.Write()

    gDigiVariance=ROOT.TGraph()
    gDigiVariance.SetName('graph_Digivariance')
    gDigidiffusion=ROOT.TGraph()
    gDigidiffusion.SetName('graph_Digidiffusion')
    gDigivarDifFrac=ROOT.TGraph()
    gDigivarDifFrac.SetName('graph_DigivarDifFrac')

    Dt=digiPar.getGas().Dt()
    for ivar in range(len(digiPerpVariance)):
        if digiNum[ivar]!=0:
            digiPerpVariance[ivar]=digiPerpVariance[ivar]/digiNum[ivar]
        else:
            digiPerpVariance[ivar]=0
        zsig=ROOT.TMath.Sqrt(digiPerpVariance[ivar])
        gDigiVariance.SetPoint(ivar,ivar+0.5,zsig)
        elecDiff=Dt*ROOT.TMath.Sqrt(ivar+0.5)
        gDigidiffusion.SetPoint(ivar,ivar+0.5,elecDiff)
        if zsig!=0:
            gDigivarDifFrac.SetPoint(ivar,ivar+0.5,elecDiff/zsig)

    outfile.cd()
    outfile.mkdir('gemSpread')
    outfile.cd('gemSpread')
    gDigiVariance.Write()
    gDigidiffusion.Write()
    gDigivarDifFrac.Write()
'''
outfile.cd()
outfile.Close()

#histsdict['pulls']['cluster x vs z']

print 'Histos are written into file:', args.outfname
print 'total runtime: {0} '.format(float((time.time()-main_start)))