<jupyter_start><jupyter_text>Sanity check to see if there is any difference between:- (1) using an analytic calculation of a baryonification correction factor to then apply to TabCorr measurement- (2) applying baryonification displacement to particles and then using TabCorr to calculate DS<jupyter_code>import numpy as np

from matplotlib import cm
from matplotlib import colors
import matplotlib.pyplot as plt

import pandas as pd

import time

from scipy.interpolate import splrep, splev
from scipy.integrate import quad

from halotools.sim_manager import UserSuppliedHaloCatalog, UserSuppliedPtclCatalog,FakeSim, CachedHaloCatalog
from halotools.mock_observables import delta_sigma, wp, return_xyz_formatted_array
from halotools.empirical_models import NFWProfile
from halotools.empirical_models import PrebuiltHodModelFactory, HodModelFactory
from halotools.empirical_models import AssembiasZheng07Cens, AssembiasZheng07Sats, TrivialPhaseSpace, NFWPhaseSpace
from halotools.utils import add_halo_hostid
from halotools.mock_observables.surface_density.surface_density_helpers import annular_area_weighted_midpoints

import baryonification as bfc
from baryonification.useful_functions import DeltaSigmas_from_density_profile, cvir_fct

from tabcorr import TabCorr

from astropy.cosmology import Planck15
redshift = (1./0.65650)-1. #z=0.523
Lbox = 1000. #Mpc/h
# Now, we're ready to calculate correlation functions for a specific model.
HOD_model = PrebuiltHodModelFactory('zheng07', redshift = redshift, cosmology=Planck15, populate=False,
                                         secondary_haloprop_key='halo_nfw_conc')


def hod_model_from_params(params, HOD_model):
    
    HOD_model.param_dict['logMmin'] = params[0]
    HOD_model.param_dict['sigma_logM'] = params[1]
    HOD_model.param_dict['logM0'] = params[2]
    HOD_model.param_dict['logM1'] = params[3]
    HOD_model.param_dict['alpha'] = params[4]
    
    return HOD_model
HOD_params = [13.031,0.38,13.27,14.08,0.76] #params from https://ui.adsabs.harvard.edu/abs/2014MNRAS.444..476R/abstract

HOD_model = hod_model_from_params(HOD_params, HOD_model)<jupyter_output><empty_output><jupyter_text>(1) using an analytic calculation of a baryonification correction factor to then apply to TabCorr measurement<jupyter_code># We could read it in like this. Thus, we can skip the previous steps in the
# future.
halotab = TabCorr.read('mdpl2_tabCorr_DS.hdf5')
baryonic_params = [2.3e13, 0.31, 4.0] #[Mc, mu, thej] Model A
# baryonic_params = [6.6e13, 0.21, 4.0] #[Mc, mu, thej] Model B

ngal, ds_baryons_in_mpc = halotab.predict_with_baryons(HOD_model,baryonic_params)
ngal, ds_in_mpc = halotab.predict(HOD_model)


length_unit_conversion = 1e6
ds_baryons_in_pc = ds_baryons_in_mpc/length_unit_conversion**2
ds_in_pc = ds_in_mpc/length_unit_conversion**2<jupyter_output>UPDATED
Normalizing power-spectrum done!<jupyter_text>(2) applying baryonification displacement to particles and then using TabCorr to calculate DS<jupyter_code># We could read it in like this. Thus, we can skip the previous steps in the
# future.
halotab_displaced = TabCorr.read('mdpl2+baryonification_tabCorr_DS.hdf5')
ngal, ds_displaced_in_mpc = halotab_displaced.predict(HOD_model)


ds_displaced_in_pc = ds_displaced_in_mpc/length_unit_conversion**2
rbin = annular_area_weighted_midpoints(np.logspace(-1, 1, 20))
plt.figure(figsize=[10,8])
plt.plot(rbin,ds_in_pc, label='MDPL2')
plt.plot(rbin, ds_baryons_in_pc, label='MDPL2+analytic correction factor')
plt.plot(rbin, ds_displaced_in_pc, label='MDPL2 displaced particles')
plt.loglog()
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel(r'$r_{\rm p} $  $\rm{[Mpc/h]}$', fontsize=25)
plt.ylabel(r'$\Delta\Sigma(R_{\rm p})$  $[h M_{\odot} / {\rm pc}^2]$', fontsize=25)
plt.legend(loc='best', fontsize=20)
plt.show()
plt.figure(figsize=[10,8])
plt.plot(rbin, ds_baryons_in_pc/ds_in_pc, label='MDPL2+analytic correction factor')
plt.plot(rbin, ds_displaced_in_pc/ds_in_pc, label='MDPL2 displaced particles')
plt.axhline(1, linestyle='--', color='k')
plt.semilogx()
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel(r'$r_{\rm p} $  $\rm{[Mpc/h]}$', fontsize=25)
plt.ylabel(r'$\Delta\Sigma_{baryons}/\Delta\Sigma_{DM}$', fontsize=25)
plt.legend(loc='best', fontsize=20)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10680287903413048, 'fraction_numerical': 0.03018342233573253, 'mean_word_length': 6.584507042253521, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '447538', 'license': 'MIT', 'n_tokens_mistral': 1691, 'n_tokens_neox': 1561, 'n_words': 332, 'path': 'TabCorr/Sanity_Check_compare_DS_with_TabCorr_for_analytic_vs_displacement.ipynb', 'repo_name': 'f-ardila/baryonic_effects', 'size': '4307'}