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

Updated

parent a9a596c3

run.jl

+138 −115
Original line number Diff line number Diff line 
using Pkg

Pkg.activate(".")



Z = 0.019422

EBV = 0.077



using Revise, Statistics, Dierckx, DataFrames

using QSFit, GFit, Gnuplot, GFitViewer, MyAstroUtils

using CL_1ES_1927p654
@@ -11,34 +8,42 @@ using Dates 
using TextParse

import StatsBase, JLD2



path = "output-" * readlines(`git branch --show-current`)[1]

mkdir("$(path)")

epoch_filenames = Vector{String}()

function read_epochs()

    out = DataFrame(id=Int[], filename=String[], date=String[], scale=Float64[], OIII_fwhm=Float64[])

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

        for file in files

            if file[end-3:end] == ".dat"

            push!(epoch_filenames, root * "/" * file)

                filename = root * "/" * file

                date = filename[27:end-4]

                push!(out, (nrow(out)+1, filename, date, NaN, NaN))

            end

        end

    end

    @assert issorted(out.date)

    return out

end



function read_spec(epoch_id; kw...)

    if epoch_id == "B03"

function read_spec(id; kw...)

    if id == "B03"

        file = "boller2003.txt"

    else

        file = epoch_filenames[epoch_id]

        irow = findfirst(epochs.id .== id)

        file = epochs[irow, :filename]

    end

    println(file)

	spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$epoch_id: $file", kw...);

    if isa(epoch_id, Number)

        spec.flux ./= all_scale[epoch_id]

        spec.err  ./= all_scale[epoch_id]

	spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$id: $file", kw...);

    if isa(id, Number)

        spec.flux ./= epochs[irow, :scale]

        spec.err  ./= epochs[irow, :scale]

    end

    return spec

end





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

function analyze_boller03()

    file = opt.path * "/boller.dat"

    if !isfile(file)

        source = QSO{q1927p654}("1ES 1927+654 (Boller03)", opt.z, ebv=opt.ebv);

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

        # source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"

        #=
@@ -50,47 +55,110 @@ spec = read_spec("B03", resolution=200.) 
        add_spec!(source, spec);

        res = qsfit(source);

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

       filename="$(path)/results_boller.html")

println(res.model[:OIII_5007].norm.val )

# 0.0002904273540659144

println(res.model[:OIII_5007].norm.val / res.model[:galaxy].norm.val )

# 60.25349678748317

               filename="$(opt.path)/results_boller.html")

        JLD2.save_object(file, res)

    else

        res = JLD2.load_object(file)

    end

    @assert res.model[:OIII_5007].norm.val == 0.00029042977121756283

    @assert res.model[:OIII_5007].norm.val / res.model[:galaxy].norm.val == 60.07373299362222

    return res

end





if !isfile("$(path)/scale_fwhm.dat")

    all_scale = fill(1., length(epoch_filenames))

    all_fwhm  = fill(NaN,length(epoch_filenames))

    for id in 1:length(epoch_filenames)

function estimate_SE_scale!()

    file = opt.path * "/" * "SE_scale.dat"

    if !isfile(file)

        for irow in 1:nrow(epochs)

            @info irow

            id = epochs[irow, :id]

            epochs[irow, :scale] = 1.

            spec = read_spec(id)

        all_scale[id] *= median(spec.flux) / 25000

            epochs[irow, :scale] = median(spec.flux) / 25000

            

        try

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

            source = QSO{q1927p654}("1ES 1927+654 ($id)", opt.z, ebv=opt.ebv);

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

            # 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);

            res = qsfit(source);

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

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

                   filename="$(opt.path)/results_$(id).html")

                

            # Update all_scale and store FWHM

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

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

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

        catch

            epochs[irow, :scale] *= res.model[:OIII_5007].norm.val

            epochs[irow, :OIII_fwhm] = res.model[:OIII_5007].fwhm.val

        end

        @assert all(epochs[.!isfinite.(epochs.OIII_fwhm), :scale] .== NaN)

        JLD2.save_object(file, epochs)

    else

        empty!(epochs)

        append!(epochs, JLD2.load_object(file))

    end

end





function analyze_job(job; Nloop = 6)

    out = fill(NaN, length(epoch_ids[job]))



    estimate_SE_scale!()

    for loop in 1:Nloop

        file = "$(opt.path)/results_$(job)_$(loop).dat"

        if !isfile(file)

            source = QSO{q1927p654}("1ES 1927+654", opt.z, ebv=opt.ebv);

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

            source.options[:wavelength_range] = [3500, 6900]

            @gp :zoom "set grid" :-

            for id in epoch_ids[job]

                spec = read_spec(id, resolution=get(resolution, job, NaN))

                add_spec!(source, spec);

                @gp :- :zoom xr=[4750,5150] spec.λ ./ (1 + source.z) spec.flux "w l t '$(id)'"

            end

    all_scale[.!isfinite.(all_fwhm)] .= NaN

    JLD2.save_object("$(path)/scale_fwhm.dat", (all_scale, all_fwhm))

            res = qsfit_multi(source);

            JLD2.save_object(file, res)

            viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], filename="$(opt.path)/results_$(job)_$(loop).html")

            viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], filename="$(opt.path)/results_$(job)_$(loop)_rebin4.html", rebin=4)

        else

    (all_scale, all_fwhm) = JLD2.load_object("$(path)/scale_fwhm.dat")

            res = JLD2.load_object(file);

        end



        OIII_norm = Float64[]

        for i in 1:length(epoch_ids[job])

            id = epoch_ids[job][i]

            push!(OIII_norm, res.multi[i][:OIII_5007].norm.val)

            epochs[epochs.id .== id, :scale] *= res.multi[i][:OIII_5007].norm.val

        end

        out = hcat(out, OIII_norm)        

    end

    out = out[:,2:end]



    @gp "set grid" :-

    for i in 1:length(epoch_ids[job])

        id = epoch_ids[job][i]

        @gp :- 1:Nloop out[i, :] "w lp t '$(id)'"

    end



    return out

end









const opt = (

z = 0.01942,

ebv = 0.077,

path = "output")



try;  mkdir(opt.path);  catch; end

const epochs = read_epochs()

analyze_boller03()

estimate_SE_scale!()



@gp :prenorm "set key bottom right" "set grid" xlab="Epoch" ylab="Value" :-

@gp :- :prenorm all_scale./1e-18 "w lp t '[OIII] norm.'" :-

@gp :- :prenorm  all_fwhm./1000 "w lp t '[OIII] FWHM (x1000 km/s)"

save(:prenorm, term="png size 800,600", output="$(path)/oiii_norm_fwhm.png")

@gp :- :prenorm epochs.scale./1e-18 "w lp t '[OIII] norm.'" :-

@gp :- :prenorm epochs.OIII_fwhm./1000 "w lp t '[OIII] FWHM (x1000 km/s)"

save(:prenorm, term="png size 800,600", output="$(opt.path)/oiii_norm_fwhm.png")



w = 10 .^(2:0.01:log10(400));

r = sqrt.(400^2 .- w.^2) ./ 2.355;
@@ -105,28 +173,29 @@ r = sqrt.(900^2 .- w.^2) ./ 2.355; 
# Ensure no sample is negative

@gp :allepochs "set grid" xr=[3e3,1e4] cbr=[1,29] :-

@gp :- :allepochs cblabel="Epoch" xlab="Rest frame wavelength[A]" ylab="Flux density [arb.units]"

for id in 1:length(epoch_filenames)

for id in epochs.id

	s = read_spec(id)

	x = s.λ;

	y = s.flux;

    y .-= median(y)

    y ./= maximum(y)

    @gp :- :allepochs x./(1+Z) y "u 1:(\$2+$id):($id) w l notit lc pal" :-

    @gp :- :allepochs x ./ (1 + opt.z) y "u 1:(\$2+$id):($id) w l notit lc pal" :-

end

@gp :- :allepochs yr=[1, 29]

save(:allepochs, term="png size 800,600", output="$(path)/all_epochs.png")

save(:allepochs, term="png size 800,600", output="$(opt.path)/all_epochs.png")

@gp :- :allepochs xr=[4900, 5100]

save(:allepochs, term="png size 800,2000", output="$(path)/all_epochs_zoom.png")

save(:allepochs, term="png size 800,2000", output="$(opt.path)/all_epochs_zoom.png")







dict_chosen_epochs = Dict(

epoch_ids = Dict(

    :vecchi   => [17, 10, 21, 23],

    :lowres   => [8,9,10,11,12,14,15,16,18,21,22,23,24,25],

    :highres  => [5, 7, 13, 17, 20],

    :test     => [7, 13],

    :limited  => [1,2,3,13],

    :all      => collect(5:26))

# deleteat!(dict_chosen_epochs[:all], 4)  #insufficient coverage in na_Hg

# deleteat!(epoch_ids[:all], 4)  #insufficient coverage in na_Hg



# Can't use resolution smaller than 150 km / s otherwise some line

# will be neglected because of spectral coverage
@@ -134,55 +203,9 @@ resolution = Dict() 
#    :lowres  => 350.,

#    :highres => 150.)



job = :all

chosen_epochs = dict_chosen_epochs[job]

Nloop = 6

(all_scale, all_fwhm) = JLD2.load_object("$(path)/scale_fwhm.dat")

for loop in 1:Nloop

    file = "$(path)/results_$(job)_$(loop).dat"

    if !isfile(file)

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

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

        source.options[:wavelength_range] = [3500, 6900]

        @gp :zoom "set grid" :-

        for id in chosen_epochs

            spec = read_spec(id, resolution=get(resolution, job, NaN))

            add_spec!(source, spec);

            @gp :- :zoom xr=[4750,5150] spec.λ ./ (1 + source.z) spec.flux "w l t '$(id)'"

        end

        res = qsfit_multi(source);

        JLD2.jldsave(file; res)

        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)



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

        for id in 1:length(chosen_epochs)

            OIII_norm[id] = res.multi[id][:OIII_5007].norm.val

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

        end

        JLD2.jldsave(file; res, all_scale, OIII_norm)

    else

        # 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));

    end

end





OIII_norm_evol = fill(NaN, length(chosen_epochs))

for loop in 1:Nloop

    file = "$(path)/results_$(job)_$(loop).dat"

    res, all_scale, OIII_norm = values(JLD2.load(file))

    OIII_norm_evol = hcat(OIII_norm_evol, OIII_norm)

    println(res.fitres.fitstat / res.fitres.dof)

end

OIII_norm_evol = OIII_norm_evol[:,2:end]

@gp "set grid" :-

for id in 1:length(chosen_epochs)

    @gp :- 1:Nloop OIII_norm_evol[id, :] "w lp t '$(chosen_epochs[id])'"

end



analyze_job(:test)

analyze_job(:highres)

ddd
@@ -256,7 +279,7 @@ for i in 1:nrow(tab) 
    tab[i, :instr] = tmp[j, :Inst]

end



f = FITS("$(path)/1ES_1927p654_results_$(job).fits", "w")

f = FITS("$(opt.path)/1ES_1927p654_results_$(job).fits", "w")

write(f, tab)

close(f)
@@ -298,7 +321,7 @@ for i in 1:nrow(tab) 
    ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"

    @gp :- domain(res.model[i])[:] res.model[i](:qso_cont) ss

end

save(:prenorm, term="png size 800,600", output="$(path)/evolution.png")

save(:prenorm, term="png size 800,600", output="$(opt.path)/evolution.png")





@gp "set grid"
@@ -310,4 +333,4 @@ for i in 1:length(chosen_epochs) 
    ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"

    @gp :- x y  .- Spline1D(x, y0)(5007.) ss

end

save(:prenorm, term="png size 800,600", output="$(path)/evolution_oiii_norm.png")
 
save(:prenorm, term="png size 800,600", output="$(opt.path)/evolution_oiii_norm.png")