Skip to content

Commit c994111

Browse files
author
llinke1
committed
Added model for baryonic effects from HMCode2020
1 parent d3a0612 commit c994111

30 files changed

+680
-10
lines changed

‎hmcode/hmcode.py‎

Lines changed: 114 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,8 @@
1616
sigma_cold_approx = False # Should the Eisenstein & Hu (1999) approximation be used for the cold transfer function?
1717

1818
def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
19-
Mmin=1e0, Mmax=1e18, nM=256, verbose=False) -> np.ndarray:
19+
Mmin=1e0, Mmax=1e18, nM=256, verbose=False,
20+
baryons=False, baryon_params=None, hmcodetweaks=True) -> np.ndarray:
2021
'''
2122
Calculates the HMcode matter-matter power spectrum
2223
Args:
@@ -26,6 +27,9 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
2627
Mmin: Minimum mass for the halo-model calculation [Msun/h]
2728
Mmax: Maximum mass for the halo-model calculation [Msun/h]
2829
nM: Number of mass bins for the halo-model calculation
30+
baryons: if true: include baryonic effects from HMCode2020; if false: Dark Matter Only
31+
baryon_params: Dictionary containing the 6 baryonic parameters from HMCode2020
32+
hmcodetweaks: if true: Use the changes to the vanilla halo model from HMCode2020; if false: Use vanilla halo model
2933
Returns:
3034
Array of matter power spectra: Pk[z, k]
3135
'''
@@ -60,6 +64,16 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
6064
print('Linear growth at z=0: {:.3}'.format(growth(1.)))
6165
print()
6266

67+
if verbose and baryons:
68+
print("Using baryonic feedback model from HMCode2020")
69+
print("B0: {:.4f}".format(baryon_params["B0"]))
70+
print("Bz: {:.4f}".format(baryon_params["Bz"]))
71+
print("Mb0: {:.1e}".format(baryon_params["Mb0"]))
72+
print("Mbz: {:.4f}".format(baryon_params["Mbz"]))
73+
print("f*0: {:.4f}".format(baryon_params["f0"]))
74+
print("f*z: {:.4f}".format(baryon_params["fz"]))
75+
print()
76+
6377
# Linear power interpolator
6478
interp, _, k_interp = CAMB_results.get_matter_power_interpolator(nonlinear=False,
6579
return_z_k=True,
@@ -96,6 +110,7 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
96110
print('Linear collapse threshold: {:.4}'.format(dc))
97111
print('Virial halo overdensity: {:.4}'.format(Dv))
98112

113+
99114
# Initialise halo model
100115
hmod = halo.model(z, Om_m, name='Sheth & Tormen (1999)', Dv=Dv, dc=dc)
101116

@@ -135,9 +150,16 @@ def power(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata,
135150
nd = 2.853 # Two-halo damping power; equation (16)
136151
ks = 0.05618*sigma8**-1.013 # One-halo damping wavenumber; equation (17)
137152
eta = 0.1281*sigma8**-0.3644 # Halo bloating parameter; equation (19)
138-
B = 5.196 # Minimum halo concentration; equation (20)
153+
B = 5.196 # Minimum halo concentration; equation (20)
139154
alpha = 1.875*(1.603)**neff # Transition smoothing; equation (23)
140-
if verbose:
155+
156+
157+
if hmcodetweaks==False: #Go to vanilla halo model if no hmcode tweaks should be used
158+
B=4
159+
eta=0
160+
alpha=1
161+
162+
if verbose and hmcodetweaks:
141163
print('Two-halo damping wavenumber; kd: {:.4} h/Mpc'.format(kd))
142164
print('Two-halo fractional damping; f: {:.4}'.format(f))
143165
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,
146168
print('Minimum halo concentration; B: {:.4}'.format(B))#, B/4.)
147169
print('Transition smoothing; alpha: {:.4}'.format(alpha))
148170
print()
149-
171+
elif verbose:
172+
print("Using vanilla halo model")
173+
174+
if baryons:
175+
B=baryon_params["B0"]*np.power(10,z*baryon_params["Bz"])
176+
Mb=baryon_params["Mb0"]*np.power(10, z*baryon_params["Mbz"])
177+
fstar=baryon_params["f0"]*np.power(10, z*baryon_params["fz"])
178+
179+
180+
181+
if verbose and baryons:
182+
print("Halo concentration parameter B:{:.4f}".format(B))
183+
print("Gas loss halo mass parameter Mb:{:1e}".format(Mb))
184+
print("Effective halo stellar mass fraction fstar:{:.4f}".format(f))
185+
186+
150187
# Halo concentration
151188
zf = _get_halo_collapse_redshifts(M, z, iz, dc, growth, CAMB_results, cold=True) # Halo formation redshift
152189
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,
157194
rv = hmod.virial_radius(M)
158195
Uk = np.ones((len(k), len(M)))
159196
for iM, (_rv, _c, _nu) in enumerate(zip(rv, c, nu)): # TODO: Remove loop for speed?
160-
Uk[:, iM] = _win_NFW(k*(_nu**eta), _rv, _c)[:, 0]
161-
profile = halo.profile.Fourier(k, M, Uk, amplitude=M*(1.-f_nu)/hmod.rhom, mass_tracer=True) # NOTE: Factor of 1-f_nu
162-
197+
if baryons:
198+
Uk[:, iM] = _win_NFW_baryons(k*(_nu**eta), _rv, _c, M[iM], Mb, fstar, Om_b, Om_m, Om_c, hmod.rhom)[:,0]
199+
200+
else:
201+
Uk[:, iM] = _win_NFW(k*(_nu**eta), _rv, _c)[:, 0]
202+
if baryons:
203+
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!
204+
205+
else:
206+
profile = halo.profile.Fourier(k, M, Uk, amplitude=M*(1.-f_nu)/hmod.rhom, mass_tracer=True) # NOTE: Factor of 1-f_nu
207+
163208
# Vanilla power spectrum calculation
164-
_, Pk_1h, _ = hmod.power_spectrum(k, Pk_lin, M, sigmaM, {'m': profile}, simple_twohalo=True)
209+
_, Pk_1h_, _ = hmod.power_spectrum(k, Pk_lin, M, sigmaM, {'m': profile}, simple_twohalo=True)
165210

166211
# HMcode tweaks
167212
P_wig = _get_Pk_wiggle(k, Pk_lin, CAMB_results) # Isolate spectral wiggle; footnote 7
168213
Pk_dwl = Pk_lin-(1.-np.exp(-(k*sigmaV)**2))*P_wig # Constuct linear spectrum with smoothed wiggle; equation (15)
169214
Pk_2h = Pk_dwl*(1.-f*(k/kd)**nd/(1.+(k/kd)**nd)) # Two-halo term; equation (16)
170-
Pk_1h = (k/ks)**4/(1.+(k/ks)**4)*Pk_1h['m-m'] # One-halo term; equation (17)
215+
Pk_1h = (k/ks)**4/(1.+(k/ks)**4)*Pk_1h_['m-m'] # One-halo term; equation (17)
171216
Pk_hm = (Pk_2h**alpha+Pk_1h**alpha)**(1./alpha) # Total prediction via smoothed sum; equation (23)
217+
218+
if hmcodetweaks==False:
219+
Pk_hm=Pk_lin+Pk_1h_['m-m']
220+
172221
Pk_HMcode[iz, :] = Pk_hm
173222

174223
# Finish
@@ -276,4 +325,59 @@ def _win_NFW(k:np.ndarray, rv:np.ndarray, c:np.ndarray) -> np.ndarray:
276325
f3 = np.sin(kv)/(ks+kv)
277326
f4 = np.log(1.+c)-c/(1.+c)
278327
Wk = (f1+f2-f3)/f4
279-
return Wk
328+
return Wk
329+
330+
331+
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:
332+
"""Normalised Fourier transform for NFW profile, including baryonic effects
333+
Is the same as Equation (25) of Mead et al. (2021)
334+
Calls _win_NFW
335+
"""
336+
Wk=_win_NFW(k, rv, c)
337+
338+
fg=(Om_b/Om_m-f)*(M/Mb)**2/(1+(M/Mb)**2) #Gas content (Eq. 24 with beta=2 specified)
339+
Wk=((Om_c)/Om_m+fg)*Wk
340+
Wk+=f
341+
342+
return Wk
343+
344+
345+
def get_Baryon_Parameters(T_AGN:float) -> dict:
346+
"""Maps 1-Param baryon feedback model from HMCode2020 to 6 baryonic parameters
347+
Uses parameters from Table 5 in Mead et al. (2021)
348+
Warning:
349+
This fit was obtained using the vanilla halo model! If the hmcode tweaks are used, different values are likely needed.
350+
351+
Args:
352+
T_AGN (float): AGN Temperature parameter
353+
354+
Returns:
355+
dict: baryonification parameters, B0, Bz, Mb0, Mbz, f0, fz
356+
"""
357+
358+
theta=np.log10(T_AGN/np.power(10,7.8))
359+
params={}
360+
params["B0"]=3.44-0.496*theta
361+
params["Bz"]=-0.0671-0.0371*theta
362+
params["Mb0"]=np.power(10, 13.87+1.81*theta) #Units: Msun/h
363+
params["Mbz"]=-0.108+0.195*theta
364+
params["f0"]=(2.01-0.3*theta)*1e-2
365+
params["fz"]=0.409+0.0224*theta
366+
return params
367+
368+
def get_Baryon_Suppression(k:np.array, zs:np.array, CAMB_results:camb.CAMBdata, T_AGN=None, hmcodetweaks=False,
369+
Mmin=1e0, Mmax=1e18, nM=256, verbose=False):
370+
"""Calculates the ratio of the powerspectrum with baryonic effect to the dark-matter-only powerspectrum for a given T_AGN
371+
Assumes the 1-param model from HMCode2020
372+
Warning:
373+
Since the fit for the baryonic effects was obtained with the vanilla halo model, it is not safe to set hmcodetweaks=True
374+
"""
375+
376+
377+
baryon_params=get_Baryon_Parameters(T_AGN)
378+
379+
Pk_dmo=power(k, zs, CAMB_results, Mmin, Mmax, nM, verbose, baryons=False, hmcodetweaks=hmcodetweaks)
380+
381+
Pk_bar=power(k, zs, CAMB_results, Mmin, Mmax, nM, verbose, baryons=True, baryon_params=baryon_params, hmcodetweaks=hmcodetweaks)
382+
383+
return Pk_bar/Pk_dmo

‎notebooks/demo_baryons.ipynb‎

Lines changed: 235 additions & 0 deletions
Large diffs are not rendered by default.

‎tests/benchmarks_baryons/cosmologies.txt‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
# Om_c, Om_b, Om_k, h, ns, sig8, w0, wa, m_nu, T_AGN
2+
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
3+
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
4+
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
5+
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
6+
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
7+
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
8+
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
9+
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
10+
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
11+
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
12+
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
13+
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
14+
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
15+
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
16+
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
17+
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
18+
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
19+
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
20+
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
21+
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
22+
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
23+
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
24+
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
25+
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
26+
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

0 commit comments

Comments
 (0)