# HiClass demo script
# prepared for Herschel KP DP workshop Dec 2009
# 
# v0.0 06 Dec 2009   CM: created
# v1.0 03 April 2010 CM: modified for HIFI PSP/SDP DP workshop, ESTEC, April 2010
#
# ---------------------------------------------------------------
# HiClass is used to save level-0.5 or higher HIFI HTP and/or 
# SpectrumDatasets into a FITS format readable by Class.
#
# hiClass will *not* export PACS or SPIRE spectra to Class.
#
#  Documentation about using HiClass can be found in Chapter 11
#  of the HIFI User Manual
# ---------------------------------------------------------------
# 
#
# Step 1. Get ObservationContext into session:
# 
# Here we get an ObservationContext, with observation identifier
# 1342180473 from a local pool called "HIFI_DP_WORKSHOP_080410_P_1". 
# Edit as appropriate to use this method.
#
 obs = getObservation("1342180473", poolName="pool_1342180473")
#
# If you do not have the observation store in a local pool, you can get it 
# from the HSA.
# Note, for this to work you must enter your HSA username and password. 
# It is convenient to save them in your ./hcss/user.props file.
#
# hcss.ia.pal.pool.hsa.haio.login_usr = your username 
# hcss.ia.pal.pool.hsa.haio.login_pwd = your password 
#
# Uncomment the line below (delete the #), and edit the observation number 
# as appropriate to use this method.
# obs = getObservation("1342180473", useHsa=True)
#
#  You can then save it to local pool with
#
# saveObservation(obs)
# 
#  be aware this can take some time and no notification message is given
# when save is complete.
#
#  Alternatively, use the Herschel Science Archive perspective in HIPE...
#
#
# Step 2.   Get the product you want to export to Class:   
#
# In this case a Level-2 HTP for the upper sideband of the WBS-H.
# Edit as appropriate for other levels or spectrometers.
# You can see the reference (obs.refs....) for each product by clicking
# on it in the ObservationContext tree viewed with the ContextViewer  
#
 HTP_WBS_H_USB = obs.refs["level2"].product.refs["WBS-H-USB"].product
#
#
# Step 3.  Simplest case - export entire HTP into Class:
#
# This method will give you a separate spectrum for each subband.
# Here, myhtp.fits will be saved in the directory HIPE is launched from, specify
# the absolute path to save elsewhere.
# See below for other methods. 
#
 HiClassTask()(product = HTP_WBS_H_USB, fileName = 'myhtp.fits')
#
# If you are exporting a map you should specify the nominal RA and dec of the 
# observation - the value you put into HSpot, in decimal degrees.
# This will not be necessary in the future when the HTP metadata contains this
# information, which is only in the ObservationContext metadata at the moment.
#
# HiClassTask()(product = HTP_WBS_H_USB,raNominal = 193.886, decNominal = -76.96, fileName = 'myhtp.fits')
#
#
# Step 4.  Read file into Class:
# Do this whatever method you use for exporting.
# Note that you need a recent version of Class, this was tested on the Oct 09 release.
# Versions of Class from the first half of 2009 are known to have a bug preventing reading of any FITS file with a long header 
#
# CLASS is not able to work directly within FITS files.  So you have to
# convert the FITS file into a CLASS file:
#   file out MyHIFISpectra.hifi mul
#   fits read myhtp.fits
# 
# Now you have a CLASS file named MyHIFISpectra.hifi (you can use whatever
# you want as an extension) you can access like you always do in CLASS:
#
#   file in MyHIFISpectra.hifi
#   find
#   get first
#   set unit f i
#   plot
#
#  The spectrometer "W" or "H" and the subband number are included in the plot title.
#
#
#   Other methods:
#
# Export single dataset(ds). Here we export the first dataset in the HTP.
# Again, this method will result in a spectrum for each segment in the dataset.
#
 ds =  HTP_WBS_H_USB.get(1)
 HiClassTask()(dataset = ds, fileName = 'myspectra.fits')
#
# 
#  Add multiple spectrum datasets to the FITS file, here two datasets, ds1 and ds2:
 ds1 = obs.refs["level1"].product.refs["WBS-H"].product.refs["3"].product["dataset"]
 ds2 = obs.refs["level1"].product.refs["WBS-H"].product.refs["4"].product["dataset"]
 h = HiClass()
 h.add(ds1)
 h.add(ds2)
 h.saveToFits("MyTwoSpectra.fits")
#
# Make a FITS file per spectrum dataset, stitching datasets so get one 
# spectrum per dataset:
 h = HiClass()
 ds1_stitched= stitch(ds = ds1) 
 h.add(ds1)
 h.saveToFits("fits1.fits")
 h = HiClass()
 ds2_stitched= stitch(ds = ds2) 
 h.add(ds2)
 h.saveToFits("fits2.fits")
#
#
#
# End of demo script