# type the following command in gnuplot, in the directory that contains the .gnu and .txt datafile, as below:
# load "shotnoise.gnu"

# use this setting for .eps, common choice for use in LaTeX
set terminal postscript eps color enhanced font "Helvetica,16pt" size 8.6cm,6cm 

# use this for .png, everyday use easier to open
#set terminal pngcairo enhanced font "Helvetica,24pt" size 1920,1080

# common formatting settings
set format y "%2.0f"
set format x "%2.0f"
set tics scale 2
set tics nomirror
set xrange [0:400]
set xtics 50
#set mxtics 5
set yrange [0:7000]
set ytics 1000
#set mytics 5
set border linewidth 2

set xlabel "Timing (seconds)"
set ylabel "Photoevents per 100ms"


# comparison countrates of Dark v Bright v Total over Position
set output 'inversion_shotnoise_countrates.png'
plot 'shotnoise.txt' using ($1*0.1):3:(sqrt($3)) w yerrorlines lc rgb "red" lw 1 title 'Higher Order Modes', 'shotnoise.txt' using ($1*0.1):2:(sqrt($2)) w yerrorlines lc rgb "blue" lw 1 title 'Common Mode', 'shotnoise.txt' using ($1*0.1):($2+$3):(sqrt($2+$3)) w yerrorlines lc rgb "black" lw 1 title 'Total'


set xlabel "Photoevents per 100ms"
set ylabel "Number of Photoevents"
# histogram distribution of data columns

set xr [107.8-300:107.8+300]
set xtics 50
set yr [0:200]
set ytics 50
binwidth=1.0087
bin(x,width)=width*floor(x/width)
#set output 'inversion_shotnoise_distribution_higher.png'
set table 'histogram_higher.txt'
plot 'shotnoise.txt' using (bin($3,binwidth)):(1.0) smooth freq with boxes lc rgb 'red' title 'Higher Modes noise distribution'
unset table
gaussh(x)=ah/(sigmah*sqrt(2*pi))*exp(-(x-muh)**2/(2*sigmah**2))
ah=4000
sigmah=11
muh=107.8
set fit errorvariables
	FIT_LIMIT=1.e-15
	FIT_MAXITER=1000
fit gaussh(x) 'histogram_higher.txt' using 1:2 via ah, sigmah, muh
set output 'inversion_shotnoise_distribution_higher.png'
set samples 10000
set xr [-300:300]
muh=0
plot 'histogram_higher.txt' using ($1-107.8):2:(sqrt($2)) w yerrorbars lw 2 ps 2 lc rgb 'red' title 'Higher Modes noise distribution', gaussh(x) w l lw 3 lc rgb 'black' title 'sigma = 11 , sqrtN = 10'

set xr [5086.5-300:5086.5+300]
set xtics 50
set yr [0:200]
set ytics 50
binwidth=7.02
bin(x,width)=width*floor(x/width)
#set output 'inversion_shotnoise_distribution_common.png'
set table 'histogram_common.txt'
plot 'shotnoise.txt' using (bin($2,binwidth)):(1.0) smooth freq with boxes lc rgb 'blue' title 'Common Mode noise distribution'
unset table
gaussc(x)=ac/(sigmac*sqrt(2*pi))*exp(-(x-muc)**2/(2*sigmac**2))
ac=28000
sigmac=75
muc=5086.5
set fit errorvariables
	FIT_LIMIT=1.e-15
	FIT_MAXITER=1000
fit gaussc(x) 'histogram_common.txt' using 1:2 via ac, sigmac, muc
set output 'inversion_shotnoise_distribution_common.png'
set samples 10000
set xr [-300:300]
muc=0
plot 'histogram_common.txt' using ($1-5086.5):2:(sqrt($2)) w yerrorbars lw 2 ps 2 lc rgb 'blue' title 'Common Mode noise distribution', gaussc(x) w l lw 3 lc rgb 'black' title 'sigma = 73 , sqrtN = 71'


set xr [5194.7-300:5194.7+300]
set xtics 50
set yr [0:200]
set ytics 50
binwidth=7.1
bin(x,width)=width*floor(x/width)
#set output 'inversion_shotnoise_distribution_total.png'
set table 'histogram_total.txt'
plot 'shotnoise.txt' using (bin($4,binwidth)):(1.0) smooth freq with boxes lc rgb 'black' title 'Total noise distribution'
unset table
gausst(x)=at/(sigmat*sqrt(2*pi))*exp(-(x-mut)**2/(2*sigmat**2))
at=28000
sigmat=75
mut=5194.7
set fit errorvariables
	FIT_LIMIT=1.e-15
	FIT_MAXITER=1000
fit gausst(x) 'histogram_total.txt' using 1:2 via at, sigmat, mut
set output 'inversion_shotnoise_distribution_total.png'
set samples 10000
set xr [-300:300]
mut=0
plot 'histogram_total.txt' using ($1-5194.7):2:(sqrt($2)) w yerrorbars lw 2 ps 2 lc rgb 'black' title 'Total noise distribution', gausst(x) w l lw 3 lc rgb 'black' title 'sigma = 74 , sqrtN = 72'


set xr [-250:250]
set yr [0:215]
set ytics 50
set xtics 100
set xlabel "Difference {/Symbol D}N from mean N (Photoevents per 100 ms)"
set ylabel "Number of Occurences"
set output 'shotnoise.eps'

plot 'histogram_higher.txt' using ($1-107.8):2:(sqrt($2)) w yerrorbars lw 2 ps 0 lc rgb 'red' title 'Asymmetric Modes', gaussh(x) w l lw 2 lc rgb 'red' notitle 'Higher sqrtN = 10', 'histogram_common.txt' using ($1-5086.5):2:(sqrt($2)) w yerrorbars lw 2 ps 0 lc rgb 'blue' title 'Symmetric Modes', gaussc(x) w l lw 2 lc rgb 'blue' notitle 'Common sqrtN = 71', 'histogram_total.txt' using ($1-5194.7):2:(sqrt($2)) w yerrorbars lw 2 ps 0 lc rgb 'black' title 'Total Modes', gausst(x) w l lw 2 lc rgb 'black' notitle 'Total sqrtN = 72'


reset