########################################################## ## ## ## DEMO SCRIPT FOR HERSCHEL ARCHIVAL WORKSHOP ## ## WORKFLOW FOR HIFI POST-PROCESSING OF POINT, MAPPING ## ## AND SPECTRAL SCAN DATA OBSERVATIONS ## ## ## ## DT - 20-10-2016 ## ## ## ########################################################## ########################################################## ## POINTED MODE OBSERVATIONS ## ########################################################## ##1. LOAD OBSERVATIONS FROM ARCHIVE OR FROM POOL obsid = 1342244603 obsFromHSA = True #set to False if working from a pool if obsFromHSA: obs = getObservation(obsid,useHsa=1) else: myPoolName = poolName #pool located in ~/.hcss/lstore here obs = getObservation(obsid,poolName=myPoolName) ''' Double click on the obs variable to visualise the products Inspect the various product levels, the quality products, OFF positions, etc. You can list the quality flags with the following command print obs.quality You can visualise the OFF positions by creating the following variable (here for the WBS-H) and open it in the spectrum explorer by double-click on it. offWBSH = obs.refs["calibration"].product.refs["pipeline-out"].product.refs["ReferenceSpectra"].product.refs["WBS-H-USB"].product ''' ############################################################################### ##2. Flag your data: use of flagTool task ## ''' - open the flagTool task to assign/remove flags manually. The flags to be set by Users are: * LINE/BRIGHT_LINE to mask lines in future tasks like baseline or fringe removal * IGNORE_DATA or SPUR_CANDIDATE for channels to be discarded at any further stage of the processing - when you are done with one polarisation, proceed with the following one. Optionally you can use the flags created for one polarisation on the other polarisation - NOTE that the flagTool can work rowflags by setting the category option to either "both" or "rowflag" ''' #Path where masks will be stored - to be changed by User yourDiskPath = "/Users/dteyssier/Work/WORKSHOP/ArchivalWS2016/" # Inspect flags set by the pipeline # Call flagtool task as a method ft = FlagTool() # Launch task on the level2 data of the WBS-H-USB ft.process(obs, 'WBS-H-USB', level=2, category="both") # When you're done, the resulting masks are copied into the # observation context. Channel and row flag masks are in two # different tables # # Save them for any further use and bookkeepeing: chanMaskPointWBSH = obs.refs["level2"].product.refs["WBS-H-USB"].product["Linemasks"] simpleFitsWriter(chanMaskPointWBSH,yourDiskPath+"/chanMaskPointWBSH.fits") ## !!! BEWARE this will not work for versions ealier than HIPE 12.1 !!! #rowMaskPointWBSH = obs.refs["level2"].product.refs["WBS-H-USB"].product["Rowmasks"] #simpleFitsWriter(rowMaskPointWBSH,yourDiskPath+"/rowMaskPointWBSH.fits") # Now apply masks created for H to the V polarisation # ''' Note: If you have exited HIPE you can reload the channel masks from disk and re-apply them to the WBS-H polarization chanMaskPointWBSH = simpleFitsReader(yourDiskPath+"/chanMaskPointWBSH.fits") ft.addMaskTables({"WBS-H-USB": chanMaskPointWBSH}) #Loads masks ft.process(obs, 'WBS-H-USB', level=2, category="both") ''' # ##Now apply WBS-H mask as starting point for WBS-V ft.addMaskTables({"WBS-V-USB": chanMaskPointWBSH}) # Proceed with adjusting flags for the V polarisation ft.process(obs, 'WBS-V-USB', level=2, category="both") # When you're done, the resulting masks are copied into the # observation context. Save them (change file path): chanMaskPointWBSV = obs.refs["level2"].product.refs["WBS-V-USB"].product["Linemasks"] simpleFitsWriter(chanMaskPointWBSV,yourDiskPath+"/chanMaskPointWBSV.fits") ## !!! BEWARE this will not work before HIPE 12.1 !!! #rowMaskPointWBSV = obs.refs["level2"].product.refs["WBS-V-USB"].product["Rowmasks"] #simpleFitsWriter(rowMaskPointWBSV,yourDiskPath+"/rowMaskPointWBSV.fits") '''Additional notes: There also exist methods to overwrite, append or delete channel and/or rowflag mask tables to a given observation. See the documentation at: herschel.esac.esa.int/hcss-doc-14.0/load/hifi_um/html/hdrg_flags_flagtool.html ''' #Save your data if you wish to saveObservation(obs,poolName=str(obsid)+"Flagged",saveCalTree=True) #Now you are ready for defringing and baseline cleaning ############################################################################### ##3. Standing wave removal ## ##3.1 Electrical standing waves: only for observations in bands 6 and 7 ''' The following is an extract from a more detailed tutorial given during the spectroscopy workshop in April 2014. See: http://herschel.esac.esa.int/twiki/bin/view/Public/OnlineSpectroscopyWorkshop2014 ''' # Read the ESW models from an external file. # From HIPE 13 onwards this will be present in the calibration product itself. # Also available from HifiCalWeb page (http://herschel.esac.esa.int/twiki/bin/view/Public/HifiCalibrationWeb) # wget http://herschel.esac.esa.int/twiki/pub/Public/HifiCalibrationWeb/hebCorrectionModels_HICAL15.fits fits = FitsArchive() catalog = fits.load( os.sep.join([yourDiskPath,'hebCorrectionModels_HICAL15.fits']) ) # Process (remove ESWs from) the observation with default input parameters # Data are corrected at Level 1 hct = hebCorrection # In this call we let the task automatically figure out where # lines may be lying and apply masks accordingly. tds = hct( obs=obs, catalog=catalog, plot=2) # # Have a look at the plots if interested plotList = hct.plotList # # Process final products from the corrected level 1.0 data obs = hifiPipeline( obs=obs, fromLevel=1.0 ) # #Optionnally you can save the corrected data saveObservation(obs,poolName=str(obsid)+"HEBCorrected",saveCalTree=True) #You can also save the table data-set containing the best correction model simpleFitsWriter(product=tds, file=os.sep.join([yourDiskPath,'tds'+str(obsid)+'.fits']) ) ''' There are a couple of things you could do alternatively to the above: 1. Indicate manually the location of line masks. This is useful for broad or complicated line profiles not optimally picked up by the automatic line masking algorithm. In HIPE 12 this has to be given in the IF scale, in HIPE 13 it will be possible to enter it in sky frequencies tds2 = hct( obs=obs, catalog=catalog, plot=2, exclude=Double1d([7400,7900]) ) 2. From HIPE 13 onwards, you will be able to short-cut the model fitting part and directly apply the solution found and stored in variable tds. This is useful if you want to apply the correction to data newly reprocessed in the archive tdsRef = simpleFitsReader(file=os.sep.join([yourDiskPath,'tds'+str(obsid)+'.fits']) ) tds = hct( obs=obs, catalog=catalog, plot=2, redo = tdsRef ) ''' ################################################################# ##3.2 Optical standing waves: use of fitHifiFringe ''' We will use the task without any baseline subtraction in order to preserve the continuum level in these absorption line observations In the treated band (7b), after electrical standing wave removal, the main optical standing wave contribution occur at periods ~100 and ~650 MHz, so we will restrict the period domain in this range and fit 2 components. The task is capable of performing an automatic line/spur identification (automask=True). Since we have spent the time on carefully flagging the data, we will not use this feature. Note that when you treat the USB spectra of a given spectrometer, the correction will be automatically applied in a similar way to the LSB spectra The standing wave components fitted by the task are stored in the observation context under e.g.: obs.refs["calibration"].product.refs["pipeline-out"].product.refs["StandingWaves_WBS-H-USB"].product and equivalent products for other polarisations and sidebands This contains for example the median continuum level that was subtracted in case you selected option subBase=True ''' #Run the fringe removal polarisation per polarisation fitHifiFringe(obs_or_htp=obs, nfringes=2, subBase=False, startPeriod=80.0, endPeriod = 1000.0, \ automask = False, product='WBS-H-USB') fitHifiFringe(obs_or_htp=obs, nfringes=2, subBase=False, startPeriod=80.0, endPeriod = 1000.0, \ automask = False, product='WBS-V-USB') ''' In case you are not interested in the continuum, you can use the subBase=1 option to remove a spline fit to the data baseline - this would supersede the action of the fitBaseline task described thereafter ''' ############################################################################### ##4. Baseline correction/subtraction ## ''' The baselines will consist of polynomial fits where the user choses the polynomial order The data treated here have continuum information we want to preserve. As such no baseline correction is necessary. Optionally in those cases thought one could decide to normalise the absorption level to the continuum level, with option basemode = "div". To subtract the continuum option basemode = "sub" should be chosen (default) The task is capable of performing an automatic line/spur identification (domask=2). Since we have spent the time on carefully flagging the data, we will not use this feature. ''' #Run the baseline correction polarisation per polarisation #Using here a polynomial order 7 fitBaseline(data=obs,doglue=True, backends=['WBS-H-USB'], order=7, domask=0, basemode='div') fitBaseline(data=obs,doglue=True, backends=['WBS-V-USB'], order=7, domask=0, basemode='div') ''' Other interesting options: - each sub-band can be treated separately (doglue=False) or in a stitched fashion - Polynom order can be adjusted for each spectrum (interactive mode), or a generic order can be used for all spectra and the "non-stop" button can be selected ''' #Optionnally you can save the corrected data - these are in principle #the final products for this point mode observation saveObservation(obs,poolName=str(obsid)+"Final",saveCalTree=True) ###################################################################### #The following covers further processing you may be interested in #5. Averaging H and V spectra ''' There are two possible tasks for this: - polarPair - accumulate (more generic) ''' #Extract the corresponding datasets - here WBS data in the USB scale htpH = obs.refs["level2"].product.refs["WBS-H-USB"].product.copy() htpV = obs.refs["level2"].product.refs["WBS-V-USB"].product.copy() wbsTotalpolarPair = polarPair(htp=htpH,htp2=htpV) #Alternatively use accumulate: this works on data-sets wbsH = obs.refs["level2"].product.refs["WBS-H-USB"].product.refs["box_001"].product["0001"].copy() wbsV = obs.refs["level2"].product.refs["WBS-V-USB"].product.refs["box_001"].product["0001"].copy() wbsTotalaccumulate = accumulate(ds=[wbsH,wbsV]) #6. Unit/scale transformations #6.1 Conversion to Tmb scale htpTmb = doMainBeamTemp(htp=wbsTotalpolarPair, cal=obs.calibration) #6.2 Conversion to velocity scale - reference frequency arbitrary here htpVelo = convertWavescale(ds=wbsTotalaccumulate,to="km/s",reference=1764.0,overwrite=False) #6.3 Conversion to flux density - BEWARE that this involves a decision on # the source size and morphology # Please refer to the Spectroscopy WS held in April 2014 for more details: # http://herschel.esac.esa.int/twiki/pub/Public/OnlineSpectroscopyWorkshop2014/HIFI_IntensityExtraction_April2014.pdf # ## Option 1: assume source is point-like, i.e. size = 0 htpJy = convertK2Jy(htp=wbsTotalpolarPair.copy(),size=0.0, hifiBeam=1, overwrite=False, cal=obs.calibration) # ## Option 2: assume source is a perfect Gaussian of HPBW 10", and use real HIFI beam htpJy = convertK2Jy(htp=wbsTotalpolarPair.copy(),shape = "gauss", size=10.0, hifiBeam=1,overwrite=False, cal=obs.calibration) # ## Option 3: provide source morphology model, and use real HIFI beam ## Here we take an observation of a semi-extended source and ## use e.g. a PACS photometer map an brightness distribution assumption obs = getObservation(1342191560,useHsa=1) obsPacs = getObservation(1342250873,useHsa=1,instrument='PACS') Source = obsPacs.refs["level2_5"].product.refs["HPPHPFMAPB"].product specK = obs.refs["level2_5"].product.refs["spectrum"].product.refs["spectrum_WBS-H-USB"].product.copy() specJy = convertK2Jy(spectrum=specK,inputShape=Source, hifiBeam=1,overwrite=False, cal=obs.calibration) # ########################################################## ## SPECTRAL SCAN OBSERVATIONS ## ########################################################## ##1. LOAD OBSERVATIONS FROM ARCHIVE OR FROM POOL obsid = 1342218527 obsFromHSA = True #set to False if working from a pool if obsFromHSA: obs = getObservation(obsid,useHsa=1) else: myPoolName = poolName #pool located in ~/.hcss/lstore here obs = getObservation(obsid,poolName=myPoolName) ''' Double click on the obs variable to visualise the products Inspect the various product levels, the quality products, OFF positions, etc. You can list the quality flags with the following command print obs.quality You can visualise the OFF positions by creating the following variable (here for the WBS-H) and open it in the spectrum explorer by double-click on it. offWBSH = obs.refs["calibration"].product.refs["pipeline-out"].product.refs["ReferenceSpectra"].product.refs["WBS-H-USB"].product ''' ############################################################################### ##2. Flag your data: use of flagTool task ## ## NOTE THAT SPECTRAL SCAN OBSERVATIONS PROCESSED WITH HIPE 13 AND LATER ## ALREADY BENEFIT FROM A DEDICATED FLAGGING AND THEREFORE THE NEED FOR ## FURTHER FLAGGING SHOULD NOT BE NECESSARY UNLESS THE USER WANTS TO REFINE ## THE ASSUMED MASK TABLES. ## ''' - open the flagTool task to assign/remove flags manually. The flags to be set by Users are: * LINE/BRIGHT_LINE to mask lines in future tasks like baseline or fringe removal * IGNORE_DATA or SPUR_CANDIDATE for channels to be discarded at any further stage of the processing - when you are done with one polarisation, proceed with the following one. Optionally you can use the flags created for one polarisation on the other polarisation - NOTE that the flagTool can work rowflags by setting the category option to either "both" or "rowflag" ''' #Path where masks will be stored - to be changed by User yourDiskPath = "/Users/dteyssier/Work/WORKSHOP/ArchivalWS2016/" # Inspect flags set by the pipeline # Call flagtool task as a method ft = FlagTool() # Launch task on the level2 data of the WBS-H-USB ft.process(obs, 'WBS-H-USB', level=2, category="both") # When you're done, the resulting masks are copied into the # observation context. Channel and row flag masks are in two # different tables # ## Save them for any further use and bookkeepeing: chanMaskScanWBSH = obs.refs["level2"].product.refs["WBS-H-USB"].product["Linemasks"] simpleFitsWriter(chanMaskScanWBSH,yourDiskPath+"/chanMaskScanWBSH.fits") ## !!! BEWARE this will not work before HIPE 12.1 !!! #rowMaskScanWBSH = obs.refs["level2"].product.refs["WBS-H-USB"].product["Rowmasks"] #simpleFitsWriter(rowMaskScanWBSH,yourDiskPath+"/rowMaskScanWBSH.fits") # Now apply masks created for H to the V polarisation # ''' Note: If you have exited HIPE you can reload the channel masks from disk and re-apply them to the WBS-H polarization chanMaskScanWBSH = simpleFitsReader(yourDiskPath+"/chanMaskScanWBSH.fits") ft.addMaskTables({"WBS-H-USB": chanMaskScanWBSH}) #Append masks to existing flag table ft.process(obs, 'WBS-H-USB', level=2, category="both") ''' # ##Now apply WBS-H mask as starting point for WBS-V ft.addMaskTables({"WBS-V-USB": chanMaskScanWBSH}) # Proceed with adjusting flags for the V polarisation ft.process(obs, 'WBS-V-USB', level=2, category="both") # When you're done, the resulting masks are copied into the # observation context. Save them (change file path): chanMaskScanWBSV = obs.refs["level2"].product.refs["WBS-V-USB"].product["Linemasks"] simpleFitsWriter(chanMaskScanWBSV,yourDiskPath+"/chanMaskScanWBSV.fits") ## !!! BEWARE this will not work before HIPE 12.1 !!! #rowMaskScanWBSV = obs.refs["level2"].product.refs["WBS-V-USB"].product["Rowmasks"] #simpleFitsWriter(rowMaskScanWBSV,yourDiskPath+"/rowMaskScanWBSV.fits") '''Additional notes: There also exist methods to overwrite, append or delete channel and/or rowflag mask tables to a given observation. See the documentation at: herschel.esac.esa.int/hcss-doc-14.0/load/hifi_um/html/hdrg_flags_flagtool.html ''' #Save your data if you wish to saveObservation(obs,poolName=str(obsid)+"Flagged",saveCalTree=True) #Now you are ready for defringing and baseline cleaning ############################################################################### ##3. Standing wave removal ## ##3.1 Electrical standing waves: only for observations in bands 6 and 7 ''' The following is an extract from a more detailed tutorial given during the spectroscopy workshop in April 2014. See: http://herschel.esac.esa.int/twiki/bin/view/Public/OnlineSpectroscopyWorkshop2014 ''' # Read the ESW models from an external file. # In HIPE 13 this will be in the calibration product. # Currently available from HifiCalWeb page (http://herschel.esac.esa.int/twiki/bin/view/Public/HifiCalibrationWeb) # wget http://herschel.esac.esa.int/twiki/pub/Public/HifiCalibrationWeb/hebCorrectionModels_HICAL15.fits fits = FitsArchive() catalog = fits.load( os.sep.join([yourDiskPath,'hebCorrectionModels_HICAL15.fits']) ) # Process (remove ESWs from) the observation with default input parameters # Data are corrected at Level 1 hct = hebCorrection # In this call we let the task automatically figure out where # lines may be lying and apply masks accordingly. tds = hct( obs=obs, catalog=catalog, plot=2) # # Have a look at the plots if interested plotList = hct.plotList # # Process final products from the corrected level 1.0 data obs = hifiPipeline( obs=obs, fromLevel=1.0 ) # #Optionnally you can save the corrected data saveObservation(obs,poolName=str(obsid)+"HEBCorrected",saveCalTree=True) #You can also save the table data-set containing the best correction model simpleFitsWriter(product=tds, file=os.sep.join([yourDiskPath,'tds'+str(obsid)+'.fits']) ) ''' There are a couple of things you could do alternatively to the above: 1. Indicate manually the location of line masks. This is useful for broad or complicated line profiles not optimally picked up by the automatic line masking algorithm. In HIPE 12 this has to be given in the IF scale, in HIPE 13 it will be possible to enter it in sky frequencies tds2 = hct( obs=obs, catalog=catalog, plot=2, exclude=Double1d([7400,7900]) ) 2. From HIPE 13 onwards, you will be able to short-cut the model fitting part and directly apply the solution found and stored in variable tds. This is useful if you want to apply the correction to data newly reprocessed in the archive tdsRef = simpleFitsReader(file=os.sep.join([yourDiskPath,'tds'+str(obsid)+'.fits']) ) tds = hct( obs=obs, catalog=catalog, plot=2, redo = tdsRef ) ''' ################################################################# ##3.2 Optical standing waves: use of fitHifiFringe ''' We will here use the task without any baseline subtraction, however if you are not interested in the continuum see the alternative call thereafter In the treated band (1a) the main optical standing wave contribution occur at periods ~100 MHz, so we will restrict the period domain in this range and fit 2 components. The task is capable of performing an automatic line/spur identification (automask=True). Since we have spent the time on carefully flagging the data, we will not use this feature. Note that when you treat the USB spectra of a given spectrometer, the correction will be automatically applied in a similar way to the LSB spectra The standing wave components fitted by the task are stored in the observation context under e.g.: obs.refs["calibration"].product.refs["pipeline-out"].product.refs["StandingWaves_WBS-H-USB"].product and equivalent products for other polarisations and sidebands This contains for example the median continuum level that was subtracted in case you selected option subBase=True ''' #Run the fringe removal polarisation per polarisation fitHifiFringe(obs_or_htp=obs, nfringes=2, subBase=False, startPeriod=80.0, endPeriod = 200.0, \ automask = False, product='WBS-H-USB') fitHifiFringe(obs_or_htp=obs, nfringes=2, subBase=False, startPeriod=80.0, endPeriod = 200.0, \ automask = False, product='WBS-V-USB') ''' In case you are not interested in the continuum, you can use the subBase=1 option to remove a spline fit to the data baseline - this would supersede the action of the fitBaseline task described thereafter ''' ############################################################################### ##4. Baseline correction/subtraction ## ''' The baselines will consist of polynomial fits where the user choses the polynomial order We assume here that we are not interested in preserving the continuum information. As such we will use the option basemode = "sub" (default) The task is capable of performing an automatic line/spur identification (domask=2). Since we have spent the time on carefully flagging the data, we will not use this feature. ''' #Run the baseline correction polarisation per polarisation #Using here a polynomial order 7 fitBaseline(data=obs,doglue=True, backends=['WBS-H-USB'], order=7, domask=0, basemode='sub') fitBaseline(data=obs,doglue=True, backends=['WBS-V-USB'], order=7, domask=0, basemode='sub') ''' Other interesting options: - each sub-band can be treated separately (doglue=False) or in a stitched fashion - Polynom order can be adjusted for each spectrum (interactive mode), or a generic order can be used for all spectra and the "non-stop" button can be selected ''' #Optionnally you can save the corrected data - these are in principle #the final products for this point mode observation saveObservation(obs,poolName=str(obsid)+"BaselineCorrected",saveCalTree=True) ###################################################################### #5. Generation of deconvolved spectra ''' Once all artifacts have been removed, the single single band solutions need to be re-generated. For this the doDeconvolution task shall be used. Deconvolution works polarisation per polarisation. Deconvolution will be very sensitive to the flagging strategy, in that any spur or unruly frequency range kept in the input data will kick in in the output data in the form of ghosts and invalid spectral information If the flagging has been done properly, doDeconvolution should be called with the following options set (default): - gain = GAIN_FIT_OFF_USE_PRESET: we will use the a priori knowledge of the sideband ratio stored in the calibration files - spur_rejection = DO_NOT_REJECT_SCAN_WITH_SPURS: we will only ignore data flagged as SPUR_CANDIDATE and IGNORE_DATA ''' #Task here called with default parameters deconH = doDeconvolution(obs=obs, spectrometer='WBS-H') deconV = doDeconvolution(obs=obs, spectrometer='WBS-V') #Optionnally you can save the corrected data - these are in principle #the final products for this point mode observation saveObservation(obs,poolName=str(obsid)+"Final",saveCalTree=True) #You can also create a deconvolved spectrum from the OFF position data deconH_OFF = doDeconvolution(obs=obs, spectrometer='WBS-H',use_reference=True) deconV_OFF = doDeconvolution(obs=obs, spectrometer='WBS-V',use_reference=True) ###################################################################### #The following covers further processing you may be interested in #6. Averaging H and V spectra ''' For spectral scans, the end product is the level 2.5 deconvolved spectra. Like for the pointed mode observations, two possible tasks here - polarPair - accumulate ''' #Extract the corresponding datasets - here WBS data in the USB scale ssbH = deconH["ssb"] ssbV = deconV["ssb"] ssbTotalpolarPair = polarPair(ds1=ssbH,ds2=ssbV) #Alternatively use accumulate: this works on data-sets ssbTotalaccumulate = accumulate(ds=[ssbH,ssbV],pointingTolerance=Double.NaN) #7. Unit/scale transformations #7.1 Conversion to Tmb scale ''' Because the doMainBeamTemp task works on HifiTimelineProduct, but the deconvolution output are so-called Spectrum1a, the task cannot be directly used on the level 2.5 Instead, the level2 products need to be transformed, and then deconvolved ''' #Copy observation context to not alter the content of the reference obsRef = obs.copy() #This will take some time if you work from archive htpH = obsRef.refs["level2"].product.refs["WBS-H-USB"].product htpV = obsRef.refs["level2"].product.refs["WBS-V-USB"].product htpTmbH = doMainBeamTemp(htp=htpH, cal=obs.calibration) htpTmbV = doMainBeamTemp(htp=htpV, cal=obs.calibration) #Deconvolve deconTmbH = doDeconvolution(obs=obsRef, spectrometer='WBS-H') deconTmbV = doDeconvolution(obs=obsRef, spectrometer='WBS-V') #Add polarizations ssbTmbTotalpolarPair = polarPair(ds1=deconTmbH["ssb"],ds2=deconTmbV["ssb"]) #7.2 Conversion to velocity scale - reference frequency arbitrary here ssbVelo = convertWavescale(ds=ssbTotalpolarPair,to="km/s",reference=950.0,overwrite=False) #7.3 Conversion to flux density - BEWARE that this involves a decision on # the source size and morphology # Please refer to the Spectroscopy WS held in April 2014 for more details: # http://herschel.esac.esa.int/twiki/pub/Public/OnlineSpectroscopyWorkshop2014/HIFI_IntensityExtraction_April2014.pdf ''' Here the same limitation as for the doMainBeamTemp task applies to the conversion to flux density has to be done at level2 prior to deconvolution ''' #Copy observation context to not alter the content of the reference obsRef = obs.copy() #This will take some time if you work from archive htpH = obsRef.refs["level2"].product.refs["WBS-H-USB"].product htpV = obsRef.refs["level2"].product.refs["WBS-V-USB"].product htpJyH = convertK2Jy(htp=htpH,size=0.0, overwrite=False, cal=obs.calibration) htpJyV = convertK2Jy(htp=htpV,size=0.0, overwrite=False, cal=obs.calibration) #Deconvolve deconJyH = doDeconvolution(obs=obsRef, spectrometer='WBS-H') deconJyV = doDeconvolution(obs=obsRef, spectrometer='WBS-V') #Add polarizations ssbJyTotalpolarPair = polarPair(ds1=deconJyH["ssb"],ds2=deconJyV["ssb"]) ############################################################################### ##8. Deconvolution of multiple obsids ## Example given on the spectral scans in bands 1a and 1b for OMC-2 FIR4 obsid1 = 1342218633 obsid2 = 1342191591 #Read observation from archive obs1a = getObservation(obsid1,useHsa=1) obs1b = getObservation(obsid2,useHsa=1) #Extract level 2.5 deconvolved products - here only for WBS-H decon1a = obs1a.refs["level2_5"].product.refs["myDecon"].product.refs["myDecon_WBS-H"].product["ssb"] decon1b = obs1b.refs["level2_5"].product.refs["myDecon"].product.refs["myDecon_WBS-H"].product["ssb"] #Create concatenated spectrum from individual SSB spectra deconsep = accumulate(ds=[decon1a,decon1b],pointingTolerance=Double.NaN) #Deconvolve using merged band 1a and 1b data - here only for WBS-H decon = doDeconvolution(obs1a, obs2_array=[obs1b]) deconboth = decon["ssb"] #Visualise results p=splot(useFrame=1) p.add(deconsep) p.add(deconboth) ############################################################################### ##10. Deconvolution of bright lines avoiding negative ghosts ## Example given on the Orion bar spectral scan in band 4a ## flagged in the previous exercise, on which fringes and ## baselines have been corrected off-line ## ## The example below is given for the H polarisation ## ## These steps are also contained in an HIFI useful script: ## Scripts > HIFI Useful Scripts > Bright Line Deconvolution #First create an SSB solution where the bright line is maintained decon_bright_notflagged = doDeconvolution(obs=obs, ignore_bright_line=0) ##Open the above result to inspect the negative ghosts #Now create a second SSB solution where the bright line is ignored decon_bright_flagged = doDeconvolution(obs=obs, ignore_bright_line=1) ##We now need to concantenate those two spectra only at the location of the bright lines #Extract SSB spectra from decon output decBright = decon_bright_notflagged["ssb"] decNoBright = decon_bright_flagged["ssb"] ##Extract flux from decBright in channels where decNoBright is NaN freqInd = decNoBright["flux"].data.where((IS_NAN(decNoBright["flux"].data))) decNoBright["flux"].data[freqInd] = decBright["flux"].data[freqInd] #The end result with bright lines and ghosts removed is in variable decNoBright ##Compare outputs in the area of negative ghosts p=splot(useFrame=1) p.add(decBright) p.add(decNoBright) ########################################################## ## SPECTRAL MAPPING OBSERVATIONS ## ########################################################## ##1. LOAD OBSERVATIONS FROM ARCHIVE OR FROM POOL obsid = 1342249579 obsFromHSA = True #set to False if working from a pool if obsFromHSA: obs = getObservation(obsid,useHsa=1) else: myPoolName = poolName #pool located in ~/.hcss/lstore here obs = getObservation(obsid,poolName=myPoolName) ''' Double click on the obs variable to visualise the products Inspect the various product levels, the quality products, OFF positions, etc. You can list the quality flags with the following command print obs.quality You can visualise the OFF positions by creating the following variable (here for the WBS-H) and open it in the spectrum explorer by double-click on it. offWBSH = obs.refs["calibration"].product.refs["pipeline-out"].product.refs["ReferenceSpectra"].product.refs["WBS-H"].product ''' ############################################################################### ##2. Flag your data: use of flagTool task ## ''' - open the flagTool task to assign/remove flags manually. The flags to be set by Users are: * LINE/BRIGHT_LINE to mask lines in future tasks like baseline or fringe removal * IGNORE_DATA or SPUR_CANDIDATE for channels to be discarded at any further stage of the processing - when you are done with one polarisation, proceed with the following one. Optionally you can use the flags created for one polarisation on the other polarisation - NOTE that the flagTool can work rowflags by setting the category option to either "both" or "rowflag" ''' #Path where masks will be stored - to be changed by User yourDiskPath = "/Users/dteyssier/Work/WORKSHOP/ArchivalWS2016/" # Inspect flags set by the pipeline # Call flagtool task as a method ft = FlagTool() # Launch task on the level2 data of the WBS-H-USB ft.process(obs, 'WBS-H-USB', level=2, category="both") ''' In this particular obsid, no noticeable spurious feature is present so the masking concerns the strong C+ line ''' # When you're done, the resulting masks are copied into the # observation context. Channel and row flag masks are in two # different tables # # Save them for any further use and bookkeepeing: chanMaskMapWBSH = obs.refs["level2"].product.refs["WBS-H-USB"].product["Linemasks"] simpleFitsWriter(chanMaskMapWBSH,yourDiskPath+"/chanMaskMapWBSH.fits") ## !!! BEWARE this will not work before HIPE 12.1 !!! #rowMaskMapWBSH = obs.refs["level2"].product.refs["WBS-H-USB"].product["Rowmasks"] #simpleFitsWriter(rowMaskMapWBSH,yourDiskPath+"/rowMaskMapWBSH.fits") # Now apply masks created for H to the V polarisation # ''' Note: If you have exited HIPE you can reload the channel masks from disk and re-apply them to the WBS-H polarization chanMaskMapWBSH = simpleFitsReader(yourDiskPath+"/chanMaskMapWBSH.fits") ft.addMaskTables({"WBS-H-USB": chanMaskMapWBSH}) #Loads masks ft.process(obs, 'WBS-H-USB', level=2, category="both") ''' # ##Now apply WBS-H mask as starting point for WBS-V ft.addMaskTables({"WBS-V-USB": chanMaskMapWBSH}) # Proceed with adjusting flags for the V polarisation ft.process(obs, 'WBS-V-USB', level=2, category="both") # When you're done, the resulting masks are copied into the # observation context. Save them (change file path): chanMaskMapWBSV = obs.refs["level2"].product.refs["WBS-V-USB"].product["Linemasks"] simpleFitsWriter(chanMaskMapWBSV,yourDiskPath+"/chanMaskMapWBSV.fits") ## !!! BEWARE this will not work before HIPE 12.1 !!! #rowMaskMapWBSV = obs.refs["level2"].product.refs["WBS-V-USB"].product["Rowmasks"] #simpleFitsWriter(rowMaskMapWBSV,yourDiskPath+"/rowMaskMapWBSV.fits") '''Additional notes: There also exist methods to overwrite, append or delete channel and/or rowflag mask tables to a given observation. See the documentation at: herschel.esac.esa.int/hcss-doc-14.0/load/hifi_um/html/hdrg_flags_flagtool.html ''' #Save your data if you wish to saveObservation(obs,poolName=str(obsid)+"Flagged",saveCalTree=True) #Now you are ready for defringing and baseline cleaning ############################################################################### ##3. Standing wave removal ## ##3.1 Electrical standing waves: only for observations in bands 6 and 7 ''' The following is an extract from a more detailed tutorial given during the spectroscopy workshop in April 2014. See: http://herschel.esac.esa.int/twiki/bin/view/Public/OnlineSpectroscopyWorkshop2014 BEWARE THAT ESW CORRECTION FOR MAP CAN BE EXTREMELY TIME AND MEMORY DEMANDING DUE TO THE LARGE AMOUNT OF SPECTRA TO BE TREATED ''' # Read the ESW models from an external file. # In HIPE 13 this will be in the calibration product. # Currently available from HifiCalWeb page (http://herschel.esac.esa.int/twiki/bin/view/Public/HifiCalibrationWeb) # wget http://herschel.esac.esa.int/twiki/pub/Public/HifiCalibrationWeb/hebCorrectionModels_HICAL15.fits fits = FitsArchive() catalog = fits.load( os.sep.join([yourDiskPath,'hebCorrectionModels_HICAL15.fits']) ) # Process (remove ESWs from) the observation with default input parameters # Data are corrected at Level 1 hct = hebCorrection # In this call we let the task automatically figure out where # lines may be lying and apply masks accordingly. tds = hct( obs=obs, catalog=catalog, plot=2) # # Have a look at the plots if interested plotList = hct.plotList # # Process final products from the corrected level 1.0 data obs = hifiPipeline( obs=obs, fromLevel=1.0 ) # #Optionnally you can save the corrected data saveObservation(obs,poolName=str(obsid)+"HEBCorrected",saveCalTree=True) #You can also save the table data-set containing the best correction model simpleFitsWriter(product=tds, file=os.sep.join([yourDiskPath,'tds'+str(obsid)+'.fits']) ) ''' There are a couple of things you could do alternatively to the above: 1. Indicate manually the location of line masks. This is useful for broad or complicated line profiles not optimally picked up by the automatic line masking algorithm. In HIPE 12 this has to be given in the IF scale, in HIPE 13 it will be possible to enter it in sky frequencies tds2 = hct( obs=obs, catalog=catalog, plot=2, exclude=Double1d([7400,7900]) ) 2. From HIPE 13 onwards, you will be able to short-cut the model fitting part and directly apply the solution found and stored in variable tds. This is useful if you want to apply the correction to data newly reprocessed in the archive tdsRef = simpleFitsReader(file=os.sep.join([yourDiskPath,'tds'+str(obsid)+'.fits']) ) tds = hct( obs=obs, catalog=catalog, plot=2, redo = tdsRef ) ''' ################################################################# ##3.2 Optical standing waves: use of fitHifiFringe ''' We will use the task without any baseline subtraction in order to preserve the continuum level in these absorption line observations In the treated band (7b), after electrical standing wave removal, the main optical standing wave contribution occur at periods ~100 and ~650 MHz, so we will restrict the period domain in this range and fit 2 components. The task is capable of performing an automatic line/spur identification (automask=True). In case the user has already careful flagged lines in an earlier step, then option automask=False can be used. For maps however an interesting option could consist in using a mask based on the average of all pixels in the map, achieving 1) a better SNR, 2) a window that accounts for possible line velocity/width changes over the map Note that when you treat the USB spectra of a given spectrometer, the correction will be automatically applied in a similar way to the LSB spectra The standing wave components fitted by the task are stored in the observation context under e.g.: obs.refs["calibration"].product.refs["pipeline-out"].product.refs["StandingWaves_WBS-H-USB"].product and equivalent products for other polarisations and sidebands This contains for example the median continuum level that was subtracted in case you selected option subBase=True ''' #Run the fringe removal polarisation per polarisation fitHifiFringe(obs_or_htp=obs, nfringes=2, subBase=False, startPeriod=80.0, endPeriod = 1000.0, \ automask = True, product='WBS-H-USB') fitHifiFringe(obs_or_htp=obs, nfringes=2, subBase=False, startPeriod=80.0, endPeriod = 1000.0, \ automask = True, product='WBS-V-USB') ''' In case you are not interested in the continuum, you can use the subBase=1 option to remove a spline fit to the data baseline - this would supersede the action of the fitBaseline task described thereafter ''' ############################################################################### ##4. Baseline correction/subtraction ## ''' The baselines will consist of polynomial fits where the user choses the polynomial order The data treated here have continuum information we want to preserve. As such no baseline correction is necessary. Optionally in those cases thought one could decide to normalise the absorption level to the continuum level, with option basemode = "div". To subtract the continuum option basemode = "sub" should be chosen (default) The task is capable of performing an automatic line/spur identification (domask=2). Since we have spent the time on carefully flagging the data, we will not use this feature. ''' #Run the baseline correction polarisation per polarisation #Using here a polynomial order 7 fitBaseline(data=obs,doglue=True, backends=['WBS-H-USB'], order=7, domask=0, basemode='sub') fitBaseline(data=obs,doglue=True, backends=['WBS-V-USB'], order=7, domask=0, basemode='sub') ''' Other interesting options: - each sub-band can be treated separately (doglue=False) or in a stitched fashion - Polynom order can be adjusted for each spectrum (interactive mode), or a generic order can be used for all spectra and the "non-stop" button can be selected ''' #Optionnally you can save the corrected data - these are in principle #the final products for this point mode observation saveObservation(obs,poolName=str(obsid)+"Final",saveCalTree=True) ###################################################################### #5. Generation of regridded cubes ''' Once all artifacts have been removed, the regridded cubes need to be re-generated. For this the doGridding task shall be used. doGridding will work per spectrometer, polarisation and spectrometer sub-band. In maps with a non-zero position angle, two categories of cubes will be created in the SPG: - cubes projected onto the equatorial frame, i.e. that will appear with their rotated orientation onto the plane of the sky - cubes rotated by the position angle, i.e. that will appear as a rectangular (or square) coverage doGridding uses a grid whose dimension are estimated based on the mapping sequence scheduled during the observation. These parameters can however be adjusted by the User See: http://herschel.esac.esa.int/hcss-doc-12.0/load/hifi_um/html/hum_dogridding_common.html ''' # The most simple call is as follows - e.g. for WBS USB cubes: htpH = obs.refs["level2"].product.refs["WBS-H-USB"].product.copy() htpV = obs.refs["level2"].product.refs["WBS-V-USB"].product.copy() cubeOutputH,cubes,cube,xPoints,yPoints,convolutionTable,grid = doGridding(htp=htpH) cubeOutputV,cubes,cube,xPoints,yPoints,convolutionTable,grid = doGridding(htp=htpV) # and the cube contexts are available as: cubesH = cubeOutputH cubesV = cubeOutputV #Optionnally you can save the corrected data - these are in principle #the final products for this point mode observation saveObservation(obs,poolName=str(obsid)+"Final",saveCalTree=True) ###################################################################### #The following covers further processing you may be interested in #6. Averaging H and V cubes ''' The general idea is to combine various HTPs into a single one that is then used as input for the gridding. This is not restricted to combining H and V but can be extended to combining sub-maps into larger maps ''' #Extract the corresponding datasets - here WBS data in the USB scale htpH = obs.refs["level2"].product.refs["WBS-H-USB"].product.copy() htpV = obs.refs["level2"].product.refs["WBS-V-USB"].product.copy() htpMerged = mergeHtps(htps=[htpH,htpV]) ##!!!! BEWARE THAT SOME BUG IS STILL AFFECTING THIS TASK IN HIPE 13 AND WILL ##!!!! REPAIRED ONLY IN HIPE 14 - PLEASE CONTACT THE HELPDESK IF YOU NEED SUPPORT BEFORE THAT cubeOutputMerged,cubes,cube,xPoints,yPoints,convolutionTable,grid = doGridding(htp=htpMerged) cubesMerged = cubeOutputMerged #7. Unit/scale transformations #7.1 Conversion to Tmb scale ''' Because the doMainBeamTemp task works on HifiTimelineProduct, the task cannot be directly used on the level 2.5 cubes Instead, the level2 products need to be transformed, and then regridded ''' #Copy observation context to not alter the content of the reference obsRef = obs.copy() #This will take some time if you work from archive htpH = obsRef.refs["level2"].product.refs["WBS-H-USB"].product htpV = obsRef.refs["level2"].product.refs["WBS-V-USB"].product htpTmbH = doMainBeamTemp(htp=htpH, cal=obs.calibration) htpTmbV = doMainBeamTemp(htp=htpV, cal=obs.calibration) #Regrid cubeOutputH,cubes,cube,xPoints,yPoints,convolutionTable,grid = doGridding(htp=htpTmbH) cubeOutputV,cubes,cube,xPoints,yPoints,convolutionTable,grid = doGridding(htp=htpTmbV) # cubesTmbH = cubeOutputH cubesTmbV = cubeOutputV #Add polarizations: see above #7.2 Conversion to velocity scale - reference frequecy set to C+ rest frequency cubeTmbVsub4 = cubesTmbV.refs['cube_WBS_V_USB_4'].product cubeVelTmbVsub4 = convertWavescale(ds=cubeTmbVsub4, to='km s-1', reference=1900.536, \ referenceUnit='GHz', overwrite=False) #7.3 Conversion to flux density - BEWARE that this involves a decision on # the source size and morphology # Please refer to the Spectroscopy WS held in April 2014 for more details: # http://herschel.esac.esa.int/twiki/pub/Public/OnlineSpectroscopyWorkshop2014/HIFI_IntensityExtraction_April2014.pdf ''' Here the same limitation as for the doMainBeamTemp task applies to the conversion to flux density has to be done at level2 prior to deconvolution ''' #Copy observation context to not alter the content of the reference obsRef = obs.copy() #This will take some time if you work from archive htpH = obsRef.refs["level2"].product.refs["WBS-H-USB"].product htpV = obsRef.refs["level2"].product.refs["WBS-V-USB"].product htpJyH = convertK2Jy(htp=htpH,size=0.0, overwrite=False, cal=obs.calibration) htpJyV = convertK2Jy(htp=htpV,size=0.0, overwrite=False, cal=obs.calibration) #Regrid cubeOutputH,cubes,cube,xPoints,yPoints,convolutionTable,grid = doGridding(htp=htpJyH) cubeOutputV,cubes,cube,xPoints,yPoints,convolutionTable,grid = doGridding(htp=htpJyV) # cubesJyH = cubeOutputH cubesJyV = cubeOutputV #Add polarizations: see above #8. Integrated intensity maps: the cube toolbox # wbsH4cube = cubesMerged.refs["cube_WBS_H_USB_4"].product.copy() # Crop around the line of interest ''' In the following, the hard-coded frequencies around 1900GHz are for obsid=1342249579 . ''' croppedCube = cropCube(cube=wbsH4cube, startWave=1900.422, endWave=1900.520) # Compute a velocity map by fitting a GAUSSIAN to each pixel # Take reference channel for velocity=0 as peak intenstity bin in the Bar [velocityMap, dispersionMap, maxFluxMap, lineIntensityMap, velCube, fittedLineCube] = computeVelocityMap(cube=croppedCube, \ reference=1900.47, \ velAlgorithm=herschel.ia.toolbox.cube.ComputeVelocityMapTask.Algorithm.GAUSSIAN_FIT) # Compute a line-integrated map from the cropped cube integratedCroppedCube = integrateSpectralMap(cube=croppedCube) ################################## ### THE END ### ##################################