From MediaWiki

Jump to: navigation, search



I don't guarantee this is good or entirely correct. I am just a beginner in all this. But it can be a good starting point.

Analysis of galaxy images

Three Phases

  • Phase 1: SExtractor
  • Phase 2: GALPHOT isophotal analysis
  • Phase 3: 2D models (GALFIT or GIM2D or BUDDA), fixed centers
    • Pure Sersic: I, B, Re, nSer
    • Sersic+Expon: Ie, Re, nSer, μ0, h
    • Sersic+Expon+NuclearComp: Ie, Re, nSer, μ0, h, Inuc, Bnuc

MGC: 2-D light profiles for all 10095 galaxies modelled using GIM2D

  • Sersic -- single component
  • De Vaucouleurs Bulge + Exponential Disk
  • Sersic Bulge + Exponential Bulge

Theory on fitting

  • Paper: Shade et al. 1996, ApJ 464, 79

The bulge component is represented by a de Vaucouleurs R1/4 law:

 I_B(r_B)=I_B(0)\ exp\left[-7.67\left(\frac{r_B}{r_e}\right)^{0.25}\right]

and the disk component by:

 I_D(r_D)=I_D(0)\ exp\left[\frac{r_D}{h}\right]

where I(0) is the central surface brightness, re is the bulge effective (half-light) radius and h is the disk scale length.

If (xc,yc) give the position of the galaxy center, then at a position (x,y), we have dx = x - xc and dy = y - yc,

 dxB = dxcosB) + dysinB)
 dyB = [ − dxsinB) + dycosB)] / arB


  r_B^2 = dx_B^2 + dy_B^2 

where θB is the position angle of the major axis of the bulge component and arB is the axial ratio (minor/major) of the bulge. A similar equation holds for the disk component. The position angles of the two component are allowed to vary independently.

Fitting using GIM2D

GIM2D needs:

1) Source extraction and corresponding segmentation image as well as a catalog of individual galaxies that is done using SExtractor. Then, using objects coordinates from this catalog, IRAF task xgal creates thumbnails to be used by GIM2D.

2) Thumbnail (postage stamps) images of all galaxies from the input image. The advantage of this is that all these images may be processed independently using different CPUs (in a cluster, for example).

3) PSF for galaxies that appear small on the images (but I use it anyway), that can be obtained either by cutting uncrowded, unsaturated star from the original image (quick and dirty) or in the usual manner with IRAF's DAOPHOT package (imexamine->pstselect->psf->seepsf). The best manual for creating PSF in IRAF is in Exercise V. Note: There are fits files available for all the exercises.

FIRST STEP: SExtractor

Once installed :), SExtractor (SE) runs from the command line:

 sex -c sex.default image.fits

Run SE from the parent dir (where it is installed), since it uses many files there (default.conv, and default.param). Or copy all these files in the working dir. Now, you must change the following lines in

                                # ASCII_VOTABLE, FITS_1.0 or FITS_LDAC
DETECT_MINAREA   10             # minimum number of pixels above threshold        CHANGE -- EXPERIMENT 
DETECT_THRESH    2              # <sigmas> or <threshold>,<ZP> in mag.arcsec-2    CHANGE -- EXPERIMENT
ANALYSIS_THRESH  2              # <sigmas> or <threshold>,<ZP> in mag.arcsec-2    CHANGE -- EXPERIMENT
DEBLEND_MINCONT  0.005          # Minimum contrast parameter for deblending       CHANGE IN CROWDED FIELDS
                                # MINIBACKGROUND, MINIBACK_RMS, -BACKGROUND,      (you can add others, but it 
                                # FILTERED, OBJECTS, -OBJECTS, SEGMENTATION,      increases the calculation time)
                                # or APERTURES
CHECKIMAGE_NAME  image_seg.fits # Filename for the check-image

In short, CATALOG_TYPE may be also FITS_LDAC, but its a binary file you cannot read, so leave ascii. DON'T use ASCII_HEAD since IRAF gets confused with the header. DETECT_MINAREA is a parameter you REALLY NEED to vary until you get satisfactory segmentation image. This parameter tells SE how bright objects to include along with DETECT_THRESH that should be less than or equal to ANALYSIS_THRESH. The combination given here is good for nearby galaxies. If you have one galaxy that dominates the image (and this is the one you need), you should change these parameters and run SE until the output says:

 ----- SExtractor 2.8.6 started on 2010-03-08 at 13:34:17 with 2 threads
 Measuring from: " "  / 352 x 502 / 0 bits FLOATING POINT data
 (M+D) Background: 72.7632    RMS: 0.540276   / Threshold: 1.08055    
 Objects: detected 1        / sextracted 1                      
 > All done (in 0 s)

which means you have extracted only one object == galaxy you need. Now, if you open the segmentation image named image_seg.fits, for example with ds9, you will see that the pixels of your galaxy differ from the rest of the image. That is because pixels of your galaxy have the value 1, and all the others 0 (check with imexamine if you don't believe). This is the thing that tells GIM2D where is the object and before that to XGAL to make a thumbnail around the center of the galaxy given in file (read the file and check - first two entries are the coordinates). Now, to include the necessary entries in the catalog file, modify default.param:


Now you know what other entries are. And NOW YOU CAN RUN SEXTRACTOR and play around.


Very good step-by-step tutorial is here.

Once installed :), run IRAF and navigate to fuzzy.gim2d. Pay attention to the following tasks:

  • gimfit2d
    in_image = ""              Name of input image
   out_image = ""              Name of output image
  (sym_image = "")             Name of symmetrized image
  (res_image = "")             Name of residual image
  (msk_image = "")             Name of pixel mask image
  (sig_image = "none")         Name of sigma image
  (logfile = "image.log")      Name of log file 
  (mdparfile = "image.mdpar")  Name of parameter value file 
  (initparam = no)      Determine parameters from input image moments (ye
      (dosym = no)             Symmetrize input image (yes/no)?
      (dobkg = no)      Recompute background parameters (yes/no)?
      (imscale = 1.)             Image scale (arcsecond/pixel) 
  (psftype = "user")           PSF type (delta|gaussian|tiny_wfpc|tiny_nic|tiny_ 
     (in_psf = "")              Name of tinytim/user PSF image 
     (c_abs = 0.)               Disk internal absorption coefficient (0-1)
  (seeing = 0.95)              Seeing FWHM (arcseconds)  
  (bkgsig = 72.73)             Background sigma (ADU)
  (nsamp = 200)                Number of parameter space samples 
  (ccdgain = 8.)               CCD Gain (electrons/DU)         
     (nframe = 1)              Number of frames combined
    (metseed = -24781)         Metropolis seed
  (osamp = 1)                  Core oversampling factor          
       (mode = "q")            
  • gscripter
 in_image = "image.fits"       Name of input image    
 in_file = ""          Name of input SExtracted file 
 (out_file1 = "image.xgal")    Name of output XGAL file  
 (out_file2 = "")      Name of output IRAF GIMFIT2D script 
    (dopyraf = no)             GIMFIT2D script to be in PYRAF form (yes/no)?
 (path = "/home/user/..where the images are/") Directory path 
     (prefix = "g")            Prefix for postage image names
 (seed = -22911)               Random generator seed     
       (rmin = 12.)            Minimum extraction radius
      (earea = 20.)            Extraction area
     (varpsf = no)             Are you using a different PSF for each object?
      (sigim = no)             Are you using sigma images?
 (psfname = "image.psf.fits")  PSF name to use when varpsf=no 
       (mode = "q")            
  • xgal
 in_image = "image.fits"      Name of input image  
 msk_image = "image_seg.fits" Name of input mask image  
 in_list = "image.xgal"       Name of input file  
   (sig_image = "none")       Name of input sigma image
       (mode = "q")            

In dark red color the lines should be changed. Mostly replace the existing with the names of your files and images. Note the file image.mdpar. We will create that one below. For now just give it the name with the extension .mdpar. BE CAREFUL with ccdgain, seeing and bkgsig -- find these parameters in the header of the original image or google it if you are using archival data. Parameter imscale is important only if you use gaussian as a PSF (try not to do that). DONT -- UNDER ANY CIRCUMSTANCES -- LET GIM2D TO DETERMINE THE BACKGROUND (use the SExtractor estimate)!!!. It will run forever, and after an hour or so you will obtain very bad results.

Run gscripter and check output files: image.xgal (the same as with the names of the images to be created with XGAL as thumbnails gx???y???.fits=galaxy and gx???y???_m.fits=segmentation image, where ??? are actual pixel coordinates of the center of the objects) and (the script you will run - check if bkgsig is equal to the 4th parameter in your file). Now, run xgal and check the images. Edit the header of gx???y???.fits image with the parameters GIM2D needs: SKYBKG, EXPTIME, GAIN, RDNOISE. For example (from IRAF):

 hedit g???y???.fits SKYBKG 72.73 add+  // the number 72.73 is again that 4th parameter in (SExtractor estimate of the skybkg) 

Now, assuming you have created PSF image, named image.psf.fits as we specified in gscripter task, you also need to create parameter file of initial guesses for all the parameters to be fitted, the magnificent image.mdpar file:

 200                                # nasamp  -- don't decrease, but an increase will increase calculation time
 500.0 500.0 1.0 1000.0 1000.0      # Flux
 0.5 0.5 0.0 1.0 1.0                # B/T  (0=pure disk) 
 2.5 2.5 0.0 5.0 5.0                # bulge Re 
 0.35 0.35 0.0 0.7 0.7              # bulge e 
 90.0 90.0 -360.0 360.0 180.0       # bulge position angle
 2.5 2.5 0.0 5.0 5.0                # exponential disk scale length (h)
 10.0 50.0 0.0 90.0 90.0            # disk inclination (0=face-on)
 90.0 90.0 -360.0 360.0 180.0       # disk position angle
 0.0 0.0 -1.0 1.0 0.1               # subpixel x offset 
 0.0 0.0 -1.0 1.0 0.1               # subpixel x offset 
 72.81 72.81 72.81 72.81 0.0        # background level db in DU
 2.0 2.0 0.2 4.0 3.8                # bulge Sersic index nb (nb = 4 for de Vaucouleurs profile)
  |   |   |   |   |
guess best min max range ... values

Without this file, GIM2D wont run. In the 1st and 2nd column you give initial guess and the best value (that should be the same) to GIM2D. The 3th and 4th column are the minimum and maximum value you want to explore and in the last column is the perurbation=variation of your minimum. The line with background level (2nd from the bottom) should have the value of YOUR bkgsig, SExtractor estimated for you (4th param in or you've found in the header or estimated with imexam. This example is the most complex GIM2D can process - you have Sersic bulge + exponential disk. You run gimfit2d in background with this command:

 cl < &>> gim.log & 

With the above deV.mdpar file, the result looks like:

Image:Galaxy.jpg Image:Galaxy_seg.jpg Image:Galaxy_r.jpg Image:Galaxy_o.jpg

  • *** original image ****** segmentation image ****** residual image ****** output model

Which is not good, since the residual image should contain only the background. The residual image is created by subtraction of the output model convolved with the psf from the original image (residual=original-model*psf). This means the model is not good. The output of this fit named gx???y???.log (large file, look at the lines pasted here):

P 241416. 226972.4 243755.1
P 0.09026347 0. 0.2217385
P 45.77557 37.7331 54.48035
P 0.3077789 0. 0.7
P -72.3067 -212.1298 4.336335
P 37.21162 33.82701 41.71941
P 8.583362 0. 22.07571
P -180.5235 -343.9721 -139.2982
P 0.699858 0.01791382 1.004733
P -3.417662 -3.781986 -3.004471
P 72.81 72.81 72.81
P 3.584051 2.408572 4.
      |        |    |
   best_value lower upper 3σ interval
CHI 0.01538232

In the same order as in .mdpar file you will find fitted values of your parameters here. In the 1st column is the best fitted value and in the 2nd and the 3th are the lower and upper 3σ errors. The lines you can and should change in the image.mdpar file are B/T and the last line to obtain different (meaningful) light profiles:

Sersic -- single component:
B/T line : 1.0 1.0 1.0 1.0 0.0  # B/T=1 fixed (no disk)
last line: 2.1 2.1 0.2 4.0 3.8
de Vaucouleurs bulge + exponential disk:
B/T line: 0.5 0.5 0.0 1.0 1.0
last line: 4.0 4.0 4.0 4.0 0.0  # you fix Sersic to exp=4.0 since the last number is the change and it is 0.0 
Sersic bulge + exponential disk:
B/T line: 0.5 0.5 0.0 1.0 1.0
last line: 2.1 2.1 0.2 4.0 3.8

For example, Millennium Galaxy Catalog is compiled with these different combinations of light profiles in GIM2D. Many "why to" you can find answered in the this short paper related to MGC catalog and in this not so short but VERY GOOD paper.

Personal tools