Merge pull request #1 from llinke1/BaryonicEffects

Added model for baryonic effects from HMCode2020

Author: alexander-mead

hmcode/hmcode.py

@@ -16,7 +16,8 @@
 sigma_cold_approx = False # Should the Eisenstein & Hu (1999) approximation be used for the cold transfer function?
 
 def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata, 
-          Mmin=1e0, Mmax=1e18, nM=256, verbose=False) -> np.ndarray:
+          Mmin=1e0, Mmax=1e18, nM=256, verbose=False, 
+          baryons=False, baryon_params=None, hmcodetweaks=True) -> np.ndarray:
     '''
     Calculates the HMcode matter-matter power spectrum
     Args:
@@ -26,6 +27,9 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
         Mmin: Minimum mass for the halo-model calculation [Msun/h]
         Mmax: Maximum mass for the halo-model calculation [Msun/h]
         nM: Number of mass bins for the halo-model calculation
+        baryons: if true: include baryonic effects from HMCode2020; if false: Dark Matter Only
+        baryon_params: Dictionary containing the 6 baryonic parameters from HMCode2020
+        hmcodetweaks: if true: Use the changes to the vanilla halo model from HMCode2020; if false: Use vanilla halo model
     Returns:
         Array of matter power spectra: Pk[z, k]
     '''
@@ -60,6 +64,16 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
         print('Linear growth at z=0: {:.3}'.format(growth(1.)))
         print()
 
+    if verbose and baryons:
+        print("Using baryonic feedback model from HMCode2020")
+        print("B0: {:.4f}".format(baryon_params["B0"]))
+        print("Bz: {:.4f}".format(baryon_params["Bz"]))
+        print("Mb0: {:.1e}".format(baryon_params["Mb0"]))
+        print("Mbz: {:.4f}".format(baryon_params["Mbz"]))
+        print("f*0: {:.4f}".format(baryon_params["f0"]))
+        print("f*z: {:.4f}".format(baryon_params["fz"]))
+        print()
+
     # Linear power interpolator
     interp, _, k_interp = CAMB_results.get_matter_power_interpolator(nonlinear=False, 
                                                                      return_z_k=True,
@@ -96,6 +110,7 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
             print('Linear collapse threshold: {:.4}'.format(dc))
             print('Virial halo overdensity: {:.4}'.format(Dv))
 
+
         # Initialise halo model
         hmod = halo.model(z, Om_m, name='Sheth & Tormen (1999)', Dv=Dv, dc=dc)
 
@@ -135,9 +150,16 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
         nd = 2.853                   # Two-halo damping power; equation (16)
         ks = 0.05618*sigma8**-1.013  # One-halo damping wavenumber; equation (17)
         eta = 0.1281*sigma8**-0.3644 # Halo bloating parameter; equation (19)
-        B = 5.196                    # Minimum halo concentration; equation (20)
+        B =   5.196                    # Minimum halo concentration; equation (20)
         alpha = 1.875*(1.603)**neff  # Transition smoothing; equation (23)
-        if verbose:
+
+        
+        if hmcodetweaks==False: #Go to vanilla halo model if no hmcode tweaks should be used
+            B=4
+            eta=0
+            alpha=1
+
+        if verbose and hmcodetweaks:
            print('Two-halo damping wavenumber; kd: {:.4} h/Mpc'.format(kd))
            print('Two-halo fractional damping; f: {:.4}'.format(f))
            print('Two-halo damping power; nd: {:.4}'.format(nd))
@@ -146,7 +168,22 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
            print('Minimum halo concentration; B: {:.4}'.format(B))#, B/4.)
            print('Transition smoothing; alpha: {:.4}'.format(alpha))
            print()
-
+        elif verbose:
+            print("Using vanilla halo model")
+
+        if baryons:
+            B=baryon_params["B0"]*np.power(10,z*baryon_params["Bz"])
+            Mb=baryon_params["Mb0"]*np.power(10, z*baryon_params["Mbz"])
+            fstar=baryon_params["f0"]*np.power(10, z*baryon_params["fz"])
+        
+
+        
+        if verbose and baryons:
+            print("Halo concentration parameter B:{:.4f}".format(B))
+            print("Gas loss halo mass parameter Mb:{:1e}".format(Mb))
+            print("Effective halo stellar mass fraction fstar:{:.4f}".format(f))
+
+        
         # Halo concentration
         zf = _get_halo_collapse_redshifts(M, z, iz, dc, growth, CAMB_results, cold=True) # Halo formation redshift
         c = B*(1+zf)/(1.+z)                                                              # Halo concentration; equation (20)
@@ -157,18 +194,30 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
         rv = hmod.virial_radius(M) 
         Uk = np.ones((len(k), len(M)))
         for iM, (_rv, _c, _nu) in enumerate(zip(rv, c, nu)): # TODO: Remove loop for speed?
-            Uk[:, iM] = _win_NFW(k*(_nu**eta), _rv, _c)[:, 0]
-        profile = halo.profile.Fourier(k, M, Uk, amplitude=M*(1.-f_nu)/hmod.rhom, mass_tracer=True) # NOTE: Factor of 1-f_nu
-
+            if baryons:
+                Uk[:, iM] = _win_NFW_baryons(k*(_nu**eta), _rv, _c, M[iM], Mb, fstar, Om_b, Om_m, Om_c, hmod.rhom)[:,0]
+        
+            else:
+                Uk[:, iM] = _win_NFW(k*(_nu**eta), _rv, _c)[:, 0]
+        if baryons:
+            profile = halo.profile.Fourier(k, M, Uk, amplitude=M/hmod.rhom, mass_tracer=True) # NOTE: No Factor of 1-f_nu for baryonic, because this effect is already included!
+    
+        else:
+            profile = halo.profile.Fourier(k, M, Uk, amplitude=M*(1.-f_nu)/hmod.rhom, mass_tracer=True) # NOTE: Factor of 1-f_nu
+    
         # Vanilla power spectrum calculation
-        _, Pk_1h, _ = hmod.power_spectrum(k, Pk_lin, M, sigmaM, {'m': profile}, simple_twohalo=True)
+        _, Pk_1h_, _ = hmod.power_spectrum(k, Pk_lin, M, sigmaM, {'m': profile}, simple_twohalo=True)
 
         # HMcode tweaks
         P_wig = _get_Pk_wiggle(k, Pk_lin, CAMB_results)   # Isolate spectral wiggle; footnote 7
         Pk_dwl = Pk_lin-(1.-np.exp(-(k*sigmaV)**2))*P_wig # Constuct linear spectrum with smoothed wiggle; equation (15)
         Pk_2h = Pk_dwl*(1.-f*(k/kd)**nd/(1.+(k/kd)**nd))  # Two-halo term; equation (16)
-        Pk_1h = (k/ks)**4/(1.+(k/ks)**4)*Pk_1h['m-m']     # One-halo term; equation (17)
+        Pk_1h = (k/ks)**4/(1.+(k/ks)**4)*Pk_1h_['m-m']     # One-halo term; equation (17)
         Pk_hm = (Pk_2h**alpha+Pk_1h**alpha)**(1./alpha)   # Total prediction via smoothed sum; equation (23)
+
+        if hmcodetweaks==False:
+            Pk_hm=Pk_lin+Pk_1h_['m-m']
+
         Pk_HMcode[iz, :] = Pk_hm
 
     # Finish
@@ -276,4 +325,59 @@ def _win_NFW(k:np.ndarray, rv:np.ndarray, c:np.ndarray) -> np.ndarray:
     f3 = np.sin(kv)/(ks+kv)
     f4 = np.log(1.+c)-c/(1.+c)
     Wk = (f1+f2-f3)/f4
-    return Wk
+    return Wk
+
+
+def _win_NFW_baryons(k:np.ndarray, rv:np.ndarray, c:np.ndarray, M:np.ndarray, Mb:float, f:float, Om_b:float, Om_m:float, Om_c:float, rhom:float) -> np.ndarray:
+    """Normalised Fourier transform for NFW profile, including baryonic effects
+    Is the same as Equation (25) of Mead et al. (2021)
+    Calls _win_NFW
+    """
+    Wk=_win_NFW(k, rv, c)
+
+    fg=(Om_b/Om_m-f)*(M/Mb)**2/(1+(M/Mb)**2) #Gas content (Eq. 24 with beta=2 specified)
+    Wk=((Om_c)/Om_m+fg)*Wk
+    Wk+=f
+ 
+    return Wk
+
+
+def get_Baryon_Parameters(T_AGN:float) -> dict:
+    """Maps 1-Param baryon feedback model from HMCode2020 to 6 baryonic parameters
+    Uses parameters from Table 5 in Mead et al. (2021)
+    Warning:
+        This fit was obtained using the vanilla halo model! If the hmcode tweaks are used, different values are likely needed.
+
+    Args:
+        T_AGN (float): AGN Temperature parameter
+
+    Returns:
+        dict: baryonification parameters, B0, Bz, Mb0, Mbz, f0, fz
+    """
+
+    theta=np.log10(T_AGN/np.power(10,7.8))
+    params={}
+    params["B0"]=3.44-0.496*theta
+    params["Bz"]=-0.0671-0.0371*theta
+    params["Mb0"]=np.power(10, 13.87+1.81*theta) #Units: Msun/h
+    params["Mbz"]=-0.108+0.195*theta
+    params["f0"]=(2.01-0.3*theta)*1e-2
+    params["fz"]=0.409+0.0224*theta
+    return params
+
+def get_Baryon_Suppression(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,  T_AGN=None, hmcodetweaks=False,
+          Mmin=1e0, Mmax=1e18, nM=256, verbose=False):
+    """Calculates the ratio of the powerspectrum with baryonic effect to the dark-matter-only powerspectrum for a given T_AGN
+    Assumes the 1-param model from HMCode2020
+    Warning:
+        Since the fit for the baryonic effects was obtained with the vanilla halo model, it is not safe to set hmcodetweaks=True
+    """
+    
+ 
+    baryon_params=get_Baryon_Parameters(T_AGN)
+ 
+    Pk_dmo=power(k, zs, CAMB_results, Mmin, Mmax, nM, verbose, baryons=False, hmcodetweaks=hmcodetweaks)
+
+    Pk_bar=power(k, zs, CAMB_results, Mmin, Mmax, nM, verbose, baryons=True, baryon_params=baryon_params, hmcodetweaks=hmcodetweaks)
+
+    return Pk_bar/Pk_dmo

notebooks/demo_baryons.ipynb

@@ -0,0 +1,26 @@
+# Om_c, Om_b, Om_k, h, ns, sig8, w0, wa, m_nu, T_AGN
+3.547912097111927121e-01 4.816635319256157288e-02 3.585979199113824289e-02 7.789472116237455834e-01 9.094177347887649754e-01 8.951244703273512071e-01 -8.433161788057881303e-01 6.360535588836309095e-01 1.281136326755458743e-01 8.227542920065280795e+07
+2.741596048465162472e-01 6.280294966545806046e-02 1.438651200806645436e-02 8.291046453083319445e-01 9.443414198827331241e-01 7.454477443569553774e-01 -9.672491277904990969e-01 -1.537910059126432305e+00 8.276311719925820709e-01 1.101958094499628246e+08
+3.516175480170747880e-01 4.563577904389605477e-02 4.706980243949032694e-02 8.572484485288791589e-01 9.778383497073761532e-01 7.389277415703934260e-01 -1.019967397763779537e+00 -1.598150664980441338e+00 1.542894920675478287e-01 1.197110574979204834e+08
+3.489524311815634383e-01 6.402529197302631037e-02 -1.741746418618480424e-02 6.481838824139475319e-01 9.469555811275808255e-01 7.378942718168570725e-01 -1.222047096798717147e+00 -2.981281720599395779e-01 2.269093490508841127e-01 1.171843988015893996e+08
+2.874303837744661694e-01 5.998034588173512677e-02 2.002651020022491735e-02 6.249466565528164486e-01 9.832259801395201171e-01 8.609528714993603948e-01 -1.067512972581895259e+00 -8.621324071699653091e-01 6.824955039749754926e-01 4.986854975460082293e+07
+2.399816404950216864e-01 3.522086809253016648e-02 2.869243775021383669e-02 7.659403426368128764e-01 9.705165378626334771e-01 8.561458062043936224e-01 -1.024650534676996072e+00 -1.808900018178993818e-02 1.397969981276574458e-01 4.788186429758129269e+07
+3.336805923580943301e-01 4.913288618429397653e-02 6.523610648118884081e-03 8.059995429664101874e-01 9.634718320000590364e-01 8.107158801315991203e-01 -9.644757035527518063e-01 -8.151102048315372883e-01 3.081783456793940612e-02 8.048263235451522470e+07
+2.429169345639058397e-01 4.725585931173909016e-02 3.534030732681661680e-02 5.935757943461362762e-01 9.058302741689066018e-01 7.562767784043993302e-01 -1.123843745339989875e+00 2.623687093279540861e-01 5.570321523412783415e-01 1.408408853234213293e+08
+3.328627080654775194e-01 4.719160584320211682e-02 3.140203846660347131e-02 5.667891679630815416e-01 9.022712073133860589e-01 7.180095721551282839e-01 -8.665843896421298043e-01 -3.397495369433241041e-01 1.612717790336017920e-01 8.927530499631467462e+07
+2.304624205426336891e-01 5.588961125233208455e-02 -5.384372442596929709e-03 6.524084904385929473e-01 9.301512089147876416e-01 8.260565186237777136e-01 -1.082912433667965857e+00 -1.466173742858535967e+00 1.180059021205153158e-01 1.876412706450518370e+08
+3.817161381415214438e-01 5.599121401432248513e-02 -2.341300385404804230e-02 8.876705509390896420e-01 9.778750903965794938e-01 8.433780378317990989e-01 -1.030383098713726797e+00 -9.105528988460713791e-01 9.639096215349929331e-02 1.705382011031487584e+08
+2.911552579667222584e-01 4.107090094385690976e-02 -1.940433758493475061e-02 7.316878275767584316e-01 9.176772782939232043e-01 8.713228568184750999e-01 -8.448882820988738995e-01 4.355834974123196091e-01 4.320930397751037155e-01 1.094249075826266259e+08
+3.168195937825471487e-01 5.449539804664460529e-02 -4.155556788601109247e-02 6.663229608682438299e-01 9.041614173861892700e-01 7.987981638489037861e-01 -1.102083272600328900e+00 -1.294982191513198710e+00 1.034029677225516419e-01 1.026481674714202583e+08
+2.341185937073772261e-01 6.275360355130391032e-02 8.106113970039498240e-03 6.387479218139348047e-01 9.590915491481416533e-01 7.045607742059394329e-01 -7.248644720551327136e-01 -2.782666548018704056e-01 7.827352272502862141e-01 4.548948289099365473e+07
+2.973316661676320694e-01 4.972120983063563204e-02 4.378264549749828760e-02 7.286912209504301519e-01 9.473489401056953696e-01 7.533951326183786490e-01 -1.101058601594468778e+00 -1.627760685606713942e-01 4.389114603050466856e-01 4.122194276090595871e+07
+3.652583848388715748e-01 6.188482315519300281e-02 -3.597509110013892564e-02 7.216144574156198033e-01 9.108575741135443993e-01 8.344480186079623607e-01 -1.131259729696594940e+00 2.548621304226241779e-01 7.269946142868826122e-01 1.374210442375697792e+08
+2.215481891911793288e-01 6.248035535412824193e-02 -2.697860091051192122e-02 5.149650224704719337e-01 9.554852469391483805e-01 7.741844567724877058e-01 -8.021261541205519929e-01 7.028369309135484855e-01 3.171388928227153459e-01 1.849394563120471835e+08
+2.581835676280237468e-01 5.045171387695143728e-02 -2.440349094323972581e-02 8.744174280195853299e-01 9.164607817582017812e-01 7.089821238784657087e-01 -1.038941763998177281e+00 9.539485714272908368e-01 7.486080194569492141e-01 1.673226536376067102e+08
+3.786893279395726619e-01 5.056575081159347257e-02 -1.840709481692070354e-02 8.088049728443952313e-01 9.661661263167761193e-01 7.747315457747420142e-01 -1.243319999163090772e+00 5.178367301605679707e-01 2.624605159228646789e-01 1.802059828031317294e+08
+2.481941150011369723e-01 3.868273797234458161e-02 3.311126721249062210e-02 5.613137266497976174e-01 9.179268308157739753e-01 8.198765583041687233e-01 -7.752627754975213170e-01 -1.138731656199134523e+00 3.103236729000947713e-01 1.393745206630210578e+08
+3.943652852121934949e-01 5.002223558607027387e-02 -3.561024974487492223e-02 5.055745150832806623e-01 9.229656029998855038e-01 7.263644435573041180e-01 -8.934047958322854921e-01 -1.363284161068121225e+00 5.063299316206330003e-01 1.218943862392111719e+08
+3.162233218441805049e-01 4.099326954980173077e-02 3.041245261822626833e-02 7.861628518463206206e-01 9.738984003915541754e-01 7.262115503114625614e-01 -1.225747717809793302e+00 1.061963278529587829e+00 3.975781938249406400e-01 6.466422608995757252e+07
+2.977168090703066761e-01 5.488592638290747117e-02 4.556232570469699594e-02 6.145784907528220042e-01 9.924808429312027558e-01 7.049718982772512632e-01 -9.668811746039052135e-01 1.782650861600660885e-01 1.058974037507532939e-01 4.991576345134069026e+07
+2.838228638632607770e-01 6.398695736429545200e-02 9.604255323437282943e-03 8.732092886400845133e-01 9.804360915612970739e-01 7.934763203110583030e-01 -8.291419304486875586e-01 -1.676311280229267009e+00 1.091439967657349408e-01 1.515653861853125691e+08
+3.593634176650323808e-01 4.197922225899301113e-02 3.076959059905348681e-03 7.424063282800044128e-01 9.867738953776010735e-01 8.206214314677451327e-01 -1.052457058435791737e+00 -6.037060293443801129e-01 4.258820863735098827e-01 1.138548978291175067e+08