applyMaskToMock now accepts build path

pipeline/applyMaskToMock.py

@@ -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)

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
+    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
+      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):
+        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",
+  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",
+  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)