A Manual for Running Gausfit, and Tempfit ========================================== Gausfit ----------- This program produces an analytic fit to emission line blends in AGN. The strongest line in the blend is fit with a sum of three Gaussians, weaker blended components are each fit with a single Gaussian. A detailed description of the algorithm used in the code is given in the Appendix section of Laor et al. (1994, ApJ 420, 110) The input spectrum is read from unit 4. The first line has 3 numbers: 1. number of data points, 2. redshift, 3. Galactic E(B-V) The second and following lines have the spectrum. Each line has: 1. Observed wavelength in A, 2. Observed flux in cgs, 3. Flux error in cgs. Note that for z<0.01 only two numbers are read from each line. An arbitrary S/N is assumed. This option is used for optical spectra (where I always used rest wavelength), where the S/N is usually not specified. The line fit parameters are read from unit 3. The first line has 6 numbers: 1. wline - rest wavelength of the strong line. 2. ma - number of fit parameters (equals 3*i where i= number of Gaussians) 3, 4. wcon1, wcon2 - wavelength range where the fit is made. Continuum level is set using the median flux at these two points 5. iblue =1 use blue half only, =0 use red half only, =2 use all the line, when calculating the Gaussian parameters initial guess 6. nfix - number of parameters to be held fixed (see below). The second line has ma number: lista() - a list of the parameter numbers: 1 = area of first Gaussian, 2 = center of first Gaussian, 3 = width of first Gaussian, 4,5,6 = as above for second Gaussian, etc... The last nfix parameters in the list are to be held fixed. The third line has nfix numbers: a() = the values of the parameters which are held fixed. The first value is read into LAST parameter in lista(). The parameters are in the following units: 1 - flux in cgs, 2 - lcenter in A, 3 - sigma (FWHM/2.3548) in A. Note that this line should exist even if nfix=0 (due to FORTRAN peculiarities). The numbers in this line are unimportant when nfix=0. The fourth line has (ma-9)/3 numbers (i.e. the Gaussians in addition to the first three. wl() = The rest wavelengths of the other components, besides the main line at wline. This line is read only if ma>9. EXAMPLE ---------- Enclosed with this message you should get a file with the spectrum of PKS 0405-123 which you can use to try out the code. Copy it to unit 4 and try the following for unit 3: 1215.67 18 1135. 1370. 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 5. 1240.15 2.e-13 1240.15 1303.49 1335.31 After running the code the output on your screen should be: ------------------------------------------------------- type in number of objects 1 <---------- You type this in 1215.67 18 1135.0 1370.0 1 0 1240.15 1303.489 1335.31 5625 data points z = .5727999 E(B-V) = 3.10000E-02 data smoothed with sigma=1.75pix Gaussian continuum windows: 1135.0 2.00803E-14 5.72006E-16 1370.0 1.62955E-14 3.07465E-16 type in iadd (-1 0 1), the continuum displacement 0 <---------- You type this in alpha_nu = -.890141 846 points left after median rejection ch^2= 1577.949 after 16 iterations lambda_line f_cont f_line EW v shift FWHM 1215.67 1.86E-14 1.35E-12 46.2 1562 15036 1215.67 1.86E-14 8.83E-13 30.2 -149 4325 1215.67 1.86E-14 3.47E-13 11.8 -211 1264 1240.15 1.82E-14 7.94E-14 2.8 -188 2101 1303.49 1.72E-14 7.21E-14 2.7 168 4518 1335.31 1.68E-14 9.72E-15 .4 -854 2276 parameters of fitted model 1214.84 1.86E-14 2.74E-12 93.7 -204 2285 ------------------------------------------------------------- In this example all parameters are free. Note that some of the numbers above are slightly different than those in Laor et al. (1994, ApJ 420, 110) where the same parameters were used. This results from the small machine dependency of the code. The earlier calculations were made on a SUN sparcstation IPC, and the current one on an HP 712/60. The parameters of the individual Gaussians are very sensitive to the initial guess, and also to the exact way each machine does the calculations. This demonstrates why the individual components parameters are not very useful. The parameters of the sum of the Gaussians are obviously much more robust since they are set by the observed spectrum. Note that one can fit more than one object at the same time. All the object spectra should be in unit 4. To get some estimate of the errors, repeat the fit with the same parameters, but respond to the input request with iadd=-1 and +1. You can get a quick view of the fit with the following SM file: -------------- lweight 2 data ftn14 read x 1 read y 2 set y=lg(y) limits x y connect x y data ftn11 read x 1 read y 3 set y=lg(y) lweight 1 histogram x y data ftn13 lweight 2 read x 1 read y 2 set y=lg(y) connect x y data ftn21 read x 1 read y 2 set y=lg(y) ptype 3 3 points x y box expand 1.1 xlabel \gl ylabel F\d\gl lweight 1 --------------- The unit names (e.g. ftn14) are the HP convention. The filled triangles (from ftn21) represent data points not included in the fit. There are other output files, their content is explained in comments at the beginning of the code. Another example for unit 3: 1215.67 18 1135. 1370. 1 3 1 2 3 4 5 6 7 8 9 13 14 15 1617 18 10 11 12 5. 1240.15 2.e-13 1240.15 1303.49 1335.31 In this case all three parameters of the Gaussian describing N V are fixed a priori. The following contains examples of unit 3 files for fitting other lines: C IV -------- 1549.05 18 1460. 1720. 2 1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 17 1486.5 1640.46 1664.15 1486.50 O I ------- 1303.5 6 1290. 1322. 2 0 1 2 3 4 5 6 1303.5 C III] -------- 1908.73 15 1830. 1980. 0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1892.03 1892.03 1857.40 Mg II -------- 2798.74 9 2690. 2910. 2 1 1 2 3 4 5 6 7 8 9 4.75 N III] -------- 1750.46 3 1725. 1767. 2 0 1 2 3 1750.46 O VI -------- 1033.816 15 950. 1110. 2 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 977.02 990.98 977.02 990.98 O IV] + Si IV -------------- 1399.605 6 1360. 1440. 2 0 1 3 4 6 5 2 1402.46 1396.75 Hbeta -------- 4861.32 12 4600. 5120. 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 4686. 4686. 4959. 5007. [O III] --------- 5006.8 6 4975. 5030 2 0 1 2 3 4 5 6 7 8 9 5006.8 ================================================================= ================================================================= Tempfit --------- This program fits a template line shape to a given emission blend. This template deblending method is likely to give realistic values for the individual components when an acceptable fit is obtained. Each template component has two free parameters: 1. flux, 2. velocity shift. Lines which are doublets are fit with two template components. The velocity difference of the two components is fixed at the doublet separation, in their flux ratio is fixed at 1:1 or 2:1. For more details see Laor et al. (1994), and Laor et al. (1995, ApJS July issue). Some of the required input files are generated by running Gausfit. These are: unit 26: The observed profile, including only data points used in the fit. unit 27: The observed profile, including all data points. unit 28: The profile which is used as a template. You should make the following file unit 33: The fit parameters: The first line has 4 numbers: 1. wline - the wavelength used for zero velocity. 2. ifix - a parameter which describes which blend is being fit. 3. ma - number of fit parameters (equals 2*i where i= number of Gaussians) 4. mfit - number of free parameters (to be found in the fit) The Second line has a list of ma numbers. Their order has the same meaning as in unit 3 input to Gausfit The next ma/2 lines have the initial values for the fit parameters. The units are: 1. line flux in cgs 2. velocity shift from wline in km/s. Next line, 1 or 2 numbers: Velocity separation between the doublet components Next line, 2 numbers: Velocity range over which chi^2 minimization is made Next line, 1 number: The objects redshift The major difference between this program and Gausfit is that here you have to specify by hand the initial guess. In many cases the code doesn't wander too much off the initial guess in parameter space, in particular if there are more than 6 free parameters (3 components). It is useful to inspect the result by eye and try a different initial guess if the fit doesn't look very good. This puts in a subjective element into the code. The code can surely be improved so that an objective, well defined, initial guess is made. Unfortunately, I didn't have time here to go beyond a code which just "works". Example: Fit to Ly alpha + N V + Si II --------------------------------------- Run Gausfit with the example given above for Ly alpha. This generates units 26, 27, and 28. Try the following for unit 33: 1215.67 2 8 6 1 3 5 7 2 4 6 8 5.20000E-13 0. 5.20000E-14 5713.00 <---- the 1238.821A component of N V 2.60000E-14 6696.00 <---- the 1242.804A component of N V 1.60000E-14 11757.0 <---- the 1263.313A Si II multiplet 983.000 <---- Velocity separation of N V -15000.0 15000.0 0.5728 After running the code the output on your screen should be: ------------------------------------------------------------ assume optically thin case (2:1) 1215.67 2 8 6 1 3 5 7 2 4 6 8 5.200E-13 .00 5.200E-14 5713.00 2.600E-14 6696.00 1.600E-14 11757.00 983.0 z= .5727999 422 -15000.0 15000.0 <-- # of data points + velocity range for fit 3 903.1705 <-- # of iteration + chi^2 5 903.1705 type in # of additional iterations 5 <------------- You type this in 3 903.1705 5 903.1705 type in # of additional iterations 0 <------------- You type this in FWHM= 2250.12 FWHM Ly alpha= 2214.443 2.27482E-12 -163.274 2.00336E-13 5534.252 1.00168E-13 6517.252 6.07335E-14 11757.0 ------------------------------------------------ The last 4 lines give the best fit parameters. Note that the velocity of Si II was not a free parameter (remained at 11757). Try the following SM file to view the fit (note that the N V line is a 2:1 sum of its components). data ftn02 read x 1 read y 2 limits 1140 1290 0 9e-14 box connect x y read y1 3 connect x y1 read y2 4 connect x y2 read y3 5 connect x y3 read y4 6 connect x y4 Another Example: Fit to C IV + He II + O III] + N IV] ----------------------------------------------------- 1. Save unit 28 in another file, say pks0405.template. 2. Run Gausfit in order to generate units 26 and 27 for the C IV blend. Try the following as an input for unit 3 1549.05 15 1460. 1720. 2 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1486.5 1640.46 1664.15 Look at the fit with SM to make sure it is OK. 3. Copy pks0405.template back to unit 28 4. Run Tempfit. Try the following for ftn33 1549.05 1 8 8 1 3 5 7 2 4 6 8 2.00000E-14 -165.97 <--- The 1548.195A component of C IV 1.00000E-14 332.7 <--- The 1550.770A component of C IV 1.00000E-15 17702.7 <--- He II 1640.46 1.00000E-15 22290. <--- O III] multiplet at 1664A 498.67 -15000.0 25000.00 0.5728 You should get the following output: ---------------------------------------- assume optically thin case (2:1) 1549.05 1 8 8 1 3 5 7 2 4 6 8 2.000E-14 -165.97 1.000E-14 332.70 1.000E-15 17702.70 1.000E-15 22290.00 498.67 z= .5727999 536 -15000.0 25000.0 3 15104.47 5 15104.47 type in # of additional iterations 0 FWHM= 2453.672 FWHM Ly alpha= 2213.845 7.60456E-13 -441.688 3.80228E-13 56.98199 9.86568E-14 16377.66 8.41043E-14 22098.88 -------------------------------------- Note the very large chi^2 (15104) for 528 (=536-8) degrees of freedom. The fit for C IV is clearly not good, as you can check by plotting the results in unit 2. The following contains examples of unit 33 files for fitting other lines: Mg II ------- 2798.74 1 4 3 1 3 2 4 1.e-14 -256.56 5e-15 513.13 769.7 -10000 10000 put here redshift O VI ----- 1033.816 1 12 6 1 3 5 7 9 11 2 4 10 8 6 12 2.00000E-14 -548.400 <--- The O VI component at 1031.926A 1.00000E-14 1102.90 <--- The O VI component at 1037.617A 1.00000E-14 -2348.5 <--- Ly beta at 1025.72A 1.00000E-15 -12430.0 <--- N III at 990.98A 1.00000E-15 -16481.0 <--- C III at 977.02A 1.00000E-15 -17782.0 <--- Ly gamma at 972.53A 1651.30 -30000.0 3000.00 put here redshift C III] ------- 1908.73 2 8 5 1 3 5 7 2 4 6 8 5.20000E-14 0. <--- C III] at 1908.73A 5.20000E-15 -8489.5 <--- Al III at 1854.716A 2.60000E-15 -7220.5 <--- Al III at 1862.789A 5.E-15 -2624.8 <--- Si III] at 1892.03A 1269.000 -15000.0 15000.0 put here redshift Si IV + O IV] ---------------- 1399.605 5 14 2 1 7 2 8 3 4 5 6 9 10 11 12 13 14 2.e-16 2.e-17 -1253.9 --- most relevant parameters are in the code 1291.4 -5000 5000 put here redshift