CL_1ES_1927p654.jl

module CL_1ES_1927p654

using QSFit, DataStructures, GFit, Statistics, Dierckx

import QSFit: DefaultRecipe,
    default_options, known_spectral_lines, LineComponent,
    add_qso_continuum!, add_patch_functs!, guess_emission_lines!, qsfit_multi

abstract type q1927p654 <: DefaultRecipe end

export q1927p654


function default_options(::Type{T}) where T <: q1927p654
    out = default_options(supertype(T))
    out[:use_host_template] = true
    #out[:host_template] = "S0"
    #out[:wavelength_range][1] = 4000
    out[:n_unk] = 0
    out[:wavelength_range] = [3500., 7200.]  # use same range at all epochs
    out[:min_spectral_coverage][:OIII_5007] = 0.5 # ensure line is taken into account
    out[:min_spectral_coverage][:OIII_4363] = 0.5
    return out
end


function QSFit.add_qso_continuum!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: q1927p654
    λ = domain(model)[:]

    comp = QSFit.powerlaw(3000)
    comp.x0.val = median(λ)
    comp.norm.val = Spline1D(λ, pspec.data.val, k=1, bc="error")(comp.x0.val)
    comp.alpha.val  = -1.5
    comp.alpha.low  = -5
    comp.alpha.high =  5

    model[:qso_cont] = comp
    push!(model[:Continuum].list, :qso_cont)
    evaluate!(model)
end

function known_spectral_lines(source::QSO{T}) where T <: q1927p654
    list = [
        # MultiCompLine(                      :Lyb                     , [BroadLine, NarrowLine]),
        # MultiCompLine(                      :Lya                     , [BroadLine, NarrowLine]),
        # NarrowLine(                         :NV_1241                ),
        # BroadLine(    custom_transition(tid=:OI_1306     , 1305.53 )),
        # BroadLine(    custom_transition(tid=:CII_1335    , 1335.31 )),
        # BroadLine(                          :SiIV_1400              ),
        # MultiCompLine(                      :CIV_1549                , [BroadLine, NarrowLine]),
        # BroadLine(                          :HeII_1640              ),
        # BroadLine(    custom_transition(tid=:OIII        , 1665.85 )),
        # BroadLine(    custom_transition(tid=:AlIII       , 1857.4  )),
        # BroadLine(                          :CIII_1909              ),
        # BroadLine(    custom_transition(tid=:CII         , 2326.0  )),
        # BroadLine(    custom_transition(tid=:F2420       , 2420.0  )),
        # MultiCompLine(                      :MgII_2798               , [BroadLine, NarrowLine]),
        # NarrowLine(   custom_transition(tid=:NeVN        , 3346.79 )),
        # NarrowLine(   custom_transition(tid=:NeVI_3426   , 3426.85 )),
        # NarrowLine(   custom_transition(tid=:OII_3727    , 3729.875)),
        # NarrowLine(   custom_transition(tid=:NeIII_3869  , 3869.81 )),
        # BroadLine(                          :Hd                     ),
        MultiCompLine(                      :Hg                      , [BroadLine, NarrowLine, BroadBaseLine]),
        NarrowLine(                         :OIII_4363              ),
        # BroadLine(                          :HeII_4686              ),
        MultiCompLine(                      :Hb                      , [BroadLine, NarrowLine, BroadBaseLine]),
        NarrowLine(                         :OIII_4959              ),
        NarrowLine(                         :OIII_5007              ),
        # AsymmTailLine(                      :OIII_5007               , :blue),
        # BroadLine(                          :HeI_5876               ),
        NarrowLine(                         :OI_6300                ),
        NarrowLine(                         :OI_6364                ),
        NarrowLine(                         :NII_6549               ),
        MultiCompLine(                      :Ha                      , [BroadLine, NarrowLine, BroadBaseLine]),
        NarrowLine(                         :NII_6583               ),
        NarrowLine(                         :SII_6716               ),
        NarrowLine(                         :SII_6731               )]
    return list
end


function QSFit.add_patch_functs!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: q1927p654
    QSFit.add_patch_functs!(parent_recipe(source), pspec, model)

    # Force Hg and Hb lines to have the same shape as Ha
    model[:Ha_br].norm.low = 1.e-10  # avoid division by zero

    @patch! begin
        model[:Hg_na].norm = model[:Hg_br].norm * (model[:Ha_na].norm / model[:Ha_br].norm)
        model[:Hg_na].voff = model[:Ha_na].voff
        model[:Hg_na].fwhm = model[:Ha_na].fwhm

        model[:Hg_bb].norm = model[:Hg_br].norm * (model[:Ha_bb].norm / model[:Ha_br].norm)
        model[:Hg_bb].voff = model[:Ha_bb].voff
        model[:Hg_bb].fwhm = model[:Ha_bb].fwhm

        model[:Hg_br].voff = model[:Ha_br].voff
        model[:Hg_br].fwhm = model[:Ha_br].fwhm
    end

    @patch! begin
        model[:Hb_na].norm = model[:Hb_br].norm * (model[:Ha_na].norm / model[:Ha_br].norm)
        model[:Hb_na].voff = model[:Ha_na].voff
        model[:Hb_na].fwhm = model[:Ha_na].fwhm

        model[:Hb_bb].norm = model[:Hb_br].norm * (model[:Ha_bb].norm / model[:Ha_br].norm)
        model[:Hb_bb].voff = model[:Ha_bb].voff
        model[:Hb_bb].fwhm = model[:Ha_bb].fwhm

        model[:Hb_br].voff = model[:Ha_br].voff
        model[:Hb_br].fwhm = model[:Ha_br].fwhm
    end
end



function guess_emission_lines!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: q1927p654
    QSFit.guess_emission_lines!(parent_recipe(source), pspec, model)
end


function qsfit_multi(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
    elapsed = time()
    @assert length(source.specs) > 1
    @assert 1 <= ref_id <= length(source.specs)
    pspecs = [PreparedSpectrum(source, id=id) for id in 1:length(source.specs)]

    multi = MultiModel()
    println(logio(source), "\nFit continuum components...")
    for id in 1:length(pspecs)
        pspec = pspecs[id]
        model = Model(pspec.domain)
        model[:Continuum] = SumReducer([])
        model[:main] = SumReducer([])
        push!(model[:main].list, :Continuum)
        select_reducer!(model, :main)
        delete!(model.revals, :default_sum)

        # TODO if source.options[:instr_broadening]
        # TODO     GFit.set_instr_response!(model[1], (l, f) -> instrumental_broadening(l, f, source.spectra[id].resolution))
        # TODO end

        QSFit.add_qso_continuum!(source, pspec, model)
        QSFit.add_host_galaxy!(source, pspec, model)
        QSFit.add_balmer_cont!(source, pspec, model)

        push!(multi, model)
        if id != ref_id
            @patch! multi[id][:galaxy].norm = multi[ref_id][:galaxy].norm
        end
    end
    fitres = fit!(source, multi, pspecs)

    for id in 1:length(pspecs)
        model = multi[id]
        pspec = pspecs[id]
        QSFit.renorm_cont!(source, pspec, model)
        freeze(model, :qso_cont)
        haskey(model, :galaxy)  &&  freeze(model, :galaxy)
        haskey(model, :balmer)  &&  freeze(model, :balmer)
        evaluate!(model)
    end
    evaluate!(multi)

    println(logio(source), "\nFit iron templates...")
    for id in 1:length(pspecs)
        model = multi[id]
        pspec = pspecs[id]
        model[:Iron] = SumReducer([])
        push!(model[:main].list, :Iron)
        QSFit.add_iron_uv!( source, pspec, model)
        QSFit.add_iron_opt!(source, pspec, model)

        if length(model[:Iron].list) > 0
            fitres = fit!(source, model, pspec)
            haskey(model, :ironuv   )  &&  freeze(model, :ironuv)
            haskey(model, :ironoptbr)  &&  freeze(model, :ironoptbr)
            haskey(model, :ironoptna)  &&  freeze(model, :ironoptna)
        end
        evaluate!(model)
    end
    evaluate!(multi)

    println(logio(source), "\nFit known emission lines...")
    for id in 1:length(pspecs)
        model = multi[id]
        pspec = pspecs[id]
        QSFit.add_emission_lines!(source, pspec, model)
        QSFit.guess_emission_lines!(source, pspec, model)
        QSFit.add_patch_functs!(source, pspec, model)

        multi[id][:OIII_5007].norm.val = 1.
        # multi[id][:OIII_5007].norm.fixed = true
    end
    fitres = fit!(source, multi, pspecs)
    for id in 1:length(pspecs)
        model = multi[id]
        pspec = pspecs[id]
        for lname in keys(pspec.lcs)
            freeze(model, lname)
        end
    end
    evaluate!(multi)

    println(logio(source), "\nFit unknown emission lines...")
    for id in 1:length(pspecs)
        model = multi[id]
        pspec = pspecs[id]
        QSFit.add_unknown_lines!(source, pspec, model)
    end
    evaluate!(multi)

    println(logio(source), "\nLast run with all parameters free...")
    for id in 1:length(pspecs)
        model = multi[id]
        pspec = pspecs[id]
        thaw(model, :qso_cont)
        haskey(model, :galaxy   )  &&  thaw(model, :galaxy)
        haskey(model, :balmer   )  &&  thaw(model, :balmer)
        haskey(model, :ironuv   )  &&  thaw(model, :ironuv)
        haskey(model, :ironoptbr)  &&  thaw(model, :ironoptbr)
        haskey(model, :ironoptna)  &&  thaw(model, :ironoptna)
        for lname in keys(pspec.lcs)
            thaw(model, lname)
        end
        for j in 1:source.options[:n_unk]
            cname = Symbol(:unk, j)
            if model[cname].norm.val > 0
                thaw(model, cname)
            else
                freeze(model, cname)
            end
        end
    end
    fitres = fit!(source, multi, pspecs)

    rerun = false
    for id in 1:length(pspecs)
        model = multi[id]
        pspec = pspecs[id]
        rerun = rerun || QSFit.neglect_weak_features!(source, pspec, model, fitres)
    end
    if rerun
        println(logio(source), "\nRe-run fit...")
        fitres = fit!(source, multi, pspecs)
    end
    println(logio(source))
    show(logio(source), fitres)

    out = QSFit.QSFitMultiResults(source, pspecs, multi, fitres)
    elapsed = time() - elapsed
    println(logio(source), "\nElapsed time: $elapsed s")
    close_logio(source)
    return out
end


end