! plot_twist_rhmu_lim_mono

device\color postscript a
set linthk -6

! use 90%CL
cl = 0.9

! sigma corresponding to this cl
! usual gaussian, 1.64@90%
sig_clg = sqrt(2)*aerf(cl)
! or one-sided limit, 1.28@90%
sig_clg1 = sqrt(2)*aerf(1 - 2*(1 - cl))

scale_xlo = 0.98
scale_xhi = 1.03
scale_ylo = 0.735
scale_yhi = 0.760
scale scale_xlo scale_xhi 5 scale_ylo scale_yhi 5
scale_ext = 0.02

scale_xlo_ext = scale_xlo - scale_ext
scale_xhi_ext = scale_xhi + scale_ext
scale_ylo_ext = scale_ylo - scale_ext
scale_yhi_ext = scale_yhi + scale_ext

x_xi = [scale_xlo;scale_xhi]
y_delta = [scale_ylo;scale_yhi]

! delta

! Balke et al., PR D37 587 (PDG)

delta_balke = 0.7486
ddelta_balke = 1.e-4*sqrt(26**2 + 28**2)
llim_delta_balke = delta_balke - sig_clg*ddelta_balke
ulim_delta_balke = delta_balke + sig_clg*ddelta_balke

! gaponenko et al.
delta_twist = 0.74964
ddelta_twist = 1.e-4*sqrt(6.6**2 + 11.2**2)
llim_delta_twist = delta_twist - sig_clg*ddelta_twist
ulim_delta_twist = delta_twist + sig_clg*ddelta_twist

!  NSERC research proposal, 2004, for 2005 result
delta_twist_05 = 0.7500
ddelta_twist_05 = 1.e-4*sqrt(5.3**2 + 2.8**2)
llim_delta_twist_05 = delta_twist_05 - sig_clg*ddelta_twist_05
ulim_delta_twist_05 = delta_twist_05 + sig_clg*ddelta_twist_05

!  NSERC research proposal, 2004, for final result
delta_twist_f = 0.7500
ddelta_twist_f = 1.e-4*sqrt(2.2**2 + 3.2**2)
llim_delta_twist_f = delta_twist_f - sig_clg*ddelta_twist_f
ulim_delta_twist_f = delta_twist_f + sig_clg*ddelta_twist_f

! pmu*xi

! Beltrami et al., PLB 194, 326 (1987)
pmuxi_bel = 1.0027
dpmuxi_bel = 1.e-4*sqrt(79**2 + 30**2)
llim_pmuxi_bel = pmuxi_bel - sig_clg*dpmuxi_bel
ulim_pmuxi_bel = pmuxi_bel + sig_clg*dpmuxi_bel
! test
! ulim_pmuxi_bel = 1.01

! pmu*xi*delta/rho

! Jodidio et al., original PRD article
pmuxidelrho_jod_o = 0.99863
dpmuxidelrho_jod_o = 1.e-4*sqrt(4.6**2 + 7.5**2)

! Jodidio et al., PRD erratum
pmuxidelrho_jod = 0.99790
dpmuxidelrho_jod = 1.e-4*sqrt(4.6**2 + 7.5**2)
! calculated using one-sided limit
llim_pmuxidelrho_jod = pmuxidelrho_jod - sig_clg1*dpmuxidelrho_jod ! 0.99677

! Stoker et al., PRL
dpmuxidelrho_sto = 1.e-4*sqrt(16**2 + 6**2 + 16**2) ! stat, correction, syst
llim_pmuxidelrho_sto = 0.9955
! work backwards to get 0.9985 as a central value
pmuxidelrho_sto_o = llim_pmuxidelrho_sto + sig_clg1*dpmuxidelrho_sto
! but Stoker et al. had same error as Jodidio, so must be modified
pmuxidelrho_sto = pmuxidelrho_sto_o*pmuxidelrho_jod/pmuxidelrho_jod_o

! lim_pmuxidelrho = 0.99682 from Jodidio (erratum), including Stoker
! find weighted sum
var_jod = dpmuxidelrho_jod**2
var_sto = dpmuxidelrho_sto**2
var_js = 1./((1./var_jod) + (1./var_sto))
pmuxidelrho_js = var_js*(pmuxidelrho_jod/var_jod + pmuxidelrho_sto/var_sto)
dpmuxidelrho_js = sqrt(var_js)
! result in paper is 0.99682; this gives 0.996828 (rounding?)
lim_pmuxidelrho_js = pmuxidelrho_js - sig_clg1*dpmuxidelrho_js

! calculate pmuxi from this limit and twist rho, delta
!pmuxi_twist_js = pmuxidelrho_js*rho_twist/delta_twist
!dpmuxi_twist_js = pmuxi_twist_js* -
!           sqrt((dpmuxidelrho_js/pmuxidelrho_js)**2 -
!                + (drho_twist/rho_twist)**2 -
!                + (ddelta_twist/delta_twist)**2 )
!
!llim_pmuxi_twist_js = pmuxi_twist_js - sig_clg*dpmuxi_twist_js
!ulim_pmuxi_twist_js = pmuxi_twist_js + sig_clg*dpmuxi_twist_js
!
! calculate xi from Q_R^mu = 0.5*(1 + xi/3 -16*xi*delta/9) >=0
!ulim_xi_qr = 1./((16.*llim_delta_twist/9.) - (1./3.))



! put in numbers from Gaponenko explicitly, rather than calculating them
llim_pmuxi_twist_js = 0.9960
ulim_pmuxi_twist_js = 1.0040

! NSERC research proposal, 2004, for 2005 result
llim_pmuxi_twist_05 = 1.0 - sig_clg*1.e-4*sqrt(7.2**2 + 6.9**2)
ulim_pmuxi_twist_05 = 1.0 + sig_clg*1.e-4*sqrt(7.2**2 + 6.9**2)

! NSERC research proposal, 2004, for final result
llim_pmuxi_twist_f = 1.0 - sig_clg*1.e-4*sqrt(3.0**2 + 3.0**2)
ulim_pmuxi_twist_f = 1.0 + sig_clg*1.e-4*sqrt(3.0**2 + 3.0**2)

!delta_vs_xi_lh = `(0.1875*(3./xi + 1.))'
delta_vs_xi_lh = `(9./(16.*xi) + 3./16.)'
xi_vs_delta_lh = `(9./(16.*delta - 3.))'

if (exist(`delta')) then destroy delta
if (exist(`xi')) then destroy xi

xi = [scale_xlo:scale_xhi:0.001]
delta = eval(delta_vs_xi_lh)
x_lhline =xi
y_lhline = delta
destroy delta
destroy xi

! xi limits

xi = llim_pmuxi_bel
yhi_l_bel = eval(delta_vs_xi_lh)
xl_pmuxi_bel = [llim_pmuxi_bel;llim_pmuxi_bel]
yl_pmuxi_bel = [scale_ylo;yhi_l_bel]

xi = ulim_pmuxi_bel
yhi_h_bel = eval(delta_vs_xi_lh)
xu_pmuxi_bel = [ulim_pmuxi_bel;ulim_pmuxi_bel]
yu_pmuxi_bel = [scale_ylo;yhi_h_bel]

xi = llim_pmuxi_twist_js
yhi_l_twist_js = eval(delta_vs_xi_lh)
xl_pmuxi_twist_js = [llim_pmuxi_twist_js;llim_pmuxi_twist_js]
yl_pmuxi_twist_js = [scale_ylo;yhi_l_twist_js]

xi = ulim_pmuxi_twist_js
yhi_h_twist_js = eval(delta_vs_xi_lh)
xu_pmuxi_twist_js = [ulim_pmuxi_twist_js;ulim_pmuxi_twist_js]
yu_pmuxi_twist_js = [scale_ylo;yhi_h_twist_js]

xi = llim_pmuxi_twist_05
yhi_l_twist_05 = eval(delta_vs_xi_lh)
xl_pmuxi_twist_05 = [llim_pmuxi_twist_05;llim_pmuxi_twist_05]
yl_pmuxi_twist_05 = [scale_ylo;yhi_l_twist_05]

xi = ulim_pmuxi_twist_05
yhi_h_twist_05 = eval(delta_vs_xi_lh)
xu_pmuxi_twist_05 = [ulim_pmuxi_twist_05;ulim_pmuxi_twist_05]
yu_pmuxi_twist_05 = [scale_ylo;yhi_h_twist_05]

xi = llim_pmuxi_twist_f
yhi_l_twist_f = eval(delta_vs_xi_lh)
xl_pmuxi_twist_f = [llim_pmuxi_twist_f;llim_pmuxi_twist_f]
yl_pmuxi_twist_f = [scale_ylo;yhi_l_twist_f]

xi = ulim_pmuxi_twist_f
yhi_h_twist_f = eval(delta_vs_xi_lh)
xu_pmuxi_twist_f = [ulim_pmuxi_twist_f;ulim_pmuxi_twist_f]
yu_pmuxi_twist_f = [scale_ylo;yhi_h_twist_f]

! delta limits

delta = llim_delta_balke
xhi_l_balke = eval(xi_vs_delta_lh)
xl_delta_balke = [scale_xlo;xhi_l_balke]
yl_delta_balke = [llim_delta_balke;llim_delta_balke]

delta = ulim_delta_balke
xhi_h_balke = eval(xi_vs_delta_lh)
xu_delta_balke = [scale_xlo;xhi_h_balke]
yu_delta_balke = [ulim_delta_balke;ulim_delta_balke]

delta = llim_delta_twist
xhi_l_twist = eval(xi_vs_delta_lh)
xl_delta_twist = [scale_xlo;xhi_l_twist]
yl_delta_twist = [llim_delta_twist;llim_delta_twist]

delta = ulim_delta_twist
xhi_h_twist = eval(xi_vs_delta_lh)
xu_delta_twist = [scale_xlo;xhi_h_twist]
yu_delta_twist = [ulim_delta_twist;ulim_delta_twist]

delta = llim_delta_twist_05
xhi_l_twist_05 = eval(xi_vs_delta_lh)
xl_delta_twist_05 = [scale_xlo;xhi_l_twist_05]
yl_delta_twist_05 = [llim_delta_twist_05;llim_delta_twist_05]

delta = ulim_delta_twist_05
xhi_h_twist_05 = eval(xi_vs_delta_lh)
xu_delta_twist_05 = [scale_xlo;xhi_h_twist_05]
yu_delta_twist_05 = [ulim_delta_twist_05;ulim_delta_twist_05]

delta = llim_delta_twist_f
xhi_l_twist_f = eval(xi_vs_delta_lh)
xl_delta_twist_f = [scale_xlo;xhi_l_twist_f]
yl_delta_twist_f = [llim_delta_twist_f;llim_delta_twist_f]

delta = ulim_delta_twist_f
xhi_h_twist_f = eval(xi_vs_delta_lh)
xu_delta_twist_f = [scale_xlo;xhi_h_twist_f]
yu_delta_twist_f = [ulim_delta_twist_f;ulim_delta_twist_f]

! set up plot shape, size, etc.

set linthk -4
enable border
orientation portrait
window 21 0 20 100 100
window 21
set xtica 90
set ytica -90
if (~exist(`letter_size')) then letter_size = 3
set %txthit letter_size
set %xlabsz letter_size
set %ylabsz letter_size
set %xnumsz letter_size
set %ynumsz letter_size
set %xlaxis 25
set %ylaxis 25


pl_xlab = `<xi>'
pl_ylab = `<delta>'
! pl_titl = `TWIST limits on muon right-handed coupling'
pl_titl = `  '

set histyp 0

stat_flag = 0

label\xaxis pl_xlab
label\yaxis pl_ylab

clear

set linthk -6

! disallowed region (Q_R^<mu> < 0)

x_disallowed =  x_lhline// -
        [scale_xhi;scale_xhi]
y_disallowed =  y_lhline// - 
        [scale_yhi;scale_yhi]

! PDG allowed region
if (exist(`delta')) then destroy delta
if (exist(`xi')) then destroy xi
xi = [min(ulim_pmuxi_bel,xhi_l_balke):max(llim_pmuxi_bel,xhi_h_balke):-0.001]
delta = eval(delta_vs_xi_lh)
x_lhline =xi
y_lhline = delta
x_pdg = [llim_pmuxi_bel; -
             min(ulim_pmuxi_bel,xhi_l_balke); -
             min(ulim_pmuxi_bel,xhi_l_balke)]//xi//-
             [max(llim_pmuxi_bel,xhi_h_balke); -
             llim_pmuxi_bel; -
             llim_pmuxi_bel]
y_pdg = [llim_delta_balke; - 
             llim_delta_balke; -
             max(llim_delta_balke,yhi_h_bel)]//delta//- 
             [min(ulim_delta_balke,yhi_l_bel); -
             min(ulim_delta_balke,yhi_l_bel); -
             llim_delta_balke]
destroy delta
destroy xi

set %txthit 0.5*letter_size

color magenta
spacevec = [.8]
set hatch 2 spacevec 120
set lintyp 299
!set lintyp 102
gra\noaxes x_pdg y_pdg
color black
set lintyp 1
gra\noaxes x_pdg y_pdg
color magenta
set lintyp 4
gra\noaxes xl_pmuxi_bel yl_pmuxi_bel
gra\noaxes xu_pmuxi_bel yu_pmuxi_bel
gra\noaxes xl_delta_balke yl_delta_balke
gra\noaxes xu_delta_balke yu_delta_balke

set cursor -1
world\percent xl_delta_balke[1] yl_delta_balke[1] xtext ytext
set %xloc xtext+1
set %yloc ytext+1
text\noconfirm `PDG limits'


! TWIST 2004 allowed region
if (exist(`delta')) then destroy delta
if (exist(`xi')) then destroy xi
xi = [min(ulim_pmuxi_twist_js,xhi_l_twist):max(llim_pmuxi_twist_js,xhi_h_twist):-0.001]
delta = eval(delta_vs_xi_lh)
x_lhline =xi
y_lhline = delta
x_twist = [llim_pmuxi_twist_js; -
             min(ulim_pmuxi_twist_js,xhi_l_twist); -
             min(ulim_pmuxi_twist_js,xhi_l_twist)]//xi//-
             [max(llim_pmuxi_twist_js,xhi_h_twist); -
             llim_pmuxi_twist_js; -
             llim_pmuxi_twist_js]
y_twist = [llim_delta_twist; - 
             llim_delta_twist; -
             max(llim_delta_twist,yhi_h_twist_js)]//delta//- 
             [min(ulim_delta_twist,yhi_l_twist_js); -
             min(ulim_delta_twist,yhi_l_twist_js); -
             llim_delta_twist]
destroy delta
destroy xi

color blue
spacevec = [.4]
set hatch 3 spacevec 30
set lintyp 277
!set lintyp 103
gra\noaxes x_twist y_twist
set lintyp 1
color black
gra\noaxes x_twist y_twist
color blue
set lintyp 3
gra\noaxes xl_pmuxi_twist_js yl_pmuxi_twist_js
gra\noaxes xu_pmuxi_twist_js yu_pmuxi_twist_js
gra\noaxes xl_delta_twist yl_delta_twist
gra\noaxes xu_delta_twist yu_delta_twist

set cursor -1
world\percent xl_delta_twist[1] yl_delta_twist[1] xtext ytext
set %xloc xtext+1
set %yloc ytext+1
text\noconfirm `TWIST (current)'


! TWIST expected final allowed region
if (exist(`delta')) then destroy delta
if (exist(`xi')) then destroy xi
xi = [min(ulim_pmuxi_twist_f,xhi_l_twist_f):max(llim_pmuxi_twist_f,xhi_h_twist_f):-0.001]
delta = eval(delta_vs_xi_lh)
x_lhline =xi
y_lhline = delta
x_twist_f = [llim_pmuxi_twist_f; -
             min(ulim_pmuxi_twist_f,xhi_l_twist_f); -
             min(ulim_pmuxi_twist_f,xhi_l_twist_f)]//xi//-
             [max(llim_pmuxi_twist_f,xhi_h_twist_f); -
             llim_pmuxi_twist_f; -
             llim_pmuxi_twist_f]
y_twist_f = [llim_delta_twist_f; - 
             llim_delta_twist_f; -
             max(llim_delta_twist_f,yhi_h_twist_f)]//delta//- 
             [min(ulim_delta_twist_f,yhi_l_twist_f); -
             min(ulim_delta_twist_f,yhi_l_twist_f); -
             llim_delta_twist_f]
destroy delta
destroy xi

color white
spacevec = [.1]
set hatch 4 spacevec 60
set lintyp 255
!set lintyp 104
gra\noaxes x_twist_f y_twist_f
color white
set lintyp 1
gra\noaxes xl_pmuxi_twist_f yl_pmuxi_twist_f
gra\noaxes xu_pmuxi_twist_f yu_pmuxi_twist_f
gra\noaxes xl_delta_twist_f yl_delta_twist_f
gra\noaxes xu_delta_twist_f yu_delta_twist_f

set cursor -1
world\percent xl_delta_twist_f[1] yl_delta_twist_f[1] xtext ytext
set %xloc xtext+1
set %yloc ytext+1
text\noconfirm `TWIST (future)'

set %txthit letter_size

! left-handed coupling line
color white
set linthk -10
graph\axesonly x_xi y_delta

! axes and labels
set linthk -10
graph\noaxes x_lhline y_lhline

vector spacevec 1
spacevec = [.4]
set hatch 1 spacevec 45
color white
!set lintyp 101
!set lintyp 299
set lintyp 1
gra\noaxes x_disallowed y_disallowed

!color white
set cursor -1
set %xloc 60
set %yloc 70
text\noconfirm `Q<_>R<^><^><Z-1.4%><mu><_><lt>0'

color white
xi_sm = [1.0]
delta_sm = [0.75]
set pchar 12
set %charsz 1
gra\noax xi_sm delta_sm

return