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())