adding in clusters notebook to test likelihood

solike/clusters/survey.py

@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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:
@@ -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

solike/clusters/sz_utils.py

@@ -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):
@@ -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'],
@@ -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