FTS Cube Analysis
In the following examples we will use the two SPIRE high resolution,
full spatially sampled spectral cubes from an observation of the Orion
Bar,
1342204919
. Tailoring the processing to create the cubes is
described
here, but in these examples we use the Naive
gridded data directly from the Herschel Science Archive:
obs = getObservation(1342204919, useHsa=True)
cubeS = obs.refs["level2"].product.refs["HR_SSW_unapodized_spectrum"].product
cubeL = obs.refs["level2"].product.refs["HR_SLW_unapodized_spectrum"].product
The two objects
cubeS
and
cubeL
should be of a class
SpectralSimpleCube.
Note that if you are interested in apodized spectra then you could use
"HR_SSW_apodized_spectrum"
instead. The extracted cubes can be used in
the Cube Spectral Analysis Toolbox (CSAT). The rest of this section
illustrate the use of command line and scripts to do some basic
manipulation and analysis of SPIRE spectral cubes.
Some utility functions, which are not yet available in HIPE, are needed for
the rest of this section. These are combined in a single jython file
CubeUtilitiesForFTS.py available
here. Save the file and loaded it in HIPE with
execfile()
:
execfile(path_to_file + "CubeUtilitiesForFTS.py")
or you can load it in a separate HIPE Editor window and run/execute it (i.e. press on the green double arrow button). Of course these utilities are simple jython
functions that can be modified if necessary.
The example shown here are available in a single jython file
cubeAnalysis_chapter.py, you may load it in HIPE and follow the examples.
Extracting spectra from a SPIRE spectral cube
To extract a spectrum from a cube pixel (also known as a
spaxel) you can use
the following method:
spec = spireProduct2SimpleSpectrum(cubeL, spaxelX=4, spaxelY=4)
The output spectrum will be of a
SimpleSpectrum class and can be used in most
of the spectral analysis tools, like the
SpectrumExplorer
and the
SpectrumFitterGUI
for example (see
Chapter 3 of the Herschel
Scripting and Data Mining Guide). The extracted spectrum will be with the
extended source calibration and the flux unit will be W m
-2
Hz
-1 sr
-1.
Alternatively, the spectrum can be extracted as a
Spectrum1d:
spec = cube.getSpectrum1d(4,4)
Extracting cube spectra from a sky position
Mapping observations produce two cubes: one for the short wavelength array
(SSW) and one for the long wavelength array (SLW). As discussed earlier, the
two cubes have different WCS and the extraction of spectra from the same sky
position is not easy to do by specifying spaxel numbers. The
spireProduct2SimpleSpectrum
method can also be used with sky positions
(the position must be within the cube):
# RA in degrees
ra = 15.0*(5.0 + 35.0/60.0 + 25.2/3600.0)
# Dec in degrees
dec = -5.0 - 24.0/60.0 - 29.0/3600.0
#
ssw = spireProduct2SimpleSpectrum(cubeS, ra=ra, dec=dec)
slw = spireProduct2SimpleSpectrum(cubeL, ra=ra, dec=dec)
The two calls will extract two spectra from the nearest spaxel from the
corresponding cube. The two spectra can be combined in a single product
following the instructions in Section 6.5.2,
"Joining SSW and SLW data". In
CubeUtilitiesForFTS.py we provide an alternative approach that combines
the SSW and SLW spectra in a single product, assigning a different segment
number to SLW and SSW, thus no assumptions of the overlap region are needed.
Here is an example on how to combine them:
ftsList = {'SSW': ssw, 'SLW': slw}
combined = combineFts(ftsList)
Double-click on
combined
in the
Variables view
will show both SSW and SLW spectra in the Spectrum Explorer, each
spectrum will have a different segment number and a different
colour. Combining spectra from SSW and SLW must be considered with
caution, because each spaxel is an average of spectra coming from
different detectors and the default SLW and SSW pixel sizes can be
quite different, i.e. for high resolution full coverage spectral
maps the default SSW pixel size is 9.5 arcsec, while the SLW pixel
size is 17.5 arcsec.
Converting an extracted spectrum to a point source calibrated one
If the spectral map contains a point source (or semi-extended
source) then it makes sense that the extracted spectrum from the
cube is converted from the default extended source flux calibration,
that provides spectra in units of W m
-2 Hz
-1 sr
-1 , to a
point source flux calibration, i.e. spectrum in units of Jy. To do
this we need the conversion factors that are provided in a
calibration file togehter with the Spectrometer beam parameters. The
way to use it is illustrated here:
# extract a spectrum from spaxel (4,4)
spec1d = spireProduct2SimpleSpectrum(cubeL,spaxelX=4,spaxelY=4)
#
cal = spireCal(pool="spire_cal_9_0")
# cal = spireCal() # to take the calibration context form the HSA
beamPar = cal.spec.beamParamList.getProduct(spec1d)
# get the point source conversion for the central detector of SLW (our
# cube is SLW).
pointConv = beamPar.getDefault()["SLWC3"]["pointConv"].data
#
# now copy and populate the output point source calibrated spectrum
pointSpec = spec1d.copy()
pointSpec.setFlux(spec1d.getFlux()*pointConv)
pointSpec.setFluxUnit(herschel.share.unit.Unit.parse("Jy"))

There are some important considerations to be noted here:
- we use the point source conversion for the central bolometer of the SLW array, namely SLWC3. The spectrum in (4,4) is an average of different detectors that may not contain SLWC3. The differences of the conversion factors are small but there are differences. The better approach would be to convert the spectra of all detectors at level-1 before the gridding part of the reprocessing. With this approach, the cube after the gridding task will have units of Jy and can be used directly to measure lines and continuum assuming a point source calibration.
- we assume that the total point source flux is contained in a spaxel (4,4). Alternatively, one can add spaxels and then perform the conversion, like this example where we add spectra in a circle with radius of 2 pixels (i.e. 35 arcsec for a default SLW cube grid), centred on (4,4):
# some necessary input classes
from herschel.ia.toolbox.cube.ExtractRegionSpectrumTask import Region, Arithmetics
# add up spectra from a circular aperture
region = extractRegionSpectrum(cube=cubeL, regionType=Region.ELLIPSE, \
arithmetics=Arithmetics.SUM, \
centerRow=4.0, centerCol=4.0, radiusRow=2.0, radiusCol=2.0)
Simultaneous fit of lines and continuum: line intensity and line velocity maps
This short script loops over all cube pixels, simultaneously
fits a list of lines and a suitable degree polynomial for the
continuum and saves an image of a selected lines. The script uses
one of the utility function in
CubeUtilitiesForFTS.py,
lineFitSpec1d()
that is a wrapper of the line fitting
useful script. This function fits simultaneously a continuum and a
number of input lines to a spectrum and outputs a number of useful
results in a
CompositeDataset
: a table with the fit
results for each line, with line centroid (useful for velocity maps),
line flux and error, and also the residual (i.e. model-data and the
fitted continuum. Within the
lineFitSpec1d()
function
the input spectrum will be corrected for the radial velocity of the
Herschel , while the correction for redshift or
source velocity are optional input parameters. In the following, we
will use the SLW cube and the example is just an illustration of the
script approach for cube analysis. And we need to prepare a line
list, with lines that we are interested in:
# Define the rest frame line frequencies (GHz) of the 12CO ladder
lineFreqs = Double1d([461.0407682,576.2679305,691.4730763,806.6518060,921.7997000,\
1036.9123930,1151.9854520,1267.0144860,1381.9951050,1496.9229090])
lineNames = String1d(["12CO(4-3)","12CO(5-4)","12CO(6-5)","12CO(7-6)","12CO(8-7)",\
"12CO(9-8)","12CO(10-9)","12CO(11-10)","12CO(12-11)","12CO(13-12)"])
Note that you can add or remove lines from this array. And it
is always good to check the line and continuum fit with a selected
spaxel, so that no strange things happen with the fit or the initial
parameters for the continuum or the lines produce acceptable
results.
slw = spireProduct2SimpleSpectrum(cubeL, spaxelX=4, spaxelY=4)
fitOutput = lineFitSpec1d(spec1d,lineFreqs,lineNames,polyOrder=4,plotIt=True)
This will perform the line and continuum fit (continuum as a
4rd order polynomial) and plot the result: the total model, the
residual and the input spectrum. If the results are acceptable then
you can continue with making the same fit on all cube
pixels.
#
# let's get the cube dimensions
#
size = cubeL.getDimensions()
nx = size[2] # X,Y convention in HIPE
ny = size[1]
# and then create the empty output arrays
lineImage = Double2d(nx,ny) # this will be the line intensity array
velImage = Double2d(nx,ny) # and this will be the line velocity map
#
# now the loop over all cube pixels
#
iLine = 2 # the index of the line we are interested in, 12CO(6-5) in this case
#
for i in range(nx):
for j in range(ny):
# extract the 1-dimensional spectrum from a spaxel
spec1d = spireProduct2SimpleSpectrum(cubeL, spaxelX=i,spaxelY=j)
# we have to take care of those pixels that are NaN and the fit will fail.
# That's why we use try: except:
try:
fitOutput = lineFitSpec1d(spec1d,lineFreqs,lineNames,polyOrder=4,\
plotIt=False) # not we don't want to plot here
lineImage[i,j] = fitOutput["FitTable"]["Flux"].data[iLine]
velImage[i,j] = fitOutput["FitTable"]["Velocity"].data[iLine]
except:
lineImage[i,j] = Double.NaN velImage[i,j] = Double.NaN
print "Fit failed for pixel (%i,%i)"%(j,i) continue
# and now, fill in the line intensity map and copy the cube WCS
co65_image = SimpleImage(description="Line intensity map for 12CO(6-5)")
co65_image.setImage(lineImage)
co65_image.setWcs(cubeL.getWcs())
disp = Display(co65_image,True,title="CO(6-5) line intensity map")
disp.zoomFit()
# and the velocity map
co65_velMap = SimpleImage("Line velocity map for 12CO(6-5)")
co65_velMap.setImage(velImage)
co65_velMap.setWcs(cubeL.getWcs())
disp2= Display(co65_velMap,True,title="CO(6-5) velocity map")
disp2.zoomFit()
Note that
fitOutput["FitTable"]
contains the
information for all fitted lines and it is not necessary to loop again
for another line, but just add a new line in the
try:
block with an additional
Double2d()
array.
Another important note is that the output line intensity image
will be in units of W m
-2 sr
-1, because the input cubes are in extended
source calibration, which means units of W m
-2 Hz
-1 sr
-1. To convert this to W
m
-2 you will need to know the beam size at the line frequency.
Continuum subtracted cube
Using the previous example and the output from the
lineFitSpec1d()
function, we can easily populate a new cube
with the same dimensions as the original one, but replacing the
spaxels with the continuum removed ones. Here is the fill
example:
#
# let's get the cube dimensions
#
size = cubeL.getDimensions()
nx = size[2] # X,Y convention in HIPE
ny = size[1]
nwave = size[0] # frequency dimension
#
# copy the input cube where we will store the continuum removed spectra
#
newCube = cubeL.copy()
#
# now the loop over all cube pixels
#
for i in range(nx):
for j in range(ny):
# extract the 1-dimensional spectrum from a spaxel
spec1d = spireProduct2SimpleSpectrum(cubeL, spaxelX=i,spaxelY=j)
# we have to take care of those pixels that are NaN and the fit will fail.
# That's why we use try: except:
try:
fitOutput = lineFitSpec1d(spec1d,lineFreqs,lineNames,polyOrder=4,\
plotIt=False) # not we don't want to plot here
newCube.setFlux(j,i,fitOutput["contRemoved"].getFlux())
except:
newCube.setFlux(j,i,Double1d(nwave,Double.NaN))
print "Fit failed for pixel (%i,%i)"%(j,i)
continue
The result of this will be a continuum removed cube.
Comparison with photometer maps
To compare a spectral map with a photometer map we have to
multiply each spectrum with the SPIRE Photometer filter
transmission. Note that this comparison is only qualitative because
there are a number of complications that have to be taken into
account: the spectrometer beam may vary significantly within the
Photometer transmission curve, the spectrometer map is using an
extended source calibration etc. The short script uses another
utility function from
CubeUtilitiesForFTS.py collection of methods,
getFtsPhot()
. The input to this function is a 1-d
spectrum and a calibration context, containing the Photometer filter
transmissions and the output is a python dictionary with three keys:
"PSW", "PMW" and "PLW" with the result. Note that running this
method on SLW spectrum will produce two numbers (for "PMW" and "PLW")
and one
NaN for "PSW" as the Photometer filter transmission coverage
is outside the spectral limits. The same way, the method will return
only one number for a spectrum from SSW.
#
# let's get the cube dimensions
#
size = cubeL.getDimensions()
nx = size[2] # X,Y convention in HIPE
ny = size[1]
#
# create output maps, only PMW and PLW as we work with SLW
s350image = Double2d(nx,ny)
s500image = Double2d(nx,ny)
#
# now the loop over all cube pixels
#
# load the calibration context from a local pool
cal = spireCal(pool="spire_cal_9_0")
# or take the cal context form the HSA
#cal = spireCal()
#
for i in range(nx):
for j in range(ny):
# extract the 1-dimensional spectrum from a spaxel
spec1d = spireProduct2SimpleSpectrum(cubeL, spaxelX=i,spaxelY=j)
# use try: except: for pixels that are NaNs
try:
photo = getFtsPhot(spec1d,cal=cal,extended=True)
# we use the extended lambda^2 weighting of the Photometer filters
s350image[j,i] = photo['PMW']
s500image[j,i] = photo['PLW']
except:
print "Skip pixel (%i,%i)"%(j,i)
# now create the actual HIPE images
map350 = SimpleImage(description="Photmetry map from FTS")
map350.setImage(s350image)
map350.setWcs(cubeL.getWcs())
disp3 = Display(map350,True,title="Photometry map at 350$\mu$m")
disp3.zoomFit()
map500 = SimpleImage(description="Photmetry map from FTS")
map500.setImage(s500image)
map500.setWcs(cubeL.getWcs())
disp4 = Display(map500,True,title="Photometry map at 500 $\mu$m")
disp4.zoomFit()
--
IvanV - 18 May 2012