diff --git a/python_tools/pipeline_source/defaults.py b/python_tools/pipeline_source/defaults.py index 2a3de54..1b9a52d 100644 --- a/python_tools/pipeline_source/defaults.py +++ b/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 }, ] diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index e5e6c87..b6d62b0 100644 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/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 != "": - print " Doing halos" +if (args.halos or args.all) and haloFileBase != "" and len(minHaloMasses) > 0: + 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,27 +911,13 @@ if (args.script or args.all) and haloFileBase != "": print " ", thisHod['name'] # estimate number of halos to get density - if haloFileDummy == '': - dataFile = catalogDir+haloFileBase+fileNums[0] + numPart = thisHod['galDensFinal'] * lbox**3 + meanDen = lbox/numPart**(1./3.) + if meanDen < 1: + minRadius = 4 else: - dataFile = catalogDir+haloFileBase.replace(haloFileDummy, - 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.))) - + minRadius = int(np.ceil(meanDen)) + setName = prefix+"hod_"+thisHod['name'] writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark", scriptDir, catalogDir, fileNums, redshifts, @@ -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" - os.system("mv %s/hod_%s.mock %s" % (catalogDir, sampleName, outFileName)) - os.system("rm %s/hod.*" % catalogDir) + # now randomly subsample the galaxies to get desired density + inFile = open(catalogDir+"/hod_"+sampleName+".mock") + 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")