################################################################ # # Title: OT2_Workshop_SPHOT_Photometry # # Scope: this script presents different ways of doing photometry # (in this case, of a point source) with SPIRE # # Luca Conversi - 21/02/2012 - ver. 1.0 # ################################################################ ######################### Demo set-up ########################## # # For this demo, you first need to download the bendoSourceFit.py script from here: # # http://herschel.esac.esa.int/twiki/pub/Public/SpireCalibrationWeb/bendoSourceFit_v0_9.py.txt # # Then you need to load into your HIPE session: # - change the extension, removing the .txt # - open into HIPE # - click on "Run all" (the double green arrows button) # ################################################################ #################### Timeline fitter set-up #################### # # First, we load an observation of Gamma Dra obs = getObservation(1342195871, useHsa=1) # Then we remove the AC component from the Level 1 scans = baselineRemovalMedian(obs.level1) # And finally we 'give' the baseline removed Level 1 scans to the fitter fitter = bendoSourceFit(scans) # ################################################################ ############# Timeline fitter on a specific source ############# # # Display PSW map on a pop-up window Display(obs.level2.refs['PSW'].product) # Now right click on the central source and from the context menu select "Get coordinates" # In the console window the pixel and world coordinates of the target will be shown # Now let's run the timeline fitter (bendoSourceFit) on the selected source raCentre = displayWorldCoords[0] decCentre = displayWorldCoords[1] # Input parameters needed by the fitter. # rPeak is the radius of the region that will include the peak of the source, in arcsec # Suggested values are 22, 30, and 42 for PSW, PMW, and PLW, respectively. arrays = ['PSW', 'PMW', 'PLW'] rPeaks = [22, 30, 42] # Run the fitter on the 3 arrays and print results for j in range(len(arrays)): param = fitter.fit(arrays[j], raCentre, decCentre, rPeaks[j]) chi = fitter.getChiSquared() print 'Source flux (in mJy) for %s array is: %f +/- %f'%(arrays[j], param[0] * 1.e3, param[3] * 1.e3) # ################################################################ ############# Timeline fitter of extracted sources ############# # # Available source extractor algorithms: # - sourceExtractorSussextractor: it does PSF fitting on map # - sourceExtractorDaophot: it does aperture photometry, so you need to correct the results. # However, the aperture correction factors are still not available, so # the user must calculate them # # The beam profiles are available in the calibration tree: # beamPSW = obs.calibration.phot.getBeamProf('PSW') # beamPMW = obs.calibration.phot.getBeamProf('PMW') # beamPLW = obs.calibration.phot.getBeamProf('PLW') # Source extractors need maps in Jy/pixel, hence they need instrument FWHM # and beam area to convert (internally) the maps. array = 'PSW' # 'PMW' 'PLW' fwhm = 17.5 # 23.9 35.2 area = 423 # 751 1587 rPeak = 22 # 30 42 # The source extractors results must be corrected for the maps pixelization, # i.e. diviving the fluxes by this factor: pixCorr = 0.941 # 0.917 0.903 # First, run sourceExtractorSussextractor on the PSW map only: this creates a source list product srcList = sourceExtractorSussextractor(image=obs.level2.refs[array].product, detThreshold=5, fwhm=fwhm, beamArea=area) # You can also drag & drop the srcList variable onto the image: it will show the sources detected # Now, given the [RA, Dec] position listed in the source list product, you can run the timeline fitter on them raCentres = srcList['sources']['ra'].data decCentres = srcList['sources']['dec'].data fluxes = srcList['sources']['flux'].data / pixCorr errors = SQRT(srcList['sources']['fluxPlusErr'].data**2 + srcList['sources']['fluxMinusErr'].data**2)/SQRT(2) chiPsw = Double1d() fluxPsw = Double1d() errPsw = Double1d() for i in range(len(raCentres)): try: param = fitter.fit(array, raCentres[i], decCentres[i], rPeak) chi = fitter.getChiSquared() except: # If the fit does not converge, a 0Jy flux will be reported in the table param = [0,0,0,0] chi = 0 # chiPsw.append(chi) fluxPsw.append(param[0] * 1.e3) errPsw.append(param[3] * 1.e3) # Finally, a create a SPIRE-only PSW catalogue spireCat = TableDataset('SPIRE 250um catalogue') spireCat.addColumn('RA', Column(raCentres)) spireCat.addColumn('Dec', Column(decCentres)) spireCat.addColumn('Sussextractor Flux [mJy]', Column(fluxes)) spireCat.addColumn('Sussextractor Error [mJy]', Column(errors)) spireCat.addColumn('Timeline Fitter Flux [mJy]', Column(fluxPsw)) spireCat.addColumn('Timeline Fitter Error [mJy]', Column(errPsw)) spireCat.addColumn('Timeline Fitter Chi2', Column(chiPsw)) # ################################################################ ############# Comparison with aperture photometry ############## # # The results must be corrected for the maps pixelization, i.e. diviving the fluxes by this factor: array = 'PSW' # 'PMW' 'PLW' pixCorr = 0.941 # 0.917 0.903 area = 423 # 751 1587 # First, convert the map to Jy/pixel dividing by the beam area mapConverted = convertImageUnit(image=obs.level2.refs[array].product, newUnit='Jy/pixel', beamArea=area) # Select/highlight mapPSWConverted in the "Variables" view and then double click on # "annularSkyAperturePhotometry" under "Applicable" in the "Tasks" view. # Click on the central source, select "Centroiding", select some values for the radii # (e.g. 30, 40, 50" for PSW) and then run the task. This is equivalent to running the following line # (centerX and centerY will change depending on where you clicked with the mouse!): result = annularSkyAperturePhotometry(image=mapConverted, centerX=91, centerY=104,\ centroid=True, fractional=1, radiusArcsec=30.0, innerArcsec=40.0, outerArcsec=50.0) # Alternatively, you can also specify a [RA, Dec] position, # e.g. the one obtained above using "Get coordinates" raCentre = str(displayWorldCoords[0]) decCentre = str(displayWorldCoords[1]) result = annularSkyAperturePhotometry(image=mapConverted, centerRA=raCentre, centerDec=decCentre,\ centroid=True, fractional=1, radiusArcsec=30.0, innerArcsec=40.0, outerArcsec=50.0) print 'Source flux (in mJy) for %s array is: %f +/- %f'%(array,\ result['Results table'][0].data[2] * 1.e3 / pixCorr,\ result['Results table'][3].data[2] * 1.e3) # ##################################################################