! plot_twist_lrs_limit


set linthk -6

! use 68.27%CL (one sigma)
!cl = 0.6827

! 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 = 200
scale_xhi = 1200
scale_ylo = -0.1
scale_yhi = 0.1
scale scale_xlo scale_xhi 5 scale_ylo scale_yhi 4
x_mass = [scale_xlo;scale_xhi]
y_zeta = [scale_ylo;scale_yhi]


!
! rho
!

! Derenzo PR 181 (1969) 1954 
rho_der = 0.7518
drho_der = 1.e-4*26.
llim_rho_der = rho_der-sig_clg*drho_der
ulim_rho_der = rho_der+sig_clg*drho_der
! lrs model, limit on zeta from Derenzo rho; rho>=0.75, so use one-sided limit
llim_rho_der = 0.75-sig_clg1*drho_der ! note this doesn't use rho_der
zeta_rho_der_1[1] = -sqrt((0.75 - llim_rho_der)*2./3.)
zeta_rho_der_1[2] = zeta_rho_der_1[1]
zeta_rho_der_2[1] = -zeta_rho_der_1[1]
zeta_rho_der_2[2] = zeta_rho_der_2[1]

! musser et al. (2005) 0.7580 pm 0.00032 pm 0.00097 pm 0.00023
rho_twist_muss = 0.75080
drho_twist_muss = 1.e-4*sqrt(3.2**2 + 9.7**2 + 2.3**2)

! macdonald et al. (2008) 0.75014 pm 0.00017 pm 0.00044 pm 0.00011
rho_twist_macd = 0.75014
drho_twist_macd = 1.e-4*sqrt(1.7**2 + 4.6**2 + 1.1**2)
llim_rho_twist_macd = rho_twist_macd-sig_clg*drho_twist_macd
ulim_rho_twist_macd = rho_twist_macd+sig_clg*drho_twist_macd
! lrs model, limit on zeta from TWIST rho; rho>=0.75, so use one-sided limit
llim_rho_twist_macd = 0.75-sig_clg1*drho_twist_macd ! note this doesn't use rho_twist
zeta_rho_twist_macd_1[1] = -sqrt((0.75 - llim_rho_twist_macd)*2./3.)
zeta_rho_twist_macd_1[2] = zeta_rho_twist_macd_1[1]
zeta_rho_twist_macd_2[1] = -zeta_rho_twist_macd_1[1]
zeta_rho_twist_macd_2[2] = zeta_rho_twist_macd_2[1]

! 2008 est, for final results
rho_twist_08fp = 0.7500
drho_twist_08fp = 1.e-4*sqrt(1.3**2 + 2.4**2 + 1.1**2)
llim_rho_twist_08fp = rho_twist_08fp-sig_clg*drho_twist_08fp
ulim_rho_twist_08fp = rho_twist_08fp+sig_clg*drho_twist_08fp
! lrs model, limit on zeta from final TWIST rho
llim_rho_twist_08fp = 0.75-sig_clg1*drho_twist_08fp
zeta_rho_twist_1_08fp[1] = -sqrt((0.75 - llim_rho_twist_08fp)*2./3.)
zeta_rho_twist_1_08fp[2] = zeta_rho_twist_1_08fp[1]
zeta_rho_twist_2_08fp[1] = -zeta_rho_twist_1_08fp[1]
zeta_rho_twist_2_08fp[2] = zeta_rho_twist_2_08fp[1]

! set up plot shape, size, etc.

enable border
orientation portrait
window 21 0 20 100 100
!window 21 0 0 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 = `(g<_>L<^>/g<_>R<^>)m<_>2<^> (GeV/c<^>2<_>)'
pl_ylab = `(g<_>R<^>/g<_>L<^>)<zeta>'
pl_titl = `General LRS model'
! pl_titl = `  '

set histyp 0

label\xaxis pl_xlab
label\yaxis pl_ylab

clear

legend frame 30 77 85 87
legend frame off

!!!
plot_rho_derenzo:

legend_text = -
`Derenzo et al. <rho>              '

!set lintyp 1
!color black
!gra\noaxes x_mass zeta_rho_der_1
!gra\noaxes x_mass zeta_rho_der_2

set lintyp 4
color white
legend on
gra\noaxes legend_text x_mass zeta_rho_der_1
legend off
gra\noaxes x_mass zeta_rho_der_2

!!!
plot_twist:

!!!
plot_rho_twist:

legend_text = -
     `TWIST <rho>                  '

set lintyp 5
color blue
legend on
gra\noaxes legend_text x_mass zeta_rho_twist_macd_1
legend off
!gra\noaxes x_mass zeta_rho_twist_macd_1
gra\noaxes x_mass zeta_rho_twist_macd_2

!!!
plot_08fp:

! plot TWIST expected final values as of 2008

legend_text = -
     `TWIST <rho> estimated final    '
set lintyp 1
color red
legend on
gra\noaxes legend_text x_mass zeta_rho_twist_1_08fp
legend off
gra\noaxes x_mass zeta_rho_twist_2_08fp

window\noconfirm 99 70 64 95 67
erasewindow
window 21

color white
set %xloc 90
set %yloc 57
set linthk -10
set cursor -11
text\noconfirm `allowed'

graph\axesonly x_mass zeta_rho_twist_2_08fp