#set term pdf fontscale 0.75 size 6in,4in
#set output 'kickpulse.pdf'
set xrange [0:500]
set yrange [-50:300]
set xlabel 'time [ns]'
set ylabel 'B_{kicker} [G]'
set term x11 0
offset = -25000+180

f=1.42
rwidth=157.
rcenter=120
length_kick=3*1.27
p=3.095
beta=7.54
eta=8.
xinf=0.077-0.0
Aperture= 0.043-0.0
k_min = (xinf - Aperture)/beta
B_to_kick = length_kick *0.3/p/1.e4*beta  #[G to T]
Bmin = k_min*beta/B_to_kick
print 'k_min = ',k_min,' Bmin=', Bmin
# set object 1 rectangle center rcenter,Bmin/2 size 155,Bmin
set arrow 1 from rcenter-rwidth/2.,0 to rcenter-rwidth/2,140 nohead dt 2
set arrow 2 from rcenter+rwidth/2.,0 to rcenter+rwidth/2,140 nohead dt 2
set arrow 3 from rcenter-rwidth/2.-20,Bmin to rcenter+rwidth/2.+20, Bmin nohead dt 2 
plot 'traces-Ir.txt' u 1:($2*f) w l t 'B_{kicker} [G]','../ps12543_mod/pulseshape_6_12543_mod.dat' u ($1+offset):($2/150) w l t 'muons turn 1', '' u ($1+offset+149):($2/150) w l lt 2 t 'muons turn 2', '' u ($1+offset+2*149):($2/150) w l lt 2 t 'muons turn 3'
unset arrow 1
unset arrow 2
unset arrow 3

pause -1
reset
set autoscale y
set xrange [0:300]
set xlabel 'time [ns]'
set term x11 1
delta_half = Aperture/eta/2.
sigma2_delta = delta_half**2/3.
set yrange [-0.6:0.6]
set ytics -0.6,0.2
plot 'traces-Ir.txt' u 1:(($2*f > Bmin) ? (xinf-$2*f*B_to_kick)/eta*100/2:1/0):($3-$2+delta_half)*100  w  errorlines pt 0 t '<{/Symbol d}> [%]'
#plot 'traces-Ir.txt' u 1:(xinf-$2*f*B_to_kick/eta*100):($3-$2+delta_half)*100  w l e t '<{/Symbol d}> [%]'
pause -1

set term x11 2
reset
set xrange [0:500]
set yrange [-50:300]
set xlabel 'time [ns]'
set ylabel 'B_{kicker} [G]'
f=1.1
rwidth=131.
rcenter=118
# set object 1 rectangle rcenter rcenter,Bmin/2 size rwidth,Bmin
set arrow 1 from rcenter-rwidth/2.,0 to rcenter-rwidth/2,140 nohead dt 2
set arrow 2 from rcenter+rwidth/2.,0 to rcenter+rwidth/2,140 nohead dt 2
set arrow 3 from rcenter-rwidth/2.-20,Bmin to rcenter+rwidth/2.+20, Bmin nohead dt 2 

plot 'traces-Ir.txt' u 1:($2*f) w l t 'B_{kicker} [G]','../ps12543_mod/pulseshape_6_12543_mod.dat' u ($1-25000+180):($2/150) w l t 'muons turn 1', '' u ($1-25000+180+149):($2/150) w l lt 2 t 'muons turn 2', '' u ($1-25000+180+2*149):($2/150) w l lt 2 t 'muons turn 3'
 unset arrow 1
 unset arrow 2
 unset arrow 3

pause -1
set term x11 3
reset
stats 'fort.13' u 2:3
f = Bmax/STATS_max_y 
stats 'fort.13' u 2:4
g= 1./STATS_max_y

print ' stats "fort.13" u 2:(($3*f > Bmin) ? $4*g:1/0)'
stats 'fort.13' u 2:(($3*f > Bmin) ? $4*g:1/0)
ave_I= STATS_mean_y
sum_I= STATS_sum_y
print 'ave_I =',ave_I,' max_y = ', STATS_max_y, 'sum_I = ', STATS_sum_y
print ' stats "fort.13" u 2:(($3*f > Bmin) ? (xinf-$3*f*B_to_kick)/eta*100/2*$4*g:1/0)'
stats 'fort.13' u 2:(($3*f > Bmin) ? (xinf-$3*f*B_to_kick)/eta*100/2*$4*g:1/0)
print '1/sum_I * sum(delta*I) = ', STATS_sum_y/sum_I,' <r> = eta<delta> = ', STATS_sum_y/sum_I * eta /100.
set xrange [0:300]
set yrange [0:250]
set ytics 0,50 nomirror
set y2range [0:1.1]
set y2tics 0,0.1
set y2label 'pulse height'
set ylabel 'B_{kicker} [G]'
set xlabel 'time [ns]'
plot 'fort.13' u 2:($3*f) w l t 'kick [G]', '' u 2:($4*g) axes x1y2 w l t 'pulse'

reset
set term x11 6
print ' stats "fort.13" u 2:(($3*f > Bmin) ? (((xinf-$3*f*B_to_kick)/eta/2)**2*$4*g):1/0)'
stats 'fort.13' u 2:(($3*f > Bmin) ? (((xinf-$3*f*B_to_kick)/eta/2)**2*$4*g):1/0)
print '1/sum_I * sum(delta**2*I) = ', STATS_sum_y/sum_I,' sqrt{<r**2>} = sqrt{(eta<delta>)**2} = ', sqrt(STATS_sum_y/sum_I) * (eta)
set xrange [0:300]
set yrange [0:250]
set ytics 0,50 nomirror
set y2range [0:1.1]
set y2tics 0,0.1
set y2label 'pulse height'
set ylabel 'B_{kicker} [G]'
set xlabel 'time [ns]'
plot 'fort.13' u 2:(($3*f > Bmin) ? (sqrt((xinf-$3*f*B_to_kick)/eta*100/2)**2*$4*g):1/0) w l t 'sigma_delta'


set term x11 4
reset
set yrange [-0.6:0.6]
set ytics -0.6,.1
set xrange [0:300]
set xlabel 'time [ns]'
set ylabel 'momentum offset [%]' 
set zeroaxis
set grid
plot 'fort.13' u 2:(($3*f > Bmin) ? (xinf-$3*f*B_to_kick)/eta*100/2:1/0):($3-$3+delta_half)*100  w  errorlines pt 0 t '<{/Symbol d}> [%]'
plot 'fort.13' u 2:(($3*f > Bmin) ? (xinf-$3*f*B_to_kick)/eta*100/2:1/0):($3-$3+delta_half)*100  w  errorlines pt 0 t '<{/Symbol d}> [%]'
pause -1
set term x11 5
splot 'fort.13' u 2:(($3*f > Bmin) ? (xinf-$3*f*B_to_kick)/eta*100/2:1/0):4

reset
set key off
set view map
set pm3d at s hidden3d 100
set surface
set dgrid3d 50, 50, 3
set view map scale 1
splot 'fort.13' u 2:(($3*f > Bmin) ? (xinf-$3*f*B_to_kick)/eta*100/2:1/0):4 w pm3d