Test cases from Simulations 140425¶
Author: Yannick Copin y.copin@ipnl.in2p3.fr
[1]:
%matplotlib notebook
import numpy as N
import matplotlib.pyplot as P
%load_ext autoreload
%autoreload 2
Simulated spectra¶
[2]:
from inspec.sim140425 import read_sim_input, fig_hist_survey
datapath = "../inspec/data/"
simpath = os.path.join(datapath, "Simulations_140425/Uncontaminated/")
filename = os.path.join(simpath, "uncontaminated_spectra.dat")
metadata = read_sim_input(filename)
colnames = metadata.dtype.names
../inspec/data/Simulations_140425/Uncontaminated/uncontaminated_spectra.dat: 7432 spectra (columns: ['idx', 'type', 'z', 'fHa', 'ewHa', 'mag'])
[3]:
fig = fig_hist_survey(metadata,
title="%s (%d entries)" % (
os.path.basename(filename), len(metadata)))
fig.tight_layout()
Adpative continuum + emission line fit¶
Read emission lines¶
[4]:
from inspec import tools
emilines = tools.parse_json(os.path.join(datapath, "emissionLines.json"))
emilines
[4]:
OrderedDict([('OVI', 1033.3),
('Lyα', 1215.67),
('NV', 1239.42),
('OI', 1305.53),
('CIV', 1545.86),
('CIII', 1908.27),
('CII', 2326.0),
('MgII', 2800.32),
('OII_', 3727.092),
('OII', 3729.874),
('Hε', 3971.195),
('Hδ', 4102.892),
('Hγ', 4341.684),
('Hβ', 4862.683),
('OIII_', 4960.295),
('OIII', 5008.24),
('NII_', 6549.86),
('Hα', 6564.613),
('NII', 6585.277),
('SII_', 6718.294),
('SII', 6732.673),
('Paα', 18756.13),
('Paβ', 12821.59),
('Paγ', 10941.09),
('Paδ', 10052.13)])
Adaptive fit¶
[5]:
from inspec.spectrum import SimSpectrum, fit_adaptive, plot_specs
sample = [10568, 15944] # ../data/Simulations_140425/
for i, idx in enumerate(sample, start=1):
print(" {}/{}: ID {} ".format(i, len(sample), idx).center(50, '-'))
# Reading
sdict = SimSpectrum.from_idx(idx, path=os.path.join(simpath, "*/"))
# Merging
merged = SimSpectrum.merge_specs(sdict, rtol=None)
subdata = metadata[N.searchsorted(metadata['idx'], idx)]
# merged.writeto("spectrum_%d_merged.fits" % idx,
# keywords=dict(zip(colnames, subdata)))
# Polynomial continuum + emission lines
cont, lines = fit_adaptive(merged, pmax=0.05, dmax=3, nmax=5, verbose=2)
# Plot
fig = plot_specs(sdict,
model=(cont, lines),
metadata=subdata,
emilines=emilines)
fig.tight_layout()
----------------- 1/2: ID 10568 ------------------
Spectrum #10568 ('../inspec/data/Simulations_140425/Uncontaminated/Survey'): 5 rolls
Ref. roll180_01: RMS=1.45766, 659 px [12032.00-18480.40] step=9.80 AA/px, mstart=7.40 AA
roll000_01: RMS=0.439268, 69 px [17321.40-17987.80] step=9.80 AA/px, mstart=4.79 AA
roll000_02: RMS=1.56681, 491 px [12509.60-17311.60] step=9.80 AA/px, mstart=4.80 AA
roll090_01: RMS=1.58077, 560 px [12509.60-17987.80] step=9.80 AA/px, mstart=4.80 AA
roll090_02: RMS=2.2463, 256 px [12509.60-15008.60] step=9.80 AA/px, mstart=4.80 AA
Common wavelength domain: 659 px [12032.00, 18480.40]
Merging 5/5 spectra
Optimization terminated successfully.
Current function value: 1311.767300
Iterations: 23
Function evaluations: 46
Testing Degree=0, Chi2=17981.67559874864, DoF=658, p-value=0
dChi2=-inf, dDoF=1, p-value=0
Optimization terminated successfully.
Current function value: 1166.464711
Iterations: 87
Function evaluations: 163
Testing Degree=1, Chi2=17040.54017360229, DoF=657, p-value=0
dChi2=-941.14, dDoF=1, p-value=1.12123e-206
Optimization terminated successfully.
Current function value: 1178.722095
Iterations: 79
Function evaluations: 154
Testing Degree=2, Chi2=17189.95358729855, DoF=656, p-value=0
dChi2=149.41, dDoF=1, p-value=1
Warning: Maximum number of iterations has been exceeded.
Testing
WARNING: Wavelength vector x is not strictly uniformly sampled: steps=[9.7988 9.8008] A
WARNING: Wavelength vector x is not strictly uniformly sampled: steps=[9.7988 9.7998 9.8008] A
WARNING: Wavelength vector x is not strictly uniformly sampled: steps=[9.7998 9.8008] A
WARNING: Model is linear in parameters; consider using linear fitting methods. [astropy.modeling.fitting]
WARNING: The fit may be unsuccessful; Maximum number of iterations reached. [astropy.modeling.optimizers]
Degree=3, Chi2=17364.00280578613, DoF=655, p-value=0
dChi2=323.46, dDoF=2, p-value=1
Converged to Degree=1, Chi2=17040.54017360229, DoF=657, p-value=0
Testing line a=13.883441261226292, m=12796.3994140625, s=7.432934115273265
Tested #lines=1, Chi2=637.4507745270935, DoF=656, p-value=0.690956
dChi2=-16403.09, dDoF=3, p-value=0
Testing line a=0.9112960900008898, m=12855.19921875, s=17.966791624381223
Tested #lines=2, Chi2=594.2620262919658, DoF=653, p-value=0.951364
dChi2=-43.19, dDoF=3, p-value=2.24407e-09
Testing line a=1.084388576651331, m=12786.599609375, s=18.27227218054737
Tested #lines=3, Chi2=593.8608459511294, DoF=650, p-value=0.943587
dChi2=-0.40, dDoF=3, p-value=0.939999
Converged to #lines=2, Chi2=594.2620262919658, DoF=653, p-value=0.951364
Model: CompoundModel
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Expression: [0] + [1]
Components:
[0]: <Gaussian1D(amplitude=13.30734502, mean=12799.74298338, stddev=28.9430113)>
[1]: <Gaussian1D(amplitude=0.90788934, mean=12858.581516, stddev=12.15996857)>
Parameters:
amplitude_0 mean_0 ... mean_1 stddev_1
------------------ ------------------ ... ------------------ ------------------
13.307345017705106 12799.742983383068 ... 12858.581515997343 12.159968571698137
----------------- 2/2: ID 15944 ------------------
Spectrum #15944 ('../inspec/data/Simulations_140425/Uncontaminated/Survey'): 4 rolls
Ref. roll180_01: RMS=1.1362, 659 px [12032.00-18480.40] step=9.80 AA/px, mstart=7.40 AA
roll000_01: RMS=0.975714, 560 px [12509.60-17987.80] step=9.80 AA/px, mstart=4.80 AA
roll090_01: RMS=1.01554, 560 px [12509.60-17987.80] step=9.80 AA/px, mstart=4.80 AA
roll090_02: RMS=1.01746, 560 px [12509.60-17987.80] step=9.80 AA/px, mstart=4.80 AA
Common wavelength domain: 659 px [12032.00, 18480.40]
Merging 4/4 spectra
Optimization terminated successfully.
Current function value: 1071.828951
Iterations: 22
Function evaluations: 44
Testing Degree=0, Chi2=14894.682597368273, DoF=658, p-value=0
dChi2=-inf, dDoF=1, p-value=0
Optimization terminated successfully.
Current function value: 1062.852972
Iterations: 93
Function evaluations: 177
Testing Degree=1, Chi2=14897.403326376869, DoF=657, p-value=0
dChi2=2.72, dDoF=1, p-value=1
Optimization terminated successfully.
Current function value: 1064.092096
Iterations: 72
Function evaluations: 137
Testing Degree=2, Chi2=14889.872785446596, DoF=656, p-value=0
dChi2=-4.81, dDoF=2, p-value=0.090274
Warning: Maximum number of iterations has been exceeded.
Testing Degree=3, Chi2=14892.02983877401, DoF=655, p-value=0
dChi2=-2.65, dDoF=3, p-value=0.448315
Converged to Degree=0, Chi2=14894.682597368273, DoF=658, p-value=0
Testing line a=11.285707393448783, m=16099.0, s=9.339131843883932
WARNING: Wavelength vector x is not strictly uniformly sampled: steps=[9.7988 9.7998 9.8008] A
WARNING: Model is linear in parameters; consider using linear fitting methods. [astropy.modeling.fitting]
WARNING: The fit may be unsuccessful; Maximum number of iterations reached. [astropy.modeling.optimizers]
Tested #lines=1, Chi2=1255.6048529758257, DoF=656, p-value=4.52575e-40
dChi2=-13639.08, dDoF=3, p-value=0
Testing line a=7.4988382652748475, m=12286.7998046875, s=9.425678019673372
Tested #lines=2, Chi2=545.3257155255861, DoF=653, p-value=0.999162
dChi2=-710.28, dDoF=3, p-value=1.23911e-153
Testing line a=0.8705252054598273, m=16157.7998046875, s=13.678936144733564
Tested #lines=3, Chi2=452.15610482457, DoF=650, p-value=1
dChi2=-93.17, dDoF=3, p-value=4.56711e-20
Testing line a=2.2060126588484517, m=12169.19921875, s=9.93845730755947
Tested #lines=4, Chi2=397.1731482319833, DoF=647, p-value=1
dChi2=-54.98, dDoF=3, p-value=6.9239e-12
Testing line a=0.3269342934967913, m=16040.2001953125, s=21.530104115246015
Tested #lines=5, Chi2=385.9701091869159, DoF=644, p-value=1
dChi2=-11.20, dDoF=3, p-value=0.0106771
Converged to #lines=5, Chi2=385.9701091869159, DoF=644, p-value=1
Model: CompoundModel
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Expression: [0] + [1] + [2] + [3] + [4]
Components:
[0]: <Gaussian1D(amplitude=11.60717741, mean=16100.14108505, stddev=19.38628273)>
[1]: <Gaussian1D(amplitude=7.84228378, mean=12288.27870412, stddev=18.97187128)>
[2]: <Gaussian1D(amplitude=1.05731839, mean=16155.2816853, stddev=10.66180227)>
[3]: <Gaussian1D(amplitude=1.84166665, mean=12175.80335579, stddev=30.62124566)>
[4]: <Gaussian1D(amplitude=0.35534719, mean=16036.24030334, stddev=10.95307651)>
Parameters:
amplitude_0 mean_0 ... mean_4 stddev_4
------------------ ----------------- ... ----------------- ----------------
11.607177406331902 16100.14108504782 ... 16036.24030334245 10.9530765083286
Synthetic magnitudes¶
[6]:
from inspec import synth
print(merged)
band = synth.Filter.bandpass("test", 14000, 16000)
print(band)
f, df = band.integrate(merged) # Flux
f /= 1e17 # Scale for the 1e-17 flux normalization
df /= 1e17
print("Flux in {!r}: {:g} +/- {:g}".format(band.name, f, df))
m, dm = band.magn(merged) # AB-magnitude
m += 17 * 2.5 # Offset for the 1e-17 flux normalization
print("AB-magn in {!r}: {:f} +/- {:f}".format(band.name, m, dm))
merged: RMS=1.12121, 659 px [12032.00-18480.40] step=9.80 AA/px, mstart=7.40 AA
Bandpass 'test' [14000.0,16000.0],
Flux in 'test': 0.00218547 +/- 0.00016646
AB-magn in 'test': 21.321959 +/- 0.082697
Continuum indices¶
[7]:
cont_indices = tools.parse_json(os.path.join(datapath, "contIndex.json"))
cont_indices
[7]:
OrderedDict([('B2640',
OrderedDict([('type', 'flambda ratio'),
('blue', [2600.0, 2630.0]),
('red', [2645.0, 2675.0]),
('reference', 'Spinrad+ 1997ApJ...484..581S')])),
('Mg_UV',
OrderedDict([('type', 'flambda normed'),
('line', [2625.0, 2725.0]),
('blue', [2525.0, 2625.0]),
('red', [2725.0, 2825.0]),
('reference', 'Daddi+ 2005ApJ...626..680D')])),
('B2900',
OrderedDict([('type', 'flambda ratio'),
('blue', [2855.0, 2885.0]),
('red', [2915.0, 2945.0]),
('reference', 'Spinrad+ 1997ApJ...484..581S')])),
('D4000 wide',
OrderedDict([('type', 'D4000'),
('blue', [3750.0, 3950.0]),
('red', [4050.0, 4250.0]),
('reference',
'Bruzual 1983ApJ...273..105B, Hamilton 1985ApJ...297..371H')])),
('D4000 narrow',
OrderedDict([('type', 'D4000'),
('blue', [3850.0, 3950.0]),
('red', [4000.0, 4100.0]),
('reference', 'Balogh+ 1999ApJ...527...54B')]))])
[8]:
from inspec import indices
indices.flambda_ratio(cont_indices['B2900'], merged, redshift=4)
[8]:
(0.9848834320449218, 0.36849523356860564)
[9]:
indices.flambda_normed(cont_indices['Mg_UV'], merged, redshift=4)
[9]:
(0.7867401890423812, 0.16767674824162523)
[10]:
indices.D4000_ratio(cont_indices['D4000 wide'], merged, redshift=3)
[10]:
(1.0179880020051273, 0.18620092550913506)
Emission line fit¶
[11]:
import astropy.modeling as AM
from inspec.lines import ComplexHa, DoubletOIII
fitter = AM.fitting.LevMarLSQFitter()
Single emission line¶
[12]:
# Compound model: [0] will be the line, [1] will be the continuum
Ha = AM.models.Gaussian1D(mean=16100, amplitude=10, stddev=20) + AM.polynomial.Legendre1D(0)
print(Ha)
Model: CompoundModel
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Expression: [0] + [1]
Components:
[0]: <Gaussian1D(amplitude=10., mean=16100., stddev=20.)>
[1]: <Legendre1D(0, c0=0.)>
Parameters:
amplitude_0 mean_0 stddev_0 c0_1
----------- ------- -------- ----
10.0 16100.0 20.0 0.0
[13]:
select = (merged.x > 15950) & (merged.x < 16250) # Adjusted spectral domain
ha = fitter(Ha, merged.x[select], merged.y[select], weights=1/merged.v[select])
print(ha)
param_cov = fitter.fit_info['param_cov'] # Adjusted parameter covariance matrix
dparams = N.sqrt(N.diagonal(param_cov))
param_corr = param_cov / N.outer(dparams, dparams)
print("Estimated redshift: {} ± {}".format(
ha[0].mean / emilines[u'Hα'] - 1,
dparams[0] / emilines[u'Hα']))
print("Correlation matrix:\n", param_corr)
Model: CompoundModel
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Expression: [0] + [1]
Components:
[0]: <Gaussian1D(amplitude=11.57675934, mean=16100.11135049, stddev=18.99123589)>
[1]: <Legendre1D(0, c0=0.26317098)>
Parameters:
amplitude_0 mean_0 stddev_0 c0_1
------------------ ------------------ ------------------ -------------------
11.576759344364786 16100.111350488007 18.991235887142437 0.26317097660655775
Estimated redshift: 1.452560623221507 ± 5.2216786456211464e-05
Correlation matrix:
[[ 1. 0.0358965 -0.53530393 -0.10418745]
[ 0.0358965 1. -0.04707842 0.0017797 ]
[-0.53530393 -0.04707842 1. -0.34731023]
[-0.10418745 0.0017797 -0.34731023 1. ]]
[14]:
fig = P.figure(figsize=(14, 6))
ax = fig.add_subplot(1, 1, 1,
xlabel=u"Wavelength [Å]",
ylabel="Flux",
title=u"Hα emission line")
l, = ax.plot(merged.x, merged.y, lw=2, label="Data")
ax.fill_between(merged.x, merged.y - N.sqrt(merged.v), merged.y + N.sqrt(merged.v),
facecolor=l.get_color(), edgecolor='none', alpha=0.3, label='_')
ax.plot(merged.x, Ha(merged.x), ls='--', label="Guess")
l, = ax.plot(merged.x, ha(merged.x), lw=2, label="Fit")
ax.plot(merged.x, ha[1](merged.x), ls=':', label="Continuum")
ax.set_xlim(ha[0].mean - 10 * ha[0].stddev, ha[0].mean + 10 * ha[0].stddev)
ax.axvspan(15950, 16250, fc='0.8', ec='none', alpha=0.3, label="Adjusted domain")
ax.legend()
[14]:
<matplotlib.legend.Legend at 0x7f68f1243550>
[NII] + Hα complex¶
[15]:
# Compound model: [0] will be the ComplexHa, [1] will be the continuum
complexHa = ComplexHa(wHa=16100, iHa=10, iNII=3, sigma=20) + AM.polynomial.Legendre1D(0)
print(complexHa)
Model: CompoundModel
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Expression: [0] + [1]
Components:
[0]: <ComplexHa(wHa=16100., iHa=10., iNII=3., sigma=20.)>
[1]: <Legendre1D(0, c0=0.)>
Parameters:
wHa_0 iHa_0 iNII_0 sigma_0 c0_1
------- ----- ------ ------- ----
16100.0 10.0 3.0 20.0 0.0
[16]:
select = (merged.x > 15950) & (merged.x < 16250) # Adjusted spectral domain
halpha = fitter(complexHa, merged.x[select], merged.y[select], weights=1/merged.v[select])
print(halpha)
param_cov = fitter.fit_info['param_cov'] # Adjusted parameter covariance matrix
dparams = N.sqrt(N.diagonal(param_cov))
param_corr = param_cov / N.outer(dparams, dparams)
print("Estimated redshift: {} ± {}".format(
halpha[0].wHa / ComplexHa.Halpha - 1,
dparams[0] / ComplexHa.Halpha))
print("Correlation matrix:\n", param_corr)
Model: CompoundModel
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Expression: [0] + [1]
Components:
[0]: <ComplexHa(wHa=16099.57175802, iHa=12.04960419, iNII=1.6035007, sigma=17.44276748)>
[1]: <Legendre1D(0, c0=0.12257132)>
Parameters:
wHa_0 iHa_0 ... sigma_0 c0_1
------------------ ------------------ ... ------------------ -------------------
16099.571758017662 12.049604192250257 ... 17.442767482400615 0.12257132382382573
Estimated redshift: 1.4524780524822423 ± 3.235282882778465e-05
Correlation matrix:
[[ 1. 0.02630667 -0.14262952 0.02388155 0.04891581]
[ 0.02630667 1. 0.23518973 -0.58280689 -0.14458737]
[-0.14262952 0.23518973 1. -0.46103498 -0.39224759]
[ 0.02388155 -0.58280689 -0.46103498 1. -0.12664899]
[ 0.04891581 -0.14458737 -0.39224759 -0.12664899 1. ]]
[17]:
fig = P.figure(figsize=(14, 6))
ax = fig.add_subplot(1, 1, 1,
xlabel=u"Wavelength [Å]",
ylabel="Flux",
title=u"[NII] + Hα complex")
l, = ax.plot(merged.x, merged.y, lw=2, label="Data")
ax.fill_between(merged.x, merged.y - N.sqrt(merged.v), merged.y + N.sqrt(merged.v),
facecolor=l.get_color(), edgecolor='none', alpha=0.3, label='_')
ax.plot(merged.x, complexHa(merged.x), ls='--', label="Guess")
l, = ax.plot(merged.x, halpha(merged.x), lw=2, label="Fit")
ax.plot(merged.x, ComplexHa.Ha_line(merged.x, *halpha.parameters[:-1]), c=l.get_color(), ls=':', label=u"Hα")
ax.plot(merged.x, ComplexHa.NII_lines(merged.x, *halpha.parameters[:-1]), c=l.get_color(), ls='-.', label="[NII]")
ax.plot(merged.x, halpha[1](merged.x), ls=':', label="Continuum")
ax.set_xlim(halpha[0].wHa - 10 * halpha[0].sigma, halpha[0].wHa + 10 * halpha[0].sigma)
ax.axvspan(15950, 16250, fc='0.8', ec='none', alpha=0.3, label="Adjusted domain")
ax.legend()
[17]:
<matplotlib.legend.Legend at 0x7f68ef17c978>
[OIII] doublet¶
[18]:
# Compound model: [0] will be the DoubletOIII, [1] will be the continuum
doubletOIII = DoubletOIII(wOIII=12300, iOIII=10, sigma=20) + AM.polynomial.Legendre1D(0)
print(doubletOIII)
Model: CompoundModel
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Expression: [0] + [1]
Components:
[0]: <DoubletOIII(wOIII=12300., iOIII=10., sigma=20.)>
[1]: <Legendre1D(0, c0=0.)>
Parameters:
wOIII_0 iOIII_0 sigma_0 c0_1
------- ------- ------- ----
12300.0 10.0 20.0 0.0
[19]:
select = (merged.x > 12100) & (merged.x < 12400)
oiii = fitter(doubletOIII, merged.x[select], merged.y[select], weights=1/merged.v[select])
print(oiii)
param_cov = fitter.fit_info['param_cov'] # Adjusted parameter covariance matrix
dparams = N.sqrt(N.diagonal(param_cov))
param_corr = param_cov / N.outer(dparams, dparams)
print("Estimated redshift: {} ± {}".format(
oiii[0].wOIII / DoubletOIII.OIII - 1,
dparams[0] / DoubletOIII.OIII))
print("Correlation matrix:\n", param_corr)
Model: CompoundModel
Inputs: ('x',)
Outputs: ('y',)
Model set size: 1
Expression: [0] + [1]
Components:
[0]: <DoubletOIII(wOIII=12288.36871635, iOIII=10.21136293, sigma=18.56625818)>
[1]: <Legendre1D(0, c0=0.27190352)>
Parameters:
wOIII_0 iOIII_0 sigma_0 c0_1
------------------ ------------------ ------------------ -------------------
12288.368716349503 10.211362934778895 18.566258178496874 0.27190351971723187
Estimated redshift: 1.45363016076496 ± 0.00012726705320362686
Correlation matrix:
[[ 1. 0.04770397 -0.12709564 -0.00205968]
[ 0.04770397 1. -0.39140655 -0.24051853]
[-0.12709564 -0.39140655 1. -0.47517549]
[-0.00205968 -0.24051853 -0.47517549 1. ]]
[20]:
fig = P.figure(figsize=(14, 6))
ax = fig.add_subplot(1, 1, 1,
xlabel=u"Wavelength [Å]",
ylabel="Flux",
title="[OIII] doublet")
l, = ax.plot(merged.x, merged.y, lw=2, label="Data")
ax.fill_between(merged.x, merged.y - N.sqrt(merged.v), merged.y + N.sqrt(merged.v),
facecolor=l.get_color(), edgecolor='none', alpha=0.3, label='_')
ax.plot(merged.x, doubletOIII(merged.x), ls='--', label="Guess")
ax.plot(merged.x, oiii(merged.x), lw=2, label="Fit")
ax.plot(merged.x, oiii[1](merged.x), ls=':', label="Continuum")
ax.set_xlim(oiii[0].wOIII - 15 * oiii[0].sigma, oiii[0].wOIII + 10 * oiii[0].sigma)
ax.axvspan(12100, 12400, fc='0.8', ec='none', alpha=0.3, label="Adjusted domain")
ax.legend()
[20]:
<matplotlib.legend.Legend at 0x7f68ef148f28>
[ ]: