Python WLanalysis.writeFits方法代码示例

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def KSmap(iinput):
	'''Input:
	i = ith zbin for zcut
	hl = 'hi' or 'lo' for higher/lower z of the zcut
	sigmaG: smoothing scale
	Wx = 1..4 of the field
	Output:
	smoothed KS map and galn map.
	'''
	Wx, sigmaG, i, hl = iinput
	print 'Wx, sigmaG, i, hl:', Wx, sigmaG, i, hl
	kmap_fn = cat_dir+'KS/W%i_KS_%s_%s_sigmaG%02d.fit'%(Wx, zbins[i],hl,sigmaG*10)
	galn_smooth_fn = cat_dir+'KS/W%i_galn_%s_%s_sigmaG%02d.fit'%(Wx, zbins[i],hl,sigmaG*10)
	
	isfile_kmap, kmap = WLanalysis.TestFitsComplete(kmap_fn, return_file = True)
	if isfile_kmap == False:
		Me1_fn = cat_dir+'Me_Mw_galn/W%i_Me1w_%s_%s.fit'%(Wx, zbins[i],hl)
		Me2_fn = cat_dir+'Me_Mw_galn/W%i_Me2w_%s_%s.fit'%(Wx, zbins[i],hl)
		Mw_fn = cat_dir+'Me_Mw_galn/W%i_Mwm_%s_%s.fit'%(Wx, zbins[i],hl)
		Me1 = WLanalysis.readFits(Me1_fn)
		Me2 = WLanalysis.readFits(Me2_fn)
		Mw = WLanalysis.readFits(Mw_fn)	
		Me1_smooth = WLanalysis.weighted_smooth(Me1, Mw, PPA=PPA512, sigmaG=sigmaG)
		Me2_smooth = WLanalysis.weighted_smooth(Me2, Mw, PPA=PPA512, sigmaG=sigmaG)
		kmap = WLanalysis.KSvw(Me1_smooth, Me2_smooth)
		WLanalysis.writeFits(kmap,kmap_fn)
	isfile_galn, galn_smooth = WLanalysis.TestFitsComplete(galn_smooth_fn, return_file = True)
	if isfile_galn == False:
		galn_fn = cat_dir+'Me_Mw_galn/W%i_galn_%s_%s.fit'%(Wx, zbins[i],hl)
		galn = WLanalysis.readFits(galn_fn)
		galn_smooth = WLanalysis.smooth(galn, sigma=sigmaG*PPA512)
		WLanalysis.writeFits(galn_smooth, galn_smooth_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def OrganizeSplitFile(ifile):
	'''read in one of the split file, pick out the redshift, and sort by fields, e2 is c2 correted, with e2-=c2'''	
	
	field = genfromtxt(split_dir+ifile,usecols=0,dtype=str)
	field = array(map(field2int,field))
	print ifile

	# generate 2 random redshift and 1 peak
	Pz = genfromtxt(split_dir+ifile,usecols=arange(14,84),dtype=str)
	Pz = (np.core.defchararray.replace(Pz,',','')).astype(float)
	
	seed(99)
	z_rand1 = array(map(DrawRedshifts,Pz)).ravel()
	seed(88)
	z_rand2 = array(map(DrawRedshifts,Pz)).ravel()
	z_peak = z_arr[argmax(Pz,axis=1)]
	z_all = concatenate([[z_peak,], [z_rand1,], [z_rand2,]]).T

	sheardata = genfromtxt(split_dir+ifile,usecols=[1,2,5,6,7,8,9,10,11,12,13,84])
	ra, dec, e1, e2, w, fitclass, r, snr, mask, m, c2, mag = sheardata.T
	e2 -= c2
	
	
	i=0
	for Wx in range(1,5):
		idx=where((field==Wx)&(mask<=1.0)&(fitclass==0)&(amin(z_all,axis=-1)>=0.2)&(amax(z_all,axis=-1)<=1.3))[0]
		print ifile, Wx, len(idx)/50000.0
		if len(idx) > 0:
			#data = (np.array([ra,dec,e1,e2,w,r,snr,m,c2,mag]).T)[idx]
			data = (np.array([ra,dec,e1,e2,w,r,snr,m,c2,mag,z_peak, z_rand1, z_rand2]).T)[idx]
			radeclist = sheardata[idx][:,[0,1]]
			xylist = list2coords(radeclist, Wx)
			xy_data = concatenate([xylist,data],axis=1)
			WLanalysis.writeFits(xy_data, W_dir(Wx)+ifile+'.fit')#,fmt=['%i','%i','%s','%s','%s','%.3f'])
		i+=1

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def Bmode(iinput):
	'''Input:
	i = ith zbin for zcut
	hl = 'hi' or 'lo' for higher/lower z of the zcut
	sigmaG: smoothing scale
	Wx = 1..4 of the field
	Output:
	smoothed KS map and galn map.
	'''
	Wx, sigmaG, i, hl = iinput
	print 'Bmode - Wx, sigmaG, i, hl:', Wx, sigmaG, i, hl
	bmap_fn = cat_dir+'KS/W%i_Bmode_%s_%s_sigmaG%02d.fit'%(Wx, zbins[i],hl,sigmaG*10)
	#galn_smooth_fn = cat_dir+'KS/W%i_galn_%s_%s_sigmaG%02d.fit'%(Wx, zbins[i],hl,sigmaG*10)
	
	isfile_kmap, bmap = WLanalysis.TestFitsComplete(bmap_fn, return_file = True)
	if isfile_kmap == False:
		Me1_fn = cat_dir+'Me_Mw_galn/W%i_Me1w_%s_%s.fit'%(Wx, zbins[i],hl)
		Me2_fn = cat_dir+'Me_Mw_galn/W%i_Me2w_%s_%s.fit'%(Wx, zbins[i],hl)
		Mw_fn = cat_dir+'Me_Mw_galn/W%i_Mwm_%s_%s.fit'%(Wx, zbins[i],hl)
		Me1 = WLanalysis.readFits(Me1_fn)
		Me2 = WLanalysis.readFits(Me2_fn)
		Mw = WLanalysis.readFits(Mw_fn)	
		Me1_smooth = WLanalysis.weighted_smooth(Me1, Mw, PPA=PPA512, sigmaG=sigmaG)
		Me2_smooth = WLanalysis.weighted_smooth(Me2, Mw, PPA=PPA512, sigmaG=sigmaG)
		### Bmode conversion is equivalent to
		### gamma1 -> gamma1' = -gamma2
		### gamma2 -> gamma2' = gamma1
		bmap = WLanalysis.KSvw(-Me2_smooth, Me1_smooth)
		WLanalysis.writeFits(bmap,bmap_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def KSmap (i):
	'''
	Smooth and KS inversion
	input: i = 1, 2, 3 ... 13
	output: nothing, write kmap, Mmask if haven't done so
	'''
	Me1, Me2, Mw, galn = fileGen(i)
	for sigmaG in sigmaG_arr:	
		print 'KSmap i, sigmaG', i, sigmaG
		KS_fn = KS_dir+'CFHT_KS_sigma%02d_subfield%02d.fits'%(sigmaG*10,i)
		mask_fn = '/scratch/02977/jialiu/KSsim/mask/CFHT_mask_ngal%i_sigma%02d_subfield%02d.fits'%(ngal_arcmin,sigmaG*10,i)
		
		if WLanalysis.TestComplete((KS_fn,mask_fn),rm=True):
			kmap = WLanalysis.readFits(KS_fn)
			Mmask = WLanalysis.readFits(mask_fn)
		else:
			Me1_smooth = WLanalysis.weighted_smooth(Me1, Mw, PPA=PPA512, sigmaG=sigmaG)
			Me2_smooth = WLanalysis.weighted_smooth(Me2, Mw, PPA=PPA512, sigmaG=sigmaG)
			galn_smooth = snd.filters.gaussian_filter(galn.astype(float),sigmaG*PPA512, mode='constant')
			## KS
			kmap = WLanalysis.KSvw(Me1_smooth, Me2_smooth)
			## mask
			maskidx = where(galn_smooth < ngal_cut) #cut at ngal=5
			Mmask = ones(shape=galn.shape)
			Mmask[maskidx]=0
			
			WLanalysis.writeFits(kmap, KS_fn)
			WLanalysis.writeFits(Mmask, mask_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def average_powspec_nonoise (cosmo):
	ps = zeros(shape=(1000,50))
	weights = (genfromtxt(KSsim_dir+'galn.txt').T[1]).astype(float)
	weights /= sum(weights)
	fn = KSsim_dir+'powspec_Mk_sum13fields/SIM_powspec_sigma05_rz1_%s_1000R.fit'%(cosmo)
	if os.path.isfile(fn):
		return WLanalysis.readFits(fn)
	else:
		for i in range(1,14):
			ips=weights[i-1]*WLanalysis.readFits(powspecMk_fn(i, cosmo))
			ps += ips
		WLanalysis.writeFits(ps,fn)
		return ps

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def Psingle_CFHT (i, sigmaG, bins, ps=0):
	if ps:
		fn = powspec_CFHT_fn(i, sigmaG)
	else:
		fn = peaks_CFHT_fn(i, sigmaG, bins)
	
	if os.path.isfile(fn):
		out = WLanalysis.readFits(fn)
	elif ps:
		kmap = WLanalysis.readFits(KSCFHT_fn(i, sigmaG))
		out = WLanalysis.PowerSpectrum(kmap, sizedeg=12.0)
		WLanalysis.writeFits(out,fn)
	else:
		kmap = WLanalysis.readFits(KSCFHT_fn(i, sigmaG))
		mask = WLanalysis.readFits(Mask_fn(i, sigmaG))
		out = WLanalysis.peaks_mask_hist(kmap, mask, bins, kmin=kmin, kmax=kmax)
		WLanalysis.writeFits(out,fn)
	return out

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def createBadFieldMask (sf):
	sf_splitfiles = os.listdir(sf_dir(sf))
	genfromtxtA = lambda fn: genfromtxt(sf_dir(sf)+fn)
	datas = map(genfromtxtA,sf_splitfiles)#3 columns: RA, DEC, GB
	datas = concatenate(datas,axis=0)
	idx = where(datas[:,-1]==1)[0]
	datas = datas[idx]
	y, x, k = datas.T
	k, galn = WLanalysis.coords2grid(x, y, array([k,]))
	for sigmaG in sigmaG_arr:
		print 'createBadFieldMask sf, sigmaG:', sf, sigmaG
		Allmask = WLanalysis.readFits(mask_fcn(sigmaG, sf))#mask for all field
		badmask_fn = badmask_fcn(sigmaG, sf)#file name for bad pointing mask, which is 75% area of Allmask
		galn_smooth = snd.filters.gaussian_filter(galn.astype(float),sigmaG*PPA512, mode='constant')
		#smooth the galn grid
		Mmask = ones(shape=galn.shape)#create mask grid
		Mmask[where(galn_smooth < ngal_cut)]=0#find the low density region in galn_smooth
		Mmask = adHocFix(Mmask,sf)
		Mmask *= Allmask#since I didn't do redshift cut in badmask, so here it takes care of it, since ALl mask has redshift cuts
		WLanalysis.writeFits(Mmask, badmask_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def average_powspec_withnoise (cosmo, sigmaG, zg='rz1', CFHT=None):
	weights = (genfromtxt(KSsim_dir+'galn.txt').T[1]).astype(float)
	weights /= sum(weights)
	
	if CFHT:
		ps = zeros(50)
		fn = KSsim_dir+'powspec_sum13fields/CFHT_powspec_sigma%02d.fit'%(sigmaG*10)
	else:
		ps = zeros(shape=(1000,50))
		fn = KSsim_dir+'powspec_sum13fields/SIM_powspec_sigma%02d_%s_%s_%04dR.fit'%(sigmaG*10, zg, cosmo, Rtol)
	if os.path.isfile(fn):
		return WLanalysis.readFits(fn)
	else:
		for i in range(1,14):
			if CFHT:
				ips=weights[i-1]*WLanalysis.readFits(powspec_CFHT_fn(i, sigmaG))[-1]
			else:
				ips=weights[i-1]*WLanalysis.readFits(powspec_fn(i, cosmo, 1000, sigmaG, zg))
			ps += ips
		WLanalysis.writeFits(ps,fn)
		return ps

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def Noise(iinput):
	'''Input: (Wx, iseed)
	Return: files of noise KS map, using randomly rotated galaxy.
	'''
	Wx, iseed = iinput
	seed(iseed)
	print 'Bmode - Wx, iseed:', Wx, iseed
	bmap_fn = cat_dir+'Noise/W%i/W%i_Noise_sigmaG10_%04d.fit'%(Wx, Wx, iseed)
	
	isfile_kmap, bmap = WLanalysis.TestFitsComplete(bmap_fn, return_file = True)
	if isfile_kmap == False:
		Me1_fn = cat_dir+'Me_Mw_galn/W%i_Me1w_1.3_lo.fit'%(Wx)
		Me2_fn = cat_dir+'Me_Mw_galn/W%i_Me2w_1.3_lo.fit'%(Wx)
		Mw_fn = cat_dir+'Me_Mw_galn/W%i_Mwm_1.3_lo.fit'%(Wx)
		Me1_init = WLanalysis.readFits(Me1_fn)
		Me2_init = WLanalysis.readFits(Me2_fn)
		#### randomly rotate Me1, Me2 ###
		Me1, Me2 = WLanalysis.rndrot(Me1_init, Me2_init)
		#################################
		Mw = WLanalysis.readFits(Mw_fn)	
		Me1_smooth = WLanalysis.weighted_smooth(Me1, Mw, PPA=PPA512, sigmaG=sigmaG)
		Me2_smooth = WLanalysis.weighted_smooth(Me2, Mw, PPA=PPA512, sigmaG=sigmaG)
		bmap = WLanalysis.KSvw(Me1_smooth, Me2_smooth)
		WLanalysis.writeFits(bmap,bmap_fn)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def Pmat (iRcosmo, Rtol=Rtol, R0 = 1):
	'''
	Input:
	iRcosmo = (i, bins, sigmaG, cosmo, zg)
	Rtol: total number of realizations, and count peaks for realizations #(R0, R0+1, .. R0+Rtotl-1)
	R0: the first realization, if not starting from 1
	if bins = 0 return power spectrum, else return peak counts.
	Return:
	A maxtrix of shape=(Rtol x bins)
	'''
	i, sigmaG, zg, bins, cosmo = iRcosmo
	if bins == 0:#powspec
		fn = powspec_fn(i, cosmo, Rtol, sigmaG, zg)
	else:#peaks
		fn = peaks_fn(i, cosmo, Rtol, sigmaG, zg, bins)
	if os.path.isfile(fn):
		mat = WLanalysis.readFits(fn)
	else:
		print 'Pmat - i, bins, sigmaG', i, bins, sigmaG
		map_fcn = Psingle (i, sigmaG, zg, bins, cosmo, pk=bins)
		#p = Pool(Rtol/4)#use multiprocessing on 1 single core
		mat = array(map(map_fcn,R_arr))
		WLanalysis.writeFits(mat, fn)
	return mat

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def fileGen(i):
	'''
	Input:
	i range from (1, 2..13)
	Return:
	Me1 = e1*w
	Me2 = (e2-c2)*w
	Mw = (1+m)*w
	galn = number of galaxies per pixel
	'''
	Me1_fn = KS_dir+'CFHT_subfield%02d_Me1.fits'%(i)
	Me2_fn = KS_dir+'CFHT_subfield%02d_Me2.fits'%(i)
	Mw_fn = KS_dir+'CFHT_subfield%02d_Mw.fits'%(i)
	galn_fn = KS_dir+'CFHT_subfield%02d_galn.fits'%(i)
	
	print 'fileGen', i
	if WLanalysis.TestComplete((Me1_fn,Me2_fn,Mw_fn,galn_fn),rm = True):
		Me1 = WLanalysis.readFits(Me1_fn)
		Me2 = WLanalysis.readFits(Me2_fn)
		Mw =  WLanalysis.readFits(Mw_fn)
		galn =WLanalysis.readFits(galn_fn)
	else:
		ifile = np.genfromtxt(full_dir+'full_subfield'+str(i) ,usecols=[0, 1, 2, 3, 4, 9, 10, 11, 16, 17])
		# cols: y, x, z_peak, z_rnd1, z_rnd2, e1, e2, w, m, c2

		#redshift cut 0.2< z <1.3
		zs = ifile[:,[2,3,4]]
		print 'zs'
		idx = np.where((amax(zs,axis=1) <= zmax) & (amin(zs,axis=1) >= zmin))[0]
		
		y, x, z_peak, z_rnd1, z_rnd2, e1, e2, w, m, c2 = ifile[idx].T

		k = array([e1*w, (e2-c2)*w, (1+m)*w])
		Ms, galn = WLanalysis.coords2grid(x, y, k)
		print 'coords2grid'
		Me1, Me2, Mw = Ms
		WLanalysis.writeFits(Me1,Me1_fn)
		WLanalysis.writeFits(Me2,Me2_fn)
		WLanalysis.writeFits(Mw,Mw_fn)
		WLanalysis.writeFits(galn,galn_fn)
	return Me1, Me2, Mw, galn

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
def SumSplitFile2Grid(Wx):
	'''For Wx field, read in each split file, 
	and create e1, e2 grid for mass construction.
	Input: Wx=1,2,3,4
	Output: (Me1, Me2, Mw, galn) split in each redshift bins'''
	isize = sizes[Wx-1]
	ishape = (len(zbins), isize, isize)
	ishape_hi = (len(zbins)-1, isize, isize)#no need to do hi for zcut=1.3 since it's everything
	Me1_hi = zeros(shape=ishape_hi)
	Me2_hi = zeros(shape=ishape_hi)#hi is for higher redshift bins, lo is lower redshift
	Mw_hi = zeros(shape=ishape_hi)
	#Mk_hi = zeros(shape=ishape_hi)
	galn_hi = zeros(shape=ishape_hi)

	Me1_lo = zeros(shape=ishape)
	Me2_lo = zeros(shape=ishape)
	Mw_lo = zeros(shape=ishape)
	#Mk_lo = zeros(shape=ishape)
	galn_lo = zeros(shape=ishape)
	
	Wfiles = os.listdir(W_dir(Wx))#get the list of split file for Wx
	for iW in Wfiles:
		datas = WLanalysis.readFits(W_dir(Wx)+iW)
		#cols: x, y, ra, dec, e1, e2, w, r, snr, m, c2, mag, z_peak, z_rand1, z_rand2
		z = datas.T[-3]#z_peak, -2 is z_rand1, -1 is z_rand2
		i = 0 #zbin count
		for zcut in zbins:
			idx0 = where(z<zcut)[0]
			idx1 = where(z>=zcut)[0]
			for idx in [idx0,idx1]:
				y, x, e1, e2, w, m = (datas[idx].T)[[0,1,4,5,6,9]]#note x, y is reversed in python
				k = array([e1*w, e2*w, (1+m)*w])
				x = radians(x)
				y = radians(y)
				print 'W'+str(Wx), iW, 'coords2grid, zbin =',zbins[i]
				A, galn = WLanalysis.coords2grid(x, y, k, size=isize)
				if len(idx)==0:#no need to calculate hi bin for zcut=1.3
					continue
				elif idx[0] == idx0[0]:
					Me1_lo[i] += A[0]
					Me2_lo[i] += A[1]
					Mw_lo[i] += A[2]
					galn_lo[i] += galn
				
				else:
					Me1_hi[i] += A[0]
					Me2_hi[i] += A[1]
					Mw_hi[i] += A[2]
					galn_hi[i] += galn
			i+=1
	print 'Done collecting small fields for W'+str(Wx)
	
	for i in range(len(zbins)):
		for hl in ('lo','hi'):
			Me1_fn = cat_dir+'Me_Mw_galn/W%i_Me1w_%s_%s.fit'%(Wx, zbins[i],hl)
			Me2_fn = cat_dir+'Me_Mw_galn/W%i_Me2w_%s_%s.fit'%(Wx, zbins[i],hl)
			Mw_fn = cat_dir+'Me_Mw_galn/W%i_Mwm_%s_%s.fit'%(Wx, zbins[i],hl)
			galn_fn = cat_dir+'Me_Mw_galn/W%i_galn_%s_%s.fit'%(Wx, zbins[i],hl)
			if hl=='hi' and i==len(zbins)-1:
				continue
			elif hl=='lo':
				WLanalysis.writeFits(Me1_lo[i],Me1_fn, rewrite = True)
				WLanalysis.writeFits(Me2_lo[i],Me2_fn, rewrite = True)
				WLanalysis.writeFits(Mw_lo[i],Mw_fn, rewrite = True)
				WLanalysis.writeFits(galn_lo[i],galn_fn, rewrite = True)
			else:
				WLanalysis.writeFits(Me1_hi[i],Me1_fn, rewrite = True)
				WLanalysis.writeFits(Me2_hi[i],Me2_fn, rewrite = True)
				WLanalysis.writeFits(Mw_hi[i],Mw_fn, rewrite = True)
				WLanalysis.writeFits(galn_hi[i],galn_fn, rewrite = True)

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
hi_m='mQ3-512b240_Om0.290_Ol0.710_w-1.000_ns0.960_si0.800'
cosmo_arr=(fidu,hi_m,hi_w,hi_s)

def galn_gen(i):
	print i
	y, x, z = WLanalysis.readFits(emucat_dir+'emulator_subfield%i_zcut0213.fit'%(i)).T
	Mz, galn = WLanalysis.coords2grid(x, y, np.array([z,]))
	WLanalysis.writeFits(galn, galn_fn)
	
# map(galn_gen,range(1,14))

# power spectrum for all of the rz1, 4 cosmo
def ps (R):
	Mk = WLanalysis.readFits(Mk_fn(i, cosmo, R))
	galn = WLanalysis.readFits(galn_fn(i))
	Mk_smooth = WLanalysis.weighted_smooth(Mk, galn, sigmaG)
	ps = WLanalysis.PowerSpectrum(Mk_smooth)
	return ps[1]
	
for i in range(1, 14):
	for cosmo in cosmo_arr:
		print i, cosmo
		p = MPIPool()
		pmat = np.array(p.map(ps, arange(1,1001)))
		pmat_fn = KS_dir + 'powspec_Mk/SIM_powspec_sigma05_subfield%i_rz1_%s_1000R.fit'%(i,cosmo)
		p.close()
		try:
			WLanalysis.writeFits(pmat,pmat_fn)
		except Exception:
			pass
print 'done'

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
			else:
				print 'generating KS'
				y, x, e1, e2, w = WLanalysis.readFits(test_dir+'yxew_subfield%i_zcut0213.fit'%(i)).T
				e1_45, e2_45 = WLanalysis.rndrot(e1, e2, deg=45)
				e1_rand, e2_rand = WLanalysis.rndrot(e1, e2, iseed=0)
				mat_e1_45,mat_e2_45,mat_e1_rand,mat_e2_rand, Mw = WLanalysis.coords2grid(x, y, array([e1_45*w, e2_45*w, e1_rand*w, e2_rand*w, w]) )[0]

				mat_e1_45_smoothed  = WLanalysis.weighted_smooth(mat_e1_45  , Mw, sigmaG=sG)
				mat_e2_45_smoothed  = WLanalysis.weighted_smooth(mat_e2_45  , Mw, sigmaG=sG)
				mat_e1_rand_smoothed= WLanalysis.weighted_smooth(mat_e1_rand, Mw, sigmaG=sG)
				mat_e2_rand_smoothed= WLanalysis.weighted_smooth(mat_e2_rand, Mw, sigmaG=sG)
				KS_45 =WLanalysis.KSvw(mat_e1_45_smoothed,mat_e2_45_smoothed)
				KS_rand=WLanalysis.KSvw(mat_e1_rand_smoothed,mat_e2_rand_smoothed)
				
				
				WLanalysis.writeFits(KS_45,KS_45_fn)
				WLanalysis.writeFits(KS_rand,KS_rand_fn)
			
			kappa_45 = KS_45.flatten()
			std_45=std(kappa_45)
			#hist_45,binedges45 = histogram(kappa_45,bins=100,range=(-5*std_45,5*std_45))
			#savetxt(test_dir+'hist45_%i_%i.ls'%(i,sG),array([hist_45,binedges45[:-1]]).T)
			
			kappa_rand = KS_rand.flatten()
			#std_rand=std(kappa_rand)
			#hist_rand,binedgesrand = histogram(kappa_rand,bins=100,range=(-5*std_rand,5*std_rand))
			#savetxt(test_dir+'histrand_%i_%i.ls'%(i,sG),array([hist_rand,binedgesrand[:-1]]).T)
			
			peaks_45 = WLanalysis.peaks_list(KS_45)
			peaks_rand=WLanalysis.peaks_list(KS_rand)
			

# 需要导入模块: import WLanalysis [as 别名]
# 或者: from WLanalysis import writeFits [as 别名]
powspec_CFHT_fn = lambda i, sigmaG: KSCFHT_dir+'CFHT_powspec_sigma%02d_subfield%02d.fits'%(sigmaG*10, i)

peaks_fn = lambda i, cosmo, Rtol, sigmaG, zg, bins: KSsim_dir+'peaks/SIM_peaks_sigma%02d_subfield%i_%s_%s_%04dR_%03dbins.fit'%(sigmaG*10, i, zg, cosmo, Rtol, bins)

peaks_sum_fn = lambda cosmo, Rtol, sigmaG, zg, bins: KSsim_dir+'peaks_sum13fields/SIM_peaks_sigma%02d_%s_%s_%04dR_%03dbins.fit'%(sigmaG*10, zg, cosmo, Rtol, bins)

for bins in bins_arr:
	for sigmaG in sigmaG_arr:
		print bins, sigmaG
		
		fn_CFHT = KSsim_dir+'peaks_sum13fields/CFHT_peaks_sigma%02d_%03dbins.fits'%(sigmaG*10, bins)
		CFHT_peaks = zeros(shape=(Rtol,bins))
		for i in i_arr:
			CFHT_peaks += WLanalysis.readFits(peaks_CFHT_fn(i, sigmaG, bins))
		try:
			WLanalysis.writeFits(CFHT_peaks, fn_CFHT)
		except Exception:
			print fn_CFHT,'already exist, but no worries'
			pass
				
		
		for zg in zg_arr:
			for cosmo in cosmo_arr:
				fn = peaks_sum_fn (cosmo, Rtol, sigmaG, zg, bins)
				peaks = zeros(shape=(Rtol,bins))
				for i in i_arr:
					peaks += WLanalysis.readFits(peaks_fn(i, cosmo, Rtol, sigmaG, zg, bins))
				try:
					WLanalysis.writeFits(peaks, fn)
				except Exception:
					fn, ' already exist, but no worries'