Monday, May 9, 2011

Going back to Working Reconstruction

I spent today trying to reproduce the result from this post when the reconstruction was working on the Mock Data.

I got it working... here is the code (also in ../logs/


import numpy as N
from pylab import *
from matplotlib import *
import crlold as crl
from correlationData import *
import galsphere as gals
import datetime as dt

workingDir = 'run2010211_1354'

# Create file names

oversample = 5.
corrBins = 40.0 #should be one more than you want it to be
mincorr = 0.01
maxcorr = 20.0
convo = 180./pi
corrspace = (maxcorr - mincorr)/corrBins
corrarray = makecorrarray(mincorr, maxcorr, corrspace)

nbins = 40
boxside = 1000.0 #in Mpc/h
xicorrmin = 0.001
xicorrmax = 20.0
nxicorr = 40

#make redshift bins
#boxside = 1000.0 #in Mpc/h
rbinmin = 0 #in Mpc/h
rbinspace = boxside/40 #in Mpc/h
rbinmax = boxside/2.0 + 0.0001 #in Mpc/h
rbins = makeRbins(rbinmin,rbinmax,rbinspace)
nrbins = rbins.size - 1

photoCD = readPhotoCDfile(photoCDfile)

#Create a matrix to input the correlation function data
# Write wps Matrix to file
#The number of redshift bins you want to skip for the
#reconstruction (this removes noisy bins from calculation)
skip = 0
makeWPSmatrix2(workingDir,rbins,corrBins, wpsMatrixFile, skip)

#Create a matrix to input the correlation function data
# Write Xi Matrix to file
makeXimatrix2(workingDir,rbins,nxicorr,xiMatrixFile, skip)

#Plot wps correlation functions
numtoplot = int(nrbins-skip) #number of correlation functions you want to display on plot
numtoplot = 10
plotWPSScreen(skip, numtoplot, workingDir, rbins)

#Plot xi correlation functions
plotXiScreen(skip, numtoplot, workingDir, rbins)

# compute the xi matrix using the spec sample,
npbins=40#the number of bins you want reconstructed
reconstructmax = 0.5 #in Gpc
reconstructionmin = 0.0 #in Gpc

phibins = linspace(start=reconstructionmin,stop=reconstructmax,num=npbins,endpoint=True)
#phibins =(linspace(2,npbins,num=npbins,endpoint=False)+0.5)*(1./npbins)

# invert to recover phi
for i in range(len(phircvd.x)):
if (phircvd.x[i] < 0):

print 'the recovered selection function is not normalized.'

outphi4 = phircvd.x
phibins4 = phibins

start = 1
end = npbins - 1

#junkbins = linspace(start=0,stop=0.5,num=100,endpoint=True)
junkbins = phibins
junk = histogram(photoCD/1000., bins=junkbins)
junkdata = junk[0]
plot(junkbins[start:end], junkdata[start:end], 'r.-', label = 'photometric data')
plot((phibins4[start:end]+1.0*phibsize), outphi4[start:end]*sum(junkdata[start:end])/sum(outphi4[start:end]), 'g-', label = 'phi')

plot(phibins4[start:end],outphi4[start:end]*sum(junkdata[start:end])/sum(outphi4[start:end]),'c.-',label='reconstruction %i'%npbins)
#plot((phibins4[2:]-1.5*phibsize)/1000, phi[2:]*size(photoCD), 'g-', label = 'phi')

xlabel('comoving distance (Gpc/h)')
ylabel('number of objects per bin')
title('Redshift Distribution')
from pylab import *


It involved reverting back to an old version of crl (I saved it as and an older version of the plotting.

Here are the correlation functions:

I am getting the same problems with slight changes to the binning, changing how well the reconstruction works. For instance the best reconstruction is here (with 40 bins):

But if I change to 39 or 41 bins, I get a very different reconstruction:

Which makes me wonder if was just dumb luck that it worked for that particular binning in the first place.

I'm also running with a more interesting photometric distribution function to see how that looks. More to come....

No comments:

Post a Comment