! Anything after an exclamation mark is a comment
!
! Read in the formatted data into vectors x and y
read\format\nomessage gauss+bkgd.asc (*) x y




! Define which variables are free to vary in the fit
scalar\vary a b c bkgd s_bkgd

! Give them initial values
a      = max(y)  ! the amplitude of the Gaussian
b      = 777     ! the mean of the Gaussian
c      =   1.2   ! the width (sigma) of the Gaussian
bkgd   = 200     ! the background level
s_bkgd =  -1     ! the slope of the background

! fit\nomess => fit without writing stuff to screen.  This doesn't include 
!               weights.  But, it is more likely to converge.  Good for 
!               a first guess.
!
! Define the fit function as a Gaussian, a*exp(-0.5*((x-b)/c)**2)/sqrt(2*pi)/c
!                         plus a sloping background, (bkgd-s_bkgd*(x-b)
!     The Gaussian is properly normalized, and the slope of the background 
!     starts at its centroid
fit\nomess y=a*exp(-0.5*((x-b)/c)**2)/sqrt(2*pi)/c+(bkgd-s_bkgd*(x-b))

! now include weights, have it put the "E1" (statistical) uncertainties in 
! the vector fit$e1, the total chi^2 in the scalar fit$chisq, the degrees of 
! freedom in the scalar fit$free, the confidence level of the fit in the 
! scalar fit$cl, and also the correlation and covariance matrices in 
! fit$corr and fit$covm respectively
fit\weights\e1\e2\chisq\free\cl\corr\cov 1/y y=a*exp(-0.5*((x-b)/c)**2)/sqrt(2*pi)/c+(bkgd-s_bkgd*(x-b))

! Let's save the result in vector yfit
yfit = a*exp(-0.5*((x-b)/c)**2)/sqrt(2*pi)/c+(bkgd-s_bkgd*(x-b))

! So we can plot the Gaussian separately, define vector ysignal
ysignal = a*exp(-0.5*((x-b)/c)**2)/sqrt(2*pi)/c
! and similarly ybkgd
ybkgd   = (bkgd-s_bkgd*(x-b))


! save the total number of signal events into scalar ysignal_sum
stat ysignal ysignal_sum\sum
! similarly for the background
stat ybkgd ybkgd_sum\sum
! can define the signal to noise this way
s2n = ysignal_sum/ybkgd_sum

! make the plot orientation portrait
orient port
! set the colour to black
col 1
! clear the graphics
cl
! set up some things
set
  histyp 0                      ! turn off histogramming
  pchar -13 .6                  ! plotting character a filled circle, radius 0.6
  %ylaxis 30                    ! the lower y-axis, % of the screen
  %yuaxis 95                    ! the upper y-axis, % of the screen
  xaxis 0                       ! turn off numbering for the x-axis
  postres 180                    ! the rest are just what I think is pretty
  linthk 4
  xlabsz 0.5
  ylabsz 0.5
  xnumsz 0.5
  ynumsz 0.5
  %xtics .5
  %ytics .5
  %xticl 1.25
  %yticl 1.25
  txthit 0.5
  nsyinc 5
  nsxinc 5
  font triumf.2
  cursor -9
  yleadz 1
  xleadz 1

! set the scales, xlow xhi ylow yhi (0 0 means autoscale)
scales 760 790 0 1700
gr\axes

set lintyp 1
set pchar -13 0.6
col 1
! Graph the data with error bars
gr\noax x y sqrt(y)


! change to a solid line
set pchar 0

! make the colour blue
col 4
! without replotting the scales, graph the background
gr\noax x ybkgd

! make the colour green
col 3
! overlay the total fit (including background)
gr\noax x yfit

! set colour back to black for residuals plot
col 1
! calculate the residuals
yresid = (yfit-y)/sqrt(yfit)

! change from points to a histogram, and y-axis upper/lower; turn x-axis on
set
  histyp 1                      ! make it a histogram
  %ylaxis 12                    ! change lower/upper yaxis so below main plot
  %yuaxis 30
  xaxis 1                       ! turn on x-axis (show numbers)

!  Vertical scales go from -3 to 3
scales 760 790 -3 2.99999

! graph the residuals
gr\axes
col 1
gr\noax x yresid
set histyp 0
col 3 ! green
gr\noax x 0*x                   ! graph a zeroline
set lintyp 9
gr\noax x 0*x+1
gr\noax x 0*x-1
set lintyp 1
col 1 ! black
set
  txthit 0.45  ! height of text
  %xloc 17     ! location of text
  %yloc 92.5
  cursor -7    ! left-aligned

! some definitions used later on
ital  = `<fitalic.2>'   ! stuff in quotes is a text variable; 
roman = `<ftriumf.2>'  ! gets evaluated when used
text ital//`F'//roman//`('//ital//`x'//roman//`) = Gaussian'

set %xloc 19
set %yloc 85
! Spit out fit result information
text ital//`x'//roman//`<_>0<^> = '//rchar(b,`f7.3')//`('//rchar(NINT(fit$e1(2)*1e3))//`)'
set %yloc 81
text `<sigma> = '//rchar(c,`f5.3')//`('//rchar(NINT(fit$e1(3)*1e3))//`)'
set %yloc 77
text `area = '//rchar(NINT(a))//`('//rchar(NINT(fit$e1(1)))//`)'
set %yloc 71
text `<chi,^>2<_>/'//rchar(fit$free)//` = '//rchar(fit$chisq/fit$free,`f6.4')
set %yloc 67
text `C.L. = '//rchar(fit$cl,`f4.1')//`%'
set %yloc 33
text ital//`y'//roman//`<_>bkgd<^> = '//rchar(bkgd,`f5.1')//`('//rchar(NINT(fit$e1(4)*10))//`) - '//rchar(s_bkgd,`f4.2')//`('//rchar(NINT(fit$e1(5)*100))//`)<times>('//ital//`x'//roman//`-'//ital//`x'//roman//`<_>0<^>)'

set %yloc 61
text `<rho,_,sigma>,area<^> = '//rchar(fit$corr(1,3)*100,`f5.1')//`%'
set %yloc 61-4
text `<rho,_,sigma>,bkgd<^> = '//rchar(fit$corr(3,4)*100,`f5.1')//`%'
set %yloc 61-4*2
text `<rho,_>area,bkgd<^> = '//rchar(fit$corr(1,4)*100,`f5.1')//`%'

fname = `fit-gauss+bkgd.eps'
ha s fname
alias\expand ttt `%/usr/local/bin/make-phys-eps-pretty.com '//fname
ttt
alias\expand ttt `%/usr/local/bin/eps2png '//fname
ttt