import os import numpy as np import re import sys import fitsio import math import time from pydl import pydlutils import maskutils import mksampleutils def fbcorrect(dd,dindx,FBopt,legacyopt,fbmatchgood,fbmatchbad,fbradarcsec=62.): """ Fiber correct galaxy with index dindx Right now only a single option defined here, may want to update this later. For instance, CMASS allowed to collide with LOWZ?? """ ## sanity check that object is really in a collision group. assert dd['IMATCH'][dindx] == 0 ### doesn't pass when I run on ICOLLIDED = 0 set. #assert dd['ICOLLIDED'][dindx] == 1 #assert dd['MULTGROUP'][dindx] > 1 ## does this work? if not dd['MULTGROUP'][dindx] > 1: return -1, -1, -1 degtorad = np.pi/180. fbrad = fbradarcsec/3600.*degtorad*1.0001 fbraddeg = fbradarcsec/3600.*1.0001 legacynbr = 0 if FBopt == 0: ## get group members. gg = np.where((dd['FINALN'] == dd['FINALN'][dindx]) & (dd['INGROUP'] == dd['INGROUP'][dindx]))[0] assert len(gg) >= 1 ## find closest CMASS. minsep = 1e8 nbrindx = -1 for g in gg: if g==dindx: continue ## you found yourself. if dd['IMATCH'][g] in fbmatchbad: continue asepdeg = mksampleutils.angsepdeg(dd['RA'][dindx],dd['DEC'][dindx], dd['RA'][g],dd['DEC'][g]) if(asepdeg < minsep): nbrindx = g minsep = asepdeg assert dd['IMATCH'][g] in fbmatchgood assert dd['IMATCH'][g] != 0 assert dd['IMATCH'][g] != 2 if asepdeg > 3.*fbraddeg: print dindx,g,'funny business??',dd['MULTGROUP'][dindx],len(gg),asepdeg/fbraddeg if nbrindx > -1: assert dd['IMATCH'][nbrindx] in fbmatchgood if nbrindx == -1: return -1, -1, -1 ## don't fiber correct this guy. if legacyopt == 1: ll = np.where(dd['IMATCH'] == 2)[0] ## in pydl version on Utah, this gives an exception: ## File "/uufs/chpc.utah.edu/sys/pkg/python/2.7.3_rhel6/lib/python2.7/site-packages/pydl-0.3.0.dev310-py2.7.egg/pydl/pydlutils/spheregroup/spherematch.py", line 48, in spherematch ## raise PydlutilsException("Change the order of the sets of coordinates!") ## pydl.pydlutils.PydlutilsException: Change the order of the sets of coordinates! #ddi, ddimatch2i, dist = pydlutils.spheregroup.spherematch(np.array([dd['RA'][dindx]]), np.array([dd['DEC'][dindx]]), dd['RA'][ll], dd['DEC'][ll],matchlength = fbraddeg,maxmatch=len(ll)+1) ddimatch2i, ddi, dist = pydlutils.spheregroup.spherematch(dd['RA'][ll], dd['DEC'][ll], np.array([dd['RA'][dindx]]), np.array([dd['DEC'][dindx]]),matchlength = fbraddeg,maxmatch=len(ll)+1) if len(dist) > 0: ## found a match to legacy object. bestl = ll[ddimatch2i[np.where(dist == dist.min())[0][0]]] assert (ddi == 0).all() angsepchk = mksampleutils.angsepdeg(dd['RA'][dindx],dd['DEC'][dindx], dd['RA'][bestl],dd['DEC'][bestl]) assert angsepchk <= fbraddeg assert (angsepchk - dist.min()) < 1.0e-3 if angsepchk < minsep: nbrindx = bestl minsep = angsepchk legacynbr = 1 return nbrindx, minsep, legacynbr else: print 'FBopt = ',FBopt,'not yet supported, aborting!' sys.exit(1) def legacycollide(dd,maskcpfits,fbseed,fbrad=62./3600.): np.random.seed(seed=fbseed) maskcp = maskutils.mask(maskcpfits) #ddi, ddzfaili, dist = pydlutils.spheregroup.spherematch(np.array([dd['RA'][zi]]), np.array([dd['DEC'][zi]]), dd['RA'][mgal], dd['DEC'][mgal], matchlength = zfraddeg, maxmatch=len(mgal)+1) #sgroup = pydlutils.spheregroup.spheregroup(dd['RA'][xx], dd['DEC'][xx], linklength = fbrad) ## index for legacy objects. li = np.where(dd['IMATCH'] == 2)[0] ## index of objects to match to: not 0, not 3. check later for 8. mi = np.where((dd['IMATCH'] != 0) & (dd['IMATCH'] != 3))[0] ddmi, ddli, dist = pydlutils.spheregroup.spherematch(dd['RA'][mi], dd['DEC'][mi], dd['RA'][li], dd['DEC'][li], matchlength = fbrad, maxmatch=len(mi)*len(li)+1) ## first randomize the order in which you go through xx = np.arange(0,len(dist),1) np.random.shuffle(xx) assert len(xx) == len(dist) dist = dist[xx] ddmi = ddmi[xx] ddli = ddli[xx] legacypairslist = [] npairstot = 0. npairsbad = 0. npairsexpected = 0. ## go through all the spherematches, remove collisions with some probability based on NTILES. for ii in range(len(ddmi)): i1 = mi[ddmi[ii]] i2 = li[ddli[ii]] if i1 == i2: continue ## same object! if dd['IMATCH'][i1] == 8: continue if dd['IMATCH'][i2] == 8: continue assert dd['IMATCH'][i2] == 2 if dd['IMATCH'][i1] == 2: ## only test each legacy-legacy pair once. if (i1,i2) in legacypairslist or (i2,i1) in legacypairslist: continue legacypairslist.append((i1,i2)) npairstot += 1 if dd['ISECT'][i1] == dd['ISECT'][i2]: pval = maskcp.weight[dd['ISECT'][i1]] else: ## choose randomly which sector. if np.random.random() < 0.5: pval = maskcp.weight[dd['ISECT'][i1]] else: pval = maskcp.weight[dd['ISECT'][i2]] npairsexpected += (1.-pval) ## probability in mask is probability of resolving pairs. So we remove pairs with probabiliy 1-pval. if np.random.random() > pval: ## randomly choose which to delete. npairsbad += 1 if np.random.random() < 0.5: dd['IMATCH'][i1] = 8 dd['Z'][i1] = dd['Z'][i2] dd['WEIGHT_CP'][i2] += dd['WEIGHT_CP'][i1] else: dd['IMATCH'][i2] = 8 dd['Z'][i2] = dd['Z'][i1] dd['WEIGHT_CP'][i1] += dd['WEIGHT_CP'][i2] print 'legacy collide finished' print npairstot, npairsbad, npairsexpected, npairsbad/float(npairstot) return 0 #def fibercollide(dd, binarymaskfits, FBopt=0, legacyopt = 0, fbmatchgood=[1,4,7],fbmatchbad=[0,2,3,5,6]): def fiberloop(dd, collidelist, FBopt=0, legacyopt = 0, fbmatchgood=[1,4,7],fbmatchbad=[0,2,3,5,6]): """ fbmatchgood are things that are allowed to receive fiber collision weight. fbmatchbad are things that are NOT allowed to receive fiber collision weight. FBopt is the fiber collision correction algorithm choice. legacyopt lets legacy galaxies receive close pair weight if they're closer than closest regular target. """ xx = collidelist ## set of indices to look for collisions for. fbnbrtot = 0 fbnonbrtot = 0 fbnbrlegacy = 0 cnt = 0 t1 = time.time() seplist = [] for x in xx: if(cnt%100 == 0): print cnt,'of',len(xx),'done.' nbrindx, minsep, legacynbr = fbcorrect(dd,x,FBopt,legacyopt,fbmatchgood,fbmatchbad) if nbrindx != -1: dd['IMATCH'][x] = 3 dd['Z'][x] = dd['Z'][nbrindx] dd['WEIGHT_CP'][nbrindx] += 1 fbnbrtot += 1 seplist.append(minsep) if legacynbr == 1: fbnbrlegacy += 1 else: fbnonbrtot += 1 cnt += 1 print 'assigned this many fb collisions: ',fbnbrtot,'this many were legacy:',fbnbrlegacy print 'did not assign nbr to this many icollided objects:',fbnonbrtot if len(seplist) > 0: print 'max separation assigned in fb units',max(seplist)/(62./3600.) hsep, xseptmp = np.histogram(seplist,bins=np.arange(0.,1.05*max(seplist),1./3600.)) xsep = 0.5*(xseptmp[1:] + xseptmp[:-1]) print 'start histogram of close pair separations' for ii in range(len(hsep)): print xsep[ii], hsep[ii] print 'end histogram of close pair separations' else: print 'No fiber collision corrections done.' t2 = time.time() print 'FB collide took',t2-t1,'seconds' def groupkey(ddi): """ Input element of group array. Output groupid key. """ return (ddi['FINALN'], ddi['INGROUP']) ## copied from match_sample.py def photokey(darr): """ SDSS magic photoid is some combo of these things, so they can be used as a unique key for each object. """ return (darr['RUN'], darr['RERUN'], darr['CAMCOL'], darr['FIELD'], darr['ID']) def calc_dzhist(ddcN, ddcS, ddlN, ddlS, binarymaskfitslowzcut): binarymask = maskutils.mask(binarymaskfitslowzcut) dd = np.concatenate((ddcN, ddcS, ddlN, ddlS)) cmassmax = len(ddcN) + len(ddcS) FBrad = 62./3600. ## now make dictionary lookup for groups. dgroups = {} maxlength = 0 for di in range(len(dd)): if dd['INGROUP'][di] == -1: continue if dd['IMATCH'][di] == 2: continue ## don't let 2's in groups, new May 25. ## only retain galaxies that are inside the binary mask. if binarymask.weight[dd['ISECT'][di]] < 0.01: continue mykey = groupkey(dd[di]) if mykey in dgroups: ## Jan 27 2015 ## only append if it's unique! there is CMASS/LOWZ targeting overlap at the boundary. ## first check that pkeynew = photokey(dd[di]) dupflag = 0 for jj in dgroups[mykey]: pkeyold = photokey(dd[jj]) if pkeyold == pkeynew: ## yes, these are duplicates. # print 'duplicate?',dd['RA'][di], dd['RA'][jj], dd['DEC'][di], dd['DEC'][jj] dupflag = 1 if dupflag == 0: dgroups[mykey].append(di) if len(dgroups[mykey]) > maxlength: maxlength = len(dgroups[mykey]) else: dgroups[mykey] = [di] ##trim dictionary to include only groups with at least two CMASS entries. print 'this many items in dict:',len(dgroups) dgroups = {key: np.array(value) for key, value in dgroups.items() if len(value) > 1} print 'after purge, this many items in dict:',len(dgroups) print 'max group size',maxlength ofpc = open('cmassFBzdist.txt','w') ofpl = open('lowzFBzdist.txt','w') for kk,groupindx in dgroups.items(): ci = np.where(groupindx < cmassmax)[0] li = np.where(groupindx >= cmassmax)[0] if len(ci) > 1: for ii in range(len(ci)): for jj in range(ii+1,len(ci)): if dd['IMATCH'][groupindx[ci[ii]]] == 1 and dd['IMATCH'][groupindx[ci[jj]]] == 1: ofpc.write('%d %d %e %e\n' % (len(ci), len(li), dd['Z'][groupindx[ci[ii]]], dd['Z'][groupindx[ci[jj]]])) if len(li) > 1: for ii in range(len(li)): for jj in range(ii+1,len(li)): if dd['IMATCH'][groupindx[li[ii]]] == 1 and dd['IMATCH'][groupindx[li[jj]]] == 1: ofpl.write('%d %d %e %e\n' % (len(ci), len(li), dd['Z'][groupindx[li[ii]]], dd['Z'][groupindx[li[jj]]])) ofpc.close() ofpl.close() def calc_cp_comp_new(ddcN, ddcS, ddlN, ddlS, binarymaskfitslowzcut, fibergood=[1,4,7]): """ copying calc_cp_comp, we want to recompute fiber collision stats on complete (non-duplicated) sample. Make collision group catalog without duplicates. """ binarymask = maskutils.mask(binarymaskfitslowzcut) print 'check binary!',(binarymask.weight*binarymask.area).sum() maskcpqpm = maskutils.mask(binarymaskfitslowzcut) maskcppair = maskutils.mask(binarymaskfitslowzcut) dd = np.concatenate((ddcN, ddcS, ddlN, ddlS)) ## these keep their order, so last cmass indx is len(ddc)-1; useful later cmassmax = len(ddcN) + len(ddcS) FBrad = 62./3600. ## histogram of all group sizes. ## now make dictionary lookup for groups. dgroups = {} maxlength = 0 for di in range(len(dd)): if dd['INGROUP'][di] == -1: continue if dd['IMATCH'][di] == 2: continue ## don't let 2's in groups, new May 25. ## only retain galaxies that are inside the binary mask. if binarymask.weight[dd['ISECT'][di]] < 0.01: continue mykey = groupkey(dd[di]) if mykey in dgroups: ## Jan 27 2015 ## only append if it's unique! there is CMASS/LOWZ targeting overlap at the boundary. ## first check that pkeynew = photokey(dd[di]) dupflag = 0 for jj in dgroups[mykey]: pkeyold = photokey(dd[jj]) if pkeyold == pkeynew: ## yes, these are duplicates. # print 'duplicate?',dd['RA'][di], dd['RA'][jj], dd['DEC'][di], dd['DEC'][jj] dupflag = 1 if dupflag == 0: dgroups[mykey].append(di) if len(dgroups[mykey]) > maxlength: maxlength = len(dgroups[mykey]) else: dgroups[mykey] = [di] ##trim dictionary to include only groups with at least two CMASS entries. print 'this many items in dict:',len(dgroups) dgroups = {key: np.array(value) for key, value in dgroups.items() if len(value) > 1} print 'after purge, this many items in dict:',len(dgroups) print 'max group size',maxlength ## compute histogram of number of gals in each size. hC = np.zeros(maxlength+1) hL = np.zeros(maxlength+1) htot = np.zeros(maxlength+1) for kk,groupindx in dgroups.items(): ci = np.where(groupindx < cmassmax)[0] li = np.where(groupindx >= cmassmax)[0] if len(ci) > 0: hC[len(ci)] += len(ci) if len(li) > 0: hL[len(li)] += len(li) htot[len(groupindx)] += len(groupindx) ## write tex for table. ofp = open('table3dr12v4.tex','w') ofp.write(r"$n_{\rm group}$ & $f_{\rm CMASS}$ & $f_{\rm LOWZ}$ & $f_{\rm C+L}$" + "\\\\\n") for i in range(2,10): ofp.write('%d & %.4f & %.4f & %.4f\\\\\n' % (i, hC[i]/hC[2:].sum(), hL[i]/hL[2:].sum(), htot[i]/htot[2:].sum())) ofp.close() ## now copying from calc_cp_comp directly. maxsect = maskcpqpm.mf['SECTOR'].max() ngroupfiber = np.zeros(maxsect + 1,dtype='int') ngroupnofiber = np.zeros(maxsect + 1,dtype='int') npairresolved = np.zeros([maxlength + 1, maxsect + 1],dtype='float') npaircollided = np.zeros([maxlength + 1, maxsect + 1],dtype='float') ntilesmax = maskcpqpm.ntiles.max() cpcompgal = np.zeros(ntilesmax + 1, dtype='float') cpcomppair = np.zeros(ntilesmax + 1, dtype='float') myhist = np.zeros(maxlength+1,dtype='int') galsnotcorrected = 0 galscorrected = 0 for k,v in dgroups.iteritems(): myhist[len(v)] += 1 assert len(v) > 1 varr = np.array(v) if 0==0: # and combinedopt == 0: ## do testing, but only correct for non-combined samples. if (dd['IMATCH'][varr] == 0).all(): galsnotcorrected += len(varr) else: galscorrected += len(varr) ## Dec 8 2014, this should also work for combined samples. myigood = np.array([]) for imatchval in fibergood: myigood = np.concatenate((myigood,varr[np.where(dd['IMATCH'][varr] == imatchval)[0]])) for iii in myigood: ngroupfiber[dd['ISECT'][iii]] += 1 myibad = varr[np.where(dd['IMATCH'][varr] == 3)[0]] for iii in myibad: ngroupnofiber[dd['ISECT'][iii]] += 1 ## iterate over all elements in case they're in different sectors. for ii in range(len(myigood)): for jj in range(ii+1,len(myigood)): npairresolved[len(v), dd['ISECT'][myigood[ii]]] += 0.5 npairresolved[len(v), dd['ISECT'][myigood[jj]]] += 0.5 for ii in range(len(myibad)): for jj in range(ii+1,len(myibad)): npaircollided[len(v), dd['ISECT'][myibad[ii]]] += 0.5 npaircollided[len(v), dd['ISECT'][myibad[jj]]] += 0.5 for ii in range(len(myigood)): for jj in range(len(myibad)): npaircollided[len(v), dd['ISECT'][myigood[ii]]] += 0.5 npaircollided[len(v), dd['ISECT'][myibad[jj]]] += 0.5 print 'total gals, total gals in groups w fibers, total gals in groups wofibers, frac:', print len(dd), ngroupfiber.sum(), ngroupnofiber.sum(), ngroupnofiber.sum()/float(ngroupfiber.sum() + ngroupnofiber.sum()+1.e-10) print 'total pairs resolved and collided, frac', npairresolved[:,:].sum(), npaircollided[:,:].sum(), float(npairresolved[:,:].sum())/float(npairresolved[:,:].sum() + npaircollided[:,:].sum() + 1.e-10) if(galsnotcorrected+galscorrected) > 0: print 'fraction uncorrected:',galsnotcorrected, galscorrected, float(galsnotcorrected)/float(galsnotcorrected+galscorrected) # else: # print 'fraction uncorrected:',galsnotcorrected, galscorrected, 0.0 print 'histogram of group sizes' s = 0 for hi in range(len(myhist)): s += myhist[hi] * hi for hi in range(len(myhist)): print hi, myhist[hi], myhist[hi]*hi, myhist[hi]*hi/float(s) ntilesmaxint = int(np.floor(ntilesmax + 0.5)) resultstable = np.zeros([ntilesmaxint+1,6]) for ntilesval in range(1,ntilesmaxint + 1): xxi = np.where(binarymask.ntiles == ntilesval)[0] i2 = ngroupfiber[xxi].sum() i1 = ngroupnofiber[xxi].sum() if i1 + i2 == 0: assert npairresolved[:,xxi].sum() == 0 assert npaircollided[:,xxi].sum() == 0 cpcompgal[ntilesval] = -1. cpcomppair[ntilesval] = -1. else: cpcompgal[ntilesval] = i1/float(i1+i2) cpcomppair[ntilesval] = float(npairresolved[:,xxi].sum())/float(npairresolved[:,xxi].sum() + npaircollided[:,xxi].sum()+1e-10) resultstable[ntilesval,0] = cpcompgal[ntilesval] resultstable[ntilesval,1] = cpcomppair[ntilesval] print 'cpcompgal results for ntiles = ',ntilesval print len(xxi), i2, i1, cpcompgal[ntilesval], cpcomppair[ntilesval] print 'here is the distribution of pair completeness as a function of group size' if maxlength >= 4: for ngroup in range(2,6): if len(xxi) == 0: continue resultstable[ntilesval,ngroup+0] = float(npairresolved[ngroup,xxi].sum())/float(npairresolved[ngroup,xxi].sum()+ npaircollided[ngroup,xxi].sum() + 1.e-10) print ngroup, npairresolved[ngroup,xxi].sum(), npaircollided[ngroup,xxi].sum(), float(npairresolved[ngroup,xxi].sum())/float(npairresolved[ngroup,xxi].sum()+ npaircollided[ngroup,xxi].sum() + 1.e-10) ## finally, fill into mask. for i in range(len(maskcpqpm.mf['WEIGHT'])): ## save pair completeness as well. maskcppair.mf['WEIGHT'][i] = cpcomppair[maskcppair.mf['NTILES'][i]] ### For QPM mocks, just a fraction of galaxies to throw out. maskcpqpm.mf['WEIGHT'][i] = cpcompgal[maskcpqpm.mf['NTILES'][i]] for s in range(len(maskcpqpm.area)): xxq = np.where(maskcpqpm.mf['SECTOR'] == s)[0] assert (maskcpqpm.mf['WEIGHT'][xxq] == maskcpqpm.mf['WEIGHT'][xxq[0]]).all() assert (maskcppair.mf['WEIGHT'][xxq] == maskcppair.mf['WEIGHT'][xxq[0]]).all() maskcpqpm.weight = maskcpqpm.mf['WEIGHT'][xxq[0]] maskcppair.weight = maskcppair.mf['WEIGHT'][xxq[0]] return myhist, resultstable ## put these in a table in the paper. def calc_cp_comp(dd, binarymaskfits, cpmaskqpmfits = None, cpmaskpairfits = None, cp_comp_opt = 0, fibergood=[1,4,7],combinedopt=0): """ Loop over sectors and compute BOSS close pair completeness. make_mask.pro does this by counting number of good pairs over total pairs in the sector. in make_mask.pro cp_comp forced to 0 if NTILES = 1. In this new code, we will loop over fiber collision groups with at least 2 CMASS targets. We will record how many objects in such groups, and what fraction got fibers. mask_cp is the fits file we will fill in with relevant probabilities. It contains NTILES inherited from boss_geom. if combinedopt = 1, remove some asserts that fail for CMASS-LOWZ combined. """ binarymask = maskutils.mask(binarymaskfits) maskcpqpm = maskutils.mask(binarymaskfits) maskcppair = maskutils.mask(binarymaskfits) dgroups = {} maxlength = 0 for di in range(len(dd)): if dd['INGROUP'][di] == -1: continue ## only retain galaxies that are inside the binary mask. if binarymask.weight[dd['ISECT'][di]] < 0.01: continue mykey = groupkey(dd[di]) if dd['IMATCH'][di] == 2: continue ## don't let 2's in groups, new May 25. if mykey in dgroups: dgroups[mykey].append(di) if len(dgroups[mykey]) > maxlength: maxlength = len(dgroups[mykey]) else: dgroups[mykey] = [di] ##trim dictionary to include only groups with at least two CMASS entries. print 'this many items in dict:',len(dgroups) dgroups = {key: value for key, value in dgroups.items() if len(value) > 1} print 'after purge, this many items in dict:',len(dgroups) maxsect = maskcpqpm.mf['SECTOR'].max() ngroupfiber = np.zeros(maxsect + 1,dtype='int') ngroupnofiber = np.zeros(maxsect + 1,dtype='int') ## pair stats. ## what we computed in fiberc_zfail.pro is really confusing, because the imatch values are changing from 0 to 3 as you go. ## we will compute the total number of pairs for which the fiber collision is resolved (both have imatch in imatchgood) and the number that are collided (have at least one imatch = 3). Add 0.5 to each objects sector (so this will count as 1 if they're both in the same sector, but no reason these won't cross sector boundaries, i think) ## we will compute total number of pairs that are collided. I will do this as a function of group size and average over all pairs, see if there's a big diff. npairresolved = np.zeros([maxlength + 1, maxsect + 1],dtype='float') npaircollided = np.zeros([maxlength + 1, maxsect + 1],dtype='float') ntilesmax = maskcpqpm.ntiles.max() cpcompgal = np.zeros(ntilesmax + 1, dtype='float') cpcomppair = np.zeros(ntilesmax + 1, dtype='float') myhist = np.zeros(maxlength+1,dtype='int') galsnotcorrected = 0 galscorrected = 0 for k,v in dgroups.iteritems(): myhist[len(v)] += 1 assert len(v) > 1 varr = np.array(v) if 0==0 and combinedopt == 0: ## do testing, but only correct for non-combined samples. assert len(varr) > 1 assert (dd['MULTGROUP'][varr] == dd['MULTGROUP'][varr[0]]).all() assert (dd['INGROUP'][varr] == dd['INGROUP'][varr[0]]).all() assert len(varr) <= dd['MULTGROUP'][varr[0]] #print dd['IMATCH'][varr] if not ((dd['IMATCH'][varr] > 0).all() or (dd['IMATCH'][varr] == 0).all()): print 'preview assert fail, turning assert off lowz failure group imatch = [0,2]!' print k print v print dd['IMATCH'][varr] assert (dd['IMATCH'][varr] > 0).all() or (dd['IMATCH'][varr] == 0).all() ## should have already been fiber corrected! if (dd['IMATCH'][varr] == 0).all(): galsnotcorrected += len(varr) else: galscorrected += len(varr) ## Dec 8 2014, this should also work for combined samples. myigood = np.array([]) for imatchval in fibergood: myigood = np.concatenate((myigood,varr[np.where(dd['IMATCH'][varr] == imatchval)[0]])) for iii in myigood: ngroupfiber[dd['ISECT'][iii]] += 1 myibad = varr[np.where(dd['IMATCH'][varr] == 3)[0]] for iii in myibad: ngroupnofiber[dd['ISECT'][iii]] += 1 ## iterate over all elements in case they're in different sectors. for ii in range(len(myigood)): for jj in range(ii+1,len(myigood)): npairresolved[len(v), dd['ISECT'][myigood[ii]]] += 0.5 npairresolved[len(v), dd['ISECT'][myigood[jj]]] += 0.5 for ii in range(len(myibad)): for jj in range(ii+1,len(myibad)): npaircollided[len(v), dd['ISECT'][myibad[ii]]] += 0.5 npaircollided[len(v), dd['ISECT'][myibad[jj]]] += 0.5 for ii in range(len(myigood)): for jj in range(len(myibad)): npaircollided[len(v), dd['ISECT'][myigood[ii]]] += 0.5 npaircollided[len(v), dd['ISECT'][myibad[jj]]] += 0.5 print 'total gals, total gals in groups w fibers, total gals in groups wofibers, frac:', print len(dd), ngroupfiber.sum(), ngroupnofiber.sum(), ngroupnofiber.sum()/float(ngroupfiber.sum() + ngroupnofiber.sum()+1.e-10) print 'total pairs resolved and collided, frac', npairresolved[:,:].sum(), npaircollided[:,:].sum(), float(npairresolved[:,:].sum())/float(npairresolved[:,:].sum() + npaircollided[:,:].sum() + 1.e-10) if(galsnotcorrected+galscorrected) > 0: print 'fraction uncorrected:',galsnotcorrected, galscorrected, float(galsnotcorrected)/float(galsnotcorrected+galscorrected) # else: # print 'fraction uncorrected:',galsnotcorrected, galscorrected, 0.0 print 'histogram of group sizes' s = 0 for hi in range(len(myhist)): s += myhist[hi] * hi for hi in range(len(myhist)): print hi, myhist[hi], myhist[hi]*hi, myhist[hi]*hi/float(s) ntilesmaxint = int(np.floor(ntilesmax + 0.5)) resultstable = np.zeros([ntilesmaxint+1,6]) for ntilesval in range(1,ntilesmaxint + 1): xxi = np.where(binarymask.ntiles == ntilesval)[0] i2 = ngroupfiber[xxi].sum() i1 = ngroupnofiber[xxi].sum() if i1 + i2 == 0: assert npairresolved[:,xxi].sum() == 0 assert npaircollided[:,xxi].sum() == 0 cpcompgal[ntilesval] = -1. cpcomppair[ntilesval] = -1. else: cpcompgal[ntilesval] = i1/float(i1+i2) cpcomppair[ntilesval] = float(npairresolved[:,xxi].sum())/float(npairresolved[:,xxi].sum() + npaircollided[:,xxi].sum()+1e-10) resultstable[ntilesval,0] = cpcompgal[ntilesval] resultstable[ntilesval,1] = cpcomppair[ntilesval] print 'cpcompgal results for ntiles = ',ntilesval print len(xxi), i2, i1, cpcompgal[ntilesval], cpcomppair[ntilesval] print 'here is the distribution of pair completeness as a function of group size' if maxlength >= 4: for ngroup in range(2,6): if len(xxi) == 0: continue resultstable[ntilesval,ngroup+0] = float(npairresolved[ngroup,xxi].sum())/float(npairresolved[ngroup,xxi].sum()+ npaircollided[ngroup,xxi].sum() + 1.e-10) print ngroup, npairresolved[ngroup,xxi].sum(), npaircollided[ngroup,xxi].sum(), float(npairresolved[ngroup,xxi].sum())/float(npairresolved[ngroup,xxi].sum()+ npaircollided[ngroup,xxi].sum() + 1.e-10) ## finally, fill into mask. for i in range(len(maskcpqpm.mf['WEIGHT'])): ## save pair completeness as well. maskcppair.mf['WEIGHT'][i] = cpcomppair[maskcppair.mf['NTILES'][i]] ### For QPM mocks, just a fraction of galaxies to throw out. maskcpqpm.mf['WEIGHT'][i] = cpcompgal[maskcpqpm.mf['NTILES'][i]] for s in range(len(maskcpqpm.area)): xxq = np.where(maskcpqpm.mf['SECTOR'] == s)[0] assert (maskcpqpm.mf['WEIGHT'][xxq] == maskcpqpm.mf['WEIGHT'][xxq[0]]).all() assert (maskcppair.mf['WEIGHT'][xxq] == maskcppair.mf['WEIGHT'][xxq[0]]).all() maskcpqpm.weight = maskcpqpm.mf['WEIGHT'][xxq[0]] maskcppair.weight = maskcppair.mf['WEIGHT'][xxq[0]] #def calc_cp_comp(dd, binarymaskfits, cpmaskqpmfits = None, cpmaskpairfits = None, cp_comp_opt = 0, fibergood=[1,4,7]): if cpmaskqpmfits is not None: fitsio.write(cpmaskqpmfits,maskcpqpm.mf,clobber=True) if cpmaskpairfits is not None: fitsio.write(cpmaskpairfits,maskcppair.mf,clobber=True) return myhist, resultstable ## put these in a table in the paper. def FBstatsdr12(ddc, ddl, binarymaskfitslowzcut): binarymask = maskutils.mask(binarymaskfitslowzcut) ## want to know total targets of each class within the binary mask. ntotC = len(np.where(binarymask.weight[ddc['ISECT']] > 0.01)[0]) ntotL = len(np.where(binarymask.weight[ddl['ISECT']] > 0.01)[0]) print 'CMASS before/after binary mask',ntotC, len(ddc) print 'LOWZ before/after binary mask',ntotL, len(ddl) print 'histogram of imatch values within lowzchunkcut mask' xx = np.where(binarymask.weight[ddc['ISECT']] > 0.01)[0] h,x = np.histogram(ddc['IMATCH'][xx],bins=np.arange(-0.5,10.5,1.0)) print 'CMASS raw histogram:',h xx = np.where(binarymask.weight[ddl['ISECT']] > 0.01)[0] h,x = np.histogram(ddl['IMATCH'][xx],bins=np.arange(-0.5,10.5,1.0)) print 'LOWZ raw histogram:',h dd = np.concatenate((ddc,ddl)) ## these keep their order, so last cmass indx is len(ddc)-1; useful later cmassmax = len(ddc) FBrad = 62./3600. ## histogram of all group sizes. ## now make dictionary lookup for groups. dgroups = {} maxlength = 0 for di in range(len(dd)): if dd['INGROUP'][di] == -1: continue if dd['IMATCH'][di] == 2: continue ## don't let 2's in groups, new May 25. ## only retain galaxies that are inside the binary mask. if binarymask.weight[dd['ISECT'][di]] < 0.01: continue mykey = groupkey(dd[di]) if mykey in dgroups: ## Jan 27 2015 ## only append if it's unique! there is CMASS/LOWZ targeting overlap at the boundary. ## first check that pkeynew = photokey(dd[di]) dupflag = 0 for jj in dgroups[mykey]: pkeyold = photokey(dd[jj]) if pkeyold == pkeynew: ## yes, these are duplicates. # print 'duplicate?',dd['RA'][di], dd['RA'][jj], dd['DEC'][di], dd['DEC'][jj] dupflag = 1 if dupflag == 0: dgroups[mykey].append(di) if len(dgroups[mykey]) > maxlength: maxlength = len(dgroups[mykey]) maxkey = mykey else: dgroups[mykey] = [di] ##trim dictionary to include only groups with at least two CMASS entries. print 'this many items in dict:',len(dgroups) dgroups = {key: np.array(value) for key, value in dgroups.items() if len(value) > 1} print 'after purge, this many items in dict:',len(dgroups) print 'max group size',maxlength ## tmp, what's the composition of the huge group? groupindx = dgroups[maxkey] ci = np.where(groupindx < cmassmax)[0] li = np.where(groupindx >= cmassmax)[0] print len(ci) print len(li) print maxkey return -1 nall = np.zeros(maxlength+1) nCall = np.zeros(maxlength+1) nLall = np.zeros(maxlength+1) ## now do some stats on the groups. lowziso = 0 cmassiso = 0 lowzisonoz = 0 cmassisonoz = 0 cmasslostall = 0 lowzlostall = 0 cmassctot = 0 lowzctot = 0 ## tmp, what's going on? ## return dgroups himatch = np.zeros(10) himatchC = np.zeros(10) himatchL = np.zeros(10) for kk, groupindx in dgroups.items(): ci = np.where(groupindx < cmassmax)[0] li = np.where(groupindx >= cmassmax)[0] if len(li) == 1: lowziso += 1 if len(ci) == 1: cmassiso += 1 if len(li) == 1 and dd['IMATCH'][groupindx[li[0]]] == 0: lowzisonoz += 1 if len(ci) == 1 and dd['IMATCH'][groupindx[ci[0]]] == 0: cmassisonoz += 1 if len(ci) > 1 and (dd['IMATCH'][groupindx[ci]] == 0).all(): cmasslostall += len(ci) if len(li) > 1 and (dd['IMATCH'][groupindx[li]] == 0).all(): lowzlostall += len(li) for tmpo in dd['IMATCH'][groupindx]: himatch[tmpo] += 1 for tmpo in dd['IMATCH'][groupindx[ci]]: himatchC[tmpo] += 1 for tmpo in dd['IMATCH'][groupindx[li]]: himatchL[tmpo] += 1 cmassctot += len(ci) lowzctot += len(li) nall[len(groupindx)] += 1 nCall[len(ci)] += 1 nLall[len(li)] += 1 print 'LOWZ numbers',lowziso,lowzisonoz,lowzlostall,lowzctot,ntotL print 'CMASS numbers',cmassiso,cmassisonoz,cmasslostall,cmassctot,ntotC print 'group size histograms' for i in range(maxlength+1): print i, nall[i], nCall[i], nLall[i] print 'histogram of imatch',himatch print 'histogram of imatch CMASS',himatchC print 'histogram of imatch LOWZ',himatchL dcoll = fitsio.read("inputfiles/bosstile-final-collated-boss2-boss38.fits") ## now scroll through things without a redshift. what do the stats look like? xx = np.where((binarymask.weight[ddc['ISECT']] > 0.01) & (ddc['IMATCH'] == 0))[0] cnt = 0 for xi in xx: mykey = groupkey(ddc[xi]) if mykey in dgroups: continue ww = np.where((dcoll['FINALN'] == ddc['FINALN'][xi]) & (dcoll['INGROUP'] == ddc['INGROUP'][xi]))[0] print mykey, ddc['MULTGROUP'][xi] print 'target file chk', dcoll['SOURCETYPE'][ww], dcoll['BOSS_TARGET1_009'][ww], dcoll['BOSS_TARGET1_009'][ww] & 2**1 == 2**1, dcoll['BOSS_TARGET1_009'][ww] & 2**0 == 2**0 cnt += 1 print 'there were ',cnt, 'CMASS weirdos' xx = np.where((binarymask.weight[ddl['ISECT']] > 0.01) & (ddl['IMATCH'] == 0))[0] cnt = 0 for xi in xx: mykey = groupkey(ddl[xi]) if mykey in dgroups: continue ww = np.where((dcoll['FINALN'] == ddl['FINALN'][xi]) & (dcoll['INGROUP'] == ddl['INGROUP'][xi]))[0] print mykey, ddl['MULTGROUP'][xi] print 'target file chk', dcoll['SOURCETYPE'][ww], dcoll['BOSS_TARGET1_009'][ww], dcoll['BOSS_TARGET1_009'][ww] & 2**1 == 2**1, dcoll['BOSS_TARGET1_009'][ww] & 2**0 == 2**0 cnt += 1 print 'there were ',cnt, 'LOWZ weirdos' print 'possibilities for hte weirdos..' print 'I didnt find their match in the spAll, but they are there.' print 'they are in chunk overlap region and the collision group stuff gets confusing, though that shouldnt leave un-fibered objects without a collision group, right?' print 'examining the printed groups, the collided objects are in the collate file and their sourcetype is LRG, but their BOSS_TARGET_009 was 0 for the two groups i checked.' #dd['IMATCH'][cmassi] == 1 def FBcmasslowzchk(ddc, ddl,binarymaskfitslowzcut,outfname,multcut=True): """ Input the output of match_sample.py for both cmass and lowz (do N and S separately). Look for CMASS-lowz pairs, plot their redshift distribution. """ dd = np.concatenate((ddc,ddl)) ## these keep their order, so last cmass indx is len(ddc)-1; useful later cmassmax = len(ddc) FBrad = 62./3600. binarymask = maskutils.mask(binarymaskfitslowzcut) ## now make dictionary lookup for groups. dgroups = {} maxlength = 0 for di in range(len(dd)): if dd['INGROUP'][di] == -1: continue if dd['IMATCH'][di] == 2: continue ## don't let 2's in groups, new May 25. ## only retain galaxies that are inside the binary mask. if binarymask.weight[dd['ISECT'][di]] < 0.01: continue mykey = groupkey(dd[di]) if mykey in dgroups: ## Jan 27 2015 ## only append if it's unique! there is CMASS/LOWZ targeting overlap at the boundary. ## first check that pkeynew = photokey(dd[di]) dupflag = 0 for jj in dgroups[mykey]: pkeyold = photokey(dd[jj]) if pkeyold == pkeynew: ## yes, these are duplicates. # print 'duplicate?',dd['RA'][di], dd['RA'][jj], dd['DEC'][di], dd['DEC'][jj] dupflag = 1 if dupflag == 0: dgroups[mykey].append(di) if len(dgroups[mykey]) > maxlength: maxlength = len(dgroups[mykey]) else: dgroups[mykey] = [di] ##trim dictionary to include only groups with at least two CMASS entries. print 'this many items in dict:',len(dgroups) dgroups = {key: value for key, value in dgroups.items() if len(value) > 1} print 'after purge, this many items in dict:',len(dgroups) zlistcmass1 = [] zlistlowz1 = [] zlistcmass2 = [] ## both objects got a riber. zlistlowz2 = [] ## both objects got a riber. ofpp = open(outfname+'.zpairs','w') cntack = 0 cntacktot = 0 mhist = np.zeros(20) ## now iterate and look for CMASS/LOWZ pairs. for kk, groupindx in dgroups.items(): ## upgrade; we want to consider groups of all sizes (???) ## complicated because we already fiber correct groups that have a multiple of the same target. ## ignore for now. if len(groupindx) != 2: continue ## MULTGROUP multiplicity. mhist[dd['MULTGROUP'][groupindx[0]]] += 1 # if dd['MULTGROUP'][groupindx[0]] != 2: # print 'whyyyyyyyyyy',dd['MULTGROUP'][np.array(groupindx)], dd['FINALN'][np.array(groupindx)], dd['INGROUP'][np.array(groupindx)] if multcut is True: if dd['MULTGROUP'][groupindx[0]] != 2: continue assert dd['FINALN'][groupindx[0]] == dd['FINALN'][groupindx[1]] assert dd['INGROUP'][groupindx[0]] == dd['INGROUP'][groupindx[1]] cmassi = -1 if (groupindx[0] < cmassmax and groupindx[1] >= cmassmax): cmassi = groupindx[0] lowzi = groupindx[1] if (groupindx[1] < cmassmax and groupindx[0] >= cmassmax): cmassi = groupindx[1] lowzi = groupindx[0] if cmassi == -1: continue angsep = mksampleutils.angsepdeg(dd['RA'][cmassi],dd['DEC'][cmassi], dd['RA'][lowzi],dd['DEC'][lowzi]) cntacktot += 1 if angsep > FBrad: # print 'mistake??',angsep/FBrad, dd['RA'][cmassi],dd['DEC'][cmassi], dd['RA'][lowzi],dd['DEC'][lowzi], dd['IMATCH'][cmassi], dd['IMATCH'][lowzi] cntack += 1 if dd['IMATCH'][cmassi] == 1 and dd['IMATCH'][lowzi] == 0: zlistcmass1.append(dd['Z'][cmassi]) assert dd['MULTGROUP'][cmassi] >= 2 if dd['IMATCH'][lowzi] == 1 and dd['IMATCH'][cmassi] == 0: zlistlowz1.append(dd['Z'][lowzi]) assert dd['MULTGROUP'][lowzi] >= 2 if dd['IMATCH'][cmassi] == 1 and dd['IMATCH'][lowzi] == 1: zlistcmass2.append(dd['Z'][cmassi]) zlistlowz2.append(dd['Z'][lowzi]) ofpp.write('%e %e\n' % (dd['Z'][cmassi], dd['Z'][lowzi])) ofpp.close() print 'here are len of lists',len(zlistcmass1), len(zlistcmass2), len(zlistlowz1), len(zlistlowz2) hedge = np.arange(0.,1.006,0.005) hcen = 0.5*(hedge[:-1] + hedge[1:]) hcmass1, xtmp = np.histogram(np.array(zlistcmass1),bins=hedge) hlowz1, xtmp = np.histogram(np.array(zlistlowz1),bins=hedge) hcmass2, xtmp = np.histogram(np.array(zlistcmass2),bins=hedge) hlowz2, xtmp = np.histogram(np.array(zlistlowz2),bins=hedge) hcmassall, xtmp = np.histogram(dd['Z'][:cmassmax],bins=hedge) hlowzall, xtmp = np.histogram(dd['Z'][cmassmax:],bins=hedge) ofp = open(outfname,'w') for i in range(len(hcen)): ofp.write('%f %d %d %d %d %d %d\n' % (hcen[i],hcmass1[i],hlowz1[i], hcmass2[i], hlowz2[i],hcmassall[i],hlowzall[i])) ofp.close() print 'most ang seps are below the FB scale',cntack,cntacktot,float(cntack)/float(cntacktot) print 'mhist' for i in range(20): print mhist[i] def fibercollide(dd, binarymaskfits, FBopt=0, legacyopt = 0, fbmatchgood=[1,4,7],fbmatchbad=[0,2,3,5,6]): binarymask = maskutils.mask(binarymaskfits) ## make sure we're doing fiber collisions before redshift failures. assert (dd['IMATCH'] != 5).all() print 'only fiber correcting galaxies inside the binary mask!!' xx = np.where((dd['ICOLLIDED'] == 1) & (dd['IMATCH'] != 0) & (binarymask.weight[dd['ISECT']] > 0.01))[0] print 'hi beth, this many galaxies are collided but imatch != 0',len(xx) #print 'new on May 16, mark those as imatch = 0 so they will be fiber corrected.' #print 'this overrides ancillary target programs that added redshifts for fiber collisions.' #dd['IMATCH'][xx] = 0 ## nevermind, let's not do that! ## see email to boss-gal-clustering regarding spAll SPECPRIMARY and ancillary targets, May 16. xx = np.where((dd['ICOLLIDED'] == 1) & (dd['IMATCH'] == 0) & (binarymask.weight[dd['ISECT']] > 0.01))[0] print 'hi beth, this many galaxies are collided but imatch == 0',len(xx) fiberloop(dd,xx,FBopt=FBopt, legacyopt = legacyopt, fbmatchgood=fbmatchgood, fbmatchbad=fbmatchbad) yy = np.where((dd['ICOLLIDED'] == 0) & (dd['IMATCH'] == 0) & (binarymask.weight[dd['ISECT']] > 0.01))[0] print 'starting second set, labelled with ICOLLIDED = 0' fiberloop(dd,yy,FBopt=FBopt, legacyopt = legacyopt, fbmatchgood=fbmatchgood, fbmatchbad=fbmatchbad) ## ttt = np.where((dd['WEIGHT_CP'] > 1) & (dd['IMATCH'] == 1))[0] s1 = (dd['WEIGHT_CP'][ttt] - 1.0).sum() ttt = np.where((dd['WEIGHT_CP'] > 1) & (dd['IMATCH'] == 2))[0] s2 = (dd['WEIGHT_CP'][ttt] - 1.0).sum() ttt = np.where((dd['WEIGHT_CP'] > 1) & (dd['IMATCH'] == 4))[0] s4 = (dd['WEIGHT_CP'][ttt] - 1.0).sum() ttt = np.where((dd['WEIGHT_CP'] > 1) & (dd['IMATCH'] == 7))[0] s7 = (dd['WEIGHT_CP'][ttt] - 1.0).sum() print 'fb weights distribution on imatch 1,2,4,7:',s1,s2,s4,s7 def fbsanity(): """ Spherematch at 62 arcsec and see if we have any surprises as a sanity check. """ pass def ifib2reweight(dd): ifibmin = 16. Nbins = 100 binsize = 0.1 ifib2bins = np.arange(ifibmin,ifibmin+(Nbins+0.5)*binsize,binsize) ifib2mag = 22.5 - 2.5*np.log10(dd['FIBER2FLUX'][:,3]) - dd['EXTINCTION'][:,3] assert (dd['FIBER2FLUX'][:,3] > 0.).all() ## indices for redshift failures, good redshifts, and upweighted galaxies. noz = np.where(dd['IMATCH'] == 5)[0] up = np.where(((dd['IMATCH'] == 1) | (dd['IMATCH'] == 4)) & (dd['WEIGHT_NOZ'] > 1.0001))[0] ## technically hup should be weighted with dd['WEIGHT_NOZ'] since it can be greater than 1, i.e., represent weight from more than one galaxy. ## noz should be weighted by WEIGHT_CP since the noz weight gets WEIGHT_CP from the redshift failure. hup, xtmp = np.histogram(ifib2mag[up], bins=ifib2bins,weights=(dd['WEIGHT_NOZ'][up]-1.)) hnoz, xtmp = np.histogram(ifib2mag[noz], bins=ifib2bins,weights=(dd['WEIGHT_CP'][noz])) hupnorm = hup/float(hup.sum()) hnoznorm = hnoz/float(hnoz.sum()) ifib2weight = hnoznorm/(hupnorm) ## set weight = 1 where there are not many galaxies of either type, to avoid crazy shot noise. xx = np.where((hup < 10) | (hnoz < 10))[0] ifib2weight[xx] = 1. ## this seems to be fine; neither distribution has many galaxies in this region of the histograms (0.6% each) for dr12v3-S. xcen = 0.5*(xtmp[1:] + xtmp[:-1]) assert (dd['WEIGHT_NOZ'][up] > 1.0001).all() weight_new = (dd['WEIGHT_NOZ'][up] - 1.0)*np.interp(ifib2mag[up],xcen,ifib2weight) weight_new = weight_new/weight_new.sum() * (dd['WEIGHT_NOZ'][up] - 1.0).sum() tmp = (dd['WEIGHT_NOZ'][up] - 1.0).sum() tmp2 = weight_new.sum() print 'sum of weights, should be equal',tmp,tmp2 assert np.fabs(tmp - tmp2) < 1.0e-3 for i in range(len(dd['WEIGHT_NOZ'][up])): dd['WEIGHT_NOZ'][up[i]] = weight_new[i] + 1.0 ## you can't do this because of the mismatch in endianness! ACK!! ##dd['WEIGHT_NOZ'][up] = weight_new + 1.0 print 'srsly',dd['WEIGHT_NOZ'][up][0:10] return xcen, hup, hnoz def wgtzfail(dd, zmatchgood=[1,4],zfailopt=0,ifib2reweightopt=0): """ zmatchgood are galaxies allowed to receive redshift failure weights. zfailopt = 0: match sector-by-sector. zfailopt = 1: match plate+mjd-by-plate+mjd ifib2reweight: apply Rita's algorithm to reweight zfails to match fiber magnitudes (1). Or skip (0). Had to make changes compared to the original code to make this work with astropy version of pydl installed on the Utah computers. """ assert zfailopt == 0 or zfailopt == 1 ## only options we have defined. zfail = np.where(dd['IMATCH'] == 7)[0] aseplistz = [] zfnbrtot = 0 zfnbrwgttot = 0 zfnonbrtot = 0 zfraddeg = 10. ## find all indices with imatch values of interest. myi = np.array([],dtype='int') for imatchval in zmatchgood: myi = np.concatenate((myi,np.where(dd['IMATCH'] == imatchval)[0])) print 'this many objects available for redshift failure weights',len(myi) for zi in zfail: if zfailopt == 0: mgal = np.intersect1d(np.where(dd['ISECT'] == dd['ISECT'][zi])[0], myi) else: mgal = np.intersect1d(np.where((dd['PLATE'] == dd['PLATE'][zi]) & (dd['MJD'] == dd['MJD'][zi]))[0], myi) ## ugh! new version of pydl fails when len(mgal) == 1. Should email Ben. ## this is a hack to make it work. if len(mgal) == 1: mgal = np.concatenate((mgal,mgal)) ## error: ## raise PydlutilsException("Change the order of the sets of coordinates!") ## pydl.pydlutils.PydlutilsException: Change the order of the sets of coordinates! #ddi, ddzfaili, dist = pydlutils.spheregroup.spherematch(np.array([dd['RA'][zi]]), np.array([dd['DEC'][zi]]), dd['RA'][mgal], dd['DEC'][mgal], matchlength = zfraddeg, maxmatch=len(mgal)+1) if len(mgal) > 0: ## fails for len(mgal) == 0 case. ddzfaili, ddi, dist = pydlutils.spheregroup.spherematch(dd['RA'][mgal], dd['DEC'][mgal], np.array([dd['RA'][zi]]), np.array([dd['DEC'][zi]]), matchlength = zfraddeg, maxmatch=len(mgal)+1) else: continue if len(ddi) > 0: ## remove self. assert (ddi == 0).all() assert zi not in mgal[ddzfaili] zfailnbr = mgal[ddzfaili[np.where(dist == dist.min())[0]]] ## more changes: because of the duplication hack, now need to select hte first elt of zfailnbr. assert dd['IMATCH'][zfailnbr[0]] in zmatchgood aseplistz.append(dist.min()) dd['IMATCH'][zi] = 5 dd['Z'][zi] = dd['Z'][zfailnbr[0]] dd['WEIGHT_NOZ'][zfailnbr[0]] += dd['WEIGHT_CP'][zi] zfnbrwgttot += dd['WEIGHT_CP'][zi] zfnbrtot += 1 else: print 'found no neighbor for ',dd['RA'][zi], dd['DEC'][zi] ## algorithm written by Rita. if ifib2reweightopt == 1: ifib2reweight(dd) print 'zfail stats:',zfnbrtot, zfnbrwgttot ## histogram aseplistz aseplistz = np.array(aseplistz) bins = np.arange(0.,aseplistz.max(),1./3600.) h,x = np.histogram(aseplistz,bins=bins) xcen = 0.5*(x[1:] + x[:-1]) return h,xcen def wgtzfailtest(dd, zmatchgood=[1,4]): """ zmatchgood are galaxies allowed to receive redshift failure weights. zfailopt = 0: match sector-by-sector. zfailopt = 1: match plate-by-plate This function compares histogram of nearest-neighbor angular separations between same-sector and same-plate/mjd choices (zfailopt = 0 vs zfailopt = 1). Needs to be updated with all the pydl changes above if I want it to work in Utah! """ zfail = np.where(dd['IMATCH'] == 7)[0] aseplistz1 = [] aseplistz2 = [] zfnbrtot = 0 zfnbrwgttot = 0 zfnonbrtot = 0 zfraddeg = 10. ## find all indices with imatch values of interest. myi = np.array([],dtype='int') for imatchval in zmatchgood: myi = np.concatenate((myi,np.where(dd['IMATCH'] == imatchval)[0])) print 'this many objects available for redshift failure weights',len(myi) c1 = np.zeros(3) cnt = 0 for zi in zfail: cnt += 1 mgal1 = np.intersect1d(np.where(dd['ISECT'] == dd['ISECT'][zi])[0], myi) mgal2 = np.intersect1d(np.where((dd['PLATE'] == dd['PLATE'][zi]) & (dd['MJD'] == dd['MJD'][zi]))[0], myi) for mgal,aseplistz in zip([mgal1, mgal2],[aseplistz1,aseplistz2]): ddi, ddzfaili, dist = pydlutils.spheregroup.spherematch(np.array([dd['RA'][zi]]), np.array([dd['DEC'][zi]]), dd['RA'][mgal], dd['DEC'][mgal], matchlength = zfraddeg, maxmatch=len(mgal)+1) if len(ddi) > 0: ## remove self. assert (ddi == 0).all() assert zi not in mgal[ddzfaili] zfailnbr = mgal[ddzfaili[np.where(dist == dist.min())[0]]] assert dd['IMATCH'][zfailnbr] in zmatchgood aseplistz.append(dist.min()) #print zfailnbr,dist.min(), else: print 'found no neighbor for ',dd['RA'][zi], dd['DEC'][zi] #print ' ' if aseplistz1[-1] == aseplistz2[-1]: c1[0] += 1 elif aseplistz1[-1] < aseplistz2[-1]: c1[1] += 1 else: c1[2] += 1 # print 'zfail stats:',zfnbrtot, zfnbrwgttot ## histogram aseplistz aseplistz1 = np.array(aseplistz1) aseplistz2 = np.array(aseplistz2) bins = np.arange(0.,max(aseplistz1.max(),aseplistz2.max()), 1./3600.) h1,x = np.histogram(aseplistz1,bins=bins) h2,x = np.histogram(aseplistz2,bins=bins) xcen = 0.5*(x[1:] + x[:-1]) assert len(h1) == len(h2) for i in range(len(h1)): print xcen[i],h1[i],h2[i] return 0 def makecptables(hlist, rlist,t1fname,t2fname, NSweights): try: ofp1 = open(t1fname,'w') ofp2 = open(t2fname,'w') except: return -1 ## histogram of group sizes. ofp1.write(r"$n_{\rm group}$ & $f_{\rm north}$ & $f_{\rm south}$" + "\n") maxl = max(len(hlist[0]), len(hlist[1])) snorth = 0 ssouth = 0 for i in range(len(hlist[0])): snorth += hlist[0][i]*i for i in range(len(hlist[1])): ssouth += hlist[1][i]*i snorth = float(snorth) ssouth = float(ssouth) for i in range(maxl): if i==0 or i==1: continue try: nval = hlist[0][i]*i/snorth except: nval = 0 try: sval = hlist[1][i]*i/ssouth except: sval = 0 ofp1.write('%d & %.4f & %.4f\\\\\n' % (i,nval,sval)) ofp1.close() ## completeness as a function of ntile and groupsize. ofp2.write(r'$n_{\rm tiles}$ & $f_{\rm fiber}$ & $c_{\rm pair}$ & $c_{\rm pair}(2)$ & $c_{\rm pair}(3)$ & $c_{\rm pair}(4)$ & $c_{\rm pair}(5)$' + '\n') for i in range(1,len(rlist[0][:,0])): ofp2.write('%d ' % (i)) for j in range(len(rlist[0][0,:])): val = NSweights[0] * rlist[0][i,j] + NSweights[1] * rlist[1][i,j] if j==0: val = 1-val ofp2.write('& %.3f ' % (val)) ofp2.write('\\\\\n') ofp2.close() if __name__ == '__main__': ## standard legacyopt = 0 ## tmp testing!! ## legacyopt = 1 #fibercollide(fitsio.read("../data/match-cmass-dr12v0-N.fits"), maskfits="../data/mask-cmass-dr12v0-N-Anderson.fits",legacyopt = legacyopt) #print 'start S' #fibercollide(fitsio.read("../data/match-cmass-dr12v0-S.fits"), maskfits="../data/mask-cmass-dr12v0-S-Anderson.fits", legacyopt = legacyopt) if 0==1: histlist = [] resultslist = [] NSweight = [] for NS in ['N','S']: dd = fitsio.read("../data/collided-cmass-dr12v3-%s-Reid.fits" % NS) binaryfitsfname = "../data/binarymask-dr12v3-%s.fits" % NS qpmfitsfname = "../data/cpgalqpm-dr12v3-%s.fits" % NS cppairfitsfname = "../data/cppair-dr12v3-%s.fits" % NS mb = maskutils.mask(binaryfitsfname) NSweight.append(len(np.where(mb.weight[dd['ISECT']] > 0.001)[0])) myhist, results = calc_cp_comp(dd, binarymaskfits = binaryfitsfname, \ cpmaskqpmfits = qpmfitsfname, cpmaskpairfits = cppairfitsfname, \ cp_comp_opt = 0, fibergood=[1,4,7]) histlist.append(myhist) resultslist.append(results) NSweight = np.array(NSweight) NSweight = NSweight/float(NSweight.sum()) makecptables(histlist,resultslist,'cmassgrouphist.txt','cmasscpntiles.txt',NSweight) ## join the binary file for final result to send to Martin. maskutils.maskNSjoin("../data/cpgalqpm-dr12v3-N.fits", "../data/cpgalqpm-dr12v3-S.fits", "../data/cpgalqpm-dr12v3.fits") maskutils.maskNSjoin("../data/cpgalqpm-dr12v3-N.fits", "../data/cpgalqpm-dr12v3-S.fits", "../data/cpgalqpm-dr12v3.fits") print 'BR switch to new geom!' sys.exit(1) maskutils.maskfitstoply(maskfitsfname="../data/cpgalqpm-dr12v3.fits",geomfitsfname="/home/bareid/boss/bosslss/trunk/geometry/boss_geometry_2014_04_05.fits") if 0==1: dNcmass = fitsio.read("../data/match-cmass-dr12v3-N.fits") dScmass = fitsio.read("../data/match-cmass-dr12v3-S.fits") dNlowz = fitsio.read("../data/match-lowz-dr12v3-N.fits") dSlowz = fitsio.read("../data/match-lowz-dr12v3-S.fits") binarymaskfitslowzcutN = "../data/binarymask-dr12v3-lowzchunkcut-N.fits" binarymaskfitslowzcutS = "../data/binarymask-dr12v3-lowzchunkcut-S.fits" FBcmasslowzchk(dNcmass, dNlowz, binarymaskfitslowzcutN,"cmasslowzFBpairzdistN.dat") FBcmasslowzchk(dScmass, dSlowz, binarymaskfitslowzcutS,"cmasslowzFBpairzdistS.dat") if 0==1: ## call calc_cp_new! ddcN = fitsio.read("../data/collided-cmass-dr12v4-N-Reid.fits") ddcS = fitsio.read("../data/collided-cmass-dr12v4-S-Reid.fits") ddlN = fitsio.read("../data/collided-lowz-dr12v4-N-Reid.fits") ddlS = fitsio.read("../data/collided-lowz-dr12v4-S-Reid.fits") myhist, results = calc_cp_comp_new(ddcN, ddcS, ddlN, ddlS, "../data/binarymask-dr12v4-lowzchunkcut.fits") ## hack to make cp tables; this function was written to combine N and S. makecptables([myhist,myhist],[results,results],'FBgrouphist-dr12v4.txt','FBgroupscpntiles.txt',np.array([0.5,0.5])) if 0==1: ## finally, compute fraction of area in different tile categories; this will go into first column of table 4. bb = maskutils.mask("../data/binarymask-dr12v4.fits") ss = (180./np.pi)**2 print 'total areas',bb.area.sum()*ss, (bb.area*bb.weight).sum()*ss ntilesmax = max(bb.ntiles) ntilesmaxint = int(np.floor(ntilesmax + 0.5)) at = np.zeros(ntilesmaxint + 1) for ntilesval in range(1,ntilesmaxint + 1): xxi = np.where(bb.ntiles == ntilesval)[0] at[ntilesval] = (bb.area[xxi]*bb.weight[xxi]).sum()*ss print ntilesval, at[ntilesval] print 'now fractions' for ntilesval in range(1,ntilesmaxint + 1): print ntilesval, at[ntilesval]/at.sum() if 0==0: ## call calc_cp_new! ddcN = fitsio.read("../data/collided-cmass-dr12v4-N-Reid.fits") ddcS = fitsio.read("../data/collided-cmass-dr12v4-S-Reid.fits") ddlN = fitsio.read("../data/collided-lowz-dr12v4-N-Reid.fits") ddlS = fitsio.read("../data/collided-lowz-dr12v4-S-Reid.fits") calc_dzhist(ddcN, ddcS, ddlN, ddlS, "../data/binarymask-dr12v4-lowzchunkcut.fits") #def calc_cp_comp_new(ddcN, ddcS, ddlN, ddlS, binarymaskfitslowzcut): #def makecptables(hlist, rlist,t1fname,t2fname, NSweights):