Different Reconstruction

Here is the reconstruction I ran last night with a different photometric distribution function, below shows the photometric and spectroscopic distributions:

I used the same number of bins and angles as yesterdays post. Here are the correlation functions:

And here are the reconstructions with several different binnings. None of them match as well as yesterdays post, but both 40 and 39 are close-ish:

I talked to Alexia this morning, she suggested that I do the binning sensitivity check and also look at how xmpl_rv (in getximat) compares to sximat. I'm doing that now. Most posts to come shortly.

Here's the code to make the following (also in logs/110509log.py):

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 = 'run201159_123'

runConstantsFile = './'+workingDir+'/runConstants.dat'

mockLRGfile,photofile,specfile,corr2dfile,corr3dfile,photo2dfile,photoCDfile,spec2dfile,spec3dfile,\

argumentFile,runConstantsFile,wpsMatrixFile,xiMatrixFile,oversample,corrBins,mincorr,maxcorr,boxside,\

xicorrmin,xicorrmax,nxicorr,rbinmin,rbinmax,nrbins,psubsize,ssubsize=readRunConstantsFileAll(runConstantsFile)

rbinspace = int((rbinmax-rbinmin)/(nrbins))

rbins = makeRbins(rbinmin,rbinmax,rbinspace)

#read in masks

minRaP,maxRaP,minDecP,maxDecP,minCDP,maxCDP,minRaS,maxRaS,minDecS,maxDecS,minCDS,maxCDS=readmasks(workingDir)

#Import the photometric redshift distribution

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

plotWPS(skip, numtoplot, workingDir, rbins)

numtoplot = 10

plotWPSScreen(skip, numtoplot, workingDir, rbins)

#Plot xi correlation functions

plotXi(skip, numtoplot, workingDir, rbins)

plotXiScreen(skip, numtoplot, workingDir, rbins)

#from commandsAlexia import *

xissfile = xiMatrixFile

wpsfile = wpsMatrixFile

nphibins = size(rbins)

# Read in the data

#spectroscopic autocorrelation function

X=N.loadtxt(xissfile,comments='#')

X=X.T

rs=X[0,:]/1000. #want box units (Gpc, not Mpc)

xispec=X[1:,:] #xi in nspecbins

#2D cross correlation between spec and photo samples

th,rlos,wps=N.loadtxt(wpsfile,comments='#',unpack=True)

th=N.unique(th)

rlos=N.unique(rlos)

#rlos=getcomovingD(rlos) #adjusts from redshift to Mpc/h

rlos=rlos/1000 #want box units (Gpc, not Mpc)

#covariance matrix of wps measurements, wps vector arranged

covar = N.identity(wps.shape[0],float)

xierr = N.identity(xispec.shape[1],float)

# compute the xi matrix using the spec sample,

npbins=39#the number of bins you want reconstructed

reconstructmax = 0.5 #in Gpc

reconstructionmin = 0.0 #in Gpc

phibsize = (reconstructmax - reconstructionmin)/npbins

phibins = linspace(start=reconstructionmin,stop=reconstructmax,num=npbins,endpoint=True)

#phibins =(linspace(2,npbins,num=npbins,endpoint=False)+0.5)*(1./npbins)

sximat=crl.getximat(th,rlos,phibins,rs,xispec,xierr)

# invert to recover phi

phircvd=crl.getCoeff(wps,covar,sximat)

for i in range(len(phircvd.x)):

if (phircvd.x[i] < 0):

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='reconstruction %i'%npbins)

xlabel('comoving distance (Gpc/h)')

ylabel('number of objects per bin')

title('Redshift Distribution')

from pylab import *

legend(loc=2)