PRO df_demo
; Demonstrate the df_ft and df_showfit routines...
;
; The fitting routine, df_fit, is called as:
;    df_fit,ngauss,x,y,a,rms, sig,yfit
; and takes the following arguments:
;          df_fit,ngauss,x,y,a,rms,sig,yfit
;          ngauss = 1,2,3 for number of gaussians
;          x = x-vector
;          y = y-vector
;          a = parameters = height,ctr,disp,...,cont
;          rms = input rms of data (single value for all y,
;                 OR vector of length y)
;          sig = output array of parameter uncertainties
;          yfit = output fit y-vector
;
; The fit-plotting routine is df_show fit and is called
; with the results of df_fit:
;    df_showfit,x,y,a,rms, sig,yfit
;
; For reference the df_common values are:
; COMMON davefits, df_continuum, df_title, df_xtitle, df_ytitle, $
;		df_ghostH, df_ghostC, df_ghostS, df_ngauss
; which allows a "shaped continuum" and plot labels
; (and a ghost gaussian - but never mind!)
;
; The "continuum" or background, df_continuum, can be:
;  - a constant (e.g., 1.0) everywhere
;  - a fixed array of the same length as x and y
; in either case the fitting process just scales it
; by a scalar which is returned as the last ("cont") parameter of a.
; - - - - - -
;  6/17/98 dd
; - - - - - -

; To have the df_common variables accessible include this line
; in procedures:
@df_common

; Create a fake data set: wiggly continuum plus two
; gaussians
; x values, 0 to 10.0
x = FLOAT(indgen(101))/10.
; y values: wiggly continuum plus two
; gaussians
continuum_shape = 1.0 + SIN(2.*!PI*x/5.)
; To make the gaussians, use the gauss1 procedure
;       F = A(0)*EXP(-Z^2/2) + A(3)
;       Z = (X-A(1))/A(2)
; so a = [gauss_amplitude, gauss_center, gauss_sigma, constant]
;   centered at 4.0 with FWHM = 1.0
a = [1.0, 4.0, 1.0/2.35, 0.0]
df_continuum = 1.0
gauss1, x, a, g1
;   centered at 6.0 with FWHM = 5.0
a = [1.0, 6.0, 5.0/2.35, 0.0]
gauss1, x, a, g2
y = 30.0*continuum_shape + 100.*g1 + 20.*g2

; Plot the EXACT data before proceeding...
; two plots coming...
!p.multi = [0,1,2]
plot, x, y, PSYM=10, TITLE = 'Noise-free Input Data'

; Now add Poisson statistics (replacing each bin's
; value with a sample from a Poisson distribution
; with the bin's value as the mean (= mu)
for ip=0,n_elements(x)-1 do y(ip) = g_poiss(y(ip))

; Plot the Statistical ("real") data too...
plot, x, y, PSYM=10, TITLE = 'Data with Poisson Statistics'
!p.multi = 0

wait, 4

; - - - - - -
; Now we have fake (binned) data (x,y) and a
; continuum shape - let's fit it...


; - - - - with two gaussians and constant "continuum" - - 
ngauss = 2
; x = x
; y = y
; initial guess for parameters
; a is three params per gaussian plus the continuum level
; the gaussian params are: height, center, sigma
;
; add a little intelligence to the initial guesses:
a = [MAX(y), 5.0, 1.0, 0.5*MAX(y), 1.0, 1.0, TOTAL(y)/n_elements(y)]
; rms error for each bin
rms = SQRT(y)  ; counting statistics
; constant continuum
df_continuum = 1.0 + 0.0*x  ; tricky way to get all 1.'s
df_title = 'Two Gaussians plus Constant Continuum'
; Fit it
df_fit,ngauss,x,y,a,rms, sig,yfit
; Show the fit result
df_showfit,x,y,a,rms, sig,yfit

; look at desperate atempt for agreememt... 
; Reduced Chi-square around 5+
wait, 4.0

; - - - - with three gaussians and constant "continuum" - - 
ngauss = 3
; x = x
; y = y
; add a little intelligence to the initial guesses:
a = [MAX(y), 4.0, 1.0, 0.5*MAX(y), 6.0, 1.0, 0.5*MAX(y), 1.0, 1.0, TOTAL(y)/n_elements(y)]
; rms error for each bin
rms = SQRT(y)  ; counting statistics
; constant continuum
df_continuum = 1.0 + 0.0*x  ; tricky way to get all 1.'s
df_title = 'Three Gaussians plus Constant Continuum'
; Fit it
df_fit,ngauss,x,y,a,rms, sig,yfit
; Show the fit result
df_showfit,x,y,a,rms, sig,yfit

; looks better...
; Reduced Chi-square around 2.0
wait, 4.0

; - - - - OK, now the two gaussians with proper shaped "continuum" - - 
ngauss = 2
; x = x
; y = y
; add a little intelligence to the initial guesses:
a = [MAX(y), 4.0, 1.0, 0.5*MAX(y), 6.0, 1.0, TOTAL(y)/n_elements(y)]
; rms error for each bin
rms = SQRT(y)  ; counting statistics
; Shaped continuum!
df_continuum = continuum_shape
df_title = 'Two Gaussians plus Shaped Continuum'
; Fit it
df_fit,ngauss,x,y,a,rms, sig,yfit
; Show the fit result
df_showfit,x,y,a,rms, sig,yfit, CHISQ = reduced_chisq

; Note that the reduced chi-square was passed out by the
; keyword CHISQ
print, ''
print, ' Reduced Chi-square = ', reduced_chisq

; In case the gaussian AREAS are needed, they can be calculated
; from the height and sigma values (and height and sigma errors):
print, ' '
print, ' The Gaussians Areas w/errors:'
; Calculate the Gaussian areas from the parameters
; Scale factor to go from Gaussian height/width to Number of counts
; assuming fitting a histogram
; assumes X is equally spaced.
XperBin = x(1)-x(0)
scale=sqrt(2.*!pi)/XperBin
; go through the gaussian components
for i=0,n_elements(a)-2,3 do begin
  areac = a(i)*a(i+2)*scale
  sigc = areac * sqrt( (sig(i)/a(i))^2 + (sig(i+2)/a(i+2))^2 )
  print, 'Gaussian #'+STRING((i+3)/3,FORMAT='(I1)')+': ', areac, '  +/-', sigc
endfor
print, ''

RETURN
END