Skip to content

Commit

Permalink
Merge pull request simonsobs#5 from nbatta/clusters
Browse files Browse the repository at this point in the history 
adding in clusters notebook to test likelihood
  • Loading branch information
@timothydmorton
timothydmorton authored Jan 13, 2020
2 parents 193161f + ed43edc commit 7fa464e
Show file tree
Hide file tree
Showing 3 changed files with 281 additions and 9 deletions.
7 changes: 0 additions & 7 deletions solike/clusters/survey.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@
from astropy.io import fits
import astropy.table as atpy


def read_clust_cat(fitsfile, qmin):
list = fits.open(fitsfile)
data = list[1].data
Expand All @@ -19,7 +18,6 @@ def read_clust_cat(fitsfile, qmin):
ind = np.where(SNR >= qmin)[0]
return z[ind], zerr[ind], Y0[ind], Y0err[ind]


def read_mock_cat(fitsfile, qmin):
list = fits.open(fitsfile)
data = list[1].data
Expand All @@ -31,7 +29,6 @@ def read_mock_cat(fitsfile, qmin):
ind = np.where(SNR >= qmin)[0]
return z[ind], zerr[ind], Y0[ind], Y0err[ind]


def read_matt_mock_cat(fitsfile, qmin):
list = fits.open(fitsfile)
data = list[1].data
Expand All @@ -46,7 +43,6 @@ def read_matt_mock_cat(fitsfile, qmin):
ind = np.where(SNR >= qmin)[0]
return z[ind], zerr[ind], Y0[ind], Y0err[ind]


def loadAreaMask(extName, DIR):
"""Loads the survey area mask (i.e., after edge-trimming and point source masking, produced by nemo).
Returns map array, wcs
Expand All @@ -58,7 +54,6 @@ def loadAreaMask(extName, DIR):

return areaMap, wcs


def loadRMSmap(extName, DIR):
"""Loads the survey RMS map (produced by nemo).
Returns map array, wcs
Expand All @@ -70,7 +65,6 @@ def loadRMSmap(extName, DIR):

return areaMap, wcs


def loadQ(source, tileNames=None):
"""Load the filter mismatch function Q as a dictionary of spline fits.
Args:
Expand Down Expand Up @@ -104,7 +98,6 @@ def loadQ(source, tileNames=None):
tckDict[tileName] = interpolate.splrep(tab["theta500Arcmin"], tab["Q"])
return tckDict


class SurveyData(object):
def __init__(self, nemoOutputDir, ClusterCat, qmin=5.6, szarMock=False, tiles=False):
self.nemodir = nemoOutputDir
Expand Down
147 changes: 145 additions & 2 deletions solike/clusters/sz_utils.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@
import numpy as np
from nemo import signals
from scipy import interpolate
from astropy.cosmology import FlatLambdaCDM

#from nemo import signals

class szutils(object):
def __init__(self,Survey):
Expand All @@ -14,7 +17,7 @@ def P_Yo(self, LgY, M, z, param_vals):

cosmoModel = FlatLambdaCDM(H0=H0, Om0=Om, Ob0=Ob, Tcmb0=2.725)

Ytilde, theta0, Qfilt = signals.y0FromLogM500(
Ytilde, theta0, Qfilt = y0FromLogM500(
np.log10(param_vals['massbias']*M / (H0/100.)),
z,
self.Survey.tckQFit['Q'],
Expand Down Expand Up @@ -69,3 +72,143 @@ def P_of_Y_per(self, LgY, MM, zz, Y_c, Y_err, param_vals):
P_Y = np.nan_to_num(self.P_Yo(LgYa2, MM, zz, param_vals))
ans = np.trapz(P_Y*P_Y_sig, LgY, np.diff(LgY), axis=1) * np.log(10)
return ans


###
'''Routines from nemo (author: Matt Hilton ) to limit dependencies'''

#------------------------------------------------------------------------------------------------------------
def calcR500Mpc(z, M500, cosmoModel):
"""Given z, M500 (in MSun), returns R500 in Mpc, with respect to critical density.
"""

if type(M500) == str:
raise Exception("M500 is a string - check M500MSun in your .yml config file: use, e.g., 1.0e+14 (not 1e14 or 1e+14)")

Ez=cosmoModel.efunc(z)
criticalDensity=cosmoModel.critical_density(z).value
criticalDensity=(criticalDensity*np.power(Mpc_in_cm, 3))/MSun_in_g
R500Mpc=np.power((3*M500)/(4*np.pi*500*criticalDensity), 1.0/3.0)

return R500Mpc

#------------------------------------------------------------------------------------------------------------
def calcTheta500Arcmin(z, M500, cosmoModel):
"""Given z, M500 (in MSun), returns angular size equivalent to R500, with respect to critical density.
"""

R500Mpc=calcR500Mpc(z, M500, cosmoModel)
theta500Arcmin=np.degrees(np.arctan(R500Mpc/cosmoModel.angular_diameter_distance(z).value))*60.0

return theta500Arcmin

#------------------------------------------------------------------------------------------------------------
def calcQ(theta500Arcmin, tck):
"""Returns Q, given theta500Arcmin, and a set of spline fit knots for (theta, Q).
"""

#Q=np.poly1d(coeffs)(theta500Arcmin)
Q=interpolate.splev(theta500Arcmin, tck)

return Q

#------------------------------------------------------------------------------------------------------------
def calcFRel(z, M500, cosmoModel, obsFreqGHz = 148.0):
"""Calculates relativistic correction to SZ effect at specified frequency, given z, M500 in MSun.
This assumes the Arnaud et al. (2005) M-T relation, and applies formulae of Itoh et al. (1998)
As for H13, we return fRel = 1 + delta_SZE (see also Marriage et al. 2011)
"""

# NOTE: we should define constants somewhere else...
h=6.63e-34
kB=1.38e-23
sigmaT=6.6524586e-29
me=9.11e-31
e=1.6e-19
c=3e8

# Using Arnaud et al. (2005) M-T to get temperature
A=3.84e14
B=1.71
#TkeV=5.*np.power(((cosmoModel.efunc(z)*M500)/A), 1/B) # HMF/Astropy
TkeV=5.*np.power(((cosmoModel.Ez(z)*M500)/A), 1/B) # Colossus
TKelvin=TkeV*((1000*e)/kB)

# Itoh et al. (1998) eqns. 2.25 - 2.30
thetae=(kB*TKelvin)/(me*c**2)
X=(h*obsFreqGHz*1e9)/(kB*TCMB)
Xtw=X*(np.cosh(X/2.)/np.sinh(X/2.))
Stw=X/np.sinh(X/2.)

Y0=-4+Xtw

Y1=-10. + (47/2.)*Xtw - (42/5.)*Xtw**2 + (7/10.)*Xtw**3 + np.power(Stw, 2)*(-(21/5.) + (7/5.)*Xtw)

Y2=-(15/2.) + (1023/8.)*Xtw - (868/5.)*Xtw**2 + (329/5.)*Xtw**3 - (44/5.)*Xtw**4 + (11/30.)*Xtw**5 \
+ np.power(Stw, 2)*(-(434/5.) + (658/5.)*Xtw - (242/5.)*Xtw**2 + (143/30.)*Xtw**3) \
+ np.power(Stw, 4)*(-(44/5.) + (187/60.)*Xtw)

Y3=(15/2.) + (2505/8.)*Xtw - (7098/5.)*Xtw**2 + (14253/10.)*Xtw**3 - (18594/35.)*Xtw**4 + (12059/140.)*Xtw**5 - (128/21.)*Xtw**6 + (16/105.)*Xtw**7 \
+ np.power(Stw, 2)*(-(7098/10.) + (14253/5.)*Xtw - (102267/35.)*Xtw**2 + (156767/140.)*Xtw**3 - (1216/7.)*Xtw**4 + (64/7.)*Xtw**5) \
+ np.power(Stw, 4)*(-(18594/35.) + (205003/280.)*Xtw - (1920/7.)*Xtw**2 + (1024/35.)*Xtw**3) \
+ np.power(Stw, 6)*(-(544/21.) + (992/105.)*Xtw)

Y4=-(135/32.) + (30375/128.)*Xtw - (62391/10.)*Xtw**2 + (614727/40.)*Xtw**3 - (124389/10.)*Xtw**4 \
+ (355703/80.)*Xtw**5 - (16568/21.)*Xtw**6 + (7516/105.)*Xtw**7 - (22/7.)*Xtw**8 + (11/210.)*Xtw**9 \
+ np.power(Stw, 2)*(-(62391/20.) + (614727/20.)*Xtw - (1368279/20.)*Xtw**2 + (4624139/80.)*Xtw**3 - (157396/7.)*Xtw**4 \
+ (30064/7.)*Xtw**5 - (2717/7.)*Xtw**6 + (2761/210.)*Xtw**7) \
+ np.power(Stw, 4)*(-(124389/10.) + (6046951/160.)*Xtw - (248520/7.)*Xtw**2 + (481024/35.)*Xtw**3 - (15972/7.)*Xtw**4 + (18689/140.)*Xtw**5) \
+ np.power(Stw, 6)*(-(70414/21.) + (465992/105.)*Xtw - (11792/7.)*Xtw**2 + (19778/105.)*Xtw**3) \
+ np.power(Stw, 8)*(-(682/7.) + (7601/210.)*Xtw)

deltaSZE=((X**3)/(np.exp(X)-1)) * ((thetae*X*np.exp(X))/(np.exp(X)-1)) * (Y0 + Y1*thetae + Y2*thetae**2 + Y3*thetae**3 + Y4*thetae**4)

fRel=1+deltaSZE

return fRel


#------------------------------------------------------------------------------------------------------------
def y0FromLogM500(log10M500, z, tckQFit, tenToA0 = 4.95e-5, B0 = 0.08, Mpivot = 3e14, sigma_int = 0.2,
cosmoModel = None, fRelWeightsDict = {148.0: 1.0}):
"""Predict y0~ given logM500 (in MSun) and redshift. Default scaling relation parameters are A10 (as in
H13).
Use cosmoModel (astropy.cosmology object) to change/specify cosmological parameters.
fRelWeightsDict is used to account for the relativistic correction when y0~ has been constructed
from multi-frequency maps. Weights should sum to 1.0; keys are observed frequency in GHz.
Returns y0~, theta500Arcmin, Q
"""

if type(Mpivot) == str:
raise Exception("Mpivot is a string - check Mpivot in your .yml config file: use, e.g., 3.0e+14 (not 3e14 or 3e+14)")

# Filtering/detection was performed with a fixed fiducial cosmology... so we don't need to recalculate Q
# We just need to recalculate theta500Arcmin and E(z) only
M500=np.power(10, log10M500)
theta500Arcmin=calcTheta500Arcmin(z, M500, cosmoModel)
Q=calcQ(theta500Arcmin, tckQFit)

# Relativistic correction: now a little more complicated, to account for fact y0~ maps are weighted sum
# of individual frequency maps, and relativistic correction size varies with frequency
fRels=[]
freqWeights=[]
for obsFreqGHz in fRelWeightsDict.keys():
fRels.append(calcFRel(z, M500, cosmoModel, obsFreqGHz = obsFreqGHz))
freqWeights.append(fRelWeightsDict[obsFreqGHz])
fRel=np.average(np.array(fRels), axis = 0, weights = freqWeights)

# UPP relation according to H13
# NOTE: m in H13 is M/Mpivot
# NOTE: this goes negative for crazy masses where the Q polynomial fit goes -ve, so ignore those
y0pred=tenToA0*np.power(cosmoModel.efunc(z), 2)*np.power(M500/Mpivot, 1+B0)*Q*fRel

return y0pred, theta500Arcmin, Q
Loading

0 comments on commit 7fa464e

Please sign in to comment.