Fit 1-D profiles to a data cube
Fit 1-D profiles to a data cube
The routine can fit complex spectra with multiple components in a data cube (actually: fit along any
axis of a hyper-cube). It is multi-threaded and capable of fitting a large number of (i.e. order a million)
spectra in a few minutes, depending of course on the number of cores available. It borrows heavily
" routine of the GIPSY package.
The type of profiles that can be fitted and have been tested are
(gaussian with asymmetric wings), and
" (peaky gaussians possibly with asymmetric
wings). See the important
" Fitting Functions
" topic for more details.
The parameters for each fitted component reside in the SMURF_FIT1D extension as cube
NDFs called COMP_1,..., COMP_N. For a gaussian the planes in these cubes correspond
to: amplitude, position, and fwhm. Further details are under topic
" Fitted Parameters
A default config file is in $SMURF_DIR/smurf_fit1d.def. A sample file for user specified values is $SMURF_DIR/smurf_fit1d_uval.def.
"(case-insensitive) or a null (!) value is supplied, a set of default configuration parameter values will be used from $SMURF_DIR/smurf_fit1d.def. See the
"topic for detailed information.
"). Instead of a comma separated string of filenames a
"file can be submitted with each filename on a separate line. Files will map to components 1..N in the order as specified. See
"for more information on parameter files.
The full specification of a parameter file name consists of three parts: an optional directory path (
" ), an optional container ndf (
" ) and plus the base name for the file (
" For convenience
" can be specified through their respective parameters, but full name(s) can also be
" . In case
" is specified, FIT1D will append a
" is given the string
" will be appended.
Note that leaving
" blank will result in a conventional
" is specified, but
" is not, all components in the container file will be
read. Note that COMP_0 will be skipped since it hold diagnostics information about the
To escape special characters (
" , etc) in the string you may need to use a set of
quotes. A null value means no component parameter files will be used. [!]
", etc) in the string you may need to use a set of singledouble quotes. A null value results in use of the current directory. [!]
".MORE.SMURF_FIT1D.COMP_#. For further details see help on PARCOMP. To escape special characters (
", etc) in the string you may need to use a set of singledouble quotes. [!]
"#.par = value, or
"#.par = [0,1] with
’being a parameter relevant for the function being fitted.
’(1..n) indicates the component profile being described. Fix indicates a parameter to be kept fixed or fitted.
" #.fid will override the default function selected in the config file. Parameter
names are described in the help item
" Fitting Functions
" comp.fid comp.par Gauss: 1 a, b, c
Gausshermite1: 2 a, b, c, h3 Gausshermite2: 3 a, b, c, h4 Voigt: 4 a, b, c, l
" #.par parameter can have a value of: 1 = fix parameter at given value, do not fit -or- 0 = use
as initial estimate, but fit. As for comp, the
’ (1..n) indicates the component.
A null value for USERVAL means that no user-supplied estimates or fixed values are to be used. [!]
": fitting a single gausshermite2 to each profile along the highest dimension of the file
". The input file is expected to be baselined and the zerolevel of the profile to be 0. The fitted profiles are stored in the file
"with a component parameter file with the gauss-hermite parameters a, b, c, h3, h4 as
"(plus ...comp_0 for the diagnostics component).
’Same as above, but fit a single gaussian instead. Alternatively, use config=
’having the following lines: $SMURF_DIR/smurf_fit1d.def function = gaussian
"initial estimate from Fit1d may be in-accurate, in the first example a fit with a gausshermite2 will be attempted regardless. By contrast, the gaussian fit to a non-gaussian profile in example 2 may have been so poor that is was rejected: in that case current gausshermite2 fit will skip the profile because there won
’t be any initial estimates.
’Fit three gaussian components to each profile using initial estimates and fixed values as defined in the file
"), e.g.: comp1.b = -5.2 fix1.b = 1 comp1.c = 6 comp2.b = 20.4 fix2.b = 1 comp2.c = 4 comp3.b = 35.3 That is: provide a user-defined fixed value for the position (parameter
") of components 1 and 2, and an initial estimate for component 3. Also provide user-defined initial estimates for the FWHM (parameter
") of components 1 and 2. Leave it to fit1d to find initial estimates for all other parameters.
". See topic
"for details. If your aim is to capture a much emission as possible e.g. in order to create a 2-D image from a 3-D cube, gausshermite2 profiles are recommmended. [
") in a Voigt fit in terms of ==PIXELS==(!). If RETRY=1 has been set, a pure gaussian fit (L=0) will be attempted in case the initial fit of H3 is out of bounds. [undef]
’functions to fit to each profile, e.g. a multi-component spectrum of maximum three gaussians. [Note: The complete fit of the gaussians is done concurrently, not iteratively starting e.g. with the strongest component]. The routine will try to find and fit ncomp functions along each profile, but may settle for less. 
". Default is to use the full extent of the axis: [undef]
") or hermite parameters (
") are out of the routine re-tries the fit with the out-of-bounds parameter(s) fixed at 0. This means that in effect the fit cascades to a simpler function: gausshermite2 - gausshermite1 - gaussian; voigt - gaussian. The result is that there are valid fits for more profiles, but the function actually fitted may vary with position. Setting the retry value to 0 prevents this from happening and may cause the fit to fail or be (very) poor. [YES]
": sort by decreasing fitted value of the amp-like parameter
": sort by decreasing fitted fwhm of the width-like parameter
": sort by increasing position along the profile axis
": sort by increasing fitted distance from the centre pixel in the profile. Sorting can be helpful, but be cautioned that it can also complicate things: if there are two components one at -10 km/s and one at 10 km/s sorting by amplitude or width can result in the parameter file for component 1 to be a mix of the -10 and 10 km/s features depending on which one was relatively stronger or wider. Similarly, sorting by position can result in low-amplitude fits of noise spikes to be mixed with stronger components. For more precise control try to run the routine iteratively with e.g. a different restricted velocity range to try pick out the different components. Default is to sort by amplitude. [
’. Estimates can be very inaccurate plus are not checked against any boundary limits until after the fit. Thus this option may not be very helpful.
1) A standard GAUSSIAN. Parameters are a = maximum, b = centre, and c = FWHM.
NOTE that if one of h3 and h4 in a gauss-hermite function is non-zero, the mean of the distribution is not the position of the maximum (Reference; Marel, P. van der, Franx, M., A new method for the identification of non-gaussian line profiles in elliptical galaxies. A.J., 407 525-539, 1993 April 20):
2) GAUSS-HERMITE1 polynomial (h3). Parameters are a (amplitude), b (position),c (width), and h3 as mentioned these are NOT the same as maximum, centre, and fwhm of the distribution as for a gaussian: maximum = [determine value and position of max from fitted profiles using e.g. collapse] centre = b h3sqrt(3) FWHM = abs( c(1-3h32) ) = c skewness = 4sqrt(3)h3
3) GAUSS-HERMITE2 polynomial (h3, h4). Same as previous, but an extra parameter h4 is included: maximum = [determine value and position of max from fitted profiles using e.g. collapse] centre = b h3sqrt(3) FWHM = abs( c(1h4sqrt(6)) ) skewness = 4sqrt(3)h3 kurtosis = 8sqrt(6)h4
4) VOIGT function. Parameters are a (area), b (centre), c (doppler FWHM), l (lorenztian FWHM), and v (area factor) with relations: maximum = [determine value of max from fitted profiles using e.g. collapse] centre = b doppler fwhm = c lorentzian fwhm = l amp = v (OUTPUT ONLY!) amplitude calculated from a (area) using the standard amp2area function for a voigt (based on the Faddeeva complex error function): amp = area / amp2area
’area 2 = position
’position 3 = fwhm
’4 = -
’5 = -
’last: function id 1 = gaussian; 2 = gausshermite1 (h3); 3 = gausshermite2 (h3,h4), 4 = voigt
1 Too many free parameters, maximum is 32.
2 No free parameters.
3 Not enough degrees of freedom.
4 Maximum number of iterations too small to obtain a solution which satisfies TOL.
5 Diagonal of matrix contains elements which are zero.
6 Determinant of the coefficient matrix is zero.
7 Square root of negative number. -10 All fitted components rejected due to minamp, minwidth, maxlorz, or range constraints.