FMT_Functions.py000644 000765 000024 00000010216 13276075100 014703 0ustar00katrinastaff000000 000000 ## For cubes from Level 2, the PACS absolute flux calibration is established
## for a perfectly extended point source (flat surface brightness). This calibration
## contains a specific "Extended Source Correction" step to account for the spatial beam
## inefficiency , i.e. to account for the fact that the spaxels don't have a perfect
## and flat response over 9.4"x9.4", but that some flux is 'lost between the spaxels'
## If you are computing your own "Extended Source Correction" in the Forward Modelling Scripts
## then this ESC should be removed before you can apply your unique correction.
def convCoord2Decimal(ra,dec):
rad = 360.*(ra[0]+(ra[1]+(ra[2]/60.))/60.)/24.
decd=ABS(dec[0]) +( (dec[1]+(dec[2]/60.))/60.)
if (dec[0] < 0): decd=(-1)*decd
return(rad,decd)
def convCoordFromDecimal(ra,dec):
ra1 = int(24.*(ra/360.))
ra2 = int(60.*(24.*(ra/360.) - ra1))
ra3 = 60*(60.*(24.*(ra/360.) - ra1) - ra2)
decc=dec
if (dec < 0.0): decc = dec * (-1)
dec1 = int(decc)
dec2 = int(60*(decc - dec1))
dec3 = 60*((60*(decc - dec1)) - dec2)
if (dec < 0.0): dec1 = (-1)*dec1
return(str(ra1)+":"+str(ra2)+":"+str(ra3),str(dec1)+":"+str(dec2)+":"+str(dec3))
from herschel.pacs.spg.spec.ExtractCentralSpectrumTask import *
calTree=getCalTree()
## For slicedRebinnedCubes
def undoExtendedSourceCorrectionSliced(slicedRebinnedCubes,calTree=calTree,verbose=1):
go = slicedRebinnedCubes.meta["extendedSourceCorrected"].value
if go:
slicedRebinnedCubesOut = slicedRebinnedCubes.copy()
for slice in range(len(slicedRebinnedCubesOut.refs)):
if verbose: print 'slice',slice
cubout = undoExtendedSourceCorrection(slicedRebinnedCubesOut.get(slice), calTree=calTree, verbose=verbose)
slicedRebinnedCubesOut.replace(slice,cubout)
slicedRebinnedCubesOut.meta["extendedSourceCorrected"].value = False
return slicedRebinnedCubesOut
else:
print "The input cubes are not extended source corrected: nothing will be done"
## For individual cubes
def undoExtendedSourceCorrectionCube(cube, calTree=None, verbose=0):
if (("extendedSourceCorrected" in cube.meta.keySet()) and (cube.meta["extendedSourceCorrected"].value == 1)):
cube=cube.copy()
if not calTree: calTree= getCalTree()
if calTree.version < 70: print "Please update your calibration tree version "+str(calTree.version) +" to version >= 70"
extLoss = calTree.spectrometer.extendedSourceLoss
extipol = LinearInterpolator(extLoss.getWavelength(), extLoss.getSignalFraction(), True)
wave, flux = cube.getWave(), cube.getFlux()
#
for x in range(5):
for y in range(5):
interpolatedFractions = extipol(wave[:,x,y])
flux[:,x,y] *= interpolatedFractions
cube.setFlux(flux)
cube.setFluxDescription("Flux")
cube.meta["extendedSourceCorrected"].value = False
else:
print "The extended source correction was not applied to these data,"
print "It will hence not be undone here"
return cube
## For mosaic cubes
def undoExtendedSourceCorrectionMosaicCube(cube, calTree=None, verbose=0):
if (("extendedSourceCorrected" in cube.meta.keySet()) and (cube.meta["extendedSourceCorrected"].value == 1)):
cube=cube.copy()
if not calTree: calTree= getCalTree()
if calTree.version < 70: print "Please update your calibration tree version "+str(calTree.version) +" to version >= 70"
extLoss = calTree.spectrometer.extendedSourceLoss
extipol = LinearInterpolator(extLoss.getWavelength(), extLoss.getSignalFraction(), True)
wave, flux = cube.getWave(), cube.getFlux()
rows = cube.dimensions[1]
cols = cube.dimensions[2]
#
for r in range(rows):
for c in range(cols):
interpolatedFractions = extipol(wave[:,r,c])
flux[:,r,c] *= interpolatedFractions
cube.setFlux(flux)
cube.setFluxDescription("Flux")
cube.meta["extendedSourceCorrected"].value = False
else:
print "The extended source correction was not applied to these data. It will hence not be undone here"
return cube
FMT_SourceModels.py000644 000765 000024 00000010066 13275053155 015347 0ustar00katrinastaff000000 000000 from herschel.pacs.spg.spec.beams import ExtendedSpectralSourceAdapter
from math import *
## A Gaussian source with a Gaussian emission line + polynomial continuum
## User-defined parameters are the FWHM of the Gaussian, its offset from the
## centre of the field (the centre of the field will be set when the model
## is used in the task pacsSpecFromModel, and is the centre of the field
## of the cube you are modelling), the flux values of the line and continuum
##
class MyGaussSpectral(ExtendedSpectralSourceAdapter):
fwhm = 0.0
offRa = 0.0
offDec = 0.0
contFlux = 0.0
contSlope = 0.0
lineAmplitude = 0.0
lineCenter = 0.0
lineWidth = 0.0
def __init__(self,fwhm,offRa,offDec,contFlux,contSlope,lineAmplitude,lineCenter,lineWidth):
self.sourceSize = fwhm/2.3548
self.offRa = offRa
self.offDec = offDec
self.contFlux = contFlux
self.contSlope = contSlope
self.lineAmplitude = lineAmplitude
self.lineCenter = lineCenter
self.lineWidth = lineWidth
# Required
def getFlux(self, relRa, relDec, wavelength):
r = sqrt((relRa-self.offRa)**2 + (relDec-self.offDec)**2)
spectrum = self.contFlux + self.contSlope * (wavelength-self.lineCenter)
spectrum += self.lineAmplitude * exp(-0.5 * ((wavelength-self.lineCenter) / self.lineWidth)**2)
flux = spectrum * exp(-0.5 * (r / self.sourceSize)**2)
return flux
def getFluxImage(self,relRa,relDec, wavelength):
r = sqrt((relRa-self.offRa)**2 + (relDec-self.offDec)**2)
spectrum = self.contFlux + self.contSlope * (wavelength-self.lineCenter)
spectrum += self.lineAmplitude * exp(-0.5 * ((wavelength-self.lineCenter) / self.lineWidth)**2)
self.fluxImage = spectrum * exp(-0.5 * (r / self.sourceSize)**2)
return self.fluxImage
# Rendering of layer on the beam grid - optional
def renderGaussLayer(self, size=109, dpix=0.5, wavelength=66.38):
self.image=Double2d(size, size)
for y in range( size ):
for x in range( size ):
xim = (-dpix)*(x-(size/2))
yim = (dpix)*(y-(size/2))
self.image.set(y,x, self.getFluxImage(xim, yim, wavelength))
return self.image
# Rendering of layer image on the beam grid - optional
def getSourceModelImage(self, size=109, dpix=0.5, wavelength=66.38):
return self.renderGaussLayer(size, dpix, wavelength)
## Model visualisation
def imageModelInPacsSpecFromModel(sourceModel,naxis, pixAngSize, wave):
sourceModelImage_01 = SimpleImage(image=sourceModel.getSourceModelImage(naxis*5, pixAngSize/5, wave))
res = Double2d(naxis,naxis)
for i in range(naxis):
for j in range(naxis):
flux = 0.0
for p in range(5):
for q in range(5):
i2=i*5+p
j2=j*5+q
flux += sourceModelImage_01.image[j2,i2]
res[j,i] = flux/25.
sourceModelImage_05 = SimpleImage(image=res)
return sourceModelImage_05
## A completely flat source
class MyFlatSourceSpectral(ExtendedSpectralSourceAdapter):
def __init__(self,contFlux,contSlope,lineCenter):
self.contFlux = contFlux
self.contSlope = contSlope
self.lineCenter = lineCenter
def getFlux(self, relRa, relDec, wavelength):
spectrum = self.contFlux + self.contSlope * (wavelength-self.lineCenter)
return spectrum
def getFluxImage(self,relRa,relDec, wavelength):
spectrum = self.contFlux + self.contSlope * (wavelength-self.lineCenter)
self.fluxImage = spectrum
return self.fluxImage
def renderGaussLayer(self, size=109, dpix=0.5, wavelength=79.0):
self.image=Double2d(size, size)
for y in range( size ):
for x in range( size ):
xim = (-dpix)*(x-(size/2))
yim = (dpix)*(y-(size/2))
self.image.set(y,x, self.getFluxImage(xim, yim, wavelength))
return self.image
def getSourceModelImage(self, size=109, dpix=0.5, wavelength=79.0):
return self.renderGaussLayer(size, dpix, wavelength)
FMT_userScript.py000644 000765 000024 00000102667 13276075061 015120 0ustar00katrinastaff000000 000000 #
# This file is part of Herschel Common Science System (HCSS).
# Copyright 2001-2011 Herschel Science Ground Segment Consortium
#
# HCSS is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as
# published by the Free Software Foundation, either version 3 of
# the License, or (at your option) any later version.
#
# HCSS is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General
# Public License along with HCSS.
# If not, see .
#
"""
This script is to explain and demonstrate the Forward Modelling that users of PACS spectroscopy
of extended sources can perform to improve the flux calibration of their data. This script and why it
is necessary are explained in the User Note "Dealing with Extended Sources", which can be found on the
Level 2 PACS section of the Herschel Explanatory Legacy Library. The plots and images created here are
also included in that user note.
https://www.cosmos.esa.int/web/herschel/legacy-documentation-pacs-level-2
https://www.cosmos.esa.int/web/herschel/pacs-overview
In this "Forward modelling script" you will estimate a flux calibration efficiency for your observation,
based on a model of the source flux distribution and its coupling to the beam efficiency for every spaxel
and in every pointing in the observation. The first is an input from the user, the second is taken from
the calibration tree. This flux calibration efficiency is included in the standard calibration applied to
all PACS spectroscopy products, but is only fully valid for flat sources. For any other source shape, the
efficiency depends on where the source lies on the detector, and hence on its shape and location in the
pointing pattern of the observation. This calibration is called the "Extended Source Calibration" (ESC).
The forward modelling requires that the user knows the shape of the source, and ideally can provide a "true"
spectrum for the source. The main correction that is computed here is to the area-integrated
flux of the source, not its full 2d distribution, but the relative 2d losses can also be estimated
following an example given at the end. If the user does not know the true spectrum of the source, then
a standard spectrum can be used and a relative correction can be computed.
In this script we will demonstrate that the process works using a simulated pointed observation of R Doradus,
which is a point source. The reason for using a point source as a demonstration is that we know what
the true spectrum of R Dor is -- via the point-source corrections that are provided for all PACS spectroscopy
observations of point sources. We can then compare the result to this known absolute.
Four scenarios are tested and demonstrated here. The first two demonstrate that the technique works and introduces
the process. The second two demonstrate the most likely ways it will be used.
- proof of concept: using a known input of a centred point source, we will show how to do the forward
modelling and that it works
- testing with an off-centred source: we show that the correction also works for off-centred sources,
that is, it can correct for the illumination pattern at different locations in the focal plane, not just
the centre
- using a standard spectrum: to demonstrate what to do if you do not know, a-priori, what the spectrum
of the source truely is (which the previous 2 cases assumed)
- showing how to run this script on a mosaic cube as if from a mapping observation (the previous 3 cases were
all for pointed observations and on a rebinned cube)
Finally, in a small bit at the end we show how to estimate the illumination pattern for your obseration to see what
the relative 2D distribution of the correction is.
(Note that despite this script using a point source as a test case, we do not recommend using this script to compute
flux corrections for off-centred point sources. The other methods recommended in the User Note on Point Sources should
be followed)
This script is based on the extensive testing work of E. Puga
K Exter. HSC. May 2018
"""
#######################
# Imports and set ups #
#######################
"""
Importing some functions written into two scripts
The first contains functions to remove the ESC from various types of cubes
The second sets up some models, to which you can add your own if you like
"""
execfile("/Users/katrina/ownCloud/HSC_work/FMT/4Users/FMT_Functions.py")
execfile("/Users/katrina/ownCloud/HSC_work/FMT/4Users/FMT_SourceModels.py")
calTree=getCalTree(version=77)
##########################################################################################################################
# Proof of concept
##########################################################################################################################
"""
We first show how the process works and also that it produces the correct result when all the inputs are also correct
"""
################
# Get the data #
################
obsid = 1342246386 # R Doradus, line scan, mapping observation (5x5 of 2.5"x2.5")
obs = getObservation(obsid, useHsa=1)
camera = 'blue'
calTree=getCalTree()
level2 = obs.level2
slicedRebinnedCubes_init = level2.getCamera(camera).rcube.product.copy()
# Select the central raster position to emulate a pointed observation. In this slice, R Dor is centred in the FoV
slicedRebinnedCubes = selectSlices(slicedRebinnedCubes_init, rasterLine=[3], rasterCol=[3])
# If you want to save for later use -
#saveSlicedCopy(slicedRebinnedCubes,"slicedRebinnedCubes")
# Read back in with -
#slicedRebinnedCubes = readSliced("slicedRebinnedCubes")
# cleaning up
del(obs,level2)
#######################
# First manipulations #
#######################
"""
Undo the ESC - since we want to compute a uniqe ESC for our source, the first step
is to remove the pipeline-standard ESC.
"""
slicedRebinnedCubes_NoESC = undoExtendedSourceCorrectionSliced(slicedRebinnedCubes, calTree=calTree)
#########################
# Create input spectrum #
#########################
"""
The input spectrum for the modelling will be the true spectrum of R Dor, which we calculate using the point-source
correction task.
"""
c1,c9,c129 = extractCentralSpectrum(slicedRebinnedCubes_NoESC.get(0), isLineScan=-1, calTree=calTree, verbose=0)
specPSC = c9
# (To know more about this task, see the User Note on Working with Point Sources)
"""
Summing up the FoV
This is what the correction computed with the forward modelling will be applied to. Extract a total spectrum of an
area of the cube where your source lies. In this case, since the point source is bright and there is nothing else in the
FoV, we will simply sum up the entire cube. If you need to extract a particular area of the cube, read the documentation
on the Cube Toolbox in the Data Analysis Guide of the help documentation.
You will only need the specSUM_noesc spectrum, we create the other here simply to show you the difference between the two
"""
specSUM_noesc = extractRegionSpectrum(cube=slicedRebinnedCubes_NoESC.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
specSUM_spg = extractRegionSpectrum(cube=slicedRebinnedCubes.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
# Plot
p=PlotXY(titleText="Comparing observed spectra")
p.addLayer(LayerXY(specPSC.getWave(),specPSC.getFlux(),name="point source corrected"))
p.addLayer(LayerXY(specSUM_noesc.getWave(),specSUM_noesc.getFlux(),name="sum of rebinned cube (noESC)"))
p.addLayer(LayerXY(specSUM_spg.getWave(),specSUM_spg.getFlux(),name="sum of rebinned cube (SPG)"))
p.xaxis.titleText="Wavelength ($\micro$m)"
p.yaxis.titleText="Flux density (Jy)"
p.legend.visible = 1
# We need to define the wavelength, peak flux, and sigma width of the line, and the continuum flux level for the model
# Fitting specPSC with the SpectrumFitterGUI will do this
# -> 66.45mum peak 442.0Jy sigma 0.0078mum, continuum 160Jy
###################################
# Coordinates of the model/source #
###################################
"""
An input to the forward modelling is the offset of the centre of the source from a reference position. This reference
position by default is the value of raNominal and decNomial in the meta data of the cube. It can also be a input directly
into the task, although in our experience this does not seem to always work.
For R Dor, what are the coordinates of the star compared to the meta data values?
For our rebinned and projected cubes we have:
>ra = slicedRebinnedCubes.get(0).getMeta()["raNominal"].value
>dec = slicedRebinnedCubes.get(0).getMeta()["decNominal"].value
>print convCoordFromDecimal(ra,dec) # ('4:36:45.59', '-62:04:37.8')
The position for R Dor from SIMBAD is 04 36 45.59127 -62 04 37.7974, i.e. the same as the meta data values.
Now you need to determine the star's actual position in the FoV. Is it actually positioned at the value
raNominal and decNominal? This is extremely difficult to determine for a pointed observation, since
the spatial sampling is so low: the PACS blue beam FWHM is about 9" and the spaxels are 9.4" in size, and this
sets a strong limit to the accuracy with which one can determine a position of a source. In this script we cheated:
the observation we took the data from is a mapping observation, so we looked at the projected cube in the Level 2
(this cube has 3" spaxels and was created as the mosaic of many raster position, and so has a good spatial sampling)
and from that could see where the star is located. Doing this, we find the following offsets:
about 0" in RA and -1.4" in Dec.
"""
######################
# Create input model #
######################
"""
To model a point source, use the Gaussian model (MyGaussSpectral) with a FWHM (sourceSize) of 0.2", as you get flux
conservation problems for smaller values.
For the flux values, it is necessary to normalise the values measured from the spectra. The fluxes we
measured are for the entire source, so in units of Jy, but the model requires Jy/square_arcsec, so you have
to normalise by the area of the source on the sky.
"""
## For the simulated pointed observation with the PSC spectrum values as input
sourceSize = 0.2
offRa = 0.0
offDec = -1.4
contFlux = 160.0/(2*Math.PI*(sourceSize/(2*SQRT(2*LOG(2))))**2)
contSlope = 0.0
lineAmplitude = 442.0/(2*Math.PI*(sourceSize/(2*SQRT(2*LOG(2))))**2)
lineCentre = 66.45
lineWidth = 0.0078
sourceModel_pscSingle = MyGaussSpectral(sourceSize,offRa,offDec,contFlux,contSlope,lineAmplitude,lineCentre,lineWidth)
#
d=Display(imageModelInPacsSpecFromModel(sourceModel_pscSingle, 109, 0.5, lineCentre), title='continuum layer Gaussian model')
######################
# Create model cubes #
######################
"""
Forward model using the source model. This takes the input model and simulates what a real observation would see,
by "folding in the beams".
Note: you can also input the centre of the cube for the model using the
parameters raSource and decSource (decimal values are required), rather than use the raNominal and
decNominal in the cube, but we found that using the default raNom and decNom is more reliable
"""
simRebinnedCubes = pacsSpecFromModel(slicedRebinnedCubes_NoESC,sourceModel_pscSingle,calTree=calTree)
#################
# Check results #
#################
"""
Before creating the correction from this result, let us check to see if the simulated cube looks like the observed cube.
Remember: we have input the true spectrum of R Dor, and the task pacsSpecFromModel simulates what PACS would see
for that source, so the simulated and observed cubes (no ESC) should look the same. We will do the comparison along the
spectral and spatial axes
"""
# Spectral - the central spaxel spectrum, the central 3x3, and the total
p1 = PlotXY(titleText='Simulated vs Observed')
spec1 = extractRegionSpectrum(cube=simRebinnedCubes.get(0), regionType=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Region.SINGLE_PIXEL, centerRow=2.0, centerCol=2.0)
spec2 = extractRegionSpectrum(cube=slicedRebinnedCubes_NoESC.get(0), regionType=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Region.SINGLE_PIXEL, centerRow=2.0, centerCol=2.0)
p1.addLayer(LayerXY(spec1.wave,spec1.flux,name="simulated cube"),0,0)
p1.addLayer(LayerXY(spec2.wave,spec2.flux,name="observed (no esc)"),0,0)
spec1 = extractRegionSpectrum(cube=simRebinnedCubes.get(0), regionType=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Region.RECTANGLE, arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM, row1=0.5, col1=0.5, row2=3.5, col2=3.5)
spec2 = extractRegionSpectrum(cube=slicedRebinnedCubes_NoESC.get(0), regionType=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Region.RECTANGLE, arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM, row1=0.5, col1=0.5, row2=3.5, col2=3.5)
p1.addLayer(LayerXY(spec1.wave,spec1.flux,name="simulated cube"),1,0)
p1.addLayer(LayerXY(spec2.wave,spec2.flux,name="observed (no esc)"),1,0)
spec1 = extractRegionSpectrum(cube=simRebinnedCubes.get(0),arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
spec2 = extractRegionSpectrum(cube=slicedRebinnedCubes_NoESC.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
p1.addLayer(LayerXY(spec1.wave,spec1.flux,name="simulated cube"),2,0)
p1.addLayer(LayerXY(spec2.wave,spec2.flux,name="observed (no esc)"),2,0)
p1.legend.visible=1
p1[0].xaxis.titleText="Wavelength ($\micro$m)"
p1[0].yaxis.titleText="Jy (central spaxel)"
p1[2].xaxis.titleText="Wavelength ($\micro$m)"
p1[2].yaxis.titleText="Jy (cental 3x3 spaxels)"
p1[4].xaxis.titleText="Wavelength ($\micro$m)"
p1[4].yaxis.titleText="Jy (entire cube)"
"""
The ratio of integrated fluxes of observed to simulated (measured using the Spectrum Fitter GUI which returns an
integrated value in units of Jy*micron) of sim/obs is
central: 0.95
3x3 : 0.99
total : 0.99
The worst match is for the central spaxel, which is almost certainly because we have not been able to
determine exactly what the offset between the star and raNominal and decNomial are. This highlights the fact that
small-area comparisons reqiure a much higher degree of accuracy in the spatial inputs to the model.
"""
# Spatial - create a 2D map of the source by collapsing along the spectral dimension
simIntegratedMap = integrateSpectralMap(cube=simRebinnedCubes.refs[0].product, startArray=Double1d([66.418]), endArray=Double1d([66.483]))
simMap = simIntegratedMap.getSimpleImage(0)
obsIntegratedMap = integrateSpectralMap(cube=slicedRebinnedCubes_NoESC.refs[0].product, startArray=Double1d([66.418]), endArray=Double1d([66.483]))
obsMap = obsIntegratedMap.getSimpleImage(0)
d1=Display(simMap)
d1.title = "simulated"
d2=Display(obsMap)
d2.title = "observed"
"""
OK, so now to the correction.
What are the various integrated fluxes?
simulated noESC SPG(with ESC) correct (PSC)
central spaxel 5.58 5.89 - -
central 3x3 7.76 7.71 - -
total field 8.03 8.13 10.58 8.64
The difference in total flux between
the noESC and the simulated is : 1.01 GOOD
the SPG and the simulated is : 1.32 This is the difference that the illumination pattern makes to an observation
the simulated and the correct is : 0.93
The unique ESC for our source is the ratio of the simulated cube's line flux to the input spectrum's line flux.
So the correction is
Observed flux: 8.13
Correction: 0.926
Corrected flux : 8.13/0.926 = 8.78 -> which is close to 8.64 (which is the correct answer)
"""
#################################################################################################################################
# Testing on an offset source
#################################################################################################################################
"""
The centre of the FoV is the best location in the illumination pattern (see the figures in the User Note). A source located here will
have the smallest correction. If the formward modelling works, then the corrected spectrum from R Dor in a cube where it is
located elsewhere in the FoV -- somewhere where the illumination is worse -- should produce the same result as the corrected
spectrum produced for a centred source.
Let us therefore test this. We will select a cube from the slicedRebinnedCube_init of this mapping observation where R Dor is
off-centred. We will then determine the correction using the specPSC as input as we did before -- because this is still the
correct spectrum of R Dor -- and then compare the result here to the result of before.
We have selected the cube slice number 2, as the source is offset by what appears to be about 1 spaxel (=9.4"). However, this offset
will have to be determined more accurately for the modelling to work.
(Unfortunately, point sources that are offset from the centre of the FoV tend to develop a skew in their spectral lines and this
adds to the difficulty of locating the source in the FoV, as you will see later. This can also affect small sources)
"""
##############
# Input cube #
##############
# using slicedRebinnedCubes_init from the beginning of this script
slicedRebinnedCubes_off = selectSlices(slicedRebinnedCubes_init, sliceNumber=[2])
slicedRebinnedCubes_off_NoESC = undoExtendedSourceCorrectionSliced(slicedRebinnedCubes_off, calTree=calTree)
specSUM_off_noesc = extractRegionSpectrum(cube=slicedRebinnedCubes_off_NoESC.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
"""
Comparing specSum_off_noesc to specSum_noesc created from the centred cube of the previous test, it is clear that
the flux from this off-centred one is lower. This is as expected.
"""
p=PlotXY(titleText="Comparing observed spectra")
p.addLayer(LayerXY(specSUM_noesc.getWave(),specSUM_noesc.getFlux(),name="centred"))
p.addLayer(LayerXY(specSUM_off_noesc.getWave(),specSUM_off_noesc.getFlux()+10.0,name="off-centred(+10)"))
p.xaxis.titleText="Wavelength ($\micro$m)"
p.yaxis.titleText="Flux density (Jy)"
p.legend.visible = 1
"""
FYI: The flux in the observed spectral line is 7.09 (as measured by the Spectrum Fitter GUI). The true flux is 8.64
(specPSC) and the line flux in the spectrum from the centred observation (the test at the top of this script) was 8.13.
"""
##########
# Offset #
##########
"""
The offset that is input into the model is of the star wrt the meta data raNominal and decNominal. In this cube, it is not the star
that is off-centred, but rather the telescope was set to place the star in an off-set position in the FoV. The meta data have not changed
- they are still the SIMBAD position - so for this test we first have to change the meta data to pretend that it is the star that is offset.
The meta data have to be the position of the centre of the star, and then we need to determine the offset of the
star from that centre.
"""
ra = MEAN(slicedRebinnedCubes_off_NoESC.get(0).getRa()[:,2,2])
dec = MEAN(slicedRebinnedCubes_off_NoESC.get(0).getDec()[:,2,2])
# Do this to the slicedRebinnedCubes and each cube therein
slicedRebinnedCubes_off_NoESC.meta["raNominal"].value = ra
slicedRebinnedCubes_off_NoESC.meta["decNominal"].value = dec
slicedRebinnedCubes_off_NoESC.get(0).meta["raNominal"].value = ra
slicedRebinnedCubes_off_NoESC.get(0).meta["decNominal"].value = dec
"""
What is the offset of the star from the central spaxel (which = raNominal and decNomial)
By making an interpolated cube of the rebinned cube,
icube = specInterpolate(slicedRebinnedCubes_off_NoESC)
icube = icube.get(0)
it is possible to get a better idea of its coordinates - open cube in the Standard Cube Viewer where you
can see the coordinates and even measure offsets with the cursor (read the Data Analysis Guide chapter on
images for instructions). Or you can just estimate the number of spaxels the source is from the centre (spaxel 10,10)
and, knowing that each spaxel is 3"in size, you get the offset value. The sky direction is indicated in the small
image view at the top right of the Cube Viewer. Another possibility is to make an image from the cube:
map = ntegrateSpectralMap(cube=icube, startArray=Double1d([66.356]), endArray=Double1d([66.569]))
map = map.getSimpleImage(0)
and you can then inspect that with the Image Viewer
A bit of trial and error is called for here - guess the offset, create a simulated rebinned cube, then
run specInterpolate to create a simulated "icube", and check the result against the icube created from the data
(from slicedRebinnedCubes_NoESC)
(An offset to the N is + in Dec, offset to W is - in RA)
Doing this, we estimate an offset of 6" to the East.
"""
sourceSize = 0.2
offRa = +6.0
offDec = 0.0
contFlux = 160.0/(2*Math.PI*(sourceSize/(2*SQRT(2*LOG(2))))**2)
contSlope = 0.0
lineAmplitude = 442.0/(2*Math.PI*(sourceSize/(2*SQRT(2*LOG(2))))**2)
lineCentre = 66.45
lineWidth = 0.0078
sourceModel = MyGaussSpectral(sourceSize,offRa,offDec,contFlux,contSlope,lineAmplitude,lineCentre,lineWidth)
d=Display(imageModelInPacsSpecFromModel(sourceModel, 109, 0.5, lineCentre), title='continuum layer Gaussian model')
simRebinnedCubes_off = pacsSpecFromModel(slicedRebinnedCubes_off_NoESC,sourceModel,calTree=calTree)
simInterpCubes = specInterpolate(simRebinnedCubes_off)
simicube = simInterpCubes.get(0)
"""
Comparing the shape of the source in simicube to icube, we see that the match is not bad but the simulated shape does not
equal the observed shape exactly - the observed star is a bit lopsided. However, for the purpose of this test we will accept this result.
"""
##############
# Simulation #
##############
"""
How does the spectrum of the simulation compare to the observation and to the input?
"""
p1 = PlotXY(titleText='Simulated vs Observed')
spec1 = extractRegionSpectrum(cube=simRebinnedCubes_off.get(0),arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
spec2 = extractRegionSpectrum(cube=slicedRebinnedCubes_off_NoESC.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
p1.addLayer(LayerXY(spec1.wave,spec1.flux,name="simulated cube"))
p1.addLayer(LayerXY(spec2.wave,spec2.flux,name="observed (no esc)"))
p1.addLayer(LayerXY(specPSC.wave, specPSC.flux, name="PSC spectrum"))
p1.legend.visible=1
p1.xaxis.titleText="Wavelength ($\micro$m)"
p1.yaxis.titleText="Jy"
"""
The correction - our unique ESC - is the difference between the flux in the simulated cube and the flux in the input spectrum.
spec1 = 7.55 (simulated)
spec2 = 7.09 (observed)
specPSC = 8.64 (true, input)
simulated/true = 0.88 -> observed/0.88 = 8.06
The observation and simulation match fairly well: the corrected results is 0.93 of the correct result, not a perfect
match but not bad, especially given the fact that our model does not 100% simulate the real source. In addition, with this amount
of offset, some of the source falls off the FoV and therefore is not included in the sum created from the observed cube or from
the simulated cube. However, our input flux is the true flux of the source which means that this input assumes a fully-centred and
fully-accounted for source.
"""
#################################################################################################################################
# Using a standard spectrum
#################################################################################################################################
"""
Next we demonstrate that this FM works even if you do not know the true spectrum of the source -- which will be the case for most
of you using this script. We will run the same modelling on the off-centred observation of previously, but as input
we will use a standard spectrum of line with a peak of 1 Jy on a continuum of 0.
"""
slicedRebinnedCubes_off = selectSlices(slicedRebinnedCubes_init, sliceNumber=[2])
slicedRebinnedCubes_off_NoESC = undoExtendedSourceCorrectionSliced(slicedRebinnedCubes_off, calTree=calTree)
# the observed spectrum
specSUM_off_noesc = extractRegionSpectrum(cube=slicedRebinnedCubes_off_NoESC.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
# integrated flux of this spectrum's line is 7.09
# adjust the coordinates as we did before
ra = MEAN(slicedRebinnedCubes_off_NoESC.get(0).getRa()[:,2,2])
dec = MEAN(slicedRebinnedCubes_off_NoESC.get(0).getDec()[:,2,2])
slicedRebinnedCubes_off_NoESC.meta["raNominal"].value = ra
slicedRebinnedCubes_off_NoESC.meta["decNominal"].value = dec
slicedRebinnedCubes_off_NoESC.get(0).meta["raNominal"].value = ra
slicedRebinnedCubes_off_NoESC.get(0).meta["decNominal"].value = dec
# create the model and the simulated cube
sourceSize = 0.2
offRa = +6.0
offDec = 0.0
contFlux = 0.0
contSlope = 0.0
lineAmplitude = 1.0/(2*Math.PI*(sourceSize/(2*SQRT(2*LOG(2))))**2)
lineCentre = 66.45
lineWidth = 0.0078
sourceModel = MyGaussSpectral(sourceSize,offRa,offDec,contFlux,contSlope,lineAmplitude,lineCentre,lineWidth)
simRebinnedCubes_off = pacsSpecFromModel(slicedRebinnedCubes_off_NoESC,sourceModel,calTree=calTree)
"""
the total flux for our input line is
1 * lineWidth * SQRT(2*java.lang.Math.PI) = 0.0196 Jy*microns
For the simulated cube we will sum up the FoV and measure line in the spectrum fitter gui, reading off the peak and width values
This is so we make a like-wise comparison to the flux of the input spectrum
"""
specSim = extractRegionSpectrum(cube=simRebinnedCubes_off.get(0),arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
# The flux (using the equation of a Gaussian)
# 0.864*0.00798*SQRT(2*java.lang.Math.PI) = 0.0173
"""
The difference, the ESC, is 0.0173/0.0196 = 0.88
This correction should be applied to the flux of the extracted source taken from the no ESC cube.
So now going back to the data and measuring the line flux of from the observed cube slicedRebinnedCubes_off_NoESC -- this being
7.09 as measured a few lines above -- the corrected spectrum here is 8.06. We know the result should be 8.64: the ratio
is 0.94, the same as the previous test.
"""
#################################################################################################################################
# Testing on a mapping observation
#################################################################################################################################
"""
The final demonstration is to show how to do the same thing with a mapping observation. For a mapping observation of an extended
source, the cube to use and to simulate will be one of the two mosaic cubes - projected or interpolated. For this example,
rather than taking a mapping observation of an extended source, we will use the same pointed observation of R Dor but treat is
as a mapping one. To this end we will create a mosaic cube from a single rebinned cube only. Those working with mapping observation
would create their mosaic cubes from all the rebinned cubes taken in the raster. But other than this difference, the procedure
outlined here is the same as you would follow: take your Level 2/2.5 rebinned cubes, remove the ESC, create a projected or
interpolated cube, and so on ...
We will create a projected _and_ interpolated cube, to demonstrate that there are differences (which all arise from the fact that
these two algorithms deal with the data differently). We will work with the raster position where R Dor is centred.
"""
# Start by creating a mosaic cube with the pipeline ESC removed
# In our example there is only one rebinned cube in slicedRebinnedCubes, but for you there will be more
slicedRebinnedCubes = selectSlices(slicedRebinnedCubes_init, rasterLine=[3], rasterCol=[3]) # centred source (from test 1)
slicedRebinnedCubes_NoESC = undoExtendedSourceCorrectionSliced(slicedRebinnedCubes, calTree=calTree)
# interpolated 3" spaxels
slicedInterpolatedCubes_NoESC = specInterpolate(slicedRebinnedCubes_NoESC)
specSUM_map_interp = extractRegionSpectrum(cube=slicedInterpolatedCubes_NoESC.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
# projected 3"spaxels
slicedProjectedCubes_NoESC = specProject(slicedRebinnedCubes_NoESC)
specSUM_map_proj = extractRegionSpectrum(cube=slicedProjectedCubes_NoESC.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
"""
specSUM_map_proj|interp is the observed spectrum of our source from our pretend mapping observation, with the
pipeline ESC removed, and taken from the projected and interpolated cube.
The two are not exactly the same:
interpolated: 7.68 Jy*microns
projected : 8.23 Jy*microns
The true flux is 8.64 Jy*micons as measured from the point-source corrected spectrum (test 1)
Now we need to do the forward modelling on the source. We know the offRa and offDec from the work done in
test 1, and we will here also work with a standard spectrum as the input to compute the correction
"""
sourceSize = 0.2
offRa = 0.0
offDec = -1.4
contFlux = 0.0/(2*Math.PI*(sourceSize/(2*SQRT(2*LOG(2))))**2)
contSlope = 0.0
lineAmplitude = 1.0/(2*Math.PI*(sourceSize/(2*SQRT(2*LOG(2))))**2)
lineCentre = 66.45
lineWidth = 0.0078
"""
As computed above, the flux of a line of this peak and width is
1 * lineWidth * SQRT(2*java.lang.Math.PI) = 0.0196 Jy*microns
"""
sourceModel = MyGaussSpectral(sourceSize,offRa,offDec,contFlux,contSlope,lineAmplitude,lineCentre,lineWidth)
simRebinnedCubes = pacsSpecFromModel(slicedRebinnedCubes_NoESC,sourceModel,calTree=calTree)
# remember, the task only produces rebinned cubes and only works on rebinned cubes
# and as we are working with an mapping observation, we now need to turn the simulated rebinned cube(s) into a mosaic cube
slicedSimInterpolatedCubes = specInterpolate(simRebinnedCubes)
slicedSimProjectedCubes = specProject(simRebinnedCubes)
# extract the spectrum from the simulated cube
specSUM_simMap_interp = extractRegionSpectrum(cube=slicedSimInterpolatedCubes.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
specSUM_simMap_proj = extractRegionSpectrum(cube=slicedSimProjectedCubes.get(0), arithmetics=herschel.ia.toolbox.cube.ExtractRegionSpectrumTask.Arithmetics.SUM)
"""
Now compute fluxes from the equation of a Gaussian
peak and width area under Gaussian
specSUM_simMap_interp: 0.872 0.00798 0.0174
specSUM_simMap_proj : 0.926 0.00798 0.0185
input : 1 0.0078 0.0195
The unique ESC is the ratio simulated/input
interpolated : 0.89
projected : 0.95
Applying the correction to the observed spectrum
interpolated : 7.68/0.89 = 8.63
projected : 8.23/0.95 = 8.66
True flux (point-source corrected value): 8.64
--> the correction works and you can also see that it does depend on which mosaic task is used
"""
#################################################################################################################################
# 2D correction map
#################################################################################################################################
"""
The correction you calculate here is to the summed spectrum of a source.
If you want to know what the distribution of the correction is over your field - that is, where the correction is higher
and where it is lower - you can do this by modelling a flat source. This will give you the illumination pattern for your
particular mapping pattern. The pattern on the resulting cube shows the relative distribution of the correction, but not the
absolute. So it can show you where in the field the source is more and where it is less affected.
"""
# for a single pointing
slicedRebinnedCubes = selectSlices(slicedRebinnedCubes_init, rasterLine=[3], rasterCol=[3])
lineCentre = 66.45
# This is just to set the wavelength range, the actual value does not matter, since what we care about
# is the relative spatial distribution of the resulting cube
contFlux = 1.
contSlope = 0.0
sourceModel = MyFlatSourceSpectral(contFlux,contSlope,lineCentre)
d=Display(imageModelInPacsSpecFromModel(sourceModel,109, 0.5, lineCentre), title='continuum layer flat model')
simRebinnedCubes_point = pacsSpecFromModel(slicedRebinnedCubes,sourceModel,calTree=calTree)
simProjectedCubes_point = specProject(simRebinnedCubes_point, outputPixelsize=0.5) # making very small spaxels makes it easier to see the result
# For a mapping observation
slicedRebinnedCubes = selectSlices(slicedRebinnedCubes_init, rasterLine=[2,3,4], rasterCol=[2,3,4])
contFlux = 1.0
contSlope = 0.0
sourceModel = MyFlatSourceSpectral(contFlux,contSlope,lineCentre)
simRebinnedCubes_map = pacsSpecFromModel(slicedRebinnedCubes,sourceModel,calTree=calTree)
simProjectedCubes_map = specProject(simRebinnedCubes_map, outputPixelsize=0.5)
"""
Inspecting the spatial distribution of the rebinned and/or projected cubes and comparing to your model, shows you where in the
field the correction is higher and where it is lower
"""
cube = simProjectedCubes_map.get(0)
map = integrateSpectralMap(cube=cube, startArray=Double1d([66.0]), endArray=Double1d([67.0])
map = map.getSimpleImage(0)
image = map.getImage()
mean = StatWithNaN(MEAN)(image)
image = image/mean
map.setImage(image)
# now you have an image scaled to the mean which you can compare to your cubes
#### EXTRAS
# plot the pipeline ESC
esc = calTree.refs["spectrometer"].product.refs["extendedSourceLoss"].product
PlotXY(esc["wavelength"].data,esc["fraction"].data)