Material prepared for the 52nd Saas-Fee Advanced Course on "The circum-galactic medium across cosmic time: an observational and modelling challenge"¶

by Alejandro Benitez-Llambay & Michele Fumagalli¶

Contact for questions or comments: alejandro.benitezllambay@unimib.it

OII Surface Brightness images from MUSE datacubes¶

In this exercise, you have access to MUSE datacubes of bright OII emitters. These have already been reduced and centered around the relevant source. Your task is:

  1. Read the files and familiarize yourself with the data.
  2. Extract the spectrum of the source and its associated uncertainty. You can do this for the entire spatial range or focus on select regions.
  3. Find the line profile and transform the spectrum's wavelength into relative velocity with respect to the peak of the OII line.
  4. Stack the cube in velocity to produce an OII surface brightness image with its associated uncertainty.
  5. Find the velocity window that maximizes the signal-to-noise ratio

Example solution¶

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from scipy.optimize import curve_fit
import glob

To carry out (1), we write a helper class that can handle the relevant MUSE datacube data

In [2]:
# We define a helper class to handle the datacube
class muse_file():
    def __init__(self, filename):
        current_cube = fits.open(filename)
        #flux in electron counts
        self.flux   = current_cube[0].data 
        #variance as reported by the intrument
        self.var    = current_cube[1].data 
        #header with useful information
        self.header = current_cube[0].header  
        #dimensions of the cube
        self.nz, self.ny, self.nx = self.flux.shape  

        #zero point for wavelength and space.
        self.l0 = self.header["CRVAL3"]
        self.x0 = self.header["CRVAL1"] 
        self.y0 = self.header["CRVAL2"]

        #Spatial sampling in arcsec
        self.dx = self.header['CD2_2'] * 60 * 60
        self.dy = self.header['CD2_2'] * 60 * 60

        #Wavelength sampling is done every 0.125 A
        self.dl = self.header["CD3_3"]

        self.l = self.l0 + np.arange(0,self.nz,1)*self.dl
        self.x = self.x0 + np.arange(0,self.nx,1)*self.dx
        self.y = self.y0 + np.arange(0,self.ny,1)*self.dy

Next, we become familiar with the data (Task 1) and extract the spectrum considering the uncertainties (Task 2), for the whole sample.

In [3]:
input_files = glob.glob('../MUSECubes/OII/*.fits')

fig = plt.figure(1,figsize=(10,10), tight_layout=True)
axes = [fig.add_subplot(int(len(input_files)/2)+1,2,i+1) for i in range(len(input_files))]


for input_file, ax in zip(input_files, axes):
    cube = muse_file(input_file)
    
    #Get the stack:
    spectrum_sum = np.nansum(cube.flux, axis=(1,2))

    #uncertainty associated with the stack
    spectrum_var = np.nansum(cube.var, axis=(1,2))
    spectrum_sigma = np.sqrt(spectrum_var)

    ax.errorbar(cube.l, spectrum_sum, yerr=spectrum_sigma, capsize=3, color='k', drawstyle='steps-mid')    
    ax.grid()
    
    ax.set_title(input_file)
    ax.set_xlabel('Wavelength / A')
    ax.set_ylabel('Flux')
    
plt.show()

in which the units of the y-axis is (10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ $A^{-1}$)

For task 3 we need to find the line profile and transform the spectrum's wavelength into relative velocity with respect to the peak of the Ly-alpha line. To show how this can be done, we consider a simple worked example: ("oii_cube1.fits")

Worked example: SB image¶

In [4]:
cube = muse_file('../MUSECubes/OII/oii_cube1.fits')

#Get the stack
spectrum_sum = np.nansum(cube.flux, axis=(1,2))

#uncertainty
spectrum_var = np.nansum(cube.var, axis=(1,2))
sigma_spec = np.sqrt(spectrum_var)

kmax = np.argmax(spectrum_sum)
lmax = cube.l[kmax]

c = 3e5 #km /s 
delta_v = c*(cube.l-lmax)/lmax

plt.errorbar(delta_v, spectrum_sum, yerr=spectrum_sigma, capsize=3, color='k', drawstyle='steps-mid')    
plt.xlabel(r'$\rm \Delta v [km / s]$')
plt.ylabel('Flux / 10$^{-20}$ erg s$^{-1}$ cm$^{-2}$ $A^{-1}$')
plt.grid()

Next, we carry out 4 and 5 simultaneously.

We aim to produce a surface brightness profile of the data cube, centred at the OII emission and maximise the S/N ratio.

If we stack the image within the entire velocity range, the signal per pixes becomes:

$f(x,y)=\displaystyle\sum_{i=0}^{N-1} f_{v_{i}}(x,y)$,

where the summation is done over the velocity range.

The uncertainty associated with the stack is thus:

$\sigma_{f} = \sqrt{\displaystyle\sum_{i=0}^{N-1} Var(v_{i})}$.

The signal to noise at the location $(x,y)$ is:

S/N = $\displaystyle\frac{f(x,y)}{\sigma_{f}}$.

If the signal is small beyond some velocity, it is clear that the S/N will deacrese. Thus, we expect the existence of an optimal velocity width that maximises the S/N. Our goal is to find a velocity window that maximises the S/N.

In [5]:
cube = muse_file('../MUSECubes/OII/oii_cube1.fits')

#Get the stack
spectrum_sum = np.nansum(cube.flux, axis=(1,2))

#uncertainty 
spectrum_var = np.nansum(cube.var, axis=(1,2))
spectrum_sigma = np.sqrt(spectrum_var)

kmax = np.argmax(spectrum_sum)
lmax = cube.l[kmax] 

c = 3e5 #km /s 
delta_v = c*(cube.l-lmax)/lmax
dv = delta_v[1]-delta_v[0]

#List to store the mean S/N of the image
SNs_mean = []
SNs_max = []
dvs = []
dvs_ids = []
for i in range(30):
    #Velocity window size
    dv_id = i

    #Stack the image in velocity within the window
    f = np.nansum(cube.flux[kmax-dv_id:kmax+dv_id+1,:,:], axis=0) * dv
    
    #Calculate the uncertainty in the pixel
    sigma_f = np.sqrt(np.nansum(cube.var[kmax-dv_id:kmax+dv_id+1,:,:], axis=0)) * dv
    
    #Store a measure of the image's S/N (e.g., the average S/N, or the maxmimum).
    SNs_mean.append(np.mean(f/sigma_f))
    SNs_max.append(np.max(f/sigma_f))
    dvs.append(2.0 * dv_id * dv)
    dvs_ids.append(dv_id)
        
idmax_SNs_mean = np.argmax(SNs_mean)
idmax_SNs_max = np.argmax(SNs_max)

plt.plot(dvs, SNs_mean/SNs_mean[idmax_SNs_mean], 'o-', label='mean S/N')
plt.plot(dvs, SNs_max/SNs_max[idmax_SNs_max], 'o-', label='max S/N')

plt.axvline(dvs[idmax_SNs_mean], color='C00', ls='--')
plt.axvline(dvs[idmax_SNs_max], color='C01', ls='--', label='Optimal velocity width')

plt.legend()

plt.xlabel('Velocity width [km/s]')
plt.ylabel('S/N / Max(S/N)')
plt.grid()

To understand what is going visually, consider the following

In [6]:
cube = muse_file('../MUSECubes/OII/oii_cube1.fits')

#Get the stack
spectrum_sum = np.nansum(cube.flux, axis=(1,2))

#uncertainty
spectrum_var = np.nansum(cube.var, axis=(1,2))
spectrum_sigma = np.sqrt(spectrum_var)

kmax = np.argmax(spectrum_sum)
lmax = cube.l[kmax] 

c = 3e5 #km /s 
delta_v = c*(cube.l-lmax)/lmax
dv = delta_v[1]-delta_v[0]

#Containers for the mean S/N of the image
SNs_mean = []
SNs_max = []
dvs = []
dvs_ids = []

fig = plt.figure(1, figsize=(10,10), tight_layout=True)
axes = [fig.add_subplot(4,3,i+1) for i in range(12)]

#hand-picked windows
ii = [0,1,2,3,4,5,6,10,14,20,22,25,28]
for i, ax in zip(ii, axes):
    #Window size
    dv_id = i

    #Stack the image in velocity within the window
    f = np.nansum(cube.flux[kmax-dv_id:kmax+dv_id+1,:,:], axis=0) * dv
    
    #Calculate the uncertainty in the pixel
    sigma_f = np.sqrt(np.nansum(cube.var[kmax-dv_id:kmax+dv_id+1,:,:], axis=0) ) * dv
    
    #Store a measure of the image's S/N (e.g., the average S/N, or the maxmimum).
    SNs_mean.append(np.mean(f/sigma_f))
    SNs_max.append(np.max(f/sigma_f))
    dvs.append(2.0 * dv_id * dv)
    dvs_ids.append(dv_id)
    
    ax.imshow(f/sigma_f, origin='lower', extent=[cube.x[0],cube.x[-1],cube.y[0],cube.y[-1]], vmin=0,vmax=50)
    ax.set_xlabel('Angular sep (arcsec)')
    ax.set_ylabel('Angular sep (arcsec)')
    ax.set_title('dv=%.2f, Max(S/N)=%.2f'%((2.0*dv_id*dv), np.max(f/sigma_f)))
plt.show()