Commit 34301d1b authored by Xavier Garrido's avatar Xavier Garrido
Browse files

add bias calculation

parent 0338a2e0
......@@ -57,6 +57,7 @@ namespace cib {
_critical_densities_.resize(nz);
cib::zeros(_hmfmz_, _mh_.size(), nz);
cib::zeros(_biasmz_, _mh_.size(), nz);
const size_t xstep = 50;
const double xmin = -6, xmax = 1;
......@@ -66,7 +67,6 @@ namespace cib {
cib::read_from_file("data/snu.txt", _snu_eff_);
cib::read_from_file("data/u_int.txt", _unfw_);
cib::read_from_file("data/Pk_int.txt", _power_);
cib::read_from_file("data/bmz_int.txt", _biasmz_);
cib::read_from_file("data/m500c.txt", _m500c_);
cib::read_from_file("data/k.txt", _k_);
......@@ -232,10 +232,8 @@ namespace cib {
void cib_interface::_update_hmf_bias_()
{
// cib::array1d & zvalues = cosmo->redshift;
auto fsigma =
[this](double z_, double omz_, double sigma_) -> double
[this](double z_, double sigma_, double delta_halo_) -> double
{
const cib::array1d delta_virs = {200, 300, 400, 600, 800, 1200, 1600, 2400, 3200};
const cib::array1d A_array = {0.185866, 0.199597, 0.211566, 0.218411, 0.248097, 0.254605, 0.260000, 0.260000, 0.260000};
......@@ -245,8 +243,7 @@ namespace cib {
const double A_exp = -0.14;
const double a_exp = 0.06;
const double delta_halo = 200;
const double dhalo = delta_halo/omz_;
const double dhalo = delta_halo_;
//interp linear
size_t idx = 0;
......@@ -267,19 +264,33 @@ namespace cib {
};
const double mean_density0 = _klass_->Omega_m() * _klass_->critical_density(0.);
const double delta_halo = 200;
for (size_t iz = 0; iz < _z_.size(); iz++) {
const double z = _z_[iz];
const double rho_z = _critical_densities_[iz];
const double omega_m_z = _omega_m_z_[iz];
// cib::Filter filt(cosmo, k, Pk[iz], zvalues[iz], rho_z, Omega_m_z, delta_halo);
for (size_t im = 0; im < _mh_.size(); im++) {
const double mass = _mh_[im];
const double radius = std::pow(3 * mass / (4*CLHEP::pi*mean_density0), 1./3.);
double sigma1 = 0.0, sigma2 = 0.0;
_compute_sigmas_(radius, _Pk_[iz], sigma1, sigma2);
const double fsig = fsigma(z, omega_m_z, sigma1);
// Halo mass function
const double fsig = fsigma(z, sigma1, delta_halo/omega_m_z);
const double dn_dm = fsig * mean_density0 * std::fabs(0.5*sigma2*1./3.)/std::pow(mass, 2);
_hmfmz_[im][iz] = mass * dn_dm * std::log(10);
// Bias
const double y = std::log10(delta_halo/omega_m_z);
const double yy = std::exp(-std::pow(4./y, 4));
const double A = 1.0 + 0.24*y*yy;
const double aa = 0.44*y - 0.88;
const double C = 0.019 + 0.107*y + 0.19*yy;
const double dc = 1.686; // critical density of the universe. Redshift evolution is small
// and neglected
const double nuu = dc/sigma1;// length of mass array
const double B = 0.183, b = 1.5, c = 2.4;
_biasmz_[im][iz] = 1. - (A*std::pow(nuu, aa)/(std::pow(nuu, aa) + std::pow(dc, aa)))
+ B*std::pow(nuu, b) + C*std::pow(nuu, c);
}
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment