set term pdf fontscale 0.5 size 6in,4in
#x-px phase space
set timestamp offset graph 1.14,0.7 rotate font ',6'
set label 'plotting scripts/2D_hist_test.gnu' at graph 1.01,0.02 rotate left font 'Verdana,6' noenhance
#the following 1-create a file named directory. 2-modify the file adding name=" at the beginning and " at the end. Then loading. Now the character string <name> is the path
!pwd > directory
!sed -n -e 's/^/name="/' -e 's/$/"/p' directory > name
load 'name'
print 'directory='.name
set label at graph 0.,1.02 name font 'Verdana,6'noenhance 

cnumber = sprintf("%1i",number)
print 'number=',number,'  cnumber=',cnumber

load '2D_hist_input_'.cnumber.'.dat'

outfile = '2D_histogram_'.cnumber.'.pdf'
print outfile
set output outfile

yf=1
xf=1
if(number == 6){
ymin=ymin*1.e6
ymax=ymax*1.e6
yf=1.e6
filename = 'bin_average_rms_6.dat'
D_filename = '2D_hist_6.dat'
load '2D_hist_fit_coef_err_6.dat'
}
if(number == 3){
xmin=xmin*1.e6
xmax=xmax*1.e6
xf=1.e6
filename = 'bin_average_rms_3.dat'
D_filename = '2D_hist_3.dat'
}
if(number == 4){
filename = 'bin_average_rms_4.dat'
D_filename = '2D_hist_4.dat'
load '2D_hist_fit_coef_err_4.dat'
}
if(number == 5){
filename = 'bin_average_rms_5.dat'
D_filename = '2D_hist_5.dat'
load '2D_hist_fit_coef_err_5.dat'
}
print filename

set key off
set view map
 set pm3d at s hidden3d 100
 set surface
 set dgrid3d 100,100,3
 set view map  scale 1
 set style increment default
 set style line 100 lt 5 lw 0.5
 unset hidden3d
 unset surf
 set pm3d interpolate 0,0
# set autoscale fix
# set contour base
# unset dgrid
set pm3d noborder
set palette defined ( 0 0.05 0.05 0.2, 0.1 0 0 1, 0.25 0.7 0.85 0.9,\
   0.4 0 0.75 0, 0.5 1 1 0, 0.7 1 0 0, 0.9 0.6 0.6 0.6,\
  1 0.95 0.95 0.95 )
    # set palette defined (0 0 0 0, 1 0 0 1, 3 0 1 0, 4 1 0 0, 6 1 1 1)
     set palette color
#     set palette defined ( 0 "white", 1 "black", 2 "blue" , 3"green")

if(number == 6){
set ylabel 'time [{/Symbol m}s]'
set xlabel '{/Symbol D}p/p [%]'
}
if(number == 3){
set xlabel 'time [{/Symbol m}s]'
set ylabel '{/Symbol D}p/p [%]'
}
if(number == 4){
set xlabel 'momentum [{/Symbol D}p/p]'
set ylabel 'x_{inf} [m]'
}
if(number == 5){
set xlabel 'momentum [{/Symbol D}p/p]'
set ylabel "x'_{inf} [rad]"
}


 xscale = (xmax-xmin)/ix_bins
 yscale=(ymax-ymin)/iy_bins
 yconvert(x) = yscale * x + (ymax+ymin)/2.
#set xrange [-0.3:0.5]
set yrange [ymin:ymax]
print D_filename
splot D_filename u ($1*xscale):(yconvert($2)):3 w pm3d
#pause -1
reset

set ylabel '<time> [{/Symbol m}s]'
set xlabel '{/Symbol D}p/p'

f(x) = offset + slope*x + q*x*x +r*x*x*x

reset
set label at graph 0.,1.02 name font 'Verdana,6'noenhance 
set timestamp offset graph 1.03,0.7 rotate font 'Verdana,6'
set label 'plotting scripts/2D_hist_test.gnu' at graph 1.015,0.02 rotate left font 'Verdana,6' noenhance

if(number == 6){
set ylabel '<time> [{/Symbol m}s]'
set xlabel '{/Symbol D}p/p'
legend  = 'd<t>/d{/Symbol d}'
units = '[{/Symbol m}s]'
units0 = " [{/Symbol m}s]"
units1 = " [{/Symbol m}s]"
units2=" [{/Symbol m}s]"
units3=" [{/Symbol m}s]"
}
if(number == 3){
set xlabel 'time [{/Symbol m}s]'
set ylabel '<{/Symbol D}p/p>'
legend  = 'd<t>/d{/Symbol d} = '
units = '[{/Symbol m}s]'
}
if(number == 4){
set ylabel '<x_{inf}> [m]'
set xlabel '{/Symbol D}p/p'
legend  = 'd<x_{inf}>/d{/Symbol d} = '
units='[m]'
units0 = " [m]"
units1 = " [m]"
units2=" [m]"
units3=" [m]"
}
if(number == 5){
set ylabel "<x'_{inf}> [rad]"
set xlabel '{/Symbol D}p/p'
legend  = "d<x'_{inf}>/d{/Symbol d} = "
units = "[rad]"
units0 = " [rad]"
units1 = " [rad]"
units2=" [rad]"
units3=" [rad]"
}

fit f(x) filename u ($1*xf):($2*yf):($5 > 2 ? (1./$5**.5):1/0) yerror via slope, offset,q, r

label0 = 'p0= '.sprintf("%5.4f",offset).'{/Symbol \261}'.sprintf("%5.4f",offset_err).units0
legends = legend.sprintf("%5.4f",slope).'{/Symbol \261}'.sprintf("%5.4f",slope_err).units
label1 = 'p1= '.sprintf("%5.4f",slope).'{/Symbol \261}'.sprintf("%5.4f",slope_err).units1
label2 = 'p2= '.sprintf("%5.4f",q).'{/Symbol \261}'.sprintf("%5.4f",q_err).units2
label3 = 'p3= '.sprintf("%5.4f",r).'{/Symbol \261}'.sprintf("%5.4f",r_err).units3

set label ' fitbinned data' at graph 0.05,0.98
set label label0 at graph 0.05,0.93
set label label1 at graph 0.05,0.88
set label label2 at graph 0.05,0.83
set label label3 at graph 0.05,0.78

set timestamp offset graph 1.14,0.7 rotate font ',6'

#plot filename u ($1*xf):($2*yf):($4*yf) w err not, f(x) t  legend.' ='.sprintf("%6.4f",slope).' [{/Symbol m}s]'
plot filename u ($1*xf):($2*yf):($4*yf) w err not, f(x) not  legends 

 
g(x) = coef0 + coef1*x + coef2*x**2 + coef3*x**3
unset label

set timestamp offset graph 1.14,0.7 rotate font ',6'
set label at graph 0.,1.02 name font 'Verdana,6'noenhance 
set label 'plotting scripts/2D_hist_test.gnu' at graph 1.015,0.02 rotate left font 'Verdana,6' noenhance

label0 = 'p0= '.sprintf("%5.4f",coef0).'{/Symbol \261}'.sprintf("%5.4f",coef0_err).units0
label1 = 'p1= '.sprintf("%5.4f",coef1).'{/Symbol \261}'.sprintf("%5.4f",coef1_err).units1
label2 = 'p2= '.sprintf("%5.4f",coef2).'{/Symbol \261}'.sprintf("%5.4f",coef2_err).units2
label3 = 'p3= '.sprintf("%5.4f",coef3).'{/Symbol \261}'.sprintf("%5.4f",coef3_err).units3

set label 'fit unbinned data' at graph 0.05, 0.98
set label label0 at graph 0.05,0.93
set label label1 at graph 0.05,0.88
set label label2 at graph 0.05,0.83
set label label3 at graph 0.05,0.78

plot filename u ($1*xf):($2*yf):($4*yf) w err not, g(x) not '{/Symbol x}^2='.sprintf("%2.4f",chisq)


#time_to_phase = t/T*2*pi = t*omega_a = t *gamma/(gamma^2-1)* omega_c = t*gamma/(gamma^2-1)*2*pi/.1492
gamma=29.3
time_to_phase_factor = gamma/(gamma**2-1)*2*pi/.1492
ttp = time_to_phase_factor

if(number == 6){
reset


set ylabel '<{/Symbol f}_a> [rad]'
set xlabel '{/Symbol D}p/p'


f(x) = offset + slope*x +q*x*x

 fit f(x) filename u ($1):($2*yf*ttp):($5>2 ? (1./$5**.5):1/0) yerror via slope, offset,q


plot filename u ($1):($2*yf*ttp):($4*yf*ttp) w err not, f(x) t  'fit- d<t>d{/Symbol d} ='.sprintf("%6.4f",slope).' [rad]'

}

print outfile