Commit b36a894b authored by Giorgio Calderone's avatar Giorgio Calderone
Browse files

Updated

parent 8a7ba810
Loading
Loading
Loading
Loading
+10 −13
Original line number Diff line number Diff line
@@ -52,9 +52,9 @@ res = fit(source);
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
       filename="$(path)/results_boller.html")
println(res.bestfit[:OIII_5007].norm.val )
# 0.00029150559312419105
# 0.0002913168348334021
println(res.bestfit[:OIII_5007].norm.val / res.bestfit[:galaxy].norm.val )
# 63.87625309927029
# 64.65011509608263



@@ -131,9 +131,9 @@ dict_chosen_epochs = Dict(

# Can't use resolution smaller than 150 km / s otherwise some line
# will be neglected because of spectral coverage
resolution = Dict(
    :lowres  => 350.,
    :highres => 150.)
resolution = Dict()
#    :lowres  => 350.,
#    :highres => 150.)

job = :all
chosen_epochs = dict_chosen_epochs[job]
@@ -151,10 +151,8 @@ for loop in 1:Nloop
            @gp :- :zoom xr=[4750,5150] spec.λ ./ (1 + source.z) spec.flux "w l t '$(id)'"
        end
        res = multi_fit(source);
        viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
               filename="$(path)/results_$(job)_$(loop).html")
        viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
               filename="$(path)/results_$(job)_$(loop)_rebin4.html", rebin=4)
        viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], filename="$(path)/results_$(job)_$(loop).html")
        viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], filename="$(path)/results_$(job)_$(loop)_rebin4.html", rebin=4)

        # Find best normalization for [OIII]
        multi2 = deepcopy(res.multi);
@@ -162,14 +160,13 @@ for loop in 1:Nloop
            for cname in collect(keys(multi2[id]))
                if cname == :OIII_5007
                    multi2[id][cname].norm.fixed = false
                    multi2[id][cname].fwhm.fixed = true  # degenerate with norm
                else
                    freeze(multi2[id], cname)
                end
            end
        end
        mzer = GFit.cmpfit()
        mzer.config.ftol = mzer.config.gtol = mzer.config.xtol = 1.e-6
        bestfit2 = fit!(multi2, source.data, minimizer=mzer);
        bestfit2 = fit!(multi2, source.data, minimizer=GFit.cmpfit())
        
        OIII_norm = fill(0., length(chosen_epochs))
        for id in 1:length(chosen_epochs)
@@ -181,7 +178,7 @@ for loop in 1:Nloop
        # res       = JLD2.load(file, "res")
        # all_scale = JLD2.load(file, "all_scale")
        # OIII_norm = JLD2.load(file, "OIII_norm")
        res, all_scale, OIII_norm = values(JLD2.load(file))
        res, all_scale, OIII_norm = values(JLD2.load(file));
    end
end

+35 −68
Original line number Diff line number Diff line
@@ -37,14 +37,11 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654

            # Split total flux between continuum and host galaxy
            vv = Spline1D(λ, source.data[id].val, k=1, bc="error")(5500.)
            model[:galaxy].norm.val  = 1/3 * vv
            model[:qso_cont].x0.val *= 2/3 * vv / Spline1D(λ, model(:qso_cont), 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.)

            if id != ref_id
                model[:galaxy].norm.fixed = true
                @patch! multi m -> begin
                    m[id][:galaxy].norm = m[ref_id][:galaxy].norm
                end
                @patch! multi[id][:galaxy].norm = multi[ref_id][:galaxy].norm
            end
        end

@@ -62,7 +59,7 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
            c.ratio.fixed = false
            c.ratio.low  = 0.1
            c.ratio.high = 1
            @patch! model m -> m[:balmer].norm *= m[:qso_cont].norm
            @patch! model[:balmer].norm *= model[:qso_cont].norm
        end
    end
    bestfit = fit!(multi, source.data, minimizer=mzer);  show(logio(source), bestfit)
@@ -204,28 +201,14 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
        end

        # Patch parameters
        if  haskey(model, :OIII_4959)  &&
            haskey(model, :OIII_5007)
            model[:OIII_4959].voff.fixed = true
            @patch! model m -> begin
                m[:OIII_4959].voff = m[:OIII_5007].voff
            end
        end
        if  haskey(model, :NII_6549)  &&
            haskey(model, :NII_6583)
            # model[:NII_6549].norm.fixed = true
            model[:NII_6549].voff.fixed = true
            @patch! model m -> begin
                # m[:NII_6549].norm = m[:NII_6583].norm / 3
                m[:NII_6549].voff = m[:NII_6583].voff
            end
        end
        if  haskey(model, :OIII_5007_bw)  &&
            haskey(model, :OIII_5007)
            @patch! model m -> begin
                m[:OIII_5007_bw].voff += m[:OIII_5007].voff
                m[:OIII_5007_bw].fwhm += m[:OIII_5007].fwhm
        @patch! model[:OIII_4959].voff = model[:OIII_5007].voff
        @patch! begin
            # model[:NII_6549].norm = model[:NII_6583].norm / 3
            model[:NII_6549].voff = model[:NII_6583].voff
        end
        @patch! begin
            model[:OIII_5007_bw].voff += model[:OIII_5007].voff
            model[:OIII_5007_bw].fwhm += model[:OIII_5007].fwhm
        end

        # Ensure luminosity at peak of the broad base component is
@@ -234,58 +217,42 @@ function multi_fit(source::QSO{TRecipe}; ref_id=1) where TRecipe <: q1927p654
            haskey(model, :bb_Hb)
            model[:bb_Hb].norm.high = 1
            model[:bb_Hb].norm.val  = 0.5
            @patch! model m -> m[:bb_Hb].norm *= m[:br_Hb].norm / m[:br_Hb].fwhm * m[:bb_Hb].fwhm
            @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 m -> m[:bb_Ha].norm *= m[:br_Ha].norm / m[:br_Ha].fwhm * m[:bb_Ha].fwhm
            @patch! model[:bb_Ha].norm *= model[:br_Ha].norm / model[:br_Ha].fwhm * model[:bb_Ha].fwhm
        end

        # Avoid division by zero
        model[:br_Ha].norm.low = 1.e-10

        model[:na_Hg].norm.fixed = 1
        model[:na_Hg].fwhm.fixed = 1
        model[:na_Hg].voff.fixed = 1
        model[:bb_Hg].norm.fixed = 1
        model[:bb_Hg].fwhm.fixed = 1
        model[:bb_Hg].voff.fixed = 1
        model[:br_Hg].fwhm.fixed = 1
        model[:br_Hg].voff.fixed = 1
        @patch! model m -> begin
            m[:na_Hg].norm = m[:br_Hg].norm * (m[:na_Ha].norm / m[:br_Ha].norm)
            m[:na_Hg].voff = m[:na_Ha].voff
            m[:na_Hg].fwhm = m[:na_Ha].fwhm

            m[:bb_Hg].norm = m[:br_Hg].norm * (m[:bb_Ha].norm / m[:br_Ha].norm)
            m[:bb_Hg].voff = m[:bb_Ha].voff
            m[:bb_Hg].fwhm = m[:bb_Ha].fwhm

            m[:br_Hg].voff = m[:br_Ha].voff
            m[:br_Hg].fwhm = m[:br_Ha].fwhm
        end

        model[:na_Hb].norm.fixed   = 1
        model[:na_Hb].fwhm.fixed   = 1
        model[:na_Hb].voff.fixed   = 1
        model[:bb_Hb].norm.fixed = 1
        model[:bb_Hb].fwhm.fixed = 1
        model[:bb_Hb].voff.fixed = 1
        model[:br_Hb].fwhm.fixed   = 1
        model[:br_Hb].voff.fixed   = 1
        @patch! model m -> begin
            m[:na_Hb].norm = m[:br_Hb].norm * (m[:na_Ha].norm / m[:br_Ha].norm)
            m[:na_Hb].voff = m[:na_Ha].voff
            m[:na_Hb].fwhm = m[:na_Ha].fwhm

            m[:bb_Hb].norm = m[:br_Hb].norm * (m[:bb_Ha].norm / m[:br_Ha].norm)
            m[:bb_Hb].voff = m[:bb_Ha].voff
            m[:bb_Hb].fwhm = m[:bb_Ha].fwhm

            m[:br_Hb].voff = m[:br_Ha].voff
            m[:br_Hb].fwhm = m[:br_Ha].fwhm
        @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

            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

            model[:br_Hg].voff = model[:br_Ha].voff
            model[:br_Hg].fwhm = model[:br_Ha].fwhm
        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

            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

            model[:br_Hb].voff = model[:br_Ha].voff
            model[:br_Hb].fwhm = model[:br_Ha].fwhm
        end

        model[:OIII_5007].norm.fixed = 1
+36 −66
Original line number Diff line number Diff line
@@ -26,8 +26,8 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654

        # Split total flux between continuum and host galaxy
        vv = Spline1D(λ, source.data[id].val, k=1, bc="error")(5500.)
        model[:galaxy].norm.val  = 1/3 * vv
        model[:qso_cont].x0.val *= 2/3 * vv / Spline1D(λ, model(:qso_cont), 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.)
    end

    # Balmer continuum and pseudo-continuum
@@ -44,7 +44,7 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
        c.ratio.fixed = false
        c.ratio.low  = 0.1
        c.ratio.high = 1
        @patch! model m -> m[:balmer].norm *= m[:qso_cont].norm
        @patch! model[:balmer].norm *= model[:qso_cont].norm
    end

    bestfit = fit!(model, source.data[id], minimizer=mzer);  show(logio(source), bestfit)
@@ -175,28 +175,14 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
    end

    # Patch parameters
    if  haskey(model, :OIII_4959)  &&
        haskey(model, :OIII_5007)
        model[:OIII_4959].voff.fixed = true
        @patch! model m -> begin
            m[:OIII_4959].voff = m[:OIII_5007].voff
        end
    end
    if  haskey(model, :NII_6549)  &&
        haskey(model, :NII_6583)
        # model[:NII_6549].norm.fixed = true
        model[:NII_6549].voff.fixed = true
        @patch! model m -> begin
            # m[:NII_6549].norm = m[:NII_6583].norm / 3
            m[:NII_6549].voff = m[:NII_6583].voff
        end
    end
    if  haskey(model, :OIII_5007_bw)  &&
        haskey(model, :OIII_5007)
        @patch! model m -> begin
            m[:OIII_5007_bw].voff += m[:OIII_5007].voff
            m[:OIII_5007_bw].fwhm += m[:OIII_5007].fwhm
    @patch! model[:OIII_4959].voff = model[:OIII_5007].voff
    @patch! begin
        # model[:NII_6549].norm = model[:NII_6583].norm / 3
        model[:NII_6549].voff = model[:NII_6583].voff
    end
    @patch! begin
        model[:OIII_5007_bw].voff += model[:OIII_5007].voff
        model[:OIII_5007_bw].fwhm += model[:OIII_5007].fwhm
    end

    # Ensure luminosity at peak of the broad base component is
@@ -205,58 +191,42 @@ function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
        haskey(model, :bb_Hb)
        model[:bb_Hb].norm.high = 1
        model[:bb_Hb].norm.val  = 0.5
        @patch! model m -> m[:bb_Hb].norm *= m[:br_Hb].norm / m[:br_Hb].fwhm * m[:bb_Hb].fwhm
        @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 m -> m[:bb_Ha].norm *= m[:br_Ha].norm / m[:br_Ha].fwhm * m[:bb_Ha].fwhm
        @patch! model[:bb_Ha].norm *= model[:br_Ha].norm / model[:br_Ha].fwhm * model[:bb_Ha].fwhm
    end

    # Avoid division by zero
    model[:br_Ha].norm.low = 1.e-10

    model[:na_Hg].norm.fixed = 1
    model[:na_Hg].fwhm.fixed = 1
    model[:na_Hg].voff.fixed = 1
    model[:bb_Hg].norm.fixed = 1
    model[:bb_Hg].fwhm.fixed = 1
    model[:bb_Hg].voff.fixed = 1
    model[:br_Hg].fwhm.fixed = 1
    model[:br_Hg].voff.fixed = 1
    @patch! model m -> begin
        m[:na_Hg].norm = m[:br_Hg].norm * (m[:na_Ha].norm / m[:br_Ha].norm)
        m[:na_Hg].voff = m[:na_Ha].voff
        m[:na_Hg].fwhm = m[:na_Ha].fwhm

        m[:bb_Hg].norm = m[:br_Hg].norm * (m[:bb_Ha].norm / m[:br_Ha].norm)
        m[:bb_Hg].voff = m[:bb_Ha].voff
        m[:bb_Hg].fwhm = m[:bb_Ha].fwhm

        m[:br_Hg].voff = m[:br_Ha].voff
        m[:br_Hg].fwhm = m[:br_Ha].fwhm
    end

    model[:na_Hb].norm.fixed = 1
    model[:na_Hb].fwhm.fixed = 1
    model[:na_Hb].voff.fixed = 1
    model[:bb_Hb].norm.fixed = 1
    model[:bb_Hb].fwhm.fixed = 1
    model[:bb_Hb].voff.fixed = 1
    model[:br_Hb].fwhm.fixed = 1
    model[:br_Hb].voff.fixed = 1
    @patch! model m -> begin
        m[:na_Hb].norm = m[:br_Hb].norm * (m[:na_Ha].norm / m[:br_Ha].norm)
        m[:na_Hb].voff = m[:na_Ha].voff
        m[:na_Hb].fwhm = m[:na_Ha].fwhm

        m[:bb_Hb].norm = m[:br_Hb].norm * (m[:bb_Ha].norm / m[:br_Ha].norm)
        m[:bb_Hb].voff = m[:bb_Ha].voff
        m[:bb_Hb].fwhm = m[:bb_Ha].fwhm

        m[:br_Hb].voff = m[:br_Ha].voff
        m[:br_Hb].fwhm = m[:br_Ha].fwhm
    @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

        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

        model[:br_Hg].voff = model[:br_Ha].voff
        model[:br_Hg].fwhm = model[:br_Ha].fwhm
    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

        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

        model[:br_Hb].voff = model[:br_Ha].voff
        model[:br_Hb].fwhm = model[:br_Ha].fwhm
    end

    bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)