Add files via upload

fileprep_coord.py

@@ -0,0 +1,99 @@
+import numpy as np
+
+
+TYPE1='vespaBC_1_gals_SFH6HSsig_dr7_sp2'
+TYPE2='LINksmhigh22_sSFRzcut_parent'
+
+raGe,decGe,zGe,jackGe=np.loadtxt('{}_jack.dat'.format(TYPE1),unpack=True)
+raGl,decGl,zGl,jackGl=np.loadtxt('{}_jack.dat'.format(TYPE2),unpack=True)
+raRe,decRe,zRe,jackRe=np.loadtxt('/uufs/astro.utah.edu/common/uuastro/astro_data/zhengzheng/mccarthy/GravBinary_Data/random_{}_jack.dat'.format(TYPE1),unpack=True)
+raRl,decRl,zRl,jackRl=np.loadtxt('random_{}_jack.dat'.format(TYPE2),unpack=True)
+
+#lam=0.73
+#mat=0.27
+
+lam=0.692885
+mat=0.307115
+
+dGe=2000*(np.sqrt(lam+mat*(1+3*zGe))-1)/mat
+dGl=2000*(np.sqrt(lam+mat*(1+3*zGl))-1)/mat
+dRe=2000*(np.sqrt(lam+mat*(1+3*zRe))-1)/mat
+dRl=2000*(np.sqrt(lam+mat*(1+3*zRl))-1)/mat
+
+datae=np.vstack((raGe,decGe,dGe,jackGe)).T
+datal=np.vstack((raGl,decGl,dGl,jackGl)).T
+
+ra_Ge=raGe*np.pi/180.
+dec_Ge=np.pi/2.-decGe*np.pi/180.
+ra_Gl=raGl*np.pi/180.
+dec_Gl=np.pi/2.-decGl*np.pi/180.
+ra_Re=raRe*np.pi/180.
+dec_Re=np.pi/2.-decRe*np.pi/180.
+ra_Rl=raRl*np.pi/180.
+dec_Rl=np.pi/2.-decRl*np.pi/180.
+
+zmaxcut=0
+if zmaxcut==1:
+	cutg_early=dGe<=270.
+	cutr_early=dRe<=270.
+        cutg_late=dGl<=270.
+        cutr_late=dRl<=270.
+	dGe=dGe[cutg_early]
+	dRe=dRe[cutr_early]
+	ra_Ge=ra_Ge[cutg_early]
+	ra_Re=ra_Re[cutr_early]
+        dec_Ge=dec_Ge[cutg_early]
+        dec_Re=dec_Re[cutr_early]
+        dGl=dGl[cutg_late]
+        dRl=dRl[cutr_late]
+        ra_Gl=ra_Gl[cutg_late]
+        ra_Rl=ra_Rl[cutr_late]
+        dec_Gl=dec_Gl[cutg_late]
+        dec_Rl=dec_Rl[cutr_late]
+	jackGe=jackGe[cutg_early]
+	jackGl=jackGl[cutg_late]
+	jackRe=jackRe[cutr_early]
+	jackRl=jackRl[cutr_late]
+
+X_Ge=dGe*np.sin(dec_Ge)*np.cos(ra_Ge)
+Y_Ge=dGe*np.sin(dec_Ge)*np.sin(ra_Ge)
+Z_Ge=dGe*np.cos(dec_Ge)
+
+X_Gl=dGl*np.sin(dec_Gl)*np.cos(ra_Gl)
+Y_Gl=dGl*np.sin(dec_Gl)*np.sin(ra_Gl)
+Z_Gl=dGl*np.cos(dec_Gl)
+
+X_Re=dRe*np.sin(dec_Re)*np.cos(ra_Re)
+Y_Re=dRe*np.sin(dec_Re)*np.sin(ra_Re)
+Z_Re=dRe*np.cos(dec_Re)
+
+X_Rl=dRl*np.sin(dec_Rl)*np.cos(ra_Rl)
+Y_Rl=dRl*np.sin(dec_Rl)*np.sin(ra_Rl)
+Z_Rl=dRl*np.cos(dec_Rl)
+
+dataGe_jack=np.vstack((X_Ge,Y_Ge,Z_Ge,jackGe)).T
+dataGl_jack=np.vstack((X_Gl,Y_Gl,Z_Gl,jackGl)).T
+dataRe_jack=np.vstack((X_Re,Y_Re,Z_Re,jackRe)).T
+dataRl_jack=np.vstack((X_Rl,Y_Rl,Z_Rl,jackRl)).T
+
+dataGe=np.vstack((X_Ge,Y_Ge,Z_Ge)).T
+dataGl=np.vstack((X_Gl,Y_Gl,Z_Gl)).T
+dataRe=np.vstack((X_Re,Y_Re,Z_Re)).T
+dataRl=np.vstack((X_Rl,Y_Rl,Z_Rl)).T
+
+
+#np.savetxt('earlyLIN_SFH_CART.dat',dataGe, fmt=['%5e','%5e','%5e'])
+np.savetxt('{}_CART.dat'.format(TYPE1),dataGe, fmt=['%5e','%5e','%5e'])
+np.savetxt('{}_CART_jack.dat'.format(TYPE1),dataGe_jack, fmt=['%5e','%5e','%5e','%5d'])
+
+#np.savetxt('lateLIN_SFH_CART.dat',dataGl, fmt=['%5e','%5e','%5e'])
+#np.savetxt('{}_CART.dat'.format(TYPE2),dataGl, fmt=['%5e','%5e','%5e'])
+#np.savetxt('{}_CART_jack.dat'.format(TYPE2),dataGl_jack, fmt=['%5e','%5e','%5e','%5d'])
+
+#np.savetxt('random_SFH_earlyLIN_CART.dat',dataRe, fmt=['%5e','%5e','%5e'])
+np.savetxt('/uufs/astro.utah.edu/common/uuastro/astro_data/zhengzheng/mccarthy/GravBinary_Data/random_{}_CART.dat'.format(TYPE1),dataRe, fmt=['%5e','%5e','%5e'])
+np.savetxt('/uufs/astro.utah.edu/common/uuastro/astro_data/zhengzheng/mccarthy/GravBinary_Data/random_{}_CART_jack.dat'.format(TYPE1),dataRe_jack, fmt=['%5e','%5e','%5e','%5d'])
+
+#np.savetxt('random_SFH_lateLIN_CART.dat',dataRl, fmt=['%5e','%5e','%5e'])
+#np.savetxt('random_{}_CART.dat'.format(TYPE2),dataRl, fmt=['%5e','%5e','%5e'])
+#np.savetxt('random_{}_CART_jack.dat'.format(TYPE2),dataRl_jack, fmt=['%5e','%5e','%5e','%5d'])

fileprep_jack_assign.py

@@ -0,0 +1,251 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+####### From https://gist.github.com/eteq/4599814
+try:
+	from scipy.spatial import cKDTree as KDT
+except ImportError:
+	from scipy.spatial import KDTree as KDT
+
+def spherematch(ra1, dec1, ra2, dec2, tol=None, nnearest=1):
+    """
+    Finds matches in one catalog to another.
+    Parameters
+    ra1 : array-like
+        Right Ascension in degrees of the first catalog
+    dec1 : array-like
+        Declination in degrees of the first catalog (shape of array must match `ra1`)
+    ra2 : array-like
+        Right Ascension in degrees of the second catalog
+    dec2 : array-like
+        Declination in degrees of the second catalog (shape of array must match `ra2`)
+    tol : float or None, optional
+        How close (in degrees) a match has to be to count as a match.  If None,
+        all nearest neighbors for the first catalog will be returned.
+    nnearest : int, optional
+        The nth neighbor to find.  E.g., 1 for the nearest nearby, 2 for the
+        second nearest neighbor, etc.  Particularly useful if you want to get
+        the nearest *non-self* neighbor of a catalog.  To do this, use:
+        ``spherematch(ra, dec, ra, dec, nnearest=2)``
+    Returns
+    -------
+    idx1 : int array
+        Indecies into the first catalog of the matches. Will never be
+        larger than `ra1`/`dec1`.
+    idx2 : int array
+        Indecies into the second catalog of the matches. Will never be
+        larger than `ra1`/`dec1`.
+    ds : float array
+        Distance (in degrees) between the matches
+    """
+
+    ra1 = np.array(ra1, copy=False)
+    dec1 = np.array(dec1, copy=False)
+    ra2 = np.array(ra2, copy=False)
+    dec2 = np.array(dec2, copy=False)
+
+    if ra1.shape != dec1.shape:
+        raise ValueError('ra1 and dec1 do not match!')
+    if ra2.shape != dec2.shape:
+        raise ValueError('ra2 and dec2 do not match!')
+
+    x1, y1, z1 = _spherical_to_cartesian(ra1.ravel(), dec1.ravel())
+
+    # this is equivalent to, but faster than just doing np.array([x1, y1, z1])
+    coords1 = np.empty((x1.size, 3))
+    coords1[:, 0] = x1
+    coords1[:, 1] = y1
+    coords1[:, 2] = z1
+
+    x2, y2, z2 = _spherical_to_cartesian(ra2.ravel(), dec2.ravel())
+
+    # this is equivalent to, but faster than just doing np.array([x1, y1, z1])
+    coords2 = np.empty((x2.size, 3))
+    coords2[:, 0] = x2
+    coords2[:, 1] = y2
+    coords2[:, 2] = z2
+
+    kdt = KDT(coords2)
+    if nnearest == 1:
+        idxs2 = kdt.query(coords1)[1]
+    elif nnearest > 1:
+        idxs2 = kdt.query(coords1, nnearest)[1][:, -1]
+    else:
+        raise ValueError('invalid nnearest ' + str(nnearest))
+
+    ds = _great_circle_distance(ra1, dec1, ra2[idxs2], dec2[idxs2])
+
+    idxs1 = np.arange(ra1.size)
+
+    if tol is not None:
+        msk = ds < tol
+        idxs1 = idxs1[msk]
+        idxs2 = idxs2[msk]
+        ds = ds[msk]
+
+    return idxs1, idxs2, ds
+
+
+def _spherical_to_cartesian(ra, dec):
+    """
+    (Private internal function)
+    Inputs in degrees.  Outputs x,y,z
+    """
+    rar = np.radians(ra)
+    decr = np.radians(dec)
+
+    x = np.cos(rar) * np.cos(decr)
+    y = np.sin(rar) * np.cos(decr)
+    z = np.sin(decr)
+
+    return x, y, z
+
+
+def _great_circle_distance(ra1, dec1, ra2, dec2):
+    """
+    (Private internal function)
+    Returns great circle distance.  Inputs in degrees.
+    Uses vicenty distance formula - a bit slower than others, but
+    numerically stable.
+    """
+    from numpy import radians, degrees, sin, cos, arctan2, hypot
+
+    # terminology from the Vicenty formula - lambda and phi and
+    # "standpoint" and "forepoint"
+    lambs = radians(ra1)
+    phis = radians(dec1)
+    lambf = radians(ra2)
+    phif = radians(dec2)
+
+    dlamb = lambf - lambs
+
+    numera = cos(phif) * sin(dlamb)
+    numerb = cos(phis) * sin(phif) - sin(phis) * cos(phif) * cos(dlamb)
+    numer = hypot(numera, numerb)
+    denom = sin(phis) * sin(phif) + cos(phis) * cos(phif) * cos(dlamb)
+    return degrees(arctan2(numer, denom))
+########################################
+
+################# From Z.Z.
+def binarySearch(alist, item):
+	first = 0
+	last = len(alist)-1
+	found = False
+	while first<=last and not found:
+	        midpoint = (first + last)//2
+	        if alist[midpoint] == item:
+	            found = True
+	        else:
+	            if item < alist[midpoint]:
+	                last = midpoint-1
+	            else:
+	                first = midpoint+1
+	if found:
+		return found,midpoint
+	else:
+		return found, 00
+##########################################
+
+
+TYPE1='vespaBC_1_gals_SFH6HSsig_dr7_sp2'
+TYPE2='LINksmhigh22_sSFRzcut_parent'
+
+# Randoms from Hong, early/late created by Kevin
+raR,decR,jackR=np.loadtxt('random.dat',unpack=True)
+data_early=np.loadtxt('{}.dat'.format(TYPE1),unpack=False)
+data_late=np.loadtxt('{}.dat'.format(TYPE2),unpack=False)
+
+#TYPE1='LINearly_sSFR'
+#TYPE2='LINlate_sSFR'
+
+raGe=data_early[:,0]
+decGe=data_early[:,1]
+zGe=data_early[:,2]
+
+raGl=data_late[:,0]
+decGl=data_late[:,1]
+zGl=data_late[:,2]
+
+# Check footprint 
+plt.scatter(raR,decR,color='y')
+plt.scatter(raGl,decGl,color='b',s=1,alpha=0.5)
+plt.show()
+
+plt.scatter(raR,decR,color='y')
+plt.scatter(raGe,decGe,color='r',s=1,alpha=0.5)
+plt.show()
+
+# Assign jack id to early/late galaxies
+indGe,indRe,ds=spherematch(raGe,decGe,raR,decR,tol=1,nnearest=1)
+indGl,indRl,ds=spherematch(raGl,decGl,raR,decR,tol=1,nnearest=1)
+
+jackGe=jackR[indRe]
+jackGl=jackR[indRl]
+
+# Assign redshift to random catalogs
+Nr=len(raR)
+
+znum=100
+histe,edgee=np.histogram(zGe, znum)
+
+randmulte=Nr/len(zGe)
+hist1e=np.around(histe*randmulte)
+binwidthe=edgee[1]-edgee[0]
+
+histl,edgel=np.histogram(zGl, znum)
+
+randmultl=Nr/len(zGl)
+hist1l=np.around(histl*randmultl)
+binwidthl=edgel[1]-edgel[0]
+
+
+zrande=[]
+for i in range(0,znum):
+        if i <znum:
+                Ne=int(hist1e[i])
+                for j in range(0,Ne):
+                        if j < Ne:
+                                ztmpe=np.random.rand()*binwidthe+edgee[i]
+                                zrande.append(ztmpe)
+
+se=len(zrande)
+raRe=raR[0:se]
+decRe=decR[0:se]
+jackRe=jackR[0:se]
+catRe=np.vstack((raRe,decRe,zrande,jackRe)).T
+#np.savetxt('/uufs/astro.utah.edu/common/uuastro/astro_data/zhengzheng/mccarthy/GravBinary_Data/random_{}_jack.dat'.format(TYPE1),catRe,fmt=['%5e','%5e','%5e','%5d'])
+
+
+zrandl=[]
+for i in range(0,znum):
+	if i <znum:
+		Nl=int(hist1l[i])
+		for j in range(0,Nl):
+			if j < Nl:
+				ztmpl=np.random.rand()*binwidthl+edgel[i]
+				zrandl.append(ztmpl)
+
+sl=len(zrandl)
+raRl=raR[0:sl]
+decRl=decR[0:sl]
+jackRl=jackR[0:sl]
+catRl=np.vstack((raRl,decRl,zrandl,jackRl)).T
+#np.savetxt('/uufs/astro.utah.edu/common/uuastro/astro_data/zhengzheng/mccarthy/GravBinary_Data/random_{}_jack.dat'.format(TYPE2),catRl,fmt=['%5e','%5e','%5e','%5d'])
+
+plt.hist(zrandl,bins=100,normed=True,label='randoms')
+plt.hist(zGl,bins=100,normed=True,alpha=0.5,label='late')
+plt.legend(loc='upper right')
+plt.show()
+
+plt.hist(zrande,bins=100,normed=True,label='randoms')
+plt.hist(zGe,bins=100,normed=True,alpha=0.5,label='early')
+plt.legend(loc='upper right')
+plt.show()
+
+
+
+earlyjack=np.vstack((raGe,decGe,zGe,jackGe)).T
+latejack=np.vstack((raGl,decGl,zGl,jackGl)).T
+#np.savetxt('{}_jack.dat'.format(TYPE1),earlyjack)
+#np.savetxt('{}_jack.dat'.format(TYPE2),latejack)