#set size 1.4,1.0
set out 'mylrs_limits_pmuxi.eps'
set term postscript enhanced eps color "Times-Roman" 16
#set grid y2 lt 18
#set grid x2 lt 18
set grid x lt 18
set grid y lt 18
set key spacing 1.2
set title 'Manifest and Non-Manifest LRS Model Limits from P_{/Symbol m}{/Symbol x}'
set ylabel '{/Symbol z}'
set xlabel 'm_2 (GeV/c^2)'

set samples 500000
set border 3
set xtics nomirror
set ytics nomirror
mw = 80.390
dpmuxi1 = 0.0073
# positive and negative zeta, for MLRS case
fp1(x) = 0.5*(-(mw/x)**2+sqrt( dpmuxi1-3.*(mw/x)**4 ) ) 
fm1(x) = 0.5*(-(mw/x)**2-sqrt( dpmuxi1-3.*(mw/x)**4 ) )
# positive and negative zeta, for non-MLRS case  
gp1(x) = ((mw/x)**2>dpmuxi1 / 2. / (2.03)**2 ) ? sqrt( dpmuxi1/2.-(mw/x)**4 ) : 0.5*(1.03*(mw/x)**2+sqrt( dpmuxi1-(2.+1.03**2)*(mw/x)**4 ) )
gm1(x) = ((mw/x)**2>dpmuxi1 / 2. / (2.03)**2 ) ? -sqrt( dpmuxi1/2.-(mw/x)**4 ) : -0.5*(1.03*(mw/x)**2+sqrt( dpmuxi1-(2.+1.03**2)*(mw/x)**4 ) )
dpmuxi2 = 0.0035
fp2(x) = 0.5*(-(mw/x)**2+sqrt( dpmuxi2-3.*(mw/x)**4 ) ) 
fm2(x) = 0.5*(-(mw/x)**2-sqrt( dpmuxi2-3.*(mw/x)**4 ) ) 
gp2(x) = ((mw/x)**2>dpmuxi2 / 2. / (2.03)**2 ) ? sqrt( dpmuxi2/2.-(mw/x)**4 ) : 0.5*(1.03*(mw/x)**2+sqrt( dpmuxi2-(2.+1.03**2)*(mw/x)**4 ) )
gm2(x) = ((mw/x)**2>dpmuxi2 / 2. / (2.03)**2 ) ? -sqrt( dpmuxi2/2.-(mw/x)**4 ) : -0.5*(1.03*(mw/x)**2+sqrt( dpmuxi2-(2.+1.03**2)*(mw/x)**4 ) )
#dpmuxi3 = 0.0009
#fp3(x) = 0.5*(-(mw/x)**2+sqrt( dpmuxi3-3.*(mw/x)**4 ) ) 
#fm3(x) = 0.5*(-(mw/x)**2-sqrt( dpmuxi3-3.*(mw/x)**4 ) )
#gp3(x) = ((mw/x)**2>dpmuxi3 / 2. / (2.03)**2 ) ? sqrt( dpmuxi3/2.-(mw/x)**4 ) : 0.5*(1.03*(mw/x)**2+sqrt( dpmuxi3-(2.+1.03**2)*(mw/x)**4 ) )
#gm3(x) = ((mw/x)**2>dpmuxi3 / 2. / (2.03)**2 ) ? -sqrt( dpmuxi3/2.-(mw/x)**4 ) : -0.5*(1.03*(mw/x)**2+sqrt( dpmuxi3-(2.+1.03**2)*(mw/x)**4 ) )
#dpmuxi4 = -0.0015
#fp4(x) = 0.5*(-(mw/x)**2+sqrt( dpmuxi4-3.*(mw/x)**4 ) ) 
#fm4(x) = 0.5*(-(mw/x)**2-sqrt( dpmuxi4-3.*(mw/x)**4 ) )
#gp4(x) = ((mw/x)**2>dpmuxi4 / 2. / (2.03)**2 ) ? sqrt( dpmuxi4/2.-(mw/x)**4 ) : 0.5*(1.03*(mw/x)**2+sqrt( dpmuxi4-(2.+1.03**2)*(mw/x)**4 ) )
#gm4(x) = ((mw/x)**2>dpmuxi4 / 2. / (2.03)**2 ) ? -sqrt( dpmuxi4/2.-(mw/x)**4 ) : -0.5*(1.03*(mw/x)**2+sqrt( dpmuxi4-(2.+1.03**2)*(mw/x)**4 ) )
#dpmuxi5 = -0.0053
#fp5(x) = 0.5*(-(mw/x)**2+sqrt( dpmuxi5-3.*(mw/x)**4 ) ) 
#fm5(x) = 0.5*(-(mw/x)**2-sqrt( dpmuxi5-3.*(mw/x)**4 ) )
#gp5(x) = ((mw/x)**2>dpmuxi5 / 2. / (2.03)**2 ) ? sqrt( dpmuxi5/2.-(mw/x)**4 ) : 0.5*(1.03*(mw/x)**2+sqrt( dpmuxi5-(2.+1.03**2)*(mw/x)**4 ) )
#gm5(x) = ((mw/x)**2>dpmuxi5 / 2. / (2.03)**2 ) ? -sqrt( dpmuxi5/2.-(mw/x)**4 ) : -0.5*(1.03*(mw/x)**2+sqrt( dpmuxi5-(2.+1.03**2)*(mw/x)**4 ) )

set yrange [-0.08:0.1]
set xtics 300,100.
set ytics -0.1,0.02
#set y2tics (-0.052,-0.037,-0.018,-0.010,0.006,0.013,0.029,0.042)
#set x2tics (352.,419.,594.,805.)
set label 'Allowed' at 1000.,-0.005
#plot [300.:1200.] gp1(x) lt 2, gm1(x) lt 2, fp1(x) lt 3, fm1(x) lt 3
plot [300.:1200.] fp1(x) lt 16 lw 2 title 'MLRS 2{/Symbol s} 1-P_{/Symbol m}{/Symbol x} = 0.0085',fm1(x) lt 16 lw 2 notitle, gp1(x) lt -1  title 'non-MLRS 2{/Symbol s} 1-P_{/Symbol m}{/Symbol x} = 0.0085', gm1(x) lt -1 notitle,\
   fp2(x) lt 4 lw 2 title 'MLRS 1{/Symbol s} 1-P_{/Symbol m}{/Symbol x} = 0.0035',fm2(x) lt 4 lw 2 notitle, gp2(x) lt 3  title 'non-MLRS 1{/Symbol s} 1-P_{/Symbol m}{/Symbol x} = 0.0035', gm2(x) lt 3 notitle
#                fp3(x) lt 3 lw 2 title '1-P_{/Symbol m}{/Symbol x} = 0.0073',fm3(x) lt 3 lw 2 notitle, gp3(x) lt 3 notitle, gm3(x) lt 3 notitle
# 
#                fp4(x) lt 4 lw 2 title 'TWIST goal {/Symbol D}P_{/Symbol m}{/Symbol x} = 0.0003',fm4(x) lt 4 lw 2 notitle, gp4(x) lt 4 notitle, gm4(x) lt 4 notitle

set out 'mylrs_limits_pmuxi.ps'
set term postscript enhanced color "Times-Roman" 16
replot