! Anything after an exclamation mark is a comment
!
! Read in the formatted data into vectors x and y
read 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 backgroun

! 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
! But first, define pi = 3.14...
pi = acos(-1)
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\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

! set the colour to black
set plotsymbolcolor 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
  %xlaxis 20                    ! the lower x-axis, % of the screen
  %xuaxis 95                    ! the upper x-axis, % of the screen
  %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
  linthk 3                      ! the rest are just what I think is pretty
  xlabsz 0.35
  ylabsz 0.35
  xnumsz 0.35
  ynumsz 0.35
  %xtics .5
  %ytics .5
  %xticl 1.25
  %yticl 1.25
  nsyinc 5
  nsxinc 5
  font triumf.2
  cursor -9
  yleadz 1
  xleadz 1
  linewidth 1

! clear the x-axis label (we'll be plotting residuals below)
set xlabel ` '
! and the y axis
set ylabel `Number of counts'
! set the scales, xlow xhi ylow yhi (0 0 means autoscale)
scales 758 792 0 1600

! Graph the data with error bars
gr x y sqrt(y)

! change to a solid line
set pchar 0
! make the colour blue
set curvecolor 4
! without replotting the scales, graph the background
gr\noax x ybkgd

! make the colour red
set curvecolor 2
! overlay the total fit (including background)
gr\noax x yfit

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

! -----------------------------------------------
! Next we plot the residuals below the main graph
! -----------------------------------------------


! some definitions of italic and roman fonts used later on
ital  = `<fitalic.2>'   ! stuff in quotes is a text variable; 
roman = `<ftriumf.2>'  ! gets evaluated when used


! 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 758 792 -3 2.999
! now label the x axis, and change y's label
set xlabel `Channel number'
set ylabel `<v0.6>(y-y<_>fit<^>)/<sigma>'
! graph the residuals
gr x yresid
set curvecolor 2 ! red
gr\noax x 0*x ! graph a zeroline
set curvecolor 1 ! black
set
  txthit 0.2  ! height of text
  %xloc 22.5     ! location of text
  %yloc 85
  cursor -7    ! left-aligned

! Spit out fit result information
text `x<_>0<^> = '//rchar(b,`%8.4f')//`('//rchar(NINT(fit$e1(2)*1e4))//`)'
set %yloc 81
text `<sigma> = '//rchar(c,`%6.4f')//`('//rchar(NINT(fit$e1(3)*1e4))//`)'
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,`%6.4f')
set %yloc 67
text `C.L. = '//rchar(fit$cl,`%4.1f')//`%'
set %yloc 61
text `<rho,_,sigma>,area<^> = '//rchar(fit$corr(1,3)*100,`%5.1f')//`%'
set %yloc 61-4
text `<rho,_,sigma>,bkgd<^> = '//rchar(fit$corr(3,4)*100,`%5.1f')//`%'
set %yloc 61-4*2
text `<rho,_>area,bkgd<^> = '//rchar(fit$corr(1,4)*100,`%5.1f')//`%'

set %yloc 33
set %xloc 30
text `y<_>bkgd<^> = '//rchar(bkgd,`%5.1f')//`('//rchar(NINT(fit$e1(4)*10))//`) - '//rchar(s_bkgd,`%4.2f')//`('//rchar(NINT(fit$e1(5)*100))//`)<times>(x-x<_>0<^>)'
hardcopy\png fit-gauss+bkgd.png