import os import numpy as np import re import sys import fitsio import mksampleutils import time 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']) class mask: def __init__(self,fitsfname): self.fitsfname = fitsfname mf = fitsio.read(fitsfname) (tmp, tlength) = mf['TILES'].shape maxsect = mf['SECTOR'].max() minus1 = np.array([-1]) ## index things by sector rather than polygon. area = np.zeros(maxsect+1,dtype='float') ntiles = np.zeros(maxsect+1,dtype='float') #tiles = np.zeros([maxsect+1,3],dtype='int') ## make this a list of numpy arrays of varying length so we can remove -1 tiles = [] ra_m90sec = np.zeros(maxsect+1) decsec = np.zeros(maxsect+1) weight = np.zeros(maxsect+1) glongpoly, glatpoly = mksampleutils.euler(mf['RAMID'],mf['DECMID']) NorSpoly = np.zeros(len(mf['RAMID']),dtype='int') -1 north = np.where(glatpoly > 0)[0] south = np.where(glatpoly <= 0)[0] NorSpoly[north] = 0 NorSpoly[south] = 1 assert ((NorSpoly == 0) | (NorSpoly == 1)).all() NorS = np.zeros(maxsect+1,dtype='int') ## is each sector N or S? ra_m90poly = np.zeros(len(mf)) ra_m90poly = mf['RAMID']-90.0 tt = np.where(ra_m90poly < 0.0)[0] ra_m90poly[tt] += 360. tt = np.where(ra_m90poly < 0.0)[0] if len(tt) != 0: print 'ramid Error!' sys.exit(1) for s in range(len(area)): xxq = np.where(mf['SECTOR'] == s)[0] assert len(xxq) > 0 assert (NorSpoly[xxq] == NorSpoly[xxq[0]]).all() NorS[s] = NorSpoly[xxq[0]] area[s] = (mf['STR'][xxq]).sum() ra_m90sec[s] = (ra_m90poly[xxq]*mf['STR'][xxq]).sum()/area[s] decsec[s] = (mf['DECMID'][xxq]*mf['STR'][xxq]).sum()/area[s] assert (mf['NTILES'][xxq] == mf['NTILES'][xxq[0]]).all() assert (mf['WEIGHT'][xxq] == mf['WEIGHT'][xxq[0]]).all() ntiles[s] = mf['NTILES'][xxq[0]] weight[s] = mf['WEIGHT'][xxq[0]] for ti in range(tlength): assert (mf['TILES'][xxq[:],ti] == mf['TILES'][xxq[0],ti]).all() # tiles[s,:] = mf['TILES'][xxq[0],:] tiles.append(np.setdiff1d(mf['TILES'][xxq[0],:],minus1)) assert len(np.setdiff1d(mf['TILES'][xxq[0],:],minus1) ) == ntiles[s] self.mf = mf self.area = area self.ntiles = ntiles self.weight = weight self.tiles = tiles self.ra_m90sec = ra_m90sec self.decsec = decsec self.NorS = NorS def maskstats(self,verbose=1): xN = np.where(self.NorS == 0)[0] xS = np.where(self.NorS == 1)[0] aN = self.area[xN].sum() aS = self.area[xS].sum() atot = self.area[:].sum() try: assert np.fabs(aN + aS - atot) < 1.e-5*atot except: 'failed north/south area check, aborting!' sys.exit(1) degperstr = (180./np.pi)**2 max1 = self.mf['SECTOR'].max() ## compute completeness mean and std dev where completeness > 0.0001 xgood = np.where((self.weight > 0.0001))[0] xNgood = np.where((self.weight > 0.0001) & (self.NorS == 0))[0] xSgood = np.where((self.weight > 0.0001) & (self.NorS == 1))[0] aa = np.zeros([3,4]) for ii, xx in zip(range(len([xNgood, xSgood,xgood])), [xNgood, xSgood,xgood]): agood = (self.area[xx]).sum() awgt = (self.area[xx]*self.weight[xx]).sum() cmean = (self.weight[xx]*self.area[xx]).sum()/agood cstd = ((self.weight[xx]**2*self.area[xx]).sum()/agood - cmean**2)**0.5 aa[ii,:] = np.array([agood,awgt,cmean,cstd]) if verbose == 1: print 'mask file has ',max1,'sectors' print 'Area [north, south, tot]: %.4f %.4f %.4f steradians.' % (aN, aS, atot) print 'Area [north, south, tot]: %.1f %.1f %.1f degrees.' % (aN*degperstr, aS*degperstr, atot*degperstr) print 'Good Area [north, south, tot]: %.4f %.4f %.4f steradians.' % (aa[0,0],aa[1,0],aa[2,0]) print 'Good Area [north, south, tot]: %.1f %.1f %.1f degrees.' % (aa[0,0]*degperstr, aa[1,0]*degperstr, aa[2,0]*degperstr) print 'Completeness-weighted Area [north, south, tot]: %.1f %.1f %.1f degrees.' % (aa[0,1]*degperstr, aa[1,1]*degperstr, aa[2,1]*degperstr) print 'Completeness mean [north, south, tot]: %.4f %.4f %.4f' % (aa[0,2],aa[1,2],aa[2,2]) print 'Completeness std [north, south, tot]: %.4f %.4f %.4f' % (aa[0,3],aa[1,3],aa[2,3]) return aa def maskfitstoply(maskfitsfname,geomfitsfname): geomplyfname = geomfitsfname.split('.fits')[0] + '.ply' fitsmaskbase = maskfitsfname.split('.fits')[0] myidlstr = """idl << EOF\n.com maskfitstoply\nmaskfitstoply,"%s","%s"\nEOF\n""" % (fitsmaskbase,geomplyfname) print 'calling idl like this:' print myidlstr os.system(myidlstr) return 0 def runransack(maskfitsfname,geomfitsfname,ransackseed,Nran): plyname = maskfitsfname.split('.fits')[0] + '.ply' try: ifp = open(plyname,'r') ifp.close() except: #ply file doesn't exist, so make it. maskfitstoply(maskfitsfname,geomfitsfname) tmpfile = maskfitsfname.split('.fits')[0] + '.'+str(ransackseed) options = '-c '+str(ransackseed) + ' -r '+str(int(np.floor(Nran))) cmd = os.environ['IDLUTILS_DIR'] + '/bin/ransack ' + options + ' ' + plyname + ' ' + tmpfile print cmd os.system(cmd) ra, dec = np.loadtxt(tmpfile,usecols=[0,1],unpack=True,skiprows=1) ipoly = np.loadtxt(tmpfile,usecols=[2],unpack=True,skiprows=1) return ra, dec, ipoly def maskNSjoin(mNfits,mSfits,mNSfitsout): mN = mask(mNfits) mS = mask(mSfits) mfinal = mask(mNfits) assert (mN.NorS == mS.NorS).all() assert (mN.mf['SECTOR'] == mS.mf['SECTOR']).all() xN = np.where(mN.NorS[mN.mf['SECTOR']] == 0)[0] xS = np.where(mN.NorS[mN.mf['SECTOR']] == 1)[0] mfinal.mf['WEIGHT'][xN] = mN.mf['WEIGHT'][xN] mfinal.mf['WEIGHT'][xS] = mS.mf['WEIGHT'][xS] assert len(xN) + len(xS) == len(mN.mf['WEIGHT']) fitsio.write(mNSfitsout,mfinal.mf,clobber=True) return 0 def makebinary(platelistfits = "/home/bareid/boss/bosslss/trunk/mksample/inputfiles/platelist.fits",\ geomfitsfname = "/home/bareid/boss/bosslss/trunk/geometry/boss_geometry_2014_05_28.fits",\ binarymaskfitsout = "/home/bareid/boss/bosslss/trunk/data/binarymask-dr12v2.fits",\ splitNSopt = 1, lowzchunkcutopt = 0): """ Input the boss geometry and platelist fits files. Output a binarymask of whether all tiles have been observed. if splitNSopt = 1, output separate binary masks for N and S. if lowzchunkcutopt = 1, cut out all sectors with a chunkvector < 7. if lowzchunkcutopt = 2, include chunk 2 and >= 7. if lowzchunkcutopt = 3, include chunks >= 3. """ pp = fitsio.read(platelistfits) goodtilelist = np.unique(pp[np.where((pp['PLATEQUALITY'] == 'good') & \ (pp['PROGRAMNAME'] == 'boss '))[0]]['TILEID']) binarymask = mask(geomfitsfname) if splitNSopt == 1: binarymaskN = mask(geomfitsfname) binarymaskS = mask(geomfitsfname) maxsect = binarymask.mf['SECTOR'].max() alltiles = np.unique(np.concatenate(binarymask.tiles)) assert (alltiles[:-1] < alltiles[1:]).all() ## slow but works? tilecheck = np.zeros(alltiles.max()+1,dtype='int') for tileval in alltiles: if tileval in goodtilelist: ## adopt list from platelist as standard for now. tilecheck[tileval] = 1 else: tilecheck[tileval] = 0 ## now loop over all sectors, see if their tile values are all observed. agood = 0. abad = 0. for s in range(maxsect+1): xxq = np.where(binarymask.mf['SECTOR'] == s)[0] assert len(xxq) > 0 if (tilecheck[binarymask.tiles[s][:]] == 1).all(): binarymask.weight[s] = 1. binarymask.mf['WEIGHT'][xxq] = 1. agood += binarymask.area[s] if splitNSopt == 1: if binarymask.NorS[s] == 0: ## this sector is in the north. binarymaskS.weight[s] = 0. binarymaskS.mf['WEIGHT'][xxq] = 0. binarymaskN.weight[s] = 1. binarymaskN.mf['WEIGHT'][xxq] = 1. else: binarymaskN.weight[s] = 0. binarymaskN.mf['WEIGHT'][xxq] = 0. binarymaskS.weight[s] = 1. binarymaskS.mf['WEIGHT'][xxq] = 1. else: binarymask.weight[s] = 0. binarymask.mf['WEIGHT'][xxq] = 0. abad += binarymask.area[s] if splitNSopt == 1: binarymaskS.weight[s] = 0. binarymaskS.mf['WEIGHT'][xxq] = 0. binarymaskN.weight[s] = 0. binarymaskN.mf['WEIGHT'][xxq] = 0. print 'sector ',s,'is missing a tile of area ', binarymask.area[s], binarymask.tiles[s], tilecheck[binarymask.tiles[s][:]] if lowzchunkcutopt == 1: xxchunkcut = np.where(((binarymask.mf['CHUNKVECTOR'][:,0] > 0) & (binarymask.mf['CHUNKVECTOR'][:,0] < 7)) |\ ((binarymask.mf['CHUNKVECTOR'][:,1] > 0) & (binarymask.mf['CHUNKVECTOR'][:,1] < 7)) |\ ((binarymask.mf['CHUNKVECTOR'][:,2] > 0) & (binarymask.mf['CHUNKVECTOR'][:,2] < 7)) |\ ((binarymask.mf['CHUNKVECTOR'][:,3] > 0) & (binarymask.mf['CHUNKVECTOR'][:,3] < 7)))[0] if lowzchunkcutopt == 2: ## cut out chunks 3-6, keep chunk 2 and chunks >= 7 xxchunkcut = np.where(((binarymask.mf['CHUNKVECTOR'][:,0] > 0) & (binarymask.mf['CHUNKVECTOR'][:,0] >= 3) & (binarymask.mf['CHUNKVECTOR'][:,0] <= 6)) |\ ((binarymask.mf['CHUNKVECTOR'][:,1] > 0) & (binarymask.mf['CHUNKVECTOR'][:,1] >= 3) & (binarymask.mf['CHUNKVECTOR'][:,1] <= 6)) |\ ((binarymask.mf['CHUNKVECTOR'][:,2] > 0) & (binarymask.mf['CHUNKVECTOR'][:,2] >= 3) & (binarymask.mf['CHUNKVECTOR'][:,2] <= 6)) |\ ((binarymask.mf['CHUNKVECTOR'][:,3] > 0) & (binarymask.mf['CHUNKVECTOR'][:,3] >= 3) & (binarymask.mf['CHUNKVECTOR'][:,3] <= 6)))[0] if lowzchunkcutopt == 3: # cut out chunk 2, keep 3-6 and >= 7 xxchunkcut = np.where(((binarymask.mf['CHUNKVECTOR'][:,0] == 2)) |\ ((binarymask.mf['CHUNKVECTOR'][:,1] == 2)) |\ ((binarymask.mf['CHUNKVECTOR'][:,2] == 2)) |\ ((binarymask.mf['CHUNKVECTOR'][:,3] == 2)))[0] if lowzchunkcutopt == 4: ## make binary mask for only chunk 2 region. ## cut everything that would be kept in regular LOWZ mask and chunks 3-6. xxchunkcut = np.where((((binarymask.mf['CHUNKVECTOR'][:,0] <= 0) | (binarymask.mf['CHUNKVECTOR'][:,0] >= 7)) &\ ((binarymask.mf['CHUNKVECTOR'][:,1] <= 0) | (binarymask.mf['CHUNKVECTOR'][:,1] >= 7)) &\ ((binarymask.mf['CHUNKVECTOR'][:,2] <= 0) | (binarymask.mf['CHUNKVECTOR'][:,2] >= 7)) &\ ((binarymask.mf['CHUNKVECTOR'][:,3] <= 0) | (binarymask.mf['CHUNKVECTOR'][:,3] >= 7))) | ((binarymask.mf['CHUNKVECTOR'][:,0] >= 3) & (binarymask.mf['CHUNKVECTOR'][:,0] <= 6)))[0] if lowzchunkcutopt == 5: ## make binary mask for only chunk 3-6 region. ## cut everything that would be kept in regular LOWZ mask and chunk 2. xxchunkcut = np.where((((binarymask.mf['CHUNKVECTOR'][:,0] <= 0) | (binarymask.mf['CHUNKVECTOR'][:,0] >= 7)) &\ ((binarymask.mf['CHUNKVECTOR'][:,1] <= 0) | (binarymask.mf['CHUNKVECTOR'][:,1] >= 7)) &\ ((binarymask.mf['CHUNKVECTOR'][:,2] <= 0) | (binarymask.mf['CHUNKVECTOR'][:,2] >= 7)) &\ ((binarymask.mf['CHUNKVECTOR'][:,3] <= 0) | (binarymask.mf['CHUNKVECTOR'][:,3] >= 7))) | ((binarymask.mf['CHUNKVECTOR'][:,0] == 2)))[0] ## test -- make sure that every sector included in binary is in binary lowz, lowz2, or lowz3-6. #makebinarylowzearlytest(): ## passed this test! if lowzchunkcutopt > 0: cutseclist = np.unique(binarymask.mf['SECTOR'][xxchunkcut]) ## compute total area lost. alost = (binarymask.mf['STR'][xxchunkcut]).sum() alostchk = binarymask.area[cutseclist].sum() assert np.fabs(alost - alostchk)/alost < 2.0e-6 print 'chunk cut removes this much area:',alost*(180./np.pi)**2 binarymask.mf['WEIGHT'][xxchunkcut] = 0. binarymask.weight[cutseclist] = 0. if splitNSopt == 1: binarymaskN.mf['WEIGHT'][xxchunkcut] = 0. binarymaskS.mf['WEIGHT'][xxchunkcut] = 0. binarymaskN.weight[cutseclist] = 0. binarymaskS.weight[cutseclist] = 0. fitsio.write(binarymaskfitsout,binarymask.mf,clobber=True) if splitNSopt == 1: binarymaskfitsoutN = binarymaskfitsout.split('.fits')[0]+'-N.fits' binarymaskfitsoutS = binarymaskfitsout.split('.fits')[0]+'-S.fits' fitsio.write(binarymaskfitsoutN,binarymaskN.mf,clobber=True) fitsio.write(binarymaskfitsoutS,binarymaskS.mf,clobber=True) print 'total good area remaining in mask:',agood print 'total rejected mask area:',abad,'= fraction',abad/(abad+agood) return binarymask def makemaskcomp(collidedfile,binarymaskfits,maskoutbase,zfailopt=-1,minzcomp=0.8,verboseopt=0,minztarg=10): """ zfailopt = -1: don't cut on sector or plate redshift completeness, keep all sectors. DEFAULT! zfailopt = 0: match sector-by-sector. zfailopt = 1: match plate-by-plate Only consider sectors/plates with more than minztarg objects with imatch !=0,3,8. this is a departure from the previous choice of at least 10 imatch = 1 objects. """ # binarymask = mask(binarymaskfits) mcomp = mask(binarymaskfits) ## completeness mask we're going to generate here. # histbins = np.arange(-0.5,len(binarymask.weight)-0.5+0.1,1.0) dd = fitsio.read(collidedfile) wnew = np.where((dd['IMATCH'] != 0) & (dd['IMATCH'] != 2))[0] wnewnocp = np.where((dd['IMATCH'] != 0) & (dd['IMATCH'] != 2) & \ (dd['IMATCH'] != 3) & (dd['IMATCH'] != 8))[0] wall = np.where(dd['IMATCH'] != 2)[0] wznom = np.where(dd['IMATCH'] == 1)[0] wzdenom = np.where((dd['IMATCH'] == 1) | (dd['IMATCH'] == 5) | (dd['IMATCH'] == 7))[0] # xx = np.where((dd['IMATCH'] == 1) | (dd['IMATCH'] == 5) | (dd['IMATCH'] == 7) | (dd['IMATCH'] == 4))[0] platemin = dd['PLATE'][xx].min() assert platemin > 0 platemax = dd['PLATE'][xx].max() platebins = np.arange(-0.5,platemax+0.6,1.0) assert platebins[-1].astype('int') == platemax # print 'Calculating sector completeness...' print 'Ntarg - Nknown:',len(wall) print 'Nobs + Ncp:', len(wnew) print 'Nobs:', len(wnewnocp) # hall, xtmp = np.histogram(dd['ISECT'][wall],bins=histbins) hnew, xtmp = np.histogram(dd['ISECT'][wnew],bins=histbins) hnewnocp, xtmp = np.histogram(dd['ISECT'][wnewnocp],bins=histbins) # hznom, xtmp = np.histogram(dd['ISECT'][wznom],bins=histbins) hzdenom, xtmp = np.histogram(dd['ISECT'][wzdenom],bins=histbins) zcompsec = hznom.astype(float)/hzdenom.astype(float) xxx = np.where(hzdenom == 0)[0] zcompsec[xxx] = 1. mysec = np.where(zcompsec < minzcomp)[0] mysec2 = np.where((zcompsec < minzcomp) & (hnewnocp > minztarg))[0] # hznomplate, ytmp = np.histogram(dd['PLATE'][wznom],bins=platebins) hzdenomplate, ytmp = np.histogram(dd['PLATE'][wzdenom],bins=platebins) zcompplate = hznomplate.astype(float)/hzdenomplate.astype(float) xxx = np.where(hzdenomplate == 0)[0] zcompplate[xxx] = 1. if verboseopt == 1: print 'this many sectors, totalling this much area dont pass redshift failure cut' print len(np.where(zcompsec < minzcomp)[0]), binarymask.area[np.where(zcompsec < minzcomp)[0]].sum() print 'sectors:',mysec print 'sectors post 10 cut',mysec2 print len(mysec2), binarymask.area[mysec2].sum() print 'this many plates dont pass redshift failure cut',len(np.where(zcompplate < minzcomp)[0]) print 'no cuts here on number of targets in sector/plate!' print 'plates',np.where(zcompplate < minzcomp)[0] print 'cut plates have this many fibers' for pp in np.where(zcompplate < minzcomp)[0]: print pp, len(np.where(dd['PLATE'][wnewnocp] == pp)[0]) print zcompplate[np.where(zcompplate < minzcomp)[0]] print 'bad sector info start here' for ss in mysec2: xx = np.where(dd['ISECT'] == ss)[0] print ss, hznom[ss], hzdenom[ss] print zcompplate[dd['PLATE'][xx]] print 'end bad sector info' xcen = 0.5*(xtmp[1:] + xtmp[:-1]) ycen = 0.5*(ytmp[1:] + ytmp[:-1]) assert np.fabs(xcen[0] - binarymask.mf['SECTOR'].min()) < 2.e-6 assert np.fabs(xcen[-1] - binarymask.mf['SECTOR'].max()) < 2.e-6 assert np.fabs(ycen[-1] - dd['PLATE'].max()) < 2.e-6 comp = hnew.astype(float)/(hall.astype(float)+1.e-10) compnocp = hnewnocp.astype(float)/(hall.astype(float)+1.e-10) xxx = np.where(hall == 0)[0] comp[xxx] = 1. compnocp[xxx] = 1. print 'std completeness ranges from ',comp.min(), comp.max() print '1+w completeness ranges from ',compnocp.min(), compnocp.max() ## now check for redshift failure sectors and/or plates. if zfailopt == -1: ## do nothing! pass elif zfailopt == 0: xxx = np.where((zcomp < minzcomp) & (hnewnocp > minztarg))[0] comp[xxx] = 0. compnocp[xxx] = 0. elif zfailopt == 1: print 'this piece of code is untested!!! Caution!!' sys.exit(1) tiletosec = {} for ss in range(len(binarymask.area)): for tval in binarymask.tile[ss]: if tval in tiletosec.keys(): tiletosec[tval].append(ss) else: tiletosec[tval] = [ss] badplatelist = np.where(zcompplate < minzcomp)[0] badtilelist = [] for pp in badplatelist: xx = np.where(dd['PLATE'] == pp)[0] yy = np.where((dd['PLATE'] == pp) & (dd['IMATCH'] != 0) & (dd['IMATCH'] != 2) & \ (dd['IMATCH'] != 3) & (dd['IMATCH'] != 8))[0] if len(yy) <= minztarg: continue ## else, passed minztarg. tilelist = np.unique(dd['SPECTILE'][xx])[0] assert len(tilelist) == 1 badtilelist.append(tilelist[0]) badseclist = [] for badtile in badtilelist: for ss in tiletosec[badtile]: badseclist.append(ss) badseclist = np.unique(np.array(badseclist)) comp[badseclist] = 0. compnocp[badseclist] = 0. else: print 'zfailopt = ',zfailopt,'not supported.' ## output new mask with completeness. ## fill new completeness in every sector, then reapply binary mask. mcomp.weight = comp mcomp.mf['WEIGHT'] = comp[mcomp.mf['SECTOR']] xs = np.where(binarymask.weight < 0.001)[0] xxq = np.where(binarymask.mf['WEIGHT'] < 0.001)[0] print 'binary mask excludes this many sectors',len(xs),'and this many polygons',len(xxq) mcomp.weight[xs] = 0. mcomp.mf['WEIGHT'][xxq] = 0. ## write out fits and ply versions. fitsio.write(maskoutbase + '.fits', mcomp.mf, clobber = True) return mcomp ## supporting functions for makecombinedmask def getcollidedtrim(collidedfile,bmask): """ Read in collided file, trim down to appropriate mask. """ cc = fitsio.read(collidedfile) xx = np.where(bmask.weight[cc['ISECT']] > 0.001)[0] cc = cc[xx] return cc def gethistcomp(dd,histbins): """ Get histograms for completeness calculation for any input data vector. """ wall = np.where(dd['IMATCH'] != 2)[0] wnew = np.where((dd['IMATCH'] != 0) & (dd['IMATCH'] != 2))[0] hall, xtmp = np.histogram(dd['ISECT'][wall],bins=histbins) hnew, xtmp = np.histogram(dd['ISECT'][wnew],bins=histbins) xcen = 0.5*(xtmp[1:] + xtmp[:-1]) return xcen, hnew, hall def getrepeats(dlong,dshort,histbins): targ_objid = {} t1 = time.time() for i in range(len(dlong)): ## native time for RERUN is '|S12', which disagrees with spAll. ## these are both derived from the collate file, so this should work. targ_objid[(dlong['RUN'][i], dlong['RERUN'][i], dlong['CAMCOL'][i], dlong['FIELD'][i], dlong['ID'][i])] = i xalllist = [] xnewlist = [] for i in range(len(dshort)): photoid = photokey(dshort[i]) if photoid in targ_objid: ilong = targ_objid[photoid] assert photokey(dlong[ilong]) == photoid assert dshort['ISECT'][i] == dlong['ISECT'][ilong] if dshort['IMATCH'][i] != 2 and dlong['IMATCH'][ilong] != 2: assert dshort['ISECT'][i] == dlong['ISECT'][ilong] xalllist.append(i) if dshort['IMATCH'][i] != 0 and dlong['IMATCH'][ilong] != 0: xnewlist.append(i) xalllist = np.array(xalllist,dtype='int') xnewlist = np.array(xnewlist,dtype='int') hnewR, xtmp = np.histogram(dshort['ISECT'][xnewlist],bins=histbins) hallR, xtmp = np.histogram(dshort['ISECT'][xalllist],bins=histbins) print 'found this many CMASS/LOWZ overlaps',len(xalllist), len(xnewlist) xcen = 0.5*(xtmp[1:] + xtmp[:-1]) return xcen, hnewR, hallR def joincombineNmasks(mLfitsfname,mL2,mL3,mLcpfitsfname,mL2cp,mL3cp,bL2cut,bL3cut,bC,mfitsOUT,mfitscpOUT): """ put this into mksamplemaster.py at the end. Name a combinedN sample with the task of concatenating the 3 N samples? """ mtot = mask(mLfitsfname) mcptot = mask(mLcpfitsfname) xx = np.where(bC.weight < 0.0001)[0] for mm in [mtot, mL2, mL3]: assert (mm.weight[xx] < 0.0001).all() ## fill in mask values for chunk 2. yy = np.where(bL2cut.weight > 0.001)[0] mtot.weight[yy] = mL2.weight[yy] mcptot.weight[yy] = mL2cp.weight[yy] for ss in yy: xxq = np.where(mtot.mf['SECTOR'] == ss)[0] assert len(xxq) > 0 mtot.mf['WEIGHT'][xxq] = mL2.mf['WEIGHT'][xxq] mcptot.mf['WEIGHT'][xxq] = mL2cp.mf['WEIGHT'][xxq] assert (mL2.mf['SECTOR'][xxq] == ss).all() ## fill in mask values for chunk 3. yy = np.where(bL3cut.weight > 0.001)[0] mtot.weight[yy] = mL3.weight[yy] mcptot.weight[yy] = mL3cp.weight[yy] for ss in yy: xxq = np.where(mtot.mf['SECTOR'] == ss)[0] assert len(xxq) > 0 mtot.mf['WEIGHT'][xxq] = mL3.mf['WEIGHT'][xxq] mcptot.mf['WEIGHT'][xxq] = mL3cp.mf['WEIGHT'][xxq] assert (mL3.mf['SECTOR'][xxq] == ss).all() fitsio.write(mfitsOUT,mtot.mf,clobber=True) fitsio.write(mfitscpOUT,mcptot.mf,clobber=True) def makecombinedmask(NorS,collidedfileCMASS,collidedfileLOWZ,collidedfileLOWZ2,collidedfileLOWZ3,binarymaskfitsCMASS,binarymaskfitsLOWZ,binarymaskfitstrimLOWZ2, binarymaskfitstrimLOWZ3, maskoutbase,zfailopt=-1,mCf=None,mLf=None,mL2f=None,mL3f=None): """ Kludgey way to make a combined CMASS/LOWZ mask. JT: make sure we don't double-count objects that are CMASS AND LOWZ. BR: this is being de-commissioned, but I will compare the results with hte new way to do this using the joined cmass-lowz samples. """ if zfailopt != -1: print 'zfailopt = ',zfailopt,'not supported!' sys.exit(1) bC = mask(binarymaskfitsCMASS) bL = mask(binarymaskfitsLOWZ) histbins = np.arange(-0.5,len(bC.weight)-0.5+0.1,1.0) dC = getcollidedtrim(collidedfileCMASS,bC) dL = getcollidedtrim(collidedfileLOWZ,bL) xcen, hnewC, hallC = gethistcomp(dC,histbins) xcen, hnewL, hallL = gethistcomp(dL,histbins) assert np.fabs(xcen[0] - bC.mf['SECTOR'].min()) < 2.e-6 assert np.fabs(xcen[-1] - bC.mf['SECTOR'].max()) < 2.e-6 hnewtot = np.zeros_like(hnewC) halltot = np.zeros_like(hallC) xcen, hnewR, hallR = getrepeats(dC,dL,histbins) if NorS == 0: bL2 = mask(binarymaskfitstrimLOWZ2) bL3 = mask(binarymaskfitstrimLOWZ3) dL2 = getcollidedtrim(collidedfileLOWZ2,bL2) dL3 = getcollidedtrim(collidedfileLOWZ3,bL3) xcen, hnewL2, hallL2 = gethistcomp(dL2,histbins) xcen, hnewL3, hallL3 = gethistcomp(dL3,histbins) xcen, hnewR2, hallR2 = getrepeats(dC,dL2,histbins) xcen, hnewR3, hallR3 = getrepeats(dC,dL3,histbins) halltot[:] = hallC[:] + hallL[:] + hallL2[:] + hallL3[:] - hallR[:] - hallR2[:] - hallR3[:] hnewtot[:] = hnewC[:] + hnewL[:] + hnewL2[:] + hnewL3[:] - hnewR[:] - hnewR2[:] - hnewR3[:] print 'total repeats:',(hallR + hallR2 + hallR3).sum(), (hnewR + hnewR2 + hnewR3).sum(), halltot.sum(), hnewtot.sum() else: halltot[:] = hallC[:] + hallL[:] - hallR[:] hnewtot[:] = hnewC[:] + hnewL[:] - hnewR[:] print 'total repeats:',(hallR).sum(), (hnewR).sum(), halltot.sum(), hnewtot.sum() ## tmp! for ss in range(len(halltot)): print ss, hallC[ss],hallL[ss],hallL2[ss],hallL3[ss],hallR[ss],hallR2[ss], hallR3[ss] if hallL[ss] > 0: assert hallL2[ss] == 0 assert hallL3[ss] == 0 if hallL2[ss] > 0: assert hallL[ss] == 0 assert hallL3[ss] == 0 if hallL3[ss] > 0: assert hallL[ss] == 0 assert hallL2[ss] == 0 comp = hnewtot.astype(float)/(halltot.astype(float)+1.e-10) xxx = np.where(halltot == 0)[0] comp[xxx] = 1. print 'area of sectors with no targets:',bC.area[xxx].sum()*(180./np.pi)**2 ## output mask. mcomp = mask(binarymaskfitsCMASS) ## check that there are no sectors that are included in lowz masks but not CMASS. xx = np.where(bC.weight < 0.001)[0] assert (bL.weight[xx] < 0.001).all() if NorS == 0: assert (bL2.weight[xx] < 0.001).all() assert (bL3.weight[xx] < 0.001).all() ## output new mask with completeness. ## fill new completeness in every sector, then reapply binary mask. mcomp.weight = comp mcomp.mf['WEIGHT'] = comp[mcomp.mf['SECTOR']] xs = np.where(bC.weight < 0.001)[0] xxq = np.where(bC.mf['WEIGHT'] < 0.001)[0] print 'binary mask excludes this many sectors',len(xs),'and this many polygons',len(xxq) mcomp.weight[xs] = 0. mcomp.mf['WEIGHT'][xxq] = 0. ## write out fits and ply versions. fitsio.write(maskoutbase + '.fits', mcomp.mf, clobber = True) ## now run some tests. ## need to read in these for this test: mCf,mLf,mL2f,mL3f if mCf is not None: ## then run tests. print 'here beth!!' mlist = [] if NorS == 0: mylist = [mCf, mLf, mL2f, mL3f] else: mylist = [mCf, mLf] for mf in mylist: print 'working on ',mf if mf is None: print 'aborting tests, returning mcomp!' return mcomp mm = mask(mf) xgoodnew = np.where(mm.weight > 0.0001)[0] if (mm.area[xgoodnew]).sum() > 0: print 'mX mean comp:',(mm.weight[xgoodnew]*mm.area[xgoodnew]).sum()/(mm.area[xgoodnew]).sum() mlist.append(mm) xx = np.where((mcomp.weight < 0.001) & (mm.weight > 0.001))[0] if len(xx) > 0: print mf,len(xx), xx[0], mcomp.weight[xx[0]], mm.weight[xx[0]] print 'this much lost area!',mcomp.area[xx].sum()*(180./np.pi)**2 wackarea = 0. for ii in range(len(mcomp.weight)): xxq = np.where(mcomp.mf['SECTOR'] == ii)[0] cv = mcomp.mf['CHUNKVECTOR'][xxq[0],0] if mcomp.weight[ii] < 0.001: continue if cv == 2: tmplist = np.array([mlist[0].weight[ii], mlist[2].weight[ii]]) elif cv >= 7: tmplist = np.array([mlist[0].weight[ii], mlist[1].weight[ii]]) else: tmplist = np.array([mlist[0].weight[ii], mlist[3].weight[ii]]) if np.fabs(mcomp.weight[ii] - 1) < 2.0e-6 and np.fabs(tmplist.min() - 1) < 2.0e-6: continue if not (mcomp.weight[ii] >= tmplist.min() and mcomp.weight[ii] <= tmplist.max()): print mcomp.weight[ii], tmplist, cv, ii, mcomp.mf['CHUNKVECTOR'][xxq[0],:] if NorS == 0: print ii, hallC[ii], hallL[ii], hallL2[ii], hallL3[ii], hallR[ii], hallR2[ii], hallR3[ii] print ii, hnewC[ii], hnewL[ii], hnewL2[ii], hnewL3[ii], hnewR[ii], hnewR2[ii], hnewR3[ii] else: print ii, hallC[ii], hallL[ii], hallR[ii] print ii, hnewC[ii], hnewL[ii], hnewR[ii] wackarea += mcomp.area[ii] ## i guess this doesn't have to be true either, if there's one of the double-counted objects in this sector. check by inspection. ##assert mcomp.weight[ii] >= tmplist.min() and mcomp.weight[ii] <= tmplist.max() ## yep, if double counted are full completeness or below min completeness, things can swing below/above. ## checked this by hand with sector 1839. so i think the new combined masks were fine. print 'wack area',wackarea*(180./np.pi)**2 print 'combined mask test passed!' return mcomp def makecombinedmasktest(mcombf,mCf,mLf,mL2f,mL3f,NorS): """ These tests have been moved up now to makecombinemask. """ mcomb = mask(mcombf) xgood = np.where(mcomb.weight > 0.0001)[0] print 'combined mean comp:',(mcomb.weight[xgood]*mcomb.area[xgood]).sum()/(mcomb.area[xgood]).sum() mlist = [] if NorS == 0: mylist = [mCf, mLf, mL2f, mL3f] else: mylist = [mCf, mLf] for mf in mylist: print 'working on ',mf mm = mask(mf) xgoodnew = np.where(mm.weight > 0.0001)[0] if (mm.area[xgoodnew]).sum() > 0: print 'mX mean comp:',(mm.weight[xgoodnew]*mm.area[xgoodnew]).sum()/(mm.area[xgoodnew]).sum() mlist.append(mm) xx = np.where((mcomb.weight < 0.001) & (mm.weight > 0.001))[0] if len(xx) > 0: print mf,len(xx), xx[0], mcomb.weight[xx[0]], mm.weight[xx[0]] print 'this much lost area!',mcomb.area[xx].sum()*(180./np.pi)**2 #assert len(xx) == 0 wackarea = 0. for ii in range(len(mcomb.weight)): xxq = np.where(mcomb.mf['SECTOR'] == ii)[0] cv = mcomb.mf['CHUNKVECTOR'][xxq[0],0] if mcomb.weight[ii] < 0.001: continue if cv == 2: tmplist = np.array([mlist[0].weight[ii], mlist[2].weight[ii]]) elif cv >= 7: tmplist = np.array([mlist[0].weight[ii], mlist[1].weight[ii]]) else: tmplist = np.array([mlist[0].weight[ii], mlist[3].weight[ii]]) if np.fabs(mcomb.weight[ii] - 1) < 2.0e-6 and np.fabs(tmplist.min() - 1) < 2.0e-6: continue if not (mcomb.weight[ii] >= tmplist.min() and mcomb.weight[ii] <= tmplist.max()): print mcomb.weight[ii], tmplist, cv, ii, mcomb.mf['CHUNKVECTOR'][xxq[0],:] wackarea += mcomb.area[ii] ## i guess this doesn't have to be true either, if there's one of the double-counted objects in this sector. check by inspection. ##assert mcomb.weight[ii] >= tmplist.min() and mcomb.weight[ii] <= tmplist.max() ## yep, if double counted are full completeness or below min completeness, things can swing below/above. print 'wack area',wackarea*(180./np.pi)**2 print 'combined mask test passed!' def checkbinary(bC,bL,bL2,bL2cut,bL3,bL3cut,NorS,inputfiles=0): if inputfiles == 1: bC = mask.maskutils(bC) bL = mask.maskutils(bL) if NorS == 0: bL2 = mask.maskutils(bL2) bL2cut = mask.maskutils(bL2cut) bL3 = mask.maskutils(bL3) bL3cut = mask.maskutils(bL3cut) ## else you're inputing the masks directly. if NorS == 0: blist = [bL, bL2, bL2cut, bL3, bL3cut] else: blist = [bL] for ss in range(0,len(bC.weight)): ## loop over the sectors. mylogic = [] if bC.weight[ss] < 0.001: for bb in blist: if not bb.weight[ss] < 0.001: cv = bC.mf['CHUNKVECTOR'][xxq[0],:] print 'wtf!?!', NorS, ss, cv, bC.weight[ss], bL.weight[ss], bL2cut.weight[ss], bL3cut.weight[ss], bL2.weight[ss], bL3.weight[ss] assert bb.weight[ss] < 0.001 else: xxq = np.where(bC.mf['SECTOR'] == ss)[0] assert len(xxq) > 0 cv = bC.mf['CHUNKVECTOR'][xxq[0],:] if NorS == 0: if np.fabs(bL.weight[ss] + bL2cut.weight[ss] + bL3cut.weight[ss] - 1.) < 2.0e-6: continue else: print NorS, ss, cv, bC.weight[ss], bL.weight[ss], bL2cut.weight[ss], bL3cut.weight[ss], bL2.weight[ss], bL3.weight[ss] assert bL.weight[ss] > 0.001 or bL2.weight[ss] > 0.001 or bL3.weight[ss] > 0.001 assert np.fabs(bL.weight[ss] + bL2.weight[ss] + bL3.weight[ss] - 1.) < 2.0e-6 else: if bL.weight[ss] <= 0.001: print NorS, ss, cv, bC.weight[ss], bL.weight[ss] assert bL.weight[ss] > 0.001 print 'passed Nors = ',NorS def makechunkmask(binaryfitsfname,geomfitsfname,chunkstart,chunkend): """ This function will generate a mask that includes sectors from binary file and also require them to be in the stated chunk range. Geomfitsfname required by maskfitstoply. """ bb = mask(binaryfitsfname) xx = np.where((bb.mf['CHUNKVECTOR'][:,0] < chunkstart) | (bb.mf['CHUNKVECTOR'][:,0] > chunkend))[0] bb.mf['WEIGHT'][xx] = 0. binarymaskfitsout = binaryfitsfname.split('.fits')[0] + 'chunk%02d-chunk%02d' % (chunkstart,chunkend) + '.fits' fitsio.write(binarymaskfitsout,bb.mf,clobber=True) maskfitstoply(binarymaskfitsout,geomfitsfname) def makechunkmaskold(geomfitsfname = "/home/bareid/boss/bosslss/trunk/geometry/boss_geometry_2014_03_12.fits",\ targfitsfname = "/home/bareid/boss/bosslss/trunk/data/trimmed-collate-cmass-dr12v1.fits",\ outfitsfname = None): """ Set outfitsfname = None if you don't want to write out the mask. This matching is approximate -- in the overlaps between chunks, there are multiple ichunk values on the same tile. Why Jeremy? Ack this is stupid! There is a chunkvector in the mask, use that instead! """ geom = mask(geomfitsfname) targ = fitsio.read(targfitsfname) maxsect = geom.mf['SECTOR'].max() alltiles = np.unique(np.concatenate(geom.tiles)) assert (alltiles[:-1] < alltiles[1:]).all() tiletochunk = {} for tt in alltiles: ## get objects only on this tile, figure out their chunk. xx = np.where((targ['TILE'][:,0] == tt) & (targ['TILE'][:,1] == -1) & (targ['TILE'][:,2] == -1))[0] yy = np.where((targ['TILE'][:,0] == tt) | (targ['TILE'][:,1] == tt) | (targ['TILE'][:,2]) == tt)[0] if len(yy) > 0: if len(xx) == 0: print tt, targ['TILE'][yy] assert len(xx) > 0 if (targ['ICHUNK'][xx] == targ['ICHUNK'][xx[0]]).all(): print 'ichunk/tile agree',tt,len(xx) tiletochunk[tt] = targ['ICHUNK'][xx[0]] else: print 'ichunk/tile disagree',tt, len(xx) print targ['ICHUNK'][xx] htmp, xtmp = np.histogram(targ['ICHUNK'][xx],bins=np.arange(targ['ICHUNK'][xx].min() - 1.5, targ['ICHUNK'][xx].max() +1.5,1.0)) xcen = 0.5*(xtmp[1:] + xtmp[:-1]) bestichunk = int(np.floor(xcen[np.where(htmp == htmp.max())[0]] + 0.5)) print 'best ichunk:',bestichunk print 'histogram:' print xcen print htmp assert bestichunk in targ['ICHUNK'][xx] #assert (targ['ICHUNK'][xx] == targ['ICHUNK'][xx[0]]).all() tiletochunk[tt] = targ['ICHUNK'][xx[0]] else: tiletochunk[tt] = -1 ichunksec = np.zeros(len(geom.area),dtype='int') ## now go through sectors and make sure all tiles in the sector agree on the chunk. for s in range(len(geom.tiles)): if len(geom.tiles[s]) == 1: if geom.tiles[s][0] in tiletochunk.keys(): ichunksec[s] = tiletochunk[geom.tiles[s][0]] else: ichunksec[s] = -1 elif len(geom.tiles[s]) == 0: ichunksec[s] = -1 else: clist = [] for tt in geom.tiles[s]: if tt in tiletochunk.keys(): clist.append(tiletochunk[tt]) if len(clist) == 0: ichunksec[s] = -1 else: #assert (np.array(clist) == clist[0]).all() #ichunksec[s] = clist[0] ichunksec[s] = -1 ## who knows! print 'ichunksec computed successfully!' for s in range(maxsect+1): xxq = np.where(geom.mf['SECTOR'] == s)[0] assert len(xxq) > 0 geom.weight[s] = ichunksec[s] geom.mf['WEIGHT'][xxq] = ichunksec[s] if outfitsfname is not None: fitsio.write(outfitsfname,geom.mf,clobber=True) return geom def comparemasksbinary(m1,m2): """ Check on if the sectors agree on what masks to include. """ if len(m1.area) != len(m2.area): #same number of sectors. print 'number of sectors not the same, aborting!!' sys.exit(1) if ((m1.area) != (m2.area)).any(): #same area in sectors. print 'area of sectors not the same, aborting!!' sys.exit(1) a = np.zeros(4) for s in range(len(m1.area)): if m1.weight[s] > 0.001 and m2.weight[s] > 0.001: a[0] += m1.area[s] elif m1.weight[s] < 0.001 and m2.weight[s] < 0.001: a[3] += m1.area[s] elif m1.weight[s] > 0.001 and m2.weight[s] < 0.001: a[1] += m1.area[s] print 'sector ',s, 'disagree: m1=%e,m2=%e.' % (m1.weight[s], m2.weight[s]),'tiles are ',m1.tiles[s] elif m1.weight[s] < 0.001 and m2.weight[s] > 0.001: a[2] += m1.area[s] print 'sector ',s, 'disagree: m1=%e,m2=%e.' % (m1.weight[s], m2.weight[s]),'tiles are ',m1.tiles[s] else: assert 0==1 ## switch to degreessqr fac = (180./np.pi)**2 a = a*fac print 'agree in mask',a[0] print 'agree not in mask',a[3] print 'disagree m1 in, m2 out',a[1] print 'disagree m1 out, m2 in',a[2] print 'fractional disagree:',(a[1]+a[2])/(a[0]+a[1]+a[2]) def comparemaskscomp(m1,m2,c1,c2): """ For sectors included in both masks, do they get the same value of completeness? c is input collided structure. """ if len(m1.area) != len(m2.area): #same number of sectors. print 'number of sectors not the same, aborting!!' sys.exit(1) if ((m1.area) != (m2.area)).all(): #same area in sectors. print 'number of sectors not the same, aborting!!' sys.exit(1) a = np.zeros(4) fac = (180./np.pi)**2 cagree = 0. cdisagree = 0. for s in range(len(m1.area)): xtmp = np.where(c1['ISECT'] == s)[0] if m1.weight[s] > 0.001 and m2.weight[s] > 0.001: a[0] += m1.area[s] if np.fabs(m1.weight[s] - m2.weight[s]) > 2.0e-6: cdisagree += m1.area[s] print 'sector disagree',s,m1.area[s]*fac, m1.weight[s], m2.weight[s] ## if they disagree, is it because of imatch2/imatch8 difference? ytmp = np.where(c1['IMATCH'][xtmp] != c2['IMATCH'][xtmp])[0] if ((c1['IMATCH'][xtmp[ytmp]] == 2) | (c2['IMATCH'][xtmp[ytmp]] == 2)).any() or ((c1['IMATCH'][xtmp[ytmp]] == 0) & (c2['IMATCH'][xtmp[ytmp]] == 3)).any() or ((c1['IMATCH'][xtmp[ytmp]] == 3) & (c2['IMATCH'][xtmp[ytmp]] == 0)).any() or (c1['IMATCH'][xtmp[ytmp]] == 8).any() or (c2['IMATCH'][xtmp[ytmp]] == 8).any(): pass ## the reasons IMATCH values could disagree are: ## a) disagree about fiber collision method (groups vs 62'' pairs) ## b) disagree object-by-object on stochastic corrections of legacy-legacy or legacy-cmass FB corrections. ## c) disagree because I switched to mask completeness calculation step before legacy pair corretion step. ## else: print 'remaining bug??',c1['IMATCH'][xtmp[ytmp]], c2['IMATCH'][xtmp[ytmp]] else: cagree += m1.area[s] elif m1.weight[s] < 0.001 and m2.weight[s] < 0.001: a[3] += m1.area[s] elif m1.weight[s] > 0.001 and m2.weight[s] < 0.001: a[1] += m1.area[s] print 'sector ',s, 'disagree: m1=%e,m2=%e.' % (m1.weight[s], m2.weight[s]),'tiles are ',m1.tiles[s] elif m1.weight[s] < 0.001 and m2.weight[s] > 0.001: a[2] += m1.area[s] print 'sector ',s, 'disagree: m1=%e,m2=%e.' % (m1.weight[s], m2.weight[s]),'tiles are ',m1.tiles[s] else: assert 0==1 a = a*fac cagree = cagree*fac cdisagree = cdisagree*fac print 'agree in mask',a[0] print 'agree not in mask',a[3] print 'disagree m1 in, m2 out',a[1] print 'disagree m1 out, m2 in',a[2] print 'perfect agreement in completeness over this area:',cagree print 'disagreement in completeness over this area:',cdisagree ## compute average value of completeness in two surveys x1 = np.where(m1.weight > 0.001)[0] x2 = np.where(m2.weight > 0.001)[0] print 'average completeness inside the two mask (area weighted): ',(m1.weight[x1]*m1.area[x1]).sum()/m1.area[x1].sum(), (m1.weight[x2]*m1.area[x2]).sum()/m2.area[x2].sum() return 0 ################################## # start functions I wrote while testing. ################################## def makebinarylowzearlytest(): ddir = "/home/bareid/boss/bosslss/trunk/data/" bt = mask(ddir+"binarymask-dr12v5.fits") b1 = mask(ddir+"binarymask-dr12v5-lowzchunkcut.fits") b2 = mask(ddir+"binarymask-dr12v5-lowzearlychunk3-6targ-chunk3-6-N.fits") b3 = mask(ddir+"binarymask-dr12v5-lowzearlychunk2targ-chunk2-N.fits") print 'read masks' btmp = mask(ddir+"binarymask-dr12v5-lowzearlychunk3-6targ-chunk3-6-S.fits") assert (btmp.mf['WEIGHT'] == 0.0).all() btmp = mask(ddir+"binarymask-dr12v5-lowzearlychunk2targ-chunk2-S.fits") assert (btmp.mf['WEIGHT'] == 0.0).all() print 'passed S test.' ## make sure for each good sector exactly one lowz mask covers it. assert (bt.mf['WEIGHT'][:] == (b1.mf['WEIGHT'][:] + b2.mf['WEIGHT'][:] + b3.mf['WEIGHT'][:])).all() print 'passed LOWZ + LOWZearly test.' def makebinarytest(spall_file = "/home/bareid/boss/bosslss/trunk/mksample/inputfiles/spAll-v5_7_0.fits", \ platelistfits = "/home/bareid/boss/bosslss/trunk/mksample/inputfiles/platelist.fits", \ targ_file = "/home/bareid/boss/bosslss/trunk/data/trimmed-collate-cmass-dr12v1.fits",\ maskfitsfname = "/home/bareid/boss/bosslss/trunk/geometry/boss_geometry_2014_03_12.fits",\ binarymaskfitsout = "/home/bareid/boss/bosslss/trunk/data/binarymaskNS-dr12v1try2.fits"): pp = fitsio.read(platelistfits) goodtilelist = np.unique(pp[np.where((pp['PLATEQUALITY'] == 'good') & \ (pp['PROGRAMNAME'] == 'boss '))[0]]['TILEID']) binarymask = mask(maskfitsfname) maxsect = binarymask.mf['SECTOR'].max() ## check which tiles are not in goodtilelist. alltiles = np.unique(np.concatenate(binarymask.tiles)) assert (alltiles[:-1] < alltiles[1:]).all() ## are all good tiles in alltiles? cnt1 = 0 cnt2 = 0 for tileval in goodtilelist: if tileval in alltiles: cnt1 += 1 else: cnt2 += 1 print tileval,'not in alltiles' print cnt1,'goodtiles in alltiles' print cnt2,'goodtiles not in alltiles' ## slow but works? tilecheck = np.zeros(alltiles.max()+1,dtype='int') for tileval in alltiles: if tileval in goodtilelist: ## adopt list from platelist as standard for now. tilecheck[tileval] = 1 else: tilecheck[tileval] = 0 ## now loop over all sectors, see if their tile values are all observed. for s in range(maxsect+1): xxq = np.where(binarymask.mf['SECTOR'] == s)[0] assert len(xxq) > 0 if (tilecheck[binarymask.tiles[s][:]] == 1).all(): binarymask.weight[s] = 1. binarymask.mf['WEIGHT'][xxq] = 1. else: binarymask.weight[s] = 0. binarymask.mf['WEIGHT'][xxq] = 0. print 'sector ',s,'is missing a tile', binarymask.tiles[s], tilecheck[binarymask.tiles[s][:]] fitsio.write(binarymaskfitsout,binarymask.mf,clobber=True) return binarymask def makebinarytestv0(spall_file = "/home/bareid/boss/bosslss/trunk/mksample/inputfiles/spAll-v5_7_0.fits", \ platelistfits = "/home/bareid/boss/bosslss/trunk/mksample/inputfiles/platelist.fits", \ targ_file = "/home/bareid/boss/bosslss/trunk/data/trimmed-collate-cmass-dr12v1.fits",\ maskfitsfname = "/home/bareid/boss/bosslss/trunk/geometry/boss_geometry_2014_03_12.fits",\ binarymaskfitsout = "/home/bareid/boss/bosslss/trunk/data/binarymaskNS-dr12v1try2.fits"): """ Upgrade this to a list of good TILE values from Stephen if it exists. """ dd_spec = fitsio.read(spall_file,columns=['RA','DEC','Z_NOQSO','ZWARNING_NOQSO','SPECPRIMARY','TILE'],iter_row_buffer=10000) dd_spec = dd_spec[np.where((dd_spec['SPECPRIMARY'] == 1) & ( dd_spec['ZWARNING_NOQSO'] == 0))[0]] print 'this many good redshifts in spall:',len(dd_spec) goodtilelist = np.setdiff1d(dd_spec['TILE'], np.array([-3,0])) goodtilelist.sort() ## tmp! ofp = open('goodtilelist1.txt','w') for i in range(len(goodtilelist)): ofp.write('%d\n' % (goodtilelist[i])) ofp.close() pp = fitsio.read(platelistfits) goodtilelist2 = np.unique(pp[np.where((pp['PLATEQUALITY'] == 'good') & (pp['PROGRAMNAME'] == 'boss '))[0]]['TILEID']) ofp = open('goodtilelist2.txt','w') for i in range(len(goodtilelist2)): ofp.write('%d\n' % (goodtilelist2[i])) ofp.close() ## those two agree I guess? will find out soon. binarymask = mask.mask(maskfitsfname) maxsect = binarymask.mf['SECTOR'].max() ## check which tiles are not in goodtilelist. alltiles = np.unique(np.concatenate(binarymask.tiles)) assert (alltiles[:-1] < alltiles[1:]).all() ## slow but works? tilecheck = np.zeros(alltiles.max()+1,dtype='int') for tileval in alltiles: if tileval in goodtilelist2: ## adopt list from platelist as standard for now. tilecheck[tileval] = 1 else: tilecheck[tileval] = 0 ## now loop over all sectors, see if their tile values are all observed. for s in range(maxsect+1): xxq = np.where(binarymask.mf['SECTOR'] == s)[0] assert len(xxq) > 0 if (tilecheck[binarymask.tiles[s][:]] == 1).all(): binarymask.weight[s] = 1. binarymask.mf['WEIGHT'][xxq] = 1. else: binarymask.weight[s] = 0. binarymask.mf['WEIGHT'][xxq] = 0. print 'sector s is missing a tile' print binarymask.tiles[s] print tilecheck[binarymask.tiles[s][:]] fitsio.write(binarymaskfitsout,binarymask.mf,clobber=True) return binarysec if __name__ == '__main__': makebinary() #makebinarytest() #chunkmask = makechunkmask()