Commit e2fde158 authored by bepinat's avatar bepinat

Ajout de create_config.py

parent adb8e04d
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
import pyfits as pf
import astropy.wcs as wcs
import astropy.coordinates as coord
import astropy.units as u
def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=True, dz=0.005, dxy=15, initw=50., wmin=30., wmax=250., dfit=100, degcont=0, sclip=10, xyclip=3, nclip=3, wsmooth=0, ssmooth=2):
'''
This function enables to create a configuration file for camel.
Parameters
----------
path: string
indicates the path where to find the cube, the catalogue and where to write the configuration files
cubefile: string
name of the input cube (fits format)
varfile: string
name of the input variance cube or spectrum (fits format)
catfile: string
name of the catalogue file: must contain a header with at least four columns: 'ID', 'ra', 'dec', 'z'
lines: string
indicates which lines to fit: 'ha' for Halpha, 'hb' for Hbeta, 'hg' for Hgamma, 'hd' for Hdelta, 'n2' for NII6548 and NII6583, 's2' for SII6716 and SII6731, 'o2' for OII, 'o3' for OIII4959 and OIII5007. The concatenation of several lines is possible to fit several lines.
suffixeout: string
suffixe for the output configartion file name
commonw: bool
specifies if a common width has to be used for the various lines
dz: float
interval of redshift around the systemic redshift (from the catalogue) to search for lines and for cutting the cube spectrally
dxyz: int
spatial extent around the object for cutting the cube (in pixels)
initw: float
initial guess for the line width (dispersion) in km/s
wmin: float
minimum width allowed in km/s
wmax: float
maximum width allowed in km/s
dfit: float
interval step for fitting the line position in km/s. Large values speed up the computation but could miss the true position.
degcont: int
degree of the polynomial used to fit the continuum
sclip: int
sigma clipping threshold (sclip * sigma) to clean the cube
xyclip: int
size of the local bos for the sigma-clipping (in pixels)
nclip: int
number of pass for the sigma-clipping
wsmooth: int
spectral smooth method
ssmooth: float
gaussian spatial smooth FWHM in pixels
'''
cat = np.genfromtxt(path + catfile, names=True)
hdul = pf.open(path + cubefile)
cube = hdul[1].data
hdr = hdul[1].header
cunit = hdr['CUNIT3']
crpix = hdr['CRPIX3']
crval = hdr['CRVAL3']
cdelt = hdr['CD3_3']
#print(cunit,crpix,crval,cdelt)
if 'CTYPE3' in hdr.ascard.keys(): del hdr['CTYPE3']
if 'CRVAL3' in hdr.ascard.keys(): del hdr['CRVAL3']
if 'CRPIX3' in hdr.ascard.keys(): del hdr['CRPIX3']
if 'CUNIT3' in hdr.ascard.keys(): del hdr['CUNIT3']
if 'CD3_3' in hdr.ascard.keys(): del hdr['CD3_3']
if 'CD3_2' in hdr.ascard.keys(): del hdr['CD3_2']
if 'CD3_1' in hdr.ascard.keys(): del hdr['CD3_1']
if 'CD2_3' in hdr.ascard.keys(): del hdr['CD2_3']
if 'CD1_3' in hdr.ascard.keys(): del hdr['CD1_3']
wlim = np.zeros(np.shape(cube)[1:])
hdu = pf.PrimaryHDU(data=wlim, header=hdr)
hdul1 = pf.HDUList(hdu)
hdr = hdul1[0].header
w = wcs.WCS(hdr[:28], hdul1)
for line in cat:
obj = int(line['ID'])
suffixe = suffixeout + '_' + str(obj) + '_' + lines
fileout = suffixe + '.config'
fo = open(path + fileout, 'w')
fo.write('FITSFILE= ' + cubefile + ' / path of the fits file containing the cube\n')
fo.write('OUTPUT = ' + suffixe + '/' + suffixe + ' / where to write the output (default: ./out.fits)\n')
fo.write('SKYFILE = ' + varfile + ' / path of the file containing the sky spectrum or cube\n')
lbda = []
if 'ha' in lines:
fo.write('HALPHA = True / do we fit Halpha (default: TRUE)\n')
lbda.append(6562.801)
else:
fo.write('HALPHA = False / do we fit Halpha (default: TRUE)\n')
if 'hb' in lines:
fo.write('HBETA = True / do we fit Hbeta4861 (default: FALSE)\n')
lbda.append(4861.363)
if 'o3' in lines:
fo.write('OIII4959= True / do we fit OIII4959 (default: FALSE)\n')
fo.write('OIII5007= True / do we fit OIII5007 (default: FALSE)\n')
lbda.append(4958.911)
lbda.append(5006.843)
if 'o2' in lines:
fo.write('OII = True / do we fit OII doublet (default: FALSE)\n')
lbda.append(3726.04)
lbda.append(3728.80)
if 'n2' in lines:
fo.write('NII6548 = True / do we fit NII6548 (default: FALSE)\n')
fo.write('NII6583 = True / do we fit NII6583 (default: FALSE)\n')
lbda.append(6548.05)
lbda.append(6583.45)
if 's2' in lines:
fo.write('SII6716 = True / do we fit SII6716 (default: FALSE)\n')
fo.write('SII6731 = True / do we fit SII6731 (default: FALSE)\n')
lbda.append(6716.44)
lbda.append(6730.82)
if 'hg' in lines:
fo.write('HGAMMA = True / do we fit Hgamma4340 (default: FALSE)\n')
lbda.append(4340.47)
if 'hd' in lines:
fo.write('HDELTA = True / do we fit Hdelta4102 (default: FALSE)\n')
lbda.append(4101.73)
fo.write('COMMW = ' + str(commonw) + ' / do we use a common width (default: FALSE)\n')
fo.write('REDSHIFT= ' + str(line['z']) + ' / initial redshift of the galaxy\n')
fo.write('REDMIN = ' + str(line['z'] - dz) + ' / minimum redshift\n')
fo.write('REDMAX = ' + str(line['z'] + dz) + ' / maximum redshift\n')
fo.write('INITW = ' + str(initw) + ' / initial linewidth in km/s during the line fitting process \n')
fo.write('WMIN = ' + str(wmin) + ' / minimum linewidth in km/s\n')
fo.write('WMAX = ' + str(wmax) + ' / maximum linewidth in km/s\n')
fo.write('DFIT = ' + str(dfit) + ' / bin for fitting the line center in km/s\n')
fo.write('DGCTNUM = ' + str(degcont) + ' / degre n for the fit of the continuum a0+...+a_n.x^n\n')
fo.write('SCLIP = ' + str(sclip) + ' / Sigma-clipping threshold to clean the cube\n')
fo.write('XYCLIP = ' + str(xyclip) + ' / Width of the box to compute the median for the sigma-clipping\n')
fo.write('NCLIP = ' + str(nclip) + ' / Number of pass for clipping (default=3)\n')
fo.write('WSMOOTH = ' + str(wsmooth) + ' / spectral smoothing: 1 for Hanning, 2 for gaussian (FWHM=3 pixels), set 0 for no smoothing (default:0)\n')
fo.write('SSMOOTH = ' + str(ssmooth) + ' / gaussian spatial smooth FWHM in pixel, set 0 for no smoothing (default:0)\n')
ra = coord.Angle(line['ra'], unit=u.degree)
dec = coord.Angle(line['dec'], unit=u.degree)
xc, yc = np.round(w.wcs_world2pix(ra, dec, 0))
minlbda = np.min(np.array(lbda)) * (1 + line['z'] - 1.5 * dz)
maxlbda = np.max(np.array(lbda)) * (1 + line['z'] + 1.5 * dz)
zmin = np.round((minlbda - crval) / cdelt + crpix - 1)
zmax = np.round((maxlbda - crval) / cdelt + crpix - 1)
fo.write('XMIN = ' + str(xc - dxy) + ' / X minimum value for spatial cut of the input cube (pixel)\n')
fo.write('XMAX = ' + str(xc + dxy) + ' / X maximum value for spatial cut of the input cube (pixel)\n')
fo.write('YMIN = ' + str(yc - dxy) + ' / Y minimum value for spatial cut of the input cube (pixel)\n')
fo.write('YMAX = ' + str(yc + dxy) + ' / Y maximum value for spatial cut of the input cube (pixel)\n')
fo.write('ZMIN = ' + str(zmin) + ' / Z minimum value for spectral cut of the input cube (spectral channel, integer)\n')
fo.write('ZMAX = ' + str(zmax) + ' / Z maximum value for spectral cut of the input cube (spectral channel, integer)')
fo.close()
def main():
gr = '116'
#gr = '83'
#gr = '28'
#gr = '32'
path = '/home/bepinat/Instruments/MUSE/groups/kinematics/gr' + gr + '/'
cubefile = 'DATACUBE_FINAL_ZAP_COSMOS-GR' + gr + '_newWCS.fits'
varfile = cubefile
catfile = 'catalog_gr' + gr + '.txt'
suffixeout = 'gr' + gr
lines = 'o3hb'
lines = 'o2'
create_config(path, cubefile, varfile, catfile, lines, suffixeout)
if __name__ == "__main__":
main()
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment