! rhmu_lim_extr.pcm

defaults
clear
destroy *

@decay_parameter_results.pcm

!scale_xlo = 0.98
!scale_xhi = 1.03
!scale_ylo = 0.735
!scale_yhi = 0.760
scale_xlo = 0.98
scale_xhi = 1.02
scale_ylo = 0.740
scale_yhi = 0.755
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_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 = evaluate(delta_vs_xi_lh)
x_lhline =xi
y_lhline = delta
destroy delta
destroy xi


! xi limits

xi = llim_pmuxi_bel
yhi_l_bel = evaluate(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 = evaluate(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_jami
yhi_l_twist_jami = evaluate(delta_vs_xi_lh)
xl_pmuxi_twist_jami = [llim_pmuxi_twist_jami;llim_pmuxi_twist_jami]
yl_pmuxi_twist_jami = [scale_ylo;yhi_l_twist_jami]

xi = ulim_pmuxi_twist_jami
yhi_h_twist_jami = evaluate(delta_vs_xi_lh)
xu_pmuxi_twist_jami = [ulim_pmuxi_twist_jami;ulim_pmuxi_twist_jami]
yu_pmuxi_twist_jami = [scale_ylo;yhi_h_twist_jami]

xi = llim_pmuxi_twist_08fp
yhi_l_twist_08fp = evaluate(delta_vs_xi_lh)
xl_pmuxi_twist_08fp = [llim_pmuxi_twist_08fp;llim_pmuxi_twist_08fp]
yl_pmuxi_twist_08fp = [scale_ylo;yhi_l_twist_08fp]

xi = ulim_pmuxi_twist_08fp
yhi_h_twist_08fp = evaluate(delta_vs_xi_lh)
xu_pmuxi_twist_08fp = [ulim_pmuxi_twist_08fp;ulim_pmuxi_twist_08fp]
yu_pmuxi_twist_08fp = [scale_ylo;yhi_h_twist_08fp]

! delta limits

delta = llim_delta_bal
xhi_l_bal = evaluate(xi_vs_delta_lh)
xl_delta_bal = [scale_xlo;xhi_l_bal]
yl_delta_bal = [llim_delta_bal;llim_delta_bal]

delta = ulim_delta_bal
xhi_h_bal = evaluate(xi_vs_delta_lh)
xu_delta_bal = [scale_xlo;xhi_h_bal]
yu_delta_bal = [ulim_delta_bal;ulim_delta_bal]

delta = llim_delta_twist_macd
xhi_l_twist_macd = evaluate(xi_vs_delta_lh)
xl_delta_twist_macd = [scale_xlo;xhi_l_twist_macd]
yl_delta_twist_macd = [llim_delta_twist_macd;llim_delta_twist_macd]

delta = ulim_delta_twist_macd
xhi_h_twist_macd = evaluate(xi_vs_delta_lh)
xu_delta_twist_macd = [scale_xlo;xhi_h_twist_macd]
yu_delta_twist_macd = [ulim_delta_twist_macd;ulim_delta_twist_macd]

delta = llim_delta_twist_08fp
xhi_l_twist_08fp = evaluate(xi_vs_delta_lh)
xl_delta_twist_08fp = [scale_xlo;xhi_l_twist_08fp]
yl_delta_twist_08fp = [llim_delta_twist_08fp;llim_delta_twist_08fp]

delta = ulim_delta_twist_08fp
xhi_h_twist_08fp = evaluate(xi_vs_delta_lh)
xu_delta_twist_08fp = [scale_xlo;xhi_h_twist_08fp]
yu_delta_twist_08fp = [ulim_delta_twist_08fp;ulim_delta_twist_08fp]

! setup for plots

set linewidth 6
set aspectratio 1  !landscape = 8.5/11

scale scale_xlo scale_xhi 4 scale_ylo scale_yhi 3

!window 21 0 20 100 100  ! not necessary if aspectratio is 1
!window 21
set xtica 90
set ytica -90
if (~exist(`letter_size')) then letter_size = 3
set %txthit letter_size
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 xlabel pl_xlab
!set ylabel pl_ylab
set xlabel ' '
set ylabel ' '

! 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_bal):max(llim_pmuxi_bel,xhi_h_bal):-0.001]
delta = evaluate(delta_vs_xi_lh)
x_lhline = xi
y_lhline = delta

x_pdg = [llim_pmuxi_bel; -
             min(ulim_pmuxi_bel,xhi_l_bal); -
             min(ulim_pmuxi_bel,xhi_l_bal)]//xi//-
             [max(llim_pmuxi_bel,xhi_h_bal); -
             llim_pmuxi_bel; -
             llim_pmuxi_bel]
y_pdg = [llim_delta_bal; - 
             llim_delta_bal; -
             max(llim_delta_bal,yhi_h_bel)]//delta//- 
             [min(ulim_delta_bal,yhi_l_bel); -
             min(ulim_delta_bal,yhi_l_bel); -
             llim_delta_bal]
destroy delta
destroy xi

set %txthit 0.7*letter_size

set color white
set areafillcolor silver
gra\overlay x_pdg y_pdg
set areafillcolor 0
set color silver
set curvelinetype 4
gra\overlay xl_pmuxi_bel yl_pmuxi_bel
gra\overlay xu_pmuxi_bel yu_pmuxi_bel
gra\overlay xl_delta_bal yl_delta_bal
gra\overlay xu_delta_bal yu_delta_bal

set color black
set textinteractive 0
world\percent xl_delta_bal[1] yl_delta_bal[1] xtext ytext
set %xtextlocation xtext+1
set %ytextlocation ytext+1
text\graphunits 'Pre-TWIST limits'


! Jamieson-MacDonald limits
if (exist(`delta')) then destroy delta
if (exist(`xi')) then destroy xi
xi = [min(ulim_pmuxi_twist_jami,xhi_l_twist_macd):max(llim_pmuxi_twist_jami,xhi_h_twist_macd):-0.001]
delta = evaluate(delta_vs_xi_lh)
x_lhline = xi
y_lhline = delta

x_twist = [llim_pmuxi_twist_jami; -
             min(ulim_pmuxi_twist_jami,xhi_l_twist_macd); -
             min(ulim_pmuxi_twist_jami,xhi_l_twist_macd)]//xi//-
             [max(llim_pmuxi_twist_jami,xhi_h_twist_macd); -
             llim_pmuxi_twist_jami; -
             llim_pmuxi_twist_jami]
y_twist = [llim_delta_twist_macd; - 
             llim_delta_twist_macd; -
             max(llim_delta_twist_macd,yhi_h_twist_jami)]//delta//- 
             [min(ulim_delta_twist_macd,yhi_l_twist_jami); -
             min(ulim_delta_twist_macd,yhi_l_twist_jami); -
             llim_delta_twist_macd]
destroy delta
destroy xi

set %txthit 0.7*letter_size

set color white
set areafillcolor blue
gra\overlay x_twist y_twist
set areafillcolor 0
set color blue
set curvelinetype 4
gra\overlay xl_pmuxi_twist_jami yl_pmuxi_twist_jami
gra\overlay xu_pmuxi_twist_jami yu_pmuxi_twist_jami
gra\overlay xl_delta_twist_macd yl_delta_twist_macd
gra\overlay xu_delta_twist_macd yu_delta_twist_macd

set color blue
set textinteractive 0
world\percent xl_delta_twist_macd[1] yl_delta_twist_macd[1] xtext ytext
set %xtextlocation xtext+1
set %ytextlocation ytext-4
text\graphunits 'current TWIST limits'


! Future expected TWIST limits as of 2008
if (exist(`delta')) then destroy delta
if (exist(`xi')) then destroy xi
xi = [min(ulim_pmuxi_twist_08fp,xhi_l_twist_08fp):max(llim_pmuxi_twist_08fp,xhi_h_twist_08fp):-0.001]
delta = evaluate(delta_vs_xi_lh)
x_lhline = xi
y_lhline = delta

x_twist = [llim_pmuxi_twist_08fp; -
             min(ulim_pmuxi_twist_08fp,xhi_l_twist_08fp); -
             min(ulim_pmuxi_twist_08fp,xhi_l_twist_08fp)]//xi//-
             [max(llim_pmuxi_twist_08fp,xhi_h_twist_08fp); -
             llim_pmuxi_twist_08fp; -
             llim_pmuxi_twist_08fp]
y_twist = [llim_delta_twist_08fp; - 
             llim_delta_twist_08fp; -
             max(llim_delta_twist_08fp,yhi_h_twist_08fp)]//delta//- 
             [min(ulim_delta_twist_08fp,yhi_l_twist_08fp); -
             min(ulim_delta_twist_08fp,yhi_l_twist_08fp); -
             llim_delta_twist_08fp]
destroy delta
destroy xi

set %txthit 0.7*letter_size

set color white
set areafillcolor red
gra\overlay x_twist y_twist
set areafillcolor 0
set color red
set curvelinetype 4
gra\overlay xl_pmuxi_twist_08fp yl_pmuxi_twist_08fp
gra\overlay xu_pmuxi_twist_08fp yu_pmuxi_twist_08fp
gra\overlay xl_delta_twist_08fp yl_delta_twist_08fp
gra\overlay xu_delta_twist_08fp yu_delta_twist_08fp

set color red
set textinteractive 0
world\percent xl_delta_twist_08fp[1] yl_delta_twist_08fp[1] xtext ytext
set %xtextlocation xtext+35
set %ytextlocation ytext+1
text\graphunits 'future TWIST limits'




!=======================
set %txthit letter_size

! axes and labels
set color black
set xaxislinewidth 3
set yaxislinewidth 3
graph\axesonly x_xi y_delta

! left-handed coupling line
if (exist(`delta')) then destroy delta
if (exist(`xi')) then destroy xi
xi = [scale_xlo:scale_xhi:0.001]
delta = evaluate(delta_vs_xi_lh)
x_lhline = xi
y_lhline = delta
destroy delta
destroy xi
set curvelinetype 1
set linewidth 2
graph\overlay x_lhline y_lhline