import os,sys,paircarsdata
from optparse import OptionParser
from casatools import msmetadata,table,measures,quanta,agentflagger,image,calibrater,ms
from casatasks import *
import logging,numpy as np,copy,glob,psutil,time,subprocess,getpass
from paircarstools.basic_func import *
from paircarstools.access_ms import *
from paircarstools.decor import *
from paircarstools.flagger import *
from astropy.io import fits
from pathlib import Path
from paircarstasks.run_paircars import *
from itertools import islice
os.system('rm -rf casa*log')

'''
Code is written by Devojyoti Kansabanik , 28 Jan, 2021
'''
global mainlog
# Logger initiating
###################
formatter = logging.Formatter('%(asctime)s %(levelname)-8s %(message)s',datefmt='%Y-%m-%d %H:%M:%S')
mainlog = logging.getLogger('paircars_main_log')
mainlog.setLevel(logging.DEBUG)
console=logging.StreamHandler(sys.stdout)
console.setFormatter(formatter)
mainlog.addHandler(console)

def main():
	usage= ' P-AIRCARS master controller for each day calibration'
	parser = OptionParser(usage=usage)
	parser.add_option('--basedir',dest="basedir",default=None,help="Name of base directory for a given day",metavar="Directory path")
	parser.add_option('--input_file',dest="input_file",default=None,help="Selfcal input file",metavar="Directory path")
	(options, args) = parser.parse_args()
	if options.basedir==None or os.path.exists(options.basedir)==False:
		print ('No base directory is given. Exiting...\n')
		os._exit(1)
	if options.input_file==None or os.path.exists(options.input_file)==False:
		print ('No input file is given.\n')
		os._exit(1)
	msg=paircars_master_controller(basedir=options.basedir,input_file=options.input_file)
	os._exit(msg)
	
def savelog(basedir=''):
	'''
	Save logger
	Parameters
	----------
	basedir : str
		Name of the basedir to save the logs
	'''
	if basedir!='':
		if os.path.isdir(basedir+'/Logs_and_Errors/Logs/')==False:
			os.makedirs(basedir+'/Logs_and_Errors/Logs/')
		filehandle=logging.FileHandler(basedir+'/Logs_and_Errors/Logs/P-AIRCARS_mainlog_'+str(os.path.basename(basedir))+'.log')
		filehandle.setFormatter(formatter)
		mainlog.addHandler(filehandle)
	mainlog.propagate = False
	return mainlog

def CPULIMIT_check():
	a=os.system('cpulimit -h > cpulimit_tmp')
	if a==256:
		os.system('rm -rf cpulimit_tmp')
		return 0
	else:
		os.system('rm -rf cpulimit_tmp')
		return 1


def paircars_master_controller(basedir='',input_file='',make_ds=True):
	'''
	P-AIRCARS master controller
	This function controls the parallelization in non-HPC platforms
	Parameters
	----------
	basedir : str
		Base directory name
	input_file : str
		Input file name
	'''
	# PAIRCARS master controller
	############################
	if basedir[-1]=='/':
		basedir=basedir[:-1]
	if basedir==None or os.path.isdir(basedir)==False:
		print ('Base directory does not exist.\n')
		return 1
	if input_file==None or os.path.exists(input_file)==False:
		print ('Input file does not exist.\n')
		return 1
	success=os.system('validating_paircars_input --basedir '+basedir+' --input_file '+input_file)
	if os.WEXITSTATUS(success)!=0:
		return 1
	os.chdir(basedir)
	sys.path.append(basedir)
	if os.path.exists(input_file)==False:
		print ('Input file does not exist.\n')
		return 1
	if input_file!=basedir+'/selfcal_inputs.py':
		os.system('cp -r '+input_file+' '+basedir+'/selfcal_inputs.py')
	import selfcal_inputs as inputs
	if os.path.isdir(inputs.basedir+'/data')==False:
		os.makedirs(inputs.basedir+'/data')
	
	bookkeep_file=inputs.basedir+'/DR_bookkeep.bkp'
	inpfil=open(basedir+'/selfcal_inputs.py','a')
	inpfil.write('bookkeeping_file=\''+str(bookkeep_file)+'\'\n')
	inpfil.close()
	cpu_list=np.load(basedir+'/cpulist.npy',allow_pickle=True).tolist()
	cpu_str=''
	for i in cpu_list:
		cpu_str+=str(i)+','
	cpu_str=cpu_str[:-1]
	max_allowed_jobs=len(cpu_list)
	
	if inputs.use_wsclean==True:
		wsclean_datadir=os.path.abspath(os.path.dirname(paircarsdata.__file__))
		wsclean_path=str(np.load(wsclean_datadir+'/wsclean_path.npy',allow_pickle=True))
		a=os.system(wsclean_path+'/wsclean > wsclean_test')
		if a!=0:
			print('WSClean is not installed. Using CASA for imaging.\n')
			use_wsclean=False
		else:
			use_wsclean=True
		os.system('rm -rf wsclean_test')
	else:
		use_wsclean=inputs.use_wsclean

	cpulimit_check=CPULIMIT_check()

	# Estimating total casa instances
	#################################
	available_paircars_instance=get_available_paircars_instance(inputs.paircars_dir,inputs.job_id,inputs.instance)

	# Removing screen logs
	######################
	if inputs.basedir[-1]=='/':
		inputs.basedir=inputs.basedir[:-1]
	paircars_base=os.path.dirname(inputs.basedir)
	if os.path.isdir(paircars_base+'/Logs_and_Errors/')==False:
		os.makedirs(paircars_base+'/Logs_and_Errors/')
	else:
		if os.path.isdir(paircars_base+'/Logs_and_Errors/Logs')==False:
			os.makedirs(paircars_base+'/Logs_and_Errors/Logs')
		else:
			log_list=glob.glob(paircars_base+'/Logs_and_Errors/Logs/*')
			if len(log_list)>0:
				for logfile in log_list:
					if 'Control_' in logfile and str(inputs.job_id) in logfile:
						pass
					else:
						os.system('rm -rf '+logfile)
		if os.path.isdir(paircars_base+'/Logs_and_Errors/Errors')==False:
			os.makedirs(paircars_base+'/Logs_and_Errors/Errors')
		else:
			error_list=glob.glob(paircars_base+'/Logs_and_Errors/Errors/*')
			if len(error_list)>0:
				for errorfile in error_list:
					if 'Control_' in errorfile and str(inputs.job_id) in errorfile:
						pass
					else:
						os.system('rm -rf '+errorfile)
	savelog(basedir=inputs.basedir)
	os.system('touch '+inputs.basedir+'/.paircars_running')
	is_internet=check_internet()

	# Deciding bandpass interval
	############################
	if inputs.quality_factor==0:
		bandpass_interval=0
		if inputs.safety_factor==0:
			calibrator_interval=90
		elif inputs.safety_factor==0:
			calibrator_interval=60
		elif inputs.safety_factor==0:
			calibrator_interval=30
	elif inputs.quality_factor==1:
		if inputs.safety_factor==0:
			bandpass_interval=7
			calibrator_interval=60
		elif inputs.safety_factor==1:
			bandpass_interval=5
			calibrator_interval=30
		elif inputs.safety_factor==2:
			bandpass_interval=3
			calibrator_interval=10
	elif inputs.quality_factor==2:
		if inputs.safety_factor==0:
			bandpass_interval=3
			calibrator_interval=20
		elif inputs.safety_factor==1:
			bandpass_interval=1
			calibrator_interval=10
		elif inputs.safety_factor==2:
			bandpass_interval=0.5
			calibrator_interval=4

	# User given timerange
	######################
	if inputs.timerange!='':
		time_list=[]
		for i in inputs.timerange.split(','):
			l=i.split('~')
			for j in l:
				time_list.append(j)
		timerange_list_mjdsecs=sorted([float("{:.2f}".format(timestamp_to_mjdsec(i,format=0))) for i in time_list])
		start_time_mjd=min(timerange_list_mjdsecs)
		end_time_mjd=max(timerange_list_mjdsecs)
	else:
		start_time_mjd=0
		end_time_mjd=0

	# Making local cal database folder
	##################################
	local_caldatabase=inputs.basedir+'/localdatabase'
	if os.path.isdir(local_caldatabase)==False:
		os.makedirs(local_caldatabase)
	mainlog.info('Local calibration database is at : '+local_caldatabase+'\n')

	# Download available cal files from paircars database #TODO
	#####################################################


	# Organising ms
	###############
	updated_MWA_obsids=0
	mainlog.info('Organising measurement sets .....\n')
	measurement_set_list=glob.glob(inputs.basedir+'/data/*.ms')
	# Making dynamic spectrum
	#########################
	'''for msname in measurement_set_list:
		ds_base_name='Dynamic_Spectrum'
		dsloc=''
		if make_ds:
			if os.path.exists(str(inputs.savedir)+'/dynamic_spectra/'+ds_base_name+'_'+os.path.basename(msname).split('.ms')[0]+'.npy')==False or \
				os.path.exists(str(inputs.savedir)+'/dynamic_spectra/'+ds_base_name+'_TB_'+os.path.basename(msname).split('.ms')[0]+'.png')==False or \
				os.path.exists(str(inputs.savedir)+'/dynamic_spectra/'+ds_base_name+'_FLUX_'+os.path.basename(msname).split('.ms')[0]+'.png')==False:
				cmd='make_solar_dynamic_spectrum --mslist '+msname+' --metafits '+metafits+' --num_threads 4 --outfile '+ds_base_name+' --outdir '+str(inputs.savedir)+\
					'/dynamic_spectra --plot True --plot_outtype png '+' --verbose '+str(verbose)+' --plotdir '+str(inputs.savedir)+'/dynamic_spectra'
				ms_mainlog.info(cmd+'\n')
				os.system(cmd)
			else:
				ms_mainlog.info('Dynamic spectra are already made.\n')
				ms_mainlog.info('Dynamic spectra plots are at : '+str(inputs.savedir)+'/dynamic_spectra/'+ds_base_name+'_TB_'+os.path.basename(msname).split('.ms')[0]+'.png and '+\
					str(inputs.savedir)+'/dynamic_spectra/'+ds_base_name+'_TB_'+os.path.basename(msname).split('.ms')[0]+'.png\n')
			dsloc=str(inputs.savedir)+'/dynamic_spectra/'+ds_base_name+'_'+os.path.basename(msname).split('.ms')[0]+'.npy'''


	# Spliting measurement set into 2 coarse channel blocks
	#######################################################
	for ms in measurement_set_list:
		if inputs.quality_factor==0:
			max_allowed_bw=5.12
		elif inputs.quality_factor==1:
			max_allowed_bw=2.56
		else:
			max_allowed_bw=1.28
		if os.path.exists(inputs.basedir+'/data/.splited_'+os.path.basename(ms).split('.ms')[0])==False:		
			AM=AccessMS(ms)
			full_bw=AM.calc_bandwidth()/10**6
			total_time=round(AM.calc_total_time(),2)
			if full_bw>max_allowed_bw:
				num_of_trials=0
				while num_of_trials<5:
					ms_split_cmd=['--msname '+ms,'--freq_interval '+str(max_allowed_bw*1000),'--freq_width '+str(max_allowed_bw*1000),'--time_width '+str(total_time),\
							'--time_interval '+str(total_time),'--force True','--savedir '+inputs.basedir+'/data','--max_flag_frac 1.0','--verbose '+str(inputs.verbose)]
					cmd='parallel_ms_split '+' '.join(ms_split_cmd)
					check_thread_availability(int(psutil.cpu_count()*inputs.cpu_frac*0.5))
					mainlog.info(cmd+'\n')
					a=os.system(cmd)
					if a==0:
						os.system('touch '+inputs.basedir+'/data/.splited_'+os.path.basename(ms).split('.ms')[0])
						splited_ms=glob.glob(inputs.basedir+'/data/splited_ms/'+os.path.basename(ms).split('.ms')[0]+'/*.ms')
						for i  in splited_ms:
							os.system('mv '+i+' '+inputs.basedir+'/data/')
						os.system('rm -rf '+inputs.basedir+'/data/splited_ms')
						os.unlink(ms)
						break
					else:
						mainlog.warning('Error orccured during spliting the ms. Try number : '+str(num_of_trials)+'\n')
						os.system('rm -rf '+inputs.basedir+'/data/splited_ms')
						num_of_trials+=1
						continue
				try:
					os.unlink(ms)
				except:
					pass
	measurement_set_list=glob.glob(inputs.basedir+'/data/*.ms')
	for ms in measurement_set_list:
		if os.path.isdir(os.path.realpath(ms))==False:
			os.system('unlink '+ms)
			measurement_set_list.remove(ms)
	msfreqs=[float(os.path.basename(a).split('.ms')[0].split('_')[-1]) for a in measurement_set_list]
	mstimes=[float(timestamp_to_mjdsec('/'.join(os.path.basename(a).split('time_')[-1].split('_freq')[0].split('_')[:3])+'/'+\
				':'.join(os.path.basename(a).split('time_')[-1].split('_freq')[0].split('_')[3:]))) for a in measurement_set_list]
	mstimes_iso=[('-'.join(os.path.basename(a).split('.ms')[0].split('_')[1:4])+' '+':'.join(os.path.basename(a).split('.ms')[0].split('_')[4:7])) for a in measurement_set_list]
	ms_OBSIDs=[]
	metafits_obsids_msdir=glob.glob(inputs.msdir+'/*.metafits')

	if len(metafits_obsids_msdir)!=0:
		for metafits in metafits_obsids_msdir:
			os.system('cp -r '+metafits+' '+inputs.basedir+'/data/'+os.path.basename(metafits))
	metafits_obsids=[int(os.path.basename(x).split('.metafits')[0]) for x in glob.glob(inputs.basedir+'/data/*.metafits')]

	if len(measurement_set_list)==0:
		mainlog.info('No measurement set present for the given time range. Please check your time range or the data directory if measurement set exists or not.\n')
		os._exit(0)
	measurement_set_list_copy=copy.deepcopy(measurement_set_list)
	for i in range(len(measurement_set_list_copy)):
		msname=measurement_set_list_copy[i]
		AMtimerange=AccessMS(msname)
		mjdstamps=AMtimerange.get_timestamps_in_mjdsecs()[0]
		start_mjdtimestamp=min(mjdstamps)
		end_mjdtimestamp=max(mjdstamps)
		allow_ms=False
		if start_time_mjd!=0 and end_time_mjd!=0:
			if start_mjdtimestamp>end_time_mjd or end_mjdtimestamp<start_time_mjd:
				allow_ms=False
			else:
				allow_ms=True
		else:
			allow_ms=True
		if allow_ms==True:
			obsid=get_OBSID_from_ms(msname)
			if obsid==0 or len(metafits_obsids)==0:
				try:
					GPStime=int(Time(mstimes_iso[i],format='iso',scale='utc').gps)
					diff_gpstime=[]
					for a in metafits_obsids:
						if abs(a-GPStime)<480:
							diff_gpstime.append(a)
					if len(diff_gpstime)!=0:
						obsid=np.min(np.array(diff_gpstime))
						ms_OBSIDs.append(obsid)
					elif is_internet==False:
						mainlog.info('No local metafits file is present and no internet connection is available to download metafits file.'+\
										' Either supply metafits file in Data Directory or connect to the internet. Exiting P-AIRCARS....\n')
						os._exit(0)
					else:
						if os.path.exists(inputs.basedir+'/data/'+str(obsid)+'.metafits')==False:
							mainlog.info('Trying to download unavailable metafits for OBS ID :'+str(obsid)+' at : '+inputs.basedir+'/data/'+str(obsid)+'.metafits.\n')
							metafits=download_metafits(msname,inputs.basedir+'/data')
						else:
							metafits=inputs.basedir+'/data/'+str(obsid)+'.metafits'
						if metafits!=None:
							mainlog.info('Metafits file downloaded at : '+metafits+'\n')
							ms_OBSIDs.append(get_OBSID_from_metafits(metafits))
						else:
							mainlog.info('Metafits file could not be downloaded for ms : '+msname+'. Removing ms from list.\n')
							measurement_set_list.remove(msname)
							msfreqs.remove(msfreqs[i])
							mstimes.remove(mstimes[i])
							mstimes_iso.remove(mstimes_iso[i]) 
				except:
					mainlog.info('Could not connect to MWA metadata server. No metafits files are found in local data directory. Exiting P-AIRCARS.....\n')
					os._exit(0)
			else:
				ms_OBSIDs.append(obsid)
		else:
			measurement_set_list.remove(msname)
	if len(measurement_set_list)==0:
		mainlog.info('No measurement set present for the given time range. Please check your time range or the data directory if measurement set exists or not.\n')
		os._exit(0)

	if inputs.calc_selfcalib_params==True:
		if inputs.quality_factor==0:
			skip_freq_pol=1280
			if inputs.safety_factor==0:
				skip_time		=	25
				skip_freq		=	1280
			elif inputs.safety_factor==1:
				skip_time		=	20
				skip_freq		=	1280
			else:
				skip_time		=	15		
				skip_freq		=	1280
		elif inputs.quality_factor==1:
			skip_freq_pol=1280
			if inputs.safety_factor==0:
				skip_time		=	20
				skip_freq		=	1280
			elif inputs.safety_factor==1:
				skip_time		=	15
				skip_freq		=	1280
			else:
				skip_time		=	10
				skip_freq		=	1280
		else:
			skip_freq_pol=1280
			if inputs.safety_factor==0:
				skip_time		=	15
				skip_freq		=	1280
			elif inputs.safety_factor==1:
				skip_time		=	10
				skip_freq		=	1280
			else:
				skip_time		=	5
				skip_freq		=   1280
	else:
		skip_time=inputs.skip_time
		skip_freq=inputs.skip_freq
		if skip_freq>1280:
			mainlog.info('Skip frequency greater than 1.28 MHz is not appropriate for bandpass self-calibration. Using 1.28 MHz for bandpass.\n')
			skip_freq=1280
			inputs.skip_freq=1280


	# Making metafits dictionary
	############################
	ms_gridpoints=[]
	ms_OBSIDs_gcal=[]
	ms_OBSIDs_bcal=[]
	ms_OBSIDs_pcal=[]
	final_image_dic={}
	if len(measurement_set_list)!=0:
		metafits_dic={}
		for i in range(len(measurement_set_list)):
			msname=measurement_set_list[i]
			metafits=str(ms_OBSIDs[i])+'.metafits'
			if os.path.isfile(inputs.basedir+'/data/'+metafits)==False: 
				mainlog.info('Trying to download unavailable metafits for OBS ID :'+str(ms_OBSIDs[i])+' at : '+inputs.basedir+'/data/'+str(ms_OBSIDs[i])+'.metafits.\n')
				downloaded_metafits=download_metafits(msname,inputs.basedir+'/data')
				if downloaded_metafits!=None:
					metafits=os.path.basename(downloaded_metafits)				
					mainlog.info('Metafits file downloaded at : '+downloaded_metafits+'\n')
					metafits_dic[msname.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]=inputs.basedir+'/data/'+metafits
					ms_gridpoints.append(fits.getheader(inputs.basedir+'/data/'+metafits)['GRIDNUM'])
				else:
					mainlog.info('Metafits file could not be downloaded for ms : '+msname+' and also not found locally at : '+inputs.basedir+'/data/'+metafits+\
								' for ms : '+msname+'\n')	
					measurement_set_list.remove(msname)
					msfreqs.remove(msfreqs[i])
					mstimes.remove(mstimes[i])
					mstimes_iso.remove(mstimes_iso[i])
					ms_OBSIDs.remove(ms_OBSIDs[i])
			else:
				mainlog.info('Metafits file found at : '+inputs.basedir+'/data/'+metafits+'\n')
				metafits_dic[msname.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]=inputs.basedir+'/data/'+metafits
				ms_gridpoints.append(fits.getheader(inputs.basedir+'/data/'+metafits)['GRIDNUM'])

		ms_list_copy=copy.deepcopy(measurement_set_list)
		ms_OBSIDs_copy=copy.deepcopy(ms_OBSIDs)
		metafits_dic_copy=copy.deepcopy(metafits_dic)
		ref_timechan_success=False

		ms_list_cals=[]
		ms_OBSIDs_cal=[]
		msfreqs_cal=[]
		mstimes_cal=[]
		for i in range(len(ms_list_copy)):
			if i==0:
				ms_list_cals.append(ms_list_copy[i])
				ms_OBSIDs_cal.append(ms_OBSIDs_copy[i])
				msfreqs_cal.append(msfreqs[i])
				mstimes_cal.append(mstimes[i])
				if inputs.do_polcal==True:
					ms=ms_list_copy[i]
					grid0=ms_gridpoints[i]
			else:
				if inputs.do_polcal==True:
					ms=ms_list_copy[i]
					meta=metafits_dic[ms.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]
					grid1=ms_gridpoints[i]
				time_diff=np.array(abs(mstimes[i]-np.array(mstimes_cal)))
				freq_diff=np.array(abs(msfreqs[i]-np.array(msfreqs_cal)))
				if (((len(time_diff)==0 and len(freq_diff)==0) or grid0!=grid1) and inputs.do_polcal==True)	or (len(time_diff)==0 and len(freq_diff)==0 and inputs.do_polcal==False):
					ms_list_cals.append(ms_list_copy[i])
					ms_OBSIDs_cal.append(ms_OBSIDs_copy[i])
					msfreqs_cal.append(msfreqs[i])
					mstimes_cal.append(mstimes[i])
				elif (((min(time_diff)>skip_time or min(freq_diff)>(skip_freq/10**3)) or grid0!=grid1) and inputs.do_polcal==True)	or\
					 (min(time_diff)>skip_time or min(freq_diff)>(skip_freq/10**3) and inputs.do_polcal==False):
					ms_list_cals.append(ms_list_copy[i])
					ms_OBSIDs_cal.append(ms_OBSIDs_copy[i])
					msfreqs_cal.append(msfreqs[i])
					mstimes_cal.append(mstimes[i])
			
		msfreqs_cal=sorted(msfreqs_cal)
		mstimes_cal=sorted(mstimes_cal)
		ms_list_copy_cals=copy.deepcopy(ms_list_cals)
		ms_OBSIDs_copy_cals=copy.deepcopy(ms_OBSIDs_cal)
		calibrator_found=False # TODO : change after done this part
		gaincal_list=[]
		bandpass_list=[]
		caltable_list=[]
		spawned_ms_jobs={}
		finished_ms=[]
		spawned_casa_instances=0
		spawned_ms=[]
		pol_ref_time=None
		if inputs.only_imaging==False: # If calibration is required
			mainlog.info('Measurement sets to perform calibrations : \n')
			for i in ms_list_cals:
				mainlog.info(i+'\n')

			screen_list=[os.path.basename(i) for i in glob.glob('/var/run/screen/S-'+str(getpass.getuser())+'/*')]
			delete_screen_list=[]
			for j in screen_list:
				if str(inputs.job_id) in j and 'P-AIRCARS_mainlog' not in j and 'screen_inputs.basedir_for' not in j:
					delete_screen_list.append(j) 
			for i in delete_screen_list:
				if 'remote' not in i and 'carta' not in i and 'jprq' not in i:
					os.system('screen -S '+i+' -X quit')

			ref_timechan_done_list=glob.glob(inputs.basedir+'/.ref_timechan_done_*')
			for i in ref_timechan_done_list:
				try:
					ref_timechan_done_msg=int(i.split('_')[-1])
					if ref_timechan_done_msg!=0:
						os.system('rm -rf '+i)
						ref_timechan_done_list.remove(i)
				except:
					os.system('rm -rf '+i)
					ref_timechan_done_list.remove(i)
			if len(ref_timechan_done_list)>=2 and os.path.exists(inputs.basedir+'/Ref_time_freq_slice_output.npy'):
				result_dic=np.load(inputs.basedir+'/Ref_time_freq_slice_output.npy',allow_pickle=True).all()
				result=list(result_dic.values())[0]
				common_ref_time=result[7]
				mainlog.info('Reference time frequency calibration has been done already.\n')
				ref_timechan_done_list=glob.glob(inputs.basedir+'/.ref_timechan_done_*')
				ref_freq_time_msname=[inputs.basedir+'/data/time_'+(a.split('time_')[-1].split('.ms_')[0])+'.ms' for a in ref_timechan_done_list]
				all_success_list=[]
				for i in ref_freq_time_msname:
					if i not in measurement_set_list:
						ref_freq_time_msname.remove(i)
				for i in measurement_set_list:
					if i in ref_freq_time_msname:
						all_success_list.append(i)
						spawned_ms.append(i)
				mainlog.info('Reference time frequency measurement set list is : '+str(ref_freq_time_msname)+'\n')
				mainlog.info('######################\n')
				ref_ms_OBSID_list=[]
				ref_time_freq_metafits=[]
				for i in ref_freq_time_msname:
					reftimefreq_ms_OBSID=get_OBSID_from_metafits(metafits_dic[i.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]])
					ref_ms_OBSID_list.append(reftimefreq_ms_OBSID)
					ref_time_freq_metafits.append(metafits_dic[i.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]])
				ref_timechan_success=True
				c=0
			############################################## # TODO : Do not apply, only make list
			# Applying solutions from calibration database




			# If no calibration is found in calibration database
			# Apply calibrator solution
			####################################################


			##########################################
			if calibrator_found==True and len(gaincal_list)!=0:
				do_fluxcal=False
			else:
				do_fluxcal=True
			while ref_timechan_success==False:
				if len(msfreqs)==0 or len(mstimes)==0:
					mainlog.info('No measurement set is present for performing reference time frequency calibration.\n')
					os._exit(1)
				if min(msfreqs_cal)>170:
					if len(msfreqs_cal)>=3:
						ref_freqstamp=[str(msfreqs_cal[0]),str(msfreqs_cal[int(len(msfreqs_cal)/2)]),str(msfreqs_cal[-1])]
					elif len(msfreqs_cal)==2:
						ref_freqstamp=[str(msfreqs_cal[0]),str(msfreqs_cal[-1])]
					else:
						ref_freqstamp=[str(msfreqs_cal[int(len(msfreqs_cal)/2)])]
				else:
					sub_msfreqs_cal=[]
					for i in msfreqs_cal:
						if i<170:
							sub_msfreqs_cal.append(i)
					if len(sub_msfreqs_cal)>=3:
						ref_freqstamp=[str(sub_msfreqs_cal[0]),str(sub_msfreqs_cal[int(len(sub_msfreqs_cal)/2)]),str(sub_msfreqs_cal[-1])]
					elif len(sub_msfreqs_cal)==2:
						ref_freqstamp=[str(sub_msfreqs_cal[0]),str(sub_msfreqs_cal[-1])]
					else:
						ref_freqstamp=[str(sub_msfreqs_cal[int(len(sub_msfreqs_cal)/2)])]
				ref_timestamp=mjdsec_to_timestamp(mstimes_cal[int(len(mstimes_cal)/2)],format=3)
				ref_timestamp_mjd=mstimes_cal[int(len(mstimes_cal)/2)]
				# Selecting reference time frequnency ms
				######################################
				ref_freq_time_msname=[]
				ref_time_freq_metafits=[]
				ref_ms_OBSID_list=[]

				for i in range(len(measurement_set_list)):
					msname=measurement_set_list[i]
					for j in ref_freqstamp:
						if j in msname and ref_timestamp in msname:
							ref_freq_time_msname.append(msname)
							ref_time_freq_metafits.append(metafits_dic[msname.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]])
				mainlog.info('Reference time frequency measurement set list is : '+str(ref_freq_time_msname)+'\n')
				mainlog.info('######################\n')

				# Self calibration for reference time frequency ms
				################################################## 
				for i in ref_freq_time_msname:
					reftimefreq_ms_OBSID=get_OBSID_from_metafits(metafits_dic[i.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]])
					ref_ms_OBSID_list.append(reftimefreq_ms_OBSID)
					AM=AccessMS(i)
					workdir=inputs.basedir+'/'+os.path.basename(i).split('.ms')[0]
					if os.path.isdir(workdir)==False:
						os.makedirs(workdir)
				
					touch_file_list=glob.glob(inputs.basedir+'/.ref_timechan_done_*'+str(reftimefreq_ms_OBSID)+'*')
					if len(touch_file_list)!=0:
						for t in touch_file_list:
							msg=t.split('_')[-1]
							try:
								msg=int(msg)
							except:
								pass
							if type(msg)==str:
								os.system('rm -rf '+t)
							else:
								msg=int(msg)
								if msg>100:
									msg-=100
								if msg!=0 and msg!=8 and msg!=9:
									os.system('rm -rf '+t)

					touch_file_list1=glob.glob(inputs.basedir+'/.Finished_*cal*'+str(reftimefreq_ms_OBSID)+'*')
					if len(touch_file_list1)!=0:
						for t in touch_file_list1:
							msg=t.split('_')[-1]
							try:
								msg=int(msg)
							except:
								pass
							if type(msg)==str:
								os.system('rm -rf '+t)
							else:
								msg=int(msg)
								if msg>100:
									msg-=100
								if msg!=0 and msg!=8 and msg!=9:
									os.system('rm -rf '+t)
				if calibrator_found==True and len(gaincal_list)!=0: # If calibration found from paircars database or calibrator
					for j in range(len(ref_freq_time_msname)):
						calibrator_OBSID=np.min(np.array([x.split('_')[0] for x in gaincal_list])) # Considering naming convention OBSID_COARSECHAN.gcal
						# Not debugged this part
						ms_gaincal_list=[]
						ms_bandpass_list=[]
						for gcal in gaincal_list:
							if int(gcal.split('_')[0])==int(calibrator_OBSID):
								ms_gaincal_list.append(gcal)
						for bcal in bandpass_list:
							if int(bcal.split('_')[0])==int(calibrator_OBSID):
								ms_bandpass_list.append(bcal)
						bandpass_freq_list=[float(get_caltable_metadata(i)['Channel 0 frequency (MHz)']) for i in ms_bandpass_list]
						bandpass_coarse_chan=np.array([freq_to_MWA_coarse(freq) for freq in bandpass_freq_list])
						reftimefreq_ms_OBSID=ref_ms_OBSID_list[int(len(ref_ms_OBSID_list)/2)]
						single_ref_freq_time_metafits=ref_time_freq_metafits[int(len(ref_time_freq_metafits)/2)]
						ref_freqstamp=ref_freqstamp[int(len(ref_freqstamp)/2)]
						single_ref_freq_time_ms=ref_freq_time_msname[j]
						AMref=AccessMS(single_ref_freq_time_ms)
						ms_freq=float(AMref.calc_meanfreq()/10**6)
						ms_coarse_chan=freq_to_MWA_coarse(ms_freq)
						ms_bandpass_list=ms_bandpass_list[np.where(bandpass_coarse_chan==ms_coarse_chan)]
						caltable_list=ms_gaincal_list+ms_bandpass_list
						#
						if abs(reftimefreq_ms_OBSID-calibrator_OBSID)>calibrator_interval*60:
							ref_time_freq=True
						else:
							ref_time_freq=False
						workdir=inputs.basedir+'/'+os.path.basename(single_ref_freq_time_ms).split('.ms')[0]
						if os.path.isdir(workdir)==False:
							if os.path.exists(workdir)==True:
								os.system('rm -rf '+workdir)
							os.makedirs(workdir)
						os.system('cp -r selfcal_inputs.py '+workdir+'/selfcal_inputs.py')
						if os.path.isfile(single_ref_freq_time_ms+'/copy.table.lock')==False:
							os.system('cp -r '+single_ref_freq_time_ms+'/table.lock '+single_ref_freq_time_ms+'/copy.table.lock')
						if os.path.isfile(single_ref_freq_time_ms+'/copy.table.info')==False:
							os.system('cp -r '+single_ref_freq_time_ms+'/table.info '+single_ref_freq_time_ms+'/copy.table.info')
						ref_time_avg=0.5
						ref_time_bin=0.5
						cpu_str=''
						for i in cpu_list:
							cpu_str+=str(i)+','
						cpu_str=cpu_str[:-1]
						cmd='run_paircars --msname '+single_ref_freq_time_ms+' --metafits '+single_ref_freq_time_metafits+' --basedir '+inputs.basedir+' --workdir '+workdir+\
							' --ref_freq_avg '+str(40)+' --ref_time_avg '+str(ref_time_avg)+' --ref_time_freq True --do_bandpass '+str(inputs.do_bandpass)+' --do_polcal '+str(inputs.do_polcal)+\
							' --cal_attenuation '+str(1.0)+' --scratch True --caltables '+str(','.join(caltable_list))+' --wsclean '+str(use_wsclean)+' --cpu_frac '+str(inputs.cpu_frac)+\
							' --do_fluxcal '+str(do_fluxcal)+' --cubical '+str(inputs.use_cubical)+' --verbose '+str(inputs.verbose)+' --input_file '+input_file+' --do_pbcor False'+\
							' --ref_time_bin '+str(ref_time_bin)+' --cpu_list '+cpu_str+' --pol_ref_time '+str(pol_ref_time)
						if inputs.ref_time!='':				
							cmd+=' --ref_time '+str(inputs.ref_time) 
						else:
							ds_base_name='Dynamic_Spectrum.npy'
							dsfile=workdir+'/'+ds_base_name
							savefile=workdir+'/ds_data'
							if j==0:
								if inputs.timerange=='':
									temp_timerange=''
								else:
									temp_timerange=inputs.timerange
								ms_mainlog.info('get_ref_time_fromDS(\''+single_ref_freq_time_ms+'\',\''+single_ref_freq_time_metafits+'\',nthreads=3,verbose='+str(inputs.verbose)+\
								',force_calculate=True,dsfile=\''+dsfile+'\',caltables=[],DRbased=False,T_limit=1,timerange=\''+str(temp_timerange)+\
								'\',docal=True,ref_ant='+str(inputs.ref_ant)+',cpu_list='+str(cpu_list)+',cpu_frac='+str(inputs.cpu_frac)+',save_file=\''+savefile+'\')\n')
								ref_time,ref_time_index,ref_timestamps,T_timeseries=get_ref_time_fromDS(single_ref_freq_time_ms,single_ref_freq_time_metafits,nthreads=3,verbose=inputs.verbose,force_calculate=True,\
									dsfile=dsfile,caltables=[],DRbased=False,T_limit=1,timerange=temp_timerange,docal=True,ref_ant=inputs.ref_ant,cpu_list=cpu_list,cpu_frac=inputs.cpu_frac,save_file=savefile)	
							else:
								ms_mainlog.info('get_ref_time_fromDS(\''+single_ref_freq_time_ms+'\',\''+single_ref_freq_time_metafits+'\',nthreads=3,verbose='+str(inputs.verbose)+\
								',input_pol_ref_time=\''+str(ref_time)+'\',force_calculate=True,dsfile=\''+dsfile+'\',caltables=[],DRbased=False,T_limit=1,timerange=\''+str(temp_timerange)+\
								'\',docal=True,ref_ant='+str(inputs.ref_ant)+',cpu_list='+str(cpu_list)+',cpu_frac='+str(inputs.cpu_frac)+',save_file=\''+savefile+'\')\n')
								ref_time,ref_time_index,ref_timestamps,T_timeseries=get_ref_time_fromDS(single_ref_freq_time_ms,single_ref_freq_time_metafits,nthreads=3,verbose=inputs.verbose,\
								input_pol_ref_time=ref_time,force_calculate=True,dsfile=dsfile,caltables=[],DRbased=False,T_limit=1,timerange=temp_timerange,docal=True,ref_ant=inputs.ref_ant,cpu_list=cpu_list,\
								cpu_frac=inputs.cpu_frac,save_file=savefile)
							cmd+=' --ref_time '+str(ref_time)
									 # TODO : Change calibrator attenuator
						screen_name=str(reftimefreq_ms_OBSID)+'_'+str(os.path.basename(ref_freq_time_msname).split('.ms')[0])+'_runpaircars_'+str(inputs.job_id)
						finished_touch_file=inputs.basedir+'/.Finished_runpaircars_'+str(reftimefFcur_cpu_listreq_ms_OBSID)+'_'+str(single_ref_freq_time_ms.split('.ms')[0])
						screen_batch_file=paircars_instance_runner(cmd,inputs.basedir,inputs.paircars_dir,screen_name,finished_touch_file,inputs.job_id)
						mainlog.info('bash '+screen_batch_file+'\n')
						os.system('bash '+screen_batch_file)
						time.sleep(1.0)
						while True:
							time.sleep(5)
							available_paircars_instance=get_available_paircars_instance(inputs.paircars_dir,inputs.job_id,inputs.instance)
							if available_paircars_instance>1:
								break
					mainlog.info('Waiting for finishing calibration of reference time frequency......\n')
					mainlog.info('Total reference ms : '+str(len(ref_freq_time_msname))+'\n')
					while True:
						all_touch_files=[]
						all_success_list=[]
						old_finished=0
						c=0
						for j in range(len(ref_freq_time_msname)):
							single_reftimefreq_ms_OBSID=ref_ms_OBSID_list[j]	
							single_ref_freq_time_ms=ref_freq_time_msname[j]
							single_ref_freqstamp=ref_freqstamp[j]
							touch_file_list=glob.glob(inputs.basedir+'/.ref_timechan_done_'+str(single_reftimefreq_ms_OBSID)+'_*')
							runtime_fail=False
							selfcal_fail=False
							other_error=False
							if len(touch_file_list)>0:
								for j in touch_file_list:
									if j not in all_touch_files:
										all_touch_files.append(j)
								i=touch_file_list[0]
								if '_selfcalerror' in i:
									selfcal_fail=True
								elif '_runtimeerror' in i:
									runtime_fail=True
								elif '_error' in i:
									other_error=True
								else:
									all_success_list.append(single_ref_freq_time_ms)
								if selfcal_fail==True:
									mainlog.info('Reference time frequency calibration failed. Trying with new reference time channel.\n')
									msfreqs_cal.remove(float(single_ref_freqstamp))
									break
								elif runtime_fail==True:
									mainlog.info('Reference time frequency calibration failed during run time. Some error occured during runtime. Contact developer to fix the problem.\n')
									msfreqs_cal.remove(float(single_ref_freqstamp))
									mstimes_cal.remove(float(ref_timestamp_mjd))
									ref_timechan_success=False
									break
								elif other_error==True:				
									mainlog.info('Reference time frequency calibration failed because some error occured. Contact developer to fix the problem.\n')
									msfreqs_cal.remove(float(single_ref_freqstamp))
									break
						if c==1:
							mainlog.info('Number of reference ms calibration finished : '+str(len(all_touch_files))+'\n')
							c=0
						if len(all_touch_files)>old_finished:
							old_finished=len(all_touch_files)
							c=1
						else:
							c=0
						if len(all_touch_files)==len(ref_freq_time_msname):
							if len(all_success_list)>=2 and len(ref_freq_time_msname)!=1:						
								mainlog.info('Reference time frequency calibration is finished for '+str(len(all_success_list))+' reference frequency slices.\n') 
								ref_timechan_success=True
								ms_OBSIDs_gcal.append(single_reftimefreq_ms_OBSID)
								if inputs.do_bandpass==True:
									ms_OBSIDs_bcal.append(single_reftimefreq_ms_OBSID)
								if inputs.do_polcal==True:
									ms_OBSIDs_pcal.append(single_reftimefreq_ms_OBSID)
								break
							elif len(all_success_list)==0:
								mainlog.info('Reference time frequency calibration has been failed for all '+str(len(all_success_list))+' reference frequency slices.\n')
								mainlog.info('Trying with new reference ms.\n') 
								mstimes_cal.remove(float(ref_timestamp_mjd))
								ref_timechan_success=False
								break
							elif len(ref_freq_time_msname)==1 and len(all_success_list)==1:
								mainlog.info('One reference frequency ms was given. Reference time frequency calibration is finished for '\
												+str(len(all_success_list))+' reference frequency slice.\n') 
								ref_timechan_success=True
								ms_OBSIDs_gcal.append(single_reftimefreq_ms_OBSID)
								if inputs.do_bandpass==True:
									ms_OBSIDs_bcal.append(single_reftimefreq_ms_OBSID)
								if inputs.do_polcal==True:
									ms_OBSIDs_pcal.append(single_reftimefreq_ms_OBSID)
								break
						else:
							time.sleep(1.0)
							if updated_MWA_obsids==0:
								obsid_file,update_msg=update_mwa_obsids(verbose=inputs.verbose)
								if update_msg==0:
									updated_MWA_obsids=1
				elif calibrator_found==False and len(touch_file_list)==0 and len(ref_freq_time_msname)>0:
					for j in range(len(ref_freq_time_msname)):
						single_ref_freq_time_ms=ref_freq_time_msname[j]
						single_ref_freq_time_metafits=ref_time_freq_metafits[j]
						single_reftimefreq_ms_OBSID=ref_ms_OBSID_list[j]
						single_ref_freqstamp=ref_freqstamp[j]
						workdir=inputs.basedir+'/'+os.path.basename(single_ref_freq_time_ms).split('.ms')[0]
						if os.path.isdir(workdir)==False:
							if os.path.exists(workdir)==True:
								os.system('rm -rf '+workdir)
							os.makedirs(workdir)
						os.system('cp -r selfcal_inputs.py '+workdir+'/selfcal_inputs.py')
						if os.path.isfile(single_ref_freq_time_ms+'/copy.table.lock')==False:
							os.system('cp -r '+single_ref_freq_time_ms+'/table.lock '+single_ref_freq_time_ms+'/copy.table.lock')
						if os.path.isfile(single_ref_freq_time_ms+'/copy.table.info')==False:
							os.system('cp -r '+single_ref_freq_time_ms+'/table.info '+single_ref_freq_time_ms+'/copy.table.info')
						AM=AccessMS(single_ref_freq_time_ms)
						ref_time_avg=ref_time_bin=AM.calc_timeres()
						cpu_str=''
						for i in cpu_list:
							cpu_str+=str(i)+','
						cpu_str=cpu_str[:-1]
						cmd='run_paircars --msname '+single_ref_freq_time_ms+' --metafits '+single_ref_freq_time_metafits+' --basedir '+inputs.basedir+' --workdir '+workdir+' --ref_freq_avg '+str(40)+\
						' --ref_time_avg '+str(ref_time_avg)+' --ref_time_freq True --do_bandpass '+str(inputs.do_bandpass)+' --do_polcal '+str(inputs.do_polcal)+' --scratch True --wsclean '+str(use_wsclean)+\
						' --cpu_frac '+str(inputs.cpu_frac)+' --do_fluxcal '+str(do_fluxcal)+' --cubical '+str(inputs.use_cubical)+' --verbose '+str(inputs.verbose)+' --input_file '+input_file+\
						' --do_pbcor False --ref_time_bin '+str(ref_time_bin)+' --cpu_list '+cpu_str+' --pol_ref_time '+str(pol_ref_time)
						if inputs.ref_time!='':				
							cmd+=' --ref_time '+str(inputs.ref_time)
						else:
							if j==0:
								ds_base_name='Dynamic_Spectrum'
								dsname=inputs.basedir+'/'+ds_base_name
								if inputs.timerange=='':
									temp_timerange=''
								else:
									temp_timerange=inputs.timerange
								ds_base_name='Dynamic_Spectrum.npy'
								dsfile=workdir+'/'+ds_base_name
								savefile=workdir+'/ds_data'
								ms_mainlog.info('get_ref_time_fromDS(\''+single_ref_freq_time_ms+'\',\''+single_ref_freq_time_metafits+'\',nthreads=3,verbose='+str(inputs.verbose)+\
								',force_calculate=True,dsfile=\''+dsfile+'\',caltables=[],DRbased=False,T_limit=1,timerange=\''+str(temp_timerange)+\
								'\',docal=True,ref_ant='+str(inputs.ref_ant)+',cpu_list='+str(cpu_list)+',cpu_frac='+str(inputs.cpu_frac)+',save_file=\''+savefile+'\')\n')
								ref_time,ref_time_index,ref_timestamps,T_timeseries=get_ref_time_fromDS(single_ref_freq_time_ms,single_ref_freq_time_metafits,nthreads=3,verbose=inputs.verbose,force_calculate=True,\
								dsfile=dsfile,caltables=[],DRbased=False,T_limit=1,timerange=temp_timerange,docal=True,ref_ant=inputs.ref_ant,cpu_list=cpu_list,cpu_frac=inputs.cpu_frac,save_file=savefile)
								common_ref_time=ref_time
							else:
								ds_base_name='Dynamic_Spectrum.npy'
								dsfile=workdir+'/'+ds_base_name
								ms_mainlog.info('get_ref_time_fromDS(\''+single_ref_freq_time_ms+'\',\''+single_ref_freq_time_metafits+'\',nthreads=3,verbose='+str(inputs.verbose)+\
								',input_pol_ref_time=\''+str(ref_time)+'\',force_calculate=True,dsfile=\''+dsfile+'\',caltables=[],DRbased=False,T_limit=1,timerange=\''+str(temp_timerange)+\
								'\',docal=True,ref_ant='+str(inputs.ref_ant)+',cpu_list='+str(cpu_list)+',cpu_frac='+str(inputs.cpu_frac)+',save_file=\''+savefile+'\')\n')
								ref_time,ref_time_index,ref_timestamps,T_timeseries=get_ref_time_fromDS(single_ref_freq_time_ms,single_ref_freq_time_metafits,nthreads=3,verbose=inputs.verbose,\
									input_pol_ref_time=ref_time,force_calculate=True,dsfile=dsfile,caltables=[],DRbased=False,T_limit=1,timerange=temp_timerange,docal=True,ref_ant=inputs.ref_ant,cpu_list=cpu_list,\
									cpu_frac=inputs.cpu_frac,save_file=savefile)
							cmd+=' --ref_time '+str(ref_time)
							cmd+=' --common_ref_time '+str(common_ref_time)
						mainlog.info('Command : '+cmd+'\n')
						screen_name=str(single_reftimefreq_ms_OBSID)+'_'+str(os.path.basename(single_ref_freq_time_ms).split('.ms')[0])+'_runpaircars_'+str(inputs.job_id)
						finished_touch_file=inputs.basedir+'/.Finished_runpaircars_'+str(reftimefreq_ms_OBSID)+'_'+str(os.path.basename(single_ref_freq_time_ms).split('.ms')[0])
						screen_batch_file=paircars_instance_runner(cmd,inputs.basedir,inputs.paircars_dir,screen_name,finished_touch_file,inputs.job_id)	
						run_batch_scripts_cpulist(screen_batch_file,cpu_list=cpu_list)	
					mainlog.info('Waiting for finishing calibration of reference time frequency......\n')
					mainlog.info('Total reference ms : '+str(len(ref_freq_time_msname))+'\n')
					while True:
						all_touch_files=[]
						all_success_list=[]
						old_finished=0
						c=0
						for j in range(len(ref_freq_time_msname)):
							single_reftimefreq_ms_OBSID=ref_ms_OBSID_list[j]	
							single_ref_freq_time_ms=ref_freq_time_msname[j]
							single_ref_freqstamp=ref_freqstamp[j]
							touch_file_list=glob.glob(inputs.basedir+'/.ref_timechan_done_'+str(single_reftimefreq_ms_OBSID)+'_*')
							runtime_fail=False
							selfcal_fail=False
							other_error=False
							if len(touch_file_list)>0:
								for j in touch_file_list:
									if j not in all_touch_files:
										all_touch_files.append(j)
								i=touch_file_list[0]
								if '_selfcalerror' in i:
									selfcal_fail=True
								elif '_runtimeerror' in i:
									runtime_fail=True
								elif '_error' in i:
									other_error=True
								else:
									all_success_list.append(single_ref_freq_time_ms)
								if selfcal_fail==True:
									mainlog.info('Reference time frequency calibration failed. Trying with new reference time channel.\n')
									msfreqs_cal.remove(float(single_ref_freqstamp))
									break
								elif runtime_fail==True:
									mainlog.info('Reference time frequency calibration failed during run time. Some error occured during runtime. Contact developer to fix the problem.\n')
									msfreqs_cal.remove(float(single_ref_freqstamp))
									mstimes_cal.remove(float(ref_timestamp_mjd))
									ref_timechan_success=False
									break
								elif other_error==True:				
									mainlog.info('Reference time frequency calibration failed because some error occured. Contact developer to fix the problem.\n')
									msfreqs_cal.remove(float(single_ref_freqstamp))
									break
						if c==1:
							mainlog.info('Number of reference ms calibration finished : '+str(len(all_touch_files))+'\n')
							c=0
						if len(all_touch_files)>old_finished:
							old_finished=len(all_touch_files)
							c=1
						else:
							c=0
						if len(all_touch_files)==len(ref_freq_time_msname):
							if len(all_success_list)>=2 and len(ref_freq_time_msname)!=1:						
								mainlog.info('Reference time frequency calibration is finished for '+str(len(all_success_list))+' reference frequency slices.\n') 
								ref_timechan_success=True
								ms_OBSIDs_gcal.append(single_reftimefreq_ms_OBSID)
								if inputs.do_bandpass==True:
									ms_OBSIDs_bcal.append(single_reftimefreq_ms_OBSID)
								if inputs.do_polcal==True:
									ms_OBSIDs_pcal.append(single_reftimefreq_ms_OBSID)
								break
							elif len(all_success_list)==0:
								mainlog.info('Reference time frequency calibration has been failed for all '+str(len(all_success_list))+' reference frequency slices.\n')
								mainlog.info('Trying with new reference ms.\n') 
								mstimes_cal.remove(float(ref_timestamp_mjd))
								ref_timechan_success=False
								break
							elif len(ref_freq_time_msname)==1 and len(all_success_list)==1:
								mainlog.info('One reference frequency ms was given. Reference time frequency calibration is finished for '\
												+str(len(all_success_list))+' reference frequency slice.\n') 
								ref_timechan_success=True
								ms_OBSIDs_gcal.append(single_reftimefreq_ms_OBSID)
								if inputs.do_bandpass==True:
									ms_OBSIDs_bcal.append(single_reftimefreq_ms_OBSID)
								if inputs.do_polcal==True:
									ms_OBSIDs_pcal.append(single_reftimefreq_ms_OBSID)
								break
						else:
							time.sleep(1.0)
							if updated_MWA_obsids==0:
								obsid_file,update_msg=update_mwa_obsids(verbose=inputs.verbose)
								if update_msg==0:
									updated_MWA_obsids=1
				else: # Checked and debugged
					if len(ref_freq_time_msname)>0:
						mainlog.info('Reference time frequency calibration is done already.\n')
						for j in range(len(ref_freq_time_msname)):
							single_ref_freq_time_ms=ref_freq_time_msname[j]
							single_ref_freq_time_metafits=ref_time_freq_metafits[j]
							single_reftimefreq_ms_OBSID=ref_ms_OBSID_list[j]
							single_ref_freqstamp=ref_freqstamp[j]
							msfreqs_cal.remove(float(single_ref_freqstamp))
							mstimes_cal.remove(float(ref_timestamp_mjd))
							while True:
								if os.path.exists(inputs.basedir+'/Ref_time_freq_slice_output.npy'):
									result_dic=np.load(inputs.basedir+'/Ref_time_freq_slice_output.npy',allow_pickle=True).all()
									break
							for i in result_dic:
								if single_ref_freq_time_ms==i:
									result=result_dic[i]
									break
							ms=list(result.keys())[0]
							values=list(result.values())
							return_msg,ref_time,ref_chan,spawned_casa_instances,ref_freq_avg,ref_time_avg,pol_ref_time,ref_time_bin,common_ref_time=values 
							spawned_ms_jobs[os.path.basename(single_ref_freq_time_ms)]=[ref_time,int(ref_chan),int(spawned_casa_instances),float(ref_freq_avg),float(ref_time_avg),\
																						str(pol_ref_time),float(ref_time_bin)]
							spawned_ms.append(single_ref_freq_time_ms)
						ref_timechan_success=True
						break
					else:
						mainlog.info('No measurement set is available for calibration.\n')
						return 1

			per_ms_cal_jobs=0
			for j in range(len(ref_freq_time_msname)):
				index=measurement_set_list.index(ref_freq_time_msname[j])
				measurement_set_list.remove(ref_freq_time_msname[j]) # Removing ref time freq ms from ms list
				del metafits_dic[ref_freq_time_msname[j].split('.ms')[0].split('_timesliced')[0].split('_ref')[0]] # Removing ref time freq metafits from list 

			for i in measurement_set_list: # Removing non-fininshed ms directories
				a=glob.glob(inputs.basedir+'/.Finished_spawned_*_'+os.path.basename(i)+'*')
				if len(a)>0:
					measurement_set_list.remove(i)
				else:
					if os.path.isdir(inputs.basedir+'/'+os.path.basename(i).split('.ms')[0]):
						os.system('rm -rf '+inputs.basedir+'/'+os.path.basename(i).split('.ms')[0])
						os.system('rm -rf '+inputs.basedir+'/.Finished_spawned_*'+os.path.basename(i).split('.ms')[0]+'*')
						os.system('rm -rf '+inputs.basedir+'/.Finished_runpaircars_*'+os.path.basename(msname).split('.ms')[0]+'*')

			while True:
				time.sleep(10)
				if os.path.exists(inputs.basedir+'/Ref_time_freq_slice_output.npy'):
					try:
						result_dic=np.load(inputs.basedir+'/Ref_time_freq_slice_output.npy',allow_pickle=True).all()
						if len(result_dic)==len(all_success_list):
							break
					except Exception as e:
						mainlog.error('Error : '+str(e)+'\n')
						pass
			keys=spawned_ms_jobs.keys()
			ref_finished_ms_list=[]
			for i in result_dic:
				ms=i
				values=result_dic[i]
				return_msg,ref_time,ref_chan,spawned_instances,ref_freq_avg,ref_time_avg,pol_ref_time,ref_time_bin,common_ref_time=values 
				ref_finished_ms_list.append(ms)
				if os.path.basename(ms) not in keys:
					spawned_ms_jobs[os.path.basename(ms)]=[ref_time,int(ref_chan),int(spawned_instances),float(ref_freq_avg),float(ref_time_avg),str(pol_ref_time),float(ref_time_bin)]
					if ms not in spawned_ms:
						spawned_ms.append(ms)
					per_ms_cal_jobs=spawned_instances
			time.sleep(5)
			for msname in measurement_set_list_copy:
				if  os.path.exists('.Cal_done_'+os.path.basename(msname)+'_success')==True:
					if msname not in spawned_ms:
						spawned_ms.append(msname)
						measurement_set_list.remove(msname)
			
			mainlog.info('Number of calibration jobs per ms : '+str(per_ms_cal_jobs)+'\n')
			mainlog.info('Calibration jobs for all reference measurement sets have been started. Reference measurement sets are : '+str(ref_finished_ms_list)+'\n')
			##########################################
			#os.system('rm -rf '+inputs.basedir+'/Nonref_time_freq_slice_output.npy')
			touch_file_list=glob.glob(inputs.basedir+'/.timechan_done_*')+glob.glob(inputs.basedir+'/.ref_timechan_done_*')
			count=len(touch_file_list)
			if len(measurement_set_list)>0:
				mainlog.info('Remaining measurement sets : '+','.join(measurement_set_list)+'\n')
				while True:
					if os.path.exists(inputs.basedir+'/Ref_time_freq_slice_output.npy'):
						result_dic=np.load(inputs.basedir+'/Ref_time_freq_slice_output.npy',allow_pickle=True).all()
						result=list(result_dic.values())[0]
						break
				ref_time=result[1]
				ref_chan=result[2]
				ref_freq_avg=float(result[4])
				ref_time_avg=float(result[5])
				pol_ref_time=result[6]
				common_ref_time=result[8]
				ref_caltable_list=[]
				ref_caltable_time_list=[]
				cal_dirs=glob.glob(inputs.basedir+'/caltables/'+str(reftimefreq_ms_OBSID)+'/*')
				common_ref_time_mjd=timestamp_to_mjdsec(common_ref_time)
				for i in cal_dirs:
					caltables=glob.glob(i+'/*.cal')
					cal_times=np.array([timestamp_to_mjdsec(get_caltable_metadata(i)['Start time']) for i in caltables])
					if len(cal_times)!=0:
						nearpos=np.argmin(np.abs(cal_times-common_ref_time_mjd))
						ref_caltable_list.append(caltables[nearpos])
				caltable_str=','.join(ref_caltable_list)
				mainlog.info('Reference caltables : '+caltable_str+'\n')
				ref_freq_list=[]
				for j in range(len(ref_freq_time_msname)):
					single_reftimefreq_ms=ref_freq_time_msname[j]
					AMref=AccessMS(single_reftimefreq_ms)
					ref_freq_list.append(AMref.calc_meanfreq()/10**6)
				ref_freq_list=np.array(ref_freq_list)
				cent_ref_freq=ref_freq_list[int(len(ref_freq_list)/2)]
				num_ms_jobs=0
				prev_num_jobs=0
				for i in range(len(measurement_set_list)):
					if os.path.exists('.Cal_done_'+os.path.basename(msname)+'_success')==False:													
						msname=measurement_set_list[i]
						obsid=ms_OBSIDs_cal[i]
						AMms=AccessMS(msname)
						ms_freq=AMms.calc_meanfreq()/10**6
						freq_diff=np.min(abs(ms_freq-ref_freq_list))
						time_diff=abs(reftimefreq_ms_OBSID-obsid)
						if time_diff>bandpass_interval*3600:
							mainlog.info('Performing bandpass as time interval is greater than : '+str(bandpass_interval)+' hr.\n')
							do_bandpass=True
						else:
							if freq_diff<(skip_freq/10**3):
								mainlog.info('Do not perform bandpass, becuase time interval is smaller than : '+str(bandpass_interval)+'.\n')
								do_bandpass=False
							else:
								mainlog.info('Performing bandpass becuase frequency interval is greater than : '+str(skip_freq)+' kHz.\n')
								do_bandpass=True
						metafits=metafits_dic[msname.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]
						if inputs.do_polcal==True:
							ref_time_freq_gridpoint=fits.getheader(ref_time_freq_metafits[0])['GRIDNUM']
							ms_gridpoint=fits.getheader(metafits)['GRIDNUM']
							if ref_time_freq_gridpoint!=ms_gridpoint:
								do_polcal=True
								mainlog.info('Beam pointing changed. Performing polarisation calibration.\n')
							elif freq_diff>(skip_freq_pol/10**3):
								do_polcal=True
								mainlog.info('Frequency difference is more than polarisation frequency interval. Performing polarisation calibration.\n')
							elif time_diff>(bandpass_interval*2)*3600:
								do_polcal=True
								mainlog.info('Time difference is more than polarisation time interval. Performing polarisation calibration.\n')
							else:
								mainlog.info('Beam pointing is same. Do not perform polarisation calibration.\n')
								do_polcal=False
						if (do_bandpass==False and do_polcal==False and time_diff>0 and np.min(np.abs(ms_freq-ref_freq_list))<=skip_freq) or (do_bandpass==True or do_polcal==True):
							if do_bandpass==False and do_polcal==False and time_diff>0 and anp.min(np.abs(ms_freq-ref_freq_list))<=skip_freq:
								mainlog.info('Performing gaincal for other times.\n')
							workdir=inputs.basedir+'/'+os.path.basename(msname).split('.ms')[0]
							if os.path.isdir(workdir)==False:
								if os.path.isfile(workdir)==True:
									os.system('rm -rf '+workdir)
								os.makedirs(workdir)
							os.system('cp -r selfcal_inputs.py '+workdir+'/selfcal_inputs.py')
							if calibrator_found==True and (len(gaincal_list)!=0 or len(bandpass_list)!=0): # Not debugged
								calibrator_OBSID=np.min(np.array([x.split('_')[0] for x in gaincal_list])) # Considering naming convention OBSID_COARSECHAN.gcal
								ms_gaincal_list=[]
								ms_bandpass_list=[]
								for gcal in gaincal_list:
									if int(gcal.split('_')[0])==int(calibrator_OBSID):
										ms_gaincal_list.append(gcal)
								for bcal in bandpass_list:
									if int(bcal.split('_')[0])==int(calibrator_OBSID):
										ms_bandpass_list.append(bcal)
								bandpass_freq_list=[float(get_caltable_metadata(i)['Channel 0 frequency (MHz)']) for i in ms_bandpass_list]
								bandpass_coarse_chan=np.array([freq_to_MWA_coarse(freq) for freq in bandpass_freq_list])
								AMref=AccessMS(msname)
								ms_freq=float(AMref.calc_meanfreq()/10**6)
								ms_coarse_chan=freq_to_MWA_coarse(ms_freq)
								ms_bandpass_list=ms_bandpass_list[np.where(bandpass_coarse_chan==ms_coarse_chan)]
								caltable_list=ms_gaincal_list+ms_bandpass_list
							else:
								caltable_list=[]
							mainlog.info('Interpolating gaincal table to : '+inputs.basedir+'/caltables/'+str(reftimefreq_ms_OBSID)+'/Interp_'+str(ms_freq)+'MHz.gcalp.\n')
							interpolated_caltable=multifreq_gaincal_interpolate(gaintables=ref_caltable_list,outputfreq=ms_freq,\
											output_gaintable=inputs.basedir+'/caltables/'+str(reftimefreq_ms_OBSID)+'/Interp_'+str(ms_freq)+'MHz.gcalp')
							if len(caltable_list)!=0:
								interpolated_caltable=interpolated_caltable+caltable_list
							ms_obsid=get_OBSID_from_metafits(metafits)
							ms_OBSIDs_gcal.append(ms_obsid)
							if do_bandpass==True:
								ms_OBSIDs_bcal.append(ms_obsid)
							if do_polcal==True:
								ms_OBSIDs_pcal.append(ms_obsid)
							cpu_str=''
							for i in cpu_list:
								cpu_str+=str(i)+','
							cpu_str=cpu_str[:-1]
							cmd='run_paircars --msname '+msname+' --metafits '+metafits+' --basedir '+inputs.basedir+' --workdir '+workdir+' --ref_freq_avg '+str(ref_freq_avg)+' --ref_time_avg '+str(ref_time_avg)+\
							' --ref_time_freq False --do_bandpass '+str(do_bandpass)+' --do_polcal '+str(do_polcal)+' --cubical '+str(inputs.use_cubical)+' --scratch False --caltables '+interpolated_caltable+' --wsclean '+\
							str(use_wsclean)+' --cpu_frac '+str(inputs.cpu_frac)+' --do_fluxcal '+str(do_fluxcal)+' --verbose '+str(inputs.verbose)+' --input_file '+input_file+' --do_pbcor False --ref_time_bin '+\
							str(ref_time_bin)+' --cpu_list '+cpu_str+' --pol_ref_time '+str(pol_ref_time)
							if inputs.ref_time!='':				
								cmd+=' --ref_time '+str(inputs.ref_time)
							else:
								try:
									cmd+=' --ref_time '+str(ref_time)
								except:
									if inputs.timerange=='':
										temp_timerange=''
									else:
										temp_timerange=inputs.timerange
									ds_base_name='Dynamic_Spectrum.npy'
									dsfile=workdir+'/'+ds_base_name
									savefile=workdir+'/ds_data'
									ms_mainlog.info('get_ref_time_fromDS(\''+single_ref_freq_time_ms+'\',\''+single_ref_freq_time_metafits+'\',nthreads=3,verbose='+str(inputs.verbose)+\
									',force_calculate=True,dsfile=\''+dsfile+'\',caltables='+str(interpolated_caltable)+',DRbased=False,T_limit=1,timerange=\''+str(temp_timerange)+',input_pol_ref_time=\''+\
									str(ref_time)+'\',ref_ant='+str(inputs.ref_ant)+',cpu_list='+str(cpu_list)+',cpu_frac='+str(inputs.cpu_frac)+',save_file=\''+savefile+'\')\n')
									ref_time,ref_time_index,ref_timestamps,T_timeseries=get_ref_time_fromDS(single_ref_freq_time_ms,single_ref_freq_time_metafits,nthreads=3,verbose=inputs.verbose,\
									force_calculate=True,dsfile=dsfile,caltables=[interpolated_caltable],DRbased=False,T_limit=1,timerange=temp_timerange,input_pol_ref_time=ref_time,\
									ref_ant=inputs.ref_ant,cpu_list=cpu_list,cpu_frac=inputs.cpu_frac,save_file=savefile)
									cmd+=' --ref_time '+str(ref_time)
							num_ms_jobs=+per_ms_cal_jobs
							mainlog.info('Command : '+cmd+'\n')
							screen_name=str(ms_obsid)+'_'+str(os.path.basename(msname).split('.ms')[0])+'_runpaircars_'+str(inputs.job_id)
							finished_touch_file=inputs.basedir+'/.Finished_runpaircars_'+str(reftimefreq_ms_OBSID)+'_'+str(os.path.basename(msname).split('.ms')[0])
							screen_batch_file=paircars_instance_runner(cmd,inputs.basedir,inputs.paircars_dir,screen_name,finished_touch_file,inputs.job_id)
							run_batch_scripts_cpulist(screen_batch_file,cpu_list=cpu_list)	
							count-=1
							spawned_ms_jobs[os.path.basename(msname)]=[ref_time,int(ref_chan),int(per_ms_cal_jobs),float(ref_freq_avg),float(ref_time_avg),str(pol_ref_time),float(ref_time_bin)]
							spawned_ms.append(msname)
							basemsdir=os.path.basename(msname).split('.ms')[0]
							c+=1
							if c>len(measurement_set_list):
								break
							if count>0:
								continue
							else:
								mainlog.info('Waiting for previous calibrations to finish....\n')
								k=0		
								prev_int_k=0
								while True:
									present_touch_file_list=glob.glob(inputs.basedir+'/.timechan_done_*')+glob.glob(inputs.basedir+'/.ref_timechan_done_*')
									if len(present_touch_file_list)-len(touch_file_list)>0 and num_ms_jobs<max_allowed_jobs:
										count=len(present_touch_file_list)-len(touch_file_list)
										num_ms_jobs-=(count*per_ms_cal_jobs)
										touch_file_list=present_touch_file_list
										break
									elif num_ms_jobs<max_allowed_jobs:
										touch_file_list=present_touch_file_list
										count+=1
										break
									else:
										time.sleep(2)
										k+=1
										if (k//10)>prev_int_k:
											mainlog.info('Calibrations are running.....\n')
											prev_int_k=k//10
						else:
							mainlog.info('Passing the ms as it is not in same coarse channel of central reference frequency.\n')		

			mainlog.info('Calibration jobs for all measurement sets have been started. Waiting to spawn all calibration instances....\n')
			ms_list=sorted(spawned_ms)
			ms_OBSIDs=[]
			metafits_dic={}
			for ms in spawned_ms:
				index=ms_list_copy.index(ms)
				ms_OBSIDs.append(ms_OBSIDs_copy[index])
				metafits_dic[ms.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]=metafits_dic_copy[ms.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]
			time.sleep(1)
			# Spawning jobs for making final caltables and images
			#####################################################
			last_imaging_ms=''
			imaging_job_cpu_list={}
			tot_cal_jobs=int(len(cpu_list)/3)
			failed_mslist=[]
			fail_count=0
			while fail_count<5:
				if len(ms_list)==0:
					mainlog.info('No measurement set remained for calibration...\n')
					break
				for i in range(len(ms_list)):
					msname=ms_list[0]
					if msname[-1]=='/':
						msname=msname[:-1]
					metafits=metafits_dic[msname.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]
					OBSID=get_OBSID_from_metafits(metafits)
					if os.path.exists('.Cal_done_'+os.path.basename(msname)+'_success')==False:
						try:
							if inputs.only_calibration==False or inputs.database_imaging==True:
								AM=AccessMS(msname)
								total_bandwidth=AM.calc_bandwidth()/1000.0
								total_time=AM.calc_total_time()
								total_chunks=int(int(total_bandwidth/inputs.image_delta_freq)*int(total_time/inputs.image_delta_time))
								if total_chunks>=max_allowed_jobs:
									cpu_str=''
									count=0
									for cpu_num in cpu_list:
										if count>3:
											break
										cpu_str+=str(cpu_num)+','
										count+=1
									cpu_str=cpu_str[:-1]
									if i>tot_cal_jobs:
										mainlog.info('Waiting for finishing the final calibration of : '+str(os.path.basename(last_imaging_ms))+'\n')
										while True:
											finished_touch_file=glob.glob(str(inputs.basedir)+'/.Cal_done_'+str(os.path.basename(last_imaging_ms))+'*')
											if len(finished_touch_file)>0:
												break
											else:
												time.sleep(10)
								elif total_chunks<len(cpu_list):
									temp_cpu_list=cpu_list[0:total_chunks]
									imaging_job_cpu_list[os.path.basename(msname)]=temp_cpu_list
									cpu_list=cpu_list[total_chunks:]
									count=0
									for cpu_num in cpu_list:
										if count>3:
											break
										cpu_str+=str(cpu_num)+','
										count+=1
									cpu_str=cpu_str[:-1]
									if len(temp_cpu_list)==0:
										image_spawned_for_ms=list(imaging_job_cpu_list.keys())
										mainlog.info('Waiting for finishing the final calibration of : \n')
										for spawned_ms in image_spawned_for_ms:
											mainlog.info(os.path.basename(spawned_ms))
										while True:
											finished_touch_file=[]
											for i in image_spawned_for_ms:
												finished_touch_file.append(glob.glob(str(inputs.basedir)+'/.Cal_done_'+str(spawned_ms)+'*')[0])
											if len(finished_touch_file)>0:
												finished_ms=[a.split('.Cal_done_')[-1].split('.ms')[0]+'.ms' for a in finished_touch_file]
												for i in finished_ms:
													temp_finished_cpu_list=imaging_job_cpu_list[i]
													for cpu in temp_finished_cpu_list:
														cpu_list.append(cpu)
													del imaging_job_cpu_list[i]
												break
											else:
												time.sleep(10)
							keys=final_image_dic.keys()
							if OBSID not in keys:
								final_image_dic[OBSID]=[msname]
							else:
								x=final_image_dic[OBSID]
								x.append(msname)
								final_image_dic[OBSID]=x
							for i in glob.glob(inputs.basedir+'/imagemodels/*'):
								if i not in ms_OBSIDs_gcal:
									ms_OBSIDs_gcal.append(int(os.path.basename(i)))
							for i in glob.glob(inputs.basedir+'/bpimagemodels/*'):
								if i not in ms_OBSIDs_bcal:
									ms_OBSIDs_bcal.append(int(os.path.basename(i)))
							for i in glob.glob(inputs.basedir+'/polimagemodels/*'):
								if i not in ms_OBSIDs_pcal:
									ms_OBSIDs_pcal.append(int(os.path.basename(i)))
							gaincal_modeldirs_list=[inputs.basedir+'/imagemodels/'+str(i) for i in ms_OBSIDs_gcal]
							bandpass_modeldirs_list=[inputs.basedir+'/bpimagemodels/'+str(i) for i in ms_OBSIDs_bcal]
							polcal_modeldirs_list=[inputs.basedir+'/polimagemodels/'+str(i) for i in ms_OBSIDs_pcal]
							gaincal_ms_list=[inputs.basedir+'/gaincalms/'+str(i) for i in ms_OBSIDs_gcal]
							bandpass_ms_list=[inputs.basedir+'/bpcalms/'+str(i) for i in ms_OBSIDs_bcal]
							polcal_ms_list=[inputs.basedir+'/polms/'+str(i) for i in ms_OBSIDs_pcal]
							for i in gaincal_modeldirs_list:
								ms_image_dirs=glob.glob(i+'/*')
								for j in ms_image_dirs:
									images=glob.glob(j+'/*.image')
									for imagename in images:
										mainlog.info('paircars_image_bookkeep(\''+imagename+'\','+str(ref_freq_avg)+','+str(ref_time_bin)+',\''+bookkeep_file+'\')\n')
										paircars_image_bookkeep(imagename,ref_freq_avg,ref_time_bin,bookkeep_file)
							for i in polcal_modeldirs_list:
								ms_image_dirs=glob.glob(i+'/*')
								for j in ms_image_dirs:
									images=glob.glob(j+'/*.image')
									for imagename in images:
										mainlog.info('paircars_image_bookkeep(\''+imagename+'\','+str(ref_freq_avg)+','+str(ref_time_bin)+',\''+bookkeep_file+'\')\n')
										paircars_image_bookkeep(imagename,ref_freq_avg,ref_time_bin,bookkeep_file)
							temp_dirs=[]
							for i in gaincal_ms_list:
								if i not in temp_dirs:
									temp_dirs.append(i)
							gaincal_ms=','.join(temp_dirs)
							temp_dirs=[]
							for i in bandpass_ms_list:
								if i not in temp_dirs:
									temp_dirs.append(i)
							bandpass_ms=','.join(temp_dirs)
							temp_dirs=[]
							for i in polcal_ms_list:
								if i not in temp_dirs:
									temp_dirs.append(i)
							polcal_ms=','.join(temp_dirs)
							mainlog.info('Making final calibration tables for : '+msname+'\n')
							last_imaging_ms=msname
							num_jobs=spawned_ms_jobs[os.path.basename(msname)][2]
							spawned_file=glob.glob(inputs.basedir+'/.Finished_spawned_*'+str(OBSID)+'*_*'+os.path.basename(msname).split('.ms')[0]+'*')
							if len(spawned_file)>0:
								spawned_file=spawned_file[0]
								num_jobs=int(spawned_file.split(os.path.basename(msname).split('.ms')[0]+'_')[-1])
							freq_avg=spawned_ms_jobs[os.path.basename(msname)][3]
							time_avg=spawned_ms_jobs[os.path.basename(msname)][4]
							screen_name=str(OBSID)+'_'+str(os.path.basename(msname).split('.ms')[0])+'_final_caltable_'+str(inputs.job_id)
							batch_file=inputs.basedir+'/'+screen_name+'.batch'
							cmd_batch_file=inputs.basedir+'/'+screen_name+'_cmd.batch'
							cmd='manage_database --msname '+msname+' --metafits '+metafits+' --num_jobs '+str(num_jobs)+' --basedir '+inputs.basedir+\
							' --gaincal_msdir '+gaincal_ms+' --bandpass_msdir '+bandpass_ms+' --polcal_msdir '+polcal_ms+\
							' --localdatabase '+local_caldatabase+' --freqavg '+str(freq_avg)+' --timeavg '+str(time_avg)+' --inputfile '+inputs.basedir+\
							'/selfcal_inputs.py --verbose '+str(inputs.verbose)+' --database_imaging '+str(inputs.database_imaging)+\
							' --onlycal '+str(inputs.only_calibration)+' --cpu_list '+cpu_str
							finished_touch_file=inputs.basedir+'/.Finished_final_caltable_'+str(OBSID)+'_'+str(os.path.basename(msname).split('.ms')[0])
							screen_batch_file=paircars_instance_runner(cmd,inputs.basedir,inputs.paircars_dir,screen_name,finished_touch_file,inputs.job_id)
							present_cpu_list=[int(a) for a in cpu_str.split(',')]
							run_batch_scripts_cpulist(screen_batch_file,cpu_list=present_cpu_list)
							ms_list.remove(msname)
							time.sleep(0.5)
						except Exception as e:
							ms_list.remove(msname)
							mainlog.warning('Error in spawning calibration for ms : '+msname+'. It will be tried again.\n')
							mainlog.error('Error message : '+str(e)+'\n')
							failed_mslist.append(msname)
				ms_list=failed_mslist
				fail_count+=1			
		
			# Waiting for finishing calibrations and smoothing coarse channel edges
			#######################################################################
			t_count=0
			ms_list=sorted(spawned_ms)
			while True:
				bptables=sorted(glob.glob(str(local_caldatabase)+'/'+str(OBSID)+'/'+str(OBSID)+'_*.bcal'))
				if len(bptables)>=len(ms_list):
					if len(ms_list)==1:
						break
					mainlog.info('smooth_bandpass_edges(bptables='+str(bptables)+')\n')
					smooth_bandpass_edges(bptables=bptables)
					break
				else:
					time.sleep(10)
					if t_count==0:
						mainlog.info('Waiting for finishing calibrations...\n')
					t_count+=1
			caltables_available=True
			bptables=glob.glob(str(local_caldatabase)+'/'+str(OBSID)+'/'+str(OBSID)+'_*.bcal')
			gtables=glob.glob(str(local_caldatabase)+'/'+str(OBSID)+'/'+str(OBSID)+'_*.gcal')	
			if len(gtables)<len(ms_list):
				mainlog.error('Gain tables are not present for all measurement set.\n')
				caltables_available=False
			elif len(bptables)<len(ms_list) and caltables_available==True:
				mainlog.error('Bandpass tables are not present for all measurement set.\n')
				caltables_available=False
			if inputs.do_polcal==True and caltables_available==True:
				poltables=glob.glob(str(local_caldatabase)+'/'+str(OBSID)+'/'+str(OBSID)+'_*.pcal')		
				if len(poltables)<len(ms_list):
					mainlog.error('Polarization tables are not present for all measurement set.\n')
					caltables_available=False
		else: # If only imaging, check if all caltables are present or not
			caltables_available=True
			ms_list=measurement_set_list
			metafits=metafits_dic[msname.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]
			OBSID=get_OBSID_from_metafits(metafits)
			bptables=glob.glob(str(local_caldatabase)+'/'+str(OBSID)+'/'+str(OBSID)+'_*.bcal')
			gtables=glob.glob(str(local_caldatabase)+'/'+str(OBSID)+'/'+str(OBSID)+'_*.gcal')	
			if len(gtables)<len(ms_list):
				mainlog.error('Gain tables are not present for all measurement set.\n')
				caltables_available=False
			elif len(bptables)<len(ms_list) and caltables_available==True:
				mainlog.error('Bandpass tables are not present for all measurement set.\n')
				caltables_available=False
			if inputs.do_polcal==True and caltables_available==True:
				poltables=glob.glob(str(local_caldatabase)+'/'+str(OBSID)+'/'+str(OBSID)+'_*.pcal')		
				if len(poltables)<len(ms_list):
					mainlog.error('Polarization tables are not present for all measurement set.\n')
					caltables_available=False
			
		# Spawn final imaging
		####################
		cpu_list_copy=copy.deepcopy(cpu_list)
		if inputs.only_calibration==False and caltables_available==True: # Only if imaging is required
			count=0
			for i in range(len(ms_list)):
				msname=ms_list[i]
				metafits=metafits_dic[msname.split('.ms')[0].split('_timesliced')[0].split('_ref')[0]]
				OBSID=get_OBSID_from_metafits(metafits)
				keys=final_image_dic.keys()
				if OBSID not in keys:
					final_image_dic[OBSID]=[msname]
				else:
					x=final_image_dic[OBSID]
					x.append(msname)
					final_image_dic[OBSID]=x
				mainlog.info('Making final imaging for : '+msname+'\n')
				AM=AccessMS(msname)
				total_bandwidth=AM.calc_bandwidth()/1000.0
				total_time=AM.calc_total_time()
				total_chunks=int(int(total_bandwidth/inputs.image_delta_freq)*int(total_time/inputs.image_delta_time))
				tot_imaging_jobs=int(len(cpu_list)/total_chunks)
				if tot_imaging_jobs==0:
					tot_imaging_jobs=1
				if (total_chunks*len(ms_list))>=max_allowed_jobs:
					cpu_str=''
					count_cpu=0
					max_cpus=int(len(cpu_list)/tot_imaging_jobs)
					while len(cpu_list)!=0:
						if count_cpu==max_cpus:
							break
						cpu_str+=str(cpu_list[0])+','
						cpu_list.remove(cpu_list[0])
						count_cpu+=1
					cpu_str=cpu_str[:-1]
					if i>0 and count>=tot_imaging_jobs:
						mainlog.info('Waiting for finishing the final imaging of : '+str(os.path.basename(last_imaging_ms))+'\n')
						while True:
							finished_touch_file=glob.glob(str(inputs.basedir)+'/.Imaging_done_'+str(os.path.basename(last_imaging_ms))+'*')
							if len(finished_touch_file)>0:
								count-=len(finished_touch_file)
								if count<0:
									count=0
								break
							else:
								time.sleep(10)
				elif total_chunks<len(cpu_list):
					temp_cpu_list=cpu_list[0:total_chunks]
					imaging_job_cpu_list[os.path.basename(msname)]=temp_cpu_list
					cpu_list=cpu_list[total_chunks:]
					for cpu_num in cpu_list:
						cpu_str+=str(cpu_num)+','
					if len(temp_cpu_list)==0:
						image_spawned_for_ms=list(imaging_job_cpu_list.keys())
						mainlog.info('Waiting for finishing the imaging of : \n')
						for spawned_ms in image_spawned_for_ms:
							mainlog.info(os.path.basename(spawned_ms))
						while True:
							finished_touch_file=[]
							for i in image_spawned_for_ms:
								finished_touch_file.append(glob.glob(str(inputs.basedir)+'/.Imaging_done_'+str(spawned_ms)+'*')[0])
							if len(finished_touch_file)>0:
								finished_ms=[a.split('.Imaging_done_')[-1].split('.ms')[0]+'.ms' for a in finished_touch_file]
								for i in finished_ms:
									temp_finished_cpu_list=imaging_job_cpu_list[i]
									for cpu in temp_finished_cpu_list:
										cpu_list.append(cpu)
									del imaging_job_cpu_list[i]
								break
							else:
								time.sleep(10)

				# Final imaging
				#################
				AM=AccessMS(str(msname))
				tot_flag_frac=AM.calc_flagfrac()
				freqs=AM.get_freqs()/10**6
				coarse_chan_0=freq_to_MWA_coarse(freqs[0])
				coarse_chan_1=freq_to_MWA_coarse(freqs[-1])
				final_caltables=glob.glob(str(local_caldatabase)+'/'+str(OBSID)+'/'+str(OBSID)+'_'+str(coarse_chan_0)+'_'+str(coarse_chan_1)+'.*cal')
				caltable_string=','.join(final_caltables)
				mainlog.info('Final caltables : '+caltable_string+'\n')
				last_imaging_ms=msname
				screen_name=str(OBSID)+'_'+str(os.path.basename(msname).split('.ms')[0])+'_final_imaging_'+str(inputs.job_id)
				batch_file=inputs.basedir+'/'+screen_name+'.batch'
				cmd_batch_file=inputs.basedir+'/'+screen_name+'_cmd.batch'
				if cpu_str[-1]==',':
					cpu_str_final=cpu_str[:-1]
				else:
					cpu_str_final=cpu_str
				cmd='spawn_final_imaging --msname '+msname+' --metafits '+metafits+' --basedir '+inputs.basedir+' --localdatabase '+local_caldatabase+\
				' --inputfile '+inputs.basedir+'/selfcal_inputs.py --verbose '+str(inputs.verbose)+' --database_imaging '+str(inputs.database_imaging)+\
				' --cpu_list '+cpu_str_final+' --caltables '+caltable_string
				finished_touch_file=inputs.basedir+'/.Finished_final_imaging_'+str(OBSID)+'_'+str(os.path.basename(msname).split('.ms')[0])
				screen_batch_file=paircars_instance_runner(cmd,inputs.basedir,inputs.paircars_dir,screen_name,finished_touch_file,inputs.job_id)
				present_cpu_list=[int(a) for a in cpu_str_final.split(',')]
				run_batch_scripts_cpulist(screen_batch_file,cpu_list=present_cpu_list)
				time.sleep(0.5)
				count+=1
				if len(cpu_list)==0:
					cpu_list=copy.deepcopy(cpu_list_copy)

				if i==0 and inputs.send_notification:
					# Tracking final imaging status
					#############################
					keys=final_image_dic.keys()
					approx_image_num=int(len(ms_list)*total_chunks*(1-tot_flag_frac))
					for OBSID in keys:
						cmd='track_final_imaging --basedir '+str(inputs.basedir)+' --savedir '+str(inputs.final_image_dir)+' --OBSID '+str(OBSID)+' --email '+str(inputs.email)+' --localdatabase '\
							+local_caldatabase+' --database_imaging '+str(inputs.database_imaging)+' --only_cal '+str(inputs.only_calibration)+' --expected_image_num '+str(approx_image_num)
						screen_name=str(OBSID)+'_track_final_imaging'
						finished_touch_file=inputs.basedir+'/.Finished_trackimage_'+str(OBSID)
						screen_batch_file=paircars_instance_runner(cmd,inputs.basedir,inputs.paircars_dir,screen_name,finished_touch_file,inputs.job_id)
						screen_cmd='bash '+screen_batch_file
						mainlog.info(screen_cmd+'\n')
						os.system(screen_cmd)
						time.sleep(1.0)

		basemsdir=os.path.basename(msname).split('.ms')[0]
		if updated_MWA_obsids==0:
			obsid_file,update_msg=update_mwa_obsids(verbose=inputs.verbose)
			if update_msg==0:
				updated_MWA_obsids=1
		if updated_MWA_obsids==0:
			obsid_file,update_msg=update_mwa_obsids(verbose=inputs.verbose)
			if update_msg==0:
				updated_MWA_obsids=1	
		return 0	
	else:
		if updated_MWA_obsids==0:
			obsid_file,update_msg=update_mwa_obsids(verbose=inputs.verbose)
			if update_msg==0:
				updated_MWA_obsids=1
		if inputs.timerange!='':
			mainlog.info('No measurement set is present in the timerange : '+inputs.timerange+'\n')
		else:
			mainlog.info('No measurment set is found.\n')
		return 0
	
if __name__=='__main__':
	result=main()
	print('#############################\n')
	print('P-AIRCARS spawns all calibration and imaging jobs.\n')
	print('#############################\n')
	os._exit(result)

