added density thresholding for halos and galaxies in preparation scripts

python_tools/pipeline_source/defaults.py

@@ -111,6 +111,9 @@ haloFileDummy = 'NNNNN'
 #   use "none" to get all halos
 minHaloMasses = [1.2e13]
 
+# density threshold for halo catalogs
+haloDenList = [0.0002]
+
 # locations of data in the halo catalog
 haloFileMCol  = 6
 haloFileXCol  = 0
@@ -135,7 +138,6 @@ lbox = 999.983 # Mpc/h
 omegaM = 0.2847979853038958
 hubble = 0.6962
 
-#galDens = 0.000225
 hodParmList = [
   {'name'       : "dr9mid", #BOSS: Manera et al. 2012, eq. 26
    'Mmin'       : 0.0,
@@ -143,7 +145,8 @@ hodParmList = [
    'sigma_logM' : 0.596,
    'alpha'      : 1.0127,
    'Mcut'       : 1.19399e13,
-   'galDens'    : 0.0002,
+   'galDens'    : 0.0002, # density passed to HOD code
+   'galDensFinal'    : 0.0002, # subsample galaxies to reach this density
   },
 ]
 

python_tools/pipeline_source/prepareCatalogs.in.py

@@ -549,7 +549,7 @@ for iSubSample in xrange(len(subSamples)):
   prevSubSample = thisSubSample
 
 # -----------------------------------------------------------------------------
-# now halos
+# now halos - filter on mass
 if (args.script or args.all) and haloFileBase != "":
 
   for minHaloMass in minHaloMasses:
@@ -612,8 +612,8 @@ if (args.script or args.all) and haloFileBase != "":
                   numSubvolumes, numSlices, True, lbox, minRadius, omegaM,
                   dataFileNameList = fileList)
 
-if (args.halos or args.all) and haloFileBase != "":
+if (args.halos or args.all) and haloFileBase != "" and len(minHaloMasses) > 0:
-  print " Doing halos"
+  print " Doing halos - mass"
   sys.stdout.flush()
 
   for minHaloMass in minHaloMasses:
@@ -664,11 +664,20 @@ if (args.halos or args.all) and haloFileBase != "":
       outFile.close()
 
       if dataFormat == "sdf":
+        idTag = "id"
+        for line in open(dataFile):
+          if "int64_t id" in line: 
+            idTag = "id"
+            break
+          if "int64_t ident" in line: 
+            idTag = "ident"
+            break
+
         SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64"
         if minHaloMass == "none":
-          command = "%s -a 200000 %s mass ident x y z vz vy vx parent_id | awk '{if ($9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, outFileName )
+          command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | awk '{if ($9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, idTag, dataFile, outFileName )
         else:
-          command = "%s -a 200000 %s mass ident x y z vz vy vx parent_id | awk '{if ($1>%g && $9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, minHaloMass, outFileName )
+          command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | awk '{if ($1>%g && $9==-1) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, idTag, dataFile, minHaloMass, outFileName )
         #os.system(command)
         subprocess.call(command, shell=True)
         outFile = open(outFileName, 'a')
@@ -697,6 +706,154 @@ if (args.halos or args.all) and haloFileBase != "":
         outFile.close()
         inFile.close()
 
+# -----------------------------------------------------------------------------
+# now halos - filter on density
+if (args.script or args.all) and haloFileBase != "":
+
+  for haloDen in haloDenList:
+    print " Doing halo script", haloDen
+    sys.stdout.flush()
+
+    # estimate number of halos to get density
+    if haloFileDummy == '':
+      dataFile = catalogDir+haloFileBase+fileNums[0]
+    else:
+      dataFile = catalogDir+haloFileBase.replace(haloFileDummy, 
+                                                 fileNums[0])
+
+    numPart = haloDen * lbox**3
+    meanDen = lbox/numPart**(1./3.)
+    if meanDen < 1: 
+      minRadius = 4
+    else:
+      minRadius = int(np.ceil(meanDen))
+  
+    setName = prefix+"halos_den"+str(haloDen)
+    fileList = []
+    for (iRedshift, redshift) in enumerate(redshifts):
+      sampleName = getSampleName(setName, redshift, False)
+      outFileName = sampleName+".dat"
+      fileList.append(outFileName) 
+
+    writeScript(setName, prefix+"halos_den"+str(haloDen)+"_z", "multidark",
+                scriptDir, catalogDir, fileNums, 
+                redshifts, 
+                numSubvolumes, numSlices, False, lbox, minRadius, omegaM,
+                dataFileNameList = fileList)
+   
+    if doPecVel:
+      writeScript(setName, prefix+"halos_den"+str(haloDen)+"_z", 
+                  "multidark",
+                  scriptDir, catalogDir, fileNums, 
+                  redshifts, 
+                  numSubvolumes, numSlices, True, lbox, minRadius, omegaM,
+                  dataFileNameList = fileList)
+
+if (args.halos or args.all) and haloFileBase != "":
+  print " Doing halos - density"
+  sys.stdout.flush()
+
+  for haloDen in haloDenList:
+    print "  halo density = ", haloDen
+    sys.stdout.flush()
+
+    for (iRedshift, redshift) in enumerate(redshifts):
+      print "   z = ", redshift
+      sys.stdout.flush()
+
+      if haloFileDummy == '':
+        dataFile = catalogDir+haloFileBase+fileNums[iRedshift]
+      else:
+        dataFile = catalogDir+haloFileBase.replace(haloFileDummy, 
+                                                   fileNums[iRedshift])
+      # check the actual density of halos
+      inFile = open(dataFile, 'r')
+      numPart = 0
+      if dataFormat == "sdf":
+        for line in inFile:
+          if "nhalos" in line:
+            numPart = int(line.split()[3].strip(';'))
+            break
+          if "npart" in line:
+            numPart = int(line.split()[3].strip(';'))
+            break
+        inFile.close()
+      else:
+        for (iLine, line) in enumerate(inFile):
+          if iLine < haloFileNumComLines: continue
+          numPart += 1
+        inFile.close()
+
+      numPartExpect = int(np.ceil(haloDen * lbox**3))
+      if numPart < numPartExpect:
+        print " ERROR: not enough halos to support that density! Maximum is %g" % (1.*numPart / lbox**3)
+        exit(-1)
+
+      sampleName = prefix+"halos_den"+str(haloDen)+"_z"+redshifts[iRedshift]
+      outFileName = catalogDir+"/"+sampleName+".dat"
+      outFile = open(outFileName, 'w')
+      outFile.write("%f\n" %(lbox))
+      outFile.write("%s\n" %(omegaM))
+      outFile.write("%s\n" %(hubble))
+      outFile.write("%s\n" %(redshift))
+      outFile.write("%d\n" %(numPartExpect))
+      outFile.close()
+
+      if dataFormat == "sdf":
+        idTag = "id"
+        for line in open(dataFile):
+          if "int64_t id" in line: 
+            idTag = "id"
+            break
+          if "int64_t ident" in line: 
+            idTag = "ident"
+            break
+
+        SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64"
+        command = "%s -a 200000 %s mass %s x y z vz vy vx parent_id | sort -rg -k 1 | awk '{if ($9==-1 && NR < %d) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, idTag, numPartExpect, outFileName )
+        #os.system(command)
+        subprocess.call(command, shell=True)
+        outFile = open(outFileName, 'a')
+        outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n")
+        outFile.close()
+      else:
+        inFile = open(dataFile, 'r')
+        haloList = []
+        for (iHalo,line) in enumerate(inFile):
+          if iHalo < haloFileNumComLines: continue
+          line = line.split(haloFileColSep)
+          x  = float(line[haloFileXCol]) * haloFilePosRescale
+          y  = float(line[haloFileYCol]) * haloFilePosRescale
+          z  = float(line[haloFileZCol]) * haloFilePosRescale
+          vz = float(line[haloFileVZCol])
+          vy = float(line[haloFileVYCol])
+          vx = float(line[haloFileVXCol])
+          mass = float(line[haloFileMCol])
+          haloList.append((x,y,z,vz,vy,vx,mass))
+        inFile.close()
+
+        haloList = np.array(haloList)
+        haloMasses = haloList[:,6]
+        haloList = haloList[np.argsort(haloMasses)[::-1]]
+      
+        outFile = open(outFileName, 'a')
+        for iHalo in xrange(numPartExpect):
+            x = haloList[iHalo,0]
+            y = haloList[iHalo,1]
+            z = haloList[iHalo,2]
+            vz = haloList[iHalo,3]
+            vy = haloList[iHalo,4]
+            vx = haloList[iHalo,5]
+            mass = haloList[iHalo,6]
+
+            # write to output file
+            outFile.write("%d %e %e %e %e %e %e %e\n" %(iHalo,x,y,z,
+                                                        vz,vy,vx,mass))
+
+        outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n")
+        outFile.close()
+
+
 # -----------------------------------------------------------------------------
 # now the HOD
 parFileText = """
@@ -754,26 +911,12 @@ if (args.script or args.all) and haloFileBase != "":
     print "  ", thisHod['name']
 
     # estimate number of halos to get density
-    if haloFileDummy == '':
+    numPart = thisHod['galDensFinal'] * lbox**3
-      dataFile = catalogDir+haloFileBase+fileNums[0]
+    meanDen = lbox/numPart**(1./3.)
+    if meanDen < 1: 
+      minRadius = 4
     else:
-      dataFile = catalogDir+haloFileBase.replace(haloFileDummy,
+      minRadius = int(np.ceil(meanDen))
-                                                 fileNums[0])
-    if dataFormat == "sdf":
-      numPart = 0
-      SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64"
-      command = "%s -a 200000 %s mass parent_id | awk '{if ($1>%g && $2==-1) print $1}' | wc" % (SDFcvt_PATH, dataFile, thisHod['Mcut'])
-      numPart = subprocess.check_output(command, shell=True)
-      numPart = int(numPart.split()[0])
-    else:
-      inFile = open(dataFile, 'r')
-      for (iHalo,line) in enumerate(inFile):
-        if iHalo < haloFileNumComLines: continue
-        line = line.split(haloFileColSep)
-        if float(line[haloFileMCol]) > thisHod['Mcut']: numPart += 1
-      inFile.close()
-
-    minRadius = int(np.ceil(float(lbox)/numPart**(1./3.)))
   
     setName = prefix+"hod_"+thisHod['name']
     writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark",
@@ -824,6 +967,7 @@ if (args.hod or args.all) and haloFileBase != "":
                                      sigma_logM=thisHod['sigma_logM'],
                                      alpha=thisHod['alpha'],
                                      Mcut=thisHod['Mcut'],
+                                     #galden=-1,
                                      galden=thisHod['galDens'],
                                      haloFile=haloFile,
                                      haloFileFormat=dataFormat,
@@ -857,9 +1001,46 @@ if (args.hod or args.all) and haloFileBase != "":
           print line
         exit(-1)
 
-      outFileName = catalogDir+"/"+sampleName+".dat"
+      # now randomly subsample the galaxies to get desired density
-      os.system("mv %s/hod_%s.mock %s" % (catalogDir, sampleName, outFileName))
+      inFile = open(catalogDir+"/hod_"+sampleName+".mock")
-      os.system("rm %s/hod.*" % catalogDir)
+      numPartExpect = thisHod['galDensFinal'] * lbox**3
+      numPartActual = 0
+      actualDen = 0
+      for (iLine,line) in enumerate(inFile):
+        numPartActual += 1
+      inFile.close()
+
+      if numPartActual < numPartExpect:
+        print " ERROR: not enough halos to support that density! Maximum is %g" % (1.*numPartActual / lbox**3)
+        exit(-1)
+
+      actualDen = 1.*numPartActual / lbox**3
+      keepFraction = float(thisHod['galDensFinal']) / actualDen
+
+      inFile = open(catalogDir+"/hod_"+sampleName+".mock")
+      outFile = open(catalogDir+"/"+sampleName+".dat", 'w')
+
+      outFile.write("%f\n" %(lbox))
+      outFile.write("%s\n" %(omegaM))
+      outFile.write("%s\n" %(hubble))
+      outFile.write("%s\n" %(redshift))
+      outFile.write("%d\n" %(numPartExpect))
+
+      numKept = 0
+      for (iLine,line) in enumerate(inFile):
+        if np.random.uniform() > keepFraction: continue
+        outFile.write(line) 
+        numKept += 1
+
+      inFile.close()
+      os.unlink(catalogDir+"/hod_"+sampleName+".mock") 
+
+      outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n")
+      outFile.close()
+      
+
+      #os.system("mv %s/hod_%s.mock %s" % (catalogDir, sampleName, outFileName))
+      #os.system("rm %s/hod.*" % catalogDir)
       os.system("rm ./hod.par")
       os.system("rm ./hod-usedvalues")