#!/usr/bin/env python

class srcList:
	#arguments are:
	#sources (string, filename of LAT source list fits file in catalog format)
	#ft1 (string, filename of event file for which the xml will be used, only used to extract ROI info)
	#out (string, name of output xml file, defaults to mymodel.xml)
	def __init__(self,sources,ft1,out='mymodel.xml',DRversion=3):
		if not fileCheck(sources): #check that file exists
			print("Error:  %s not found."%sources)
			return
		if fileCheck(out):
			print('Warning: %s already exists, file will be overwritten if you proceed with makeModel.'%out)
		if DRversion not in [1,2,3]:
			print('Invalid choice, %i, of data release version.  Must be 1, 2, or 3.\nWill produce file assuming DR3')
			self.DR=3
			self.VarLim=24.725
		else:
			self.DR=DRversion
			self.VarLim=(24.725 if DRversion==3 else 21.666 if DRversion==2 else 18.475)
		self.srcs=sources
		self.out=out
		self.roi=getPos(ft1)
	
	#define a quick print function to make sure everything looks irght
	def Print(self):
		print('Source list file: ',self.srcs)
		print('Output file name: ',self.out)
		print('Selecting %s degrees around (ra,dec)=(%s,%s)'%(self.roi[2],self.ra,self.dec))
	
	#make the xml file
	#arguments are:
	#GDfile (str) -- optional, location and name of Galactic diffuse model to use
	#GDname (str) -- optional, name of Galactic diffuse component to use in xml model
	#ISOfile (str) -- optional, location and name of Isotropic diffuse template to use
	#ISOname (str) -- optional, name of Isotropic diffuse component to use in xml model
	#normsOnly (bool) -- optional, flag to only set normalizations parameters free
	#extDir (str) -- optional, directory with extended source templates
	#radLim (float) -- optional, radius in degrees from center of ROI beyond which source parameters are fixed
	#maxRad (float) -- optional, absolute maximum radius beyond which sources are fixed, this may be necessary when doing binned analysis and a variable source beyond radLim would be set free but this source is beyond the boundaries of the square region used for the binned likelihood
	#ExtraRad (float) -- optional, radius beyond ROI radius in event file out to which sources will be included with fixed parameters, defaul tof 10 is good for analyses starting around 100 MeV, but for higher energy fits this can be decreased
	#sigFree (float) -- optional, average significance, using FITS catalog file, below which source parameters are fixed, even if within radLim.  This corresponds to TS_value if using the XML catalog file.
	#varFree (float) -- optional, variability index above which source parameters are free, if beyond radLim and/or below sigFree only the normalization parameter is set free
	#psForce (bool) -- optional, flag to force extended sources to be point sources
	#E2CAT (bool) -- optional, flag to force use catalog names for extended sources (only matters if using catalog FITS file)
	#makeRegion (bool) -- optional, flag to also generate ds9 region file
	#GIndexFree (bool) -- optional, the Galactic diffuse is given a power-law spectral shape but the by default the index is frozen, setting this flag to True allows that to be free for additional freedom in diffuse fit
	#oldNames (bool) -- optional, flag to use the naming convention from make1FGLxml.py and make2FGLxml.py with a leading underscore and no spaces
	def makeModel(self,GDfile="$(FERMI_DIR)/refdata/fermi/galdiffuse/gll_iem_v07.fits",GDname='gll_iem_v07',ISOfile="$(FERMI_DIR)/refdata/fermi/galdiffuse/iso_P8R3_SOURCE_V3_v1.txt",ISOname='iso_P8R3_SOURCE_V3_v1',normsOnly=False,extDir='',radLim=-1,maxRad=None,ExtraRad=10,sigFree=5,varFree=True,psForce=False,E2CAT=False,makeRegion=True,GIndexFree=False,wd='',oldNames=False):
		self.radLim=(self.roi[2] if radLim<=0 else radLim)
		self.maxRad=(self.radLim if maxRad==None else maxRad)
		if self.maxRad<self.radLim:
			print("NOTE: maxRad (%.1f deg) is less than radLim (%.1f deg), meaning maxRad parameter is useless"%(self.maxRad,self.radLim))
		self.var=varFree
		self.psF=psForce
		self.E2C=E2CAT
		self.nO=normsOnly
		extDir=(extDir if extDir!='' else '$(FERMI_DIR)/data/pyBurstAnalysisGUI/templates')#make sure the default for FSSC STs is correct
		self.extD=(extDir if extDir[-1]=='/' else extDir+'/')
		self.ER=ExtraRad
		self.sig=sigFree
		self.reg=makeRegion
		self.GIF=GIndexFree
		if makeRegion:
			rhold=self.out.split('.')[:-1]
			rhold=rhold[0].split('/')[-1]#just in case the user has specified the full path for the output XML model file
			wd=(os.getcwd() if wd=='' else wd)
			self.regFile=wd+'/ROI_'
			for r in rhold:
				self.regFile+=r
			self.regFile+='.reg'
		print('Creating file and adding sources from 4FGL')
		#want ability to use either the FITS or xml versions of the catalog
		#need to tweak FITS version to have new functionality and then work out xml version
		if self.srcs.split('.')[-1]=='xml':
			addSrcsXML(self,GDfile,GDname,ISOfile,ISOname,oldNames)
		else:
			addSrcsFITS(self,GDfile,GDname,ISOfile,ISOname,oldNames)
	
import astropy.io.fits as pyfits
import os
from xml.dom import minidom
from xml.dom.minidom import parseString as pS
from numpy import floor,log10,cos,sin,arccos,pi,array,log,exp
acos=arccos
print("This is make4FGLDR3xml version 01r09.")
print("The default diffuse model files and names are for P8R3, 4FGL DRs, and the default fermitools directory locations.")
d2r=pi/180.
####varValue=24.725#21.666 for dr2,18.475 for 4fgl 

def addSrcsXML(sL,GD,GDn,ISO,ISOn,oldNames=False):
	varValue=sL.VarLim
	inputXml=minidom.parse(sL.srcs)
	outputXml=minidom.getDOMImplementation().createDocument(None,'source_library',None)
	outputXml.documentElement.setAttribute('title','source library')
	catalog=inputXml.getElementsByTagName('source')
	Sources={}
	ptSrcNum=0
	extSrcNum=0
	for src in catalog:
		if src.getAttribute('type')=='PointSource':
			for p in src.getElementsByTagName('spatialModel')[0].getElementsByTagName('parameter'):
				if p.getAttribute('name')=='RA':
					srcRA=float(p.getAttribute('value'))
				if p.getAttribute('name')=='DEC':
					srcDEC=float(p.getAttribute('value'))
		else:
			try:
			  srcDEC=float(src.getAttribute('DEC'))
			  srcRA=float(src.getAttribute('RA'))
			except:
			  for p in src.getElementsByTagName('spatialModel')[0].getElementsByTagName('parameter'):
			    if p.getAttribute('name')=='RA':
			      srcRA=float(p.getAttribute('value'))
			    if p.getAttribute('name')=='DEC':
			      srcDEC=float(p.getAttribute('value'))
		dist=angsep(sL.roi[0],sL.roi[1],srcRA,srcDEC) #check that source is within ROI radius + 10 degress of ROI center
		if srcRA==sL.roi[0] and srcDEC==sL.roi[1]:
			dist=0.0
		if dist<=sL.roi[2]+sL.ER:
			spec=src.getElementsByTagName('spectrum')
			specType=spec[0].getAttribute('type')
			specPars=spec[0].getElementsByTagName('parameter')
			Ext=(True if (src.getAttribute('type')=='DiffuseSource' and not sL.psF) else False)
			sname=src.getAttribute('name')
			
			if oldNames:#if you want the same naming convention as in make1FGLxml.py and make2FGLxml.py, e.g., preceeded by an underscore and no spaces
				sn='_'
				for N in str(sname).split(' '):
					sn+=N
			varIdx=float(src.getAttribute('Variability_Index'))
			#varidx=-1#hold over from old XMLs having and then not having and then now having Variability Index
			Sources[sname]={'ra':srcRA,'dec':srcDEC,'E':Ext,'stype':str(specType)}
			specOut=outputXml.createElement('spectrum')
			if str(specType)=='PLSuperExpCutoff2':
			  specOut.setAttribute('type','PLSuperExpCutoff')
			else:
			  specOut.setAttribute('type',specType)
			spatialOut=outputXml.createElement('spatialModel')
			srcOut=outputXml.createElement('source')
			srcOut.setAttribute('name',sname)
			srcOut.setAttribute('ROI_Center_Distance',"%.2f"%dist)
			if dist>=sL.roi[2] or dist>=sL.maxRad:
				Sources[sname]['free']=False
				for p in specPars:
				  specOut.appendChild(parameter_element("0","%s"%str(p.getAttribute('name')),"%s"%str(p.getAttribute('max')),"%s"%str(p.getAttribute('min')),"%s"%str(p.getAttribute('scale')),"%s"%str(p.getAttribute('value'))))
			elif dist>sL.radLim:
				if sL.var and varIdx>=varValue:
					Sources[sname]['free']=True
					for p in specPars:
					  freeFlag=("1" if p.getAttribute('name')==spec[0].getAttribute('normPar') else "0")
					  specOut.appendChild(parameter_element("%s"%freeFlag,"%s"%str(p.getAttribute('name')),"%s"%str(p.getAttribute('max')),"%s"%str(p.getAttribute('min')),"%s"%str(p.getAttribute('scale')),"%s"%str(p.getAttribute('value'))))
				else:
					Sources[sname]['free']=False
					for p in specPars:
					  specOut.appendChild(parameter_element("0","%s"%str(p.getAttribute('name')),"%s"%str(p.getAttribute('max')),"%s"%str(p.getAttribute('min')),"%s"%str(p.getAttribute('scale')),"%s"%str(p.getAttribute('value'))))
			elif float(src.getAttribute('TS_value'))>=sL.sig:
				Sources[sname]['free']=True
				#specOut.setAttribute('apply_edisp',ed)
				for p in specPars:
				  freeFlag=("1" if p.getAttribute('name')==spec[0].getAttribute('normPar') or (not sL.nO and p.getAttribute('free')=="1") else "0")
				  specOut.appendChild(parameter_element("%s"%freeFlag,"%s"%str(p.getAttribute('name')),"%s"%str(p.getAttribute('max')),"%s"%str(p.getAttribute('min')),"%s"%str(p.getAttribute('scale')),"%s"%str(p.getAttribute('value'))))
			else:
				if sL.var and varIdx>=varValue:
					Sources[sname]['free']=True
					#specOut.setAttribute('apply_edisp',ed)
					for p in specPars:
					  freeFlag=("1" if p.getAttribute('name')==spec[0].getAttribute('normPar') else "0")
					  specOut.appendChild(parameter_element("%s"%freeFlag,"%s"%str(p.getAttribute('name')),"%s"%str(p.getAttribute('max')),"%s"%str(p.getAttribute('min')),"%s"%str(p.getAttribute('scale')),"%s"%str(p.getAttribute('value'))))
				else:
					Sources[sname]['free']=False
					#specOut.setAttribute('apply_edisp','false')
					for p in specPars:
					  specOut.appendChild(parameter_element("0","%s"%str(p.getAttribute('name')),"%s"%str(p.getAttribute('max')),"%s"%str(p.getAttribute('min')),"%s"%str(p.getAttribute('scale')),"%s"%str(p.getAttribute('value'))))
			if Ext:
				spatial=src.getElementsByTagName('spatialModel')
				spatType=spatial[0].getAttribute('type')
				spatPars=spatial[0].getElementsByTagName('parameter')
				if str(spatType)=='SpatialMap':
				  spatialOut.setAttribute('type','SpatialMap')
				  spatialOut.setAttribute('map_based_integral','true')
				  efile=sL.extD+spatial[0].getAttribute('file').split('/')[-1]
				  spatialOut.setAttribute('file',efile)
				  print('Extended source %s in ROI, make sure %s is the correct path to the extended template.'%(sname,efile))
				else:#have to do above to get correct extended source template file localtion
				  spatialOut.setAttribute('type',str(spatType))
				  for p in spatPars:#for radial disks and gaussians, can just do the following
				    spatialOut.appendChild(parameter_element("0","%s"%str(p.getAttribute('name')),"%s"%str(p.getAttribute('max')),"%s"%str(p.getAttribute('min')),"%s"%str(p.getAttribute('scale')),"%s"%str(p.getAttribute('value'))))
				    print('Extended source %s in ROI, with %s spatial model.'%(sname,str(spatType)))
				srcOut.setAttribute('type','DiffuseSource')
				extSrcNum+=1
				#print('Extended source %s in ROI, make sure %s is the correct path to the extended template.'%(sname,efile))
			else:
				spatialOut.setAttribute('type','SkyDirFunction')
				spatialOut.appendChild(parameter_element("0","RA","360.0","-360.0","1.0","%.4f"%srcRA))
				spatialOut.appendChild(parameter_element("0","DEC","360.0","-360.0","1.0","%.4f"%srcDEC))
				srcOut.setAttribute('type','PointSource')
				ptSrcNum+=1
			srcOut.appendChild(specOut)
			srcOut.appendChild(spatialOut)
			outputXml.documentElement.appendChild(srcOut)
	gal=outputXml.createElement('source')
	gal.setAttribute('name',GDn)
	gal.setAttribute('type','DiffuseSource')
	galspec=outputXml.createElement('spectrum')
	galspec.setAttribute('type','PowerLaw')
	#galspec.setAttribute('apply_edisp','false')
	galspec.appendChild(parameter_element("1","Prefactor","10","0","1","1"))
	if sL.GIF:
		galspec.appendChild(parameter_element("1","Index","1","-1","1","0"))
	else:
		galspec.appendChild(parameter_element("0","Index","1","-1","1","0"))
	galspec.appendChild(parameter_element("0","Scale","1e6","2e1","1","100"))
	galspatial=outputXml.createElement('spatialModel')
	galspatial.setAttribute('type','MapCubeFunction')
	galspatial.setAttribute('file',GD)
	galspatial.appendChild(parameter_element("0","Normalization","1e3","1e-3","1","1"))
	gal.appendChild(galspec)
	gal.appendChild(galspatial)
	outputXml.documentElement.appendChild(gal)
	iso=outputXml.createElement('source')
	iso.setAttribute('name',ISOn)
	iso.setAttribute('type','DiffuseSource')
	isospec=outputXml.createElement('spectrum')
	isospec.setAttribute('type','FileFunction')
	isospec.setAttribute('file',ISO)
	isospec.setAttribute('apply_edisp','false')
	isospec.appendChild(parameter_element("1","Normalization","10","0.01","1","1"))
	isospatial=outputXml.createElement('spatialModel')
	isospatial.setAttribute('type','ConstantValue')
	isospatial.appendChild(parameter_element("0","Value","10","0","1","1"))
	iso.appendChild(isospec)
	iso.appendChild(isospatial)
	outputXml.documentElement.appendChild(iso)
	xmlStr=outputXml.toprettyxml(' ').splitlines(True)
	outStr=filter(lambda xmlStr: len(xmlStr) and not xmlStr.isspace(),xmlStr)
	outfile=open(sL.out,'w')
	outfile.write(''.join(outStr))
	outfile.close()
	if not sL.psF:
		print('Added %i point sources and %i extended sources'%(ptSrcNum,extSrcNum))
		if extSrcNum>0:
			print('If using unbinned likelihood you will need to rerun gtdiffrsp for the extended sources or rerun the makeModel function with optional argument psForce=True')
	else:
		print('Added %i point sources, note that any extended sources in ROI were modeled as point sources becaue psForce option was set to True'%ptSrcNum)
	if sL.reg:
		BuildRegion(sL,Sources)
	return

#function to cycle through the source list and add point source entries
def addSrcsFITS(sL,GD,GDn,ISO,ISOn,oldNames):
	model=open(sL.out,'w') #open file in write mode, overwrites other files of same name
	catfile=pyfits.open(sL.srcs) #open source list file and access necessary fields, requires LAT source catalog definitions and names
	data=catfile['LAT_Point_Source_Catalog'].data
	extendedinfo=catfile['ExtendedSources'].data
	extName=extendedinfo.field('Source_Name')	
	extFile=extendedinfo.field('Spatial_Filename')
	extFunc=extendedinfo.field('Spatial_Function')
	extSize=extendedinfo.field('Model_SemiMajor')
	extRa=extendedinfo.field('RAJ2000')
	extDec=extendedinfo.field('DEJ2000')
	name=data.field('Source_Name')
	Sigvals=data.field('Signif_Avg')
	VarIdx=data.field('Variability_Index')
	EName=data.field('Extended_Source_Name')
	ra=data.field('RAJ2000')
	dec=data.field('DEJ2000')
	plflux=data.field('PL_Flux_Density')
	lpflux=data.field('LP_Flux_Density')
	coflux=data.field('PLEC_Flux_Density')
	pivot=data.field('Pivot_Energy')
	plIndex=data.field('PL_Index')
	lpIndex=data.field('LP_Index')
	lpbeta=data.field('LP_beta')
	if sL.DR==3:
		plecIndex=data.field('PLEC_IndexS')
		plecexpFact=data.field('PLEC_ExpfactorS')
		plecexpIndex=data.field('PLEC_Exp_Index')
	else:
		plecIndex=data.field('PLEC_Index')
		plecexpFact=data.field('PLEC_Expfactor')
		plecexpIndex=data.field('PLEC_Exp_Index')
	spectype=data.field('SpectrumType')
	model.write('<?xml version="1.0" ?>\n')
	model.write('<source_library title="source library">\n')
	model.write('\n<!-- Point Sources -->\n')
	step=(sL.roi[2]+sL.ER)/5. #divide ROI radius plus ExtraRadius degrees into 5 steps for ordering of sources
	i=1
	radii=[]
	ptSrcNum=0
	extSrcNum=0
	Sources={}#dictionary for sources, useful for creating region file later.
	while i<6:
		if i*step<=sL.roi[2]+sL.ER:
			radii+=[step*i]
		else:
			radii+=[sL.roi[2]+sL.ER] #just in case of rounding errors
		i+=1
	for x in radii:
		if x==sL.roi[2]+sL.ER:
			model.write('\n<!-- Sources between [%s,%s] degrees of ROI center -->\n' %(x-step,x))
		else:
			model.write('\n<!-- Sources between [%s,%s) degrees of ROI center -->\n' %(x-step,x))
		for n,plf,lpf,cof,r,d,p,pli,lpi,lpb,pleci,plecef,plecei,t,TS,En,vi in zip(name,plflux,lpflux,coflux,ra,dec,pivot,plIndex,lpIndex,lpbeta,plecIndex,plecexpFact,plecexpIndex,spectype,Sigvals,EName,VarIdx):
			E=(True if n[-1]=='e' else False)
			dist=angsep(sL.roi[0],sL.roi[1],r,d) #check that source is within ROI radius + 10 degress of ROI center
			if r==sL.roi[0] and d==sL.roi[1]:
				dist=0.0
			if (dist<x and dist>=x-step) or (x==sL.roi[2]+10. and dist==x):
				if E and not sL.psF:
					Sources[En]={'ra':r,'dec':d,'stype':t,'E':E}
					extSrcNum+=1
					Name='<source ROI_Center_Distance="%.3f" name="%s" type="DiffuseSource">\n' %(dist,En)
				else:
					if E and not sL.E2C:#even if forcing all to point sources, use extended name except if E2CAT flag is set
						Sources[En]={'ra':r,'dec':d,'stype':t,'E':E}
						Name='<source ROI_Center_Distance="%.3f" name="%s" type="PointSource">\n' %(dist,En)
					else:
						Sources[n]={'ra':r,'dec':d,'stype':t,'E':E}
						if oldNames:
							srcname='_'
							for N in n.split(' '):
								srcname+=N
							Name='<source ROI_Center_Distance="%.3f" name="%s" type="PointSource">\n' %(dist,srcname)
						else:
							Name='<source ROI_Center_Distance="%.3f" name="%s" type="PointSource">\n' %(dist,n)
					ptSrcNum+=1
				if t=='PowerLaw':
					#uncomment out the two lines immediately following later
					spec,free=PLspec(sL,plf,pli,p,dist,TS,vi)
				elif t=='LogParabola':
					fixAll=(True if n=='4FGL J0534.5+2201i' else False)
					spec,free=LPspec(sL,lpf,lpi,p,lpb,dist,TS,vi,fixAll)
				else:
					if sL.DR==3:
						spec,free=CO4spec(sL,cof,pleci,p,plecef,plecei,dist,TS,vi)
					else:
						spec,free=CO2spec(sL,cof,pleci,p,plecef,plecei,dist,TS,vi)
				if E and not sL.E2C:
					Sources[En]['free']=free
				else:
					Sources[n]['free']=free
				if E and not sL.psF:
					efile=None
					efunc=None
					eSize=None
					eR=None
					eD=None
					for EXTNAME,EXTFILE,EXTFUNC,EXTSIZE,EXTRA,EXTDEC in zip(extName,extFile,extFunc,extSize,extRa,extDec):
						if En==EXTNAME:
							efunc=EXTFUNC
							efunc=('RadialGaussian' if efunc=='RadialGauss' else efunc)
							if efunc=='SpatialMap':
							  efile=sL.extD+EXTFILE
							else:
							  eSize=(EXTSIZE/(-2*log(0.32))**0.5 if efunc=='RadialGaussian' else EXTSIZE)
							  eR=EXTRA
							  eD=EXTDEC
					if efunc=='SpatialMap':
					  if efile==None:
					    print('could not find a match for',En,'in the list:')
					    print(extName)
					    efile=''
					  skydir='\t<spatialModel file="%s" map_based_integral="true" type="SpatialMap">\n'%(efile)
					  print('Extended source %s in ROI, make sure %s is the correct path to the extended template.'%(En,efile))
					  skydir+='\t\t<parameter free="0" max="1000" min="0.001" name="Prefactor" scale="1" value="1"/>\n'
					  skydir+='\t</spatialModel>\n'
					else:
					  skydir='\t<spatialModel type="%s">\n'%efunc
					  skydir+='\t<parameter free="0" max="360" min="-360" name="RA" scale="1" value="%s"/>\n'%eR
					  skydir+='\t<parameter free="0" max="90" min="-90" name="DEC" scale="1" value="%s"/>\n'%eD
					  if efunc=='RadialDisk':
					    skydir+='\t<parameter free="0" max="10" min="0" name="Radius" scale="1" value="%s"/>\n'%eSize
					  else:
					    skydir+='\t<parameter free="0" max="10" min="0" name="Sigma" scale="1" value="%s"/>\n'%eSize
					  skydir+='\t</spatialModel>\n'
					  print('Extended source %s in ROI with %s spatial model.'%(En,efunc))
				else:
					skydir='\t<spatialModel type="SkyDirFunction">\n'
					skydir+='\t\t<parameter free="0" max="360.0" min="-360.0" name="RA" scale="1.0" value="%s"/>\n' %r
					skydir+='\t\t<parameter free="0" max="90" min="-90" name="DEC" scale="1.0" value="%s"/>\n' %d
					skydir+='\t</spatialModel>\n'
				skydir+='</source>'
				(src,)=(Name+spec+skydir,)
				ptsrc=pS(src).getElementsByTagName('source')[0]
				ptsrc.writexml(model)
				model.write('\n')
	catfile.close() #close file
	if not sL.psF:
		print('Added %i point sources and %i extended sources'%(ptSrcNum,extSrcNum))
		if extSrcNum>0:
			print('If using unbinned likelihood you will need to rerun gtdiffrsp for the extended sources or rerun the makeModel function with optional argument psForce=True')
	else:
		print('Added %i point sources, note that any extended sources in ROI were modeled as point sources becaue psForce option was set to True'%ptSrcNum)
	#add galactic diffuse with PL spectrum, fix index to zero for general use, those who want it to be free can unfreeze parameter manually
	model.write('\n<!-- Diffuse Sources -->\n')
	Name='\n<source name="%s" type="DiffuseSource">\n' %GDn
	spec='\t<spectrum type="PowerLaw">\n'
	spec+='\t\t<parameter free="1" max="10" min="0" name="Prefactor" scale="1" value="1"/>\n'
	if sL.GIF:
		spec+='\t\t<parameter free="1" max="1" min="-1" name="Index" scale="1.0" value="0"/>\n'
	else:
		spec+='\t\t<parameter free="0" max="1" min="-1" name="Index" scale="1.0" value="0"/>\n'
	spec+='\t\t<parameter free="0" max="2e2" min="5e1" name="Scale" scale="1.0" value="1e2"/>\n'
	spec+='\t</spectrum>\n'
	skydir='\t<spatialModel file="%s" type="MapCubeFunction">\n' %GD
	skydir+='\t\t<parameter free="0" max="1e3" min="1e-3" name="Normalization" scale="1.0" value="1.0"/>\n'
	skydir+='\t</spatialModel>\n'
	skydir+='</source>'
	(src,)=(Name+spec+skydir,)
	galdiff=pS(src).getElementsByTagName('source')[0]
	galdiff.writexml(model)
	model.write('\n')
	Name='<source name="%s" type="DiffuseSource">\n' %ISOn
	spec='\t<spectrum type="FileFunction" file="%s"  apply_edisp="false">\n' %ISO
	spec+='\t\t<parameter free="1" max="10" min="1e-2" name="Normalization" scale="1" value="1"/>\n'
	spec+='\t</spectrum>\n'
	skydir='\t<spatialModel type="ConstantValue">\n'
	skydir+='\t\t<parameter free="0" max="10.0" min="0.0" name="Value" scale="1.0" value="1.0"/>\n'
	skydir+='\t</spatialModel>\n'
	skydir+='</source>'
	(src,)=(Name+spec+skydir,)
	iso=pS(src).getElementsByTagName('source')[0]
	iso.writexml(model)	
	model.write('\n</source_library>')
	model.close()
	if sL.reg:
		BuildRegion(sL,Sources)
	return

def BuildRegion(sL,Sources):
	myreg=open(sL.regFile,'w')#note that this will overwrite previous region files of the same name
	myreg.write('# Region File format: DS9 version ?')#I don't actually know, but I think it's one of the later ones, need to verify
	myreg.write('\n# Created by make3FGLxml.py')
	myreg.write('\nglobal font="roman 10 normal" move =0')
	for k in Sources.keys():
		src=Sources[k]
		#get color based on if the source is free or not
		color=('green' if src['free'] else 'magenta')
		if src['E']:#if the source is extended, always have the point be a "big" box
			myreg.write('\nJ2000;point(%.3f,%.3f) # point = box 18 color = %s text={%s}'%(src['ra'],src['dec'],color,k))
		else:#if the source is a point source, choose the point type based on spectral model
			ptype=('cross' if (src['stype']=='PLSuperExpCutoff' or src['stype']=='PLSuperExpCutoff2') else 'diamond' if src['stype']=='LogParabola' else 'circle')
			myreg.write('\nJ2000;point(%.3f,%.3f) # point = %s 15 color = %s text={%s}'%(src['ra'],src['dec'],ptype,color,k))
	myreg.close()
	return

def PLspec(sL,f,i,p,dist,TS,vi):
	varValue=sL.VarLim
	fscale=int(floor(log10(f)))
	spec='\t<spectrum type="PowerLaw">\n'
	spec+='\t<!-- Source is %s degrees away from ROI center -->\n' %dist
	if dist>sL.roi[2]: #if beyond ROI, shouldn't attempt to fit parameters
		spec+='\t<!-- Source is outside ROI, all parameters should remain fixed -->\n'
		spec+='\t\t<parameter free="0" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
		free=False
	elif(dist>sL.radLim):
		if dist>sL.maxRad:
			spec+='\t<!-- Source is outside specified radius limit of %s -->\n'%sL.radLim
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		elif vi<varValue or not sL.var:
			spec+='\t<!-- Source is outside specified radius limit of %s -->\n'%sL.radLim
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		else:
			spec+='\t<!-- Source is outside specified radius limit of %s but variability index %.2f is greater than %.2f and varFree set to True-->\n'%(sL.radLim,vi,varValue)
			free=True
			spec+='\t\t<parameter free="1" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
	elif(TS<sL.sig):
		if vi<varValue or not sL.var:
			spec+='\t<!-- Source signficance %.1f is less than specified minimum for a free source of %s -->\n'%(TS,sL.sig)
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		else:
			spec+='\t<!-- Source significane %.1f is less than specified minimum for a free source of %s but variability index %.2f is greater than %.2f and varFree is set to True -->\n'%(TS,sL.sig,vi,varValue)
			spec+='\t\t<parameter free="1" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=True
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
	else:
		spec+='\t\t<parameter free="1" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		free=True
		if sL.nO:
			spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
		else:
			spec+='\t\t<parameter free="1" max="10.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
	spec+='\t\t<parameter free="0" max="5e5" min="30" name="Scale" scale="1.0" value="%f"/>\n' %p
	spec+='\t</spectrum>\n'
	return spec,free

def CO2spec(sL,f,i,p,a,ei,dist,TS,vi):
	varValue=sL.VarLim
	f*=exp(a*(p**ei))
	fscale=int(floor(log10(f)))
	spec='\t<spectrum type="PLSuperExpCutoff2">\n'
	spec+='\t<!-- Source is %s degrees away from ROI center -->\n' %dist
	#i=(i if i>=0 else 2.)#some pulsars with index1 < 0 assuming standard convention, means rising counts spectrum at low E, odd
	if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
		spec+='\t<!-- Source is outside ROI, all parameters should remain fixed -->\n'
		free=False
		spec+='\t\t<parameter free="0" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="Index1" scale="-1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="100" min="-1" name="Expfactor" scale="0.001" value="%f"/>\n'%(a*1000.)
	elif(dist>sL.radLim):
		if dist>sL.maxRad:
			spec+='\t<!-- Source is outside specified radius limit of %s -->\n'%sL.radLim
			spec+='\t\t<parameter free="0" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		elif vi<varValue or not sL.var:
			spec+='\t<!-- Source is outside specified radius limit of %s -->\n'%sL.radLim
			spec+='\t\t<parameter free="0" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		else:
			spec+='\t<!-- Source is outside specified radius limit of %s but variability index %.2f is greater than %.2f and varFree is set to True -->\n'%(sL.radLim,vi,varValue)
			spec+='\t\t<parameter free="1" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=True
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="Index1" scale="-1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="100" min="-1" name="Expfactor" scale="0.001" value="%f"/>\n'%(a*1000.)
	elif(TS<sL.sig):
		if vi<varValue or not sL.var:
			spec+='\t<!-- Source significance %.1f is less than specified minimum for a free source of %s -->\n'%(TS,sL.sig)
			spec+='\t\t<parameter free="0" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		else:
			spec+='\t<!-- Source significance %.1f is less than specified minimum for a free source of %s but variability index %.2f is greater than %.2f and varFree is set to True -->\n'%(TS,sL.sig,vi,varValue)
			spec+='\t\t<parameter free="1" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=True
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="Index1" scale="-1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="100" min="-1" name="Expfactor" scale="0.001" value="%f"/>\n'%(a*1000.)
	else:
		free=True
		spec+='\t\t<parameter free="1" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		if sL.nO:
			spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="Index1" scale="-1.0" value="%s"/>\n' %i
			spec+='\t\t<parameter free="0" max="100" min="-1" name="Expfactor" scale="0.001" value="%f"/>\n'%(a*1000.)
		else:
			spec+='\t\t<parameter free="1" max="10.0" min="0.0" name="Index1" scale="-1.0" value="%s"/>\n' %i
			spec+='\t\t<parameter free="1" max="100" min="-1" name="Expfactor" scale="0.001" value="%f"/>\n'%(a*1000.)
	spec+='\t\t<parameter free="0" max="5e5" min="30" name="Scale" scale="1.0" value="%f"/>\n'%p
	spec+='\t\t<parameter free="0" max="5" min="0" name="Index2" scale="1.0" value="%f"/>\n'%ei
	spec+='\t</spectrum>\n'
	return spec,free

def CO4spec(sL,f,i,p,a,ei,dist,TS,vi):
	varValue=sL.VarLim
	fscale=int(floor(log10(f)))
	spec='\t<spectrum type="PLSuperExpCutoff4">\n'
	spec+='\t<!-- Source is %s degrees away from ROI center -->\n' %dist
	if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
		spec+='\t<!-- Source is outside ROI, all parameters should remain fixed -->\n'
		free=False
		spec+='\t\t<parameter free="0" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="IndexS" scale="-1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="100" min="-10" name="ExpfactorS" scale="0.1" value="%f"/>\n'%(a*10.)
	elif(dist>sL.radLim):
		if dist>sL.maxRad:
			spec+='\t<!-- Source is outside specified radius limit of %s -->\n'%sL.radLim
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		elif vi<varValue or not sL.var:
			spec+='\t<!-- Source is outside specified radius limit of %s -->\n'%sL.radLim
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		else:
			spec+='\t<!-- Source is outside specified radius limit of %s but variability index %.2f is greater than %.2f and varFree is set to True -->\n'%(sL.radLim,vi,varValue)
			spec+='\t\t<parameter free="1" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=True
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="IndexS" scale="-1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="100" min="-10" name="ExpfactorS" scale="0.1" value="%f"/>\n'%(a*10.)
	elif(TS<sL.sig):
		if vi<varValue or not sL.var:
			spec+='\t<!-- Source significance %.1f is less than specified minimum for a free source of %s -->\n'%(TS,sL.sig)
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		else:
			spec+='\t<!-- Source significance %.1f is less than specified minimum for a free source of %s but variability index %.2f is greater than %.2f and varFree is set to True -->\n'%(TS,sL.sig,vi,varValue)
			spec+='\t\t<parameter free="1" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=True
		spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="IndexS" scale="-1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="100" min="-10" name="ExpfactorS" scale="0.1" value="%f"/>\n'%(a*10.)
	else:
		free=True
		spec+='\t\t<parameter free="1" max="1e6" min="0" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		if sL.nO:
			spec+='\t\t<parameter free="0" max="10.0" min="0.0" name="IndexS" scale="-1.0" value="%s"/>\n' %i
			spec+='\t\t<parameter free="0" max="100" min="-10" name="ExpfactorS" scale="0.1" value="%f"/>\n'%(a*10.)
		else:
			spec+='\t\t<parameter free="1" max="10.0" min="0.0" name="IndexS" scale="-1.0" value="%s"/>\n' %i
			spec+='\t\t<parameter free="1" max="100" min="-10" name="ExpfactorS" scale="0.1" value="%f"/>\n'%(a*10.)
	spec+='\t\t<parameter free="0" max="5e5" min="30" name="Scale" scale="1.0" value="%f"/>\n'%p
	spec+='\t\t<parameter free="0" max="5" min="0" name="Index2" scale="1.0" value="%f"/>\n'%ei
	spec+='\t</spectrum>\n'
	return spec,free

def LPspec(sL,f,i,p,b,dist,TS,vi,fixAll=False):
	varValue=sL.VarLim
	fscale=int(floor(log10(f)))
	spec='\t<spectrum type="LogParabola">\n'
	spec+='\t<!-- Source is %s degrees away from ROI center -->\n' %dist
	if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
		spec+='\t<!-- Source is outside ROI, all parameters should remain fixed -->\n'
		spec+='\t\t<parameter free="0" max="1e6" min="0" name="norm" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		spec+='\t\t<parameter free="0" max="5.0" min="0.0" name="alpha" scale="1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="10.0" min="-5.0" name="beta" scale="1.0" value="%s"/>\n'%b
		free=False
	elif fixAll: #if set to fix all, should be only for crab IC component
		spec+='\t<!-- Source is kept fixed in 4FGL-DR3 analysis -->\n'
		spec+='\t\t<parameter free="0" max="1e6" min="0" name="norm" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		spec+='\t\t<parameter free="0" max="5.0" min="0.0" name="alpha" scale="1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="10.0" min="-5.0" name="beta" scale="1.0" value="%s"/>\n'%b
		free=False
	elif(dist>sL.radLim):
		if dist>sL.maxRad:
			spec+='\t<!-- Source is outside specified radius limit of %s -->\n'%sL.radLim
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="norm" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		elif vi<varValue or not sL.var:
			spec+='\t<!-- Source is outside specified radius limit of %s -->\n'%sL.radLim
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="norm" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		else:
			spec+='\t<!-- Source is outside specified radius limit of %s but variability index %.2f is greater than %.2f and varFree is set to True -->\n'%(sL.radLim,vi,varValue)
			spec+='\t\t<parameter free="1" max="1e6" min="0" name="norm" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=True
		spec+='\t\t<parameter free="0" max="5.0" min="0.0" name="alpha" scale="1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="10.0" min="-5.0" name="beta" scale="1.0" value="%s"/>\n'%b
	elif(TS<sL.sig):
		if vi<varValue or not sL.var:
			spec+='\t<!-- Source significance %.1f is less than specified minimum for a free source of %s -->\n'%(TS,sL.sig)
			spec+='\t\t<parameter free="0" max="1e6" min="0" name="norm" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=False
		else:
			spec+='\t<!-- Source significance %.1f is less than specified minimum for a free source of %s but variability index %.2f is greater than %.2f and varFree is set to True -->\n'%(TS,sL.sig,vi,varValue)
			spec+='\t\t<parameter free="1" max="1e6" min="0" name="norm" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
			free=True
		spec+='\t\t<parameter free="0" max="5.0" min="0.0" name="alpha" scale="1.0" value="%s"/>\n' %i
		spec+='\t\t<parameter free="0" max="10.0" min="-5.0" name="beta" scale="1.0" value="%s"/>\n'%b
	else:
		free=True
		spec+='\t\t<parameter free="1" max="1e6" min="0" name="norm" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
		if sL.nO:
			spec+='\t\t<parameter free="0" max="5.0" min="0.0" name="alpha" scale="1.0" value="%s"/>\n' %i
			spec+='\t\t<parameter free="0" max="10.0" min="-5.0" name="beta" scale="1.0" value="%s"/>\n'%b
		else:
			spec+='\t\t<parameter free="1" max="5.0" min="0.0" name="alpha" scale="1.0" value="%s"/>\n' %i
			spec+='\t\t<parameter free="1" max="10.0" min="-5.0" name="beta" scale="1.0" value="%s"/>\n'%b
	spec+='\t\t<parameter free="0" max="5e5" min="30" name="Eb" scale="1.0" value="%s"/>\n'%p
	spec+='\t</spectrum>\n'
	return spec,free

#this function searches the header of the ft1 to find the Position keyword and extract the ra and dec values
def getPos(ft1):
	file=pyfits.open(ft1)
	num=file[1].header['NDSKEYS']
	header=file[1].header
	right='POS(RA,DEC)'
	i=1
	keynum=0
	while i<=num:  #this step is necessary since it is not clear that the POS key word will have the same number always
		word='DSTYP%i' %i
		test=file[1].header[word]
		if(test==right):
			keynum=i
			i=num
		i+=1
	if(keynum==0):  #DSKEYS start numbering at 1, if this value hasn't been updated, KEYword doesn't exist
		print('Error: No position keyword found in fits header (assuming position is RA and DEC.  Exiting...')
		exit()
	keyword='DSVAL%i' %keynum
	try:
		ra,dec,rad=header[keyword].strip('CIRCLE()').split(',') #gets rid of the circle and parenthesis part and splits around the comma
		float(ra)
	except:
		ra,dec,rad=header[keyword].strip('circle()').split(',')
	file.close()
	return float(ra),float(dec),float(rad)
	
#calculates the angular separation between two points on the sky
def angsep(ra1,dec1,ra2,dec2):
	ra1*=d2r
	dec1*=d2r
	ra2*=d2r
	dec2*=d2r
	diffCosine=cos(dec1)*cos(dec2)*cos(ra1-ra2)+sin(dec1)*sin(dec2)
	dC='%.10f'%diffCosine#when the source is right at the center of the roi python sometimes adds extraneous digits at the end of the value i.e. instead of 1.0
	#it returns 1.0000000000000024, which throws an error with the acos function
	return acos(float(dC))/d2r #returns values between 0 and pi radians

#Check if a given file exists or not
def fileCheck(file):
	if (not os.access(file,os.F_OK)):
		return 0
	return 1

#copied from Damien's macro
def parameter_element(free, name, maximum, minimum, scale, value):
    """Create an XML document parameter description element"""
    impl = minidom.getDOMImplementation()
    xmldoc_out = impl.createDocument(None,None,None)
    parameter = xmldoc_out.createElement('parameter')
    parameter.setAttribute('free', str(free))
    parameter.setAttribute('name', str(name))
    parameter.setAttribute('max', str(maximum))
    parameter.setAttribute('min', str(minimum))
    parameter.setAttribute('scale', str(scale))
    parameter.setAttribute('value', str(value))
    return parameter

def mybool(Input):
	return {'True':True,'False':False,'T':True,'F':False,'t':True,'f':False,'TRUE':True,'FALSE':False,"true":True,"false":False,"1":True,"0":False}.get(Input)

def cli():
	import argparse
	
	helpString="Creates an xml model from the 3FGL catalog (FITS or xml version) for a specific ROI,\
		    coordinates of the ROI center are taken from an input event file,\
		    the radius for including sources is 10 degrees beyond the extraction radius used in the event file,\
		    sources with free parameters within the original extraction radius are chosen based on nearness to center, significance, and variability."
	parser=argparse.ArgumentParser(description=helpString)
	parser.add_argument("catalog",type=str,help="Catalog file to use, can be FITS or xml.")
	parser.add_argument("ev",type=str,help="Event file with ROI information in header.")
	parser.add_argument("-o","--outputxml",type=str,default='mymodel.xml',help="Name of output xml file, is set to overwrite files of same name.")
	parser.add_argument("-G","--galfile",type=str,default='$(FERMI_DIR)/refdata/fermi/galdiffuse/gll_iem_v07.fits',help="Name and location of Galactic diffuse model to use, will default to gll_iem_v06.fits.")
	parser.add_argument("-g","--galname",type=str,default='gll_iem_v07',help="Name of Galactic diffuse component in output model, will default to gll_iem_v06.")
	parser.add_argument("-I","--isofile",type=str,default='$(FERMI_DIR)/refdata/fermi/galdiffuse/iso_P8R3_SOURCE_V3_v1.txt',help="Name of isotropic diffuse template for output model, will default to P8R3 SOURCE class model.")
	parser.add_argument("-i","--isoname",type=str,default='iso_P8R3_SOURCE_V3_v1',help="Name of isotropic diffuse component in output model, default is for P8R3 SOURCE class.")
	parser.add_argument("-N","--normsonly",type=mybool,default=False,help="Flag to only let the normalizations of parameters be free, default is False.",nargs="?",const=True,choices=['True','False','T','F','t','f','TRUE','FALSE','true','false',1,0])
	parser.add_argument("-e","--extDir",type=str,default='',help="Path to directory with LAT extended source templates, will default to STs default.")#need to figure out what that is
	parser.add_argument("-r","--radLim",type=float,default=-1.,help="Radius, in degrees, from ROI center beyond which all source parameters should be fixed, will default to selection radius.")
	parser.add_argument("-R","--maxRad",type=float,default=None,help="Absolute maximum radius, in degrees, from ROI center beyond which all source parameters should be fixed, even variable sources will not be freed beyond this radius, defaults to radLim value.")
	parser.add_argument("-ER","--ExtraRad",type=float,default=10.,help="Radius beyond event file ROI out to which sources will be included in the model with all parameters fixed, default is 10, good for analyses starting around a few hundred MeV, can be decreased for high energy only fits.")
	parser.add_argument("-s","--sigFree",type=float,default=5.,help="Average significance below which all source parameters are fixed, defaults to 5.  Note, if using the 3FGL catalog xml file as input, this is actually a cut on TS, so adjust accordingly.")
	parser.add_argument("-v","--varFree",type=mybool,default=True,help="Flag to set normalization of significantly variable sources, even if source is beyond radius limit or below TS limit, default is True.",choices=['True','False','T','F','t','f','TRUE','FALSE','true','false',1,0])
	parser.add_argument("-p","--psForce",type=mybool,default=False,help="Flag to cast extended sources as point sources, default is False.",nargs="?",const=True,choices=['True','False','T','F','t','f','TRUE','FALSE','true','false',1,0])
	parser.add_argument("-E2C","--E2CAT",type=mybool,default=False,help="Flag to use catalog names for extended sources, default is False.",nargs="?",const=True,choices=['True','False','T','F','t','f','TRUE','FALSE','true','false',1,0])
	parser.add_argument("-m","--makeRegion",type=mybool,default=True,help="Flag to create ds9 region file as well as the xml model, default is True.",choices=['True','False','T','F','t','f','TRUE','FALSE','true','false',1,0])
	parser.add_argument("-GIF","--GIndexFree",type=mybool,default=False,help="Flag to use a power-law modification to the Galactic diffuse model spectrum and have the index be free, default is False.",nargs="?",const=True,choices=['True','False','T','F','t','f','TRUE','FALSE','true','false',1,0])
	parser.add_argument("-wd","--writeDir",type=str,default='',help="Directory to write the output ds9 region file in if not the current working directory or if you are specifying the full path to the newly made XML file.")
	parser.add_argument("-ON","--oldNames",type=mybool,default=False,help="Flag to use the make2FLGxml style naming convention, underscore before name and no spaces, default is False.",nargs="?",const=True,choices=['True','False','T','F','t','f','TRUE','FALSE','true','false',1,0])
	parser.add_argument('-DR','--DataRelease',type=int,default=3,help='Choice of data release, 3 for 4FGL-DR3, 2 for 4FGL-DR2, and 1 for 4FGL.',choices=[1,2,3])
	
	args=parser.parse_args()

	sL=srcList(args.catalog,args.ev,args.outputxml,args.DataRelease)
	
	sL.makeModel(args.galfile,args.galname,args.isofile,args.isoname,args.normsonly,args.extDir,args.radLim,args.maxRad,args.ExtraRad,args.sigFree,args.varFree,args.psForce,args.E2CAT,args.makeRegion,args.GIndexFree,args.writeDir,args.oldNames)
	
	
if __name__=='__main__': cli()