Public Web>SpireCalibrationWeb>FtsCubeSection9 (revision 1)~~EditAttach~~

FTS Cube Analysis

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