;+ ;Project: ; SDAC ;Name: ; RUN_CURVEFIT ; ;PURPOSE: ; This procedure interfaces the least squares fit procedure with ; the data and response functions for SPEX. ; ;Category: ; SPEX, ANALYSIS, FITTING ; ;Explanation: ; This procedure controls the input to the functional minimization ; program, usually CURVFIT. In the standard technique there is ; an inner and outer loop. In the outer loop the conversion ; factors are recompute with the new values of the parameters ; obtained from CURVFIT. The inner loop counter is used within ; CURVFIT. ; If there are no degrees of freedom for two free parameters and two ; data values, and one of the free parameters is a normalization parameter, ; while the other is a shape parameter directly after it, ; the procedure, SPECTRAL_RATIO, is used to obtain the "fitted" values. ; ; ;Keyword Inputs: ; E_IN - Mean of m input energy bins, geometric mean is best normally, or 2Xm edges ; W_IN - Width of input energy bins. ; COUNT_2_FLUX - conversion vector from incident counts to count flux ; for hard x-ray data normally the product of energy bin width and area ; LIE - Livetime of observing interval in seconds (or unit of choice) ; DRM - Detector response matrix dimensioned n x m where n is the number ; of output energy bins and m is the number of input energy bins. ; F_MODEL - Name of functional model (string) ; OBSERVE (OBS) - Source counting rate spectrum ; EOBSERVE (EOBS) - Uncertainties on OBS calculated from Gaussian stat. on observed ; USE_MOD_SIG - If set, use Gaussian statistics on expected rate ; BACKGROUND (BACK) - Backgrount counting rate spectrum ; EBACK (EBACK) - Uncertainties on BCK ; UNCERT_MIN - minimum relative uncertainty used to calculate weighting ; E_OUT - Mean of n output energy bins, or 2XN edges. ; W_USE - Indices to use in E_OUT and OBS ; APAR - Parameters used in call to F_MODEL. ; FREE - set (1) if parameter is free, reset (0) if it is fixed. ; RANGE - Min and Max values allowed for each of APAR, ; dimensioned 2,n where n is the number of parameters ; ITER - Number of times in outside loop where conversion factors are ; recomputed. ; If ITER is set to 0, then fitting is done in count rate space, ; and there are no conversion factors used. ; NITER - Number of times in loop inside Curvefit where parameter ; increments are calculated ; QUIET - If set then informational print statements are suppressed. ; NOLAMBDA - see curvfit ; PCONV - pileup conversion factor if set, the ratio of a single photon pulse-ht ; spectrum divided by the pileup pulse height spectrum ; These corrections can be obtained using [sdac.spex]convolve.pro ; NOFIT - Make all computations but don't run curvefit if set. ;KEYWORD Outputs: ; APAR - Input parameter values are adjusted to new values in CURVEFIT ; CHISQR - Calculated CHISQR based on OUTPUT APAR and input values ; SIGMAA - Uncertainties on the APAR from CURVEFIT, somewhat crude ; CONVERSION - model count spectrum divided by model photon spectrum ; NRESID - Residuals on fitted points normalized by uncertainties used, same ; length as wuse ;USAGE: run_curvefit, e_in=e_in, w_in=w_in, drm=drm, f_model=f_model, $ ; e_out=e_out, live=live, count_2_flux=count_2_flux, wuse=wuse, observe=obs, $ ; uncert_min=uncert_min, eobserve=eobs, background=bck, eback=ebck, $ ; use_mod_sig=use_mod_sig, iter=iter, niter=niter, quiet=quiet, free=free, range=range, $ ; chisqr=chisqr, apar=apar, sigmaa=sigmaa, conversion=conversion, $ ; nolambda=nolambda ;CALLS: ; GET_CONVERSION, PRINTX, CHECKVAR, F_DIV, FCOUNT_RATE ;History: ; Written by RAS 92/9/26, ; Mod 11-feb-94 to use model uncertainties ; 28-oct-94, added pconv keyword to correct for pileup effects ; 23-mar-95, ras, e_out can be mid-points or edges (2xm) inside curvfit ; 23-jan-96, added nresid output and nofit keyword ; Version 6, richard.schwartz@gsfc.nasa.gov, 18-aug-1997, added section ; to fit to observed counts, not photons if iter is set to le 0. ; Version 7, richard.schwartz@gsfc.nasa.gov, 7-oct-1997, added keyword nofit ; Version 8, richard.schwartz@gsfc.nasa.gov, 28-jan-1998 added interface to spectral_ratio(). ; 5-may-2003, richard.schwartz@gsfc.nasa.gov, prevent crate_sig crash with iter=0 ;- pro run_curvefit, e_in=e_in, w_in=w_in, drm=drm, f_model=f_model, e_out=e_out, $ live=live, count_2_flux=count_2_flux, wuse=wuse, observe=obs, uncert_min=uncert_min,$ eobserve=eobs, background=back, eback=eback, use_mod_sig=use_mod_sig, $ iter=iter, niter=niter, quiet=quiet, free=free, range=range, $ chisqr=chisqr, apar=apar, sigmaa=sigmaa, conversion=conversion, pconv=pconv, $ nolambda=nolambda, nresid=nresid, nofit=nofit PRINTX, 'RUN_CURVEFIT printx, 'Input parameters for '+f_model+' on entering loop: ' printx, apar checkvar, iter, 9 ;External loop where response function is used. ; ; If there are no degrees of freedom, try to obtain parameters using a ratio solution. ; norm_par = model_components( f_model ) wfree = where( free, nfree ) wnorm = where( wfree(0) eq norm_par, first_par_norm) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; If there are only two free parameters and no degrees of freedom, and the ; first free parameter normalizes the function with the non-linear parameter, ; then interpolate on array of channel ratios to obtain the non-linear parameters value. ; Obtaining that, reconcile the normalization to one of the observed count-rates. ; if n_elements(wuse) eq 2 $ ; two data points and nfree eq 2 $ ; two free parameters and total(free(norm_par)) eq 1 $ ; 1 and only 1 of free parameters is a normalization and (wfree(1 < (n_elements(wfree)-1)) - wfree(0)) eq 1 and first_par_norm eq 1 then begin ;first free is normalization. eflux = sqrt( eobs^2 + eback^2 ) apar_in =(1.+fltarr(2,1))#apar ;*** apar_in(*, wfree(1)) = range(wfree(1),*) ;*** ;*** apar(wfree) = spectral_ratio( obs-back, ecflux=eflux, epar=sigmaa, $ ;*** drm=drm, model=f_model, photon_edges=e_in, wuse=wuse, apar_in = apar_in) ;*** chisqr = 0.0 ;*** nresid = obs * 0.0 ;*** get_conversion, e_in=e_in, w_in=w_in, drm=drm, e_out=e_out, $ ;*** f_model=f_model, apar=apar, model_out=model_out, conversion=conversion ;*** conversion = f_div( conversion, fcheck(pconv, 1.0) ) ;*** ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; endif else begin if (size(e_out))(0) eq 1 then x = e_out(wuse) else x = e_out(*,wuse) ;begin ;independent variable if iter gt 0 then begin for i=0,iter+1 do begin get_conversion, e_in=e_in, w_in=w_in, drm=drm, e_out=e_out, $ f_model=f_model, apar=apar, model_out=model_out, conversion=conversion conversion = f_div( conversion, fcheck(pconv, 1.0) ) y = f_div( obs(wuse)-back(wuse),conversion(wuse) ) use_mod_sig = fcheck(use_mod_sig, 0) ; ras, temp, 5-jun-02 ,1) if use_mod_sig then begin model_crate = model_out * conversion + back crate_sig = f_div( model_crate, live *count_2_flux) wnz = where(crate_sig ge 0., nnz) if nnz ge 1 then crate_sig(wnz) = crate_sig(wnz)^0.5 eobs_2use= crate_sig endif else begin eobs_2use = eobs wzero = where(eobs eq 0.0,nzero) ;June 5 2002, ras, use 1 normalized to the time, area, and energy if nzero ge 1 then eobs_2use(wzero) = 1./(live*count_2_flux)[wzero] ;crate_sig(wzero) endelse ey = f_div( (sqrt( eobs_2use^2 + eback^2))(wuse), conversion(wuse) ) ey = ey > (fcheck(uncert_min,0.0) * (y * (1-use_mod_sig) + use_mod_sig*model_out(wuse)) ) wt = 1/ey^2 ;USE CURVFIT.PRO TO AVOID CONFUSION WITH USERLIB VERSION CURVEFIT if i lt (iter+1) and not fcheck(NOFIT, 0) then $ yfit=curvfit( x,y,wt,apar,sigmaa,niter=niter, funct=f_model, $ quiet=quiet,free=free,range=range,chisqr=chisqr,$ nolambda=nolambda) endfor nresid = f_div( y-model_out(wuse), ey ) chisqr = total((y-model_out(wuse))^2*wt) / (n_elements(wuse)-total(free)) ; print,'Using new parameters and conversion factors CHISQR= ',CHISQR endif else begin ; ; Here all the iteration is done within CURVFIT. The fitting is done on the ; actual counts, and not the modified (divided by conversion factor) counts. ; The fear here is that this will be substantially slower. ; for i=0,2 do begin y = obs(wuse)-back(wuse) y_mod = fcount_rate( e_out, apar, e_in, drm, wuse=wuse ) use_mod_sig = fcheck(use_mod_sig,0) if use_mod_sig then begin model_crate = y_mod + back crate_sig = f_div( model_crate, live *count_2_flux) wnz = where(crate_sig ge 0., nnz) if nnz ge 1 then crate_sig(wnz) = crate_sig(wnz)^0.5 eobs_2use= crate_sig endif else begin eobs_2use = eobs wzero = where(eobs eq 0.0,nzero) if nzero ge 1 then begin model_crate = y_mod + back crate_sig = f_div( model_crate, live *count_2_flux) eobs_2use(wzero) = crate_sig(wzero) endif endelse ey = (sqrt( eobs_2use^2 + eback^2))(wuse) ey = ey > (fcheck(uncert_min,0.0) * (y * (1-use_mod_sig) + use_mod_sig*y_mod(wuse)) ) wt = 1/ey^2 ;USE CURVFIT.PRO TO AVOID CONFUSION WITH USERLIB VERSION CURVEFIT if i lt 2 and not fcheck(NOFIT, 0) then $ yfit=curvfit( x,y,wt,apar,sigmaa,niter=niter, funct='fcount_rate', $ quiet=quiet,free=free,range=range,chisqr=chisqr,$ nolambda=nolambda) endfor y_mod = fcount_rate( e_out, apar) nresid = f_div( y-y_mod, ey ) ;wuse, e_in, and drm are still set. chisqr = total((y- y_mod)^2*wt) / (n_elements(wuse)-total(free)) ; print,'Using new parameters and conversion factors CHISQR= ',CHISQR get_conversion, e_in=e_in, w_in=w_in, drm=drm, e_out=e_out, $ f_model=f_model, apar=apar, model_out=model_out, conversion=conversion conversion = f_div( conversion, fcheck(pconv, 1.0) ) endelse endelse end