consistent handling of comoving coords in sims

c_tools/stacking/pruneVoids.cpp

@@ -163,7 +163,7 @@ int main(int argc, char **argv) {
   gsl_function expanF;
   expanF.function = &expanFun;
   struct my_expan_params expanParams;
-  expanParams.Om = 0.27;
+  expanParams.Om = args.omegaM_arg;
   expanParams.w0 = -1.0;
   expanParams.wa = 0.0;
   expanF.params = &expanParams;
@@ -741,7 +741,12 @@ int main(int argc, char **argv) {
     } else {
 
       voids[iVoid].redshiftInMpc = voids[iVoid].barycenter[2];
+      if (args.useComoving_flag) {
+        voids[iVoid].redshift = gsl_interp_eval(interp, dL, redshifts,
+                 voids[iVoid].redshiftInMpc, acc);
+      } else {
         voids[iVoid].redshift = voids[iVoid].barycenter[2]/LIGHT_SPEED*100.;
+      }
 
       nearestEdge = 1.e99;
      

c_tools/stacking/pruneVoids.ggo

@@ -25,6 +25,8 @@ option "isObservation" - "We are working with observational data" flag off
 
 option "useComoving" - "Void positions are in comoving coordinates" flag off
 
+option "omegaM" - "Omega_M for redshift convertion" double optional default="0.27"
+
 option "zMin" - "Minimum redshift of sample" double optional default="0.0"
 option "zMax" - "Maximum redshift of sample" double optional default="10.0"
 

python_tools/pipeline_source/prepareCatalogs.in.py

@@ -230,7 +230,8 @@ newSample = Sample(dataFile = "{dataFile}",
                    profileBinSize = "auto",
                    includeInHubble = True,
                    partOfCombo = False,
-                   {autoStack}
+                   {autoStack},
+                   numAPSlices = {numAPSlices},
                    isCombo = False,
                    boxLen = {boxLen},
                    usePecVel = {usePecVel},
@@ -372,6 +373,7 @@ newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, Fal
                                          omegaM=Om,
                                          boxLen=lbox,
                                          autoStack=autoStack,
+                                         numAPSlices=numAPSlices,
                                          usePecVel=useVel,
                                          minRadius=minRadius,
                                          numSubvolumes=numSubvolumes,

python_tools/void_python_tools/backend/classes.py

@@ -76,11 +76,11 @@ class Sample:
   profileBinSize = 2 # Mpc
   autoNumInStack = -1 # set to >0 to automatically generate stacks of size N
   autoPartInStack = -1 # set to >0 to automatically generate stacks with N particles
+  numAPSlices = 1
   volumeLimited = True
   includeInHubble = True
   partOfCombo  = False
   isCombo = False
-  useComoving = False # if True, convert population to comoving coordinates
   
   comboList = []
 
@@ -88,7 +88,6 @@ class Sample:
   boxLen = 1024 # Mpc/h
   usePecVel = False
   subsample = 1.0
-  useLightCone = True
   numSubvolumes = 1
   mySubvolume = 1
 
@@ -98,13 +97,14 @@ class Sample:
                nickName="", maskFile="", selFunFile="",
                zBoundary=(), zRange=(), zBoundaryMpc=(),
                minVoidRadius=0, fakeDensity=0.01, volumeLimited=True,
+               numAPSlices=1,
                includeInHubble=True, partOfCombo=False, isCombo=False, 
                comboList=(), profileBinSize=2.0, skyFraction=0.19,
                boxLen=1024, usePecVel=False, omegaM=0.27, 
                numSubvolumes=1, mySubvolume=1, dataFormat="sdss",
-               useComoving=False,
+               useComoving=True,
                dataType="observation",
-               subsample=1.0, useLightCone=True, autoNumInStack=-1,
+               subsample=1.0, useLightCone=False, autoNumInStack=-1,
                autoPartInStack=-1):
     self.dataFile = dataFile
     self.fullName = fullName
@@ -121,6 +121,7 @@ class Sample:
     self.partOfCombo = partOfCombo 
     self.isCombo = isCombo
     self.comboList = comboList
+    self.numAPSlices = numAPSlices
     self.zobovDir = None
     self.profileBinSize = profileBinSize
     self.skyFraction = skyFraction

python_tools/void_python_tools/backend/launchers.py

@@ -459,6 +459,7 @@ def launchPrune(sample, binPath,
     cmd += observationLine
     cmd += periodicLine
     cmd += useComovingFlag
+    cmd += " --omegaM=" + str(sample.omegaM)
     cmd += " --outputDir=" + zobovDir
     cmd += " --sampleName=" + str(sampleName)
     #cmd += " &> " + logFile
@@ -594,8 +595,8 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None,
   if useComoving or not useLightCone:
     zMin = LIGHT_SPEED/100.*vp.angularDiameter(zMin, Om=omegaM)
     zMax = LIGHT_SPEED/100.*vp.angularDiameter(zMax, Om=omegaM)
-    print min(sample.zRange[1],stack.zMax)*LIGHT_SPEED/100., zMax
+    print min(sample.zRange[1],stack.zMax)*LIGHT_SPEED/100., zMax, sample.zBoundaryMpc[1]
-    print max(sample.zRange[0],stack.zMin)*LIGHT_SPEED/100., zMin
+    print max(sample.zRange[0],stack.zMin)*LIGHT_SPEED/100., zMin, sample.zBoundaryMpc[0]
   else:
     zMin *= LIGHT_SPEED/100.
     zMax *= LIGHT_SPEED/100.
@@ -1185,7 +1186,8 @@ def launchProfile(sample, stack, voidDir=None, logFile=None, continueRun=None,
 
 
 # -----------------------------------------------------------------------------
-def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
+def launchFit(sample, stack, useComoving=None, 
+              logFile=None, voidDir=None, figDir=None,
               continueRun=None, thisDataPortion=None):
 
   sampleName = sample.fullName
@@ -1259,6 +1261,11 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None,
       ret,fits,args = vp.compute_inertia(voidDir, stack.rMax, mode="2d", nBootstraps=500, rMaxInertia=0.7)
       badChain = False
 
+    ## rescale to expected stretch if using comoving coords
+    #if useComoving or not sample.useLightCone:
+    #  exp = args[1]
+    #  expError = args[2]
+
     #np.save(voidDir+"/chain.npy", ret)
     np.savetxt(voidDir+"/fits.out", fits)
 
@@ -1378,11 +1385,8 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
       #   (not not fully accurate) plots
       for (iZBin, stack) in enumerate(sample.getUniqueZBins()):
 
-        if sample.dataType == "observation":
+        aveDist = vp.aveStretch(sample, stack.zMin, stack.zMax, 
-          aveDist = vp.aveStretchCone(stack.zMin, stack.zMax, 
+                                Om=sample.omegaM)
-                                      skyFrac = sample.skyFraction)
-        else:
-          aveDist = vp.aveStretch(stack.zMin, stack.zMax, Om=sample.omegaM)
 
         aveDistList[iZBin, 0] = stack.zMin
         aveDistList[iZBin, 1] = aveDist
@@ -1427,11 +1431,8 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
           #                            skyFrac = sample.skyFraction, 
           #                            voidRedshifts=voidRedshifts)
 
-          if sample.dataType == "observation":
+          aveDist = vp.aveStretch(sample, zBin.zMin, zBin.zMax, 
-            aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+                                  Om=sample.omegaM)
-                                        skyFrac = sample.skyFraction)
-          else:
-            aveDist = vp.aveStretch(zBin.zMin, zBin.zMax, Om=sample.omegaM)
 
           expList[0, iR, iZBin, 2] = aveDist
 
@@ -1439,11 +1440,8 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
             for (iRCheck,rBinCheck) in enumerate(allRBins):
               if zBin.zMin == zBinCheck.zMin and zBin.zMax == zBinCheck.zMax:
                 if rBin.rMin == rBinCheck.rMin and rBin.rMax == rBinCheck.rMax:
-                  if sample.dataType == "observation":
+                  aveDist = vp.aveStretch(sample, zBin.zMin, zBin.zMax,
-                    aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+                                          Om=sample.omegaM)
-                                                skyFrac = sample.skyFraction)
-                  else:
-                    aveDist = vp.aveStretch(zBin.zMin, zBin.zMax)
 
                   allExpList[iSample, iRCheck, iZCheck, 0] = exp
                   allExpList[iSample, iRCheck, iZCheck, 1] = expError
@@ -1470,8 +1468,8 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
         #sys.stderr = open(logFile, 'a')
         plotTitle = "Sample: "+sample.nickName+", "+thisDataPortion+" voids"
         if doPlot:
-          #vp.do_all_obs(zbase, expList, workDir+"/avedistortion_",
+          #vp.do_all_obs(sample, zbase, expList, workDir+"/avedistortion_",
-          vp.do_all_obs(zbase, expList, aveDistList,
+          vp.do_all_obs(sample, zbase, expList, aveDistList,
                         rlist, plotTitle=plotTitle, plotAve=True)
           figure(1).savefig(figDir+"/hubble_"+setName+"_"+sampleName+"_"+\
                             thisDataPortion+\
@@ -1540,11 +1538,8 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
         if zBin.zMaxPlot > plotZmax: plotZmax = zBin.zMaxPlot
         if zBin.zMinPlot < plotZmin: plotZmin = zBin.zMinPlot
 
-        if sample.dataType == "observation":
+        aveDist = vp.aveStretch(sample, zBin.zMin, zBin.zMax,
-          aveDist = vp.aveStretchCone(zBin.zMin, zBin.zMax, 
+                                Om=sample.omegaM)
-                                      skyFrac = sample.skyFraction)
-        else:
-          aveDist = vp.aveStretch(zBin.zMin, zBin.zMax)
 
         aveDistList[iZ, 0] = zBin.zMin
         aveDistList[iZ, 1] = aveDist
@@ -1567,10 +1562,11 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
           #plotTitle = "all samples, "+thisDataPortion+" voids"
           plotTitle = setName
           #plotTitle = ''
-        vp.do_all_obs(zbase, allExpList, aveDistList,
+        vp.do_all_obs(sample, zbase, allExpList, aveDistList,
                       rlist, plotTitle=plotTitle, sampleNames=shortSampleNames,
                       plotAve=True, mulfac = 1.0, biasLine = 1.16, 
-                      plotZmin=plotZmin, plotZmax=plotZmax)
+                      plotZmin=plotZmin, plotZmax=plotZmax,
+                      plotOm1=True)
         figure(1).savefig(figDir+"/hubble_combined_"+setName+"_"+ \
                           thisDataPortion+\
                           ".eps",bbox_inches='tight')
@@ -1588,10 +1584,11 @@ def launchHubble(dataPortions=None, dataSampleList=None, logDir=None,
           #plotTitle = "all samples, "+thisDataPortion+\
           #            " voids (systematics corrected)"
           plotTitle = setName + "(systematics corrected)"
-        vp.do_all_obs(zbase, allExpList, aveDistList,
+        vp.do_all_obs(sample, zbase, allExpList, aveDistList,
                       rlist, plotTitle=plotTitle, sampleNames=shortSampleNames,
                       plotAve=True, mulfac = 1.16, 
-                      plotZmin=plotZmin, plotZmax=plotZmax)
+                      plotZmin=plotZmin, plotZmax=plotZmax,
+                      plotOm1=True)
         figure(1).savefig(figDir+"/hubble_combined_"+setName+"_"+\
                           thisDataPortion+\
                           "_debiased.eps",bbox_inches='tight')