## BACKGROUND SUBTRACTION ########################## # Function for Background Subtraction # -------------------------------------------- def backgroundSubtract(obs,selectOffs,subtractDark,darkObsid,kernelWidth,res): # Parameters: # 1. Interactive selection of off-axis detectors? selectOffs = 0/1 # 2. Subtract a dark sky observation? subtractDark = 0/1 # 3. Set Kernel width (GHz) for smoothing: Example: kernelWidth = 21 # 4. Spectral resolution of the input data: res = 'HR','LR','MR' productName = 'spectrum_point' # 9. Define the off-axis detectors in a PyDictionary # This list is used for the off-axis subtraction. It can be edited to remove # detectors or add other unvignetted off-axis detectors. # SSWE2 and SSWF2 are excluded from the SSW list as habitual outliers offAxisDets = {'SSW':['SSWB2', 'SSWB3', 'SSWB4', 'SSWC2', 'SSWC5', 'SSWD2', \ 'SSWD6', 'SSWE5', 'SSWF3'], \ 'SLW':['SLWD3', 'SLWC4', 'SLWB3', 'SLWB2', 'SLWC2', 'SLWD2']} ## READ THE DATA ################# # Update the format of the observation context to match HIPE v13 obs = updateObsContext(obs) # Get the point source calibrated data pointSpec = obs.level2.getProduct('%s_%s'%(res, productName)) # Get the dark sky if requested if subtractDark == 1: darkObs = getObservation(darkObsid, useHsa=True) darkObs = updateObsContext(darkObs) darkSkyPointSpec = darkObs.level2.getProduct('%s_%s'%(res, productName)) # Define the centre detectors centreDetectors = ['SLWC3', 'SSWD4'] # Check for any missing off-axis detectors offAxisDetsCheck = {'SSW':[],'SLW':[]} for key in offAxisDets.keys(): for det in offAxisDets[key]: if pointSpec['0000'][det]: offAxisDetsCheck[key].append(det) offAxisDets = offAxisDetsCheck ## SUBTRACT THE DARK SKY ######################### # A dark sky observation is made on every observational day. This can be # downloaded from the archive and simply subtracted. if subtractDark == 1: try: darkSubtracted = pointSpec.copy() - darkSkyPointSpec except: darkSubtracted = pointSpec.copy() for det in pointSpec[0].channelNames: try: darkSubtracted[0][det] = pointSpec[0][det].copy() \ - darkSkyPointSpec[0][det] except: pass ## SUBTRACT THE MEAN CO-ALIGNED OFF-AXIS DETECTORS #################################################### # The off-axis detectors from the same observation can be used to correct # the spectral shape, assuming there is no structured background emission. # Ideally, the off-axis detectors should be examined by hand and the best # ones selected for averaging. # Note that no error propagation is performed for the smoothed subtraction. pointOffAxisSubSm = pointSpec.copy() # Smooth the off-axis detectors offAxisSm = smooth(ds=pointSpec, \ filter='Gaussian', width=kernelWidth, unit='GHz', \ edge='SHRINK_WIDTH', variant='wave-flag-weight', \ overwrite=0) # Create a variable to keep track of the off-axis detectors used in the subtraction detectorsUsed = [] # Select which off-axis detectors to include in the mean sum to subtract # from the central detectors. # To skip the visual inspection, set selectOffs = 0 above and all listed # off-axis detectors will be included in the median sum. plotSmooth = 1 if selectOffs: for det in centreDetectors: count = len(offAxisDets[det[:3]]) averageOffSm = 0 plotAll = plotDetectors(offAxisDets[det[:3]], offAxisSm, highlight='0') check = [''] check = raw_input('Exclude any detectors?: y/n') plotAll.close() if check == None: # print 'Cancelling off-axis selection and subtraction ....' plotSmooth = 0 break if check == 'y': count = 0 for offDet in offAxisDets[det[:3]]: plotAll = plotDetectors(offAxisDets[det[:3]], offAxisSm, highlight=offDet) ok = [''] ok = raw_input('Exclude %s from the mean sum? Type "y" to exclude or any other key to keep: '%offDet) if ok == None: plotAll.close() # print 'Cancelling off-axis selection and subtraction ....' plotSmooth = 0 break elif ok != 'y': # print 'Including %s'%offDet averageOffSm = averageOffSm + offAxisSm['0000'][offDet] detectorsUsed.append(offDet) count += 1 else: # print 'Excluding %s'%offDet continue plotAll.close() if count == 0: plotSmooth = 0 break else: for offDet in offAxisDets[det[:3]]: averageOffSm = averageOffSm + offAxisSm['0000'][offDet] detectorsUsed.append(offDet) # Find the mean of the off-axis sum averageOffSm = averageOffSm/float(count) # Subtract the smoothed off-axis sum from the input data pointOffAxisSubSm['0000'][det] = pointSpec['0000'][det] - averageOffSm # Else, run without visual inspection and take the MEDIAN of all six off-axis # detectors in the list for each centre detector. # Note: the median is taken to replace the visual inspection to remove outlying # detectors else: for det in centreDetectors: medianOffSmList = [] for offDet in offAxisDets[det[:3]]: #in case any detectors are missing try: medianOffSmList.append(offAxisSm['0000'][offDet].flux) detectorsUsed.append(offDet) except: pass medianOffSm = findMedianSpectrum(medianOffSmList) # Subtract the smoothed median off-axis sum from the input data pointOffAxisSubSm['0000'][det].flux = pointSpec['0000'][det].flux - medianOffSm ## PLOT THE RESULTS ##################### plot1 = PlotXY() inLegend = 1 for det in centreDetectors: # The original data direct from the observation (no subtraction) plot1.addLayer(LayerXY(pointSpec['0000'][det].wave, pointSpec['0000'][det].flux, \ color=java.awt.Color.RED, name='Original data', inLegend=inLegend)) if subtractDark: # The dark subtracted data plot1.addLayer(LayerXY(pointSpec['0000'][det].wave, darkSubtracted['0000'][det].flux, \ color=java.awt.Color.BLACK, name='Dark subtracted', inLegend=inLegend)) if plotSmooth: # Averaged smoothed off-axis detector subtracted data plot1.addLayer(LayerXY(pointSpec['0000'][det].wave, pointOffAxisSubSm['0000'][det].flux, \ color=java.awt.Color.GREEN, name='Off-axis subtracted', \ inLegend=inLegend)) inLegend = 0 plot1.xtitle = '%s [%s]'%(pointSpec['0000'][det].waveDescription, \ pointSpec['0000'][det].waveUnit.dialogName) plot1.ytitle = '%s [%s]'%(pointSpec['0000'][det].fluxDescription, \ pointSpec['0000'][det].fluxUnit.dialogName) plot1.xrange = [MIN(pointSpec['0000']['SLWC3'].wave), MAX(pointSpec['0000']['SSWD4'].wave)] plot1.titleText = '%i OD%s'%(pointSpec.meta['obsid'].value,\ pointSpec.meta['odNumber'].value) plot1.subtitleText = '%s, %s reps'%(pointSpec.meta['object'].value, pointSpec.meta['numRepetitions'].value) plot1.subtitle.fontSize = 10 plot1.legend.fontSize = 10 plot1.legend.visible = 1 # Print the off-axis detectors that were included in the off-axis subtraction # print '\n Off-axis detectors used for subtraction:\n', detectorsUsed plot1.saveAsPNG('/home/salvarez/backgroundSubtractionPlot_%s.png'%(pointSpec.meta['obsid'].value)) # plot1.close() ## RETURN THE RESULTS ###################### return pointOffAxisSubSm