print 'input: filename=dir/all_EndOfTracking_phase_space.dat; events=100; time_off=35e-6'
set terminal pdf enhanced fontscale 0.75 size 6.0in, 4.0in
outputfile = filename.'_hist.pdf' 
print outputfile
set output outputfile
# inputs required: filename =<EndOfTracking_phase_space.dat> (phase space file), events = 'number of events in histogram'
set key right inside samplen 1
#stats filename u 10
Y_MIN=0
nturns = 100

set timestamp offset graph 1.08,0.65 rotate font 'Verdana,6'
set label 3 'plotting_scripts/hist_any_phase_space_file_mp.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
set label 4 filename at graph 0.2,1.03 font 'Verdana,8' noenhance
print 'histogram mean radial position next'
n=100 #number of intervals
scale=events/n
n=events
max=50. #max value
min=-50 #min value
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 boxwidth width*0.9
set style fill solid 0.5 #fillstyle
set tics out nomirror
set xlabel  "<radial offset> [mm]"
stats filename 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 filename 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)
set label 1 'x_e = '.ave_x.' mm' at graph 0.05,0.95
set label 2 '{/Symbol s}_{x_e} = '.rms_x.' mm' at graph 0.7,0.95
set label 8 '<x_e^2> = '.ave_xe.' mm^2' at graph 0.05,0.9
set label 7 'start = '.start_time.'{/Symbol m}s' at graph 0.05,0.81
plot "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u (hist(($24)*1.e3,width)):(1.0) smooth freq w boxes lc rgb"green" not

print 'histogram of x at end of track'
set xlabel  "radial offset [mm]"
set timestamp offset graph 1.08,0.8 rotate font 'Verdana,6'
stats filename u (($4)*1.e3):($10 > time_off ? ($4*1.e3):1/0)
rms_x=sprintf("%3.2f",STATS_stddev_y)
ave_x=sprintf("%3.2f",STATS_mean_y)
set label 1 '<x> = '.ave_x.' mm' at graph 0.05,0.95
set label 2 '{/Symbol s}_{x} = '.rms_x.' mm' at graph 0.7,0.95

plot  "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u (hist(($4)*1.e3,width)):(1.0) smooth freq w boxes lc rgb"green" not

#print GPVAL_Y_MAX, scale
#Y_MAX = (0+floor(GPVAL_Y_MAX/scale)) * scale
#set ytics 0,(Y_MAX-Y_MIN)/5.,GPVAL_Y_MAX 
#clear
#replot
set timestamp offset graph 1.08,0.8 rotate font 'Verdana,6'


print 'histogram of average y next'
max=50. #max value
min=-50. #min value
width=(max-min)/n #interval width
set xrange [min:max]
#set yrange [0:10]
#to put an empty boundary around the
#data inside an autoscaled graph.
#set offset graph 0.05,0.05,0.05,0.0
set xtics min,(max-min)/10,max
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
set tics out nomirror
set xlabel "<y> [mm]" 
stats filename u ($26*1.e3):($10 > time_off ? ($26*1.e3):1/0)
rms_x=sprintf("%3.2f",STATS_stddev_x)
ave_x=sprintf("%3.2f",STATS_mean_x)
set label 1 '<y> = '.ave_x.' mm' at graph 0.05,0.95
set label 2 '{/Symbol s}_{<y>} = '.rms_x.' mm' at graph 0.7,0.95
plot  "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u (hist(($26)*1.e3,width)):(1.0) smooth freq w boxes lc rgb"green" not "y[mm]"

print 'histogram of y at end of tracking next'
set xlabel "y [mm]" 
stats filename u (($6*1.e3)):($10>time_off ? (($6*1.e3)**2):1/0)
ave_y2=sprintf("%3.2f",STATS_mean_y)
stats filename u (($6*1.e3)):($10>time_off ? ($6*1.e3):1/0)
rms_y=sprintf("%3.2f",STATS_stddev_y)
ave_y=sprintf("%3.2f",STATS_mean_y)
print 'ave_y2= '.ave_y2,' rms_y ='.rms_y, ' ave_y ='.ave_y
set label 1 '<y> = '.ave_y.' mm' at graph 0.05,0.95
set label 2 '{/Symbol s}_y = '.rms_y.' mm' at graph 0.7,0.95
set label 9 '<y^2> = '.ave_y2.' mm^2' at graph 0.7,0.9
plot  "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u (hist(($6)*1.e3,width)):(1.0) smooth freq w boxes lc rgb"green" not "y[mm]"

print 'histogram momentum next'
max=.5  #max value
min =-0.5 #min value
width=(max-min)/n #interval width
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 xlabel  "dp/p[%]"
stats filename u ($9*100):($10 > time_off ? ($9*100):1/0)
sig_p = sprintf("%3.2f",STATS_stddev_y)
p_ave = sprintf("%3.2f",STATS_mean_y)
set label 1 '{/Symbol s}_p = '.sig_p.'[%]' at graph 0.05,0.95
set label 2 '<{/Symbol D}p/p> = '.p_ave.' [%]' at graph 0.65,0.95

plot "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u (hist(($9)*1.e2,width)):(1.0) smooth freq w boxes lc rgb"green" not "dp/p[%]"

print 'histogram time'
set autoscale
width=1
stats filename u ($10*1.e6):($10> time_off ? ($10*1.e6):1/0)
max=floor(STATS_max_x -STATS_mean_x +1) #max value
min=floor(STATS_min_x-STATS_mean_x) #min value
width=(max-min)/n*2 #interval width
set xrange [0:]
print ' min= ', min, ' max= ', max, ' width = ', width,'  STATS_min_x = ', STATS_min_x,' STATS_max_x = ', STATS_max_x,' STATS_mean_x = ', STATS_mean_x
#set xtics min,(max-min)/10
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
set xtics 0,64
set xlabel "time [{/Symbol m}s]" 
#stats filename u ($10)*1.e6
rms_x=sprintf("%3.2f",STATS_stddev_y)
ave_x=sprintf("%3.2f",STATS_mean_y-time_off*1.e6)
Nmuons=sprintf("%5i",STATS_records)
print  'Number of muons = ',Nmuons
unset label
set label 1 '<t> = '.ave_x.' {/Symbol m}s' at graph 0.05,0.95
set label 2 '{/Symbol s}_t = '.rms_x.' {/Symbol m}s' at graph 0.7,0.95
unset label 2
set label 3 Nmuons.' muons' at graph 0.05,0.9
#plot "< awk '{if ($10 > time_off){print $0}}' ".filename u (hist(($10*1.e6-STATS_mean_x),width)):(1.0) smooth freq w boxes lc rgb "green" not "t[ns]"
#plot "< awk '{if ($10 > time_off){print $0}}' ".filename u (hist(($10*1.e6),width)):(1.0) smooth freq w boxes lc rgb "green" not "t[ns]"
plot filename u (hist(($10*1.e6),width)):(1.0) smooth freq w boxes lc rgb "green" not "t[ns]"

#
reset
set timestamp offset graph 1.08,0.8 rotate font 'Verdana,6'
set label 3 'plotting_scripts/hist_any_phase_space_file_mp.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
set label 4 filename at graph 0.2,1.03 font 'Verdana,8' noenhance
max=6.75 #max value
min=6.65 #min value
width=(max-min)/n #interval width
print 'width = ', width
set autoscale
set xrange [min:max]
set xtics min,0.01
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
set xlabel "cyclotron frequency [MHz] " 
#stats filename u ($23)*1.e-6
stats filename u ($23*1.e-6):($10 > time_off ? ($23*1.e-6):1/0)
rms_x=sprintf("%3.2f",STATS_stddev_x)
ave_x=sprintf("%3.2f",STATS_mean_x)
set label 1 '<f> = '.ave_x.' [MHz]' at graph 0.05,0.95
set label 2 '{/Symbol s}_f = '.rms_x.' MHz' at graph 0.7,0.95
set label 7 'start = '.start_time.'{/Symbol m}s' at graph 0.05,0.90
set grid
set mxtics 5
plot "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u (hist(($23)*1.e-6,width)):(1.0) smooth freq w boxes lc rgb "green" not
#clear
xdif = (GPVAL_DATA_X_MAX-GPVAL_DATA_X_MIN)/5
#
set autoscale
set timestamp offset graph 1.06,0.8 rotate font 'Verdana,6'
set label 3 'plotting_scripts/hist_any_phase_space_file_mp.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
set label 4 filename at graph 0.1,1.03 font 'Verdana,8' noenhance
set xlabel '{/Symbol D}p/p [%]'
set ylabel 'Cyclotron Frequency [MHz]'
set xrange [-0.4:0.4]
set yrange [6.65:6.75]
set ytics 6.65,0.02
set xtics -0.4,0.2
plot "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u ($9*100):($23/1.e6)  w p pt 5 ps .25 not

#
set autoscale
reset
set timestamp offset graph 1.06,0.8 rotate font 'Verdana,6'
set label 3 'plotting_scripts/hist_any_phase_space_file_mp.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
set label 4 filename at graph 0.1,1.03 font 'Verdana,8' noenhance
set xlabel "Amplitude ((x-{/Symbol hd})^2/{/Symbol b}+x'^2{/Symbol b}) [mm-mrad]"
set ylabel 'Hits'
width=0.5
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
beta = 7.6969
eta= 8.045
plot "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u (hist((($4-$9*eta)**2/beta+$5**2*beta)*1.e6,width)):(1.0) smooth freq w boxes lc rgb "green" not

#
set autoscale
reset
set timestamp offset graph 1.06,0.8 rotate font 'Verdana,6'
set label 3 'plotting_scripts/hist_any_phase_space_file_mp.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
set label 4 filename at graph 0.1,1.03 font 'Verdana,8' noenhance
set xlabel "Amplitude ((x-{/Symbol hd})^2/{/Symbol b}+x'^2{/Symbol b}) [mm-mrad]"
set ylabel 'Hits'
width=0.5
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
beta = 7.6969
eta= 8.045
plot "< awk '{if (($10 > time_off)&&($9<0.001)&&($9>-0.001)){print $0}}' ".filename u (hist((($4-$9*eta)**2/beta+$5**2*beta)*1.e6,width)):(1.0) smooth freq w boxes lc rgb "green" t '|{/Symbol D}p/p|<0.001'

#
set autoscale
reset
set timestamp offset graph 1.06,0.8 rotate font 'Verdana,6'
set label 3 'plotting_scripts/hist_any_phase_space_file_mp.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
set label 4 filename at graph 0.1,1.03 font 'Verdana,8' noenhance
set xlabel "Amplitude (x^2/{/Symbol b}+x'^2{/Symbol b}) [mm-mrad]"
set ylabel 'Hits'
width=width
set boxwidth width*0.9
set style fill solid 0.5 #fillstyle
beta = 7.6969
eta= 8.045
plot "< awk '{if (($10 > time_off)&&($9<0.01)&&($9>-0.01)){print $0}}' ".filename u (hist((($4)**2/beta+$5**2*beta)*1.e6,width)):(1.0) smooth freq w boxes lc rgb "green" not

reset
set timestamp offset graph 1.06,0.8 rotate font 'Verdana,6'
set label 3 'plotting_scripts/hist_any_phase_space_file_mp.gnu' at graph 1.02,0.02 rotate left font 'Verdana,6' noenhance
set label 4 filename at graph 0.1,1.03 font 'Verdana,8' noenhance
set ylabel "Amplitude ((x-{/Symbol hd})^2/{/Symbol b}+x'^2{/Symbol b}) [mm-mrad]"
set xlabel '<x_e> [mm]'
plot filename u ($24*1000):($10>time_off ? (($4-$9*eta)**2/beta+$5**2*beta)*1.e6:1/0) pt 5 ps .5 not

print outputfile