---
# System prepended metadata

title: HowToLens with GLEE

---

# Lens Modeling: How to



## 0. Handling FITS Files and Generating Needed Data


### FITS Format
`.fits` is the most used digital file format in astronomy. These files are opened with several different software, but we'll use **[SAOImageDS9](https://sites.google.com/cfa.harvard.edu/saoimageds9)**. The main functionalities we'll need are the viewing of data and the creation of regions. See: [DS9Guide](https://www.jb.man.ac.uk/~gbendo/Sci/Pict/DS9guide.pdf)

A `.fits` file consists of 2 main parts, the `header` and the `data`. 

A typical `header` looks like this

```
IMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -32 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                 2362                                                  
NAXIS2  =                 2423                                                  
EXTEND  =                    T                                                  
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
FILENAME= '../reduced_data/WG0214-2105_F160W_drz_sci.fits' / name of file       
FILETYPE= 'SCI      '          / type of data found in data file                
                                                                                
TELESCOP= 'HST'                / telescope used to acquire data                 
INSTRUME= 'WFC3  '             / identifier for instrument used to acquire data 
EQUINOX =               2000.0 / equinox of celestial coord. system             
                                                                                
              / DATA DESCRIPTION KEYWORDS                                       
                                                                                
ROOTNAME= '../reduced_data/WG0214-2105_F160W' / rootname of the observation set 
IMAGETYP= 'EXT               ' / type of exposure identifier                    
PRIMESI = 'WFC3  '             / instrument designated as prime                 
                                                                                
              / TARGET INFORMATION                                              
                                                                                
TARGNAME= 'WG0214-2105                    ' / proposer's target name            
RA_TARG =   3.356812500000E+01 / right ascension of the target (deg) (J2000)    
DEC_TARG=  -2.109304722222E+01 / declination of the target (deg) (J2000)        
                                                                                
              / PROPOSAL INFORMATION                                            
                                                                                
PROPOSID=                15652 / PEP proposal identifier                        
LINENUM = 'Z3.009         '    / proposal logsheet line number                  
PR_INV_L= 'Treu                          ' / last name of principal investigator
PR_INV_F= 'Tommaso             ' / first name of principal investigator         
PR_INV_M= 'L.                  ' / middle name / initial of principal investigator
               
                                 [...]
                                 [...]
                                 [...]                                 
```


It includes the metadata corresponding to the observation done. Among others you can find useful information about exposure time, pixel size, WCS, etc.

On the other side, the `data` is multidimensional array, where each entry corresponds to a pixel In essence, the array 
```
data = np.array([[1, 0, 0], [0, 1, 0], [0,0,1]])
```       
represents the following image:
<center>
    
![](https://hackmd.io/_uploads/B1-SGsJ6h.png =400x400)
</center>

### Cutout

Observation data includes a lot of information which is not needed, so in order to model we would like to have a cutout of the system we're trying to model.

In order to do this we use the `Cutout2D` function from `astropy`

```
class Cutout2D(
    data: data,
    position: tuple[int, int],
    size: int,
    wcs: Any | None = None,
    mode: str = "trim",
    fill_value: float = np.nan,
    copy: bool = False
)
```

where `data` is the data from the initial `.fits` file. Note: It is good practice to update the header (particularly the WCS) accordingly. An example of the script to do a cutout can be seen here:

```
import astropy
from astropy.io import fits
import matplotlib.pyplot as plt
from astropy.visualization import astropy_mpl_style
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
import tkinter, tkinter.filedialog
import matplotlib.pyplot as plt
tkinter.Tk().withdraw()


# Open the SCI file
print("Select the SCI File to be cropped")
file_path = tkinter.filedialog.askopenfilename()
hdulist = fits.open(file_path)

#Initialize variables for header, data and WCS
header = hdulist['PRIMARY'].header
data = hdulist['PRIMARY'].data
wcs = WCS(header)

# Read the coordinates for the center of the cutout
c0x =int(input('Cutout center, x-axis, pixels: \n'))
c0y =int(input('Cutout center, y-axis, pixels: \n'))
center = (c0x, c0y)

# Read size of cutout
size=int(input('Size of the (square) cutout, pixels: \n'))

#Use Cutout2D to generate a cutout
cutout = Cutout2D(data, center, size, wcs=wcs)


#Keep the header but modified for the cutout - this is good practice
cheader = cutout.wcs.to_header()
primaryhdu = fits.PrimaryHDU(cutout.data, cheader)
hdulist = fits.HDUList([primaryhdu])

#plot the file
plt.figure()
plt.style.use(astropy_mpl_style)
plt.imshow(cutout.data, cmap='gray',origin= 'lower')
plt.colorbar()
plt.show()

# Save the file
hdulist.writeto('Cutout.fits', overwrite=True)

```

![](https://hackmd.io/_uploads/H1kSQnyan.png)



### Background Reduction

Similarly, we'll use `astropy` to reduce the background of the image. We will subtract a floor that leaves a very small (close to zero, they can be negative!) value in every pixel. To obtain the floor value we make a region in ds9 on an empty patch. In the region statistics we look at the computed median and use this as our floor.

```
import astropy
import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
from astropy.visualization import astropy_mpl_style
from astropy.wcs import WCS
from astropy.nddata import Cutout2D
import tkinter, tkinter.filedialog
import matplotlib.pyplot as plt
tkinter.Tk().withdraw()


# Open the SCI file
print("Select the SCI File to reduce background from")
file_path = tkinter.filedialog.askopenfilename()
hdulist = fits.open(file_path)
 
#Initialize variables for header, data and WCS
header = hdulist['PRIMARY'].header
data = hdulist['PRIMARY'].data
wcs = WCS(header)

# Reduce the value in each pixel
bg= float(input("Input the background value: \n"))
newdata = data - bg


#Keep the header
primaryhdu = fits.PrimaryHDU(newdata, header)

#Create fits file
hdulist = fits.HDUList([primaryhdu])

# Save the file
hdulist.writeto('Noisereduced.fits', overwrite=True)

```


### PSF

The Point Spread Function (PSF) is the geometric reaction an optical system has to a point source. In Hubble it is characterized by the diffraction spikes on star. In JWST, newer, it is characterized by the diamont shape that stars have.

![](https://hackmd.io/_uploads/rkuua0D-a.png)

 It affects the observation of the system, as spikes and *airy* rings are visible around the quasars. This airy rings are not to be  confused with the arc generated by the gravitational lenses. The latter will be discused further down this guide.


To generate the PSF please refer to ``starred``, a python package to generate psf from stars in the sci file by Michalewicz et al. 

Files needed for the reconstruction of the PSF are:
- Cutout of the stars in the SCI File. (Make them all the same size).
- Patching the cutout of the stars, using a random normal distribution with a mean and RMS obtained from the BGR Sci file. Make sure the BGR is used.
- Error Map of each star. A basic $\sigma$ map of the star cutout. Use the weight files for this. 
 
In the case of WG0214 in IR band (160), four stars with a decent amount of spikes/rings were used.

![](https://hackmd.io/_uploads/BJB1A0wWT.png)


After feeding these images, error maps, and extra model information, the final PSF reconstruction looks like this:

![](https://hackmd.io/_uploads/B1VmyJd-T.png)

In which the Airy rings around the point source are quite visible. The spikes are also visible if the contrast is adjusted.

**Important:** A PSF generated for one set of observations is NOT compatible with another set of obs. The PSF is wavelenght and construction dependent. Since Hubble is in LEO, the effects of the earth atmosphere, solar induced contraction and expansion, this varies. It also varies in other optical systems, such as the ELT, VLA, JWST, etc. It is always better to obtain the psf from the observed data rather than other sources.

[GitHub Repo](https://gitlab.com/cosmograil/starred)

### Masks

We initially need to do 3 masks. We first make the region files for this.

#### 1-Pixel mask

This is just a mask with all but one single corner pixel set to 0.

#### Arc Mask 

This one covers both, the arc and the light from the source galaxy - leaving the lens galaxy uncovered. Make sure the region file is saved in the image coordinate system (pixels), DS9 format, and double check there is no other figures in the file.

<center>

 ![](https://hackmd.io/_uploads/HJgKkhkT3.png =300x300)  ![](https://hackmd.io/_uploads/BkCF12kph.png =300x300)
</center>

#### Object Mask 

This one covers all the light objects not belonging to the system to be modelled. Save each region individually. Make sure each region file is saved with the image coordinate system (pixels), DS9 format and double check there is no other figures in the file.

<center>

 ![](https://hackmd.io/_uploads/H1Y3g3kTh.png =300x300)  ![](https://hackmd.io/_uploads/BkRRlhyah.png =308x308)
</center>

Now, we use the script to generate the fits files for the masks.

![](https://hackmd.io/_uploads/B1GRVny62.png)

### Error Map

For this we also need a region, which is more or less "empty", so that we can get a sample of the noise in the whole image.

<center>

 ![](https://hackmd.io/_uploads/HJevShJ6h.png =300x300)  ![](https://hackmd.io/_uploads/rk-wBhJa2.png =300x300)
</center>

Then we use this region (Once again, make sure the region file is saved in the image coordinate system, DS9 format, and double check there is no other figures in the file) and the errormap script to generate an errormap.

<center>
    
![](https://hackmd.io/_uploads/Hk6FkNMah.png =300x300)

</center>


## 1. Mass Modelling

### The configfile

We begin with a `configfile` this text file (extensionless) contains the data of the model. Our initial version looks like this:

```
chi2type            1
minimiser       siman
seed                1

siman_iter          1
siman_nT         1000
siman_dS           10
siman_Sf          1.1
siman_k             2
siman_Ti           10
siman_Tf          1.1
siman_Tmin          1

lenses_vary 2

    spemd
    4.80  #x-coord   flat:4.7,4.9  step:0.01
    4.80  #y-coord   flat:4.7,4.9  step:0.01
    0.90  #b/a       flat:0.8,1  step:0.01
    0.0  #theta     flat:0,6.283:  step:0.01
    0.828  #theta_e   flat:0.5,1.2  scale:  step:0.01
    0.0000001E-4  #r_core    exact:
    0.500522  #gam       flat:0.3,0.7  step:0.01

    shear
    0.0000 flat:0,0.3 step:0.01 #magnitude
    0.0000 noprior: step:0.01  #theta

sources 1

    Dds/Ds    1.000000  exact:
    source weighted
    srcx     4.82  exact:
    srcy     4.90  exact:
    4
   4.3200000     3.9200000  error:0.04
   5.6000000     4.7200000  error:0.04
   5.1200000     5.6000000  error:0.04
   4.24000000    5.3600000  error:0.04
```


From top to bottom:

This first part

```
chi2type            1
minimiser       siman
seed                1

siman_iter          1
siman_nT         1000
siman_dS           10
siman_Sf          1.1
siman_k             2
siman_Ti           10
siman_Tf          1.1
siman_Tmin          1
```

contains information on how the modelling is gonna take place. Since we're initially modelling first the `chi2type` we'll use is 1, which corresponds to the source plane positions. We'll be optimizing that first. We'll be initially the using simulated annealing (siman) agorithm with the parameters, Temperature, Step Size, etc., as defined above. 

Now, this part 
```
lenses_vary 2

    spemd
    4.80  #x-coord   flat:4.7,4.9  step:0.01
    4.80  #y-coord   flat:4.7,4.9  step:0.01
    0.90  #b/a       flat:0.8,1  step:0.01
    0.0  #theta     flat:0,6.283:  step:0.01
    0.828  #theta_e   flat:0.5,1.2  scale:  step:0.01
    0.0000001E-4  #r_core    exact:
    0.500522  #gam       flat:0.3,0.7  step:0.01

    shear
    0.0000 flat:0,0.3 step:0.01 #magnitude
    0.0000 noprior: step:0.01  #theta
```

defines the initial values for the lens galaxy. The mass model assumed in the above case is a single power elliptical mass distribution (SPEMD). `x/y-coord` are the coordinates (in arcseconds) of the center of the lens, `b/a` is the eccentricity of the galaxy, `theta` is the angle of the galaxy measured in radians from the +x-axis, `theta_e` is the einstein radius in `arcsecongd`, `r_core` is the radius of the galaxy core. Usually a very small number., lastly `gamm` is a factor in the mass distribution (see literature).

Shear is the distortion caused due to the large scale structures of the universe (0.1% to 1%) and causes some minor effects, which are factored into the model by defining it.

We define the priors to be the interval inbetween which the actual value will be in, not everything has to have a prior, but it is easier with one. `flat:` means the probability of every value is the same, but other options like `gaussian:` are possible (see glee wiki). 


Lastly, the final part of the initial config file

```
sources 1

    Dds/Ds    1.000000  exact:
    source weighted
    srcx     4.82  exact:
    srcy     4.90  exact:
    4
   4.3200000     3.9200000  error:0.04
   5.6000000     4.7200000  error:0.04
   5.1200000     5.6000000  error:0.04
   4.24000000    5.3600000  error:0.04
```

This correspond to the sources. `srcx` and `srcy` are the weighted averages if the position of the images, followed by the number of images (this case 4) and their position.


### Source Plane Position Optimization (chi2type=1) 

We run `glee -m -i configfile` and this returns a `configfile.001` file. (For organization reasons, I'll refer to it as `configfile01`). This new updated file contains the first optimized version of the file. On this new file, you can run `glee -c configfile01` to see the current `chi2` value. This should be samll (less than TBC)


### Image Plane Position Optimization (chi2type=2)

We change the `chi2type` to 2 in `configfile01` and rerun an annealing (`glee -m -i configfile1`). This returns `configfile02`. Again, check that the chi2 is below the value.

### Visualization

We can see the results of this optimization running the `gleeview.py` tool. Before running this, we need to generate the `.fits` cube with the required data (caustics, critical curves, contours, ...), we do this with the `glee -s configfile02`. The result should look like this. The sources (green diamonds) should be bundled up together. The caustics are the red lines, the critical curve is the blue line, and the images are the circle. The center of the lens is the blue square.

![](https://hackmd.io/_uploads/B1AvLpJ6n.png)








## 2. Lens Plane Modelling  

### Lens Light 
To begin modeling the light of the lens, we have to modify the mass-optimized `configfile`. We first set the `chi2type` to 16. Then we set all the mass parameters to `exact:` (leave the `scale:` in `theta_e`). Now we add the sources and we do so by adding the `esource` (extended sources) section. We also add a new algorithm, the MCMC Chains. The file now looks like:

```
chi2type           16
minimiser       siman
seed                1

siman_iter          1
siman_nT         1000
siman_dS           10
siman_Sf          1.1
siman_k             2
siman_Ti           10
siman_Tf          1.1
siman_Tmin          1

mcmc_n         100000
mcmc_dS		0.04375
mcmc_dSini          0
mcmc_k              2

lenses_vary 2

    spemd
    4.794076  #x-coord   exact:
    4.785908  #y-coord   exact:
    0.809780  #b/a       exact:
    0.105977  #theta     exact:
    0.881267  #theta_e   exact:  scale:
  1.0000e-11  #r_core    exact:
    0.540466  #gam       exact:

    shear
    0.083024  #magnitude  exact:
    0.744317  #theta     exact:

sources 1

    Dds/Ds    1.000000  exact:
    source weighted
    srcx      4.335643  exact:
    srcy      4.377210  exact:
    4
    4.320000      3.920000  error:0.04
    5.600000      4.720000  error:0.04
    5.120000      5.600000  error:0.04
    4.240000      5.360000  error:0.04

esources 1

    Dds/Ds    1.000000  exact:
    esource_ngy          20
    esource_ngx          120
    esource_dx           0.080000
    esource_data         /afs/mpa/home/allansch/Lens/WG0214/WG0214_F160W_sci_CO_NR.fits
    esource_err          /afs/mpa/home/allansch/Lens/WG0214/ErrorMap.fits
    esource_arcmask      /afs/mpa/home/allansch/Lens/WG0214/Mask_Arc1.fits
    esource_lensmask     /afs/mpa/home/allansch/Lens/WG0214/Mask_Lens1.fits
    esource_mod_light    LensOnly
    esource_psf          /afs/mpa/home/allansch/Lens/WG0214/PSF.fits
    esource_sub_agn_psf  /afs/mpa/home/allansch/Lens/WG0214/PSF_agn.fits
    esource_sub_agn_psf_factor     3
    esource_sub_esr_psf  /afs/mpa/home/allansch/Lens/WG0214/PSF_esr.fits
    esource_sub_esr_psf_factor     3
    esource_regopt       SpecRegPrecSigFigOnce
    esource_reglampre    1
    esource_reglamnup    1000
    esource_regtype      curv
    esource_reglam       1
    esource_reglamlo     1e-05
    esource_reglamhi     100000
    esource_light  2
        sersic
            4.800000  #x-coord   flat:4.72,4.88  label:sersic_x  step:0.01
            4.800000  #y-coord   flat:4.72,4.88  label:sersic_y  step:0.01
            0.800000  #q         flat:0.6,1  step:0.01
            3.000000  #PA        noprior:  step:0.01
            0.300000  #amp       noprior: min:0  step:0.01
            0.900000  #r_eff     flat:0,1.5  step:0.01
            0.700000  #n_sersic  flat:0.5,9 step:0.01
        sersic
            4.810615  #x-coord   link:sersic_x  a:0,1,1
            4.791608  #y-coord   link:sersic_y  a:0,1,1
            0.800000  #q         flat:0.6,1  step:0.01
            3.000000  #PA        noprior:  step:0.01
            0.300000  #amp       noprior: min:0  step:0.01
            0.900000  #r_eff     flat:0,1.5  step:0.01
            0.700000  #n_sersic  flat:0.5,9 step:0.01 
    esource_end
```

Specific details on what everything means in the GLEE wiki. Here we’re modelling the lens galaxy using 2 sersic profiles. Notice how the centroids of these profiles are linked to each other. The triplet following `a: a1,a2,a3` means the value will be transformed according to: $y = a_1+a_2\cdot x^{a_3}$.

Once everything is set up, you run one siman to jump into parameter space. This is followed by several MCMC chains (done with `glee -M -i configfile0x`). Pay special care that the acceptance has a value of 0.25±0.05. You can adjust this by changing the value of `mcmc_dS`. After each chain (generally 100k-300k iterations), save the covariance matrix using the command `glee -h configfile0x.mcmc` and add it to the next `configfile` by adding the lines

```
sampling_f     gaussian
sampling_cov   /afs/mpa/home/allansch/Lens/WG0214/WG0214_0x.cov
```

after the `mcmc_k` line. Then repeat.

To visualize the returned file, `configfile0x.mcmc` we use `gleeviewsample.py`. here we can see the sampling and it's `log_evidence` value, which tells us how good the fit is and also gives us a corner plot of how the parameters vary with each other. 
<center>

![](https://hackmd.io/_uploads/BJuPMmxp3.png)
</center>

after several chains:

![](https://hackmd.io/_uploads/SkrjSVGph.png)




Furthermore, by running `glee -f 2 configfile0x` we can export `fits` cube that contains (1) the original image, (2) the model, (3) residuals, (4) normalized residuals.



### Lens Light + Quasar Light

Now we set all the parameters regarding the light of the lens to `exact:`, we add the number of images (in this case 4) to the number of profiles in `esource_light`. Also, we now remove the `arcmask` and replace it with a mask of the exact same dimensions, but in which the data is just zeros except for one pixel in any corner, set to 1. This because glee requires a mask, so we put in this minimally invasive mask to make up for it.

The file now looks like
```
chi2type           16
minimiser       siman
seed                1

siman_iter          1
siman_nT         1000
siman_dS           10
siman_Sf          1.1
siman_k             2
siman_Ti           10
siman_Tf          1.1
siman_Tmin          1

mcmc_n         100000
mcmc_dS           0.3
mcmc_dSini          0
mcmc_k              2

lenses_vary 2

    spemd
    4.794076  #x-coord   exact:
    4.785908  #y-coord   exact:
    0.809780  #b/a       exact:
    0.105977  #theta     exact:
    0.881267  #theta_e   exact:  scale:
  1.0000e-11  #r_core    exact:
    0.540466  #gam       exact:

    shear
    0.083024  #magnitude  exact:
    0.744317  #theta     exact:

sources 1

    Dds/Ds    1.000000  exact:
    source weighted
    srcx      4.335643  exact:
    srcy      4.377210  exact:
    4
    4.320000      3.920000  error:0.04
    5.600000      4.720000  error:0.04
    5.120000      5.600000  error:0.04
    4.240000      5.360000  error:0.04

esources 1

    Dds/Ds    1.000000  exact:
    esource_ngy          20
    esource_ngx          120
    esource_dx           0.080000
    esource_data         /afs/mpa/home/allansch/Lens/WG0214/WG0214_F160W_sci_CO_NR2.fits
    esource_err          /afs/mpa/home/allansch/Lens/WG0214/ErrorMap.fits
    esource_arcmask      /afs/mpa/home/allansch/Lens/WG0214/Mask_1pix.fits
    esource_lensmask     /afs/mpa/home/allansch/Lens/WG0214/Mask_Lens1.fits
    esource_mod_light    LensOnly
    esource_psf          /afs/mpa/home/allansch/Lens/WG0214/PSF.fits
    esource_sub_agn_psf  /afs/mpa/home/allansch/Lens/WG0214/PSF_agn.fits
    esource_sub_agn_psf_factor     3
    esource_sub_esr_psf  /afs/mpa/home/allansch/Lens/WG0214/PSF_esr.fits
    esource_sub_esr_psf_factor     3
    esource_regopt       SpecRegPrecSigFigOnce
    esource_reglampre    1
    esource_reglamnup    1000
    esource_regtype      curv
    esource_reglam       9.99288
    esource_reglamlo     1e-05
    esource_reglamhi     100000
    esource_light  6
                    
                                ...
                    
        psf
            5.205020  #x-coord   noprior:  step:0.01
            5.495106  #y-coord   noprior:  step:0.01
          170.869832  #amp       noprior:  min:0  step:2.61399
        psf
            5.580431  #x-coord   noprior:  step:0.01
            4.644163  #y-coord   noprior:  step:0.01
          225.230848  #amp       noprior:  min:0  step:2.88811
        psf
            4.386284  #x-coord   noprior:  step:0.01
            3.922756  #y-coord   noprior:  step:0.01
          351.470253  #amp       noprior:  min:0  step:2.63407
        psf
            4.193880  #x-coord   noprior:  step:0.01
            5.309741  #y-coord   noprior:  step:0.01
          149.845956  #amp       noprior:  min:0  step:1.54551
    esource_end
```