# THIS SCRIPT PROVIDES: # BACKGROUND SUBTRACTION OF OFF-AXIS DETECTORS SELECTED USING 3SIGMA-CLIPPING # - PNG file with two plots: spectra before and after subtraction and off-axis # detectors included in the median that was subtracted # - Fits file with the spectra after the subtraction # For processing a complete list of observation Ids # Parameters for complete CSV table obsidS = String1d() targetS = String1d() meansSLWS = Float1d() meansSSWS = Float1d() meansSSW2S = Float1d() gapS = Float1d() relGapS = Float1d() radecoffS = Float1d() # Import list of obs ids and selection tags pagenumber = '7' obsidList = asciiTableReader(file='/home/salvarez/Desktop/EXCEL/newPage%sobsids.txt'%(pagenumber), tableType='CSV') obsidList = asciiTableReader(file='/home/salvarez/Desktop/EXCEL/lowresObsid.csv', tableType='CSV') columnSelections = obsidList['c0'].data[1:] columnObsIds = obsidList['c1'].data[1:] # for val,data in enumerate(columnSelections): if data == "Good": obsid = columnObsIds[val] print obsid # Run the good script [target,radecoff,meanSLW,meanSSW,meanSSW2,gap,relGap,res] = good(obsid) # Data for CSV table obsidS = obsidS.append(str(obsid)) targetS = targetS.append(target) radecoffS = radecoffS.append(radecoff) meansSLWS = meansSLWS.append(meanSLW) meansSSWS = meansSSWS.append(meanSSW) meansSSW2S = meansSSW2S.append(meanSSW2) gapS = gapS.append(gap) relGapS = relGapS.append(relGap) # # Creating table tds = TableDataset() tds.addColumn("ObservationId", Column(obsidS)) tds.addColumn("Object", Column(targetS)) tds.addColumn("MedianSLW", Column(meansSLWS)) tds.addColumn("MedianSSW", Column(meansSSWS)) tds.addColumn("MedianSSW2", Column(meansSSW2S)) tds.addColumn("Gap", Column(gapS)) tds.addColumn("RelGap", Column(relGapS)) tds.addColumn("raDecOffset", Column(radecoffS)) from herschel.share.unit import Scalar from herschel.share.unit import Angle from herschel.share.unit import FluxDensity tds["MedianSLW"].unit = Scalar.PERCENT tds["MedianSSW"].unit = Scalar.PERCENT tds["MedianSSW2"].unit = Scalar.PERCENT tds["Gap"].unit = FluxDensity.JANSKYS tds["RelGap"].unit = Scalar.PERCENT tds["raDecOffset"].unit = Angle.SECONDS_ARC #asciiTableWriter(table=tds,file = '/home/salvarez/Desktop/HPDP/BKG/GOOD/page%s_median3goods.csv'%(pagenumber)) # #-------------------------------------------------- def good(obsid): obs = getObservation(obsid,useHsa=True) # res = obs.meta['commandedResolution'].getValue() if res =='MR': res = 'LR' if res =='CR': res = 'HR' if res =='H+LR': res = 'HR' # print res target = obs.refs["level2"].product.refs["%s_spectrum_point"%(res)].product.meta['object'].value radecoff = obs.refs["level2"].product.refs["%s_spectrum_point"%(res)].product.meta['raDecOffset'].value if res == 'LR': # Ratios before background subtraction [meanSLW,meanSSW,meanSSW2] = ratiosFirstRings2(obs.refs["level2"].product.refs["%s_spectrum_point"%(res)].product) # # Subtract Background pointSpec = backgroundSubtract(obs,1,0,21,'%s'%(res)) simpleFitsWriter(product=pointSpec, file='/home/salvarez/Desktop/HPDP/BKG/GOOD/%s_BackSubt_median3sigma.fits'%(obsid), warn=False) # # Calculate gap [gap, relGap] = calcGap(pointSpec,res) # else: meanSLW =0 meanSSW = 0 meanSSW2 = 0 gap = 0 relGap = 0 return target,radecoff,meanSLW,meanSSW,meanSSW2,gap,relGap,res # ------------------------------------------ def backgroundSubtract(obs,selectOffs,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' offAxisDets = {'SSW':['SSWB2', 'SSWB3', 'SSWB4', 'SSWC2', 'SSWC5', 'SSWD2', \ 'SSWD6', 'SSWE5', 'SSWF3','SSWE2','SSWF2'], \ '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)) # 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 MEDIAN 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 = [] plotSmooth = 1 plotFinal = PlotXY() plotFinal.plotSize = (6,2) topLeft = SubPlot(SubPlotGridConstraints(0,1)) noSkipDetect = [] for det in centreDetectors: medianOffSmList = [] myArray = [] for val,offDet in enumerate(offAxisDets[det[:3]]): if det == 'SLWC3': # For HR: #myArray.append(offAxisSm['0000'][offDet].fluxColumn.data[1678]) # For LR: myArray.append(offAxisSm['0000'][offDet].fluxColumn.data[67]) elif det== 'SSWD4': # For HR: #myArray.append(offAxisSm['0000'][offDet].fluxColumn.data[1521]) # For LR: myArray.append(offAxisSm['0000'][offDet].fluxColumn.data[61]) med = MEDIAN(myArray) mad = MedianAbsoluteDeviation(med)(myArray) for val,offDet in enumerate(offAxisDets[det[:3]]): #in case any detectors are missing try: if ABS(myArray[val]-med)<= 3*mad: medianOffSmList.append(offAxisSm['0000'][offDet].flux) detectorsUsed.append(offDet) noSkipDetect.append(offDet) continue except: pass if det == 'SLWC3': medianOffSmSLWC3 = findMedianSpectrum(medianOffSmList) layerMedianSLWC3 = LayerXY(pointOffAxisSubSm['0000'][det].wave,medianOffSmSLWC3, line=Style.DASHED, \ color=java.awt.Color.BLACK, name='Median', stroke=2) # Subtract the smoothed median off-axis sum from the input data pointOffAxisSubSm['0000'][det].flux = pointSpec['0000'][det].flux - medianOffSmSLWC3 elif det == 'SSWD4': medianOffSmSSWD4 = findMedianSpectrum(medianOffSmList) layerMedianSSWD4 = LayerXY(pointOffAxisSubSm['0000'][det].wave,medianOffSmSSWD4, line=Style.DASHED, \ color=java.awt.Color.BLACK, name='Median', stroke=2) # Subtract the smoothed median off-axis sum from the input data pointOffAxisSubSm['0000'][det].flux = pointSpec['0000'][det].flux - medianOffSmSSWD4 color = [226,46,46,114,21,21,242,61,245,123,127,202,81,141,245,81,221,245,55,186,156,55,122,25,124,149,58,229,42,9,225,186,72,149,237,48,202,211,29,245,146,27,164,164,163,185,179,179,0,204,0] i=0 j=1 k=2 for offDet in noSkipDetect: x=color[i] y=color[j] z=color[k] plotDetectors(topLeft,offDet,x,y,z,pointSpec,det,offAxisSm) i=i+3 j=j+3 k=k+3 topLeft.addLayer(layerMedianSLWC3) topLeft.addLayer(layerMedianSSWD4) layerMedianSSWD4.inLegend = 0 ## PLOT THE RESULTS ##################### bottomLeft = SubPlot(SubPlotGridConstraints(0,0)) inLegend = 1 layer1 = LayerXY(Double1d(0), Double1d(0)) for det in centreDetectors: # The original data direct from the observation (no subtraction) layer1 = LayerXY(pointSpec['0000'][det].wave, pointSpec['0000'][det].flux, \ color=java.awt.Color.RED, name='Original data', inLegend=inLegend) bottomLeft.addLayer(layer1) if plotSmooth: # Averaged smoothed off-axis detector subtracted data bottomLeft.addLayer(LayerXY(pointSpec['0000'][det].wave, pointOffAxisSubSm['0000'][det].flux, \ color=java.awt.Color.GREEN, name='Off-axis subtracted', \ inLegend=inLegend)) inLegend = 0 for subP in [bottomLeft,topLeft]: # Remove all axes labels to start with (add needed ones later) subP.baseLayerXY.xaxis.tick.labelVisible = 0 subP.baseLayerXY.xaxis.title.visible = 0 subP.baseLayerXY.yaxis.tick.labelVisible = 0 subP.baseLayerXY.yaxis.title.visible = 0 # Set axis ranges to be the same for all sub-plots subP.baseLayerXY.xaxis.range = [400.0, 1600.0] # subP.baseLayerXY.yaxis.range = [-0.1, 1.1] # Set the ticks for the xaxis to be at nice intervals # subP.baseLayerXY.xaxis.tick.setFixedValues(Double1d([500, 1000, 1500]), \ # Double1d(range(400, 1600, 100))) plotFinal.addSubPlot(subP) # Set the gap between the sub-plots to be zero (i.e. plots touching each other) plotFinal.gridLayout.setGap(0, 0) # Set the titles for the left hand axes of left hand plots # and make the tick labels visible on these axes topLeft.baseLayerXY.yaxis.title.visible = 1 topLeft.baseLayerXY.yaxis.titleText = "Flux Density" topLeft.baseLayerXY.yaxis.tick.labelVisible = 1 bottomLeft.baseLayerXY.yaxis.title.visible = 1 bottomLeft.baseLayerXY.yaxis.titleText = "Flux Density" bottomLeft.baseLayerXY.yaxis.tick.labelVisible = 1 # Make the tick labels visible on the lower xaxis bottomLeft.baseLayerXY.xaxis.tick.labelVisible = 0 topLeft.baseLayerXY.xaxis.tick.labelVisible = 1 topLeft.baseLayerXY.xaxis.title.visible = 1 topLeft.baseLayerXY.xaxis.titleText = "Frequency (GHz)" # Use the upper xaxis for a title for each column # topLeft.baseLayerXY.xaxis.getAuxAxis(0).titleText = '%i OD%s'%(pointSpec.meta['obsid'].value,\ # pointSpec.meta['odNumber'].value) topLeft.baseLayerXY.xaxis.getAuxAxis(0).title.visible = 0 # Use the plot subtitle to contain a single label for all 3 x-axes # plotFinal.subtitle.text = "Frequency (GHz)" # plotFinal.subtitle.fontSize = plotFinal.title.fontSize # from herschel.ia.gui.plot import PlotTitle # plotFinal.subtitle.position = PlotTitle.BOTTOMCENTER # plot1.baseLayerXY.xaxis.titleText = '%s [%s]'%(pointSpec['0000'][det].waveDescription, \ # pointSpec['0000'][det].waveUnit.dialogName) # plot1.baseLayerXY.yaxis.titleText = '%s [%s]'%(pointSpec['0000'][det].fluxDescription, \ # pointSpec['0000'][det].fluxUnit.dialogName) # plotFinal.addSubPlot(plot1) # # plotAll.baseLayerXY.xaxis.titleText ='%s [%s]'%(pointSpec['0000'][det].waveDescription, \ # pointSpec['0000'][det].waveUnit.dialogName) # plotAll.baseLayerXY.xaxis.fontSize = 10 # plotAll.baseLayerXY.yaxis.titleText = '%s [%s]'%(pointSpec['0000'][det].fluxDescription, \ # pointSpec['0000'][det].fluxUnit.dialogName) # plotAll.baseLayerXY.yaxis.fontSize = 10 # plotFinal.addSubPlot(plotAll) plotFinal.titleText = '%i OD%s'%(pointSpec.meta['obsid'].value,\ pointSpec.meta['odNumber'].value) plotFinal.subtitleText = '%s, %s reps'%(pointSpec.meta['object'].value, pointSpec.meta['numRepetitions'].value) plotFinal.subtitle.fontSize = 10 plotFinal.legend.fontSize = 10 plotFinal.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 plotFinal.saveAsPNG('/home/salvarez/Desktop/changes/%s_BackSubt_median3sigma.png'%(pointSpec.meta['obsid'].value)) # plotFinal.close() ## RETURN THE RESULTS return pointOffAxisSubSm # -------------------------------------------- # Function for plotting the off-axis detectors. def plotDetectors(plotAll,offDet,x,y,z,pointSpec,det,offAxisSm): plotAll.addLayer(LayerXY(offAxisSm['0000'][offDet].wave, offAxisSm['0000'][offDet].flux,\ name=offDet, color=java.awt.Color(x,y,z), stroke=1)) return plotAll # -------------------------------------------- # Function for taking the median across a Jython List def findMedianSpectrum(list): medianOut = Double1d() # Take the median across the python list for each # spectral element for i in range (len(list[0])): temp = Double1d() for listBit in list: temp.append(listBit[i]) medianOut.append(MEDIAN(temp)) return medianOut # ---------------------------------------------- # Calculate gap def calcGap(pointSpec,res): spectrumSLW = pointSpec['0000']['SLWC3'].flux spectrumSSW = pointSpec['0000']['SSWD4'].flux freqSLW = pointSpec['0000']['SLWC3'].wave freqSSW = pointSpec['0000']['SSWD4'].wave if res=='HR': # Long wavelength indListSLW = freqSLW.where(freqSLW>944) indSLW = 1658 overLapSpectSLW = spectrumSLW[indSLW:len(freqSLW)] # Short wavelength indListSSW = freqSSW.where(freqSSW<1018) indSSW = 247 overLapSpectSSW = spectrumSSW[0:indSSW] else: # Long wavelength indListSLW = freqSLW.where(freqSLW>944) indSLW = 66 overLapSpectSLW = spectrumSLW[indSLW:len(freqSLW)] # Short wavelength indListSSW = freqSSW.where(freqSSW<1018) indSSW = 10 overLapSpectSSW = spectrumSSW[0:indSSW] # Calculate gap gap = MEDIAN((overLapSpectSLW-overLapSpectSSW)) # Calculate gap percentage relative to long wavelength spectra relGap = MEDIAN((overLapSpectSLW-overLapSpectSSW)/overLapSpectSLW)*100 return gap, relGap # -------------------------------------------------- def ratiosFirstRings2(pointSpec): # ring1SLW = ['SLWB2','SLWB3','SLWC2','SLWC4','SLWD2','SLWD3'] ring1SSW = ['SSWD3','SSWC3','SSWC4','SSWE4','SSWE3'] #D5 is a dead bolometer # ring2SSW = ['SSWB2', 'SSWB3', 'SSWB4', 'SSWC2', 'SSWC5', 'SSWD2','SSWD6', 'SSWE5', 'SSWF3','SSWE2','SSWF2'] # offAxisSm = pointSpec # central detectors data SLW freq = offAxisSm['0000']['SLWC3'].wave flux = offAxisSm['0000']['SLWC3'].flux ratiosSLW = Double2d(len(ring1SLW),len(flux)) medianratiosSLW = Double1d() k=0 # For every detector in the first ring of SLW: for detector in ring1SLW: ratios = Double1d() fluxdetector = offAxisSm['0000'][detector].flux for i in range(0,len(fluxdetector)): ratios.append((fluxdetector[i]/flux[i])*100) ratiosSLW[k,:] = ratios k += 1 # Median of the different ratios for each detector # I take out the dependency on frequency medianratiosSLW = medianratiosSLW.append(MEDIAN(ratios)) ratio = ratios #This is the average for all detectors averageSLW = MEDIAN(medianratiosSLW) # central detectors data SSW freq = offAxisSm['0000']['SSWD4'].wave flux = offAxisSm['0000']['SSWD4'].flux ratiosSSW = Double2d(len(ring1SSW),len(flux)) medianratiosSSW = Double1d() k=0 for detector in ring1SSW: ratios = Double1d() fluxdetector = offAxisSm['0000'][detector].flux for i in range(0,len(fluxdetector)): ratios.append((fluxdetector[i]/flux[i])*100) ratiosSSW[k,:] = ratios k += 1 # Median of the different ratios for each detector # I get an array dependent on frequency medianratiosSSW = medianratiosSSW.append(MEDIAN(ratios)) ratio = ratios #This is the average in frequency by taking the Median averageSSW = MEDIAN(medianratiosSSW) #SECOND RING OF SSW ratiosSSW2 = Double2d(len(ring2SSW),len(flux)) medianratiosSSW2 = Double1d() k=0 for detector in ring2SSW: ratios = Double1d() fluxdetector = offAxisSm['0000'][detector].flux for i in range(0,len(fluxdetector)): ratios.append((fluxdetector[i]/flux[i])*100) ratiosSSW2[k,:] = ratios k += 1 # Median of the different ratios for each detector # I take out the dependency on frequency medianratiosSSW2 = medianratiosSSW2.append(MEDIAN(ratios)) ratio = ratios #This is the average for all detectors averageSSW2 = MEDIAN(medianratiosSSW2) return averageSLW, averageSSW, averageSSW2