#   findSpurs_flagData_demo.py  :  Prepared for the ESAC HIFI OT1 DPWS
#   v0.0 11/Mar/2011 CM: created
#

# --- Find spurs, displaying flagged spectra ---
#
obs=getObservation("1342181161",poolName="1342181161_sscan_dbs_1b_r4")

#  In the ObservationContext tree inspect the datasets of one of the level 2
#  HTPs, you can scroll through them using the arrow down key.
#  When you find one with a spur in it, right click to open in Spectrum Explorer
#  Note that while spectrum explorer is open, you can display regions that are 
#  already flagged by right clicking on the window, then select "VIEW->FLAGS".
#  For this particular observation, you will note that the automatic spur finder 
#  already flagged these as bad.


# --- Flagging data ---
#
#  Documentation:  Chapter 5 of the HIFI Users Manual ("Flags in HIFI Data")
#
#  Data can be flagged using the Flagging task in the SpectrumExplorer or in the
#  command line. When working on the command line, the flagged data needs to be 
#  set back into the ObservationContext in order that other tasks can 
#  recognise the flags.


#  Extract one of the datasets and view it in SpectrumExplorer.
#  To be able to change the ObservationContext, we need to extract the HTP
#  first and then the dataset from that.

#  Extract the HTP for the WBS-H backend
htpwbsh=obs.refs["level2"].product.refs["WBS-H-USB"].product
# 3. Now get datset 32
ds = htpwbsh.get(32)

#  Double click on the 'ds' variable to open in SpectrumExplorer

#  To set flags in SpectrumExplorer:
#  - Right click on the plot and select Tools -> SelectPoints
#  - Draw a box around the region you want to flag (left click, drag).
#     Circles will appear around the selected points
#  - Right click on the selected points choose flag from the PointSelection menu
#  - Choose a manual flag, enter the flag value you want (try 21 to ignore data)
#  - Click OK
#  - If you now saved the observation context (right click on obs and send to
#    local store), the data would be stored flagged


# In the command line
#get the fourth sub-band
subband = 4
segment=ds.getPointSpectrum(0).getSegment(subband)
# define frequency range (in GHz) to be masked
freq1=589.6 
freq2=590.0
#make index of those frequencies
flag=segment.getFlag() # fetch the flags
freq=segment.getWave() # fetch the frequencies
index=freq.where((freq>=freq1).and(freq<=freq2))
#assign flag=1 (meaning bad pixel) to those indices
flagPixels(ds=ds, mask={(0,subband):index}, flag=21)
#check that it worked. Also one can View Flags in SpectrumExplorer
print ds.getFlag(4)

#  You should see that the flagged pixels have value 21, the start and end
#  pixels have value 4 and the reminder value 0


#  Now put the flagged dataset back into the HTP
htpwbsh.set("32", ds)
#  Put the HTP back into the Observation Context
obs.level["level2"].setProduct("WBS-H-USB",htpwbsh)