-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathPlotCand_v2.py
156 lines (132 loc) · 6.77 KB
/
PlotCand_v2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
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
150
151
152
153
154
155
#!/usr/bin/env python
# A separate function to extract and plot
# heimdall candidate
# This script is a modified version of the heimdall plotting scipt 'trans_freq_time.py'
#
import os,sys,math
import numpy as np
import glob
from itertools import chain
from os.path import basename
from itertools import tee, izip
def pairwise(iterable):
"s -> (s0,s1), (s1,s2), (s2, s3), ..."
a, b = tee(iterable)
next(b, None)
return izip(a, b)
def plotParaCalc(snr,filter,dm,fl,fh,tint):
#Extract block factor plot in seconds
extimefact = 1.0
# Total extract time Calc
# Extract according to the DM delay
cmd = 'dmsmear -f %f -b %f -n 2048 -d ' % (fl+(fh-fl)/2,fh-fl) + str(dm) + " -q 2>&1 "
p = os.popen(cmd)
cand_band_smear = p.readline().strip()
p.close()
extime = extimefact/2 + extimefact*float(cand_band_smear)
if extime < 1.0: extime = 1.0
# Tbin calc
# For Filter widths startting from 2^0 to 2^12=4096
#widths = [2048,2048,2048,1024,1024,512,512,256,256,128,128,64,32]
#tbin = widths[filter]
bin_width = tint * (2 ** filter)
#So that we have at least 4 bins on pulse
if filter <= 4 and snr > 20:
tbin = 4*int(extime / bin_width)
else:
tbin = 2*int(extime / bin_width)
if tbin < 16:
tbin = 16
#Fbin Calc
fbin = int(round(math.pow(float(snr)/4.0,2)))
if fbin < 16:
fbin = 16
fbin_base2 = int(round(math.log(fbin,2)))
fbin = pow(2,fbin_base2)
if fbin > 512:
fbin = 512
# Fraction of extraction to plot each time calc
if tbin>512:
frac = np.linspace(0,1,np.ceil(tbin/512.0))
else:
frac = np.array([0,1])
return tbin,fbin,extime,frac
def extractPlotCand(fil_file,frb_cands,noplot,fl,fh,tint,Ttot,kill_time_range,kill_chans,source_name):
# Half of this time will be subtracted from the Heimdall candidate time
extimeplot = 1.0
if(noplot is not True):
if(frb_cands.size >= 1):
if(frb_cands.size>1):
frb_cands = np.sort(frb_cands)
frb_cands[:] = frb_cands[::-1]
if(frb_cands.size==1): frb_cands = [frb_cands]
for indx,frb in enumerate(frb_cands):
time = frb['time']
dm = frb['dm']
filter = frb['filter']
width = tint * (2 ** filter)*(10**3) # Width in msec
snr = frb['snr']
tbin,fbin,extime,frac=plotParaCalc(snr,filter,dm,fl,fh,tint)
#print tbin,fbin,extime,frac
stime = time-(extimeplot*0.1) # Go 10% back data
if(stime<0): stime = 0
if(stime+extime>=Ttot): extime=Ttot-stime
#if(any(l<=stime<=u for (l,u) in kill_time_ranges)):
if(any(l<=time<=u for (l,u) in kill_time_range)):
print "Candidate inside bad-time range"
else:
if(indx<1000):
candname = '%04d' % (indx) + "_" + '%.3f' % (time) + "sec_DM" + '%.2f' % (dm)
cmd = "dspsr -cepoch=start -N %s" % (source_name) + \
" -b " + str(tbin) + \
" -S " + str(stime) + \
" -c " + str(extime) + \
" -T " + str(extime) + \
" -D " + str(dm) + \
" -O " + candname + " -e ar " + \
fil_file
print cmd
os.system(cmd)
# If kill_chans, do first manual and then an automatic smoothing for remaining channels
temp = ""
if kill_chans:
for k in kill_chans:
if(k!=2048): temp = temp +" "+str(k)
cmd = "paz -z \"" + temp + "\" -m %s.ar" % (candname)
print cmd
os.system(cmd)
cmd = "paz -r -b -L -m %s.ar" % (candname)
os.system(cmd)
# Correct the variable baseline, this script writes out .norm files
cmd = "running_mean_sub %s.ar" % (candname)
os.system(cmd)
for i,j in pairwise(frac):
cmd = "psrplot -p F -j 'D, F %d' " % (int(fbin)) + \
" -c 'flux:below:l = SNR: %.2f'" % (float(snr)) + \
" -c 'flux:below:r = Wid: %.2f ms'" % (float(width)) + \
" -c x:unit=ms -c above:c='' " + \
" -c 'x:range=(%f,%f)' " % (i,j) + \
" -c 'freq:cmap:map=heat' " + \
" -c 'freq:crop=0.9' " + \
" -D %s_%.2f.ps/cps %s.norm" % (candname,i,candname)
print cmd
os.system(cmd)
cmd = "gs -sDEVICE=pdfwrite -dNOPAUSE -dBATCH -dSAFER -sOutputFile=%s_frb_cand.pdf *sec*DM*.ps" % (source_name)
print cmd
os.system(cmd)
else:
print "No candidate found"
return
if __name__ == "__main__":
fil_file = str(sys.argv[1]) # Filterbank file
FinalList = str(sys.argv[2]) # Final list of candidate (output of FRB_detector_Fast.py)
frb_cands = np.loadtxt(FinalList,dtype={'names': ('snr','time','samp_idx','dm','filter','prim_beam'),'formats': ('f4', 'f4', 'i4','f4','i4','i4')})
fl = 1100
fh = 1500
noplot=0
tint=0.000128
Ttot = 20 # Total length of the file
kill_time_range=[]
kill_chans=[]
source_name="Fake"
extractPlotCand(fil_file,frb_cands,noplot,fl,fh,tint,Ttot,kill_time_range,kill_chans,source_name)