update halo mass function

src/camel/Astro/cib_interface.cc

@@ -47,23 +47,31 @@ namespace cib {
 
     // Redshift & comoving distances
     cib::read_from_file("data/z.txt", _z_);
-    _comoving_distances_.resize(_z_.size());
-    _Ez_.resize(_z_.size());
-    _dVcdz_.resize(_z_.size());
-    _dchidz_.resize(_z_.size());
+    const size_t nz = _z_.size();
+    _comoving_distances_.resize(nz);
+    _Ez_.resize(nz);
+    _dVcdz_.resize(nz);
+    _dchidz_.resize(nz);
+    _omega_m_z_.resize(nz);
+    _omega_b_z_.resize(nz);
+    _critical_densities_.resize(nz);
+
+    cib::zeros(_hmfmz_, _mh_.size(), nz);
 
     const size_t xstep = 50;
     const double xmin = -6, xmax = 1;
     const double dx = (xmax - xmin)/(xstep-1);
     for (size_t i = 0; i < xstep; i++) _x_.push_back(std::pow(10, xmin + i*dx));
 
-    cib::read_from_file("data/hmf_red.txt", _hmfmz_);
     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_);
+    cib::read_from_file("data/pk.txt", _Pk_);
+
     _initialized_ = true;
   }
 
@@ -83,20 +91,8 @@ namespace cib {
       _gnu_[i] = x*((std::exp(x) + 1)/(std::exp(x) - 1)) - 4;
     }
 
-    // Update comoving distances for each z given cosmo. parameters
-    const double Om0 = _klass_->Omega_m();
-    const double Ode0 = _klass_->Omega_L();
-    const double H0 = _klass_->get_H0();  //Km.s-1.Mpc-1
-
-    for (size_t i = 0; i < _z_.size(); i++) {
-      const double z = _z_[i];
-      _comoving_distances_[i] = _klass_->com_distance(z);
-
-      const double ez = std::sqrt(Om0*std::pow(1+z, 3) + Ode0);
-      _Ez_[i] = ez;
-      _dVcdz_[i] = _c_*std::pow(_comoving_distances_[i], 2) / (H0 * ez);
-      _dchidz_[i] = _c_*1e-3/ (H0 * ez);
-    }
+    _update_cosmology_();
+    _update_hmf_bias_();
   }
 
   void cib_interface::compute_cls(cib::array3d & cl_cib_,
@@ -213,6 +209,103 @@ namespace cib {
 
 
   // Private methods
+  void cib_interface::_update_cosmology_()
+  {
+    // Update comoving distances for each z given cosmo. parameters
+    const double Om0 = _klass_->Omega_m();
+    const double Ode0 = _klass_->Omega_L();
+    const double H0 = _klass_->get_H0();  //Km.s-1.Mpc-1
+
+    for (size_t i = 0; i < _z_.size(); i++) {
+      const double z = _z_[i];
+      _comoving_distances_[i] = _klass_->com_distance(z);
+      _omega_b_z_[i] = _klass_->get_Ob(z);
+      _omega_m_z_[i] = _klass_->get_Om(z);
+      _critical_densities_[i] = _klass_->critical_density(z);
+
+      const double ez = std::sqrt(Om0*std::pow(1+z, 3) + Ode0);
+      _Ez_[i] = ez;
+      _dVcdz_[i] = _c_*std::pow(_comoving_distances_[i], 2) / (H0 * ez);
+      _dchidz_[i] = _c_*1e-3/ (H0 * ez);
+    }
+  }
+
+  void cib_interface::_update_hmf_bias_()
+  {
+    // cib::array1d & zvalues = cosmo->redshift;
+
+    auto fsigma =
+      [this](double z_, double omz_, double sigma_) -> 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};
+        const cib::array1d a_array = {1.466904, 1.521782, 1.559186, 1.614585, 1.869936, 2.128056, 2.301275, 2.529241, 2.661983};
+        const cib::array1d b_array = {2.571104, 2.254217, 2.048674, 1.869559, 1.588649, 1.507134, 1.464374, 1.436827, 1.405210};
+        const cib::array1d c_array = {1.193958, 1.270316, 1.335191, 1.446266, 1.581345, 1.795050, 1.965613, 2.237466, 2.439729};
+
+        const double A_exp = -0.14;
+        const double a_exp = 0.06;
+        const double delta_halo = 200;
+        const double dhalo = delta_halo/omz_;
+
+        //interp linear
+        size_t idx = 0;
+        while (delta_virs[idx+1] < dhalo) idx++;
+        if (idx >= delta_virs.size()) idx = delta_virs.size()-1;
+        const double dx = double(dhalo - delta_virs[idx])/double(delta_virs[idx+1]-delta_virs[idx]);
+        const double A0 = A_array[idx] + (A_array[idx+1]-A_array[idx])*dx;
+        const double a0 = a_array[idx] + (a_array[idx+1]-a_array[idx])*dx;
+        const double b0 = b_array[idx] + (b_array[idx+1]-b_array[idx])*dx;
+        const double c0 = c_array[idx] + (c_array[idx+1]-c_array[idx])*dx;
+
+        const double s = sigma_;
+        const double A = A0 * std::pow(1.+z_, A_exp);
+        const double a = a0 * std::pow(1.+z_, a_exp);
+        const double alpha = std::pow(10, -std::pow(0.75/std::log10(dhalo/75.), 1.2));
+        const double b = b0 * std::pow(1.+z_, alpha);
+        return A * (std::pow(s/b, -a) + 1)*std::exp(-c0/std::pow(s, 2));
+      };
+
+    const double mean_density0 = _klass_->Omega_m() * _klass_->critical_density(0.);
+    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);
+        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);
+      }
+    }
+  }
+
+  void cib_interface::_compute_sigmas_(const double rad_, const cib::array1d & Pk_,
+                                       double & sigma1_, double & sigma2_) const
+  {
+    cib::array1d integ1;
+    cib::array1d integ2;
+    cib::zeros(integ1, _k_.size());
+    cib::zeros(integ2, _k_.size());
+    for (size_t i = 0; i < _k_.size(); i++) {
+      const double rk = rad_*_k_[i];
+      double wrk1 = 1.0;
+      if (rk > 1.4e-6) wrk1 = 3*(std::sin(rk)-rk*std::cos(rk)) / std::pow(rk, 3);
+      integ1[i] = Pk_[i] * std::pow(_k_[i],3) * std::pow(wrk1, 2);
+      double wrk2 = 0.0;
+      if (rk > 1e-3) wrk2 = (9*rk*std::cos(rk)+3*std::sin(rk)*(rk*rk-3)) / std::pow(rk, 3);
+      integ2[i] = Pk_[i] * std::pow(_k_[i],3) * wrk1 * wrk2;
+    }
+
+    // dk = 1e-7
+    const double dlnk = std::log(_k_[1] / _k_[0]);
+    sigma1_ = std::sqrt(0.5 / std::pow(CLHEP::pi, 2) * cib::integration_simps(integ1, dlnk));
+    sigma2_ = 1.0 / std::pow(CLHEP::pi*sigma1_, 2) * cib::integration_simps(integ2, dlnk);
+  }
 
   void cib_interface::_compute_djsub_dlnMh_(cib::array3d & dj_sub_) const
   {
@@ -281,8 +374,8 @@ namespace cib {
     cib::zeros(sfr_, mhalo_.size(), _z_.size());
     for (size_t i = 0; i < mhalo_.size(); i++) {
       for (size_t j = 0; j < _z_.size(); j++) {
-        const double ob = _klass_->get_Ob(_z_[i]);
-        const double om = _klass_->get_Om(_z_[i]);
+        const double ob = _omega_b_z_[j];
+        const double om = _omega_m_z_[j];
         sfr_[i][j] = mhdot[i][j] * ob/om * sfr_mhdot[i][j];
       }
     }
@@ -416,7 +509,7 @@ namespace cib {
     const double delta_h_tsz = 500;
     for (size_t m = 0; m < _mh_.size(); m++) {
       for (size_t z = 0; z < _z_.size(); z++) {
-        const double rho_crit = _klass_->critical_density(_z_[z]);
+        const double rho_crit = _critical_densities_[z];
         const double r3 = 3*_m500c_[m][z]/(4*CLHEP::pi*delta_h_tsz*rho_crit);
         r_delta_[m][z] = std::pow(r3, 1./3.);
       }

src/camel/Astro/cib_interface.h

@@ -41,6 +41,12 @@ namespace cib {
 
   private:
 
+    /// Update cosmplogy quantities
+    void _update_cosmology_();
+
+    /// Update halo mass function and bias
+    void _update_hmf_bias_();
+
     /// Compute djc_dlnMh
     void _compute_djc_dlnMh_(cib::array3d & dj_cen_) const;
 
@@ -64,6 +70,8 @@ namespace cib {
 
     void _compute_P_e_(cib::array3d & pe_) const;
 
+    void _compute_sigmas_(const double rad_, const cib::array1d & Pk_, double & sigma1_, double & sigma2_) const;
+
   private:
 
     /// Initalized flag
@@ -101,6 +109,9 @@ namespace cib {
 
     /// Comoving distances in Mpc
     cib::array1d _comoving_distances_;
+    cib::array1d _omega_m_z_;
+    cib::array1d _omega_b_z_;
+    cib::array1d _critical_densities_;
 
     ///
     cib::array1d _dVcdz_;
@@ -114,6 +125,10 @@ namespace cib {
     cib::array2d _biasmz_;
     cib::array2d _m500c_;
     cib::array1d _x_;
+
+    // Internal k, pk
+    cib::array1d _k_;
+    cib::array2d _Pk_;
   };
 
 } // end of namespace