#!/usr/bin/python
import os,sys
from math import cos
from scipy.optimize import brentq  #recommended ROOT finder

## __________ TWIST inputs __________
Pmuxi_90pc = 0.99922  #90% conf. from Glen's script, which refs Carl's posting

## __________ PDG inputs __________
m_W1   = 80.398   #GeV/c^2; from Glen's script
V_ud_L = 0.97     #from Glen's script

## __________ magic __________
def f(zeta_g,t,t_theta,alpha,omega):
    #following are Eqs. 1.26 and 1.27 of JB's thesis

    #debugging this
#    Pmu = 1. - 2.*t_theta*t_theta \
#          - 2.*zeta_g*zeta_g\
#          - 4.*t_theta*zeta_g*cos(alpha+omega)
#    xi  = 1. - 2.*(t*t+zeta_g*zeta_g)
#    return (Pmu*xi - Pmuxi_90pc)  #will be zero for root
#    print zeta_g,t,t_theta,alpha,omega
    Pmuxi =1-4*t*t-4*zeta_g*zeta_g-4*t*zeta_g
    return Pmuxi-Pmuxi_90pc

def mixing_angle(gL_times_m2_over_gR):

    x = gL_times_m2_over_gR    #convenience
    t = (m_W1 * m_W1) / (x*x)  #Eq. 1.29 in JB's thesis

    V_ud_R = V_ud_L  #TEMP
    t_theta = t*abs(V_ud_R/V_ud_L)  #Eq. 1.30 in JB's thesis

    alpha = 0. #TEMP
    omega = 0. #TEMP

    try:
        plus_angle  = brentq(f, 0.  , 0.1, (t,t_theta,alpha,omega))
        minus_angle = brentq(f, -0.1, 0. , (t,t_theta,alpha,omega))
    except ValueError: #no solution in this region
        return None,None

    return plus_angle,minus_angle

if __name__=='__main__':

#    print f(0.0161,0.0100997,0.0100997,0,0)
#    plus_angle,minus_angle = mixing_angle(800.)
#    print f(plus_angle,0.0100997,0.0100997,0,0)
#    sys.exit(1)

#    plus_angle,minus_angle = mixing_angle(800.)
#    print 800,plus_angle,minus_angle
#    sys.exit(1)
 
    points = []
    for m2 in range(1,1200,1):
        plus_angle, minus_angle = mixing_angle(m2)
        if plus_angle != None:
            points.append((m2,plus_angle))
        if minus_angle != None:
            points.append((m2,minus_angle))

## __________ from here on, it's just plotting __________
#at this stage you can "print points" if you want to use a different application for printing

## _____ lipstick on pig _____
from ROOT import *
from clean_canvas import *
gStyle = makegstyle(gStyle)

## _____ make plot _____
g = TGraph()
for i,(x,y) in enumerate(points):
    g.SetPoint(i,x,y)
g.Draw('AP')
g.SetMarkerStyle(20)
au = raw_input('>')