diff --git a/pipeline/applyMaskToMock.py b/pipeline/applyMaskToMock.py deleted file mode 100755 index af7723c..0000000 --- a/pipeline/applyMaskToMock.py +++ /dev/null @@ -1,392 +0,0 @@ -#!/usr/bin/env python - -# applies a mask to a given dataset - -import numpy as np -import os -import sys -import void_python_tools as vp -import argparse -import imp -import healpy as hp - -# ------------------------------------------------------------------------------ - -def my_import(name): - mod = __import__(name) - components = name.split('.') - for comp in components[1:]: - mod = getattr(mod, comp) - return mod - -# ----------------------------------------------------------------------------- - -LIGHT_SPEED = 299792.458 - -parser = argparse.ArgumentParser(description='options') -parser.add_argument('--scripts', dest='script', action='store_const', - const=True, default=False, - help='write scripts') -parser.add_argument('--parmFile', dest='parmFile', - default="", - help='path to parameter file') -parser.add_argument('--subsamples', dest='subsample', action='store_const', - const=True, default=False, - help='write subsamples') -parser.add_argument('--halos', dest='halos', action='store_const', - const=True, default=False, - help='write halos') -parser.add_argument('--hod', dest='hod', action='store_const', - const=True, default=False, - help='write hod') -parser.add_argument('--all', dest='all', action='store_const', - const=True, default=False, - help='write everything') -args = parser.parse_args() - - -filename = args.parmFile -print " Loading parameters from", filename -if not os.access(filename, os.F_OK): - print " Cannot find parameter file %s!" % filename - exit(-1) -parms = imp.load_source("name", filename) -globals().update(vars(parms)) - - -#------------------------------------------------------------------------------ -def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): - - sampleName = setName - - sampleName += "_z" + redshift - - if iVol != -1: sampleName += "_d" + iVol - - return sampleName - -#------------------------------------------------------------------------------ -# for given dataset parameters, outputs a script for use with analyzeVoids -def writeObservationScript(setName, dataFileNameBase, maskFileName, - scriptDir, catalogDir, fileNums, redshifts, - useVel, minRadius, omegaM, - suffix=".dat"): - - - if useVel: setName += "_pv" - - scriptFileName = scriptDir + "/" + setName + ".py" - scriptFile = open(scriptFileName, 'w') - - scriptFile.write("""#!/usr/bin/env/python -import os -from void_python_tools.backend.classes import * - -continueRun = True # set to True to enable restarting aborted jobs -startCatalogStage = 1 -endCatalogStage = 4 - -startAPStage = 1 -endAPStage = 7 - -ZOBOV_PATH = os.getenv("PWD")+"/../zobov/" -CTOOLS_PATH = os.getenv("PWD")+"/../c_tools/" -freshStack = True -errorBars = "CALCULATED" -numIncoherentRuns = 100 -ranSeed = 101010 -useLCDM = False -bias = 1.16 - -dataPortions = ["central", "all"] -dataSampleList = [] -""") - - - dataInfo = """ -setName = "{setName}" - -workDir = "{voidOutputDir}/{setName}/" -inputDataDir = "{inputDataDir}" -figDir = "{figDir}/{setName}/" -logDir = "{logDir}/{setName}/" - -numZobovDivisions = {numZobovDivisions} -numZobovThreads = {numZobovThreads} - """ - scriptFile.write(dataInfo.format(setName=setName, figDir=figDir, - logDir=logDir, voidOutputDir=voidOutputDir, - inputDataDir=catalogDir, - numZobovDivisions=numZobovDivisions, - numZobovThreads=numZobovThreads)) - - sampleInfo = """ -newSample = Sample(dataFile = "{dataFile}", - fullName = "{sampleName}", - nickName = "{sampleName}", - dataType = "observation", - maskFile = "{maskFileName}", - zBoundary = ({zMin}, {zMax}), - zRange = ({zMin}, {zMax}), - minVoidRadius = {minRadius}, - includeInHubble = True, - partOfCombo = False, - isCombo = False, - usePecVel = {usePecVel}) -dataSampleList.append(newSample) -newSample.addStack({zMin}, {zMax}, {minRadius} , {minRadius}+2, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius} , {minRadius}+4, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius}+2, {minRadius}+6, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius}+6, {minRadius}+10, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius}+10, {minRadius}+18, True, False) -newSample.addStack({zMin}, {zMax}, {minRadius}+18, {minRadius}+24, True, False) - """ - for (iFile, redshift) in enumerate(redshifts): - fileNum = fileNums[iFile] - zBox = float(redshift) - - zMin = zBox - zMax = zBox + 0.1 - zMin, zMax = getMaskedZRange(lbox, zBox, zMin, zMax, omegaM) - dataFileName = dataFileNameBase + fileNum + suffix - - sampleName = getSampleName(setName, redshift, useVel, - iSlice=-1, iVol=-1) - - scriptFile.write(sampleInfo.format(dataFile=dataFileName, - dataFormat=dataFormat, - sampleName=sampleName, - maskFileName=maskFileName, - zMin=zMin, - zMax=zMax, - usePecVel=useVel, - minRadius=minRadius)) - - scriptFile.close() - return - - -#------------------------------------------------------------------------------ -# now place these particles on a lightcone, restrict redshift range, apply mask -def applyMask(inFileName, outFileName, maskFileName, lbox, zBox, zMin, zMax, omegaM): - mask = hp.read_map(maskFileName) - nside = hp.get_nside(mask) - - inFile = open(inFileName, 'r') - outFile = open(outFileName, 'w') - - Om = float(omegaM) - - # converter from redshift to comoving distance - zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000) - zVsDX = np.zeros(len(zVsDY)) - for i in xrange(len(zVsDY)): - zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om) - - zBoxInMpc = vp.angularDiameter(zBox, Om=Om) - - for (iLine,line) in enumerate(inFile): - if iLine < 5: - print >> outFile, line.rstrip() - continue - - line = line.split(' ') - uniqueID = int(line[0]) - x = float(line[1]) - lbox/2. - y = float(line[2]) - lbox/2. - z = float(line[3]) - lbox/2. - vz = float(line[4]) - - # TODO if usePecvel - - redshift = np.sqrt(x*x + y*y + z*z) - redshift = np.interp(zBoxInMpc+100./LIGHT_SPEED*redshift, zVsDX, zVsDY) - - if redshift < zMin or redshift > zMax: continue - - vec = np.array((x,y,z)) - theta, phi = hp.vec2ang(vec) - theta = theta[0] - phi = phi[0] - RA = phi*180./np.pi - Dec = 180./np.pi*(np.pi/2.-theta) - - pix = hp.vec2pix(nside, x, y, z) - if mask[pix] <= 0.2: continue - - print >> outFile, RA, Dec, redshift*LIGHT_SPEED, uniqueID, x, y, z - - inFile.close() - outFile.close() - -#------------------------------------------------------------------------------ -def getMaskedZRange(lbox, zBox, zMin, zMax, omegaM): - - if zMin < zBox: zMin = zBox - Om = float(omegaM) - - # converter from redshift to comoving distance - zVsDY = np.linspace(0., zBox+8*lbox*100./LIGHT_SPEED, 10000) - zVsDX = np.zeros(len(zVsDY)) - for i in xrange(len(zVsDY)): - zVsDX[i] = vp.angularDiameter(zVsDY[i], Om=Om) - - zBoxInMpc = vp.angularDiameter(zBox, Om=Om) - - boxWidthZ = np.interp(vp.angularDiameter(zBox,Om=Om)+100. / \ - LIGHT_SPEED*lbox/2., zVsDX, zVsDY)-zBox - - print "RANGE", zBox, zBox+boxWidthZ - if zMax > zBox+boxWidthZ: zMax = zBox+boxWidthZ - - return zMin, zMax - -#------------------------------------------------------------------------------ -#------------------------------------------------------------------------------ - -inCatalogDir = catalogDir - -prefix = "masked_" + prefix -voidOutputDir += "masked/" -figDir += "masked/" -logDir += "masked/" -scriptDir += "masked/" -catalogDir += "masked/" -dataType = "observation" - -if not os.access(scriptDir, os.F_OK): os.mkdir(scriptDir) -if not os.access(catalogDir, os.F_OK): os.mkdir(catalogDir) - -#------------------------------------------------------------------------------ -# first the directly downsampled runs -# Note: ss0.002 ~ SDSS DR7 dim2 -# ss0.0004 ~ SDSS DR9 mid -baseResolution = float(numPart)/lbox/lbox/lbox # particles/Mpc^3 -for thisSubSample in subSamples: - - keepFraction = float(thisSubSample) / baseResolution - maxKeep = keepFraction * numPart - minRadius = int(np.ceil(lbox/maxKeep**(1./3))) - - if args.script or args.all: - print " Doing subsample", thisSubSample, " scripts" - setName = prefix+"ss"+str(thisSubSample) - writeObservationScript(setName, "masked_md.ss"+str(thisSubSample)+"_z", - maskFileName, - scriptDir, catalogDir, fileNums, - redshifts, False, minRadius, omegaM) - writeObservationScript(setName, "masked_md.ss"+str(thisSubSample)+"_z", - maskFileName, - scriptDir, catalogDir, fileNums, - redshifts, True, minRadius, omegaM) - - if args.subsample or args.all: - print " Doing subsample", thisSubSample - - for (iRedshift, redshift) in enumerate(redshifts): - print " redshift", redshift - - zMin = float(redshift) - zMax = zMin + 0.1 - zMin, zMax = getMaskedZRange(lbox, float(redshift), zMin, zMax, omegaM) - - if dataFormat == "multidark" or dataFormat == "random": - inFileName = inCatalogDir+"/md.ss"+str(thisSubSample)+"_z"+redshift+".dat" - outFileName = catalogDir+"/masked_md.ss"+str(thisSubSample)+"_z"+redshift+".dat" - - applyMask(inFileName, outFileName, maskFileName, lbox, float(redshift), - zMin, zMax, omegaM) - - -# ----------------------------------------------------------------------------- -# now halos -if (args.script or args.all) and dataFormat == "multidark": - print " Doing halo scripts" - - for minHaloMass in minHaloMasses: - - if dataFormat == "multidark": - setName = prefix+"halos_min"+str(minHaloMass) - writeScript(setName, "md.halos_min"+str(minHaloMass)+"_z", - scriptDir, catalogDir, fileNums, - redshifts, - numSubvolumes, numSlices, False, lbox, minRadius, omegaM) - writeScript(setName, "md.halos_min"+str(minHaloMass)+"_z", - scriptDir, catalogDir, fileNums, - redshifts, - numSubvolumes, numSlices, True, lbox, minRadius, omegaM) - -if args.halos or args.all: - print " Doing halos" - - for minHaloMass in minHaloMasses: - print " min halo mass = ", minHaloMass - - for (iRedshift, redshift) in enumerate(redshifts): - print " z = ", redshift - - dataFile = catalogDir+haloFileBase+fileNums[iRedshift] - inFile = open(dataFile, 'r') - numPart = 0 - for line in inFile: - line = line.split(',') - if minHaloMass == "none" or float(line[6]) > minHaloMass: - numPart += 1 - inFile.close() - - sampleName = "md.halos_min"+str(minHaloMass)+"_z"+redshift - outFile = open(catalogDir+"/"+sampleName+".dat", 'w') - -# ----------------------------------------------------------------------------- -# now the SDSS HOD - -if (args.script or args.all) and dataFormat == "multidark": - print " Doing DR7 HOD scripts" - if dataFormat == "multidark": - setName = prefix+"hod_dr72dim2" - writeScript(setName, "md.hod_dr72dim2_z", - scriptDir, catalogDir, fileNums, redshifts, - numSubvolumes, numSlices, False, lbox, 5, omegaM) - writeScript(setName, "md.hod_dr72dim2_z", - scriptDir, catalogDir, fileNums, redshifts, - numSubvolumes, numSlices, True, lbox, 5, omegaM) - -if args.hod or args.all: - print " Doing DR7 HOD" - for (iRedshift, redshift) in enumerate(redshifts): - print " z = ", redshift - - sampleName = getSampleName("md.hod_dr72dim2", redshift, False) - outFileName = catalogDir+"/"+sampleName+".dat" - -# ----------------------------------------------------------------------------- -# now the BOSS HOD -if (args.script or args.all) and dataFormat == "multidark": - print " Doing DR9 HOD scripts" - if dataFormat == "multidark": - setName = prefix+"hod_dr9mid" - writeScript(setName, "md.hod_dr9mid_z", - scriptDir, catalogDir, fileNums, redshifts, - numSubvolumes, numSlices, False, lbox, 15, omegaM) - writeScript(setName, "md.hod_dr9mid_z", - scriptDir, catalogDir, fileNums, redshifts, - numSubvolumes, numSlices, True, lbox, 15, omegaM) - -if args.hod or args.all: - print " Doing DR9 HOD" - for (iRedshift, redshift) in enumerate(redshifts): - print " z = ", redshift - - outFileName = catalogDir+"/"+sampleName+".dat" - - if dataFormat == "multidark" or dataFormat == "random": - sampleName = getSampleName("md.hod_dr9mid", redshift, False) - inFileName = inCatalogDir+"/"+sampleName+".dat" - outFileName = catalogDir+"/masked_"+sampleName+".dat" - - zMin = float(redshift) - zMax = zMin + 0.1 - zMin, zMax = getMaskedZRange(lbox, float(redshift), zMin, zMax, omegaM) - - applyMask(inFileName, outFileName, maskFileName, lbox, float(redshift), - zMin, zMax, omegaM) diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index ca300bc..551f298 100755 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/python_tools/pipeline_source/prepareCatalogs.in.py @@ -236,7 +236,7 @@ for thisSubSample in subSamples: fileNums, redshifts, numSubvolumes, numSlices, True, lbox, minRadius, omegaM, subsample=1.0) - elif dataFormat == "gadget": + elif dataFormat == "gadget" or dataFormat == "lanl": writeScript(setName, particleFileBase, scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, False, lbox, minRadius, omegaM, @@ -311,7 +311,7 @@ for thisSubSample in subSamples: # ----------------------------------------------------------------------------- # now halos -if (args.script or args.all) and dataFormat == "multidark": +if (args.script or args.all) and (dataFormat == "multidark" or dataFormat == "lanl"): print " Doing halo scripts" for minHaloMass in minHaloMasses: @@ -319,10 +319,16 @@ if (args.script or args.all) and dataFormat == "multidark": dataFile = catalogDir+haloFileBase+fileNums[0] inFile = open(dataFile, 'r') numPart = 0 - for line in inFile: - line = line.split(',') - if minHaloMass == "none" or float(line[6]) > minHaloMass: - numPart += 1 + if dataFormat == "multidark": + for line in inFile: + line = line.split(',') + if minHaloMass == "none" or float(line[6]) > minHaloMass: + numPart += 1 + elif dataFormat == "lanl": + for line in inFile: + line = line.split(' ') + if minHaloMass == "none" or float(line[0]) > minHaloMass: + numPart += 1 inFile.close() minRadius = 2*int(np.ceil(lbox/numPart**(1./3.))) @@ -349,10 +355,15 @@ if args.halos or args.all: dataFile = catalogDir+haloFileBase+fileNums[iRedshift] inFile = open(dataFile, 'r') numPart = 0 - for line in inFile: - line = line.split(',') - if minHaloMass == "none" or float(line[6]) > minHaloMass: - numPart += 1 + if dataFormat == "multidark": + for line in inFile: + line = line.split(',') + if minHaloMass == "none" or float(line[6]) > minHaloMass: + numPart += 1 + elif dataFormat == "lanl": + line = line.split(' ') + if minHaloMass == "none" or float(line[0]) > minHaloMass: + numPart += 1 inFile.close() sampleName = prefix+"halos_min"+str(minHaloMass)+"_z"+redshift @@ -366,12 +377,20 @@ if args.halos or args.all: inFile = open(dataFile, 'r') for (iHalo,line) in enumerate(inFile): - line = line.split(',') - if minHaloMass == "none" or float(line[6]) > minHaloMass: - x = float(line[0]) - y = float(line[1]) - z = float(line[2]) - vz = float(line[5]) + if dataFormat == "multidark": + line = line.split(',') + if minHaloMass == "none" or float(line[6]) > minHaloMass: + x = float(line[0]) + y = float(line[1]) + z = float(line[2]) + vz = float(line[5]) + elif dataFormat == "lanl": + line = line.split(' ') + if minHaloMass == "none" or float(line[0]) > minHaloMass: + x = float(line[1]) + y = float(line[2]) + z = float(line[3]) + vz = float(line[4]) # write to output file outFile.write("%d %e %e %e %e\n" %(iHalo,x,y,z,vz)) @@ -422,14 +441,13 @@ BOX_SIZE {boxSize} root_filename hod """ -if (args.script or args.all) and dataFormat == "multidark": +if (args.script or args.all) and (dataFormat == "multidark" or dataFormat == "lanl"): print " Doing DR7 HOD scripts" - if dataFormat == "multidark": - setName = prefix+"hod_dr72dim2" - writeScript(setName, "md.hod_dr72dim2_z", + setName = prefix+"hod_dr72dim2" + writeScript(setName, prefix+"hod_dr72dim2_z", scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, False, lbox, 5, omegaM) - writeScript(setName, "md.hod_dr72dim2_z", + writeScript(setName, prefix+"hod_dr72dim2_z", scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, True, lbox, 5, omegaM) @@ -440,7 +458,7 @@ if args.hod or args.all: parFileName = "./hod.par" parFile = open(parFileName, 'w') - haloFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift] + haloFile = catalogDir+haloFileBase+fileNums[iRedshift] parFile.write(parFileText.format(omegaM=omegaM, hubble=hubble, redshift=redshift, @@ -456,7 +474,7 @@ if args.hod or args.all: os.system(hodPath+" "+parFileName+">& /dev/null") - sampleName = getSampleName("md.hod_dr72dim2", redshift, False) + sampleName = getSampleName(prefix+"hod_dr72dim2", redshift, False) outFileName = catalogDir+"/"+sampleName+".dat" os.system("mv hod.mock" + " " + outFileName) @@ -464,14 +482,13 @@ if args.hod or args.all: # ----------------------------------------------------------------------------- # now the BOSS HOD -if (args.script or args.all) and dataFormat == "multidark": +if (args.script or args.all) and (dataFormat == "multidark" or dataFormat == "lanl"): print " Doing DR9 HOD scripts" - if dataFormat == "multidark": - setName = prefix+"hod_dr9mid" - writeScript(setName, "md.hod_dr9mid_z", + setName = prefix+"hod_dr9mid" + writeScript(setName, prefix+"hod_dr9mid_z", scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, False, lbox, 15, omegaM) - writeScript(setName, "md.hod_dr9mid_z", + writeScript(setName, prefix+"hod_dr9mid_z", scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, True, lbox, 15, omegaM) @@ -483,7 +500,7 @@ if args.hod or args.all: # these parameters come from Manera et al 2012, eq. 26 parFileName = "./hod.par" parFile = open(parFileName, 'w') - haloFile = catalogDir+"/mdr1_halos_z"+fileNums[iRedshift] + haloFile = catalogDir+haloFileBase+fileNums[iRedshift] parFile.write(parFileText.format(omegaM=omegaM, hubble=hubble, redshift=redshift, @@ -499,7 +516,7 @@ if args.hod or args.all: os.system(hodPath+" "+parFileName+">& /dev/null") - sampleName = getSampleName("md.hod_dr9mid", redshift, False) + sampleName = getSampleName(prefix+"hod_dr9mid", redshift, False) outFileName = catalogDir+"/"+sampleName+".dat" os.system("mv hod.mock" + " " + outFileName)