#Demonstration on how combine PACS and SPIRE spectra of a semi-extended source to compare the total source flux density applying: # 1._The semi-extended source correction to a PACS SED spectrum of a semi-extended source # 2._SECT tool to the SPIRE spectrum (large Gaussian reference beam) # #Spectrometer Workshop: PACS #23 April 2014 #Created by E. Puga # # 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 . # """ For HIPE Track 12 and onwards This script shows you how to push together PACS and SPIRE SED spectra into a single spectral product, i.e. to "combine" PACS and SPIRE spectra. This is intended for point sources, and furthermore, for PACS it is intended for spectra covering the full spectral range ("SED" mode), although it will work for smaller r?ange/line scans also. For SPIRE the entire spectral range is always covered. You begin with a single spectrum, point-source loss corrected where necessary, of your source, for all the bands you have. You will end with a Spectrum1d product that contains all the individual PACS and SPIRE spectra, not mathematically merged but rather simply contained in one product. The basic operation is simple, and there are some extras you can do. All you need to do this is contained in this script, but for the SPIRE part you will need to be working on a SPIRE buildof HIPE, and for the PACS part, on a PACS build. The final part can be done on either. WHAT THIS SCRIPT DOES: 1- Takes the PACS and SPIRE point source (i.e. single) spectra--normally this would be a single spectrum per band and for an object that is, obviously, a point source. This will usually mean 4 PACS spectra and 2 SPIRE spectra. For PACS this means running the bit at the end of the pipeline scripts, where a spectrum is extracted from the cube and point source loss corrected. For SPIRE you work from (a properly reduced) ObservationContext, grabbing the appropriate Level 2 product (see the SPIRE data reduction guide to learn more). **This part of the job has to be done on a PACS HIPE build for PACS and a SPIRE build for SPIRE** 2- The PACS and the SPIRE spectra are separately pushed into an "allPACS" Spectrum1d and an "allSPIRE" SimpleSpectrum **this having to be done on a SPIRE build for SPIRE, but can be done on either build for PACS** (after that you can work on any type of HIPE build). This requires the use of functions included at the top of this script (prepPacsSpectra and prepSireSpectra). Load those functions into HIPE memory by clicking the single green arrow when lined up on the "def..." line of each. 3- A final function (mergePacsNSpireSpec) then pushes the spectra into a single product--a Spectrum1d--with different "segment" numbers for the different spectra. At the same time it converts the flux and wavelength units to the same value for all the spectra. WHAT THIS SCRIPT DOES NOT DO: o- It does not do the PACS point source correction: that has to be done to the separate spectra yourself (see the demo below) o- It does not merge the spectra into a single one--there is no flux scaling to match the PACS and the SPIRE, and no wavelength regridding WHAT IS THE POINT OF THIS SCRIPT? o- This script is provided for ease of viewing an SED of an object over the full PACS to SPIRE range. By having all the spectra in a single product you can open it with the Spectrum Explorer and view the separate segments--the separate spectra--individually and all together o- We do not fully merge the spectra because (i) there may be offsets between the continuum of PACS and SPIRE, and to adjust for this the astronomer needs to decide what is best to do, and, more importantly (ii) the spectral resolution and spectral sampling differ so much across the PACS to SPIRE range that a simple merging does not have a strong scientific case WHAT IS NECESSARY TO DO BEFORE YOU RUN THIS SCRIPT: o- PACS: to have pipeline-reduced your data (or to use SPG-reduced data). At the end of all the PACS pipeline scripts you are offered a task to extract a single spectrum of your point source: you take the spectrum from a single spaxel (that with most of the source in it, usually the central spaxel) and apply the point source flux correction. PACS provides, in HIPE Track 9 and beyond, two tasks that do this and it is always documented in the pipeline scripts and the PDRG. The recommended task is "extractCentralSpectrum", but if your source is not located in the central spaxel you will need to use "extractSpaxelSpectrum" instead: this task is explained in Scripts->Pacs Useful Scripts->"Spectroscopy: Point source loss correction (any spaxel)". o- SPIRE: you can use data directly from the Herschel Science Archive, provided the observation was processed with SPG 8 or above. There are two possible ways of extracting SPIRE point source spectra from the two wavelength arrays, one for sparse mode and one for mapping observations; however that for mapping has not been completed so you can only work with sparse mode observations. FURTHER READING: If you have a random spectrum, e.g. from a single spaxel or a combined spectrum taken from parts of a cube, or some other Observatory's spectrum, and you want to turn that into a Spectrum1d and at the same time combine several of them, the steps you need to follow are easier. Much of the bulk in the prepPacsSpec script (especially) is concerned with creating Meta data in the output product, and error catches; making the output, merged, Spectrum1d is conceptually very easy: create the flux, wavelength, weights, flags, and segment arrays from the input spectra, and push them into a Spectrum1d; add units; ...basta. How to create and work with Spectrum1d (as well as Spectrum2d and cubes and SimpleSpectrum) products can be found in the Scripting Guide of HIPE, currently Chapter 3 is about spectra and cubes. AUTHORS Katrina Exter (PACS); Ivan Valtchanov (SPIRE) MODIFICATION HISTORY: v1 - May 2012 v2 - May 2013: simplifying the SPIRE part, now only ine def prepSpireSpec() is needed, its output is a SimpleSpectrum. v3 - May 2013: and then adjusted PACS part to account for changes in the output of extractCentralSpectrum v4 - Sept 2013: and again for the same reason, change the example target, clean up the meta data a bit, and add weights to the PACS part as these now exist in the input spectra """ # imports from herschel.share.unit import Angle import re from herschel.pacs.spg.spec.beams import ExtendedSource from herschel.ia.gui.plot.renderer.PlotLegendEngine import Position from math import * ##################### ### THE FUNCTIONS ### ##################### # Load these into memory (place your cursor to the left of "def..." line for each function, # and click the single green arrow at the top of HIPE def ls(directory="./",pattern="",case=1,verbose=1): """ls(directory="./",pattern="",case=1) Returns the files in "directory" (relative path allowed) with filename containing "pattern" (incl. wildcards) case : 0/1 for case (in)sensitive listing (prefer case = 1 for pattern filtering) """ import sys,os,fnmatch if case: allfiles = SORT(String1d(os.listdir(directory))) else: allfiles = String1d([file.lower() for file in SORT(String1d(os.listdir(directory)))]) files = allfiles if pattern: files = [file for file in allfiles if fnmatch.fnmatch(file, pattern)] if verbose: for file in files: print file return files ############### #### PACS ##### ############### #Gaussian Model definition class MySource(ExtendedSource): sigma = 0.0 offRa = 0.0 offDec = 0.0 def __init__(self,diameter, offRa, offDec): self.sigma = diameter/2.3548 self.offRa = offRa self.offDec = offDec def getFlux(self,relRa,relDec): flux = exp(-(0.5*((relRa-self.offRa)/self.sigma)**2+0.5*((relDec-self.offDec)/self.sigma)**2)) return flux def renderGauss(self, size=109, dpix=0.5): self.image=Double2d(size, size) for y in range( size ): for x in range( size ): xim = (-dpix)*(x-(size/2+0.5)) yim = (dpix)*(y-(size/2+0.5)) self.image.set(y,x, source.getFlux(xim, yim)) return self.image def getSourceModelImage(self, size=109, dpix=0.5): return self.renderGauss(size, dpix) def prepPacsSpec(spectra=[], goods=None): """ For HIPE Track 12 and beyond A function to combine several PACS spectra into a single product. "Combine" means pushing the spectra into a single Spectrum1d object--it does not merge the data mathematically. Each separate spectrum placed in the Spectrum1d has a unique segment number (they increment from 0) and hence they can be distinguished from each other. In this function you are also offered the option of "flagging" the input spectra. The flag dataset takes on value of 0 in the "good" spectral ranges and 1 elsewhere. This is merely a cosmetic move--in case you want to record this information within the spectrum itself. The data themselves are not changed. >Copy the flags from the input product, default behaviour and you need do nothing >Set all flags to 0 by setting goods=[0] NOTE THIS CHANGE FROM PREVIOUS VERSIONS >If you want to set the hardwired flags (see below), set goods=[1] >If you want to set your own flags, then you need to define one list of left and right limit(s) per spectrum, and hold these together in a list that you input as "goods" (see the example given below) For the hardwired flags we have adopted spectral ranges that are appropriate for SEDs using the spectral leakage plots of the PACS Observers Manual (Sec. 4.8). So this is not just masking out the ends of the spectra where the response takes a dive, but also where there could be spectral leakage. NOTE: April 2012 the combined allPacs Spectrum1d (and the combined allSpire) has an error column in it. This has been added because SPIRE have an error column, and the two Spectrum1d need to have the same layout so they can be combined into an allPacsSpire. However, at present PACS do not propagate errors from the PacsRebinnedCube to the extracted and point source loss corrected specrum, so this column is filled with 0. NOTE: Sept 2013. The input point source spectrum from PACS now has a filled weights and flag dataset, and so the weights are now copied over to the output and the users is offered the option of copying over the flags also TASK CALL allPacsSpec = prepPacsSpec(spectra=spectra [,goods=mygoods]) INPUT spectra = python list of all the spectra, of class SimpleSpectrum, to be combined. They can be in any order. If you have 4 spectra (centralSpectrumB1,B2,R1,R2) you create this list via: spectra=[] spectra.append(centralSpectrumB1) spectra.append(centralSpectrumB2) spectra.append(centralSpectrumR1) spectra.append(centralSpectrumR2) If you have only one, then only one goes in the list. The wavelength unit for PACS data is micrometer and the flux unit is Jy/pixel; we copy over the flux and wavelength unit of the first spectrum in the input list to the output, just in case they are different to this. Note that the flux and wave unit of the have to be the same for all spectra in the output Spectrum1d. MANDATORY goods = python list of lists of good regions, that is areas to have a flag of 0. The inner lists are for each spectrum and the outer list simply holds them together so they can be set via a single parameter in this task. NOTE THAT THIS HAVE CHANGED SINCE THE PREVIOUS VERSION: the default is now to use the input flags - The default behaviour is to adopt the flags from the input (you do not need to specify anything). - To set goods yourself: if you have 4 input spectra which are, in this order, band B2A, B2B, R1 and R1, then you can do the following: good1=[53,72.5] good2=[71,98] good3=[100,190] good4=[100,190] mygoods=[good1,good2,good3,good4] - These values are also the hardwired ones, and if you wish to use these you can just use: mygoods=[1]. - Or, set all flags to 0 with: mygoods=[0] If you define limits for one spectrum, you must for all (note: setting the limits to e.g. [40,240] is the same as saying, for any spectral band, that "the whole range is good"). The order of good1, good2,... should follow the order of the spectra in "spectra" above. Good should always have an even number of entries (left limit, right limit, left limit, right limit...), and there can be more than 1 couple defined. OPTIONAL OUTPUT: speco = Spectrum1d, where the different input spectra in the "spectra" list are distinguished by "segment" number MODIFICATION HISTORY: v1 - Apr 2011 v2 - Apr 2012: added error column to match SPIRE (it is filled with 0s) v3 - Apr 2013: changed ra and dec bit v4 - Sepc 2013: add new option for flags (and set to the default) and tidied up meta data AUTHOR: Katrina Exter (PACS, KU Leuven) following interactions with Ivan Valtchanov (SPIRE, ESAC) """ # error catches sey=str(spectra.__class__) sey=str(re.search("list",sey)) if (sey == "None"): print "Error: input spectra are not in a python list.\n" return if (len(spectra) ==0): print "Error: you have no any spectra in your input list" return if (goods != None): if ((goods[0] != 0)&(goods[0]!=1)): if (len(goods) != len(spectra)): print "Error: the length of good and spectra is not the same" return # Make the following a try/except because it is not crucial and some older data do not have # the object name in the meta data try: object=spectra[0].meta["object"].string except: object=None if (object): for i in range(1,len(spectra)): try: objecti=spectra[i].meta["object"].string if (objecti != object): print "Warning: The target name for spectrum",i,"in your list is not the same as that for the first spectrum in your list" # only a warning because the spelling could be different except: object=None if (object == None): object="Unknown" # set up the Spectrum1d. This needs to be filled with long arrays where the segment number changes # with spectrum number. Start the process by setting up the array, initially from the first spectrum # and then from rest, and then push them into the Spectrum1d for i in range(len(spectra)): if (i == 0): sey=str(spectra[i].class) sey1=str(re.search("SimpleSpectrum",sey)) sey2=str(re.search("Spectrum1d",sey)) poo1="none" poo2="none" if (sey1 != "None"): poo1="ss" if (sey2 != "None"): poo2="s1d" if ((sey1 == "None")&(sey2=="None")): print "Error: entry",i," is not a SimpleSpectrum or Spectrum1d.\n" return spec = spectra[i].copy() wav = spec.getWave() flux = spec.getFlux() weight = spec.getWeight() seg = Int1d(len(wav),i) err = SQRT(1/weight) meta = [] meta.append(spec["dataset"].meta) # units - take that of the first input spectrum and to be used below if (poo1 == "ss"): fluxunit = spec["dataset"].getFluxUnit() waveunit = spec["dataset"].getWaveUnit() else: fluxunit = spec.getFluxUnit() waveunit = spec.getWaveUnit() band=spectra[i].getMeta()["band"].string # set the flags if (goods == None): sey=None print "All wavelengths for all spectra get a copy of the input flags" print "...spectrum",i,"is of band",band flag = spec.getFlag() elif ((len(goods) == 1) & (goods[0]==1)): print "The hardwired flags will be set for all spectra" print "...spectrum",i,"is of band",band flag = Int1d(len(wav),0) flag=flagSet(band,goods[0],flag,wav) if (flag==None): return elif ((len(goods) == 1) & (goods[0]==0)): print "The flags are all set to 0" print "...spectrum",i,"is of band",band flag = Int1d(len(wav),0) else: print "Setting user-defined flags" print "...spectrum",i,"is of band",band flag = Int1d(len(wav),0) flag=flagSet(band,goods[i],flag,wav) if (flag==None): return else: sey=str(spectra[i].class) sey1=str(re.search("SimpleSpectrum",sey)) sey2=str(re.search("Spectrum1d",sey)) poo1="none" poo2="none" if (sey1 != "None"): poo1="ss" if (sey2 != "None"): poo2="s1d" if ((sey1 == "None")&(sey2=="None")): print "Error: entry",i," is not a SimpleSpectrum or Spectrum1d.\n" return # now fill the arrays with data from the other spectra: (.include method does not seem to work so have to do this manually) spec = spectra[i].copy() wav2 = spec.getWave() flux2 = spec.getFlux() weight2 = spec.getWeight() seg2 = Int1d(len(wav2),i) err2 = SQRT(1/weight2) band=spectra[i].getMeta()["band"].string # set the flags if (goods == None): sey=None print "...spectrum",i,"is of band",band flag2 = spec.getFlag() elif ((len(goods) == 1) & (goods[0]==1)): print "...spectrum",i,"is of band",band flag2 = Int1d(len(wav2),0) flag2=flagSet(band,goods[0],flag2,wav2) if (flag==None): return elif ((len(goods) == 1) & (goods[0]==0)): print "...spectrum",i,"is of band",band flag2 = Int1d(len(wav2),0) else: print "...spectrum",i,"is of band",band flag2 = Int1d(len(wav2),0) flag2=flagSet(band,goods[i],flag2,wav2) if (flag==None): return wav.append(wav2) flux.append(flux2) weight.append(weight2) flag.append(flag2) seg.append(seg2) err.append(err2) meta.append(spec["dataset"].meta) # fill the Spectrum1d speco = Spectrum1d(description="PACS Spectrum1d, containing different input spectra") speco.setWave(wav) speco.setWeight(weight) speco.setFlux(flux) speco.setFlag(flag) speco.setSegment(seg) # to add the error I have to do some massaging unfortunately speco.addColumn("error",Column(err,unit=speco.getFluxUnit())) speco.setWaveUnit(waveunit) speco.setFluxUnit(fluxunit) speco.setWaveDescription("Wavelength") speco.setFluxDescription("Flux") ras=[] # to calculate the mean ra and dec decs=[] for i in range(len(spectra)): meta1=meta[i] sey1="spaxelRa"+str(i) sey2="spaxelDec"+str(i) ras.append(meta1["spaxelRa"].double) decs.append(meta1["spaxelDec"].double) sey3="spaxelRa for segment "+str(i) sey4="spaxelDec for segment "+str(i) speco.meta.set(sey1,DoubleParameter(meta1["spaxelRa"].double)) speco.meta.set(sey2,DoubleParameter(meta1["spaxelDec"].double)) speco.setMetaUnit(sey1,Angle.DEGREES) speco.setMetaUnit(sey2,Angle.DEGREES) speco.setMetaDescription(sey1,sey3) speco.setMetaDescription(sey2,sey4) sey1="band"+str(i) sey3="band for segment "+str(i) speco.meta.set(sey1,StringParameter(spectra[i].getMeta()["band"].string)) speco.setMetaDescription(sey1,sey3) sey1="obsid"+str(i) sey3="obsid for segment "+str(i) speco.meta.set(sey1,LongParameter(spectra[i].getMeta()["obsid"].long)) speco.setMetaDescription(sey1,sey3) speco.meta.set("ra_mean",DoubleParameter(MEAN(ras))) speco.meta.set("dec_mean",DoubleParameter(MEAN(decs))) speco.setMetaDescription("ra_mean","mean ra of all segements") speco.setMetaDescription("dec_mean","mean dec of all segements") speco.setMetaUnit("ra_mean",Angle.DEGREES) speco.setMetaUnit("dec_mean",Angle.DEGREES) sey1="object" sey3="Target name" speco.meta.set(sey1,StringParameter(object)) speco.setMetaDescription(sey1,sey3) speco.meta.set("combiningComment",StringParameter("n/a")) speco.setMetaDescription("combiningComment","A combination of PACS spectra: see originals for full MetaData") return(speco) def flagSet(band,gd,flag,wve): flg=flag.copy() if (gd==1): if (band == "B3A"): sey=[56.0,70.0] elif (band == "B2A"): sey=[53.0,72.5] elif (band == "B2B"): sey=[71.0,98.0] else: sey=[100,190] idx=wve.where(wve < sey[0]) flg[idx]= 1.0 idx=wve.where(wve > sey[1]) flg[idx]= 1.0 else: sey=ABS(len(gd)/2.) - len(gd)/2 if ((len(gd) > 0) &(sey != 0.0)): print "Error: you have an odd number of ""good"" entries" return flg[:]= 1.0 for i in range(0,len(gd),2): idx=wve.where((wve >= gd[i])&(wve <= gd[i+1])) if (idx.length() == 0): print "Error: your limits",gd[i],gd[i+1],"are outside the range of the spectrum; start again" flg=None else: flg[idx]= 0.0 return(flg) ############### #### SPIRE #### ############### def prepSpireSpec(obsContext,toWavelength=False,apod=False): """ For HIPE Track 11 and beyond NAME: prepSpireSpec PURPOSE: combine the two FTS arrays into one Spectrum1d object, only the central detectors are used each FTS array will have a different segment number USAGE: result = prepSpireSpec(obsContext,apod="unapodized",toWavelength=False) INPUTS: obsContext - a SPIRE FTS observation context, sparse spatial sampling only toWavelength - if set to True then the output spectrum will be in wavelength in microns apod - if the apodized version of the final level-2 spectra is to be used (not recommended). OUTPUT: object of class SimpleSpectrum with the combined SSW and SLW spectra, each with a different segment number HISTORY: Created: Apr 2011, Ivan Valtchanov, HSC, ESAC, ESA, following interactions with Katrina Exter MODIFICATION HISTORY: v1.0 - Apr 2011, created v2.0 - Apr 2012, modifications: change waveUnit to a boolean toWavelength as the default wave unit is now frequency in GHz; also added an "error" column to the output as SPIRE works with errors and not weights v3.0 - Apr 2012, refactored to accept input of a list with SSW and SLW spectra v4.0 - Apr 2013, refactored to use spireProduct2SimpleSpectrum() """ # if (obsContext.meta["instMode"].string != "SOF1"): print "Only sparse mode observations can be used in prepSpireSpec() function" return None # res = obsContext.meta["commandedResolution"].string if res=="CR": res = "HR" apodName = "unapodized" if (apod): apodName = "apodized" # error catches sey=str(obsContext.__class__) sey=str(re.search("Context",sey)) if (sey == "None"): spec = spireProduct2SimpleSpectrum(obsContext, detector=["SLWC3", "SSWD4"], scan='0000') else: spss = obsContext.level2.getProduct("%s_%s_spectrum"%(res,apodName)) # hardcoded for the central detectors, can be changed to use any pair of detectors (with HIPE 11 on) # in previous versions level-2 spectra for sparse observations onply included the central detectors spec = spireProduct2SimpleSpectrum(spss, detector=["SLWC3", "SSWD4"]) # check if conversion of frequency to wavelength is needed if (toWavelength): # convert to wavelength (microns) spec = convertWavescale(ds=spec,to="micrometer", overwrite=False) waveDescription = "Wavelength" return spec ############## #### Both #### ############## def mergePacsNSpireSpec(pacs,spire,waveunit="micrometer"): """ For Track 12 and beyond Function to merge the PACS and SPIRE spectra, contained in a PACS Spectrum1d and a SPIRE SimpleSpectrum, into a final Spectrum1d object. This is aimed at the point-source corrected spectra of a single target observed with both instruments, and moreover with the full SED range observed by PACS. Each spectrum in the final product--in most cases this will be 4 PACS spectra (contained in the PACS input Spectrum1d) and 2 SPIRE spectra (contained in the SPIRE SimpleSpectrum input)-- will have a unique (and successive) segment number. In this task you read in the PACS and SPIRE spectra and specify what wavelength units you want the result to be in, the choices being: micrometer, GHz The flux unit is assumed to be Jansky (and the flux unit, which has to be the same for all spectra within a Spectrum1d, will be copied over from the PACS Spectrum1d *only*). INPUT pacs = PACS Spectrum1d MANDATORY spire = SPIRE SimpleSpectrum MANDATORY waveunit = desired wave unit ("micrometer", "GHz") DEFAULT micrometer OUTPUT Spectrum1d that is the combination of the PACS and SPIRE. MODIFICATION HISTORY: v1 - Apr 2011 v2 - Apr 2012 v3 - Apr 2013 (IV modifications) v4 - Sept 2013 added dealing with weights and errors and flags from PACS, and weights from SPIRE, and clean up Meta Data AUTHOR: Katrina Exter (PACS, KU Leuven) following interactions with Ivan Valtchanov (SPIRE, ESAC) """ # error catching sey=str(pacs.class) sey=str(re.search("Spectrum1d",sey)) if (not isinstance(pacs,herschel.ia.dataset.spectrum.Spectrum1d)): print "Error: the PACS spectrum is not a Spectrum1d.\n" return if (not isinstance(spire,herschel.ia.dataset.spectrum.SimpleSpectrum)): print "Error: the SPIRE spectrum is not a SimpleSpectrum.\n" return # make this a try/except as sometimes this does not work try: objectp=pacs.meta["object"].string except: objectp="Unknown" try: objects=spire.meta["object"].string except: objects="Unknown" if (objectp != objects): print "Warning: The target name is not the same for PACS and SPIRE:",objectp, "and", objects # ony a warning because the spelling could be different # set up the output spectrum tSpectrum=Spectrum1d() # get the input datasets but first make the wavelengths the same if (waveunit=="micrometer"): pacsnew = convertWavescale(ds=pacs,to="micrometer", overwrite=False) spirenew = convertWavescale(ds=spire,to="micrometer", overwrite=False) elif (waveunit =="GHz"): pacsnew = convertWavescale(ds=pacs,to="GHz", overwrite=False) spirenew = convertWavescale(ds=spire,to="GHz", overwrite=False) else: print "Error: allowed values for waveunit are only: ""micrometer"" and ""GHz""" return # PACS is a Spectrum1d and SPIRE a SimpleSpectrum, so slight differences in # syntax will be necessary pwave=pacsnew.getWave() pweight=pacsnew.getWeight() pflag=pacsnew.getFlag() pseg=pacsnew.getSegment() pflux=pacsnew.getFlux() perr= pacsnew["error"].data # sflux=spirenew.getFlux() swave=spirenew.getWave() sweight=spirenew.getWeight() sflag=spirenew.getFlag() serr= SQRT(1/sweight) sflag = spirenew["dataset"]["mask"].data # for SPIRE, increment from the PACS...so if the SPIRE segment numbers are not 0,1 there will be a jump sey=pseg[-1]+1 sseg=spire.getSegment() + sey # set up the total spectrum datasets of flux and wave twave = Double1d() tflux = Double1d() tweight = Double1d() tflag = Int1d() tseg = Int1d() terr = Double1d() twave.append(pwave) twave.append(swave) tflux.append(pflux) tflux.append(sflux) tweight.append(pweight) tweight.append(sweight) tflag.append(pflag) tflag.append(sflag) tseg.append(pseg) tseg.append(sseg) terr.append(perr) terr.append(serr) tSpectrum.setWave(twave) tSpectrum.setWeight(tweight) tSpectrum.setFlux(tflux) tSpectrum.setFlag(tflag) tSpectrum.setSegment(tseg) # to add the error I have to do some massaging tSpectrum.addColumn("error",Column(terr,unit=herschel.share.unit.FluxDensity.JANSKYS)) tSpectrum.setFluxUnit(herschel.share.unit.FluxDensity.JANSKYS) tSpectrum.setFluxDescription("Flux") if (waveunit == "micrometer"): tSpectrum.setWaveDescription("Wavelength") tSpectrum.setWaveUnit(herschel.share.unit.Length.MICROMETERS) else: tSpectrum.setWaveDescription("frequency") tSpectrum.setWaveUnit(herschel.share.unit.Frequency.GIGAHERTZ) # if (pacs.getFluxUnit().isEquivalent(herschel.share.unit.Unit.parse("Jy"))==False): print "Warning: the PACS input does not have flux units of Jy. Proceeding anyway, " print " and adopting the PACS flux unit *only* for the final result" tSpectrum.setFluxUnit(pacs.getFluxUnit()) tSpectrum.setFluxDescription(pacs.getFluxDescription()) tSpectrum.addColumn("error",Column(terr,pacs.getFluxUnit())) # error unit also change if (pacs.getFluxUnit().isEquivalent(herschel.share.unit.Unit.parse("Jy"))==False): print "Warning: the SPIRE input does not have flux units of Jy. Proceeding anyway," print " and adopting the PACS flux unit *only* for the final result" tSpectrum.setFluxUnit(pacs.getFluxUnit()) tSpectrum.setFluxDescription(pacs.getFluxDescription()) tSpectrum.addColumn("error",Column(terr,pacs.getFluxUnit())) pra=pacs.meta["ra_mean"].double pdec=pacs.meta["dec_mean"].double sra=spire.meta["raNominal"].double sdec=spire.meta["decNominal"].double tSpectrum.meta.set("combiningComment",StringParameter("n/a")) tSpectrum.setMetaDescription("combiningComment","A combination of PACS and SPIRE spectra: see originals for full MetaData") tSpectrum.meta.set("ra",DoubleParameter(MEAN([pra,sra]))) tSpectrum.meta.set("dec",DoubleParameter(MEAN([pdec,sdec]))) tSpectrum.setMetaUnit("ra",Angle.DEGREES) tSpectrum.setMetaUnit("dec",Angle.DEGREES) tSpectrum.setMetaDescription("ra","mean ra of PACS and SPIRE") tSpectrum.setMetaDescription("dec","mean dec of PACS and SPIRE") tSpectrum.meta.set("object",StringParameter(objectp)) tSpectrum.setMetaDescription("object","Target name (taken from PACS)") tSpectrum.meta.set("instrument",StringParameter("PACS+SPIRE")) # Now copying over the meta data from the input pacs and spire product, appending with # a PACS or SPIRE so those that are the same do not overwrite each other for key in pacs.meta.keySet(): key1="PACS:"+key tSpectrum.meta[key1] = pacs.meta[key].copy() for key in spire["dataset"].meta.keySet(): key1="SPIRE:"+key # for SPIRE copying over the meta data of the "dataset" rather than the product tSpectrum.meta[key1] = spire["dataset"].getMeta()[key].copy() return(tSpectrum) ########################## #### Now Run the code #### ########################## ##################### ### THE PACS PART ### ##################### path= '/home/epuga/PACS_ICC/PACS_TelBackNorm/calset64/' sourceName='NGC7027' model='Gaussian' diamP, offRa, offDec=12.0, 0., 0. #in arcsecs #Model Visualization source = MySource(diamP, offRa, offDec) #obsidList=[1342186968,1342186969]#NGC7027 According to Pierre, this is the real SED, but it was so early in time (OD182), that there is no orderSel # Chose the obsids obsids = {} obsids[1342186968] = "SEDA" obsids[1342186969] = "SEDB" #Each pool has been separated in red and blue ################################################################### calTree=getCalTree() #Load FinalCubes and concatenate them cubes=[] for obsid in obsids.keys(): for camera in ['blue','red']: pattern='*'+ str(obsid) +'*'+ camera+'*'+ '_FinalCubes' files=ls(directory=path,pattern = pattern, case=1, verbose=0) for file in files: print file cubes.append(readSliced(file, poolLocation=path, verbose=0)) cubes=concatenateSliced(cubes) spectra = [] spectraCorr = [] mygoods=[] mygoodsCorr=[] bandBorders = {"B2A":[55.,73.], "B3A":[55.,70.], "B2B":[70.,95.], "R1":[102.,190.], "":[50,230]} verbose=0 # for extractCentralSpectrum smoothing = 'wavelet' # or: smoothing = 'filter' nLowFreq = 4 for slice in range(len(cubes.refs)): rcube = cubes.get(slice) band = rcube.meta['band'].value c1,c9,c129 = extractCentralSpectrum(rcube, slice=slice, smoothing=smoothing, nLowFreq=nLowFreq, calTree=calTree, verbose=1) if (model == 'Gaussian'): # in the case of a Gaussian, it is the FWHM source = MySource(diamP, offRa, offDec)# in the case of a Gaussian, it is the FWHM central3x3=False e2p1 = specExtendedToPointCorrection(c1,calTree=calTree,userSource=source,central3x3=central3x3, extra=True) central3x3=True e2p9 = specExtendedToPointCorrection(c9,calTree=calTree,userSource=source,central3x3=central3x3, extra=True) else: #specExtendedToPointCorrection this assumes a disc top hat central3x3=False e2p1 = specExtendedToPointCorrection(c1,calTree=calTree,srcRadX=diam/2.,srcRadY=diam/2.,srcAngle=0.0,raOffset=0.0,decOffset=0.0,central3x3=central3x3,extra=True) central3x3=True e2p9 = specExtendedToPointCorrection(c9,calTree=calTree,srcRadX=diam/2.,srcRadY=diam/2.,srcAngle=0.0,raOffset=0.0,decOffset=0.0,central3x3=central3x3,extra=True) c1e2pCorr = divide(ds1=c1, ds2=e2p1, overwrite=False) c9e2pCorr = divide(ds1=c9, ds2=e2p9, overwrite=False) # Concatenate all bands (slices) to one list for central spaxel spectra.append(c1) spectraCorr.append(SimpleSpectrum(c1e2pCorr)) mygoods.append(bandBorders[band]) mygoodsCorr.append(bandBorders[band]) allPacs=prepPacsSpec(spectra,goods=mygoods) # to set your own good/bad flags allPacsCorr=prepPacsSpec(spectraCorr,goods=mygoodsCorr) # to set your own good/bad flags ###################### ### THE SPIRE PART ### ###################### cal = spireCal(pool="spire_cal_12_1") ### 1. get your observation, e.g. to test you can get the public NGC7027 data in the HSA #obs = getObservation(1342240019,useHsa=True) path = '/home/epuga/CROSS/SPIRE_analysis/Peter/DATA_Ros/cal_obs_hipe11/' obsid=1342259592 #NGC7027 filename = '%s/0x%8.8X_level2_HIPE11_ALPHA.fits'%(path,obsid) beamRef=60. diamS=12. #NGC7027 sds = FitsArchive().load(filename) ### 2. combine the two central detectors into one product allSpire = prepSpireSpec(sds,toWavelength=True) # basic with wavelength set to m-um ### 3. Apply SECT correction sds= semiExtendedCorrector(spectrum=sds, calibration=cal, gaussRefBeamDiam=beamRef, sourceModel=SemiExtendedSourceModelShape('gaussian', diamS, 0.0, 0.0, 0.0, 0.0, 1.0)) ### 4. combine the two central detectors into one product allSpireCorr = prepSpireSpec(sds,toWavelength=True) # basic with wavelength set to m-um ####################################### ### COMBINE THE PACS WITH THE SPIRE ### ####################################### # After having created allSpire and allPacs following the instructions here, you can # combine the two inputs into a final Spectrum1d. You can chose from a selection of # wavelength units but the flux units should be "Jy". The PACS spectra are added # first, taking on their segment numbers, and the SPIRE then added, segment numbers # following on from the last PACS one. In this way the spectra can still be distinguished # from each other. Appropriate meta data are added (you can view these by using the # Dataset viewer on the Spectrum1d, the input and output ones) # ### Examples of use allPacsNSpire1 = mergePacsNSpireSpec(allPacs,allSpire) # default allPacsNSpire2 = mergePacsNSpireSpec(allPacs,allSpire,waveunit="GHz") allPacsNSpireCorr1 = mergePacsNSpireSpec(allPacsCorr,allSpireCorr) # default allPacsNSpireCorr2 = mergePacsNSpireSpec(allPacsCorr,allSpireCorr,waveunit="GHz") pp=PlotXY() colorSpOff=java.awt.Color.black nameSpOff='Point Source' lineStOff=Style.SOLID colorSpOn=java.awt.Color.red nameSpOn='Semi-Extended Correction' lineStOn=Style.SOLID #this is the basic: PACS can have central3x3 True or False for i in range(allPacsNSpire2.getSegmentCount()+1): try: freq,flux, flag = allPacsNSpire2.getSegment(i).wave, allPacsNSpire2.getSegment(i).flux, allPacsNSpire2.getSegment(i).flag idx= flag.where(flag == 0) if i == 0: pp.addLayer(LayerXY(freq[idx], flux[idx], line=lineStOff, color=colorSpOff, name=nameSpOff, inLegend=True)) else: pp.addLayer(LayerXY(freq[idx], flux[idx], line=lineStOff, color=colorSpOff, name=nameSpOff, inLegend=False)) except: pass #this is the spectrum with the SE correction for i in range(allPacsNSpireCorr2.getSegmentCount()+1): try: freq,flux, flag = allPacsNSpireCorr2.getSegment(i).wave, allPacsNSpireCorr2.getSegment(i).flux, allPacsNSpireCorr2.getSegment(i).flag idx= flag.where(flag == 0) if i == 0: pp.addLayer(LayerXY(freq[idx], flux[idx], line=lineStOn, color=colorSpOn, name=nameSpOn, inLegend=True)) else: pp.addLayer(LayerXY(freq[idx], flux[idx], line=lineStOn, color=colorSpOn, name=nameSpOn, inLegend=False)) except: pass pp.legend.visible=1 pp.legend.position = Position.RIGHTMIDDLE pp.titleText=allPacsNSpire2.meta['object'].value #pp.setXtitle("Wavelength ($\micro$m)") pp.setXtitle("Frequency (GHz)") pp.setYtitle("Observed Flux (Jy)") pp.xaxis.tick.gridLines = 1 pp.yaxis.tick.gridLines = 1 pp.xaxis.range=[400., 5500.] pp.setPlotSize(8, 3.0) pp.xaxis.type=Axis.LOG pp.yaxis.type=Axis.LOG c = 299792458.0 # m/s xaux = ReciprocalAuxAxis(1.0e-3*c) pp.getXaxis().removeAuxAxis(0) pp.getXaxis().addAuxAxis(xaux) xaux.setTitleText("Wavelength ($\micro$m)") #Tools>>ExternalTools>>VOSpec #Source Query>>Observational Spectra Services >> The ISO Data Archive InterOperability System Query >> ISO LWS01 f=convertUnits(data=p_VOSpec.wave, sourceUnit='microm', targetUnit='GHz', quantity='Frequency') pp.addLayer(LayerXY(Double1d(f), Double1d(p_VOSpec.flux), line=1, color=java.awt.Color.gray, name='ISO LWS'))