#!/usr/bin/python
import os,sys
from math import sqrt
from ROOT import *
from clean_canvas import *
gStyle = makegstyle(gStyle)

me = 0.511
mmu = 105.658
Wem = (mmu**2 + me**2)/(2*mmu)
x0 = me/Wem

## __________
def frange(start,end,nintervals):
    step = (float(end)-float(start))/float(nintervals-1)
    return [start+step*i for i in range(nintervals)]

## __________
def F_IS(x,rho,eta):
    return (x*(1-x)
            + (2./9.)*rho*(4*x*x - 3*x - x0**2)
            + eta*x0*(1-x))
## __________
def F_AS(x,xi,delta):
    return ((1./3.)*(xi*sqrt(x**2-x0**2))*
            (1-x+(2./3.)*delta*
             (4*x - 3. + (sqrt(1-x0**2)-1))))
## __________
def k(x): return sqrt(x**2-x0**2)

## __________
def git(vx,vy,ti):
    gr = TGraph()
    gr.SetTitle(ti)
    for i,(px,py) in enumerate(zip(vx,vy)):
        gr.SetPoint(i,px,py)
    return gr


## ____________________
if __name__=='__main__':

    Pmu, xi, delta, rho, eta = 1., 1., 0.75, 0.75, 0.

    a, b = -0.2, 0.2

    vs = {'xs':[], 'Nus':[], 'Nds':[], 'sum':[], 'diff':[], 'asy':[]}
    for x in frange(0.5,1.,10):
        Nu = k(x) * ( F_IS(x,rho,eta)*(b-a) + Pmu * F_AS(x,xi,delta) * (b**2-a**2) )
        Nd = k(x) * ( F_IS(x,rho,eta)*((-a)-(-b)) + Pmu * F_AS(x,xi,delta) * (a**2-b**2) )
        print 'x=',x,'Nu=',Nu,'Nd=',Nd
        vs['xs'].append(x)
        vs['Nus'].append(Nu)
        vs['Nds'].append(Nd)
        vs['sum'].append(Nu+Nd)
        vs['diff'].append(Nu-Nd)
        vs['asy'].append( (Nu-Nd)/(Nu+Nd) )

    grs = [git(vs['xs'],vs['Nus'],'Nus'),
           git(vs['xs'],vs['Nds'],'Nds'),
           git(vs['xs'],vs['sum'],'sum'),
           git(vs['xs'],vs['diff'],'diff'),
           git(vs['xs'],vs['asy'],'asy')]

    cs = []
    for i,gr in enumerate(grs):
        cs.append(TCanvas(gr.GetTitle(),gr.GetTitle()))
        gr.Draw('AC')
    au = raw_input('>')