Commits
Laszlo Szucs authored 600a7638580
1 1 | import re |
2 + | import pdb |
2 3 | |
3 4 | import pipeline.infrastructure as infrastructure |
4 5 | import pipeline.infrastructure.filenamer as filenamer |
6 + | import pipeline.infrastructure.casatools as casatools |
7 + | from pipeline.infrastructure import casa_tasks |
8 + | import pipeline.domain.measures as measures |
5 9 | |
6 10 | from .imageparams_base import ImageParamsHeuristics |
7 11 | |
8 12 | LOG = infrastructure.get_logger(__name__) |
9 13 | |
10 14 | |
11 15 | class ImageParamsHeuristicsVLA(ImageParamsHeuristics): |
12 16 | |
13 17 | def __init__(self, vislist, spw, observing_run, imagename_prefix='', proj_params=None, contfile=None, |
14 18 | linesfile=None, imaging_params={}): |
15 19 | ImageParamsHeuristics.__init__(self, vislist, spw, observing_run, imagename_prefix, proj_params, contfile, |
16 20 | linesfile, imaging_params) |
17 21 | self.imaging_mode = 'VLA' |
18 22 | |
19 23 | def robust(self): |
20 24 | """See PIPE-680 and CASR-543""" |
21 25 | return 0.5 |
22 26 | |
23 27 | def uvtaper(self, beam_natural=None, protect_long=None): |
24 28 | return [] |
25 29 | |
30 + | def uvrange(self): |
31 + | """ |
32 + | Restrict uvrange in case of very extended emission. |
33 + | |
34 + | If the amplitude of the shortest 5 per cent of the covered baselines |
35 + | is more than 2 times that of the 50-55 per cent baselines, then exclude |
36 + | the shortest 5 per cent of the baselines. |
37 + | |
38 + | See PIPE-681 and CASR-543. |
39 + | |
40 + | :return: None or string in the form of '> {x}klambda', where |
41 + | {x}=0.05*max(baseline) |
42 + | """ |
43 + | def get_mean_amplitude(vis, uvrange=None, axis='amplitude', spw=None): |
44 + | stat_arg = {'vis': vis, 'uvrange': uvrange, 'axis': axis, |
45 + | 'useflags': True, 'spw': spw} |
46 + | job = casa_tasks.visstat(**stat_arg) |
47 + | stats = job.execute(dry_run=False) # returns stat in meter |
48 + | |
49 + | # loop through keys to determine average from all spws |
50 + | mean = 0.0 |
51 + | for k, v in stats.items(): |
52 + | mean = mean + v['mean'] |
53 + | |
54 + | return mean / float(len(stats)) |
55 + | # |
56 + | qa = casatools.quanta |
57 + | |
58 + | # Can it be that more than one visibility (ms file) is used? |
59 + | vis = self.vislist[0] |
60 + | ms = self.observing_run.get_ms(vis) |
61 + | spw_str = ','.join([str(spw) for spw in self.spwids]) |
62 + | |
63 + | # Determine the largest covered baseline in klambda. Generally this |
64 + | # is the maximum baseline at the highest frequency. |
65 + | light_speed = qa.getvalue(qa.convert(qa.constants('c'), 'm/s'))[0] |
66 + | max_bl = 0.0 |
67 + | for spwid in self.spwids: |
68 + | real_spwid = self.observing_run.virtual2real_spw_id(spwid, ms) |
69 + | spw = ms.get_spectral_window(real_spwid) |
70 + | mean_freq_Hz = spw.mean_frequency.to_units(measures.FrequencyUnits.HERTZ) |
71 + | mean_wave_m = light_speed / float(mean_freq_Hz) # in meter |
72 + | # Get max baseline |
73 + | job = casa_tasks.visstat(vis=vis, spw=str(spwid), axis='uvrange', useflags=True) |
74 + | uv_stat = job.execute(dry_run=False) # returns stat in meter |
75 + | max_bl = max(max_bl, uv_stat['DATA_DESC_ID=%s' % spwid]['max'] / mean_wave_m) |
76 + | |
77 + | # Define bin for lowest 5% of baselines (in wavelength units) |
78 + | uvll = 0.0 |
79 + | uvul = 0.05 * max_bl |
80 + | |
81 + | uvrange = '{:0.1f}~{:0.1f}klambda'.format(uvll / 1000.0, uvul / 1000.0) |
82 + | mean_SBL = get_mean_amplitude(vis=vis, uvrange=uvrange, spw=spw_str) |
83 + | |
84 + | # Range for 50-55% bin |
85 + | uvll = 0.5 * max_bl |
86 + | uvul = 0.5 * max_bl + 0.05 * max_bl |
87 + | uvrange = '{:0.1f}~{:0.1f}klambda'.format(uvll / 1000.0, uvul / 1000.0) |
88 + | mean_MBL = get_mean_amplitude(vis=vis, uvrange=uvrange, spw=spw_str) |
89 + | |
90 + | # Compare amplitudes and decide on return value |
91 + | meanratio = mean_SBL / mean_MBL |
92 + | |
93 + | if meanratio > 2.0: |
94 + | LOG.info('Selecting uvrange>{:0.1f}klambda to avoid very extended emission.'.format(0.05 * max_bl / 1000.0)) |
95 + | return ">{:0.1f}klambda".format(0.05 * max_bl / 1000.0) |
96 + | else: |
97 + | return None |
98 + | |
26 99 | def pblimits(self, pb): |
27 100 | """ |
28 101 | PB gain level at which to cut off normalizations (tclean parameter). |
29 102 | |
30 103 | See PIPE-674 and CASR-543 |
31 104 | """ |
32 105 | pblimit_image = -0.1 |
33 106 | pblimit_cleanmask = -0.1 # not used at the moment, negative values are untested (see PIPE-674) |
34 107 | |
35 108 | return pblimit_image, pblimit_cleanmask |