#!/usr/bin/python
import os,sys
from ROOT import *
from clean_canvas import *
from nice_colours import *
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')
f = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum68-3-2_zone3_intercept30_raw.root')
#f = TFile('/home/e614/analysis_2006_2007/treesum/pass1/Sum2/PsPact/sum68-3-2_zone4_intercept30_raw.root')
hi = f.Get('PACT/pc6vs5_11_after')

c = TCanvas()
c.Divide(2,2)
c.cd(1)
hi.SetTitle('; PC5 after rotation (arb units); PC6 after rotation (arb units)')
hi.Draw('colz')
c.cd(2)
px = hi.ProjectionX().Clone()
px.Clone().Draw()
c.cd(3)

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

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

#print 'nentries in 2D histo = ',hi.GetEntries()
nbx = px.GetNbinsX()
t = 0.
for ix in range(1,nbx+1):
    t += px.GetBinContent(ix)
print 'total content 1D projection = ',t
print 'integral for fit function from 129 to 500  = ',ff.Integral(129.,500.)
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')


au = raw_input('>')