Commits

Takeshi Nakazato authored 21ae7faa13b Merge
Merge branch 'main' into PIPE-2278-hsd_skycal-error-indexerror-index-0-is-out-of-bounds-for-axis-0-with-size-0
No tags

pipeline/extern/almarenorm.py

Modified
95 95 from casatools import measures as metool
96 96 from casatools import calibrater as cbtool
97 97
98 98 mytb=tbtool()
99 99 myms=mstool()
100 100 mycb=cbtool()
101 101
102 102
103 103 # Interface function for use by ALMA Pipeline renormalization task (PIPE-1687).
104 104 def alma_renorm(vis: str, spw: List[int], create_cal_table: bool, threshold: Union[None, float], excludechan: Dict,
105 - correct_atm: bool, atm_auto_exclude: bool, bwthreshspw: Dict, caltable: str) -> tuple: ## PIPE-1469 for bwthreshspw
105 + atm_auto_exclude: bool, bwthreshspw: Dict, caltable: str) -> tuple: ## PIPE-1469 for bwthreshspw
106 106 """
107 107 Interface function for ALMA Pipeline: this runs the ALMA renormalization
108 108 heuristic for input vis, and returns statistics and metadata required by
109 109 the Pipeline renormalization task.
110 110
111 111 Args:
112 112 vis: path to the measurement set to renormalize.
113 113 spw: Spectral Window IDs to process.
114 114 create_cal_table: save renorm corrections to a calibration table.
115 115 threshold: the threshold above which the scaling must exceed for the
155 155 atmExcludeCmd = {}
156 156 calTableCreated = False
157 157
158 158 # Store "all TDM" information before trying the renormalization so
159 159 # that the weblog is rendered properly.
160 160 alltdm = rn.tdm_only
161 161
162 162 # Only attempt the correction if the data are not all TDM.
163 163 if not alltdm:
164 164 # Run the renormalisation.
165 - rn.renormalize(createCalTable=create_cal_table, fillTable=True, docorrThresh=threshold, correctATM=correct_atm, spws=spw, excludechan=excludechan,
166 - atmAutoExclude=atm_auto_exclude, bwthreshspw=bwthreshspw, datacolumn='DATA') # add bwthreshspw
165 + rn.renormalize(createCalTable=create_cal_table, fillTable=create_cal_table, docorrThresh=threshold, spws=spw, excludechan=excludechan,
166 + atmAutoExclude=atm_auto_exclude, bwthreshspw=bwthreshspw, datacolumn='DATA')
167 167
168 168 # Create spectra plots.
169 169 rn.plotSpectra(includeSummary=False)
170 170
171 171 # Check if the renormalisation has been saved to a cal table.
172 172 calTableCreated = os.path.exists(rn.rntable)
173 173
174 174 # Get stats (dictionary) indexed by source, spw
175 175 stats = rn.rnpipestats
176 176
1188 1188
1189 1189 # Need to keep track of the number of times we write to the renorm table.
1190 1190 rntb_iter = 0
1191 1191
1192 1192 if fillTable:
1193 1193 casalog.post('All targets, scans, spws will be filled in resulting calibration table. This may greatly enlarge the size of the final table!')
1194 1194 casalog.post(f'TDM windows {self.tdmspws} will also be included in {self.rntable}.')
1195 1195
1196 1196 else:
1197 1197 casalog.post('No corrections will be saved (createCalTable=False)!')
1198 + fillTable=False
1198 1199
1199 1200 # Check if docorrThresh is set correctly
1200 1201 if docorrThresh is not None:
1201 1202 if type(docorrThresh) is not float:
1202 1203 casalog.post('Correction of CORRECTED_DATA requested, but docorrThresh is set incorrectly! Cannot procede.')
1203 1204 casalog.post(' set to None for automatic thresholding during apply, or input a float to use')
1204 1205 casalog.post('*** Terminating renormalization run ***', 'INFO', 'ReNormalize')
1205 1206 raise TypeError('Correction of CORRECTED_DATA requested, but docorrThresh is set incorrectly. Use None or float.')
1206 1207 if docorrThresh > 1.5:
1207 1208 casalog.post('WARNING: Correction of CORRECTED_DATA requested, but docorrThresh is set very high')
2153 2154 if self.atmWarning[str(fldname)][str(ispw)]:
2154 2155 if atmAutoExclude:
2155 2156 casalog.post('Atmospheric features above the threshold have been mitigated by atmAutoExclude.')
2156 2157 casalog.post('Equivalent manual call: '+exclude_cmd)
2157 2158 else:
2158 2159 casalog.post('ATM features may be falsely triggering renorm!')
2159 2160 casalog.post('Suggested channel exclusion: '+exclude_cmd)
2160 2161 elif exclude_cmd != "No flagging suggested.":
2161 2162 casalog.post('Atmospheric features mitigated.')
2162 2163 casalog.post('Equivalent manual call: '+exclude_cmd)
2163 -
2164 - # If fillTable is selected as an option, then we need to also fill in
2165 - # all the science TDM windows as well (if any) for this target.
2166 - if len(self.tdmspws) > 0 and fillTable:
2167 - casalog.post(f'\nFilling in TDM spws in calibration table {self.rntable} with value of 1.0')
2168 - for ispw in self.tdmspws:
2169 - #casalog.post(f' Writing filler solutions for spw {ispw}')
2170 - # Not all targets are in all scans, we need to iterate over only those scans containing the target
2171 - target_scans = np.intersect1d(self.msmeta.scansforintent('*TARGET*'), self.msmeta.scansforfield(target))
2172 -
2173 - # Make an additional cut to catch only those scans which contain the current spw (usually only relevant
2174 - # for spectral scan datasets)
2175 - target_scans = np.intersect1d(target_scans, self.msmeta.scansforspw(ispw))
2176 -
2177 - # If user input list of scans to use, cross check those with the list of all scans on targets to make
2178 - # sure it's necessary to perform this loop.
2179 - target_scans = np.intersect1d(target_scans, targscans)
2180 -
2181 - # global list of scans with the current spw
2182 - spwscans=list(self.msmeta.scansforspw(ispw))
2183 2164
2184 - for iscan in target_scans:
2185 - #casalog.post(f' Scan {iscan}')
2186 - # get the fields to process in this scan - i.e. mosaics have many fields per scan
2187 - Tarfld = list(self.msmeta.fieldsforscan(iscan))
2165 + if createCalTable:
2166 + # If fillTable is selected as an option, then we need to also fill in
2167 + # all the science TDM windows as well (if any) for this target.
2168 + if len(self.tdmspws) > 0 and fillTable:
2169 + casalog.post(f'\nFilling in TDM spws in calibration table {self.rntable} with value of 1.0')
2170 + for ispw in self.tdmspws:
2171 + #casalog.post(f' Writing filler solutions for spw {ispw}')
2172 + # Not all targets are in all scans, we need to iterate over only those scans containing the target
2173 + target_scans = np.intersect1d(self.msmeta.scansforintent('*TARGET*'), self.msmeta.scansforfield(target))
2188 2174
2189 - for ifld in Tarfld:
2190 - # Skip over target scans that don't have the current spw - can happen for spectral scans?
2191 - if spwscans.count(iscan) == 0:
2192 - continue
2193 - casalog.post(f' Writing filler solutions for spw {ispw}, scan {iscan}, field {ifld}')
2194 - N_fill = np.ones((self.num_corrs, self.msmeta.nchan(ispw), self.nAnt))
2195 - self.writeCalTable(ispw, ifld, iscan, N_fill, rntb_iter)
2196 - rntb_iter += 1
2175 + # Make an additional cut to catch only those scans which contain the current spw (usually only relevant
2176 + # for spectral scan datasets)
2177 + target_scans = np.intersect1d(target_scans, self.msmeta.scansforspw(ispw))
2178 +
2179 + # If user input list of scans to use, cross check those with the list of all scans on targets to make
2180 + # sure it's necessary to perform this loop.
2181 + target_scans = np.intersect1d(target_scans, targscans)
2182 +
2183 + # global list of scans with the current spw
2184 + spwscans=list(self.msmeta.scansforspw(ispw))
2185 +
2186 + for iscan in target_scans:
2187 + #casalog.post(f' Scan {iscan}')
2188 + # get the fields to process in this scan - i.e. mosaics have many fields per scan
2189 + Tarfld = list(self.msmeta.fieldsforscan(iscan))
2190 +
2191 + for ifld in Tarfld:
2192 + # Skip over target scans that don't have the current spw - can happen for spectral scans?
2193 + if spwscans.count(iscan) == 0:
2194 + continue
2195 + casalog.post(f' Writing filler solutions for spw {ispw}, scan {iscan}, field {ifld}')
2196 + N_fill = np.ones((self.num_corrs, self.msmeta.nchan(ispw), self.nAnt))
2197 + self.writeCalTable(ispw, ifld, iscan, N_fill, rntb_iter)
2198 + rntb_iter += 1
2197 2199
2198 2200
2199 2201
2200 2202 # PIPE 1168 (3)
2201 2203 # Loops through the scalingValue dict and populates the pipeline needed dictionary
2202 2204 self.rnpipestats = {}
2203 2205 target_field_ids = self.msmeta.fieldsforintent('*TARGET*')
2204 2206 target_fields = np.unique(self.msmeta.namesforfields(target_field_ids))
2205 2207 for trg in self.rnstats['N'].keys(): #target_fields:
2206 2208 self.rnpipestats[trg] = {}
2283 2285 A book-keeping value that keeps track of how many times we
2284 2286 have written to the table to avoid over-writing data.
2285 2287
2286 2288 Output:
2287 2289 None, the calibration table initialized by self.renormalize() is
2288 2290 updated with the contents of input N and all related metadata.
2289 2291 """
2290 2292 # We need the number of antennas from the array shape as this will be
2291 2293 # the number of rows we'll be adding to the table currently.
2292 2294 nPol, nChan, nAnt = N.shape
2295 +
2296 + # For the single pol case, fill the other pol with ones so that the
2297 + # Pipeline can handle it.
2298 + if nPol == 1:
2299 + N_fill = np.ones(N.shape)
2300 + N = np.concatenate((N,N_fill))
2293 2301
2294 2302 # Antenna based corrections so we need 1 row for each antenna
2295 2303 self.rntb.addrows(nrow = nAnt)
2296 2304
2297 2305 # Placing the renorm values themselves into the table
2298 2306 self.rntb.putcol('FPARAM', N, startrow=iter_n*nAnt, nrow=nAnt)
2299 2307
2300 2308 # Following the convention of the Tsys table, filling all values with 0.1 for the error.
2301 2309 self.rntb.putcol('PARAMERR', N*0.0+0.1, startrow=iter_n*nAnt, nrow=nAnt)
2302 2310

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut