TWiki> Public Web>HipeUserDefs (revision 6)EditAttach

HIPE Contributed Software

This page contains contributed prototype scripts, which may one day enter in a HIPE build.

Combining PACS and SPIRE spectra

We are working on tasks to combine the spectra from PACS and SPIRE into a single spectrum. We have begun with the case of a point source, where the SPIRE and PACS data have been separately reduced through their respective pipelines and a single, point-source loss corrected spectrum has been extracted from the cubes. Combining data from these instruments makes most sense if you are dealing with wide wavelength range data, i.e. SEDs for PACS (for SPIRE you always get an SED). For PACS to cover the full SED you need two separate observations, and these would need to be reduced separately and combined.

The approach taken has been to combine all the PACS bands you have into a single Spectrum1d class product, where the bands are distinguished by segment number. The same is done for the SPIRE data. For the PACS data one can in addition "clean" up spectral regions that are affected by the light leak (see the PACS Observer's Manual to learn more about these order leaks). The combining of spectra (PACS to PACS, SPIRE to SPIRE) is done in separate scripts that one must run at the end of the respective pipelines. The two Spectrum1d products (PACS all and SPIRE all) are then combined in a third script. All of these Spectrum1d products can be opened in the Spectrum Explorer, so one can inspect the combined spectrum and the input spectra in the same environment.

At present the scripts to do this are jython functions. We do not yet mathematically "stitch" the spectra, rather they are simply pushed together. We allow for units of microns, GHz, or wavenumber on the X axis, but at present only Jy on the Y axis.

In the future these scripts will be improved, extended, and converted to HIPE tasks. We will also address the question of extended sources, which due to the different spaxel and beam sizes is a more complex question. But if you want to combine PACS and SPIRE spectra of point sources right now, you can download the scripts, read their instructions, and run them in HIPE.

This is all also explained in the PACS Pipeline Reference Manual (Chp. 3.6) and the SPIRE Data Reduction Guide (probably Chp. 6.3), which documentation you will find in the very last documentation set of the Track 7 continuous integration build of HIPE (and any build after that, of course).

The PACS part

First you must reduce your data through the pipeline up to the task "specAddNodCubes", where you end up with a single PacsRebinnedCube ListContext. From this you extract the spectrum of your point source, which most likely will be located in the central spaxel (if it is not well centred in the central spaxel, please note, your resulting spectrum will be wrong--see the PACS Data Reduction Guide Chap. 3 for more information). Then you apply the point-source loss correction to it. You do this for all the data of your source you want to combine: so, if you have two observations that cover the entire PACS range, that will give you two red and two blue spectra. You will run the following pipeline tasks (you can find these in the PACS pipeline scripts and documented in the PACS Data Reduction Guide):

slice = 0
spaxelX, spaxelY = 2,2
# get the red spectrum from the first obsid you have run the pipeline on
centralSpectrumR1 = extractSpaxelSpectrum(slicedRebinnedCubeR1, slice=slice, spaxelX=spaxelX, spaxelY=spaxelY)
centralSpectrumR1 = pointSourceLossCorrection(centralSpectrumR1, calTree=calTree)
# get the blue spectrum  from the first obsid you have run the pipeline on
centralSpectrumB1 = extractSpaxelSpectrum(slicedRebinnedCubeB1, slice=slice, spaxelX=spaxelX, spaxelY=spaxelY)
centralSpectrumB1 = pointSourceLossCorrection(centralSpectrumB1, calTree=calTree)
# get the red spectrum from the first obsid you have run the pipeline on
centralSpectrumR2 = extractSpaxelSpectrum(slicedRebinnedCubeR2, slice=slice, spaxelX=spaxelX, spaxelY=spaxelY)
centralSpectrumR2 = pointSourceLossCorrection(centralSpectrumR2, calTree=calTree)
# get the blue spectrum  from the first obsid you have run the pipeline on
centralSpectrumB2 = extractSpaxelSpectrum(slicedRebinnedCubeB2, slice=slice, spaxelX=spaxelX, spaxelY=spaxelY)
centralSpectrumB2 = pointSourceLossCorrection(centralSpectrumB2, calTree=calTree)

Next you can chose (i.e. it is not compulsory) to clean up any regions you want not to be included in the combining (e.g. the ends of the spectra where the filter bandpasses take a dive, or the ends where the light leak spectral regions are located--see the PACS Observer's Manual for more information). Finally you then combine the spectra into one Spectrum1d class product:


The scripts themselves are provided as a single attachment to this wiki page (the script: remove the .txt from the filename, which this wiki pages adds by default). You need to download the file, and then either open it in HIPE and press the double green arrow (to run the script and put the functions in to memory) or you can load them memory bu executing the file:

Two functions are in here: cleanPacsSpectrum, pacsSpectraIntoS1d.


outSpec = cleanPacsSpectrum(spectrum=inSpec, good=[], vis=False)
Where good is a list of good regions you want to be included in the eventual combination with the SPIRE spectra. An example of good is good=[60.,65.,70.,80.]. There must be an even number of good values, these being the short and long ranges of your good regions. By default the light leak regions and the very ends of the spectra are considered to be bad. Bad regions are given flags in the Spectrum1d, a flag of 1 meaning bad and a flag of 0 meaning good. This is to avoid actually cutting out data. Setting the vis parameter to True will simply produce a PlotXY of your input spectrum with the flagged-out regions indicated. Note, "outSpec" can be the same as "inSpec". "spectrum" must be of class SimpleSpectrum (as will be the case if you follow the instructions here). Various Meta data are copied from the input spectra to the output spectrum.


where "spectra" is, as the instructions above make clear, a simple list of all the SimpleSpectra that you want to combine. The output is a Spectrum1d, where the first spectrum in the list is the first segment of the Spectrum1d, the second spectrum is the second segment, and so on. It does not matter in what order your spectra are as you add them. Various Meta data are copied from the input spectra to the output spectrum. Note that it is expected that the PACS spectra are in units of micrometer vs. Jy.

The SPIRE FTS part

The script provides the interface to combine the two SPIRE FTS spectra from the central detectors SSWD4 and SLWC3. The output of this procedure is one Spectrum1d with subsequent segment numbers for SSW (194-313 μm) and SLW (303-671 μm) arrays.

Here is the full description of the script:

        combine the two FTS arrays into one Spectrum1d object, 
            each FTS array will be with a different segment number
        result = prepSpireSpec(obsContext,startSegment=0,apod="unapodized",waveUnit=None,shiftS=0.0,shiftL=0.0)
        obsContext - a SPIRE FTS observation context, sparse spatial sampling only
        startSegment - the starting index of the segment, default is 0 => SSW will be segment 0, SLW will be 1.
        apod - one of "unapodized" (default) or "apodized", which variant of the spectra to use
        waveUnit - if None then use the default unit (wavenumbers in cm-1)
                 can be "micrometer" or "GHz" for the unit of the output spectrum.
        shiftS - a flux shift to be aplied to the SSW spectrum, i.e. to match with the photometry, in the same units as the flux
        shiftL - a flux shift to be applied to the SLW spectrum, i.e. to match with the photometry
        object of class Spectrum1d with the combined SSW and SLW spectra, each with a different segment number
        Created: Apr 2011, Ivan Valtchanov, HSC, ESAC, ESA, following interactions wth Katrina Exter
        v1.0 - Apr 2011

Here is the simplest call of the script

obs = getObservation(1342189124,useHsa=True) # public NGC7027 data in the HSA
allSpireSpectra = prepSpireSpec(obsContext,waveUnit="micrometer")


The script to combine the data is run in this way:

combinedSpec = combinePacsSpireSingleSpectrum(allPacsSpectra, allSpireSpectra)
Where allPacsSpectra and allSpireSpectra here are the Spectrum1d products that are the result of the instructions given above. The script will simply push all the SPIRE and PACS spectral segments into a single Spectrum1d. Various Meta data are copied from the input spectra to the output spectrum.

You can get the script to do this from this wiki page (the script: remove the .txt from the filename, which this wiki pages adds by default). To load it into HIPE memory you either open it in HIPE and press the double green arrow (to run the script and put the function in to memory) or you can put it into memory by executing it:


The parameters of the task are:

combinedSpec = combinePacsSpireSingleSpectrum(pacs=allPacsSpectra, spire=allSpireSpectra, fluxunit="Jansky", waveunit="micrometer")
Where the "fluxunit" can take on values of "Jansky" only, and the "waveunit" can take on values of "micrometer", "GHz", and "cm-1". Note the spelling, the words must be exactly as written here.

SPIRE Photometry

Source extraction using timeline data

Tests have demonstrated that a source fitter working on the detectors' timeline works better than the map-based, such as sourceExtractorDaophot or sourceExtractorSussex. The algorithm will be included in future HIPE releases in the form of a task.

For the time being, you can use the DP/jython script written by G. Bendo: it will fit a Gaussian function to the baseline-subtracted SPIRE timelines in a SpireListContext.

Example usage

This example is based on fitting the peak of Gamma Dra in ObsID 0x50005984 in the PSW band.

fitter=bendoSourceFit(inputContext)"PSW", 269.1515617, 51.488894, 200)"PSW", 269.1515617, 51.488894, 22, 300, 350)

The first line defines an instance of the fitter object.

The second line calls a method in which the data within a 200 arcsec circle centerd on RA=269.1515617 and Dec=51.488894 is fit with a Gaussian function. The default is to fit an elliptical Gaussian function with a variable background. The first parameter will be the peak flux density.

The third line calls a methods in which a background is measured within an annulus between radii of 300 and 350 arcsec and then a Gaussian function is fit to both the central 22 arcsec and the background annulus. The default function, an elliptical Gaussian function with a variable background, is still used in this case.

See the comments at the beginning of the code to learn how to select optional functions, set parameters for the fits, or get additional data based on the resulting fits (e.g. uncertainties in the best fitting parameters).

-- KatrinaExter - 01 Apr 2011 -- IvanV - 28 Mar 2011

Topic attachments
I Attachment History Action Size Date Who Comment
Texttxt r1 manage 4.9 K 2011-04-01 - 16:56 KatrinaExter script containing function to combine PACS and SPIRE point source spectra
Texttxt r1 manage 7.6 K 2011-04-01 - 16:46 KatrinaExter script containing functions to "clean" and combine PACS spectra
Texttxt r1 manage 4.4 K 2011-04-12 - 14:22 IvanV SPIRE FTS preparation of spectra
Edit | Attach | Watch | Print version | History: r8 < r7 < r6 < r5 < r4 | Backlinks | Raw View | Raw edit | More topic actions...
Topic revision: r6 - 2011-04-14 - IvanV
This site is powered by the TWiki collaboration platform Powered by Perl