import sys
sys.path.append('/home/vwalbr/fopiroot/fopiroot_dev/macro/tpc/FOPI/vwalbre')
import ROOT,glob,math,argparse
from anaFile import anaFile
from functions import *

def RooFit(histo,outfname,c):
    c.Divide(2,1)
    x=ROOT.RooRealVar("x", "x", -10,10)
    varset=ROOT.RooArgList('set')
    varset.add(x)
    #PreFit
    DataSetPre=ROOT.RooDataHist("DataSet", "Dataset with x",varset, histo)

    meanLPre=ROOT.RooRealVar('meanLPre','Mean of  gaussian L Pre',0,-0.5 ,0.5)
    sigLPre=ROOT.RooRealVar('sigLPre','Sigma of gaussian L Pre',2.,1.,3.)

    width=ROOT.RooRealVar('width','width',3.,2.,4.)

    meanRPre=ROOT.RooRealVar('meanRPre','Mean of  gaussian R Pre',0,-0.5,0.5)
    sigRPre=ROOT.RooRealVar('sigRPre','Sigma of gaussian R Pre',0.5,0.3,0.8)
    width2=ROOT.RooRealVar('width2','width2',0.6,0.3,.9)

    gausLPre=ROOT.RooVoigtian("gausLPre","gausLPre p.d.f.",x,meanLPre,width,sigLPre)
 
    gausRPre=ROOT.RooVoigtian("gausRPRe","gausR p.d.f.",x,meanRPre,width2,sigRPre) 

    c0 = ROOT.RooRealVar("c0","cof 0", 1,-1.,1.)
    c1 = ROOT.RooRealVar("c1","cof 1", 0.1,-1.,1.)
    c2 = ROOT.RooRealVar("c2","cof 2", -0.1,-1.,1.)
    bkg=ROOT.RooChebychev("bkg","bkg",x,ROOT.RooArgList(c0,c1,c2))

    gm1frac=ROOT.RooRealVar("gm1frac","fraction of gm1",0.5,0.,1.0)
    gm2frac=ROOT.RooRealVar("gm2frac","fraction of gm2",0.5,0,1.0)
    Sig=ROOT.RooAddModel("Sig","Sig",ROOT.RooArgList(gausRPre,gausLPre),ROOT.RooArgList(gm1frac))
    DGaus=ROOT.RooAddModel("DGaus","DoubleGaus",ROOT.RooArgList(gausLPre,bkg),ROOT.RooArgList(gm2frac))
    
    #RelDiff1=gausLPre.fitTo(DataSetPre,ROOT.RooFit.Range(-2,0))
    #RelDiff1=gausRPre.fitTo(DataSetPre,ROOT.RooFit.Range(-10,10))
    #RelDiff2=gausRPre.fitTo(DataSetPre,ROOT.RooFit.Range(0 ,2))

    #DoppleFit
    #DataSet=ROOT.RooDataHist("DataSet", "Dataset with x",varset, histo)
    #meanL=ROOT.RooRealVar('meanL','Mean of  gaussian L',meanLPre.getVal(),meanLPre.getVal()-0.4,meanLPre.getVal()+0.4)
    #sigL=ROOT.RooRealVar('sigL','Sigma of gaussian L',sigLPre.getVal(),sigLPre.getVal()-0.4,sigLPre.getVal()+0.4)
    #meanR=ROOT.RooRealVar('meanR','Mean of  gaussian R',meanRPre.getVal(),meanRPre.getVal()-0.4,meanRPre.getVal()+0.4)
    #sigR=ROOT.RooRealVar('sigR','Sigma of gaussian R',sigRPre.getVal(),sigRPre.getVal()-0.4,sigRPre.getVal()+0.4)
    #gausL=ROOT.RooGaussian("gausL","gausL p.d.f.",x,meanL,sigL) 
    #gausR=ROOT.RooGaussian("gausR","gausR p.d.f.",x,meanR,sigR) 

    #gm1frac=ROOT.RooRealVar("gm1frac","fraction of gm1",0.5,0,1.)
    #DGaus=ROOT.RooAddModel("DGaus","DoubleGaus",ROOT.RooArgList(gausL,gausR),ROOT.RooArgList(gm1frac))
    #print meanLPre.getVal(), meanRPre.getVal()
    #test= float(round(meanRPre.getVal(),2))
    #test2= float(abs(round(meanLPre.getVal(),2)))
    #if test<0.30 or test2<0.3:
    #    LRange=-sigLPre.getVal()-0.4
    #    RRange=sigRPre.getVal()+0.4
    #elif test>=1.0:
    #    LRange=-sigLPre.getVal()-1.
    #    RRange=sigRPre.getVal()+1.
    #else:
    #    LRange=-sigLPre.getVal()-0.6
    #    RRange=sigRPre.getVal()+0.6
    #RelDiffL=DGaus.fitTo(DataSet,ROOT.RooFit.Range(LRange,RRange))
    c.cd()
    xframe=x.frame()
    DataSetPre.plotOn(xframe)
    #gausLPre.plotOn(xframe,ROOT.RooFit.LineColor(2))
    #gausL.plotOn(xframe,ROOT.RooFit.LineColor(2))
    #gausR.plotOn(xframe,ROOT.RooFit.LineColor(2))
    #gausLPre.paramOn(xframe,ROOT.RooFit.Layout(0.65,0.65,0.89))
    xframe.getAttText().SetTextSize(0.02)
    xframe.Draw()
    #gausR.paramOn(xframe)
    c.Update()
    FitPar=[]
    #FitPar.append(meanL.getVal())
    #FitPar.append(sigL.getVal())
    #FitPar.append(meanR.getVal())
    #FitPar.append(sigR.getVal())
    #FitPar.append(gm1frac.getVal())
    #print test
    #print meanLPre.getVal(), sigLPre.getVal(), meanRPre.getVal(), sigRPre.getVal()
    #print LRange, RRange
    return FitPar, c

def prod_reldiff(hdiff, reldiff, reldiff1d, xyres, xy):
    reldiff.Reset()
    for ibi in range(hdiff.GetNbinsX()+1):
        for jbi in range(hdiff.GetNbinsY()+1):
            #if xy[0].GetBinContent(ibi,jbi)== 0 or xy[1].GetBinContent(ibi,jbi)==0:
            #    continue
            cont=xyres[0].GetBinContent(ibi,jbi)-xyres[1].GetBinContent(ibi,jbi)
            hdiff.SetBinContent(ibi,jbi,cont)
            divider=xyres[1].GetBinContent(ibi,jbi)
            newcont=-9999999999
            if divider > -99999 and cont > -99999:
                newcont=cont/divider
                reldiff1d.Fill(newcont)
                #newcont=abs(newcont)
                reldiff.SetBinContent(ibi,jbi,newcont)
            
    histos=[]
    histos.append(reldiff)
    histos.append(reldiff1d)
    
    return histos 

def prod_reldiff_slice(counter, hdiffbin, xyslice, reldiffbin, reldiffbin1d):
    reldiffbin[-1].Reset()
    for ibi in range(hdiffbin[-1].GetNbinsX()+1):
        for jbi in range(hdiffbin[-1].GetNbinsY()+1):
            cont=hdiffbin[-1].GetBinContent(ibi,jbi)
            divider=xyslice[1][counter].GetBinContent(ibi,jbi)
            newcont=-9999999999
            #if divider==0:
            #    continue
            if divider > -999 and cont > -999:
                newcont=cont/divider
                reldiffbin1d[-1].Fill(newcont)
                #if newcont==1:
                #    print ibi,jbi,cont,divider
                    #newcont=abs(newcont)
                reldiffbin[-1].SetBinContent(ibi,jbi,newcont)
    histos_slice=[]
    histos_slice.append(reldiffbin)
    histos_slice.append(reldiffbin1d)

    return histos_slice

def set_titles(xyres):
    xyres.GetXaxis().SetTitle("X Pos (cm)")
    xyres.GetYaxis().SetTitle("Y Pos (cm)")
    xyres.GetZaxis().SetTitle("Residual (#mu m)")
    xyres.GetZaxis().SetTitleOffset(1.6)
    xyres.SetStats(0)
    return xyres

def set_titles3D(xyzres):
    xyzres.GetXaxis().SetTitle("X Pos (cm)")
    xyzres.GetYaxis().SetTitle("Y Pos (cm)")
    xyzres.GetZaxis().SetTitle("Z Pos (cm)")
    xyzres.GetZaxis().SetTitleOffset(1.6)
    xyzres.SetStats(0)
    return xyzres

def set_titles_hdiff(hdiff):
    hdiff.GetZaxis().SetRangeUser(-1000,1000)
    hdiff.GetXaxis().SetTitle("X Pos (cm)")
    hdiff.GetYaxis().SetTitle("Y Pos (cm)")
    hdiff.GetZaxis().SetTitle("Difference (#mu m)")
    hdiff.GetZaxis().SetTitleOffset(1.6)
    hdiff.SetStats(0)
    return hdiff

def set_titles_hdiffyz(hdiff):
    hdiff.GetZaxis().SetRangeUser(-1000,1000)
    hdiff.GetXaxis().SetTitle("Z Pos (cm)")
    hdiff.GetYaxis().SetTitle("Y Pos (cm)")
    hdiff.GetZaxis().SetTitle("Difference (#mu m)")
    hdiff.GetZaxis().SetTitleOffset(1.6)
    hdiff.SetStats(0)
    return hdiff

def set_titles_hdiff3D(hdiff3D):
    hdiff3D.GetXaxis().SetTitle("X Pos (cm)")
    hdiff3D.GetYaxis().SetTitle("Y Pos (cm)")
    hdiff3D.GetZaxis().SetTitle("Difference (#mu m)")
    hdiff3D.GetZaxis().SetTitleOffset(1.6)
    return hdiff3D

def prod_reldiff3D(hdiff, reldiff,reldiff1d, xyzres,xyz):
    reldiff.Reset()

    hdiff2=hdiff.Clone("hdiff2")

    RelDiffyx=ROOT.TH2D("RelDiffyx","RelDiff proj yx ",32,-16.,16.,32,-16.,16.)
    binyx=0
    countyx=0

    RelDiffCUTyx=ROOT.TH2D("RelDiffCUTyx","RelDiff CUT proj yz ",32,-16.,16.,32,-16.,16.)
    binyxCUT=0
    countyxCUT=0

    Zero=0
    
    for ibi in range(hdiff.GetNbinsX()+1):        
        for jbi in range(hdiff.GetNbinsY()+1):
            for kbi in range(hdiff.GetNbinsZ()+1):
                if xyz[0].GetBinContent(ibi,jbi,kbi)== 0 or xyz[1].GetBinContent(ibi,jbi,kbi)==0:
                    continue
                if kbi == 30:
                #    print xyz[1].GetBinContent(ibi,jbi,kbi)
                    RelDiffyx.SetBinContent(ibi,jbi,xyz[1].GetBinContent(ibi,jbi,kbi))
                    RelDiffCUTyx.SetBinContent(ibi,jbi,xyz[0].GetBinContent(ibi,jbi,kbi))
                cont=xyzres[0].GetBinContent(ibi,jbi,kbi)-xyzres[1].GetBinContent(ibi,jbi,kbi)
                hdiff.SetBinContent(ibi,jbi,cont)          
                divider=xyzres[1].GetBinContent(ibi,jbi,kbi)              
                newcont=-9999999999
                #if divider ==0: 
                #    continue
                if divider > -999 and cont > -999:
                    newcont=cont/divider
                    if newcont==-1.0:
                        print ibi, kbi, jbi
                    reldiff1d.Fill(newcont)
                    binyx+=cont
                    countyx+=1
                    if abs(newcont)>2:
                        binyxCUT+=cont
                        countyxCUT+=1
                    reldiff.SetBinContent(ibi,jbi,kbi,cont)
                    #if newcont==-1.0:
                    #    if countyx!=0:
                            #bini=binyx/countyx
                            #RelDiffyx.SetBinContent(ibi,jbi,newcont)
            if countyxCUT!=0:
                biniCUT=binyxCUT/countyxCUT
                #RelDiffCUTyx.SetBinContent(ibi,jbi,biniCUT)
            binyx=0
            countyx=0
            binyxCUT=0
            countyxCUT=0

    print "**********", Zero

    RelDiffyz=ROOT.TH2D("RelDiffyz","RelDiff proj yz ",39,-4,74,32,-16.,16.)
    binyz=0
    countyz=0
    
    RelDiffCUTyz=ROOT.TH2D("RelDiffCUTyz","RelDiff CUT proj xy ",39,-4,74,32,-16.,16.)
    binyzCUT=0
    countyzCUT=0
    
    for ibi in range(hdiff2.GetNbinsZ()+1):        
        for jbi in range(hdiff2.GetNbinsY()+1):
            for kbi in range(hdiff2.GetNbinsX()+1):
                #if xyz[0].GetBinContent(kbi,jbi,ibi)== 0 or xyz[1].GetBinContent(kbi,jbi,ibi)==0:
                #    continue
                cont=xyzres[0].GetBinContent(kbi,jbi,ibi)-xyzres[1].GetBinContent(kbi,jbi,ibi)       
                divider=xyzres[1].GetBinContent(kbi,jbi,ibi)      
                if divider==0:
                    continue
                newcont=-9999999999
                if divider > -999 and cont > -999:
                    newcont=cont/divider
                    binyz+=cont
                    countyz+=1
                    if abs(newcont)>2:
                        binyzCUT+=cont
                        countyzCUT+=1
            if countyz!=0:
                bini=binyz/countyz
                RelDiffyz.SetBinContent(ibi,jbi,bini)
            if countyzCUT!=0:
                biniCUT=binyzCUT/countyzCUT
                RelDiffCUTyz.SetBinContent(ibi,jbi,biniCUT)
            binyz=0
            countyz=0
            binyzCUT=0
            countyzCUT=0

    histos=[]
    histos.append(reldiff)
    histos.append(reldiff1d)


    return histos, RelDiffyx, RelDiffCUTyx, RelDiffyz, RelDiffCUTyz