################################################################
#
# 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)
#
##################################################################