Personal tools
You are here: Home 工作空间 Meetings and Conferences Hands on workshop during ISM2010 POF5 demo script
Document Actions

POF5 demo script

An error has been corrected. The variable poolName should be "0x50001833_1.2.0RC5".

Click here to get the file

Size 13.9 kB - File type text/python-source

File contents

#
# $Id: POF5_pipeline.py,v 1.70 2009/12/18 09:30:23 pasquale Exp $
#
# Official POF5 pipeline, i.e. the pipeline script for the 
# SPIRE Large Map observing mode.
#
# Author: Pasquale Panuzzo, CEA Saclay, France, pasquale.panuzzo@cea.fr
#
# History:
# 27 Jun 2006: First version
# 13 Jul 2006: (SG) Version for 1st DP End-to-End test
# 15 Sep 2006: (PP) Use EngConversion instead of DataProcessing
# 15 Nov 2006: (PP) Compatible with new calibration context
# 30 Nov 2006: (PP) Rename bolometer module into detresponse
# 16 Apr 2007: (PP) Moved in main
# 30 Apr 2007: add more comments
#  9 May 2007: engConversion output is now a EdpContext. Add map making (experimental)
# 29 May 2007: adapt to new Product API
# 30 May 2007: use new detector response tasks
# 20 Jun 2007: add obsLoader usage
# 12 Nov 2007: add level2 context
# 29 Jan 2008: remove detector response usage
# 20 Feb 2008: add non linearity correction task
# 17 Mar 2008: use the engConversion to process the complete level 0 context
# 10 Jun 2008: add SIAM product
# 30 Jun 2008: add ProductSink usage in non-interactive mode
# 12 Jul 2008: Fix usage of ProductSink (SPR-0772)
# 23 Jul 2008: Unload product reference (SPR-0745)
#  4 Sep 2008: Fix SPR-0858
#  7 Oct 2008: Introduce new photometer flux convertion task (SCR-0728). Simplify the syntax.
#  8 Oct 2008: Introduce electrical crosstalk correction (SCR-0886).
#  9 Oct 2008: Change bsmConverter in calcBsmAngles (SCR-0929).
#              Make pipeline script resistant to missing data (SPR-0930, SPR-0860, SPR-0866)
# 10 Oct 2008: New Level1Context (SCR-0770)
# 14 Oct 2008: Introduce optical crosstalk correction (SCR-0886). Set modelName of level1 and scan context (SPR-0598)
# 17 Oct 2008: Promote observations (SPR-0941)
# 21 Oct 2008: Introduce Electrical Filter Response correction (SCR-0957)
# 10 Nov 2008: Introduce Bolometer Time Response correction
# 12 Nov 2008: remove the usage of deprecated auxiliary context methods
# 16 Jan 2009: Better handling of broken observation (SPR-1039)
# 27 Jan 2009: Change usage of ObsLoader for SCR-1166. Fix imports. Add dummy browse product generation.
#  4 Mar 2009: Remove usage of ScanContext (SCR-1218). Set useSink to true (SCR-1229).
#              Use addRef(ProductRef) to attach saved products to Level1Context (SPR-1074)
#  6 Mar 2009: Fix Level1Context creation
# 10 Mar 2009: Add "level" parameter in pipeline signature to control the staring 
#              point of the pipeline [SCR-1261]
# 13 Mar 2009: Pass calibration products to corrElecFilterResponse and corrBolTimeResponse (SCR-1275, 1276)
# 21 Apr 2009: Use herschel.ia.task.MissingDataException [SPR-1353].
# 15 May 2009: Add progress setting [SCR-1260]
# 25 May 2009: Set deglitching values [SCR-1339]
# 29 May 2009: Use PhotOptCrossCorrection [SCR-1271]
# 12 Jun 2009: Select calibration products on bias and time basis. Pass ChanNoise to MADMap
# 31 Jul 2009: Extract the calibration and auxiliary products at the start of the script [SPIRE-SCR-1611]
#              Use CreateSpirePointingTask to generate pointing [SPIRE-SCR-1648]
#              Change name of Electrical crosstalk correction task [SPIRE-1649]
#  7 Aug 2009: Update task and parameters names of deglithcing.
# 11 Aug 2009: Change name of useSink into tempStorage [SPIRE-1646]
# 28 Aug 2009: Update parameter name of photOptCrossCorrection [SPIRE-1827]
#  5 Sep 2009: Update parameter name for photFluxConversion [SPIRE-1826] and temperatureDriftcorrection [SPIRE-1826]
#  8 Sep 2009: Add tasks to use turnaround data [SPIRE-1399]
# 10 Sep 2009: Fix turnaround retrival. Add removeBaseline task [SPIRE-1925]
#  5 Oct 2009: Add check to exclude BB gluing when a PCAL BB is in the middle [SPIRE-2020]
# 19 Oct 2009: Generate a proper browse product and browse image [SPIRE-2055]
# 18 Nov 2009: Fix browse product type.
# 27 Nov 2009: corrBolTimeResponse -> bolometerResponseCorrection
# 11 Dic 2009: Improve progress message [SPIRE-2231]
# 15 Dec 2009: Initialize the ProductSink with a TemporalPool [SPIRE-2095]. Update some comments.
# 18 Dec 2009: Fix import [SPIRE-2301]

# Import all needed classes
from herschel.spire.all import *
from herschel.ia.all import *
from herschel.ia.task.mode import *
from herschel.ia.pg import ProductSink
from java.lang import *
from java.util import *
from herschel.ia.obs.util import ObsParameter
from herschel.ia.pal.pool.lstore.util import TemporalPool

# Import the script tasks.py that contains the task definitions
from herschel.spire.ia.pipeline.scripts.POF5.POF5_tasks import *

# Input definition
from herschel.spire.ia.pipeline.scripts.POF5.POF5_input import *

# Import the script obsLoader.py that allows to load an ObservationContext from a storage.
from herschel.spire.ia.scripts.tools.obsLoader import *

# Open the input dialog to enter inputs
#inputs.openDialog()
#
import time
from time import localtime
print time.asctime(localtime())

poolName = "0x50001833_1.2.0RC5"
myobsid = 0x50001833   
storage = ProductStorage(poolName)
refs=storage.select(Query(ObservationContext, "obsid="+hex(myobsid)+"L"))
obs = refs[0].product  

#create a logger for the pipeline
logger=TaskModeManager.getMode().getLogger()

# Open a dialog to load the ObservationContext if "obs" is not defined.
try:
	obsid=obs.obsid
except NameError:
	loader=ObsLoader()
	obs=loader.getObs().product
pass

# Shift of time origin for plots
t0=obs.startDate.microsecondsSince1958()*1e-6

# get and print the obsid
obsid=obs.obsid
print "processing OBSID=", obsid,"("+hex(obsid)+")"

# Extract from the ObservationContext the calibration products that
# will be used in the script
bsmPos=obs.calibration.phot.bsmPos
lpfPar=obs.calibration.phot.lpfPar
detAngOff=obs.calibration.phot.detAngOff
elecCross=obs.calibration.phot.elecCross
optCross=obs.calibration.phot.optCross
chanTimeConst=obs.calibration.phot.chanTimeConst
chanNum=obs.calibration.phot.chanNum
fluxConvList=obs.calibration.phot.fluxConvList
tempDriftCorrList=obs.calibration.phot.tempDriftCorrList

# Extract from the observation context the auxiliary products that
# will be used in the script
hpp=obs.auxiliary.pointing
siam=obs.auxiliary.siam


# Set this to FALSE if you don't want to use the ProductSink
# and do all the processing in memory
tempStorage=Boolean.FALSE

# Initialize the ProductSink with a TemporalPool that will be removed when the
# HIPE session is closed, in case of interactive mode.
# The TemporalPool is created in a directory starting from the path defined by the
# var.hcss.workdir property. If this directory is inaccessible or not convenient, please
# change this property to a proper value.
if TaskModeManager.getType().toString() == "INTERACTIVE" and tempStorage:
    pname="tmp"+hex(System.currentTimeMillis())[2:-1]
    tmppool=TemporalPool.createTmpPool(pname,TemporalPool.CloseMode.DELETE_ON_CLOSE)
    ProductSink.getInstance().productStorage=ProductStorage(tmppool)


# From Level 0 to Level 0.5
if inputs.level=="level0":
	# Make Engineering conversion of level 0 products
	level0_5= engConversion(obs.level0,cal=obs.calibration, tempStorage=tempStorage)
	# Add the result to the observation in level 0.5
	obs.level0_5=level0_5
else:
	level0_5=obs.level0_5
pass

# set the progress
inputs.progress=20
# counter for computing progress
count=0

# From Level 0.5 to Level 1
if inputs.level=="level0" or inputs.level=="level0_5":
	# Create Level1 context
	level1=Level1Context(obsid)
	bbids=level0_5.getBbids(0xa103)
	nlines=len(bbids)
	print "number of scan lines:",nlines
	#
	# Loop over scan lines
	for bbid in bbids:
		block=level0_5.get(bbid)
		print "processing BBID="+hex(bbid)
		# Now move to engineering data products
		pdt  = block.pdt
		nhkt = block.nhkt
		if pdt == None:
			logger.severe("Building block "+hex(bbid)+" doesn't contain a PDT. Cannot process this building block.")
			print "Building block "+hex(bbid)+" doesn't contain a PDT. Cannot process this building block."
			continue
		if nhkt == None:
			logger.severe("Building block "+hex(bbid)+" doesn't contain a NHKT. Cannot process this building block.")
			print "Building block "+hex(bbid)+" doesn't contain a NHKT. Cannot process this building block."
			continue
		#
		# access and attach turnaround data to the nominal scan line
		bbCount=bbid & 0xFFFF
		pdtFollow=None
		nhktFollow=None
		pdtTrail=None
		nhktTrail=None
		if bbid < MAX(Long1d(bbids)):
			blockFollow=level0_5.get(0xaf000000L+bbCount)
			pdtFollow=blockFollow.pdt
			nhktFollow=blockFollow.nhkt
			if pdtFollow != None and pdtFollow.sampleTime[0] > pdt.sampleTime[-1]+3.0:
				pdtFollow=None
				nhktFollow=None
		if bbCount >1:
			blockTrail=level0_5.get(0xaf000000L+bbCount-1)
			pdtTrail=blockTrail.pdt
			nhktTrail=blockTrail.nhkt
			if pdtTrail != None and pdtTrail.sampleTime[-1] < pdt.sampleTime[0]-3.0:
				pdtTrail=None
				nhktTrail=None
		pdt=joinPhotDetTimelines(pdt,pdtTrail,pdtFollow)
		nhkt=joinNhkTimelines(nhkt,nhktTrail,nhktFollow)
		#
		# calculate BSM angles
		bat=calcBsmAngles(nhkt,bsmPos=bsmPos)
		#
		# create the SpirePointingProduct
		spp=createSpirePointing(detAngOff=detAngOff,bat=bat,hpp=hpp,siam=siam)
		#
		# run electrical crosstalk correction
		pdt=elecCrossCorrection(pdt,elecCross=elecCross)
		#
		# run the deglitch
		pdt=deglitchTimeline(pdt, scaleMin=1.0, scaleMax=8.0, scaleInterval=5, holderMin=-1.9,\
			holderMax=-0.3, correlationThreshold=0.69)
		#
		# run electrical Low Pass Filter response correction
		pdt=lpfResponseCorrection(pdt,lpfPar=lpfPar)
		#
		# run the flux conversion
		fluxConv=fluxConvList.getProduct(pdt.meta["biasMode"].value,pdt.startDate)
		pdt=photFluxConversion(pdt,fluxConv=fluxConv)
		#
		# run the temeperature drift correction
		tempDriftCorr=tempDriftCorrList.getProduct(pdt.meta["biasMode"].value,pdt.startDate)
		pdt=temperatureDriftCorrection(pdt,tempDriftCorr=tempDriftCorr)
		#
		# run bolometer time response correction
		pdt=bolometerResponseCorrection(pdt,chanTimeConst=chanTimeConst)
		#
		# run optical crosstalk correction
		pdt=photOptCrossCorrection(pdt,optCross=optCross)
		#
		# add pointing
		psp=associateSkyPosition(pdt,spp=spp)
		#
		# cut the timeline back to scan line range.
		# If you want include turnaround data in map making, call the following
		# task with the option "extend=True"
		psp=cutPhotDetTimelines(psp)
		#
		if inputs.plot:
			# Let's plot the signal of the central detector of PSW array.
			name="PSWE8"
			x=psp.sampleTime-t0  
			y=psp.getSignal(name)
			p=PlotXY(x,y,xtitle="Time [s]",ytitle="Signal [Jy]",titleText="Detector signal timeline")
			p[0].name=name
		# Store Photometer Scan Product in Level 1 product storage
		if tempStorage:
			ref=ProductSink.getInstance().save(psp)
			level1.addRef(ref)
		else:
			level1.addProduct(psp)
		#
		print "Completed BBID=0x%x (%i/%i)"%(bbid,count+1,nlines)
		# set the progress
		count=count+1
		inputs.progress = 20+(60*count)/nlines
	#
	if level1.count == 0:
		logger.severe("No scan line processed due to missing data. This observation CANNOT be processed!")
		print "No scan line processed due to missing data. This observation CANNOT be processed!"
		raise MissingDataException("No scan line processed due to missing data. This observation CANNOT be processed!")
	#
	obs.level1=level1
	# promote to LEVEL1_PROCESSED
	obs.obsState = ObservationContext.OBS_STATE_LEVEL1_PROCESSED
else:
	level1=obs.level1
pass

# Flag to switch on and off the baseline removal
useRemoveBaseline=True

# Create a SpireListContext to be used as input of map making
scans=SpireListContext()

# Run baseline removal and populate the map making input
for i in range(level1.count):
	if useRemoveBaseline:
		psp=level1.getProduct(i)
		psp=removeBaseline(psp,chanNum=chanNum)
		if tempStorage:
			ref=ProductSink.getInstance().save(psp)
			scans.addRef(ref)
		else:
			scans.addProduct(psp)
	else:
		scans.addRef(level1.refs[i])
pass

# Run mapmaking
if inputs.mapping == 'naive':
	mapPlw=naiveScanMapper(scans, array="PLW")
	inputs.progress=85
	mapPmw=naiveScanMapper(scans, array="PMW")
	inputs.progress=90
	mapPsw=naiveScanMapper(scans, array="PSW")
else:
	chanNoise=obs.calibration.phot.chanNoiseList.getProduct(level1.getProduct(0).meta["biasMode"].value,\
		level1.getProduct(0).startDate)
	mapPlw=madScanMapper(scans, array="PLW",chanNoise=chanNoise)
	inputs.progress=85
	mapPmw=madScanMapper(scans, array="PMW",chanNoise=chanNoise)
	inputs.progress=90
	mapPsw=madScanMapper(scans, array="PSW",chanNoise=chanNoise)
pass

# Create a context with level 2 products (maps) and attach it to the observation context
level2=MapContext()
level2.refs.put("PLW",ProductRef(mapPlw))
level2.refs.put("PMW",ProductRef(mapPmw))
level2.refs.put("PSW",ProductRef(mapPsw))
obs.level2=level2

# promote to LEVEL2_PROCESSED
obs.obsState = ObservationContext.OBS_STATE_LEVEL2_PROCESSED

# Create browse product and image
createRgbImage=CreateRgbImageTask()
browseProduct=createRgbImage(red=mapPlw,green=mapPmw,blue=mapPsw,percent=98.0,redFactor=1.0,\
	greenFactor=1.0,blueFactor=1.0)

# Populate metadata of the browse product
for par in ObsParameter.values():
	if obs.meta.containsKey(par.key) and par.key != "fileName":
		browseProduct.meta[par.key]=obs.meta[par.key].copy()
pass
browseProduct.startDate=obs.startDate
browseProduct.endDate=obs.endDate
browseProduct.instrument=obs.instrument
browseProduct.modelName=obs.modelName
browseProduct.description="Browse Product"
browseProduct.type="BROWSE"

# Attach the browse product to the ObservationContext
obs.browseProduct=browseProduct

# Generate the browse image
d=Display(True)
d.setImage(browseProduct)
obs.browseProductImage=d.renderedImage.asBufferedImage

print "completed OBSID=",obsid,"("+hex(obsid)+")"

try:
	loader.saveObs(obs)
except NameError:
	pass
	# do nothing
pass

inputs.progress=100

print time.asctime(localtime())

#### End of the script ####




Powered by Plone CMS, the Open Source Content Management System

This site conforms to the following standards: