Newer
Older
Z = 0.019422
EBV = 0.077
using Revise, Statistics, Serialization, Dierckx
using QSFit, GFit, Gnuplot, GFitViewer
using CL_1ES_1927p654
function fitline(x, y, err=[]; center=NaN, sigma=NaN)
domain = Domain(x)
isfinite(center) || (center = mean(x))
isfinite(sigma) || (sigma = (maximum(x) - minimum(x)) / 10.)
(length(err) == 0) && (err = (y .* 0.1) .+ maximum(y) / 10)
model = Model(domain, Reducer(sum, [:continuum, :line]),
:continuum => GFit.OffsetSlope(0., mean(x), 0.),
:line => GFit.Gaussian(1., center, sigma))
bestfit = fit!(model, Measures(y, err))
return (model, bestfit)
end
for (root, dirs, files) in walkdir("AT2018zf")
for file in files
if file[end-3:end] == ".dat"
global_offset = 0.
all_offset = Vector{Float64}()
all_scale = Vector{Float64}()
all_norm = Vector{Float64}()
all_fwhm = Vector{Float64}()
for i in 1:length(epoch_filenames)
s = Spectrum(Val(:ASCII), epoch_filenames[i], columns=[1,2]);
x = s.λ;
j = findall(5060 .< x .< 5150);
x = x[j];
y = s.flux[j];
offset = 0. # -median(y);
y .+= offset;
#scale = 1 / maximum(y[5050 .< x .< 5150]) # quantile(y, 0.95);
scale = 1 / Spline1D(x ./ (1+Z), y)(5100.)
y .*= scale;
j = findall(isfinite.(x));
x = x[j];
y = y[j];
model, bestfit = fitline(x, y)
push!(all_offset, offset)
push!(all_scale, scale)
push!(all_norm, bestfit[:line].norm.val)
push!(all_fwhm, bestfit[:line].sigma.val / 5007 * 3e5 * 2.355)
@gp tit=string(i) x y "w lp notit" x model() "w l notit"
(readline() == "q") && break
end
@gp :prenorm "set key bottom right" "set grid" xlab="Epoch" ylab="Value" :-
@gp :- :prenorm all_norm "u (\$0+1):1 w lp t '[OIII] norm.'" all_fwhm./100 "u (\$0+1):1 w lp t '[OIII] FWHM (x100 km/s)"
save(:prenorm, term="png size 800,600", output="output/oiii_norm_fwhm.png")
w = 10 .^(2:0.01:log10(400));
r = sqrt.(400^2 .- w.^2) ./ 2.355;
@gp w r "w l t '400 km/s'"
w = 10 .^(2:0.01:log10(900));
r = sqrt.(900^2 .- w.^2) ./ 2.355;
@gp :- w r "w l t '900 km/s'"
@gp :- "set grid" xlab="Actual FWHM" ylab="Instr. resolution"
if epoch_id == "B03"
file = "boller2003.txt"
else
file = epoch_filenames[epoch_id]
end
println(file)
spec = Spectrum(Val(:ASCII), file, columns=[1,2]; kw...);
if isa(epoch_id, Number)
spec.flux .+= all_offset[epoch_id]
spec.flux .*= all_scale[epoch_id]
spec.flux .+= global_offset
end
spec.err .= 0.05 .* spec.flux;
return spec
end
# Ensure no sample is negative, or adjust global_offset.
@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)
s = read_spec(id)
x = s.λ;
y = s.flux;
@gp :- :allepochs x./(1+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="output/all_epochs.png")
@gp :- :allepochs xr=[4900, 5100]
save(:allepochs, term="png size 800,2000", output="output/all_epochs_zoom.png")
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
s = read_spec("B03")
x = s.λ;
j = findall(5060 .< x .< 5150);
x = x[j];
y = s.flux[j];
model, bestfit = fitline(x, y)
@gp x y "w l"
@gp :- x model() "w l"
bestfit[:line].sigma.val / 5007 * 3e5 * 2.355
source = QSO{q1927p654}("1ES 1927+654 (Boller03)", Z, ebv=EBV);
source.options[:min_spectral_coverage][:OIII_5007] = 0.35
source.options[:oiii_galaxy] = NaN
# 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);
(model, bestfit) = fit(source);
viewer(model, source.data, bestfit, selcomps=[:qso_cont, :galaxy, :balmer])
println( bestfit[:OIII_5007].norm.val )
# 0.0002909406810630153
println( bestfit[:OIII_5007].norm.val / bestfit[:galaxy].norm.val )
# 60.6503967542268
file = epoch_filenames[chosen_epoch]
source = QSO{q1927p654}("1ES 1927+654 ($chosen_epoch)", Z, ebv=EBV);
# source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"
(model, bestfit) = fit(source);
viewer(model, source.data, bestfit, selcomps=[:qso_cont, :galaxy, :balmer])
dict_chosen_epochs = 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])
# Can't use resolution smaller than 150 km / s otherwise some line
# will be neglected because of spectral coverage
for job in [:highres, :lowres]
chosen_epochs = dict_chosen_epochs[job]
source = QSO{q1927p654}("1ES 1927+654", Z, ebv=EBV);
source.options[:min_spectral_coverage][:OIII_5007] = 0.5
@gp :zoom "set grid" :-
for id in 1:length(chosen_epochs)
spec = read_spec(chosen_epochs[id], resolution=resolution[job])
add_spec!(source, spec);
@gp :- :zoom xr=[4750,5150] spec.λ ./ (1 + source.z) spec.flux "w l t '$(chosen_epochs[id])'"
end
(model, bestfit) = multi_fit(source);
viewer(model, source.data, bestfit, selcomps=[:qso_cont, :galaxy, :balmer])
run(`cp /tmp/gfitviewer.html output/results_$(job).html`)
serialize("output/results_$(job).dat", bestfit)
@gp( [bestfit[id][:qso_cont].alpha.val for id in 1:length(chosen_epochs)],
[bestfit[id][:br_Ha].norm.patched for id in 1:length(chosen_epochs)] ./ all_scale[chosen_epochs],
"w l")
@gp(:-,
[bestfit[id][:qso_cont].alpha.val for id in 1:length(chosen_epochs)],
[bestfit[id][:br_Hb].norm.patched for id in 1:length(chosen_epochs)],
"w l")
@gp( [bestfit[id][:qso_cont].alpha.val for id in 1:length(chosen_epochs)],
[bestfit[id][:qso_cont].norm.patched for id in 1:length(chosen_epochs)],
"w l")
@gp( [bestfit[id][:qso_cont].norm.val for id in 1:length(chosen_epochs)],
[bestfit[id][:br_Ha].norm.patched for id in 1:length(chosen_epochs)],
"w l")