Skip to content
Commits on Source (2)
Hide whitespace changes
Inline Side-by-side
run.jl
View file @ 4d11b546
using Pkg using Pkg
Pkg.activate(".") Pkg.activate(".")
Z = 0.019422
EBV = 0.077
using Revise, Statistics, Dierckx, DataFrames using Revise, Statistics, Dierckx, DataFrames
using QSFit, GFit, Gnuplot, GFitViewer, MyAstroUtils using QSFit, GFit, Gnuplot, GFitViewer, MyAstroUtils
using CL_1ES_1927p654 using CL_1ES_1927p654
...@@ -11,86 +8,157 @@ using Dates ...@@ -11,86 +8,157 @@ using Dates
using TextParse using TextParse
import StatsBase, JLD2 import StatsBase, JLD2
path = "output-" * readlines(`git branch --show-current`)[1] function read_epochs()
mkdir("$(path)") out = DataFrame(id=Int[], filename=String[], date=String[], scale=Float64[], OIII_fwhm=Float64[])
epoch_filenames = Vector{String}() for (root, dirs, files) in walkdir("AT2018zf")
for (root, dirs, files) in walkdir("AT2018zf") for file in files
for file in files if file[end-3:end] == ".dat"
if file[end-3:end] == ".dat" filename = root * "/" * file
push!(epoch_filenames, root * "/" * file) date = filename[27:end-4]
push!(out, (nrow(out)+1, filename, date, NaN, NaN))
end
end end
end end
@assert issorted(out.date)
return out
end end
function read_spec(epoch_id; kw...) function read_spec(id; kw...)
if epoch_id == "B03" if id == "B03"
file = "boller2003.txt" file = "boller2003.txt"
else else
file = epoch_filenames[epoch_id] irow = findfirst(epochs.id .== id)
file = epochs[irow, :filename]
end end
println(file) println(file)
spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$epoch_id: $file", kw...); spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$id: $file", kw...);
if isa(epoch_id, Number) if isa(id, Number)
spec.flux ./= all_scale[epoch_id] spec.flux ./= epochs[irow, :scale]
spec.err ./= all_scale[epoch_id] spec.err ./= epochs[irow, :scale]
end end
return spec return spec
end end
source = QSO{q1927p654}("1ES 1927+654 (Boller03)", Z, ebv=EBV); function analyze_boller03()
source.options[:min_spectral_coverage][:OIII_5007] = 0.35 file = opt.path * "/boller.dat"
# source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51" if !isfile(file)
#= source = QSO{q1927p654}("1ES 1927+654 (Boller03)", opt.z, ebv=opt.ebv);
Resolution: from Sect 3 it is 6 / 5007 * 3e5 ~ 360 km/s source.options[:min_spectral_coverage][:OIII_5007] = 0.35
however with this value the [OIII] FWHM is 100 km/s (limit) # source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"
I tried 180, 200 and 250, and the [OIII] luminosity do not change #=
=# Resolution: from Sect 3 it is 6 / 5007 * 3e5 ~ 360 km/s
spec = read_spec("B03", resolution=200.) however with this value the [OIII] FWHM is 100 km/s (limit)
add_spec!(source, spec); I tried 180, 200 and 250, and the [OIII] luminosity do not change
res = qsfit(source); =#
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], spec = read_spec("B03", resolution=200.)
filename="$(path)/results_boller.html") add_spec!(source, spec);
println(res.model[:OIII_5007].norm.val ) res = qsfit(source);
# 0.0002904273540659144 viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
println(res.model[:OIII_5007].norm.val / res.model[:galaxy].norm.val ) filename="$(opt.path)/results_boller.html")
# 60.25349678748317 JLD2.save_object(file, res)
else
res = JLD2.load_object(file)
if !isfile("$(path)/scale_fwhm.dat") end
all_scale = fill(1., length(epoch_filenames)) @assert res.model[:OIII_5007].norm.val == 0.00029042977121756283
all_fwhm = fill(NaN,length(epoch_filenames)) @assert res.model[:OIII_5007].norm.val / res.model[:galaxy].norm.val == 60.07373299362222
for id in 1:length(epoch_filenames) return res
spec = read_spec(id) end
all_scale[id] *= median(spec.flux) / 25000
try function estimate_SE_scale!()
source = QSO{q1927p654}("1ES 1927+654 ($id)", Z, ebv=EBV); 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)
epochs[irow, :scale] = median(spec.flux) / 25000
source = QSO{q1927p654}("1ES 1927+654 ($id)", opt.z, ebv=opt.ebv);
source.options[:min_spectral_coverage][:OIII_5007] = 0.5 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" # source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"
spec = read_spec(id) spec = read_spec(id)
add_spec!(source, spec); add_spec!(source, spec);
res = qsfit(source); res = qsfit(source);
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
filename="$(path)/results_$(id).html") filename="$(opt.path)/results_$(id).html")
# Update all_scale and store FWHM # Update all_scale and store FWHM
@info res.model[:OIII_5007].norm.val @info res.model[:OIII_5007].norm.val
all_scale[id] *= res.model[:OIII_5007].norm.val epochs[irow, :scale] *= res.model[:OIII_5007].norm.val
all_fwhm[ id] = res.model[:OIII_5007].fwhm.val epochs[irow, :OIII_fwhm] = res.model[:OIII_5007].fwhm.val
catch 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
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
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 end
out = hcat(out, OIII_norm)
end end
all_scale[.!isfinite.(all_fwhm)] .= NaN out = out[:,2:end]
JLD2.save_object("$(path)/scale_fwhm.dat", (all_scale, all_fwhm))
else @gp "set grid" :-
(all_scale, all_fwhm) = JLD2.load_object("$(path)/scale_fwhm.dat") 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 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 "set key bottom right" "set grid" xlab="Epoch" ylab="Value" :-
@gp :- :prenorm all_scale./1e-18 "w lp t '[OIII] norm.'" :- @gp :- :prenorm epochs.scale./1e-18 "w lp t '[OIII] norm.'" :-
@gp :- :prenorm all_fwhm./1000 "w lp t '[OIII] FWHM (x1000 km/s)" @gp :- :prenorm epochs.OIII_fwhm./1000 "w lp t '[OIII] FWHM (x1000 km/s)"
save(:prenorm, term="png size 800,600", output="$(path)/oiii_norm_fwhm.png") save(:prenorm, term="png size 800,600", output="$(opt.path)/oiii_norm_fwhm.png")
w = 10 .^(2:0.01:log10(400)); w = 10 .^(2:0.01:log10(400));
r = sqrt.(400^2 .- w.^2) ./ 2.355; r = sqrt.(400^2 .- w.^2) ./ 2.355;
...@@ -105,28 +173,29 @@ r = sqrt.(900^2 .- w.^2) ./ 2.355; ...@@ -105,28 +173,29 @@ r = sqrt.(900^2 .- w.^2) ./ 2.355;
# Ensure no sample is negative # Ensure no sample is negative
@gp :allepochs "set grid" xr=[3e3,1e4] cbr=[1,29] :- @gp :allepochs "set grid" xr=[3e3,1e4] cbr=[1,29] :-
@gp :- :allepochs cblabel="Epoch" xlab="Rest frame wavelength[A]" ylab="Flux density [arb.units]" @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) s = read_spec(id)
x = s.λ; x = s.λ;
y = s.flux; y = s.flux;
y .-= median(y) y .-= median(y)
y ./= maximum(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 end
@gp :- :allepochs yr=[1, 29] @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] @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], :vecchi => [17, 10, 21, 23],
:lowres => [8,9,10,11,12,14,15,16,18,21,22,23,24,25], :lowres => [8,9,10,11,12,14,15,16,18,21,22,23,24,25],
:highres => [5, 7, 13, 17, 20], :highres => [5, 7, 13, 17, 20],
:test => [7, 13],
:limited => [1,2,3,13], :limited => [1,2,3,13],
:all => collect(5:26)) :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 # Can't use resolution smaller than 150 km / s otherwise some line
# will be neglected because of spectral coverage # will be neglected because of spectral coverage
...@@ -134,55 +203,10 @@ resolution = Dict() ...@@ -134,55 +203,10 @@ resolution = Dict()
# :lowres => 350., # :lowres => 350.,
# :highres => 150.) # :highres => 150.)
job = :all analyze_job(:test)
chosen_epochs = dict_chosen_epochs[job] analyze_job(:highres)
Nloop = 6 analyze_job(:all)
(all_scale, all_fwhm) = JLD2.load_object("$(path)/scale_fwhm.dat") ddd
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
...@@ -256,7 +280,7 @@ for i in 1:nrow(tab) ...@@ -256,7 +280,7 @@ for i in 1:nrow(tab)
tab[i, :instr] = tmp[j, :Inst] tab[i, :instr] = tmp[j, :Inst]
end end
f = FITS("$(path)/1ES_1927p654_results_$(job).fits", "w") f = FITS("$(opt.path)/1ES_1927p654_results_$(job).fits", "w")
write(f, tab) write(f, tab)
close(f) close(f)
...@@ -298,7 +322,7 @@ for i in 1:nrow(tab) ...@@ -298,7 +322,7 @@ for i in 1:nrow(tab)
ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2" ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"
@gp :- domain(res.model[i])[:] res.model[i](:qso_cont) ss @gp :- domain(res.model[i])[:] res.model[i](:qso_cont) ss
end 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" @gp "set grid"
...@@ -310,4 +334,4 @@ for i in 1:length(chosen_epochs) ...@@ -310,4 +334,4 @@ for i in 1:length(chosen_epochs)
ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2" ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"
@gp :- x y .- Spline1D(x, y0)(5007.) ss @gp :- x y .- Spline1D(x, y0)(5007.) ss
end 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")