rcube = fitsReader(file = '/Users/kexter/PacsRebinnedCubeHCSS14987.fits')
#
# Or if you want to test on a projectedCube, just to make sure it works you can make one from
# your rebinned cube.
#
scube=specProject(rcube,outputPixelsize=6.0) # just for testing, making large spaxels so the script runs not too slowely
# Chose the cube and set the parameters
cube=rcube
ModelSpecList=[] # a list to hold the model fit spectra for each line you run the loop on,
ParamsList=[] # and a list to hold the parameters of the fits
gaussP1=[157.7,145.5] # the lines you want to fit
wide=4 # the width about the line centre that covers the line+local continuum; from this range the line peak and continuum level are estimated
gaussP2=[0.1,0.1] # the approximate width of the lines
porder = 1 # order of the poly I want to fit to the continuum
# need to change the line "M1 = sf.addModel("poly",[porder],[cont,0])"
# if increase the order from 1, as then there will be more parameters in the final []
c = herschel.share.unit.Constant.SPEED_OF_LIGHT.value
#
# Weights. For masking IN or OUT regions from the fitting
# You can hardwire the weights (i), ignore the whole thing, or set them to be some wavelength range around the line (ii)
# Chose option (i), (ii) or leave them commented-out. Which you do will change some of the script in the loop, so read the rest of the
# script to find out where to comment in/out the dependent lines
# (i)
weights = [[157.1,158.3],[144.9,146.1]]
# (ii)
#weights = [0.6,0.6]
#
doplot = 1 # (1) yes, I do want to see a PlotXY of the fits as they are done, or not (0)
convunits = 0 # (1) convert to W/m2 for integrated line flux, or not (0)
for i in range(len(gaussP1)):
cnt=0 # just for the plotting
sz=cube.dimensions[0] # the length of the wavelength array of the cube
Totmodel=Double4d(cube.dimensions[1],cube.dimensions[2],3,sz) # to hold 3 model fit spectra (as 5x5 "cubes") [wavelength, M1+M2, M2: see later]
Params=Double4d(cube.dimensions[1],cube.dimensions[2],2,3) # to hold the parameters of the models defined (2): M1 and M2
tit=str(gaussP1[i])
if (doplot==1):
p=PlotXY(titleText=tit)
print "Doing line",gaussP1[i]
print " spaxel"
for x in range(cube.dimensions[1]):
for y in range(cube.dimensions[2]):
print x,y
#
# A catch for NaN spectra
# Note: if a spectrum is completely NaN it will be because it is an edge spectrum from a projected cube; in these cubes
# many spaxels at the edge of the FoV are actually empty (there was no observation in that piece of sky) but the
# spaxel is present to help make the FoV of the projected cube a regular, rectangular grid. For these spaxels we do
# not want to try to fit the spectrum, so we will skip those. For these spaxels, we will fill the Totmodel and Params with
# NaN results, so you know when looking at the output that these are not "bad fits", but are "no fits"
# Adding in a catch for case where spectra at edges are all 0 rather than NaN
#
flx = cube.getFlux()[:,x,y]
idx1= flx.where(IS_FINITE(flx))
idx = flx.where(IS_NAN(flx)) # find all array indices corresponding to flux values of NaN
if ((len(flx) == len(flx[idx])) | (SUM(flx[idx1]) == 0.0)):
Totmodel[x,y,:,:]=Float.NaN
Params[x,y,:,:]=Float.NaN
else:
#
# For the fitting it helps to have a decent guess at the peak flux of your line.
# If your continuum is significant, then this should be included in the first guess of the Gaussian's flux, as the fitter
# needs to know the value of the peak as it would be above a continuum of 0.
# To do both we will simply look for the max (peak) and min (continuum) values in a range that is
# +/- wide of the fwhm. If you have fat lines or a spectrum crowded with lines, this could give the wrong first-guess
# result. You can change this value of course
#
wv=cube.getWave(x,y)
flx=cube.getFlux(x,y)
mx=gaussP1[i]+(gaussP2[i]*wide)
mn=gaussP1[i]-(gaussP2[i]*wide)
idx=wv.where((wv < mx) & (wv > mn) & IS_FINITE(flx)) # locating the array index of the appropriate wavelengths
pk=MAX(flx[idx]) # calculating the mean flux from the flux array values in this index range
cont=MIN(flx[idx])
#
# Initialise and set up the parameters of the fitting
#
sf=SpectrumFitter(cube,x,y,False)
sf.useFitter(1)
M1 = sf.addModel("poly",[porder],[cont,0]) # model, order, [parameter_0, parameter_1,...]
M2 = sf.addModel("gauss",[pk-cont,gaussP1[i],gaussP2[i]]) # model, [P0, P1, P2] being amplitude, wave, fwhm for the ith spectral line
#
# Set the weights. Comment in the line that is appropriate to the choice you made above
# (i)
sf.setMask(weights[i][0],weights[i][1],1.0)
# (ii)
#sf.setMask(gaussP1[i]-weights[i],gaussP1[i]+weights[i],1.0)
#
# Do the fitting
# If for some reason the fitting does not work, we need to put in a catch so that the
# entire for loop does not die
#
try:
sf.doGlobalFit()
sf.residual()
#
# Display the results
#
if (doplot==1):
tit="spaxel "+str(x)+","+str(y)
p.addLayer(LayerXY(wv,flx,line=1,color=java.awt.Color.red,xtitle="wavelength",ytitle=tit),x,y)
p.addLayer(LayerXY(wv,M1.synthetic+M2.synthetic,line=1,color=java.awt.Color.blue),x,y)
p.addLayer(LayerXY(wv,sf.getResidual(),line=1,color=java.awt.Color.green),x,y)
p[cnt].setXrange([gaussP1[i]-(gaussP2[i]*wide*2),gaussP1[i]+(gaussP2[i]*wide*2)])
p[cnt+1].setXrange([gaussP1[i]-(gaussP2[i]*wide*2),gaussP1[i]+(gaussP2[i]*wide*2)])
p[cnt+2].setXrange([gaussP1[i]-(gaussP2[i]*wide*2),gaussP1[i]+(gaussP2[i]*wide*2)])
cnt+=3 # a counter so the various spectra are plotted to the line you are fitting only
#
# Saving the models and parameters so they are accessible after you have run this script
# Parameters I am chosing to save here are: Gauss integrated (not peak!) flux, wavelength, fwhm
# Model I am saving are the total model (Gauss + Poly) and also just the Gauss. The first is
# so one can compare the model fit to the spectrum later, the second is for making line-flux maps
#
Totmodel[x,y,0,:]=wv
Totmodel[x,y,1,:]=M1.synthetic+M2.synthetic # Gauss + Poly
Totmodel[x,y,2,:]=M2.synthetic # Gauss
#
# Get the global model results-they are in the order you defined them in (M1 then M2)
# Extract the wave, fwhm and their errors for M2 (Gauss) from this. If you want to, you
# can also extract the parameters of the poly--here we do not
# The integrated flux of the Gauss is more astronomically relevant than the peak flux
# The peak flux and its error is gotten in a similar way to the fwhm and wave, but the
# integrated flux is taken from the models directly (hence a different syntax for getting
# them)
# NOTE: the SpectrumFitter "integrated" flux is calculated via a trapezium--integration
# under your curve, using the datapoint of your spectrum. Hence, if your spectrum
# is sparsely sampled, the flux intflx=M2.getIntegral() returns will be a slight
# underestimate.
# You could alternatively numberically calculate the integrated flux. For a Gaussian
# the equation would be, instead of M2.getIntegral(),
# intflx = amp * sig * SQRT(2.* java.lang.Math.PI) * hzPerMicron * 1.e-26 # in W/m2
# intflx = amp * sig * SQRT(2.* java.lang.Math.PI) # in Jy
#
sey=sf.getGlobalResult()[1] # get the results of the Gauss
amp = sey["Parameters"].data[0] # peak flux of Gauss
centre=sey["Parameters"].data[1] # wavelength of Gauss
sig=sey["Parameters"].data[2] # fwhm of the Gauss
amp_err= sey["StdDev"].data[0] # and their respective errors
centre_err=sey["StdDev"].data[1]
sig_err=sey["StdDev"].data[2]
intflx=M2.getIntegral() # the integrated flux for the Gauss, or...
if (convunits == 1):
hzPerMicron = c / (centre*1.e-6) / (centre*1.e-6) / 1.e6
intflx = amp * sig * SQRT(2.* java.lang.Math.PI) * hzPerMicron * 1.e-26 # in W/m2
intflx_err = SQRT( (amp_err/amp)**2+(sig_err/sig)**2 ) * intflx
Params[x,y,0,0]=intflx
Params[x,y,1,0]=intflx_err
Params[x,y,0,1]=centre
Params[x,y,1,1]=centre_err
Params[x,y,0,2]=sig
Params[x,y,1,2]=sig_err
except:
# if the fitting fails, fill all the results with NaN
Totmodel[x,y,:,:]=Float.NaN
Params[x,y,:,:]=Float.NaN
ModelSpecList.append(Totmodel)
ParamsList.append(Params)
# Inspect the results..........
#
# To take the model fits out of the array they are held in, to e.g. plot then with PlotXY as used
# within the for loop:
x=2
y=2
Totmodel=ModelSpecList[0] # first line you worked on
totfit=Totmodel[x,y,1,:]
wave=Totmodel[x,y,0,:]
gaussfit=Totmodel[x,y,2,:]
p=PlotXY(wave,totfit,line=1)
spec=cube.getFlux()[:,x,y]
p.addLayer(LayerXY(wave,spec))
# To plot the flux vs error for all spaxels:
Params=ParamsList[0] # first line you worked on
flx=Params[:,:,0,0]
err=Params[:,:,1,0]
PlotXY(RESHAPE(flx),RESHAPE(err/flx),line=0)
# To image the fluxes (you can easily modify this to make a velocity map, ...)
# NOTE: if you were fitting on a projectedCube (created by specProject) made from a dithered observation
# set, taken to allow you to sample the PACS spectral PSF with at least a Nyquist value, then the
# positions of the spaxels in the cube are a true representation of the sky footprint of the instrument.
# Hence, the simple map we show you how to make here, from the results of your fitting will also
# have the correct sky footprint, i.e. the correct WCS.
#
# However, if you are working on a single pointed PacsRebinnedCube (as created by specAddNodCubes,
# for example), the beam is not properly sampled and, as we have said in many places in the PACS
# documentation, you should not run the cube through specProject. Instead, your science measurements
# should be made on the PacsRebinnedCube itself.
# Now, the WCS of these cubes is *not* correct, it is only approximately correct. This is
# because the sky footprint of the PacsProjectedCube is not regular (see the PDRG or remember the
# lectures from earlier in this workshop to learn why). Hence, a (simple) map made from your fits
# to a PacProjectedCube will not have the correct sky footprint/WCS. The spaxels of the cube do
# have in them their correct RA and Dec coordinates, and hence you can, with some scripting,
# create a map from your fits. However the maps made as we show do not have the correct sky
# footprint. (The WCS of the PacsRebinnedCube is an approximation based on the correct RA, Dec
# information in each spaxel but roughly interpreted into a regular grid; this was done to allow
# the cube to be opened in the various display and spectral fitting tasks of HIPE.
line=0 # first line
Params=ParamsList[line] # first line you worked on
map1=SimpleImage(image=Params[:,:,0,0])
map1.setWcs(cube.getWcs())