Commit a74e92c9 authored by bepinat's avatar bepinat
Browse files

merge of versions

parent 4a492681
......@@ -40,6 +40,26 @@ logging.basicConfig(level=logging.DEBUG)
logger = logging.getLogger('CAMEL')
def dz_from_dv(dv, z):
'''
This function enables to compute the redshift interval corresponding to a given velocity interval
Parameters
----------
dv: float
velocity interval in km/s
z: float
mean redshift of the source
Returns
-------
dz: float
corresponding redshift interval
'''
dz = dv * 1e3 / ct.c * (1 + z)
return dz
def compute_rv(z, zcosmo=None):
"""This function computes the velocity from local and global redshifts
......@@ -145,7 +165,7 @@ def readconfig(filename, config):
# XXX or give an error message but in that case, must test that the keywords have correct values
# XXX more keywords? add default values?
needed_keys = ['FITSFILE', 'OUTPUT', 'SKYFILE', 'HALPHA', 'NII6548', 'NII6583', 'SII6716', 'SII6731', 'OIII4959', 'OIII5007', 'HBETA', 'OII3729', 'OII3726', 'EXTRAL', 'COMMW', 'REDSHIFT', 'REDMIN', 'REDMAX', 'INITW', 'WMIN', 'WMAX', 'DFIT', 'DGCTNUM', 'MFIT', 'SCLIP', 'XYCLIP', 'NCLIP', 'SPSF', 'WSMOOTH', 'SSMOOTH', 'THRES', 'MEDIAN', 'FITSPSF', 'XMIN', 'YMIN', 'ZMIN', 'XMAX', 'YMAX', 'ZMAX']
needed_keys = ['FITSFILE', 'OUTPUT', 'SKYFILE', 'HALPHA', 'NII6548', 'NII6583', 'SII6716', 'SII6731', 'OIII4959', 'OIII5007', 'HBETA', 'OII3729', 'OII3726', 'EXTRAL', 'COMMW', 'REDSHIFT', 'REDMIN', 'REDMAX', 'ZRMIN', 'ZRMAX', 'INITW', 'WMIN', 'WMAX', 'DFIT', 'DGCTNUM', 'MFIT', 'SCLIP', 'XYCLIP', 'NCLIP', 'SPSF', 'WSMOOTH', 'SSMOOTH', 'THRES', 'MEDIAN', 'FITSPSF', 'XMIN', 'YMIN', 'ZMIN', 'XMAX', 'YMAX', 'ZMAX']
for key in needed_keys:
if not(key in conf.keys()):
conf[key] = None
......@@ -157,6 +177,12 @@ def readconfig(filename, config):
if conf['COMMW'] is None:
conf['COMMW'] = False
# XXX or give an error message but in that case, must test that the keywords have correct values
if conf['ZRMIN'] is None:
deltared = conf['REDMAX'] - conf['REDMIN']
conf['ZRMIN'] = conf['REDMIN'] - deltared / 2.
if conf['ZRMAX'] is None:
deltared = conf['REDMAX'] - conf['REDMIN']
conf['ZRMAX'] = conf['REDMAX'] + deltared / 2.
return conf
......@@ -600,7 +626,6 @@ def build_parinfo_lines(els, linename, conf, free_ratio, comind, multi_comp=1):
zcosmo = conf['REDSHIFT']
zmax = conf['REDMAX']
zmin = conf['REDMIN']
#xxx ICI redefinir zmin et zmax avec le nouveau mot clé pour définir les limites et idem utiliser ce nouveau zmin, zmax pour explorer en z les parametres initiaux
parbase = {'value': 0., 'fixed': 1, 'limited': [0, 0], 'limits': [0., 0.], 'tied': '', 'mpmaxstep': 0, 'mpminstep': 0, 'parname': ''}
parinfo = []
......@@ -743,11 +768,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
#Condition on the wavelength to be within the redshift cuts
#XXX try zmin > zmax?
#xxx CES zmin, zmax sont bien ceux à utiliser pour les bornes en lbda
zcosmo = conf['REDSHIFT']
zmax = conf['REDMAX']
zmin = conf['REDMIN']
zrange = zmax - zmin
zmax = conf['ZRMAX']
zmin = conf['ZRMIN']
argsorted = np.argsort(np.array([els.lines[i].index for i in els.lines]))
lineindex = np.array([els.lines[i].index for i in els.lines])[argsorted]
......@@ -762,8 +784,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
for i in range(linefit.size):
if (not linefit[i]):
continue
linf = np.max([(zmin - zrange / 2. + 1) * linewave[i], np.min(wave0)])
lsup = np.min([(zmax + zrange / 2. + 1) * linewave[i], np.max(wave0)])
linf = np.max([(zmin + 1) * linewave[i], np.min(wave0)])
lsup = np.min([(zmax + 1) * linewave[i], np.max(wave0)])
condo = np.where((wave0 >= linf) & (wave0 <= lsup))
if np.size(condo) != 0:
firstind0[i] = condo[0][0]
......@@ -775,8 +797,8 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
for i in range(linefit.size):
if (not linefit[i]):
continue
linf = np.max([(zmin - zrange / 2. + 1) * linewave[i], np.min(wave0)])
lsup = np.min([(zmax + zrange / 2. + 1) * linewave[i], np.max(wave0)])
linf = np.max([(zmin + 1) * linewave[i], np.min(wave0)])
lsup = np.min([(zmax + 1) * linewave[i], np.max(wave0)])
condo = np.where((wave >= linf) & (wave <= lsup))
if np.size(condo) != 0:
firstind[i] = condo[0][0]
......@@ -819,14 +841,19 @@ def buildmaps(conf, cube, hdr, var=None, plot=False, debug=False, free_ratio=Fal
if var.ndim == 3:
weight = ((1. / np.sqrt(var)) / np.sum(1. / np.sqrt(var), axis=0)) * var.shape[0]
#XXX try zmin > zmax?
zcosmo = conf['REDSHIFT']
zmax = conf['REDMAX']
zmin = conf['REDMIN']
#Redshift step
zstep = (conf['DFIT'] * 1e3 / ct.c + zcosmo) / (1 - conf['DFIT'] * 1e3 / ct.c) - zcosmo
zstep = dz_from_dv(conf['DFIT'], zcosmo)
#zstep = (conf['DFIT'] * 1e3 / ct.c + zcosmo) / (1 - conf['DFIT'] * 1e3 / ct.c) - zcosmo
#if zstep >= (zmax - zmin) / 2: zstep = (zmax - zmin) / 2.
zmean = (zmax + zmin) / 2.
zzmax = zmean + np.floor((zmax - zmean) / zstep + 1) * zstep
zzmin = zmean + np.ceil((zmin - zmean) / zstep) * zstep
#xxx ici utiliser le nouveau zmin, zmax (nouveau mot clé), car cela définit les bornes de l'exploration en z
#Parameters contraints and informations: each line has z, sigma, intensity and ratio constraint
......@@ -1386,7 +1413,9 @@ def main(argv):
parser.add_argument('--REDSHIFT', action="store", dest="REDSHIFT", default=None, type=float, help="mean redshift of the object")
parser.add_argument('--REDMIN', action="store", dest="REDMIN", default=None, type=float, help="minimum redshift allowed for fitting")
parser.add_argument('--REDMAX', action="store", dest="REDMAX", default=None, type=float, help="maximum redshift allowed for fitting")
parser.add_argument('--INITW', action="store", dest="INITW", default=None, type=float, help="initial linewidth in km/s")
parser.add_argument('--ZRMIN', action="store", dest="ZRMIN", default=None, type=float, help="minimum redshift for spectral range considered around lines")
parser.add_argument('--ZRMAX', action="store", dest="ZRMAX", default=None, type=float, help="maximum redshift for spectral range considered around lines")
parser.add_argument('--INITW', action="store", dest="INITW", default=None, type=float, help="initial line width in km/s")
parser.add_argument('--WMIN', action="store", dest="WMIN", default=None, type=float, help="minimum width allowed for fitting in km/s")
parser.add_argument('--WMAX', action="store", dest="WMAX", default=None, type=float, help="maximum width allowed for fitting in km/s")
parser.add_argument('--DFIT', action="store", dest="DFIT", default=None, type=float, help="bin in km/s for searching the position of the lines")
......
This diff is collapsed.
......@@ -17,7 +17,28 @@ import logging
logging.basicConfig(level=logging.DEBUG)
logger = logging.getLogger('CAMEL')
def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=True, dv=1000., dxy=15, deltav=500., initw=50., wmin=30., wmax=250., dfit=100., degcont=0, sclip=10, xyclip=3, nclip=3, wsmooth=0, ssmooth=2.):
def dz_from_dv(dv, z):
'''
This function enables to compute the redshift interval corresponding to a given velocity interval
Parameters
----------
dv: float
velocity interval in km/s
z: float
mean redshift of the source
Returns
-------
dz: float
corresponding redshift interval
'''
dz = dv / ct.c.to('km/s').value * (1 + z)
return dz
def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=True, dv=500., dxy=15, deltav=2000., 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.
......@@ -38,11 +59,11 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
commonw: bool
specifies if a common width has to be used for the various lines
dv: float
interval of velocity around the systemic redshift (from the catalogue) to search for lines and for cutting the cube spectrally
velocity range in km/s to determine redshift range for fit
dxy: int
spatial extent around the object for cutting the cube (in pixels)
deltav: float
velocity range in km/s to determine redshift range for fit
interval of velocity around the systemic redshift (from the catalogue) to search for lines and for cutting the cube spectrally
initw: float
initial guess for the line width (dispersion) in km/s
wmin: float
......@@ -174,24 +195,33 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
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')
if dv is not None:
dz = dv / ct.c.to('km/s').value * (1 + line['Z'])
fo.write('REDMIN = ' + str(line['Z'] - dz) + ' / minimum redshift\n')
fo.write('REDMAX = ' + str(line['Z'] + dz) + ' / maximum redshift\n')
dz = dz_from_dv(dv, line['Z'])
fo.write('REDMIN = ' + str(line['Z'] - dz) + ' / minimum redshift for line fit\n')
fo.write('REDMAX = ' + str(line['Z'] + dz) + ' / maximum redshift for line fit\n')
else:
fo.write('REDMIN = ' + str(line['Z'] - 0.1) + ' / minimum redshift\n')
fo.write('REDMAX = ' + str(line['Z'] + 0.1) + ' / maximum redshift\n')
logger.error('dv must be not None and non-zero')
if deltav is not None:
deltaz = dz_from_dv(deltav, line['Z'])
fo.write('ZRMIN = ' + str(line['Z'] - deltaz) + ' / minimum redshift for spectral range around lines \n')
fo.write('ZRMAX = ' + str(line['Z'] + deltaz) + ' / maximum redshift for spectral range around lines \n')
else:
if dv is not None:
deltaz = 2 * dz
else:
deltaz = None
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('INITW = ' + str(initw) + ' / initial line width in km/s during the line fitting process\n')
fo.write('WMIN = ' + str(wmin) + ' / minimum line width in km/s\n')
fo.write('WMAX = ' + str(wmax) + ' / maximum line width 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('DGCTNUM = ' + str(degcont) + ' / degree 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')
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)
......@@ -211,9 +241,9 @@ def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=T
fo.write('XMAX = None / X maximum value for spatial cut of the input cube (pixel)\n')
fo.write('YMIN = None / Y minimum value for spatial cut of the input cube (pixel)\n')
fo.write('YMAX = None / Y maximum value for spatial cut of the input cube (pixel)\n')
if dv is not None:
minlbda = np.min(np.array(lbda)) * (1 + line['Z'] - 1.5 * dz) # XXX VERIFIER QUE CAMEL FAIT CA POUR DETERMINER LE LBDA RANGE AUTOURS DES RAIES: PAS TOUT A FAIT. CAMEL FAIT ~ 2 x dz. Il faudrait faire la même chose. En fait idéalmement, zmin, zmax devraient donner directement le bon résusltat, sans facteur 1.5 ou 2...
maxlbda = np.max(np.array(lbda)) * (1 + line['Z'] + 1.5 * dz)
if deltaz is not None:
minlbda = np.min(np.array(lbda)) * (1 + line['Z'] - deltaz)
maxlbda = np.max(np.array(lbda)) * (1 + line['Z'] + deltaz)
zmin = np.round((minlbda - crval) / cdelt + crpix - 1)
zmax = np.round((maxlbda - crval) / cdelt + crpix - 1)
fo.write('ZMIN = ' + str(zmin) + ' / Z minimum value for spectral cut of the input cube (spectral channel, integer)\n')
......@@ -254,9 +284,9 @@ def main(argv):
#parser.add_argument('--dxy', action="store", dest="dxy", default='15', help="size in pixels around the center")
parser.add_argument('--commonw', action="store", dest="commonw", default='True', help="parameter to use a common width in velocity for the various lines")
parser.add_argument('--dv', action="store", dest="dv", default='1000.', help="velocity difference in km/s to determine spectral range considered around lines")
parser.add_argument('--dv', action="store", dest="dv", default='500.', help="velocity difference in km/s to determine spectral range considered around lines")
parser.add_argument('--dxy', action="store", dest="dxy", default='15', help="size in pixels around the center")
parser.add_argument('--deltav', action="store", dest="deltav", default='400.', help="velocity range in km/s to determine redshift range for fit")
parser.add_argument('--deltav', action="store", dest="deltav", default='2000.', help="velocity range in km/s to determine redshift range for fit")
parser.add_argument('--initw', action="store", dest="initw", default='50.', help="initial guess for the line width (dispersion) in km/s")
parser.add_argument('--wmin', action="store", dest="wmin", default='30.', help="minimum width allowed in km/s")
parser.add_argument('--wmax', action="store", dest="wmax", default='250.', help="maximum width allowed in km/s")
......@@ -280,7 +310,6 @@ def main(argv):
if opts['varfile'] is None:
opts['varfile'] = opts['cubefile']
#create_config(opts['pathb'], opts['cubefile'], opts['varfile'], opts['catfile'], opts['lines'], opts['suffixeout'], dv=float(opts['dv']), dxy=int(opts['dxy']))
create_config(opts['pathb'], opts['cubefile'], opts['varfile'], opts['catfile'], opts['lines'], opts['suffixeout'],
dv=float(opts['dv']), dxy=int(opts['dxy']),
deltav=float(opts['deltav']), dfit=float(opts['dfit']), commonw=(opts['commonw'] == 'True'),
......@@ -290,9 +319,10 @@ def main(argv):
wsmooth=int(opts['wsmooth']), ssmooth=float(opts['ssmooth'])
)
# Faire un REDMIN, REDMAX (VELMIN/VELMAX, VLOTH/VUPTH, ...) pour les pixels spectraux à considérer autours des raies et utiliser un VMIN, VMAX ou DV pour les bornes du fit. Du coup, y aurait plus besoin de faire la conversion DV -> Redmin, Redmax dans create_config.
#Eventuellement faire une condition dans CAMEL: si on a qu'un des deux paramètres (REDMIN/REDMAX ou VMIN/VMAX ou peu importe le choix de nom final) faire une condition telle qu'actuellement
# Du coup, il vaut probablement mieux laisser un REDMIN/REDMAX relié au DV (pour que ça marche encore avec les anciennes versions des programmes/config files) et définir un nouveau VELMIN, VELMAX pour le bornes spectrales
#create_config(opts['pathb'], opts['cubefile'], opts['varfile'], opts['catfile'], opts['lines'], opts['suffixeout'], dv=float(opts['dv']), dxy=int(opts['dxy']))
# XXX Faire un REDMIN, REDMAX (VELMIN/VELMAX, VLOTH/VUPTH, ...) pour les pixels spectraux à considérer autours des raies et utiliser un VMIN, VMAX ou DV pour les bornes du fit. Du coup, y aurait plus besoin de faire la conversion DV -> Redmin, Redmax dans create_config.
# XXX Eventuellement faire une condition dans CAMEL: si on a qu'un des deux paramètres (REDMIN/REDMAX ou VMIN/VMAX ou peu importe le choix de nom final) faire une condition telle qu'actuellement
# XXX Du coup, il vaut probablement mieux laisser un REDMIN/REDMAX relié au DV (pour que ça marche encore avec les anciennes versions des programmes/config files) et définir un nouveau ZRMIN, ZRMAX pour le bornes spectrales
if __name__ == "__main__":
main(sys.argv[1:])
......
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import os, sys
import argparse
import ipdb
import numpy as np
#import pyfits as pf
import astropy.io.fits as fits
import astropy.constants as ct
import astropy.wcs as wcs
import astropy.coordinates as coord
import astropy.units as u
import logging
logging.basicConfig(level=logging.DEBUG)
logger = logging.getLogger('CAMEL')
def dz_from_dv(dv, z):
'''
This function enables to compute the redshift interval corresponding to a given velocity interval
Parameters
----------
dv: float
velocity interval in km/s
z: float
mean redshift of the source
Returns
-------
dz: float
corresponding redshift interval
'''
dz = dv / ct.c.to('km/s').value * (1 + z)
return dz
def create_config(path, cubefile, varfile, catfile, lines, suffixeout, commonw=True, dv=500., dxy=15, deltav=2000., 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 configuration file name
commonw: bool
specifies if a common width has to be used for the various lines
dv: float
velocity range in km/s to determine redshift range for fit
dxy: int
spatial extent around the object for cutting the cube (in pixels)
deltav: float
interval of velocity around the systemic redshift (from the catalogue) to search for lines and for cutting the cube spectrally
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 box for the sigma-clipping (in pixels)
nclip: int
number of passes for the sigma-clipping
wsmooth: int
spectral smooth method
ssmooth: float
Gaussian spatial smooth FWHM in pixels
'''
#ipdb.set_trace()
cat = np.genfromtxt(path + catfile, names=True, case_sensitive=False)
try:
hdul = fits.open(path + cubefile)
logger.info('Using cube %s ' % (path + cubefile) )
except:
logger.info('Not able to read cube %s !' % (path + cubefile) )
return
try:
cube = hdul[1].data
hdr = hdul[1].header
logger.info('Using extension 1 of cube %s ' % (path + cubefile))
except:
try:
cube = hdul[0].data
hdr = hdul[0].header
logger.info('Using extension 0 of cube %s ' % (path + cubefile))
except:
logger.info('Problem in data stored in extensions 0 and 1 of cube %s ' % (path + cubefile))
cunit = hdr['CUNIT3']
crpix = hdr['CRPIX3']
crval = hdr['CRVAL3']
cdelt = hdr['CD3_3']
try:
cdelt = hdr['CD3_3']
except:
cdelt = hdr['CDELT3']
logger.debug(" wcs check: %s %d %.4f %.6f " % (cunit,crpix,crval,cdelt))
if 'CTYPE3' in hdr.keys(): del hdr['CTYPE3']
if 'CRVAL3' in hdr.keys(): del hdr['CRVAL3']
if 'CRPIX3' in hdr.keys(): del hdr['CRPIX3']
if 'CUNIT3' in hdr.keys(): del hdr['CUNIT3']
if 'CDELT3' in hdr.keys(): del hdr['CDELT3']
if 'CD3_3' in hdr.keys(): del hdr['CD3_3']
if 'CD3_2' in hdr.keys(): del hdr['CD3_2']
if 'CD3_1' in hdr.keys(): del hdr['CD3_1']
if 'CD2_3' in hdr.keys(): del hdr['CD2_3']
if 'CD1_3' in hdr.keys(): del hdr['CD1_3']
if '' in hdr.keys(): del hdr['']
if 'COMMENT' in hdr.keys(): del hdr['COMMENT']
if 'HISTORY' in hdr.keys(): del hdr['HISTORY']
wlim = np.zeros(np.shape(cube)[1:])
hdu = fits.PrimaryHDU(data=wlim, header=hdr)
hdul1 = fits.HDUList(hdu)
hdr = hdul1[0].header
w = wcs.WCS(hdr, hdul1)
#w = wcs.WCS(hdr[:28], hdul1)
logger.debug(' Naxis %d ' % (w.naxis))
lc = np.size(cat)
if lc == 1: cat = [cat]
logger.debug(" catatalog length %s " % (str(lc)) )
for line in cat:
obj = int(line['ID'])
suffixe = suffixeout + '_' + str(obj) + '_' + lines
fileout = suffixe + '/' + suffixe + '.config'
logger.info('Writing config ' + fileout)
print(path + fileout)
fo = open(path + fileout, 'w')
fo.write('FITSFILE= ' + '../' + cubefile + ' / path of the fits file containing the cube\n')
fo.write('OUTPUT = ' + 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')
if dv is not None:
dz = dz_from_dv(dv, line['Z'])
fo.write('REDMIN = ' + str(line['Z'] - dz) + ' / minimum redshift for line fit\n')
fo.write('REDMAX = ' + str(line['Z'] + dz) + ' / maximum redshift for line fit\n')
else:
fo.write('REDMIN = ' + str(line['Z'] - 0.1) + ' / minimum redshift\n')
fo.write('REDMAX = ' + str(line['Z'] + 0.1) + ' / maximum redshift\n')
logger.error('dv must be not None and non-zero')
if deltav is not None:
deltaz = dz_from_dv(deltav, line['Z'])
fo.write('ZRMIN = ' + str(line['Z'] - deltaz) + ' / minimum redshift for spectral range around lines \n')
fo.write('ZRMAX = ' + str(line['Z'] + deltaz) + ' / maximum redshift for spectral range around lines \n')
else:
if dv is not None:
deltaz = 2 * dz
else:
deltaz = None
fo.write('INITW = ' + str(initw) + ' / initial line width in km/s during the line fitting process\n')
fo.write('WMIN = ' + str(wmin) + ' / minimum line width in km/s\n')
fo.write('WMAX = ' + str(wmax) + ' / maximum line width in km/s\n')
fo.write('DFIT = ' + str(dfit) + ' / bin for fitting the line center in km/s\n')
fo.write('DGCTNUM = ' + str(degcont) + ' / degree 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)
logger.debug(' converting ra,dec %f,%f ' % (line['RA'],line['DEC']) )
xc, yc = np.round(w.wcs_world2pix(ra, dec, 0))
#xc, yc, _ = np.round(w.wcs_world2pix(ra, dec, 0, 0))
logger.debug(' converting x,y %f %f ' % (xc,yc) )
if dxy is not None:
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')
else:
fo.write('XMIN = None / X minimum value for spatial cut of the input cube (pixel)\n')
fo.write('XMAX = None / X maximum value for spatial cut of the input cube (pixel)\n')
fo.write('YMIN = None / Y minimum value for spatial cut of the input cube (pixel)\n')
fo.write('YMAX = None / Y maximum value for spatial cut of the input cube (pixel)\n')
if deltaz is not None:
minlbda = np.min(np.array(lbda)) * (1 + line['Z'] - deltaz)
maxlbda = np.max(np.array(lbda)) * (1 + line['Z'] + deltaz)
zmin = np.round((minlbda - crval) / cdelt + crpix - 1)
zmax = np.round((maxlbda - crval) / cdelt + crpix - 1)
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)')
else:
fo.write('ZMIN = None / Z minimum value for spectral cut of the input cube (spectral channel, integer)\n')
fo.write('ZMAX = None / Z maximum value for spectral cut of the input cube (spectral channel, integer)')
fo.close()
def main_byhand():
gr = '116'
#gr = '83'
#gr = '28'
#gr = '32'
lines = 'o3hb'
lines = 'o2'
suffixeout = 'gr' + gr
path = '/home/bepinat/Instruments/MUSE/groups/kinematics/gr' + gr + '/' + lines + '/'
cubefile = '../DATACUBE_FINAL_ZAP_COSMOS-GR' + gr + '_newWCS.fits'
varfile = cubefile
catfile = '../catalog_gr' + gr + '.txt'
create_config(path, cubefile, varfile, catfile, lines, suffixeout)
def main(argv):
parser = argparse.ArgumentParser()
parser.add_argument('--lines', '-l', action="store", dest="lines", default="", help="name of the lines")
parser.add_argument('--suffixe', '-s', action="store", dest="suffixeout", default="", help="suffixe")
parser.add_argument('--path', '-p', action="store", dest="pathb", default="", help="initial path")
parser.add_argument('--cat', '-c', action="store", dest="catfile", default="", help="path of the catalogue with respect to initial path")
parser.add_argument('--data', '-d', action="store", dest="cubefile", default="", help="path of the cube with respect to initial path")
parser.add_argument('--var', '-v', action="store", dest="varfile", default=None, help="path of the variance with respect to initial path")
#parser.add_argument('--dv', action="store", dest="dv", default='1000.', help="velocity difference in km/s to determine redshift range for fit")
#parser.add_argument('--dxy', action="store", dest="dxy", default='15', help="size in pixels around the center")
parser.add_argument('--commonw', action="store", dest="commonw", default='True', help="parameter to use a common width in velocity for the various lines")
parser.add_argument('--dv', action="store", dest="dv", default='500.', help="velocity difference in km/s to determine spectral range considered around lines")
parser.add_argument('--dxy', action="store", dest="dxy", default='15', help="size in pixels around the center")
parser.add_argument('--deltav', action="store", dest="deltav", default='2000.', help="velocity range in km/s to determine redshift range for fit")
parser.add_argument('--initw', action="store", dest="initw", default='50.', help="initial guess for the line width (dispersion) in km/s")
parser.add_argument('--wmin', action="store", dest="wmin", default='30.', help="minimum width allowed in km/s")
parser.add_argument('--wmax', action="store", dest="wmax", default='250.', help="maximum width allowed in km/s")
parser.add_argument('--dfit', action="store", dest="dfit", default='100.', help="interval step for fitting the line position in km/s. Large values speed up the computation but could miss the true position")
parser.add_argument('--degcont', action="store", dest="degcont", default='0', help="degree of the polynomial used to fit the continuum")
parser.add_argument('--sclip', action="store", dest="sclip", default='10', help="sigma clipping threshold (sclip * sigma) to clean the cube")
parser.add_argument('--xyclip', action="store", dest="xyclip", default='3', help="size of the local box for the sigma-clipping (in pixels)")
parser.add_argument('--nclip', action="store", dest="nclip", default='3', help="number of passes for the sigma-clipping")
parser.add_argument('--wsmooth', action="store", dest="wsmooth", default='0', help="spectral smooth method")
parser.add_argument('--ssmooth', action="store", dest="ssmooth", default='2.', help="Gaussian spatial smooth FWHM in pixels")
try:
options = parser.parse_args()
except:
sys.exit(2)
opts = vars(options) # this is a dictionnary
#print(opts)
if opts['varfile'] is None:
opts['varfile'] = opts['cubefile']
create_config(opts['pathb'], opts['cubefile'], opts['varfile'], opts['catfile'], opts['lines'], opts['suffixeout'],
dv=float(opts['dv']), dxy=int(opts['dxy']),
deltav=float(opts['deltav']), dfit=float(opts['dfit']), commonw=(opts['commonw'] == 'True'),
initw=float(opts['initw']), wmin=float(opts['wmin']), wmax=float(opts['wmax']),
degcont=int(opts['degcont']),
sclip=int(opts['sclip']), xyclip=int(opts['xyclip']), nclip=int(opts['nclip']),
wsmooth=int(opts['wsmooth']), ssmooth=float(opts['ssmooth'])
)
#create_config(opts['pathb'], opts['cubefile'], opts['varfile'], opts['catfile'], opts['lines'], opts['suffixeout'], dv=float(opts['dv']), dxy=int(opts['dxy']))
# XXX Faire un REDMIN, REDMAX (VELMIN/VELMAX, VLOTH/VUPTH, ...) pour les pixels spectraux à considérer autours des raies et utiliser un VMIN, VMAX ou DV pour les bornes du fit. Du coup, y aurait plus besoin de faire la conversion DV -> Redmin, Redmax dans create_config.
# XXX Eventuellement faire une condition dans CAMEL: si on a qu'un des deux paramètres (REDMIN/REDMAX ou VMIN/VMAX ou peu importe le choix de nom final) faire une condition telle qu'actuellement