In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from astropy.constants import c
from astropy.convolution import Gaussian1DKernel, convolve
import sys , os
sys.path.append("../../src/lines")
import lineTools as lt
### working dir. and files
wd = "/home/stephane/Science/RadioGalaxy/ALMA/absorptions/analysis/a/"
os.chdir(wd)
datadir = "dataSpecAll/"
dbline = "lineAll.db"
transfile = "splatalogue.csv"
dirplot = "plots/"
In [2]:
trans = pd.read_csv(transfile, sep=":")
nline = len(trans)
class color:
PURPLE = '\033[95m'
CYAN = '\033[96m'
DARKCYAN = '\033[36m'
BLUE = '\033[94m'
GREEN = '\033[92m'
YELLOW = '\033[93m'
RED = '\033[91m'
BOLD = '\033[1m'
UNDERLINE = '\033[4m'
END = '\033[0m'
print("%s %d %s transition lines from splatalogue"%(color.BOLD, nline, color.END))
In [3]:
al = lt.analysisLines(dbline)
res = al.getInfoDB(verbose = False)
## Frequency range to check for the line frequency at rest applying the redshift (one source can have several z)
##
df = 0.1
ndetected = 0
snmin = 3.0
sourcewithnoredshift =[]
linesFound = []
sumtable = ""
print("## Searching for lines in %s.."%(transfile))
print("## Tolerance : %3.3f MHz"%(df*1000.))
print("## S/N min : %3.2f"%(snmin))
for i in range(nline):
freqLine = trans.loc[i,'Freq']
nameTrans = trans.loc[i,'Species']
transLow = freqLine - df
transHigh = freqLine + df
for s in res:
lines = al.findLinesSource(s[0])
if len(lines) > 0 and not isinstance(s[1], int):
for z in s[1]:
if not np.isnan(z):
for l in lines:
frest = l[4] * (1. + z)
if frest > transLow and frest < transHigh and l[2] >= snmin:
resDict = {}
print("## Line found..")
print("### Transition : %s %s %s - %3.5f GHz"%(color.BOLD, nameTrans, color.END, freqLine))
print("### Source : %s %s %s "%(color.BOLD, s[0], color.END))
print("### Redshift : %f"%(z))
print("### Detection : %3.3f GHz (rest) - %3.3f GHz (sky)"%( frest, l[4]))
print("### S/N = %s%3.2f%s - I = %s%3.3f%s mJy"%(color.BOLD, l[2], color.END, color.BOLD, 1000*l[3],color.END))
print("### lineid : %s%d%s -- dataid : %s%d%s"%(color.BOLD, l[0], color.END, color.BOLD, l[1], color.END))
print("###")
ndetected += 1
## formatted
vel = c.value * 1e-3 * (float(frest) - freqLine) / freqLine
sumtable +="| %s | %2.4f | %s | %3.3f | %3.3f | %3.2f | %2.1f | %d | %d | | \n"%(s[0], z, nameTrans, freqLine, vel, 1000.*l[3], l[2] , l[1], l[0])
resDict['species'] = nameTrans
resDict['source'] = s[0]
resDict['frequency'] = freqLine
resDict['freqsky'] = l[4]
resDict['redshift'] = float(z)
resDict['freqatrest'] = float(frest)
resDict['peak'] = l[3]
resDict['snr'] = l[2]
resDict['dataid'] = l[1]
resDict['lineid'] = l[0]
## query the datafile
cmdsql = "select filedata FROM dataset WHERE dataid = '%d'"%(l[1])
resdb = al.query(cmdsql)
resDict['filedata'] = resdb[0][0]
linesFound.append(resDict)
print("## %s %d %s tentative detections. "%(color.BOLD, ndetected, color.END) )
print("## Done..")
print(sumtable)
In [4]:
def showText(posx , pys, text, ywidth = 1.):
posy = pys
dy = ywidth / len(text)
for t in text:
plt.text(posx,posy,t)
posy += dy
def plotSpectra(det , dv , amp , vwidth = 250. , lineplot = True):
"Plot the spectra with line detection and transition name"
plt.figure(figsize=(8.0,4.5))
plt.plot(dv,amp)
x1 = max(min(dv),-1. * vwidth)
x2 = min(max(dv),vwidth)
ind = np.where(np.abs(dv) < 200 )
ymin = min(amp[ind])
ymax = max(amp[ind])
y1 = ymin - (ymax-ymin)*0.25
y2 = ymax + (ymax-ymin)*0.1
dx = x2-x1
dy = y2-y1
plt.xlim(x1,x2)
plt.ylim(y1,y2)
txt = []
txt.append("Source : %s"%(det['source']))
txt.append("Redshift : %s"%(det['redshift']))
txt.append("Species : %s"%(det['species']))
txt.append("Freq. Sky : %3.3f GHz"%(det['freqsky']))
vel = c.value * 1e-3 * (det['freqatrest']- det['frequency']) / det['frequency']
txt.append("dv : %3.3f km/s"%(vel))
txt.append("Lineid : %d"%(det['lineid']))
txt.append("Dataid : %d"%(det['dataid']))
showText(x1+dx*0.01, y1+dy*0.01 , txt, (y2-y1)*0.4)
# xline= [det['frequency'], det['frequency']]
xline = [0., 0.]
xdetect = [vel, vel]
yline = [0., 2.]
if x1 <0 and x2>0 and lineplot:
plt.plot(xdetect, yline, "g--")
plt.plot(xline, yline, "r--")
plt.text((x1-x2)*0.02, y1 + (y2-y1) * 0.3, det['species'] , rotation = 90. )
figname = "%s%s-%s-%s.png"%(dirplot,det['source'],det['species'],det['lineid'])
plt.savefig(figname)
plt.show()
def plotSpectraPublication(vel, amp , info):
"Plot published-like plot of the spectra with txt array info"
plt.figure(figsize=(8.0,4.5))
plt.step(dv,amp, "b-" , linewidth = 1.0)
plt.xlabel ("v (km/s)")
plt.ylabel(r"$S^*$")
vwidth = 300.
x1 = info['velline'] - vwidth/ 2.0
x2 = info['velline'] + vwidth/ 2.0
ind = np.where(np.abs(dv) < 200 )
ymin = min(amp[ind])
ymax = max(amp[ind])
y1 = ymin - (ymax-ymin)*0.3
y2 = ymax + (ymax-ymin)*0.1
dx = x2-x1
dy = y2-y1
xdetect = info['velline']
xline = [xdetect,xdetect]
yline = [y1,y1 + dy * 0.15]
plt.plot(xline,yline,"g--")
xtxt = (x2+x1)/2. + dx*0.3
ytxt = y1 + dy * 0.05
plt.text(xtxt,ytxt,"%s"%(info['source']))
plt.text(xdetect-dx*0.03 , y1 + dy*0.18,info["species"], rotation = 90)
plt.xlim(x1,x2)
plt.ylim(y1,y2)
figname = "%s%s-%s-%s.pdf"%(dirplot,info['source'],info['species'],info['lineid'])
print(figname)
#plt.savefig(figname)
plt.show()
In [5]:
pl = lt.plotLines("fake", "fake", "fake")
for det in linesFound:
datafile = datadir + det['filedata'][2:]
freq , amp = pl.extractData(datafile)
amp = amp / np.mean(amp)
freq = freq * (1. + det['redshift'])
dv = c.value * 1e-3 * (freq - det['frequency']) / det['frequency']
## smoothing to 1 km/s
# Create kernel
ddv = abs(dv[1]-dv[0])
kdv = 0.5 / ddv ## 1km/s
g = Gaussian1DKernel(stddev=kdv)
z = convolve(amp, g)
plotSpectra(det,dv,z)
The detection are for following after checking each one.
| Source | redshift | Species | Frequency | Vel. | Peak | S/N | dataid | lineid | Comments |
|---|---|---|---|---|---|---|---|---|---|
| --- | --- | --- | GHz | km/s | mJy | --- | --- | --- | --- |
| J0219+0120 | 1.6230 | HDO | 255.050 | 59.075 | 32.59 | 4.8 | 4567 | 7042 | emission line, TBC |
| J0418+380 | 0.0485 | HNCv=0 | 362.630 | -29.982 | -967.54 | 7.4 | 1676 | 3989 | TBC |
| J0538-440 | 0.8940 | HCNv=0 | 443.116 | 33.434 | 12.72 | 3.2 | 6140 | 8744 | abs. line! but telluric?? |
| J051002+180041 | 0.4160 | H2Ov=0 | 474.689 | 35.129 | -26.21 | 6.6 | 167 | 167 | |
| J1159-2142 | 0.6170 | HDO | 157.404 | 7.946 | 22.36 | 3.9 | 6085 | 8645 | ?? |
| J1159-2142 | 0.6170 | HDO | 157.404 | -36.446 | -24.55 | 4.1 | 6085 | 8646 | ?? |
| J1337-1257 | 0.5390 | HDO | 364.907 | 32.025 | 18.54 | 5.6 | 6163 | 8805 | 1 line in 3 parts |
| J1337-1257 | 0.5390 | HDO | 364.907 | 49.853 | -67.40 | 23.2 | 6163 | 8806 | 1 line in 3 parts |
| J1337-1257 | 0.5390 | HDO | 364.907 | 68.831 | 13.73 | 5.0 | 6163 | 8807 | 1 line in 3 parts |
| J1832-2039 | 0.1033 | HDO | 263.832 | -43.877 | -175.91 | 3.4 | 5904 | 8242 | ?? |
| J1832-2039 | 0.1033 | HDO | 263.832 | -46.916 | -187.31 | 5.3 | 5920 | 8283 | good |
| J1832-2039 | 0.1033 | HDO | 263.832 | -60.033 | 97.83 | 3.2 | 5920 | 8284 | emission similar to previous |
| J1833-2103 | 0.1926 | C17O | 337.061 | 29.409 | -23.42 | 3.3 | 5712 | 7682 | TBC |
| J1833-2103 | 0.1926 | C17O | 337.061 | 29.276 | -22.96 | 3.9 | 5728 | 7743 | similar to previous only vel. |
| J1833-2103 | 0.8858 | SiOv=0 | 607.599 | -8.859 | -85.52 | 4.6 | 5777 | 7922 | |
| J1833-2103 | 0.8858 | SiOv=0 | 607.599 | -9.594 | -96.12 | 7.0 | 5793 | 7995 | same line as above |
| J1924-292 | 0.3526 | HDO | 143.727 | -21.000 | 2.33 | 3.1 | 5966 | 8451 | emission line |
| J2148+0657 | 0.7000 | NH3v=0 | 172.832 | 170.928 | 8.29 | 4.9 | 4204 | 6680 | TBC |
| J2148+0657 | 0.8490 | SiOv=0 | 173.688 | 61.019 | -10.53 | 4.3 | 4762 | 7206 | TBC , same line as HCO+ |
| J2148+0657 | 0.8479 | SiOv=0 | 173.688 | -117.369 | -10.53 | 4.3 | 4762 | 7206 | TBC , different z |
| J2148+0657 | 0.8983 | HCO+v=0 | 178.375 | -32.588 | -10.53 | 4.3 | 4762 | 7206 | TBC , good candidate, z! |
| J2148+0657 | 0.9765 | HDO | 185.709 | -9.668 | -10.53 | 4.3 | 4762 | 7206 | TBC , good candidate , z! |
| J2148+0657 | 0.8648 | NH3v=0 | 189.725 | -47.892 | 8.29 | 4.9 | 4204 | 6680 | looks an abs. |
| J2148+0657 | 0.8983 | C18O | 439.089 | -21.928 | 9.49 | 4.2 | 2041 | 4869 | same detectino withe 2 species |
| J2148+0657 | 0.8983 | H2Ov=0 | 439.151 | -64.270 | 9.49 | 4.2 | 2041 | 4869 |
In [6]:
## plot the selected lines (publishable form)
lines = pd.read_csv("lineSelectedEG.csv", sep="|")
for i in range(len(lines)):
txtinfo = {}
cmdsql = "select filedata FROM dataset WHERE dataid = '%s'"%(lines.loc[i,'dataid'])
resdb = al.query(cmdsql)
datafile = resdb[0][0]
datafile = datadir + datafile
freq , amp = pl.extractData(datafile)
amp = amp / np.mean(amp)
freq = freq * (1. + lines.loc[i,'redshift'])
dv = c.value * 1e-3 * (freq - lines.loc[i,"Frequency"]) / lines.loc[i,'Frequency']
txtinfo['source'] = lines.loc[i,'Source']
txtinfo['species'] = lines.loc[i,'Species']
txtinfo['frequency'] = lines.loc[i,'Frequency']
txtinfo['velline'] = lines.loc[i,'Vel']
txtinfo['lineid'] = lines.loc[i,'lineid']
## smoothing to 1 km/s
# Create kernel
ddv = abs(dv[1]-dv[0])
kdv = 0.1 / ddv ## 1km/s
g = Gaussian1DKernel(stddev=kdv)
z = convolve(amp, g)
plotSpectraPublication(dv, z , txtinfo)
In [7]:
al = lt.analysisLines(dbline)
res = al.getInfoDB(verbose = False)
sourcewithoutz =[]
for s in res:
redshiftFound = False
if not isinstance(s[1], int):
for z in s[1]:
if not np.isnan(z):
redshiftFound = True
else:
redshiftFound = True
if not redshiftFound:
sourcewithoutz.append(s[0])
print("## Sources with no redshift known:")
print(sourcewithoutz)
## lines of these sources.
for s in sourcewithoutz:
lines = al.findLinesSource(s)
print("## Source: %s"%(s))
print("### Lines: %d"%(len(lines)))
In [8]:
snrmin = 4.0
linesFound = []
sumtable = ""
for s in sourcewithoutz :
lines = al.findLinesSource(s)
if len(lines) > 0:
print("## Source: %s"%(s))
for l in lines:
if l[2] >= snrmin:
z = 0.0
nameTrans = "unknown"
frest = l[4] * (1. + z)
freqLine = frest
resDict = {}
vel = c.value * 1e-3 * (float(frest) - freqLine) / freqLine
sumtable +="| %s | %2.4f | %s | %3.3f | %3.3f | %3.2f | %2.1f | %d | %d | | \n"%(s, z, nameTrans, freqLine, vel, 1000.*l[3], l[2] , l[1], l[0])
resDict['species'] = nameTrans
resDict['source'] = s
resDict['frequency'] = freqLine
resDict['freqsky'] = l[4]
resDict['redshift'] = float(z)
resDict['freqatrest'] = float(frest)
resDict['peak'] = l[3]
resDict['snr'] = l[2]
resDict['dataid'] = l[1]
resDict['lineid'] = l[0]
## query the datafile
cmdsql = "select filedata FROM dataset WHERE dataid = '%d'"%(l[1])
resdb = al.query(cmdsql)
resDict['filedata'] = resdb[0][0]
linesFound.append(resDict)
#print(sumtable)
In [9]:
pl = lt.plotLines("fake", "fake", "fake")
for det in linesFound:
datafile = datadir + det['filedata'][2:]
freq , amp = pl.extractData(datafile)
amp = amp / np.mean(amp)
freq = freq * (1. + det['redshift'])
dv = c.value * 1e-3 * (freq - det['frequency']) / det['frequency']
## smoothing to 1 km/s
# Create kernel
ddv = abs(dv[1]-dv[0])
kdv = 0.5 / ddv ## 1km/s
g = Gaussian1DKernel(stddev=kdv)
z = convolve(amp, g)
plotSpectra(det,dv,z , lineplot = False)