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) .+ 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 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 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) #scale /= bestfit[:line].norm.val show(bestfit) 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") function read_spec(epoch_id; kw...) file = epoch_filenames[epoch_id] println(file) spec = Spectrum(Val(:ASCII), file, columns=[1,2]; kw...); spec.flux .+= all_offset[epoch_id] spec.flux .*= all_scale[epoch_id] spec.flux .+= global_offset 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") 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[: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]) resolution = Dict( :lowres => 300., :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[:use_host_template] = false @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