# Title: OT1_Workshop_SPHOT_combine_obs
#
# Scope: this script loads multiple observations (in this particular case,
#        2 cross-scans parallel mode obs.), "gluing" level 1 product from each
#        obs. and then run the map-maker: the output will be a map for each
#        SPIRE channel containing data from all observations
#
# Luca Conversi - 2011/03/10


# List of observations IDs
obsids  = [1342183046, 1342183047]
pools   = ['spirephot_ngc6946_nom_rp', 'spirephot_ngc6946_orth_rp']
scanCon = 0

# Start cycle on observations
for i in range(len(obsids)):
	print "Processing ObsID =", obsids[i],"("+hex(obsids[i])+")"
	#
	# Load observations into variable obs
	obs = getObservation(obsids[i], poolName=pools[i])
	#
	# Diaply PSW map of the input obseravtions
	Display(obs.refs["level2"].product.refs["PLW"].product, title=hex(obs.obsid) + ' PLW Map')
	#
	if scanCon == 0:
		#
		# Copy level 1 product from 1st obs.
		scanCon = obs.level1
	else:
		for ref in obs.level1.refs:
			#
			# Attach level 1 product from other obs.
			scanCon.addRef(ref)


# Create a SpireListContext to be used as input of map-making
chanNum = obs.calibration.phot.chanNum
scans   = SpireListContext()

# Run baseline removal and populate the map-making input
for i in range(scanCon.count):
	psp = scanCon.getProduct(i)
	psp = removeBaseline(psp, chanNum=chanNum)
	scans.addProduct(psp)
	print "%2.2f" %(100*(i+1)/float(scanCon.count)) + "% done"

# Create combined map
plw_combined = naiveScanMapper(scans, array="PLW")
pmw_combined = naiveScanMapper(scans, array="PMW")
psw_combined = naiveScanMapper(scans, array="PSW")


# Display combined map
Display(plw_combined, title='Combined PLW Map')