run.jl

using Pkg
Pkg.activate(".")

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


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)
        end
    end
end


# Estimate norm. and FWHM of [OIII] in all epochs, and normalize all
# spectra at 5100A.
all_scale  = 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];
    scale = median(y)
    y ./= scale
    j = findall(isfinite.(x));
    x = x[j];
    y = y[j];
    model, bestfit = fitline(x, y)
    scale *= bestfit[:line].norm.val
    y ./=    bestfit[:line].norm.val
    model, bestfit = fitline(x, y)
    show(bestfit)

    push!(all_scale, scale)
    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_scale./1e-14 "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"





function read_spec(epoch_id; kw...)
    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_scale[epoch_id]
    end
    spec.err .= 0.05 .* spec.flux;
    return spec
end


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




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
# 657.1651386497541


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


chosen_epoch = 17
file = epoch_filenames[chosen_epoch]
source = QSO{q1927p654}("1ES 1927+654 ($chosen_epoch)", Z, ebv=EBV);
source.options[:min_spectral_coverage][:OIII_5007] = 0.5
source.options[:oiii_galaxy] = 60.65
# source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"
spec = read_spec(chosen_epoch)
add_spec!(source, spec);
(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
resolution = Dict(
    :lowres  => 350.,
    :highres => 150.)

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
    source.options[:oiii_galaxy] = 60.65
    @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)
end



@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")



        if  haskey(model[id], :galaxy)  &&
            haskey(model[id], :OIII_5007)
            model[id][:galaxy].norm.fixed = true
            patch!(model[id]) do m
                m[id][:galaxy].norm = m[id][:OIII_5007].norm / 59.08
            end
        end