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

Updated

parent aefc458a

run.jl

+43 −61
Original line number Diff line number Diff line
@@ -8,8 +8,9 @@ using Revise, Statistics, Serialization, Dierckx, DataFrames 
using QSFit, GFit, Gnuplot, GFitViewer, MyAstroUtils

using CL_1ES_1927p654

using Dates

using FileIO





mkdir("output")

epoch_filenames = Vector{String}()

for (root, dirs, files) in walkdir("AT2018zf")

    for file in files
@@ -19,8 +20,6 @@ for (root, dirs, files) in walkdir("AT2018zf") 
    end

end







function read_spec(epoch_id; kw...)

    if epoch_id == "B03"

        file = "boller2003.txt"
@@ -70,14 +69,14 @@ if !isfile("output/scale_fwhm.dat") 
            # source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"

            spec = read_spec(id)

            add_spec!(source, spec);

            (model, bestfit) = fit(source);

            viewer(model, source, bestfit, showcomps=[:qso_cont, :galaxy, :balmer],

            res = fit(source);

            viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],

                   filename="output/results_$(id).html")



            # Update all_scale and store FWHM

            @info bestfit[:OIII_5007].norm.val

            all_scale[id] *= bestfit[:OIII_5007].norm.val

            all_fwhm[ id]  = bestfit[:OIII_5007].fwhm.val

            @info res.bestfit[:OIII_5007].norm.val

            all_scale[id] *= res.bestfit[:OIII_5007].norm.val

            all_fwhm[ id]  = res.bestfit[:OIII_5007].fwhm.val

        catch

        end

    end
@@ -134,11 +133,12 @@ resolution = Dict( 
    :lowres  => 350.,

    :highres => 150.)





job = :all

chosen_epochs = dict_chosen_epochs[job]

Nloop = 6

(all_scale, all_fwhm) = deserialize("output/scale_fwhm.dat")

for loop in 1:Nloop

    if !isfile("output/results_$(job)_$(loop).dat")

        source = QSO{q1927p654}("1ES 1927+654", Z, ebv=EBV);

        source.options[:min_spectral_coverage][:OIII_5007] = 0.5

        @gp :zoom "set grid" :-
@@ -153,24 +153,9 @@ for loop in 1:Nloop 
        viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],

               filename="output/results_$(job)_$(loop)_rebin4.html", rebin=4)



    models = Dict()

        # Find best normalization for [OIII]

        model2 = deepcopy(res.model);

        for id in 1:length(chosen_epochs)

        models[(id, :x)] = domain(model[id])[:]

        models[(id, :y)] = model[id]()

        models[(id, :l5100)] = 5100 .* Spline1D(domain(model[id])[:], model[id]())(5100.)



        for cname in [:qso_cont, :OIII_5007,

                      :br_Ha, :na_Ha, :bb_Ha,

                      :br_Hb, :na_Hb, :bb_Hb]

            models[(id, cname)] = model[id](cname)

        end

    end

    serialize("output/results_$(job)_$(loop).dat", res)

    FileIO.save("output/results_$(job)_$(loop).jld2", "res", res)



    model2 = deepcopy(res.model)

    for id in 1:length(chosen_epochs)

        @info id

            for cname in collect(keys(model2[id]))

                if cname == :OIII_5007

                    model2[id][cname].norm.fixed = false
@@ -181,20 +166,17 @@ for loop in 1:Nloop 
        end

        mzer = GFit.cmpfit()

        mzer.config.ftol = mzer.config.gtol = mzer.config.xtol = 1.e-6

    bestfit2 = fit!(model2, source.data, minimizer=mzer)

    #=

    @gp([bestfit[ id][:OIII_5007].norm.val for id in 1:length(chosen_epochs)],

        [bestfit2[id][:OIII_5007].norm.val for id in 1:length(chosen_epochs)])

    @gp([bestfit2[id][:OIII_5007].norm.val for id in 1:length(chosen_epochs)],

        ([bestfit[ id][:OIII_5007].fwhm.val for id in 1:length(chosen_epochs)] .-

         [bestfit2[id][:OIII_5007].fwhm.val for id in 1:length(chosen_epochs)]) ./

        [bestfit[ id][:OIII_5007].fwhm.val for id in 1:length(chosen_epochs)])

    =#

        bestfit2 = fit!(model2, source.data, minimizer=mzer);

        

        OIII_norm = fill(0., length(chosen_epochs))

        for id in 1:length(chosen_epochs)

            OIII_norm[id] = bestfit2[id][:OIII_5007].norm.val

            all_scale[chosen_epochs[id]] *= bestfit2[id][:OIII_5007].norm.val

        end

    serialize("output/scale_$(job)_$(loop).dat", all_scale)

        serialize("output/results_$(job)_$(loop).dat", (res, all_scale, OIII_norm))

    else

        (res, all_scale, OIII_norm) = deserialize("output/results_$(job)_$(loop).dat")

    end

end