In [1]:
#Requirements :
# (a) Make sure that Suman's code directory is available
# (b) Make sure that interfacecosmology is available
#Data Files:
# T
#Comparison with Suman's code
#Example with m000 to check that things are working
In [48]:
SumanCode = "/Users/rbiswas/src/newHalopy"
import sys
sys.path.append(SumanCode)
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.ticker as ticker
import matplotlib.pyplot as plt
import camb_utils.cambio as cio
import psutils as psu
import fsigma as fs
import MF as MF
from interfaces import FCPL
import massfunctions as mf
import utils.plotutils as pu
import halomass as hmf
import sys
majorformatter = ticker.ScalarFormatter(useOffset =False)
In [3]:
# INPUT FILES (CAMB outputs)
dirn = "./"
#M000 (LCDM)
tkfile = dirn +"example_data/M000/m000n0_transfer_fin_out.dat"
pkfile = dirn + "example_data/M000/m000n0_matterpower_fin.dat"
#M000n1 (same LCDM, with some fraction of CDM replaced by massive neutrinos)
ntkfile = dirn +"example_data/M000n1/m000n1_transfer_fin_out.dat"
npkfile = dirn +"example_data/M000n1/m000n1_matterpower_fin.dat"
In [4]:
#Cosmological Models:
# Same cosmology represented differently
M000 = FCPL(H0 = 71. ,Om0 = 0.26479, Ob0 = 0.044793, ns = 0.963, As = 2.16e-9)
M000s = FCPL(H0 = 71., Om0 = 0.26479, Ob0 = 0.044793, ns = 0.963, sigma8 = 0.8)
#
M000n1 = FCPL(H0 = 71., Om0 = 0.26479, Ob0 = 0.044793,Neff= 0.046, ns = 0.963, sigma8 = 0.8, sigmamnu=0.94)
In [5]:
#Suman's reading of the file
Tks = cio.loadtransfers(filename = tkfile)
k = Tks[:,0]
Tk = Tks[:,-1]
sumanrhoc = 2.77536627e11
sumannorm = fs.Pk_norm(sigma8 = 0.8, ns = 0.963, k = k , Tk=Tk, hubble= 0.71)
#working power sp
psfrompk = psu.powerspectrum(koverh = None, asciifile=pkfile)
psfromtkcbscb = psu.powerspectrum(koverh = None, pstype ="cb",sigma8ype= "cb", asciifile =tkfile, cosmo = M000s, z =0., method="CAMBoutgrowth")
#psfromtkcbscb = psu.powerspectrum(koverh = None, pstype ="cb", asciifile =tkfile, cosmo = M000s)
In [6]:
#Masses in units of solar masses
Masses = np.logspace(8,16,100)
Rad = np.linspace(.1,20,100)
MassesinMsunOverh = Masses* M000.h
sumansigmar = np.zeros(len(Rad))
sumandlsidlm = np.zeros(len(Masses))
sumandndlnM = np.zeros(len(Masses))
sumandndlnM1 = np.zeros(len(Masses))
sumandndlnM2 = np.zeros(len(Masses))
sumansigmam = np.zeros(len(Masses))
sumanfsigma = np.zeros(len(Masses))
sumanfsigma1 = np.zeros(len(Masses))
sumanfsigma2 = np.zeros(len(Masses))
for i, mass in enumerate(MassesinMsunOverh) :
sumansigmaM = fs.sigmam ( Om = M000.Om0, ns = 0.963, M= mass, N=sumannorm , k = k , Tk = Tk , hubble = 0.71 )
sumanlogsigmaM = fs.logsigm ( Om = M000.Om0 , ns = 0.963, M= mass, N = sumannorm, k = k, Tk = Tk, sigmam = sumansigmaM, hubble = 0.71 )
sumansigmam[i] = sumansigmaM
sumanfsigma[i] = MF.MF_fit(sumansigmaM, z= 0) #mass function as a function of halo mass
#mass function as a function of halo mass
sumandndlnM[i] = sumanfsigma[i] *(sumanrhoc*M000.Om0*sumanlogsigmaM)/mass
sumanfsigma1[i] = MF.MF_fit(sumansigmaM*M000.growth(z = 1.0)[0], z= 1) #mass function as a function of halo mass
#mass function as a function of halo mass
sumandndlnM1[i] = sumanfsigma1[i] *(sumanrhoc*M000.Om0*sumanlogsigmaM)/mass
sumanfsigma2[i] = MF.MF_fit(sumansigmaM*M000.growth(z = 2.0)[0], z= 2) #mass function as a function of halo mass
#mass function as a function of halo mass
sumandndlnM2[i] = sumanfsigma2[i] *(sumanrhoc*M000.Om0*sumanlogsigmaM)/mass
In [7]:
MF.MF_fit(sumansigmaM*M000.growth(z = 1.0)[0], z= 0)
Out[7]:
In [8]:
fig = plt.figure()
ax = plt.subplot()
ax.plot(MassesinMsunOverh, sumanfsigma)
ax3 = ax.twiny()
xt = ax.get_xticks()[1::2]
axtlabs= ["%1.1f" % z for z in (1.0/psu.sigmaM(M = xt/M000.h, ps = psfromtkcbscb,cosmo = M000s))]
ax3.set_xticks(xt)
ax3.set_xticklabels(axtlabs)
#ax2.set_xscale('log')
#axt.set_xscale('log')
ax.set_xscale('log')
ax.grid(True)
#ax2.set_xscale('log')
In [55]:
plt.plot(MassesinMsunOverh, hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s , deltac = 1.686,fittingform = "MICE")/hmf.dndlnM0(M= Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686),'k-',label = "z =0")
plt.plot(MassesinMsunOverh, hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s , deltac = 1.686,fittingform = "MICE", z = 1.0)/hmf.dndlnM0(M= Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z= 1.0),"r--",label ="z = 1")
plt.xscale('log')
plt.legend(loc="best")
plt.ylim(0.9,1.1)
plt.xlabel("Mass ($h^{-1} M_\odot$)")
plt.ylabel("Crocce/Bhattacharya")
plt.grid(True)
In [9]:
#fig = plt.figure()
fig,ax , ax1 = pu.settwopanel(setdifflimits=[0.999,1.001])
#ax = fig.axes
ax.loglog(MassesinMsunOverh, - sumandndlnM, 'k-')
ax.loglog(MassesinMsunOverh, - sumandndlnM1, 'r-')
ax.loglog(MassesinMsunOverh, - sumandndlnM2,'b-')
ax.loglog(MassesinMsunOverh, hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.674, z = 0.)/M000s.h**3, 'k--', lw =4)
#ax.loglog(MassesinMsunOverh, - sumandndlnM1)
ax.loglog(MassesinMsunOverh, hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.674, z =1.0)/M000s.h**3, 'r--', lw =4)
ax.loglog(MassesinMsunOverh, hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.674, z =2.0)/M000s.h**3, 'b--', lw =4)
ax1.plot(MassesinMsunOverh, - sumandndlnM/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.674, z = 0.)*M000s.h**3, "k-", lw =1.0)
ax1.plot(MassesinMsunOverh, - sumandndlnM1/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.684, z = 1.)*M000s.h**3, "r-.", lw =4.0)
ax1.plot(MassesinMsunOverh, - sumandndlnM2/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.686, z = 2.)*M000s.h**3, "b--", lw =4.0)
fig.savefig("Comparison_Sumanvsmine.pdf")
In [26]:
MM = (1.1752e13)/0.71;
xx = hmf.dndlnM0(M = MM, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =2.03625)/0.71**3.0
print xx
In [24]:
ss1, ss2 = 1.0*(psu.sigmaM(M =MM, ps = psfromtkcbscb , cosmo = M000s, deltac =1.686, z = 0.0) *M000s.growth(z= 2.0362537644)[0]),1.0*(psu.sigmaM(M =MM, ps = psfromtkcbscb , cosmo = M000s, deltac =1.686, z = 0.0) *M000s.growth(z= 2.0)[0])
In [27]:
print ss1, xx, MM
In [13]:
print MF.MF_fit(ss1, z= 2) :
In [47]:
plt.plot(Masses*0.71, hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =2.018)/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =2.0))
plt.plot(Masses*0.71, hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =1.006)/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =1.0))
plt.plot(Masses*0.71, hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =0.0)/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =0.0))
plt.xscale('log')
plt.grid(True)
plt.xlabel(r'M ($h^{-1} M_\odot$)')
plt.ylabel("MF(zstep)/MF(zround)", fontsize = 14)
plt.xlim(1.e12,1.5e15)
plt.tight_layout()
plt.savefig("Universalfunction_redshfit.pdf")
In [29]:
rati = hmf.dndlnM0(M = Masses", ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =2.0) /hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac =1.686, z =2.0)
In [22]:
print rati
In [14]:
print mf.__fsigmaBhattacharya(sigma = ss1, deltac = 1.686, z =2.03625377644) ~
In [15]:
hmf.dndlnM0(M = MM , ps = psfromtkcbscb, z =2.03625,cosmo = M000s, deltac=1.686)
Out[15]:
In [16]:
#ax1 = plt.subplot()
plt.plot(MassesinMsunOverh, - sumandndlnM/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.674, z = 0.)*M000s.h**3, "k--")
plt.plot(MassesinMsunOverh, - sumandndlnM1/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.684, z = 1.00)*M000s.h**3, "r--")
plt.plot(MassesinMsunOverh, - sumandndlnM2/hmf.dndlnM0(M = Masses, ps = psfromtkcbscb, cosmo = M000s, deltac=1.686
, z = 2.)l*M000s.h**3, "b--")
#plt.ylim(0.999,1.001)
In [ ]:
sumanfsigma1/ mf.__fsigmaBhattacharya(sigma = M000s.growth(z = 1.0)[0]*psu.sigmaM(M = Masses, ps = psfromtkcbscb,cosmo = M000s),deltac=1.684, z =1.)
In [ ]:
xt =ax.get_xticks()[1::2]
print xt
In [ ]:
axtlabs = "{}}.format(1.1)(1.0/psu.sigmaM(M = xt/M000.h, ps = psfromtkcbscb,cosmo = M000s)
In [ ]:
1.0/psu.sigmaM(M = 10**15, ps = psfromtkcbscb,cosmo = M000s)
In [ ]:
otf = "/Users/rbiswas/doc/projects/NeutrinoMF/data/CAMB_files/cmbM000.tf";
ops = psu.powerspectrum(koverh = None, asciifile=otf, cosmo = M000s)
In [53]:
plt.loglog(psfromtkcbscb[0],psfromtkcbscb[1],'k-')
plt.loglog(ops[0],ops[1],'r--')
plt.loglog(psfrompk[0],psfrompk[1],'g-')
Out[53]:
In [48]:
opsi = psu.powerspectrum(koverh=psfromtkcbscb[0], asciifile= otf, cosmo = M000s)
ax = plt.subplot()
ax.plot(opsi[0],opsi[1]/psfromtkcbscb[1],label = 'wrt CAMB TF')
ax.plot(opsi[0],opsi[1]/ psfrompk[1],'r-',label='wrt CAMB PS')
ax.set_xscale('log')
ax.set_ylim(0.995,1.005)
ax.yaxis.set_major_formatter(majorformatter)
ax.axvline(x=min(ops[0]))
ax.axvline(x=max(ops[0]))
ax.axvline(x=min(psfrompk[0]),color='k')
ax.axvline(x=max(psfrompk[0]),color='k')
ax.legend(loc="best")
ax.grid(True)
In [51]:
plt.plot(ops[0],ops[1])
Out[51]:
In [25]:
yy = hmf.dndlnM0(M = MM, ps = ops, cosmo = M000s, deltac =1.686, z =2.0)/0.71**3.0
print yy
In [21]:
sumandndlnM2/hmf.dndlnM0(M = Masses,ps = psfromtkcbscb, cosmo = M000s, z= 2.0,deltac = 1.686)*M000s.h**3
Out[21]:
In [39]:
M000s.growth(z= 2.0)[0]
Out[39]:
In [41]:
3.3753/3.4152
Out[41]:
In [ ]: