##########################################################
## 					     		##
##  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           ###
##################################