Line: 1 to 1 | ||||||||
---|---|---|---|---|---|---|---|---|
FTS Cube Analysis | ||||||||
Line: 7 to 7 | ||||||||
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 | ||||||||
Changed: | ||||||||
< < | described here, but in these examples we use the Naive | |||||||
> > | described here (Sect. 6.4.3 of the SPIRE Data Reduction Guide), but in these examples we use the Naive | |||||||
gridded data directly from the Herschel Science Archive:
obs = getObservation(1342204919, useHsa=True) | ||||||||
Line: 44 to 44 | ||||||||
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 | ||||||||
Changed: | ||||||||
< < | SpectrumFitterGUI for example (see Chapter 3 of the Herschel | |||||||
> > | 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. | ||||||||
Line: 347 to 347 | ||||||||
| ||||||||
Changed: | ||||||||
< < |
| |||||||
> > |
|
Line: 1 to 1 | ||||||||
---|---|---|---|---|---|---|---|---|
Added: | ||||||||
> > |
FTS Cube AnalysisIn 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"].productThe 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 cubeTo 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 positionMapping 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. ThespireProduct2SimpleSpectrum 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 oneIf 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")) ![]()
# 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 mapsThis 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 cubeUsing the previous example and the output from thelineFitSpec1d() 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) continueThe result of this will be a continuum removed cube. Comparison with photometer mapsTo 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
|