#!/usr/bin/python
import os,sys
from ROOT import *
#from clean_canvas import *
from nice_colours import *

def makegstyle(cStyle):
    cStyle.SetPadBorderMode(0)
    cStyle.SetOptStat(0)	   ## Hide stats
    cStyle.SetOptTitle(1)	   ## Display histogram titles.
    cStyle.SetCanvasColor(10)	   ## Background color (white)
    cStyle.SetPadColor(10)	   ## Background color (white)
    cStyle.SetFillColor(10)	   ## Background color (white)
    cStyle.SetTitleFillColor(10)   ## Background color (white)
    cStyle.SetStatColor(10)	   ## Background color (white)
    cStyle.SetPalette(1)	   ## Nicer colour scale for 2D histograms.  
    cStyle.SetTitleBorderSize(0)   ##duh
    cStyle.SetOptFit(1) #added
    return cStyle

gStyle = makegstyle(gStyle)
gStyle = nice_colours(gStyle)

## __________
def gauss_with_exp_tails(x,par):
    ## copied from /home/rpm1/Root/tailgaussianfit.C
    ## Arguments: x is the array of variables (only one element for a 1D
    ## function like this one, so we always refer to x[0]), and par is
    ## the array of fit parameters.
    ## The five fit parameters are, in order:
    ## N, lambdaneg, lambdapos, sigma, mean
    if ((x[0] - par[4]) < -par[1]*par[3]):
        ## Exponential tail
        return par[0]*exp(par[1]*par[1]/2.0)*exp(par[1]*(x[0] - par[4])/par[3])
    elif ((x[0] - par[4]) > par[2]*par[3]):
        ## Exponential tail (notice the sign flip!)
        return par[0]*exp(par[2]*par[2]/2.0)*exp(-par[2]*(x[0] - par[4])/par[3])
    else:
        ## Gaussian peak
        return par[0]*exp(-(x[0] - par[4])*(x[0] - par[4])/(2.0*par[3]*par[3]));

#f = TFile('/home/e614/jbueno/psPACT_systematic_rotated_PC56/results/sum68-3-2_clk_raw.root')

f4 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum68-3-2_zone4_intercept30_raw.root')
f3 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum68-3-2_zone3_intercept30_raw.root')
#f4 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum74-4-2_zone4_lower_intercept_raw.root')
#f3 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum74-4-2_zone3_lower_intercept_raw.root')
#f4 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum83-4-2_zone4_lower_intercept_raw.root')
#f3 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum83-4-2_zone3_lower_intercept_raw.root')
#f4 = TFile('/home/e614/jbueno/PC5678_tests/s83_PC6_Q4_raw.root')  #NOT ROTATED, no adjustedment to intercept
#f3 = TFile('/home/e614/jbueno/PC5678_tests/s83_PC6_Q3_raw.root')  #NOT ROTATED, no adjustedment to intercept
#f4 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum87-2-2_zone4_lower_intercept_raw.root')
#f3 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum87-2-2_zone3_lower_intercept_raw.root')

#f4 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/s84_9_2_Q4_raw.root')
#f3 = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/s84_9_2_Q3_raw.root')

h4 = f4.Get('PACT/pc6vs5_11_after')
h3 = f3.Get('PACT/pc6vs5_11_after')

h0 = f3.Get('PACT/pc6vs5_11_before')

h4.SetTitle('ZONE 4 - rotated; PC5 after rotation (arb units); PC6 after rotation (arb units)')
h3.SetTitle('ZONE 3 - rotated; PC5 after rotation (arb units); PC6 after rotation (arb units)')
h0.SetTitle('All zones, BEFORE rotation')

for h in [h3,h4]:
    h.GetXaxis().SetRangeUser(0,300.)
    h.GetYaxis().SetRangeUser(0,300.)

c = TCanvas()
c.Divide(2,2)

#c.cd(1)
#h0.Draw('colz')

c.cd(1)
h4.Draw('colz')
c.cd(2)
h3.Draw('colz')


px4 = h4.ProjectionX().Clone()
gROOT.cd()
px3 = h3.ProjectionX().Clone()
gROOT.cd()

px4.SetLineColor(2)
px3.SetLineColor(4)

px3.GetXaxis().SetRangeUser(0,300.)
px3.SetTitle('projection of zones 4 and 3')

c.cd(3)

ff = TF1('gauss_with_exp_tails',gauss_with_exp_tails, 127.,200.,5)
ff.SetParNames('maxval','lambda_LHS','lambda_RHS','sigma','mean')
ff.SetParameters(px3.GetMaximum(), 1.0, 1.0, px3.GetRMS(), 150)
px3.Fit(ff,'EMR','',127,200)

ff.SetLineColor(1)
ff.SetLineStyle(3)

ff2 = TF1('gauss_with_exp_tails2',gauss_with_exp_tails, 0.,127.,5)
ff2.SetParameters(ff.GetParameter(0),ff.GetParameter(1),
                  ff.GetParameter(2),ff.GetParameter(3),
                  ff.GetParameter(4))
ff2.SetLineStyle(7)

px3.Clone().Draw()
px4.Clone().Draw('same')


#ff2.SetLineWidth(0.8)

ff.Draw('same')
ff2.Draw('same')

c.cd(4)

ff3 = TF1('gaus1','gaus',127.,200.)
px3.Fit(ff3,'EMR','',127,200)

ff4 = TF1('gaus2','gaus',0.,127.)
ff4.SetParameters(ff3.GetParameter(0),ff3.GetParameter(1),ff3.GetParameter(2))
ff4.SetLineStyle(7)

px3.Clone().Draw()
px4.Clone().Draw('same')


#ff3.SetLineWidth(0.8)
#ff4.SetLineWidth(0.8)

ff3.Draw('same')
ff4.Draw('same')


print '_____ gaus with exp tails _____'
print 'integral for fit function from 127 to 500   = ',ff.Integral(127.,500.)
print 'integral for fit function from -500 to 127  = ',ff.Integral(-500.,129.)

print '_____ gaus _____'
print 'integral for fit function from 127 to 500   = ',ff3.Integral(127.,500.)
print 'integral for fit function from -500 to 127  = ',ff3.Integral(-500.,129.)

##print 'integral for fit function from -500 to 129 = ',ff2.Integral(-500.,129.)
##
###print ff2.Integral(0.,129)
##
##c.cd(4)
##px.Draw()
##ff2.Draw('same')


def onegaus(x,par):
    g1 = par[0]*exp(-0.5*((x[0]-par[1])/par[2])**2)
    return g1
def twogaus(x,par):
    g1 = par[0]*exp(-0.5*((x[0]-par[1])/par[2])**2)
    g2 = par[3]*exp(-0.5*((x[0]-par[4])/par[5])**2)
    return g1+g2

c2 = TCanvas()
c2.Divide(1,2)

c2.cd(1)

h43 = h4.Clone()
h43.Add(h3,1.)
h43.GetXaxis().SetRangeUser(60.,300.)
h43.SetTitle('zones 4 and 3 combined')
h43.Draw('colz')

c2.cd(2)

p43 = h43.ProjectionX().Clone()
p43.SetTitle('zones 4 and 3 combined')
p43.GetXaxis().SetRangeUser(60.,300.)

## _____ CHANGE ME _____
ff5 = TF1('twogaus1',twogaus,127.,200.,6)   #s68
#ff5 = TF1('twogaus1',twogaus,85.,175.,6)    #s74
#ff5 = TF1('twogaus1',twogaus,100.,200.,6)   #s83, s87

gStyle.SetOptFit(1)

ff5.SetParNames('scale1','mean1','sigma1','scale2','mean2','sigma2')

## _____ CHANGE ME _____
ff5.SetParameters(0.1*p43.GetMaximum(),100.,15.,0.7*p43.GetMaximum(),150.,15.)  #s68, s74
#ff5.SetParameters(0.3*p43.GetMaximum(),125.,15.,0.7*p43.GetMaximum(),175.,15.)  #s83, s87

## _____ CHANGE ME _____
p43.Fit(ff5,'EMR','',85.,175.) #s68, s74
#p43.Fit(ff5,'EMR','',100.,200.) #s83, s87

p43.Draw()

ff6 = TF1('onegaus_gas',onegaus,100,200,3)
ff6.SetParameters(ff5.GetParameter(3),ff5.GetParameter(4),ff5.GetParameter(5))
ff6.SetLineStyle(3)

## _____ CHANGE ME _____
low = 127  #s68,s74
#low = 153  #s83, s87
#low = 136  #s84 test

l = TLine(low,0,low,13500)
l.SetLineWidth(3)
l.SetLineStyle(7)
l.Draw('same')

ff6.Draw('same')

i1 = ff6.Integral(-500.,low)
i2 = ff6.Integral(low,500.)

print 'integral for single gaus from -500 to',low,' =', i1
print 'integral for single gaus from',low,'to 500 =', i2
print 'contamination = ',i1/(i1+i2)

c3 = TCanvas()
p43.Draw()
ff6.Draw('same')

au = raw_input('>')