Newer
Older
function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
Giorgio Calderone
committed
Nspec = length(source.domain)
elapsed = time()
mzer = GFit.cmpfit()
Giorgio Calderone
committed
# Initialize components and guess initial values
Giorgio Calderone
committed
for id in 1:Nspec
λ = source.domain[id][:]
model = Model(source.domain[id],
:qso_cont => QSFit.qso_cont_component(TRecipe),
:Continuum => SumReducer([:qso_cont]))
Giorgio Calderone
committed
# TODO if source.options[:instr_broadening]
# TODO GFit.set_instr_response!(model, (l, f) -> instrumental_broadening(l, f, source.spectra[id].resolution))
# TODO end
Giorgio Calderone
committed
Giorgio Calderone
committed
c.norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(c.x0.val)
Giorgio Calderone
committed
end
for id in 1:Nspec
Giorgio Calderone
committed
λ = source.domain[id][:]
# Host galaxy template
if source.options[:use_host_template] &&
(minimum(λ) .< 5500 .< maximum(λ))
model[:galaxy] = QSFit.hostgalaxy(source.options[:host_template])
model[:Continuum] = SumReducer([:qso_cont, :galaxy])
# Split total flux between continuum and host galaxy
vv = Spline1D(λ, source.data[id].val, k=1, bc="error")(5500.)
model[:galaxy].norm.val = 1/2 * vv
model[:qso_cont].x0.val *= 1/2 * vv / Spline1D(λ, model(:qso_cont), k=1, bc="error")(5500.)
@patch! multi[id][:galaxy].norm = multi[ref_id][:galaxy].norm
Giorgio Calderone
committed
end
# Balmer continuum and pseudo-continuum
if source.options[:use_balmer]
tmp = [:qso_cont, :balmer]
(:galaxy in keys(model)) && push!(tmp, :galaxy)
model[:balmer] = QSFit.balmercont(0.1, 0.5)
model[:Continuum] = SumReducer(tmp)
c = model[:balmer]
Giorgio Calderone
committed
c.norm.val = 0.1
c.norm.fixed = false
c.norm.high = 0.5
c.ratio.val = 0.5
c.ratio.fixed = false
Giorgio Calderone
committed
c.ratio.high = 1
Giorgio Calderone
committed
end
end
bestfit = fit!(multi, source.data, minimizer=mzer); show(logio(source), bestfit)
Giorgio Calderone
committed
# QSO continuum renormalization
for id in 1:Nspec
model = multi[id]
freeze(model, :qso_cont)
c = model[:qso_cont]
Giorgio Calderone
committed
if c.norm.val > 0
println(logio(source), "$id: Cont. norm. (before): ", c.norm.val)
Giorgio Calderone
committed
while true
residuals = (model() - source.data[id].val) ./ source.data[id].unc
Giorgio Calderone
committed
(count(residuals .< 0) / length(residuals) > 0.9) && break
(c.norm.val < initialnorm / 5) && break # give up
Giorgio Calderone
committed
c.norm.val *= 0.99
Giorgio Calderone
committed
end
println(logio(source), "$id : Cont. norm. (after) : ", c.norm.val)
Giorgio Calderone
committed
else
println(logio(source), "$id: Skipping cont. renormalization")
Giorgio Calderone
committed
end
freeze(model, :qso_cont)
(:galaxy in keys(model)) && freeze(model, :galaxy)
(:balmer in keys(model)) && freeze(model, :balmer)
Giorgio Calderone
committed
end
Giorgio Calderone
committed
# Fit iron templates
Giorgio Calderone
committed
for id in 1:Nspec
Giorgio Calderone
committed
λ = source.domain[id][:]
iron_components = Vector{Symbol}()
if source.options[:use_ironuv]
fwhm = 3000.
source.options[:instr_broadening] || (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
comp = QSFit.ironuv(fwhm)
(_1, _2, coverage) = spectral_coverage(λ, source.spectra[id].resolution, comp)
threshold = get(source.options[:min_spectral_coverage], :ironuv, source.options[:min_spectral_coverage][:default])
if coverage >= threshold
Giorgio Calderone
committed
push!(iron_components, :ironuv)
else
println(logio(source), "Ignoring ironuv component on prediction $id (threshold: $threshold)")
Giorgio Calderone
committed
end
end
if source.options[:use_ironopt]
fwhm = 3000.
source.options[:instr_broadening] || (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
comp = QSFit.ironopt_broad(fwhm)
(_1, _2, coverage) = spectral_coverage(λ, source.spectra[id].resolution, comp)
threshold = get(source.options[:min_spectral_coverage], :ironopt, source.options[:min_spectral_coverage][:default])
if coverage >= threshold
fwhm = 500.
source.options[:instr_broadening] || (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
model[:ironoptbr] = comp
model[:ironoptna] = QSFit.ironopt_narrow(fwhm)
model[:ironoptbr].norm.val = 0.1 # TODO: guess a sensible value
model[:ironoptna].norm.val = 0.0
freeze(model, :ironoptna) # will be freed during last run
Giorgio Calderone
committed
push!(iron_components, :ironoptbr, :ironoptna)
else
println(logio(source), "Ignoring ironopt component on prediction $id (threshold: $threshold)")
Giorgio Calderone
committed
end
end
if length(iron_components) > 0
model[:Iron] = SumReducer(iron_components)
model[:main] = SumReducer([:Continuum, :Iron])
bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
Giorgio Calderone
committed
else
model[:Iron] = @expr m -> [0.]
model[:main] = SumReducer([:Continuum, :Iron])
Giorgio Calderone
committed
end
(:ironuv in keys(model)) && freeze(model, :ironuv)
(:ironoptbr in keys(model)) && freeze(model, :ironoptbr)
(:ironoptna in keys(model)) && freeze(model, :ironoptna)
Giorgio Calderone
committed
end
Giorgio Calderone
committed
# Add emission lines
line_names = [collect(keys(source.line_names[id])) for id in 1:Nspec]
line_groups = [unique(collect(values(source.line_names[id]))) for id in 1:Nspec]
Giorgio Calderone
committed
for id in 1:Nspec
Giorgio Calderone
committed
λ = source.domain[id][:]
resid = source.data[id].val - model() # will be used to guess line normalization
for (cname, comp) in source.line_comps[id]
model[cname] = comp
end
Giorgio Calderone
committed
for (group, lnames) in QSFit.invert_dictionary(source.line_names[id])
Giorgio Calderone
committed
end
model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]...])
Giorgio Calderone
committed
if haskey(model, :MgII_2798)
model[:MgII_2798].voff.low = -1000
model[:MgII_2798].voff.high = 1000
Giorgio Calderone
committed
end
if haskey(model, :OIII_5007_bw)
model[:OIII_5007_bw].fwhm.val = 500
model[:OIII_5007_bw].fwhm.low = 1e2
model[:OIII_5007_bw].fwhm.high = 1e3
model[:OIII_5007_bw].voff.low = 0
model[:OIII_5007_bw].voff.high = 2e3
Giorgio Calderone
committed
end
for cname in line_names[id]
model[cname].norm_integrated = source.options[:norm_integrated]
Giorgio Calderone
committed
end
# Guess values
evaluate!(model)
for cname in line_names[id]
Giorgio Calderone
committed
resid_at_line = Spline1D(λ, resid, k=1, bc="nearest")(c.center.val)
Giorgio Calderone
committed
# If instrumental broadening is not used and the line profile
# is a Gaussian one take spectral resolution into account.
# This is significantly faster than convolving with an
# instrument response but has some limitations:
# - works only with Gaussian profiles;
# - all components must be additive (i.e. no absorptions)
# - further narrow components (besides known emission lines)
# will not be corrected for instrumental resolution
if !source.options[:instr_broadening]
Giorgio Calderone
committed
c.spec_res_kms = source.spectra[id].resolution
else
println(logio(source), "Line $cname is not a Gaussian profile: Can't take spectral resolution into account")
Giorgio Calderone
committed
end
end
end
# Patch parameters
# model[:OIII_4959].norm = model[:OIII_5007].norm / 3
model[:OIII_4959].voff = model[:OIII_5007].voff
@patch! begin
model[:OIII_5007_bw].voff += model[:OIII_5007].voff
model[:OIII_5007_bw].fwhm += model[:OIII_5007].fwhm
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
@patch! begin
# model[:OI_6300].norm = model[:OI_6364].norm / 3
model[:OI_6300].voff = model[:OI_6364].voff
end
@patch! begin
# model[:NII_6549].norm = model[:NII_6583].norm / 3
model[:NII_6549].voff = model[:NII_6583].voff
end
@patch! begin
# model[:SII_6716].norm = model[:SII_6731].norm / 1.5
model[:SII_6716].voff = model[:SII_6731].voff
end
@patch! model[:na_Hb].voff = model[:na_Ha].voff
# The following are required to avoid degeneracy with iron
# template
@patch! begin
model[:Hg].voff = model[:br_Hb].voff
model[:Hg].fwhm = model[:br_Hb].fwhm
end
@patch! begin
model[:br_Hg].voff = model[:br_Hb].voff
model[:br_Hg].fwhm = model[:br_Hb].fwhm
end
@patch! begin
model[:na_Hg].voff = model[:na_Hb].voff
model[:na_Hg].fwhm = model[:na_Hb].fwhm
end
# Ensure luminosity at peak of the broad base component is
# smaller than the associated broad component:
if haskey(model, :br_Hb) &&
haskey(model, :bb_Hb)
model[:bb_Hb].norm.high = 1
model[:bb_Hb].norm.val = 0.5
@patch! model[:bb_Hb].norm *= model[:br_Hb].norm / model[:br_Hb].fwhm * model[:bb_Hb].fwhm
end
if haskey(model, :br_Ha) &&
haskey(model, :bb_Ha)
model[:bb_Ha].norm.high = 1
model[:bb_Ha].norm.val = 0.5
@patch! model[:bb_Ha].norm *= model[:br_Ha].norm / model[:br_Ha].fwhm * model[:bb_Ha].fwhm
# Force Hg and Hb lines to have the same shape as Ha
model[:br_Ha].norm.low = 1.e-10 # avoid division by zero
@patch! begin
model[:na_Hg].norm = model[:br_Hg].norm * (model[:na_Ha].norm / model[:br_Ha].norm)
model[:na_Hg].voff = model[:na_Ha].voff
model[:na_Hg].fwhm = model[:na_Ha].fwhm
Giorgio Calderone
committed
model[:bb_Hg].norm = model[:br_Hg].norm * (model[:bb_Ha].norm / model[:br_Ha].norm)
model[:bb_Hg].voff = model[:bb_Ha].voff
model[:bb_Hg].fwhm = model[:bb_Ha].fwhm
Giorgio Calderone
committed
model[:br_Hg].voff = model[:br_Ha].voff
model[:br_Hg].fwhm = model[:br_Ha].fwhm
Giorgio Calderone
committed
end
@patch! begin
model[:na_Hb].norm = model[:br_Hb].norm * (model[:na_Ha].norm / model[:br_Ha].norm)
model[:na_Hb].voff = model[:na_Ha].voff
model[:na_Hb].fwhm = model[:na_Ha].fwhm
Giorgio Calderone
committed
model[:bb_Hb].norm = model[:br_Hb].norm * (model[:bb_Ha].norm / model[:br_Ha].norm)
model[:bb_Hb].voff = model[:bb_Ha].voff
model[:bb_Hb].fwhm = model[:bb_Ha].fwhm
Giorgio Calderone
committed
model[:br_Hb].voff = model[:br_Ha].voff
model[:br_Hb].fwhm = model[:br_Ha].fwhm
Giorgio Calderone
committed
end
model[:OIII_5007].norm.fixed = 1
model[:OIII_5007].norm.val = 1.
bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
Giorgio Calderone
committed
for lname in line_names[id]
Giorgio Calderone
committed
end
end
# Add unknown lines
Giorgio Calderone
committed
if source.options[:n_unk] > 0
for id in 1:Nspec
tmp = OrderedDict{Symbol, GFit.AbstractComponent}()
Giorgio Calderone
committed
for j in 1:source.options[:n_unk]
tmp[Symbol(:unk, j)] = line_component(TRecipe, QSFit.UnkLine(5e3))
Giorgio Calderone
committed
tmp[Symbol(:unk, j)].norm_integrated = source.options[:norm_integrated]
end
for (cname, comp) in tmp
model[cname] = comp
end
model[:UnkLines] = SumReducer(collect(keys(tmp)))
model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]..., :UnkLines])
Giorgio Calderone
committed
for j in 1:source.options[:n_unk]
Giorgio Calderone
committed
end
end
else
# Here we need a :UnkLines reducer, even when n_unk is 0
for id in 1:Nspec
model = multi[id]
model[:UnkLines] = @expr m -> [0.]
model[:main] = SumReducer([:Continuum, :Iron, line_groups[id]...])
Giorgio Calderone
committed
end
end
Giorgio Calderone
committed
# Set "unknown" line center wavelength where there is a maximum in
# the fit residuals, and re-run a fit.
for id in 1:Nspec
Giorgio Calderone
committed
λ = source.domain[id][:]
λunk = Vector{Float64}()
while true
(length(λunk) >= source.options[:n_unk]) && break
Giorgio Calderone
committed
# Avoid considering again the same region (within 1A) TODO: within resolution
for l in λunk
Δ[findall(abs.(l .- λ) .< 1)] .= 0.
end
# Avoidance regions
for rr in source.options[:unk_avoid]
Δ[findall(rr[1] .< λ .< rr[2])] .= 0.
end
# Do not add lines close to from the edges since these may
# affect qso_cont fitting
Δ[findall((λ .< minimum(λ)*1.02) .|
(λ .> maximum(λ)*0.98))] .= 0.
iadd = argmax(Δ)
(Δ[iadd] <= 0) && break # No residual is greater than 0, skip further residuals....
push!(λunk, λ[iadd])
cname = Symbol(:unk, length(λunk))
model[cname].norm.val = 1.
model[cname].center.val = λ[iadd]
model[cname].center.low = λ[iadd] - λ[iadd]/10. # allow to shift 10%
model[cname].center.high = λ[iadd] + λ[iadd]/10.
thaw(model, cname)
bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
Giorgio Calderone
committed
end
end
Giorgio Calderone
committed
# ----------------------------------------------------------------
# Last run with all parameters free
println(logio(source), "\nLast run with all parameters free...")
Giorgio Calderone
committed
for id in 1:Nspec
model = multi[id]
thaw(model, :qso_cont)
(:galaxy in keys(model)) && thaw(model, :galaxy)
(:balmer in keys(model)) && thaw(model, :balmer)
(:ironuv in keys(model)) && thaw(model, :ironuv)
(:ironoptbr in keys(model)) && thaw(model, :ironoptbr)
(:ironoptna in keys(model)) && thaw(model, :ironoptna)
Giorgio Calderone
committed
for lname in line_names[id]
Giorgio Calderone
committed
end
for j in 1:source.options[:n_unk]
cname = Symbol(:unk, j)
Giorgio Calderone
committed
else
Giorgio Calderone
committed
end
end
end
Giorgio Calderone
committed
# Disable "unknown" lines whose normalization uncertainty is larger
# than 3 times the normalization
needs_fitting = false
for id in 1:Nspec
Giorgio Calderone
committed
for ii in 1:source.options[:n_unk]
cname = Symbol(:unk, ii)
Giorgio Calderone
committed
if bestfit[id][cname].norm.val == 0.
Giorgio Calderone
committed
needs_fitting = true
Giorgio Calderone
committed
elseif bestfit[id][cname].norm.unc / bestfit[id][cname].norm.val > 3
Giorgio Calderone
committed
needs_fitting = true
println(logio(source), "Disabling $cname (unc. / norm. > 3)")
Giorgio Calderone
committed
end
end
end
if needs_fitting
Giorgio Calderone
committed
end
Giorgio Calderone
committed
Giorgio Calderone
committed
elapsed = time() - elapsed
println(logio(source), "\nElapsed time: $elapsed s")
QSFit.close_logio(source)
return out