#================================================================
# various settings
#================================================================
INPUT  = "acc_part.dat"
OUTPUT = "E_hist"
OUTPUTA = sprintf("%s.eps", OUTPUT);
OUTPUTB = sprintf("%s.pdf", OUTPUT);


iGun = 500   # [mA]
NPS  = 50000  # macro particles in GPT

G2E  = 0.511
dE   = 0.015      # [MeV]

N2I = iGun/(NPS*dE) # [mA/(Num MeV)]

Xmin = 0            # [MeV]
Xmax = 10.5	    # [MeV]
Xsamples = (Xmax-Xmin)/dE+1

#================================================================
# Statistics analysis
#================================================================

#------ find max and peak energy (kernel density estimation) ------

DENSITY = "density.txt"
set table DENSITY
set samples Xsamples
plot INPUT using (($7-1)*G2E):(dE):(dE) smooth kdensity
stats DENSITY using ($1):($2)
unset table

MaxEnergy = STATS_max_x
PekEnergy = STATS_min_x+(STATS_max_x - STATS_min_x)/(STATS_records-1)*STATS_index_max_y

#----- Analysis all energy -----------

set xrange [0.0: 1.1*MaxEnergy]
set yrange [0.0: 1.1*MaxEnergy]

stats INPUT using (($7-1)*G2E)
N100     = STATS_records
meanE100 = STATS_mean
Ib100   = iGun*N100/NPS
SigmaE100 = STATS_stddev

#----- Analysis over 70% of peak energy -----------

set xrange [0.8*PekEnergy: 1.1*MaxEnergy]
set yrange [0.8*PekEnergy: 1.1*MaxEnergy]

stats INPUT using (($7-1)*G2E)
N080     = STATS_records
meanE080 = STATS_mean
Ib080   = iGun*N080/NPS
SigmaE080 = STATS_stddev

#----- save data to file ------------------------
set print "analysis.txt"
result100 = sprintf("%.3f\t%.3f\t%.3f\t%.3f\t%.2f\t", MaxEnergy, PekEnergy, meanE100, SigmaE100, Ib100)
result080 = sprintf("%.3f\t%.3f\t%.2f", meanE080, SigmaE080, Ib080)

print "-------------------------------------------------------"
print "E range 100%\t\t\t\trange 80%"
print "Max E\tPeak E\tEav\tSigmaE\tIb\tEav\tSigmaE\tIb"
print "[MeV]\t[MeV]\t[MeV]\t[MeV]\t[mA]\t[MeV]\t[MeV]\t[mA]"
print "======================================================="
print result100, result080
print "-------------------------------------------------------"
set print

#================================================================
#  screen plot
#================================================================
#---- etc -------
set terminal wxt enhanced
set key off

#---- x axis -------
set xrange [0:Xmax]
set xtics border mirror in scale 3,1.5 font ",22" 0,1
set format x "%.0f"
set mxtics 5
set xlabel "Energy [MeV]" font ",24"

#---- y axis -------
set yrange [0:2000]
set ytics border mirror in scale 3,1.5 font ",22" 0,500
set mytics 5
set format y "%.0f"
set ylabel "Intensity [mA/MeV]" font ",24"


#---- for histgram ----
bin(x) = dE*floor((x+0.5*dE)/dE)


#---- label --------
labelEmax = sprintf("E_{max} : %.3f [MeV]", MaxEnergy)
labelEpek = sprintf("E_{peak} : %.3f [MeV]", PekEnergy)
labelEav8 = sprintf("<E>_{080} : %.3f [MeV]", meanE080)
labelS080 = sprintf("{/Symbol s}_{E080} : %.3f [MeV]", SigmaE080)
labelI100 = sprintf("I_{100} : %.2f [mA]", Ib100)
labelI080 = sprintf("I_{080} : %.2f [mA]", Ib080)

set label labelEmax at graph 0.67,0.90 tc rgb "black" font ",20"
set label labelEpek at graph 0.67,0.85 tc rgb "black" font ",20"
set label labelEav8 at graph 0.67,0.80 tc rgb "black" font ",20"
set label labelS080 at graph 0.67,0.75 tc rgb "black" font ",20"
set label labelI100 at graph 0.67,0.70 tc rgb "black" font ",20"
set label labelI080 at graph 0.67,0.65 tc rgb "black" font ",20"

#-------- plot --------------------
plot INPUT using (bin(($7-1)*G2E)):(N2I) smooth frequency with histeps lw 1.5 lt 7 lc rgb "red" notitle

pause 3 "press [Enter] key to save and quit"

#----------------------------------------------------------
set terminal push
set terminal postscript eps color solid enhanced font "Helvetica,18"
set output OUTPUTA
replot
#set terminal pop
unset output

#---------- fix bouding box ----
command = sprintf("epstopdf %s", OUTPUTA)
system command