#!/usr/bin/gnuplot
infile1 = '2017-04-20/histo.dat'
infile2 = '2017-04-20/bi_tx_2017_04_20'
infile3 = '2017-04-20/blue_tx_2017_04_20'
infile4 = '2017-04-20/red_tx_2017_04_20'

set terminal postscript eps color enhanced solid font "Helvetica,18" size 10cm, 20cm
set output "hist_spec_raw.eps"

###### STATS FOR HISTOGRAM NORMALIZATION ##############
stats infile1 u 2 nooutput
sum2= STATS_sum
stats infile1 u 3 nooutput
sum3= STATS_sum
stats infile1 u 4 nooutput
sum4= STATS_sum
stats infile1 u 5 nooutput
sum5= STATS_sum
stats infile1 u 6 nooutput
sum6= STATS_sum
stats infile1 u 7 nooutput
sum7= STATS_sum
stats infile1 u 8 nooutput
sum8= STATS_sum
stats infile1 u 9 nooutput
sum9= STATS_sum

###### FITTING ########################################
lorentz_dis_b(x) = (1 + (4*lambda_b**2 - 4*lambda_b*cos(phi_b))*gamma_b**2/ ( 4*(x0_b-x)**2 + (gamma_b)**2 ) - ( 8*lambda_b*sin(phi_b)*(x-x0_b)*gamma_b ) / ( 4*(x0_b-x)**2 + gamma_b**2 ) ) 
gamma_b = 6.9
x0_b = 40
phi_b = 0.01
lambda_b = 0.047

lorentz_dis_r(x) = (1 + (4*lambda_r**2 - 4*lambda_r*cos(phi_r))*gamma_r**2/ ( 4*(x0_r-x)**2 + (gamma_r)**2 ) - ( 8*lambda_r*sin(phi_r)*(x-x0_r)*gamma_r ) / ( 4*(x0_r-x)**2 + gamma_r**2 ) ) 
gamma_r = 6.9
x0_r = 40
phi_r = 0.01
lambda_r = 0.047

lorentz_dis_bs(x) = (1 + (4*lambda_bs**2 - 4*lambda_bs*cos(phi_bs))*gamma_bs**2/ ( 4*(x0_bs-x)**2 + (gamma_bs)**2 ) - ( 8*lambda_bs*sin(phi_bs)*(x-x0_bs)*gamma_bs ) / ( 4*(x0_bs-x)**2 + gamma_bs**2 ) ) 
gamma_bs = 6.9
x0_bs = 40
phi_bs = 0.01
lambda_bs = 0.047

lorentz_dis_rs(x) = (1 + (4*lambda_rs**2 - 4*lambda_rs*cos(phi_rs))*gamma_rs**2/ ( 4*(x0_rs-x)**2 + (gamma_rs)**2 ) - ( 8*lambda_rs*sin(phi_rs)*(x-x0_rs)*gamma_rs ) / ( 4*(x0_rs-x)**2 + gamma_rs**2 ) ) 
gamma_rs = 6.9
x0_rs = 40
phi_rs = 0.01
lambda_rs = 0.047

lorentz_dis_sum(x) = (1 + (4*lambda_sum**2 - 4*lambda_sum*cos(phi_sum))*gamma_sum**2/ ( 4*(x0_sum-x)**2 + (gamma_sum)**2 ) - ( 8*lambda_sum*sin(phi_sum)*(x-x0_sum)*gamma_sum ) / ( 4*(x0_sum-x)**2 + gamma_sum**2 ) ) 
gamma_sum = 6.9
x0_sum = 40
phi_sum = 0.01
lambda_sum = 0.047

fit lorentz_dis_sum(x) infile2 u (($1-190000)/1000*2):11:12 via  x0_sum, gamma_sum, lambda_sum, phi_sum
fit lorentz_dis_b(x) infile2 u (($1-190000)/1000*2):7:8 via  x0_b, gamma_b, lambda_b, phi_b
fit lorentz_dis_r(x) infile2 u (($1-190000)/1000*2):9:10 via  x0_r, gamma_r, lambda_r, phi_r
fit lorentz_dis_rs(x) infile4 u (($1-190000)/1000*2):4:5 via   x0_rs, gamma_rs, lambda_rs, phi_rs
fit lorentz_dis_bs(x) infile3 u (($1-190000)/1000*2):2:3 via   x0_bs, gamma_bs, lambda_bs, phi_bs

#################################
#######
####GENERAL SETTINGS####
set border linewidth 1.5
set samples 5000
set multiplot 
#set bar 2
unset bar
unset key
set style fill solid 0.3 noborder
set lmargin at screen 0.15
set rmargin at screen 0.9

###########################
#####      SPEC       #####
set xrange [23:58]
set yrange [58:105]
#set ylabel 'Transmission (%)'

#BI-SPEC
set tmargin at screen 0.55
set bmargin at screen 0.25

set ytics 10 nomirror 
set xtics 10 nomirror
set xlabel 'Probe detuning (MHz)'

#plot infile2  u (($1-190000)/1000*2):($11*100):($12*100) w yerr ps 0.8 pt 9 lc rgb '#4B0082' lw 2 not, lorentz_dis_sum(x)*100 lt 1 lc rgb '#4B0082' lw 2 not, infile2  u (($1-190000)/1000*2):($7*100):($8*100) w yerr ps 0.8 pt 7 lc 3 lw 2 not, lorentz_dis_b(x)*100 lc 3 lw 2 not, infile2  u (($1-190000)/1000*2):($9*100):($10*100) w yerr ps 0.8 pt 5 lc 1 lw 2 not, lorentz_dis_r(x)*100 lc 1 lw 2 not

plot lorentz_dis_r(x)*100 lc 1 lw 2 not, lorentz_dis_b(x)*100 lc 3 lw 2 not, lorentz_dis_sum(x)*100 lt 1 lc rgb '#4B0082' lw 2 not, infile2  u (($1-190000)/1000*2):($9*100):($10*100) w yerr ps 0.8 pt 5 lc 1 lw 2 not, infile2  u (($1-190000)/1000*2):($7*100):($8*100) w yerr ps 0.8 pt 7 lc 3 lw 2 not, infile2  u (($1-190000)/1000*2):($11*100):($12*100) w yerr ps 0.8 pt 9 lc rgb '#4B0082' lw 2 not

#BLUE-SPEC
unset xtics
unset xlabel
set border 14
set ytics 80,10,100
set tmargin at screen 0.95
set bmargin at screen 0.65

plot infile3 u (($1-190000)/1000*2):($2*100):($3*100) w yerr ps 0.8 pt 7 lw 2 lc 3 not,  lorentz_dis_bs(x)*100 lc 3 lw 2 not

#RED-SPEC
set ylabel 'Transmission (%)'
set border 14
set tmargin at screen 0.75
set bmargin at screen 0.45

plot infile4 u (($1-190000)/1000*2):($4*100):($5*100) w yerr ps 0.8 pt 5 lw 2 lc 1 not,  lorentz_dis_rs(x)*100 lc 1 lw 2 not

set border 15
###########################
#####      HISTO      #####
set lmargin at screen 0.7
set rmargin at screen 0.95
set xtics nomirror
set ytics nomirror
unset ylabel
#set ylabel 'norm. frequency' #offset 1.5
set xlabel 'photons' offset 0,0.5

#BLUE HISTO
set xrange [0:110]
set yrange [0:0.18]
set ytics 0.1 offset 0.5
set xtics 50

set tmargin at screen 0.85
set bmargin at screen 0.8



plot infile1  u 1:($4/sum4) lc -1 with boxes , '' u 1:($3/sum3) lc 3 lw 3 with steps

#RED HISTO
set tmargin at screen 0.65
set bmargin at screen 0.6
plot infile1  u 1:($7/sum7) lc -1 with boxes , '' u 1:($6/sum6) lc 1 lw 3 with steps

#BI HISTO
set xrange [0:200]
set xtics 100

set tmargin at screen 0.35
set bmargin at screen 0.3
plot infile1  u 1:($9/sum9) lc -1 with boxes , '' u 1:($8/sum8) lc rgb '#4B0082' lw 3 with steps

unset xlabel
unset xtics
#blue
set tmargin at screen 0.45
set bmargin at screen 0.4
plot infile1  u 1:($8/sum8) lc -1 with boxes , '' u 1:($2/sum2) lc 3 lw 3 with steps

#red
set tmargin at screen 0.4
set bmargin at screen 0.35
plot infile1  u 1:($8/sum8) lc -1 with boxes , '' u 1:($5/sum5) lc 1 lw 3 with steps




unset multiplot   
unset output