Added 2019 FLOYDS epochs. Updated plots (5231b9a0) · Commits · Giorgio Calderone / CL_1ES_1927p654

funcs.jl

0 → 100644

+198 −0

Original line number	Diff line number	Diff line

		function read_epochs()
		out = DataFrame(filename=String[], date=String[], day=Int[],
		telescope=String[], instr=String[],
		scale=Float64[], scale_rel_unc=Float64[],
		OIII_fwhm=Float64[], OIII_fwhm_unc=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, (filename, date, 0, "", "", NaN, NaN, NaN, NaN))
		end
		end
		end
		@assert issorted(out.date)

		dd = Date.(out.date, Ref("yyyymmdd"))
		out[!, :day] .= getproperty.(dd - dd[1], :value)
		@assert issorted(out.day)

		(d, c) = csvread("tabula-2019-Trakhtenbrot-mod.csv", ',', header_exists=true)
		tmp = DataFrame(collect(d), Symbol.(c))
		tmp[!, :Date] .= string.(tmp.Date)

		for i in 1:nrow(out)
		j = findfirst(tmp.Date .== out[i, :date])
		if isnothing(j)
		@warn "Can't find epoch $(out[i, :date]) in tabula-2019-Trakhtenbrot"
		continue
		end
		out[i, :telescope] = tmp[j, :Telescope]
		out[i, :instr] = tmp[j, :Inst]
		end

		out.day .+= 72 # 2018-03-06 corresponds to day 72

		return out
		end


		# function read_spec(epoch::DataFrameRow; kw...)
		# irow = findfirst(epochs.date .== date)
		# @assert !isnothing(irow)
		# file = epochs[irow, :filename]
		# spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$date: $file", kw...);
		# spec.flux ./= epochs[irow, :scale]
		# spec.err ./= epochs[irow, :scale]
		# return spec
		# end

		function read_spec(epoch::DataFrameRow; kw...)
		file = epoch.filename
		spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$(epoch.date): $file", kw...);
		spec.flux ./= epoch.scale
		spec.err ./= epoch.scale
		return spec
		end


		function boller03_analysis(opt)
		file = opt.path * "/boller.dat"
		if !isfile(file)
		source = QSO{q1927p654}("1ES 1927+654 (Boller+03)", opt.z, ebv=opt.ebv);
		#=
		Resolution: from Sect 3 it is 6 / 5007 * 3e5 ~ 360 km/s
		however with this value the [OIII] FWHM hits the limit at 100 km/s
		I tried 180, 200 and 250, and the [OIII] luminosity do not change
		=#
		spec = Spectrum(Val(:ASCII), "boller2003.txt", columns=[1,2]; label="Boller03", resolution=200.)
		source.options[:min_spectral_coverage][:OIII_5007] = 0.35
		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_single_epoch_scale!(opt, epochs)
		file = opt.path * "/" * "SE_scale.dat"
		if !isfile(file)
		# Estimate first guess calibration factor
		for irow in 1:nrow(epochs)
		epochs[irow, :scale] = 1.
		spec = read_spec(epochs[irow, :])
		epochs[irow, :scale] = median(spec.flux)
		end
		epochs[:, :scale] .= median(epochs.scale) / 25_000
		@assert all(epochs[:, :scale] .== 5.7866e-20)

		for irow in 1:nrow(epochs)
		date = epochs[irow, :date]
		source = QSO{q1927p654}("1ES 1927+654 ($(date))", opt.z, ebv=opt.ebv);
		source.options[:min_spectral_coverage][:OIII_5007] = 0.35
		spec = read_spec(epochs[irow, :])
		add_spec!(source, spec);
		res = qsfit(source);
		viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
		filename="$(opt.path)/results_$(date).html")

		# Update scale and OIII_fwhm
		epochs[irow, :scale] *= res.model[:OIII_5007].norm.val
		epochs[irow, :scale_rel_unc] = res.model[:OIII_5007].norm.unc / res.model[:OIII_5007].norm.val
		epochs[irow, :OIII_fwhm] = res.model[:OIII_5007].fwhm.val
		epochs[irow, :OIII_fwhm_unc] = res.model[:OIII_5007].fwhm.unc
		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 multi_epoch_analysis!(opt, _epochs, job; Nloop = 6)
		estimate_single_epoch_scale!(opt, _epochs)
		epochs = _epochs[opt.subsets[job], :]

		out = fill(NaN, nrow(epochs))
		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 irow in 1:nrow(epochs)
		spec = read_spec(epochs[irow, :])
		add_spec!(source, spec);
		@gp :- :zoom xr=[4750,5150] spec.λ ./ (1 + source.z) spec.flux "w l t '$(epochs[irow, :date])'"
		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
		@info file
		res = JLD2.load_object(file);
		println(loop, " ", res.fitres.gofstat / res.fitres.dof, " ", res.fitres.elapsed / 3600.)
		end

		OIII_norm = Float64[]
		for i in 1:nrow(epochs)
		push!(OIII_norm, res.multi[i][:OIII_5007].norm.val)
		epochs[i, :scale] *= ((1 + res.multi[i][:OIII_5007].norm.val) / 2)
		end
		out = hcat(out, OIII_norm)

		f = open("$(opt.path)/scales_$(job)_$(loop).txt", "w")
		show(f, epochs)
		close(f)
		end
		out = out[:,2:end]

		@gp "set grid" "set key outside top horiz" xlab="Iteration" ylab="[OIII]λ5007 norm. (arb. units)" :-
		for i in 1:nrow(epochs)
		@gp :- 1:Nloop out[i, :] "w lp t '$(epochs[i, :date])'" :-
		end
		@gp :- xr=[1,6]
		save(term="pngcairo size 1000,700", output="$(opt.path)/calib_$(job).png")

		@gp :- yr=[0.96, 1.04]
		save(term="pngcairo size 1000,700", output="$(opt.path)/calib_$(job)_zoom.png")

		_epochs[opt.subsets[job], :scale] .= epochs.scale
		return out
		end


		function cont_at(bestfit, λ)
		A1 = bestfit[:qso_cont].norm.val - bestfit[:qso_cont].norm.unc
		A = bestfit[:qso_cont].norm.val
		A2 = bestfit[:qso_cont].norm.val + bestfit[:qso_cont].norm.unc
		if A1 > A2
		A1, A2 = A2, A1
		end
		@assert A1 <= A2

		B1 = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val - bestfit[:qso_cont].alpha.unc)
		B = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val)
		B2 = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val + bestfit[:qso_cont].alpha.unc)
		if B1 > B2
		B1, B2 = B2, B1
		end
		@assert B1 <= B2

		dd = extrema([A1B1, A1B2, A2B1, A2B2])

		return λ * A * B, λ * (dd[2]-dd[1])
		end

tabula-2019-Trakhtenbrot-mod.csv

0 → 100644

+29 −0

Original line number	Diff line number	Diff line
		Date,MJD start,dt,Telescope,Inst,Exp
		20180306,58183.5,72,MDM Hiltner 2.4m,OSMOS,600
		20180308,58185.5,74,MDM Hiltner 2.4m,OSMOS,1800
		20180309,58186.5,75,MDM Hiltner 2.4m,OSMOS,2700
		20180313,58190.0,79,Liverpool Telescope 2m,SPRAT,600
		20180323,58200.5,89,APO 3.5m,DIS,300
		20180423,58231.6,119,Las Cumbres 2m,FLOYDS,1800
		20180424,58232.5,120,Shane 3m,Kast,1500
		20180507,58245.6,133,Las Cumbres 2m,FLOYDS,1800
		20180514,58252.5,140,Las Cumbres 2m,FLOYDS,1800
		20180528,58266.5,153,Las Cumbres 2m,FLOYDS,1800
		20180603,58272.6,159,Las Cumbres 2m,FLOYDS,1800
		20180611,58280.6,167,Las Cumbres 2m,FLOYDS,1800
		20180613,58282.5,169,APO 3.5m,DIS,600
		20180624,58293.3,180,Las Cumbres 2m,FLOYDS,1800
		20180628,58297.5,184,Las Cumbres 2m,FLOYDS,3600
		20180706,58305.5,192,Las Cumbres 2m,FLOYDS,3600
		20180716,58315.5,201,Keck,LRIS,900
		20180717,58316.5,202,Las Cumbres 2m,FLOYDS,3600
		20180727,58326.5,212,Las Cumbres 2m,FLOYDS,3600
		20180811,58341.0,227,Hale200 inch,DBSP,1200
		20180812,58342.4,228,Las Cumbres 2m,FLOYDS,3600
		20180908,58369.3,254,Las Cumbres 2m,FLOYDS,3600
		20181113,58435.2,319,Las Cumbres 2m,FLOYDS,3600
		20180828,58358.4,244,HST,COS,2866
		20180828,58358.5,244,HST,STIS,2740
		20190319,0,119,Las Cumbres 2m,FLOYDS,1800
		20190406,0,119,Las Cumbres 2m,FLOYDS,1800
		20190519,0,119,Las Cumbres 2m,FLOYDS,1800