Skip to content
Commits on Source (2)
Hide whitespace changes
Inline Side-by-side
run.jl
View file @ 4d11b546
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,86 +8,157 @@ using Dates
using TextParse
import StatsBase, JLD2
path = "output-" * readlines(`git branch --show-current`)[1]
mkdir("$(path)")
epoch_filenames = Vector{String}()
for (root, dirs, files) in walkdir("AT2018zf")
for file in files
if file[end-3:end] == ".dat"
push!(epoch_filenames, root * "/" * file)
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"
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);
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"
#=
Resolution: from Sect 3 it is 6 / 5007 * 3e5 ~ 360 km/s
however with this value the [OIII] FWHM is 100 km/s (limit)
I tried 180, 200 and 250, and the [OIII] luminosity do not change
=#
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
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)
spec = read_spec(id)
all_scale[id] *= median(spec.flux) / 25000
try
source = QSO{q1927p654}("1ES 1927+654 ($id)", 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"
#=
Resolution: from Sect 3 it is 6 / 5007 * 3e5 ~ 360 km/s
however with this value the [OIII] FWHM is 100 km/s (limit)
I tried 180, 200 and 250, and the [OIII] luminosity do not change
=#
spec = read_spec("B03", resolution=200.)
add_spec!(source, spec);
res = qsfit(source);
viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
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
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)
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[: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
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
out = hcat(out, OIII_norm)
end
all_scale[.!isfinite.(all_fwhm)] .= NaN
JLD2.save_object("$(path)/scale_fwhm.dat", (all_scale, all_fwhm))
else
(all_scale, all_fwhm) = JLD2.load_object("$(path)/scale_fwhm.dat")
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,10 @@ 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)
analyze_job(:all)
ddd
......@@ -256,7 +280,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 +322,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 +334,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")