set term pdf fontscale 0.75 size 6in,4in
outfile = 'BetaCrossE.pdf'
set output outfile
set timestamp offset graph 1.04,0.7 rotate font 'Verdana,6'
set label 5 at graph 0.06,0.05 dir noenhance font ',10'
set label 3 'plotting_scripts/BetaCrossE.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
# s-px phase space
set xlabel 'x [mm]'
set ylabel 'p_x [mrad]'
file = dir.'/all_EndOfTracking_phase_space.dat'
print 'file = ', file
print 'x-px phase space next'
stats file u ($4*1000):($10 > time_off ? ($5*1000):1/0)
rms_x=sprintf("%3.3f",STATS_stddev_x)
rms_xp=sprintf("%3.3f",STATS_stddev_y)
set label 1 '{/Symbol s}_x = '.rms_x.' mm' at graph 0.05,0.95
set label 2 '{/Symbol s}_{/Symbol q} = '.rms_xp.' mrad' at graph 0.7,0.95
print 'rms_x = ', rms_x
#plot file u ($4*1000):($10 > time_off ? ($5*1000):1/0) pt 5 ps .25 not
#unset timestamp

# stats average energy and std dev energy
max=0.5 #max value microseconds
min=-0.5
width=(max-min)/n #interval width
#function used to map a value to the intervals
hist(x,width)=width*floor(x/width)+width/2.0
set xrange [min:max]
set xtics min,(max-min)/10,max
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
set tics out nomirror
set label 5 at graph 0.06,0.95 dir noenhance font ',8'
set xlabel  "{/Symbol D}p/p [%]"
set ylabel "Frequency"
#set ytics 50
#count and plot
stats file u 9:($10 > time_off ? ($9):1/0)
avee =sprintf("%5.4f",STATS_mean_y*100)
set label 1  '<{/Symbol D}p/p>= '.avee.'%' at graph 0.65,0.9
stddev =sprintf("%4.3f",STATS_stddev_y*100)
set label  2 '{/Symbol s}_p/p= '.stddev.'%' at graph 0.65,0.8
#set title 'Energy at End of Tracking' noenhance
print '<dp/p>= ',avee
#histogram energy
set timestamp offset graph 1.04,0.7 rotate font 'Verdana,6'
#plot file u (hist(($9*100),width)):($10 > time_off ? (1.0):1/0) smooth freq w boxes lc rgb"green" not
plot "< awk '{if ($10 > time_off) {print $0}}' ".file u (hist(($9*100),width)):(1.0) smooth freq w boxes lc rgb"green" not 

# C_e
f = 2/3.e8/1.45*1.e6
betax = 7.8
gamma = 1./betax
eta = 8.8
etap=0.

# stats average C_e
reset
stats file u (($4-$9*eta)**2*gamma + ($5-$9*etap)**2*betax)*1.e6:($10 > time_off ?($17/$19*$9)*f:1/0)
avece =sprintf("%3.3f",STATS_mean_y)
rmsce =sprintf("%3.3f",STATS_stddev_y)
set label 1 '<C_e>= '.avece.' ppm' at graph 0.6,0.9
set label 2 '{/Symbol s}_e= '.rmsce.' ppm' at graph 0.6,0.85
print '<C_e> =',avece
#unset label 2

#C_e histogram
set timestamp offset graph 1.07,  0.7 rotate font 'Verdana,6'
set label 3 'plotting_scripts/BetaCrossE.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
set label 5 at graph 0.55,0.05 dir noenhance font ',8'
set xrange [0:4]
set xtics 0,0.5
set xlabel 'C_e [ppm]'
max = 4 #max value microseconds
min=0
width=(max-min)/(1.0*n) #interval width
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
plot "< awk '{if ($10 > time_off) {print $0}}' ".file u (hist(($17/$19*$9*f),width)):(1.0) smooth freq w boxes lc rgb"green" not

set label 5 at graph 0.55,0.1 dir noenhance font ',8'
set xrange [0:0.5]
set xtics 0.,0.1
max=0.5
width=(max-min)/(1.0*n)
set boxwidth width*0.9
plot "< awk '{if ($10 > time_off) {print $0}}' ".file u (hist(($17/$19*$9*f),width)):(1.0) smooth freq w boxes lc rgb"green" not

#C_e vs amplitude
unset xtics
set xtics auto
set autoscale
set label 5 at graph 0.06,0.05 dir noenhance font ',10'
set xlabel "amplitude [mm-mrad]"
set ylabel "C_e [ppm]"
set label 5 at graph 0.06,0.05 dir noenhance font ',10'
plot file u (($4-$9*eta)**2*gamma + ($5-$9*etap)**2*betax)*1.e6:($10 > time_off ? ($17/$19*$9)*f:1/0) pt 5 ps .25 not

#C_e vs dp/p
set xlabel '{/Symbol D}p/p [%]'
plot file u ($9*100):($10 > time_off ? ($17/$19*$9)*f:1/0) pt 5 ps .25 not

# vertical phase space
stats file u ($6*1000):($10 > time_off ?($7*1000):1/0)
rms_y=sprintf("%3.3f",STATS_stddev_x)
rms_yp=sprintf("%3.3f",STATS_stddev_y)
set xlabel 'y [mm]'
set ylabel 'p_y [mrad]'
set label 1 '{/Symbol s}_y = '.rms_y.' mm' at graph 0.05,0.95
set label 2 '{/Symbol s}_{/Symbol q} = '.rms_yp.' mrad' at graph 0.7,0.95
#plot file u ($6*1000):($10 > time_off ? ($7*1000):1/0) pt 5 ps .25 not
print 'rms_y =', rms_y

# C_p vs amplitude
reset
unset xtics
set xtics auto
set autoscale
set xlabel "amplitude [mm-mrad]"
set ylabel "C_p [ppm]"
betay = 19.5
gamma = 1./betay
eta = 8.8
etap=0.
f = 1/1.45*1.e6/2.
p0c  = 3.094353005E9
mass=105.6583715E6

beta(x) = (x+1)*p0c/((x+1)**2*p0c**2+mass**2)**0.5
beta2(x) = ((x+1)*p0c)**2/((x+1)**2*p0c**2+mass**2)
ga(x) = (p0c**2*(x+1)**2 + mass**2)**.5/mass
bet(x) = (1.-1./ga(x)**2)**.5

# stats average C_p
stats file u (($6)**2*gamma + ($7)**2*betay)*1.e6:($10 > time_off ?($21/beta2($9)/$19)*f:1/0)
avecp =sprintf("%3.3f",STATS_mean_y)
set label 1 '<C_p>= '.avecp.' ppm' at graph 0.6,0.9
unset label 2
set xlabel 'vertical amplitude [mm-mrad]'
print '<C_p> = ',avecp
#C_p vs amplitude
set label 1 '<C_p>= '.avecp.' ppm' at graph 0.6,0.1
set label 5 at graph 0.55,0.05 dir noenhance font ',10'
plot file u (($6)**2*gamma + ($7)**2*betay)*1.e6:($10 > time_off ? ($21/beta2($9)/$19)*f:1/0) pt 5 ps .25 not

#C_p vs y'^2
set xlabel '{/Symbol q}^2 [mm^2]
plot file u ($7**2)*1.e6:($10 > time_off ? ($21/beta2($9)/$19)*f:1/0) pt 5 ps .25 not


#C_p histogram
set label 1 '<C_p>= '.avecp.' ppm' at graph 0.6,0.95
set label 5 at graph 0.55,-0.1 dir noenhance font ',8'
set xlabel 'C_p [ppm]'
set ylabel 'dN/dC_p'
max=1.4 #max value microseconds
min=0
width=(max-min)/(1.0*n) #interval width
set xrange [min:max]
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
set xlabel 'C_p [ppm]'
plot "< awk '{if ($10 > time_off) {print $0}}' ".file u (hist(($21/beta2($9)/$19)*f,width)):(1.0) smooth freq w boxes lc rgb"green" not

#C_b vs dp/p
set xlabel '{/Symbol D}p/p [%]'
set ylabel '{/Symbol D}B [ppm]'
set xrange [-0.5:0.5]
stats file u($9*100):($10 > time_off ? ($32/$19/1.45*1.e6):1/0)
avedb =sprintf("%3.3f",STATS_mean_y)
rmsdb =sprintf("%3.3f",STATS_stddev_y)
set label 1 '<{/Symbol D}B> ='.avedb at graph 0.6,0.1
set label 2 '{/Symbol s}_{B}= '.rmsdb.' ppm' at graph 0.6,0.15
set label 5 at graph 0.06,0.05 dir noenhance font ',10'
plot file u ($9*100):($10 > time_off ?($32/$19)*1.e6/1.45:1/0) pt 5 ps .25 not


# C_b vs horizontal amplitude

set autoscale x
set xlabel 'Horizontal amplitude [mm-mrad]'
set xrange [0:400]
plot file u (($4-$9*eta)**2*gamma + ($5-$9*etap)**2*betax)*1.e6:($10 > time_off ? ($32/$19)*1.e6/1.45:1/0) pt 5 ps .25 not

reset
stats file u (($24)*1.e3):($10 > time_off ? ($24*1.e3):1/0)
rms_x=sprintf("%3.2f",STATS_stddev_y)
ave_x=sprintf("%3.2f",STATS_mean_y)
stats file u (($24)*1.e3):($10 > time_off ? ($24**2*1.e6):1/0)
ave_xe=sprintf("%3.2f",STATS_mean_y)
start_time = sprintf("%2.0f",time_off*1.e6)
print '<x_e> =', ave_x,' <x_e^2> = ',ave_xe,' start time = ', start_time 
print outfile