### This script fetches SPIRE imaging data, using a list generated by the
### archive (in the CSV format), processes them to level 1 and saves them
### to disk in the correct format for later use by Scanamorphos.
### See important instructions below.


#######################################################

### modified pipeline script for the SPIRE Large Map observing mode.
### author: Pasquale Panuzzo, CEA Saclay, France, pasquale.panuzzo@cea.fr
### modified by Helene Roussel to be included in the Scanamorphos package

###  This file has been adapted from the POF5_pipeline.py script
###  distributed as part of the Herschel Common Science System (HCSS).
###  Copyright 2001-2010 Herschel Science Ground Segment Consortium

###  HCSS is free software: you can redistribute it and/or modify
###  it under the terms of the GNU Lesser General Public License as
###  published by the Free Software Foundation, either version 3 of
###  the License, or (at your option) any later version.

#######################################################

## Import classes and definitions:
from herschel.spire.all import *
from herschel.spire.util 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
from herschel.spire.ia.pipeline.scripts.POF5.POF5_tasks import *
from herschel.spire.ia.pipeline.scripts.POF5.POF5_input import *

fits = FitsArchive()

#######################################################


## local settings:

dir_root = "/ricci5nb/tshimizu/"
path = dir_root+"L1_HIPE13/"

### number of observations:
n_obs = 313

use_calib_tree = "spire_cal_13_1"
cal = spireCal(pool="spire_cal_13_1")
## observations:

# Load in all of the OBSIDs for the BAT AGN
f_obsid = open('/ricci9nb/tshimizu/Github/spire-catalog-photometry/bat_agn_spire_obsids.txt', 'r')
lines = f_obsid.readlines()[0]
lines = lines.split('\r')
list_obsids = [int(x.split('\t')[-1]) for x in lines]
list_names = [x.split('\t')[0] for x in lines]
#list_obsids = [1342245154]
#list_names = ['MCG-01-24-012']

source_prev = 'none'
vscan_prev = 0.
for i_obs in range(65, n_obs):
	##
	num_obsid = list_obsids[i_obs]
	source = list_names[i_obs]
	#source = str.lower(str(source))
	##
	if (source == source_prev):
		first_scan_num = first_scan_num + nscans
	else:
		first_scan_num = 1
	##
	##
	print ""
	print "Downloading obsid " + `num_obsid`
	obs = getObservation(num_obsid, useHsa=True, instrument="SPIRE", verbose=True)
	if obs.meta["obsMode"].value != "Parallel Mode":
		nrep = obs.meta["numRepetitions"].value
	else:
		nrep = 0
        vscan = obs.meta["scanSpeed"].value
	print source
	print "vscan = "+`vscan`
	print ""
	if ((source == source_prev) and (vscan != vscan_prev)): continue
	##
	path_processed = path+source+"_processed_subscans_"
	source_prev = source
	vscan_prev = vscan
	##
	##
	##
	calVersion = obs.calibration.version
	print "calibration version : "+calVersion
	if (calVersion != use_calib_tree):
		obs.calibration.update(cal)
	##
	##
	## Extract the calibration products from the observation context:
	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
	chanRelGains = obs.calibration.phot.chanRelGain
	fluxConvList = obs.calibration.phot.fluxConvList
	##
	## Extract the auxiliary products from the observation context:
	hpp=obs.auxiliary.pointing
	siam=obs.auxiliary.siam
	timeCorr = obs.auxiliary.timeCorrelation
	##
	##
	## start from level0:
	level0_5 = engConversion(obs.level0, cal=obs.calibration, tempStorage=Boolean.FALSE)
	##
	##
	##
	bbids = level0_5.getBbids(0xa103)
	nsubscans = len(bbids)
	if nrep == 0:
		nscans = 1
		nlegsperscan = nsubscans
	else:
		nscans = 2*nrep
		nlegsperscan = nsubscans/nscans
	print "number of scans: ", nscans
	print "number of subscans: ", nsubscans
	print "number of legs per scan = ", nlegsperscan
	##
	##
	## Create Level1 context:
	level1_allscans = Level1Context(num_obsid)
	##
	##
	##
	## Loop over scan legs:
	for bbid in bbids:
		##
		block=level0_5.get(bbid)
		## 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
		##
		## attach turnaround data:
		bbCount=bbid & 0xFFFF
		pdtLead=None
		nhktLead=None
		pdtTrail=None
		nhktTrail=None
		if bbCount > 1:
			blockLead=level0_5.get(0xaf000000L+bbCount-1)
			pdtLead=blockLead.pdt
			nhktLead=blockLead.nhkt
			if pdtLead != None and pdtLead.sampleTime[-1] < pdt.sampleTime[0]-3.0:
				pdtLead=None
				nhktLead=None
		if bbid < MAX(Long1d(bbids)):
			blockTrail=level0_5.get(0xaf000000L+bbCount)
			pdtTrail=blockTrail.pdt
			nhktTrail=blockTrail.nhkt
			if pdtTrail != None and pdtTrail.sampleTime[0] > pdt.sampleTime[-1]+3.0:
				pdtTrail=None
				nhktTrail=None
		pdt=joinPhotDetTimelines(pdt,pdtLead,pdtTrail)
		nhkt=joinNhkTimelines(nhkt,nhktLead,nhktTrail)
		##
		## BSM angles:
		bat=calcBsmAngles(nhkt,bsmPos=bsmPos)
		##
		## pointing product:
		spp=createSpirePointing(detAngOff=detAngOff,bat=bat,hpp=hpp,siam=siam)
		##
		## electrical crosstalk correction:
		pdt=elecCrossCorrection(pdt,elecCross=elecCross)
		##
		## deglitching:
		pdt=waveletDeglitcher(pdt, scaleMin=1.0, scaleMax=8.0, scaleInterval=5, holderMin=-1.9, holderMax=-0.3, correlationThreshold=0.69)
		##
		## electrical filter response correction:
		pdt=lpfResponseCorrection(pdt,lpfPar=lpfPar)
		##
		## flux conversion:
		fluxConv=fluxConvList.getProduct(pdt.meta["biasMode"].value,pdt.startDate)
		pdt=photFluxConversion(pdt,fluxConv=fluxConv)
		##
		## temperature drift correction:
		###	tempDriftCorr=tempDriftCorrList.getProduct(pdt.meta["biasMode"].value,pdt.startDate)
		###	pdt=temperatureDriftCorrection(pdt,tempDriftCorr=tempDriftCorr)
		##
		## bolometer time response correction:
		pdt=bolometerResponseCorrection(pdt,chanTimeConst=chanTimeConst)
		##
		## optical crosstalk correction:
		pdt=photOptCrossCorrection(pdt,optCross=optCross)
		##
		## add pointing:
		psp=associateSkyPosition(pdt,spp=spp)
		##
		## remove duplicate turnaround data:
		psp=cutPhotDetTimelines(psp,extend=True)
		##
        	## correct for any drifting of the onboard clock:
		psp=timeCorrelation(psp,timeCorr)
		##
        	## record the calibration tree version:
		psp.calVersion=calVersion
		##
		## store product in level 1 product storage:
        	level1_allscans.addProduct(psp)
		##
		print "Completed BBID="+hex(bbid)
		##
		##
		#######################################################
		##
		##
	print
	print "Apply relative gains for bolometers for better extended maps"
	for i in range(level1_allscans.getCount()):
		level1_allscans.getRefs().set(i,ProductRef(applyRelativeGains(level1_allscans.getProduct(i), gains = chanRelGains)))
	print "Finished applying relative gains"
	print
	## Save level 1 data as fits files:
	##
	if nlegsperscan > 99:
		prefix_1="00"
		prefix_2="0"
	else:
		prefix_1="0"
		prefix_2=""

	for iscan in range(nscans):
		scan = first_scan_num + iscan
		if scan < 10:
			str_scan = "0"+str(scan)
		else:
			str_scan = str(scan)
		root_filename = path_processed+"level1_scan"+str_scan
#
		for ileg in range(nlegsperscan):
			bbid = ileg+iscan*nlegsperscan
			product = level1_allscans.getProduct(bbid)
			subs = bbid+1
			if subs < 10:
				str_subs = prefix_1+str(subs)
			else:
				if subs < 100:
					str_subs = prefix_2+str(subs)
				else:
					str_subs = str(subs)
			filename = root_filename+"_leg"+str_subs+".fits"
			fits.save(filename,product)
			print filename


### END OF SCRIPT

#######################################################
