; $Id: minuitfit.pro,v 1.2 2002/01/14 20:59:26 amunter Exp $
FUNCTION minuitfit,x,y,w,sigma=sigma,minimizer=minimize
;+
; NAME:
;       MINUITFIT
; PURPOSE:
;       fit y=f(x) where:
;       f(x) = a0+a1*x+a2*x^2+arbitrary number of gaussians+lorentzians
; CATEGORY:
;       ?? - fitting
; CALLING SEQUENCE:
;       yfit = minuitfit(x,y,w)
; INPUTS:
;       x = independent variable, must be a vector.
;       y = dependent variable, must have the same number of points
;               as x.
;       quiet = set to inhibit printing curfit iterations.
; OUTPUTS:
;       yfit = fitted function.
; OPTIONAL OUTPUT PARAMETERS:
;       a = coefficients. a three element vector as described above.
;
; COMMON BLOCKS:
;       None.
; SIDE EFFECTS:
;       
; RESTRICTIONS:
;       Need to have the minlink.sl library in your $IDL_PATH
; PROCEDURE:
;       Uses MINUIT, the fortran callable CERN library to minimize
;       Chi-squared.
; MODIFICATION HISTORY:
;       Paul Schleger, Riso National Lab. Oct 1993
;-
;
common fitp,ngau,nlor,a,vary,pnames,gw,ufunct,nufun
COMMON fitinfo_com, fit
COMMON com_minchoice, minchoice, minpath
;The dimension of the convolution, and more
COMMON com_convolution, condim, epar, resoly, tilt, ghorder  
COMMON preferences, pref

IF NOT keyword_set(minimizer) THEN minimizer='MFit'
;Check for valid minimizer choices
IF (minimizer NE 'MFit') AND (minimizer NE 'MFit+') AND (minimizer NE 'Simp') THEN $
        minimizer='MFit' 
n=n_elements(y)
w=replicate(0.,n)+w
index=where(vary,count)
IF NOT keyword_set(sigma) THEN sigma=fltarr(n_elements(a))
IF (nufun NE 0) THEN BEGIN
        message, 'Cannot do a minuit fit on user functions. Returning....'
        return, 0.0
END
;Initialize MINUIT
llll=get_kbrd(1)
result=CALL_EXTERNAL(minpath,'_mninit_')
;print,!ERROR


;Set up the parameters
FOR i=0,n_elements(a)-1 DO BEGIN
        name=getpname(i,ngau,nlor,pnames)
        IF (STRPOS(name,'I') NE -1) THEN upperbound=10.0*max(y) ELSE $
        IF (STRPOS(name,'W') NE -1) THEN upperbound=10.0*(max(x)-min(x)) ELSE upperbound=0.0
        ;breakpoint
        result=CALL_EXTERNAL(minpath,'_mnparm_',i+1,name,$
     
        
        a(i),(ABS(a(i))+1.0)/100.0,0.0,FLOAT(upperbound))
 
END
chisquared=0.0                                  ;Must initialize chisquared to a float before 

;Fix the fixed parameters
fixem=WHERE(vary NE 1,count)                    ;Our first par number is 0
IF count NE 0 THEN BEGIN
 fixem=fltarr(n_elements(fixem))
 fixem(*)=FLOAT(WHERE(vary NE 1)+1.0)           ;Minuit says that the first par is 1
 result=CALL_EXTERNAL(minpath,'_mnexcm_','FIX',fixem,FIX(n_elements(fixem)))
END

;Send the data, and relevant parameters to define the fitting function
result=CALL_EXTERNAL(minpath,'send_pars',x,y,w,FIX(n_elements(x)),$
        ngau,nlor,gw,a,chisquared)

;Do the 2D convolution if condim == 2
IF (pref.vert NE 0) THEN BEGIN
result=CALL_EXTERNAL(minpath,'do_2d_convolution',FIX(ghorder),resoly,tilt,epar)
END
        
;Set the MINUIT printout level to a minimum of stuff
result=CALL_EXTERNAL(minpath,'_mnexcm_','SET PRI',[0.0],1)        

;Minimize
IF minimizer EQ 'Simp' THEN $
        result=CALL_EXTERNAL(minpath,'_mnexcm_','SIMPLEX',FLTARR(1),0)     
IF (minimizer EQ 'MFit') OR (minimizer EQ 'MFit+') THEN $
         result=CALL_EXTERNAL(minpath,'_mnexcm_','MINIMIZE',FLTARR(1),0)

;Get the parameter values
result=CALL_EXTERNAL(minpath,'get_pars',a,sigma,FIX(n_elements(a)),chisquared)
;Note that the fitting is done on Chi-squared, NOT the reduced chi-squared, so we've got to 
;divide chisquared by the number of degrees of freedom, to make it close to 1
chisquared = chisquared/(n_elements(x)-TOTAL(vary))

;Now do an error analysis
IF (minimizer EQ 'MFit+') THEN $
        result=CALL_EXTERNAL(minpath,'_mnexcm_','MINOS',FLTARR(1),0)

;redefine the fit information structure
IF TOTAL(vary) GT 0 THEN covar=FLTARR(TOTAL(vary),TOTAL(vary)) ELSE covar = 0.
fit={a:a,eplus:FLTARR(n_elements(a)),eminus:FLTARR(n_elements(a)),eparab:FLTARR(n_elements(a)),$
        globcc:FLTARR(n_elements(a)),covar:covar,chisquared:chisquared,$
        fmin:0.0,fedm:0.0,errdef:0.0,istat:0,fitstat:' ',npari:0,nparx:0}
        
;Get the parameter errors.
;NOTE the stupid way we must get the values. This seems to be the only way to reliably
;pass parameters _FROM_ a C program _TO_ IDL; i.e. don't use structures or array elements.
;Just use either a simple INT, DBL or FLT array, or a simple variable
FOR i=0,n_elements(a)-1 DO BEGIN
        eplus=0.0 & eminus=0.0 & eparab=0.0 & globcc=0.0
        result=CALL_EXTERNAL(minpath,'_mnerrs_',FIX(i+1),eplus,eminus,$
                eparab,globcc)
        fit.eplus(i)=eplus & fit.eminus(i)=eminus & fit.eparab(i)=eparab & fit.globcc(i)=globcc
;If the user wanted an MINOS analysis then say that the sigmas are the average of the
;positive and negative MINOS errors
        IF (minimizer EQ 'MFit+') THEN sigma(i) = 0.5*(ABS(eplus)+ABS(eminus))
END

                
;Get the covariance matrix. Again in a stupid way.
result=CALL_EXTERNAL(minpath,'_mnemat_',covar,FIX(TOTAL(vary)));
fit.covar=covar

;Get the status of the fit
fmin=0.0 & fedm=0.0 & errdef = 0.0 & npari=0 & nparx=0 & istat=0
result=CALL_EXTERNAL(minpath,'_mnstat_',fmin,fedm,errdef,npari,nparx,istat)
fit.fmin=fmin & fit.fedm=fedm & fit.errdef=errdef 
fit.npari=npari & fit.nparx=nparx & fit.istat=istat
CASE fit.istat OF
        0: fit.fitstat='Covariance Matrix: Not calculated'
        1: fit.fitstat='Covariance Matrix: Diagonal approximation only, not accurate'
        2: fit.fitstat='Covariance Matrix: Full matrix, but forced positive definite'
        3: fit.fitstat='Covariance Matrix: Full accurate matrix. Indicates normal convergence'
END
;Send the data, and relevant parameters to define the fitting function, one more time
result=CALL_EXTERNAL(minpath,'send_pars',x,y,w,FIX(n_elements(x)),$
        ngau,nlor,gw,a,chisquared)
RETURN, chisquared
END
