Skip to content

Commit

Permalink
Initial Commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@nfshermy
nfshermy committed Jun 5, 2018
0 parents commit d294e09
Show file tree
Hide file tree
Showing 62 changed files with 11,979 additions and 0 deletions.
493 changes: 493 additions & 0 deletions HTML.py

Large diffs are not rendered by default.

120 changes: 120 additions & 0 deletions HostMatchAlgorithms_Ravi.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
#!/usr/bin/env python
"""
Chris D'Andrea, Ravi Gupta
21 July 2015
The routines hereing are called by the desHostMatch.py script, which the directional
light radius (DLR) to a galaxy from an RA and DEC given either
a. Second moments (get_DLR)
b. A_IMAGE, B_IMAGE, THETA_IMAGE (get_DLR_ABT)
v0.1: July 2015 -- Split off host-match functions into a separate script to be imported.
Makes swapping out the modules simpler.
v0.2: June 2016 -- Add function to compute HOST_CONFUSION from DLR (Ravi)
"""

import sys
import os
import numpy as np
import math as m

##################### DEFINE GLOBAL VARIABLES #################################################

rad = np.pi/180 # convert deg to rad
pix_arcsec = 0.264 # pixel scale (arcsec per pixel)
pix2_arcsec2 = 0.264**2 # pix^2 to arcsec^2 conversion factor
pix2_deg2 = pix2_arcsec2/(3600**2) # pix^2 to deg^2 conversion factor


####################### Function to compute DLR distance from 2nd moments #########################

def get_DLR(RA_SN, DEC_SN, RA, DEC, X2, Y2, XY, angsep):
# inputs are arrays

global numFailed
rPHI = np.empty_like(angsep)
d_DLR = np.empty_like(angsep)

# convert from IMAGE units (pixels^2) to WORLD (arcsec^2)
X2_ARCSEC = X2*pix2_arcsec2
Y2_ARCSEC = Y2*pix2_arcsec2
XY_ARCSEC = XY*pix2_arcsec2

dX2Y2 = X2_ARCSEC - Y2_ARCSEC # difference
sX2Y2 = X2_ARCSEC + Y2_ARCSEC # sum

# ensures not all 2nd moments are same & arctan is not undefined later
good = ( (dX2Y2 + sX2Y2 -2*XY_ARCSEC != 0) &
((np.fabs(XY_ARCSEC) > 1.e-6) | (np.fabs(dX2Y2) > 1.e-6)) )
bad = np.invert(good)

# spatial rms of object profile along semi-major axis. If bad, set to 0.0
A_ARCSEC = np.where(good, np.sqrt(0.5*sX2Y2 + np.sqrt((0.5*dX2Y2)**2 + XY_ARCSEC**2)), 0.0)
# spatial rms of object profile along semi-minor axis. If bad, set to 0.0
B_ARCSEC = np.where(good, np.sqrt(0.5*sX2Y2 - np.sqrt((0.5*dX2Y2)**2 + XY_ARCSEC**2)), 0.0)

# if bad, the default THETA is 0.0
THETA = np.where(good, 0.5*np.arctan2(2.0*XY_ARCSEC, dX2Y2), 0.0)
# angle between RA-axis and SN-host vector
GAMMA = np.arctan((DEC_SN - DEC)/(np.cos(DEC_SN*rad)*(RA_SN - RA)))
# angle between semi-major axis of host and SN-host vector
PHI = THETA + GAMMA

rPHI = np.where(good, A_ARCSEC*B_ARCSEC/np.sqrt((A_ARCSEC*np.sin(PHI))**2 +
(B_ARCSEC*np.cos(PHI))**2), 0.0)

# directional light radius
# where 2nd moments are bad, set d_DLR = 99.99
d_DLR = np.where(good, angsep/rPHI, 99.99)

return [d_DLR, A_ARCSEC, B_ARCSEC, rPHI]

######################## Function to compute DLR distance from ####################################
######################## A_IMAGE, B_IMAGE, THETA_IMAGE ############################################

def get_DLR_ABT(RA_SN, DEC_SN, RA, DEC, A_IMAGE, B_IMAGE, THETA_IMAGE, angsep):
# inputs are arrays

global numFailed
rPHI = np.empty_like(angsep)
d_DLR = np.empty_like(angsep)

# convert from IMAGE units (pixels) to WORLD (arcsec^2)
A_ARCSEC = A_IMAGE*pix_arcsec
B_ARCSEC = B_IMAGE*pix_arcsec

# angle between RA-axis and SN-host vector
GAMMA = np.arctan((DEC_SN - DEC)/(np.cos(DEC_SN*rad)*(RA_SN - RA)))
# angle between semi-major axis of host and SN-host vector
PHI = THETA_IMAGE + GAMMA # angle between semi-major axis of host and SN-host vector

rPHI = A_ARCSEC*B_ARCSEC/np.sqrt((A_ARCSEC*np.sin(PHI))**2 +
(B_ARCSEC*np.cos(PHI))**2)

# directional light radius
# where 2nd moments are bad, set d_DLR = 99.99
d_DLR = angsep/rPHI

return [d_DLR, A_ARCSEC, B_ARCSEC, rPHI]

################### Function to compute host confusion from DLR ##################################

def compute_HC(DLR):
dlr = np.sort(DLR) # first sort DLR array
e = 1e-5 # define epsilon, some small number
if len(dlr)==1: # only 1 object in radius
HC = -99.0
else:
delta12 = dlr[1] - dlr[0] + e
D1D2 = dlr[0]**2/dlr[1] + e
Dsum = 0
for i in range(0, len(dlr)):
for j in range(i+1, len(dlr)):
didj = dlr[i]/dlr[j] + e
delta_d = dlr[j] - dlr[i] + e
Dsum += didj/((i+1)**2*delta_d)
HC = m.log10(D1D2*Dsum/delta12)
return HC

############################### MAIN ##########################################################

228 changes: 228 additions & 0 deletions WSTemp_check.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,228 @@
import numpy as np
import pandas as pd
import easyaccess as ea
from glob import glob
import os
import sys
from collections import Counter
import fnmatch

base = '/pnfs/des/persistent/gw/exp'

ls = '/data/des41.b/data/rbutler/sb/bench/PostProcessing/event2_trim.list'
biglist = '/data/des41.b/data/rbutler/desgw/event2test/exposures.list'
colnames = ['expnum','nite','mjd_obs','radeg','decdeg','band','exptime','propid','obstype','object']


#df = pd.read_csv(biglist,delim_whitespace=True)#,header=0,names=colnames)

#expnum = np.genfromtxt(biglist,delimiter=[6],dtype=['int'],skip_header=1,usecols=(0))

#print expnum[:10]
#print nite[:10]

#print df[:10]
#print df.ix[:10,'band']
#print list(df)

#sys.exit()

if os.path.exists(ls):
f = open(ls,'r')
exps = f.readlines()
f.close()
exps = map(lambda s: int(s.strip()), exps)
else:
exps = [497767,498267,498268]

exps = exps[:]
#exps = [506423,506424,506425,506426,506427,506428,506429,506430,506431,506432]
#exps = [506430]
season = 208
bands = 'i'

if bands=='all':
bands = None
else:
bands = bands.split(',')

print exps

connection = ea.connect('desoper')
cursor = connection.cursor()

if len(exps)>1:
exps = str(tuple(exps))
query = 'select expnum, band, nite from exposure where expnum in '+exps
else:
query = 'select expnum, band, nite from exposure where expnum='+str(exps[0])

df = connection.query_to_pandas(query)

df = df.loc[df['BAND'].isin(bands)]
df = df.sort_values(by='EXPNUM')
df = df.reset_index(drop=True)

length = len(df)
bad = []
total = 0
notemp = 0
fznotfound = 0
starnotfound = 0
other = 0

chips = range(1,63)
chips.remove(2),chips.remove(31),chips.remove(61)

for ind in range(len(df)):
expnum = str(df.get_value(ind,'EXPNUM'))
nite = str(df.get_value(ind,'NITE'))
band = str(df.get_value(ind,'BAND'))
path = os.path.join(base,nite)
path = os.path.join(path,expnum)
path = os.path.join(path,'dp'+str(season))
print
print str(ind+1)+'/'+str(length)+' - '+path
print
for c in chips:
c = '%02d' % c
cpath = os.path.join(path,band+'_'+c)
#print cpath
if os.path.isdir(cpath):
ffail = 'RUN03_expose_mkWStemplate.FAIL'
allfiles = os.listdir(cpath)
for fp in allfiles:
if fnmatch.fnmatch(fp,'WSTemplate_*_GW*GWV1*.stdout'):
stdout = fp
if ffail in allfiles:
#if 1==1:
f = open(os.path.join(cpath,fp))
lines = f.readlines()
f.close()
#print expnum,c,lines[-1]
if '***FATAL' in lines[-1]:
total+=1
if '[.fz] not found' in lines[-1]:
fznotfound+=1
#print 'fz',os.path.join(cpath,stdout)
elif '*/ not found' in lines[-1]:
print expnum,c
for lin in lines:
if 'DEBUG: exposure' in lin:

starnotfound+=1
#print '*/',os.path.join(cpath,stdout)
elif 'No good templates.' in lines[-1]:
print expnum,c
notemp+=1
#print 'NGT',os.path.join(cpath,stdout)
else:
for i in lines[-3:]:
print i
other+=1


# for ind in range(len(df)):
# expnum = str(df.get_value(ind,'EXPNUM'))
# nite = str(df.get_value(ind,'NITE'))
# band = str(df.get_value(ind,'BAND'))
# path = os.path.join(base,nite)
# path = os.path.join(path,expnum)
# path = os.path.join(path,'dp'+str(season))
# print
# print str(ind+1)+'/'+str(length)+' - '+path
# print
# for c in chips:
# c = '%02d' % c
# cpath = os.path.join(path,band+'_'+c)
# #print cpath
# if os.path.isdir(cpath):
# globout = glob(os.path.join(cpath,'WSTemplate_*_GW*GWV1*.stdout'))
# pathfail = os.path.join(cpath,'RUN03_expose_mkWStemplate.FAIL')
# if len(globout)==1 and os.path.isfile(pathfail):
# #print globout[0]
# f = open(globout[0],'r')
# lines = f.readlines()
# f.close()
# #print lines[-1]
# if '***FATAL' in lines[-1]:
# total+=1
# if '[.fz] not found' in lines[-1]:
# fznotfound+=1
# print 'fz',globout[0]
# #for g in lines[-1:]:
# #print g
# #print '-'*45
# elif '*/ not found' in lines[-1]:
# starnotfound+=1
# print '*/',globout[0]
# elif 'No good templates.' in lines[-1]:
# notemp+=1
# else:
# other+=1
# #print e, c
# #print lines[-1]
# #for l in lines[-3:]:
# # print l
# #print
# #bad.append(int(lines[-12].split('_')[1]))
# #bad.append(lines[-12])

# for e in exps:
# print str(exps.index(e)+1)+'/'+str(length)+' - '+str(e)
# path = os.path.join(base,str(e))
# path = os.path.join(path,'dp'+str(season))
# for c in range(1,63):
# ccd = '%02d' % c
# globout = glob(path+'/*'+ccd+'/WSTemplate_*_GW*GWV1*.stdout')
# globfail = glob(path+'/*'+ccd+'/RUN03_expose_mkWStemplate.FAIL')
# #print globout
# if len(globout)==1 and len(globfail)==1:
# #print globout[0]
# f = open(globout[0],'r')
# lines = f.readlines()
# f.close()
# #print lines[-1]
# if '***FATAL' in lines[-1]:
# total+=1
# if 'not found' in lines[-1]:
# notfound+=1
# print globout[0]
# for g in lines[-1:]:
# print g
# print '-'*45
# elif 'No good templates.' in lines[-1]:
# notemp+=1
# else:
# other+=1
# #print e, c
# #print lines[-1]
# #for l in lines[-3:]:
# # print l
# #print
# #bad.append(int(lines[-12].split('_')[1]))
# #bad.append(lines[-12])

print
print '---'
print

print total
print '[.fz] not found: %d, %.1f%%' % (fznotfound,((float(fznotfound)/float(total))*100))
print '...*/ not found: %d, %.1f%%' % (starnotfound,((float(starnotfound)/float(total))*100))
print 'no good templates: %d, %.1f%%' % (notemp,((float(notemp)/float(total))*100))
print 'other: %d, %.1f%%' % (other,((float(other)/float(total))*100))

sys.exit()

counts = Counter(bad)

dcounts = dict(counts)
#print dcounts

for i in dcounts:
print i,':',dcounts[i]

total = len(bad)
print
print 'total :',total
Loading

0 comments on commit d294e09

Please sign in to comment.