#hifidemo_scripting1_apr2010.py: demo of how to write scripts to
#    manipulate HIFI data in HIPE. Baseline subtraction using the 
#    median signal is taken as
#    an example, since HIPE does not provide a convenient
#    way to do this, presently.
#
#
# v0.0 04/Apr/2010 ACAB: created for April 2010 HIFI DP Workshop
# v0.1 07/Apr/2010 ACAB: use getObservation
#
#

#Basics of Jython scripting is described in the "Scripting and 
#Data Mining Manual", accessible from the Hipe Help menu.

#load band 4B SScan LChop observation 
#from ~/.hcss/lstore/
obs=getObservation("1342180554",poolName="pool_1342180554")

#extract WBS-H-USB product from ObservationContext
htp=obs.refs["level2"].product.refs["WBS-H-USB"].product

#This htp contains a collection of level 2 "spectrumdatasets".
#A single level 2 spectrumdataset, say number 60, can be selected as follows.
#(use the context viewer to find out which spectrumdatasets are available)
sds=htp.refs["60"].product["dataset"]

#A spectrumdataset can consist of multiple scans (read-outs), each of
#which can contain multiple 'segments' (WBS sub-bands or HRS modules).

#    -number of scans
scans=sds.rowCount
print scans
#    -array of segment indices, e.g. [1,2,3,4] for WBS bands 1-5.
segments=sds.getSegmentIndices()
print segments

#REST OF SCRIPT: As a first order baseline correction,
#                subtract median of flux over entire IF from
#                individual sub-bands.

#The example spectrumdataset is a WBS spectrum with 4 sub-bands.
#If we want to work on the overall spectrum, 'stitch' them first.
#(If you don't how to do this, use GUI under Applicable Tasks and copy
# equivalent command from console) 
sds_s = stitch(ds=sds,variant="crossoverPoints",edgeTolerance=0.01,stepsize=0.0)

#Get most basic spectrum to work on: a single segment
#Since we 'stitched' the sub-bands above, this consists
#of one spectrum covering the entire IF.
segments_s=sds_s.getSegmentIndices()
print segments_s
spectrum_s = sds_s.getPointSpectrum(0).getSegment(segments_s[0])

#subtract median of flux from stitched spectrum and inspect result
#in SpectrumExplorer
#note that as we modified the flux column in 'spectrum',
#it gets modified in 'sds_s' as well, since they are bound
#to the same object
flux=spectrum_s.getFlux()
mflux=MEDIAN(flux)
spectrum_s.setFlux(flux-mflux)

#To subtract this median flux from the un-stitched data,
#loop over segments the (i.e. WBS sub-bands) of this scan in this dataset
segments=sds.getSegmentIndices()
for p in segments:
 fluxp=sds.getPointSpectrum(0).getSegment(p).getFlux()-mflux
 sds.getPointSpectrum(0).getSegment(p).setFlux(fluxp)

#put modified dataset back in htp
htp.setProduct("60", DatasetWrapper(sds))

#The user should add a loop over the available SpectrumDatasets
#for all available data to be median-subtracted.

#some tips:

#Since 'htp' is bound to 'obs', 'obs' has been changed as well.
#If 'htp' were generated with a deep copy, e.g. htp1=htp.copy(),
#one would put it back in 'obs' as follows.
#obs.level["level2"].setProduct("WBS-H-USB",htp)

#In case a certain statement results in an error,
#a crash can be avoided using try/except
try:
 htp=obs.refs["level2"].product.refs["HRS-H-USB"].product
except:
 print "No HRS data available"

#scripts can be used in batch mode on the command line without
#starting the GUI, using 'jylaunch':
#$PATH/apps/jylaunch scriptname