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'):
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
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
#GDfile is galactic diffuse model and ISOfile is optional Isotropic template file
#radLim is an angular radius limit beyond which all the spectral parameters of all sources will be fixed
#tslim, for 2FGL only, sources with TS below this limit will have spectral parameters fixed
#signif, for 1FGL only, sources with signiv_avg below this value will have spectral parameters fixed
#psForce, if true forces extended sources to be case as point sources
def makeModel(self,GDfile='$(GLAST_EXT)/extFiles/v0r9/galdiffuse/gll_iem_v02.fit',GDname='GAL_v02',ISOfile=None,ISOname='Extragalactic Diffuse',extDir='',radLim=-1,signif=4,psForce=False): #the version number may differ, but default syntax from ModelEditor didn't work
#quick check to see if 1FGL or 2FGL fits file
self.radLim=(self.roi[2] if radLim<=0 else radLim)
self.psF=psForce
self.extD=extDir
mycheck=pyfits.open(self.srcs)
try:
test=mycheck[1].data.field('Cutoff')
mycheck.close()
print 'Creating file and adding sources for 2FGL'
use=2
except:
mycheck.close()
print 'Creating file and adding sources for 1FGL'
use=1
if use==1:
addSrcs1(self,GDfile,GDname,ISOfile,ISOname,signif)
else:
addSrcs2(self,GDfile,GDname,ISOfile,ISOname,signif)
import pyfits
import os
from xml.dom.minidom import parseString as pS
from numpy import floor,log10,cos,sin,arccos,pi,array,log
acos=arccos
#import ROOT #note that this is only done to turn tab completion on for functions and filenames
print "This is make2FGLxml version 04r1."
print "For use with the gll_psc_v02.fit and gll_psc_v05.fit and later LAT catalog files."
#print "NOTE: You must have run gtselect on the event file you use as input."
d2r=pi/180.
#function to cycle through the source list and add point source entries
def addSrcs2(sL,GD,GDn,ISO,ISOn,signif):
model=open(sL.out,'w') #open file in write mode, overwrites other files of same name
file=pyfits.open(sL.srcs) #open source list file and access necessary fields, requires LAT source catalog definitions and names
data=file[1].data
extendedinfo=file[4].data
extName=extendedinfo.field('Source_Name')
extFile=extendedinfo.field('Spatial_Filename')
name=data.field('Source_Name')
ExtName=data.field('Extended_Source_Name')
ra=data.field('RAJ2000')
dec=data.field('DEJ2000')
flux=data.field('Flux_Density')
pivot=data.field('Pivot_Energy')
index=data.field('Spectral_Index')
cutoff=data.field('Cutoff')
spectype=data.field('SpectrumType')
beta=data.field('beta')
#extended=data.field('Extended')
sigma=data.field('Signif_Avg')
model.write('\n')
model.write('\n')
model.write('\n\n')
step=(sL.roi[2]+5)/5 #divide ROI radius plus 5 degrees into 5 steps for ordering of sources, get next highest int for ease
i=1
radii=[]
ptSrcNum=0
extSrcNum=0
while i<6:
if i*step<=sL.roi[2]+5:
radii+=[step*i]
else:
radii+=[sL.roi[2]+5] #just in case of rounding errors
i+=1
for x in radii:
if x==sL.roi[2]+5.:
model.write('\n\n' %(x-step,x))
else:
model.write('\n\n' %(x-step,x))
for n,f,i,r,d,p,c,t,b,S,EN in zip(name,flux,index,ra,dec,pivot,cutoff,spectype,beta,sigma,ExtName):
dist=angsep(sL.roi[0],sL.roi[1],r,d) #check that source is within ROI radius + 5 degress of ROI center
if r==sL.roi[0] and d==sL.roi[1]:
dist=0.0
if (dist=x-step) or (x==sL.roi[2]+5. and dist==x):
sn=''
for N in n.split(' '):
sn+=N
if EN!='' and not sL.psF:
extSrcNum+=1
if EN[0].isdigit():#use the extended source name to set these apart even further
Name='\n' %EN
else:
Name='\n' %EN
else:
ptSrcNum+=1
if sn[0].isdigit():
Name='\n' %sn
else:
Name='\n' %sn
if EN!='Vela X' and EN!='MSH 15-52':
if t=='PowerLaw':
spec=PLspec(sL,f,i,p,dist,S,signif)
elif t=='PowerLaw2':#no value for flux from 100 MeV to 100 GeV in fits file0
if i!=1.:#so calculate it by integrating PowerLaw spectral model
F=f*p**i/(-i+1)*(1.e5**(-i+1)-1.e2**(-1+1))
else:
F=f*p*log(1.e3)
spec=PL2spec(sL,F,i,dist,S,signif)
elif t=='LogParabola':
spec=LPspec(sL,f,i,p,b,dist,S,signif)
else:
spec=COspec(sL,f,i,p,c,dist,S,signif)
else:
if EN=='Vela X':
spec=VXspec(sL,i,dist)
else:
spec=MSHspec(sL,i,dist)
if EN!='' and not sL.psF:
eFile=None
for exN,exf in zip(extName,extFile):
if exN==EN:
if len(sL.extD)>0:
if sL.extD[-1]=='/':
eFile=sL.extD+exf
else:
eFile=sL.extD+'/'+exf
else:
eFile=exf
break
if eFile==None:
skydir='\t\n'%sL.extD
print 'Extended source %s (%s) in ROI, could not find name of matching template file.'%(EN,sn)
else:
skydir='\t\n'%eFile
print 'Extended source %s (%s) in ROI, verify that %s is the correct path.'%(EN,sn,eFile)
skydir+='\t\t\n'
skydir+='\t\n'
else:
skydir='\t\n'
skydir+='\t\t\n' %r
skydir+='\t\t\n' %d
skydir+='\t\n'
skydir+=''
(src,)=(Name+spec+skydir,)
ptsrc=pS(src).getElementsByTagName('source')[0]
ptsrc.writexml(model)
model.write('\n')
file.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\n')
Name='\n\n' %GDn
spec='\t\n'
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\n'
skydir='\t\n' %GD
skydir+='\t\t\n'
skydir+='\t\n'
skydir+=''
(src,)=(Name+spec+skydir,)
galdiff=pS(src).getElementsByTagName('source')[0]
galdiff.writexml(model)
model.write('\n')
if ISO==None: #null means no file so assume user wants an isotropic power law component
Name='\n' %ISOn
spec='\t\n'
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\n'
else: #if a file is given assume it is an isotropic template
Name='\n' %ISOn
spec='\t\n' %ISO
spec+='\t\t\n'
spec+='\t\n'
#both of the above options use the same spatial model
skydir='\t\n'
skydir+='\t\t\n'
skydir+='\t\n'
skydir+=''
(src,)=(Name+spec+skydir,)
iso=pS(src).getElementsByTagName('source')[0]
iso.writexml(model)
model.write('\n')
model.close()
return
def PLspec(sL,f,i,p,dist,sig,siglim):
fscale=int(floor(log10(f)))
spec='\t\n'
spec+='\t\n' %dist
if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
spec+='\t\n'
spec+='\t\t\n' %(fscale,f/10**fscale)
spec+='\t\t\n' %i
elif(dist>sL.radLim):
spec+='\t\n'%sL.radLim
spec+='\t\t\n' %(fscale,f/10**fscale)
spec+='\t\t\n' %i
elif sigsL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
spec+='\t\n'
spec+='\t\t\n'%(fscale,F/10**fscale)
spec+='\t\t\n' %i
elif(dist>sL.radLim):
spec+='\t\n'%sL.radLim
spec+='\t\t\n'%(fscale,F/10**fscale)
spec+='\t\t\n' %i
elif sigsL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
spec+='\t\n'
spec+='\t\t\n'
spec+='\t\t\n' %i
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\n'
return spec
def MSHspec(sL,i,dist):
spec='\t\n'
spec+='\t\n' %dist
if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
spec+='\t\n'
spec+='\t\t\n'
spec+='\t\t\n' %i
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\n'
return spec
def COspec(sL,f,i,p,c,dist,sig,siglim):
fscale=int(floor(log10(f)))
spec='\t\n'
spec+='\t\n' %dist
if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
spec+='\t\n'
spec+='\t\t\n' %(fscale,f/10**fscale)
spec+='\t\t\n' %i
if c<=1e5:
spec+='\t\t\n'%c
else:
spec+='\t\t\n'%(2.*c,c)
elif(dist>sL.radLim):
spec+='\t\n'%sL.radLim
spec+='\t\t\n' %(fscale,f/10**fscale)
spec+='\t\t\n' %i
if c<=1e5:
spec+='\t\t\n'%c
else:
spec+='\t\t\n'%(2.*c,c)
elif sigsL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
spec+='\t\n'
spec+='\t\t\n' %(fscale,f/10**fscale)
spec+='\t\t\n' %i
spec+='\t\t\n'%b
elif(dist>sL.radLim):
spec+='\t\n'%sL.radLim
spec+='\t\t\n' %(fscale,f/10**fscale)
spec+='\t\t\n' %i
spec+='\t\t\n'%b
elif sig\n')
model.write('\n')
model.write('\n\n')
step=(sL.roi[2]+5)/5 #divide ROI radius plus 5 degrees into 5 steps for ordering of sources, get next highest int for ease
i=1
radii=[]
while i<6:
if i*step<=sL.roi[2]+5:
radii+=[step*i]
else:
radii+=[sL.roi[2]+5] #just in case of rounding errors
i+=1
for x in radii:
if x==sL.roi[2]+5.:
model.write('\n\n' %(x-step,x))
else:
model.write('\n\n' %(x-step,x))
for n,f,i,r,d,p,s in zip(name,flux,index,ra,dec,pivot,sigma):
dist=angsep(sL.roi[0],sL.roi[1],r,d) #check that source is within ROI radius + 5 degress of ROI center
if (dist=x-step) or (x==sL.roi[2]+5. and dist==x):
fscale=int(floor(log10(f)))
sn=''
for N in n.split(' '):
sn+=N
if sn[0].isdigit():
Name='\n' %sn
else:
Name='\n' %sn
spec='\t\n'
spec+='\t\n' %dist
if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
spec+='\t\n'
spec+='\t\t\n' %(fscale,f/10**fscale)
spec+='\t\t\n' %i
elif(dist>sL.radLim):
spec+='\t\n'%sL.radLim
spec+='\t\t\n' %(fscale,f/10**fscale)
spec+='\t\t\n' %i
elif s\n')
Name='\n\n' %GDn
spec='\t\n'
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\n'
skydir='\t\n' %GD
skydir+='\t\t\n'
skydir+='\t\n'
skydir+=''
(src,)=(Name+spec+skydir,)
galdiff=pS(src).getElementsByTagName('source')[0]
galdiff.writexml(model)
model.write('\n')
if ISO==None: #null means no file so assume user wants an isotropic power law component
Name='\n' %ISOn
spec='\t\n'
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\t\n'
spec+='\t\n'
else: #if a file is given assume it is an isotropic template
Name='\n' %ISOn
spec='\t\n' %ISO
spec+='\t\t\n'
spec+='\t\n'
#both of the above options use the same spatial model
skydir='\t\n'
skydir+='\t\t\n'
skydir+='\t\n'
skydir+=''
(src,)=(Name+spec+skydir,)
iso=pS(src).getElementsByTagName('source')[0]
iso.writexml(model)
model.write('\n')
model.close()
return
#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