CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
MassFunction.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2010 by Federico Marulli *
3  * federico.marulli3@unibo.it *
4  * *
5  * This program is free software; you can redistribute it and/or *
6  * modify it under the terms of the GNU General Public License as *
7  * published by the Free Software Foundation; either version 2 of *
8  * the License, or (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public *
16  * License along with this program; if not, write to the Free *
17  * Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ********************************************************************/
20 
35 #include "Cosmology.h"
36 
37 using namespace std;
38 
39 using namespace cbl;
40 
41 
42 // =====================================================================================
43 
44 
45 double cbl::cosmology::Cosmology::mass_function (const double Mass, const double redshift, const std::string model_MF, const std::string method_SS, const bool store_output, const std::string output_root, const double Delta, const std::string interpType, const int norm, const double k_min, const double k_max, const double prec, const std::string input_file, const bool is_parameter_file, const bool default_delta, const double delta_t)
46 {
47  double fact = (m_unit) ? 1 : m_hh;
48  double MASS = Mass*fact;
49 
50  const bool unit1 = true;
51 
52  double SSS = sigma2M(MASS, method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file, unit1);
53  double Sigma = sqrt(SSS);
54  double Dln_Sigma = dnsigma2M(1, MASS, method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file, unit1)*(MASS/(2.*SSS));
55 
56  double MF = m_MF_generator(MASS, Sigma, Dln_Sigma, redshift, model_MF, Delta, default_delta, delta_t)*pow(fact, 4.);
57 
58  if (m_fNL!=0 && default_delta) MF *= MF_correction(MASS, redshift, model_MF, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
59  else if (m_fNL!=0 && !default_delta) ErrorCBL("Non Gaussianity still not available for user-defined density contrast threshold!", "mass_function", "MassFunction.cpp");
60 
61  return MF;
62 }
63 
64 // =====================================================================================
65 
66 
67 double cbl::cosmology::Cosmology::mass_function_fR (const double Mass, const double redshift, const std::string model_MF, const double f_R0, const bool store_output, const double Delta, const std::string interpType, const int norm, const double k_min, const double k_max, const double prec, const std::string input_file, const bool is_parameter_file, const bool default_delta, const double delta_t)
68 {
69  double fact = (m_unit) ? 1 : m_hh;
70  double MASS = Mass*fact;
71 
72  const bool unit1 = true;
73 
74  set_RhoZero(rho_m(0., true));
75 
76  cbl::Path path;
77  string pathToMGCAMB = path.DirCosmo()+"/External/MGCAMB/";
78 
79  if (system(("cp "+pathToMGCAMB+"params_MG_backup.ini "+pathToMGCAMB+"params_MG.ini").c_str())) {}
80  int MG_flag = (f_R0==0.) ? 0 : 3;
81  string cosmoTag = (f_R0==0.) ? "LCDM" : "fR"+to_string(f_R0);
82 
83  string sed;
84  sed = "sed -i 's/MG_flag = 0/MG_flag = "+cbl::conv(MG_flag, cbl::par::fINT)+"/g' "+pathToMGCAMB+"params_MG.ini";
85  if (system(sed.c_str())) {}
86  sed = "sed -i 's/F_R0 = 0.001d0/F_R0 = "+cbl::conv(f_R0, cbl::par::fDP8)+"d0/g' "+pathToMGCAMB+"params_MG.ini";
87  if (system(sed.c_str())) {}
88 
89  double SSS = sigma2M(MASS, "MGCAMB", redshift, store_output, cosmoTag, interpType, k_max, input_file, is_parameter_file, unit1);
90  double Sigma = sqrt(SSS);
91  double Dln_Sigma = dnsigma2M(1, MASS, "MGCAMB", redshift, store_output, cosmoTag, interpType, k_max, input_file, is_parameter_file, unit1)*(MASS/(2.*SSS));
92 
93  double MF = m_MF_generator(MASS, Sigma, Dln_Sigma, redshift, 1., model_MF, Delta, default_delta, delta_t)*pow(fact, 4.);
94 
95  if (m_fNL!=0 && default_delta) MF *= MF_correction(MASS, redshift, model_MF, store_output, cosmoTag, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
96  else if (m_fNL!=0 && !default_delta) ErrorCBL("Non Gaussianity still not available for user-defined density contrast threshold!", "mass_function", "MassFunction.cpp");
97 
98  return MF;
99 }
100 
101 
102 // =====================================================================================
103 
104 
105 double cbl::cosmology::Cosmology::m_mass_function (const double Mass, std::shared_ptr<void> mass_function_params)
106 {
107  std::shared_ptr<glob::STR_MF> pp = static_pointer_cast<glob::STR_MF>(mass_function_params);
108 
109  double fact = (m_unit) ? 1 : m_hh;
110  double MASS = Mass*fact;
111 
112  const bool unit1 = true;
113 
114  double SSS = sigma2M(MASS, pp->method_SS, 0., pp->store_output, pp->output_root, pp->interpType, pp->k_max, pp->input_file, pp->is_parameter_file, unit1);
115  double Sigma = sqrt(SSS);
116  double Dln_Sigma = dnsigma2M(1, MASS, pp->method_SS, 0., pp->store_output, pp->output_root, pp->interpType, pp->k_max, pp->input_file, pp->is_parameter_file, unit1)*(MASS/(2.*SSS));
117 
118  double MF = m_MF_generator(MASS, Sigma, Dln_Sigma, pp->redshift, pp->model_MF, pp->Delta, pp-> default_delta, pp->delta_t)*pow(fact, 4.);
119 
120  if (m_fNL!=0 && pp->default_delta) MF *= MF_correction(MASS, pp->redshift, pp->model_MF, pp->store_output, pp->output_root, pp->interpType, pp->norm, pp->k_min, pp->k_max, pp->prec, pp->input_file, pp->is_parameter_file);
121  else if (m_fNL!=0 && !pp->default_delta) ErrorCBL("Non Gaussianity still not available for user-defined density contrast threshold!", "m_mass_function", "MassFunction.cpp");
122 
123  return MF;
124 }
125 
126 
127 // =====================================================================================
128 
129 
130 double cbl::cosmology::Cosmology::mass_function_fast (const double Mass, const double redshift, const std::string model_MF, const std::string method_SS, const bool store_output, const std::string output_root, const double Delta, const std::string interpType, const int norm, const double k_min, const double k_max, const double prec, const std::string input_file, const bool is_parameter_file)
131 {
132  double fact = (m_unit) ? 1 : m_hh;
133  double MASS = Mass*fact;
134 
135 
136  // ---------- read the grid file ----------
137 
138  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root, interpType, k_max);
139 
140  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
141 
142  double MMass, Sigma, Dln_Sigma;
143  vector<double> mass, sigma, dln_sigma;
144  while (fin >>MMass>>Sigma>>Dln_Sigma) {
145  mass.push_back(MMass);
146  sigma.push_back(Sigma);
147  dln_sigma.push_back(Dln_Sigma);
148  }
149 
150 
151  // ---------- compute the MF ----------
152 
153  double sig = interpolated(MASS, mass, sigma, "Steffen");
154  double dlsig = interpolated(MASS, mass, dln_sigma, "Steffen");
155 
156  double MF = m_MF_generator(MASS, sig, dlsig, redshift, model_MF, Delta)*pow(fact, 4.);
157  if (std::isnan(MF)) ErrorCBL("MF = " + conv(MF,par::fDP3) + "!", "", "MassFunction.cpp");
158 
159  if (m_fNL!=0) MF *= MF_correction(MASS, redshift, method_SS, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
160 
161  return MF;
162 }
163 
164 
165 // =====================================================================================
166 
167 
168 double cbl::cosmology::Cosmology::mass_function (const double Mass, const double Sigma, const double Dln_Sigma, const double redshift, const std::string model_MF, const bool store_output, const std::string output_root, const double Delta, const std::string interpType, const int norm, const double k_min, const double k_max, const double prec, const std::string method_SS, const std::string input_file, const bool is_parameter_file)
169 {
170  double fact = (m_unit) ? 1 : m_hh;
171  double MASS = Mass*fact;
172 
173  double MF = m_MF_generator(MASS, Sigma, Dln_Sigma, redshift, model_MF, Delta)*pow(fact, 4.);
174 
175  if (m_fNL!=0) MF *= MF_correction(MASS, redshift, method_SS, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
176 
177  return MF;
178 }
179 
180 
181 // =====================================================================================
182 
183 
184 double cbl::cosmology::Cosmology::mass_function (const double Mass, const double Sigma, const double Dln_Sigma, const double redshift, const double D_N, const std::string model_MF, const bool store_output, const std::string output_root, const double Delta, const std::string interpType, const int norm, const double k_min, const double k_max, const double prec, const std::string method_SS, const std::string input_file, const bool is_parameter_file)
185 {
186  double fact = (m_unit) ? 1 : m_hh;
187  double MASS = Mass*fact;
188 
189  double MF = m_MF_generator(MASS, Sigma, Dln_Sigma, redshift, D_N, model_MF, Delta)*pow(fact, 4.);
190 
191  if (m_fNL!=0) MF *= MF_correction(MASS, redshift, method_SS, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
192 
193  return MF;
194 }
195 
196 
197 // =====================================================================================
198 
199 
200 double cbl::cosmology::Cosmology::m_MF_generator (const double Mass, const double Sigma, const double Dln_Sigma, const double redshift, const std::string model_MF, const double Delta, const bool default_delta, const double delta_t)
201 {
202  double D_N = DN(redshift);
203 
204  return m_MF_generator(Mass, Sigma, Dln_Sigma, redshift, D_N, model_MF, Delta, default_delta, delta_t);
205 }
206 
207 
208 // =====================================================================================
209 
210 
211 double cbl::cosmology::Cosmology::m_MF_generator (const double Mass, const double Sigma, const double Dln_Sigma, const double redshift, const double D_N, const std::string model_MF, const double Delta, const bool default_delta, const double delta_t)
212 {
213  const double deltacz = (default_delta) ? deltac(redshift) : fabs(delta_t*D_N);
214  const double sigmaz = Sigma*D_N;
215 
216  const double RHO = rho_m(0., true);
217 
218  if (model_MF=="PS") // Press & Schechter
219  return sqrt(2./par::pi)*RHO/(Mass*Mass)*deltacz/sigmaz*fabs(Dln_Sigma)*exp(-(deltacz*deltacz)*0.5/(sigmaz*sigmaz));
220 
221  else if (model_MF=="ST") { // Sheth & Tormen
222  const double aa = 0.707;
223  const double qq = 0.3;
224  const double AA = 0.3222;
225  return AA*sqrt(2.*aa/par::pi)*RHO/(Mass*Mass)*deltacz/sigmaz*(1+pow(sigmaz/(sqrt(aa)*deltacz),2.*qq))*fabs(Dln_Sigma)*exp(-(aa*deltacz*deltacz)*0.5/(sigmaz*sigmaz));
226  }
227 
228  else if (model_MF=="Jenkins") // Jenkins et al. (2001)
229  return 0.315*exp(-pow(fabs(log(1./sigmaz)+0.61),3.8))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
230 
231  else if (model_MF=="Warren") { // Warren et al. (2006)
232  const double AA = 0.7234;
233  const double aa = 1.625;
234  const double bb = 0.2538;
235  const double cc = 1.1982;
236  return AA*(pow(sigmaz,-aa)+bb)*exp(-cc/(sigmaz*sigmaz))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
237  }
238 
239  else if (model_MF=="Reed") { // Reed et al. (2007)
240  const double cc = 1.08;
241  const double aa = 0.764/cc;
242  const double pp = 0.3;
243  const double AA = 0.3222;
244  const double Fatt = AA*sqrt(2.*aa/par::pi);
245  const double n_eff = -6.*Dln_Sigma-3.;
246  const double G1 = exp(-pow(log(1./sigmaz)-0.4,2)/0.72);
247  const double G2 = exp(-pow(log(1./sigmaz)-0.75,2)/0.08);
248  return Fatt*RHO/(Mass*Mass)*deltacz/sigmaz*(1+pow(sigmaz*sigmaz/(aa*deltacz*deltacz),pp)+0.6*G1+0.4*G2)*fabs(Dln_Sigma)*exp(-(cc*aa*deltacz*deltacz)*0.5/(sigmaz*sigmaz)-(0.03/pow(n_eff+3.,2)*pow(deltacz/sigmaz,0.6)));
249  }
250 
251  else if (model_MF=="Pan") { // Pan (2007)
252  const double alpha = 0.435; // check!!!
253  return 4.*alpha/sqrt(2.*par::pi)*RHO/(Mass*Mass)*deltacz/pow(sigmaz,2.*alpha)*fabs(Dln_Sigma)*exp(-(deltacz*deltacz)*0.5/pow(sigmaz,4.*alpha));
254  }
255 
256  else if (model_MF=="ShenH") { // halo MF, Shen et al. (2006) // check!!!
257  const double alpha = -0.55;
258  const double beta = -0.56;
259  const double ni = pow(deltacz/sigmaz,2);
260  const double ni_fni = sqrt(ni/(2.*par::pi))*exp(-ni*pow(1.-beta*pow(ni,alpha),2)*0.5)*(1.-beta/pow(ni,-alpha)*(1+alpha+(alpha*(alpha+1.))*0.5));
261  return ni_fni*RHO/(Mass*Mass)*fabs(Dln_Sigma)*2;
262  }
263 
264  else if (model_MF=="ShenF") { // filament MF, Shen et al. (2006) // check!!!
265  const double alpha = -0.28;
266  const double beta = -0.012;
267  const double ni = pow(deltacz/sigmaz,2);
268  const double ni_fni = sqrt(ni/(2.*par::pi))*exp(-ni*pow(1.-beta*pow(ni,alpha),2)*0.5)*(1.-beta/pow(ni,-alpha)*(1+alpha+(alpha*(alpha+1.))*0.5));
269  return ni_fni*RHO/(Mass*Mass)*fabs(Dln_Sigma)*2;
270  }
271 
272  else if (model_MF=="ShenS") { // sheet MF, Shen et al. (2006) // check!!!
273  const double alpha = -0.61;
274  const double beta = 0.45;
275  const double ni = pow(deltacz/sigmaz,2);
276  const double ni_fni = sqrt(ni/(2.*par::pi))*exp(-ni*pow(1.-beta*pow(ni,alpha),2)*0.5)*(1.-beta/pow(ni,-alpha)*(1+alpha+(alpha*(alpha+1.))*0.5));
277  return fabs(ni_fni*RHO/(Mass*Mass)*fabs(Dln_Sigma)*2); // check!!!
278  }
279 
280  else if (model_MF=="Tinker") { // Tinker et al. (2008)
281 
282  //if (redshift>2) WarningMsgCBL("the Tinker mass function has been tested for z<~2!", "m_MF_generator", "MassFunction.cpp");
283 
284  double A0, a0, b0, c0;
285 
286  if (Delta==200) {A0 = 0.186; a0 = 1.47; b0 = 2.57; c0 = 1.19;}
287  else if (Delta==300) {A0 = 0.200; a0 = 1.52; b0 = 2.25; c0 = 1.27;}
288  else if (Delta==400) {A0 = 0.212; a0 = 1.56; b0 = 2.05; c0 = 1.34;}
289  else if (Delta==600) {A0 = 0.218; a0 = 1.61; b0 = 1.87; c0 = 1.45;}
290  else if (Delta==800) {A0 = 0.248; a0 = 1.87; b0 = 1.59; c0 = 1.58;}
291  else if (Delta==1200) {A0 = 0.255; a0 = 2.13; b0 = 1.51; c0 = 1.80;}
292  else if (Delta==1600) {A0 = 0.260; a0 = 2.30; b0 = 1.46; c0 = 1.97;}
293  else if (Delta==2400) {A0 = 0.260; a0 = 2.53; b0 = 1.44; c0 = 2.24;}
294  else if (Delta==3200) {A0 = 0.260; a0 = 2.66; b0 = 1.41; c0 = 2.44;}
295  else {
296  A0 = (Delta<1600) ? 0.1*log10(Delta)-0.05 : 0.26;
297  a0 = 1.43+pow(max(log10(Delta)-2.3,0.),1.5); // check!!!
298  b0 = 1.+pow(log10(Delta)-1.6,-1.5);
299  c0 = 1.2+pow(max(log10(Delta)-2.35,0.),1.6); // check!!!
300  }
301 
302  const double alpha = pow(10.,-pow(max(0.75/log10(Delta/75.),0.),1.2)); // check!!!
303  const double AA = A0*pow(1.+redshift,-0.14);
304  const double aa = a0*pow(1.+redshift,-0.06);
305  const double bb = b0*pow(1.+redshift,-alpha);
306  const double cc = c0;
307 
308  return AA*(pow(sigmaz/bb,-aa)+1.)*exp(-cc/(sigmaz*sigmaz))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
309  }
310 
311  else if (model_MF=="Tinker08_Interp") { // Tinker et al. (2008)
312 
313  //if (redshift>2) WarningMsgCBL("the Tinker mass function has been tested for z<~2!", "m_MF_generator", "MassFunction.cpp");
314 
315  vector<double> _Delta = {200, 300, 400, 600, 800, 1200, 1600, 2400, 3200};
316  vector<double> _A0 = {1.858659e-01, 1.995973e-01, 2.115659e-01, 2.184113e-01, 2.480968e-01, 2.546053e-01, 2.600000e-01, 2.600000e-01, 2.600000e-01};
317  vector<double> _a0 = {1.466904, 1.521782, 1.559186, 1.614585, 1.869936, 2.128056, 2.301275, 2.529241, 2.661983};
318  vector<double> _b0 = {2.571104, 2.254217, 2.048674, 1.869559, 1.588649, 1.507134, 1.464374, 1.436827, 1.405210};
319  vector<double> _c0 = {1.193958, 1.270316, 1.335191, 1.446266, 1.581345, 1.795050, 1.965613, 2.237466, 2.439729};
320 
321  cbl::glob::FuncGrid interp_A0(_Delta, _A0, "Spline");
322  cbl::glob::FuncGrid interp_a0(_Delta, _a0, "Spline");
323  cbl::glob::FuncGrid interp_b0(_Delta, _b0, "Spline");
324  cbl::glob::FuncGrid interp_c0(_Delta, _c0, "Spline");
325 
326  const double alpha = pow(10.,-pow(max(0.75/log10(Delta/75.),0.),1.2)); // check!!!
327  const double AA = interp_A0(Delta)*pow(1.+redshift,-0.14);
328  const double aa = interp_a0(Delta)*pow(1.+redshift,-0.06);
329  const double bb = interp_b0(Delta)*pow(1.+redshift,-alpha);
330  const double cc = interp_c0(Delta);
331 
332  return AA*(pow(sigmaz/bb,-aa)+1.)*exp(-cc/(sigmaz*sigmaz))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
333  }
334 
335 
336  else if (model_MF=="Crocce") { // Crocce et al. (2010)
337  const double AA = 0.58*pow(1.+redshift,-0.13);
338  const double aa = 1.37*pow(1.+redshift,-0.15);
339  const double bb = 0.3*pow(1.+redshift,-0.084);
340  const double cc = 1.036*pow(1.+redshift,-0.024);
341  return AA*(pow(sigmaz,-aa)+bb)*exp(-cc*pow(sigmaz,-2))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
342  }
343 
344  else if (model_MF=="Angulo_FOF") // FoF MF, Angulo et al. (2012)
345  return 0.201*(pow(2.08/sigmaz,1.7)+1)*exp(-1.172/(sigmaz*sigmaz))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
346 
347  else if (model_MF=="Angulo_Sub") // SUBFIND MF, Angulo et al. (2012)
348  return 0.265*(pow(1.675/sigmaz,1.9)+1)*exp(-1.4/(sigmaz*sigmaz))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
349 
350  else if (model_MF=="Bhattacharya") { // Bhattacharya et al. (2011) // check!!!
351  const double aa = 0.788/pow(1+redshift,0.01);
352  const double qq = 1.795;
353  const double AA = 0.333/pow(1+redshift,0.11);
354  const double pp = 0.807;
355  const double nu = deltacz/sigmaz;
356  //if (redshift>=2) aa = 0.799/pow(1+redshift,0.024);
357  return RHO/(Mass*Mass)*fabs(Dln_Sigma)*AA*sqrt(2/par::pi)*exp(-nu*nu*0.5*aa)*(1+pow(1/(nu*nu*aa),pp))*pow(nu*sqrt(aa),qq);
358  }
359 
360  else if (model_MF=="Courtin") { // Courtin et al. (2010)
361  const double aa = 0.695;
362  const double qq = 0.1;
363  const double AA = 0.348;
364  return AA*sqrt(2.*aa/par::pi)*RHO/(Mass*Mass)*deltacz/sigmaz*(1+pow(sigmaz/(sqrt(aa)*deltacz),2.*qq))*fabs(Dln_Sigma)*exp(-(aa*deltacz*deltacz)*0.5/(sigmaz*sigmaz));
365  }
366 
367  else if (model_MF=="Manera") { // Manera et al. (2010)
368 
369  double aa, qq;
370  if (fabs(redshift-0.)<1.e-30) {
371  aa = 0.709;
372  qq = 0.248;
373  }
374  else if (fabs(redshift-0.5)<1.e-30) {
375  aa = 0.724;
376  qq = 0.241;
377  }
378  else
379  ErrorCBL("the Manera et al. mass function has been tested only for z=0 and 0.5!", "m_MF_generator", "MassFunction.cpp");
380 
381  return 0.3222*sqrt(2.*aa/par::pi)*RHO/(Mass*Mass)*deltacz/sigmaz*(1+pow(sigmaz/(sqrt(aa)*deltacz),2.*qq))*fabs(Dln_Sigma)*exp(-(aa*deltacz*deltacz)*0.5/(sigmaz*sigmaz));
382  }
383 
384  else if (model_MF=="Watson_FOF") { // FoF MF, Watson et al. (2012)
385  const double a = 0.282;
386  const double b = 1.406;
387  const double c = 2.163;
388  const double d = 1.210;
389  return a*(pow((b/sigmaz),c)+1)*exp(-d/(sigmaz*sigmaz))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
390  }
391 
392  else if (model_MF=="Watson_SOH") { // Spherical Overdensity halo MF, Watson et al. (2012)
393  const double Omega = OmegaM(redshift);
394  const double C = exp(0.023*(Delta/178-1));
395  const double d = -0.456*Omega-0.139;
396  const double p = 0.072;
397  const double q = 2.13;
398  const double Gamma = C*pow((Delta/178),d)*exp(p*(1-Delta/178)/pow(sigmaz,q));
399  double A = 0.194;
400  double alpha = 2.267;
401  double beta = 1.805;
402  double gamma = 1.287;
403  if (redshift>0 && redshift<6) {
404  A = Omega*(1.097*pow((1+redshift),-3.216)+0.074);
405  alpha = Omega*(3.136*pow((1+redshift),-3.056)+2.349);
406  beta = Omega*(5.907*pow((1+redshift),-3.599)+2.344);
407  gamma = 1.318;
408  }
409  if (redshift>=6) {
410  A = 0.563;
411  alpha = 0.874;
412  beta = 3.810;
413  gamma = 1.453;
414  }
415  return RHO/(Mass*Mass)*fabs(Dln_Sigma)*Gamma*A*(pow((beta/sigmaz),alpha)+1.)*exp(-gamma/(sigmaz*sigmaz));
416  }
417 
418  else if (model_MF=="Peacock") { // Peacock et al. (2007)
419  const double aa = 1.529;
420  const double bb = 0.704;
421  const double cc = 0.412;
422  const double nu = deltacz/sigmaz;
423  return nu*exp(-cc*nu*nu)/pow(1+aa*pow(nu,bb),2)*(bb*aa*pow(nu,bb-1)+2*cc*nu*(1+aa*pow(nu,bb)))*RHO/(Mass*Mass)*fabs(Dln_Sigma);
424  }
425 
426  else if (model_MF=="Despali_Z0") { // Despali et al. 2016
427  const double aa = 0.794; // 0.794 ± 0.005, z=0 fit
428  const double pp = 0.247; // 0.247 ± 0.009, z=0 fit
429  const double AA = 0.333; // 0.333 ± 0.001, z=0 fit
430  const double nu = pow(deltacz/sigmaz, 2);
431  const double nup = aa*nu;
432  return 2.*AA*(1.+pow(nup,-pp))*sqrt(nup/(2.*par::pi))*exp(-0.5*nup)*RHO/(Mass*Mass)*fabs(Dln_Sigma);
433  }
434 
435  else if (model_MF=="Despali_AllZ") { // Despali et al. 2016
436  const double aa = 0.7663; // 0.7663 ± 0.0013, all z fit
437  const double pp = 0.2579; // 0.2579 ± 0.0026, all z fit
438  const double AA = 0.3298; // 0.3298 ± 0.0003, all z fit
439  const double nu = pow(deltacz/sigmaz, 2);
440  const double nup = aa*nu;
441  return 2.*AA*(1.+pow(nup,-pp))*sqrt(nup/(2.*par::pi))*exp(-0.5*nup)*RHO/(Mass*Mass)*fabs(Dln_Sigma);
442  }
443 
444  else if (model_MF=="Despali_AllZAllCosmo") { // Despali et al. 2016
445  const double aa = 0.7689; // 0.7689 ± 0.0011, all z fit
446  const double pp = 0.2536; // 0.2536 ± 0.0026, all z fit
447  const double AA = 0.3295; // 0.3295 ± 0.0003, all z fit
448  const double nu = pow(deltacz/sigmaz, 2);
449  const double nup = aa*nu;
450  return 2.*AA*(1.+pow(nup,-pp))*sqrt(nup/(2.*par::pi))*exp(-0.5*nup)*RHO/(Mass*Mass)*fabs(Dln_Sigma);
451  }
452 
453  else if (model_MF=="Despali_HighM") { // Despali et al. 2016
454  const double aa = 0.8199; // 0.8199 ± 0.0010, M>3e13
455  const double pp = 0.0; // 0, M>3e13
456  const double AA = 0.3141; // 0.3141 ± 0.0006, M>3e13
457  const double nu = pow(deltacz/sigmaz, 2);
458  const double nup = aa*nu;
459  return 2.*AA*(1.+pow(nup,-pp))*sqrt(nup/(2.*par::pi))*exp(-0.5*nup)*RHO/(Mass*Mass)*fabs(Dln_Sigma);
460  }
461 
462  else return ErrorCBL("model_MF not allowed!", "m_MF_generator", "MassFunction.cpp");
463 }
464 
465 
466 // =====================================================================================
467 
468 
469 double cbl::cosmology::Cosmology::n_haloes (const double Mass_min, const double Mass_max, const double z_min, const double z_max, const bool angle_rad, const std::string model_MF, const std::string method_SS, const bool store_output, const std::string output_root, const double Delta, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file)
470 {
471 
472  // ---------- read the grid file ----------
473 
474  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root,interpType, k_max, input_file, is_parameter_file);
475 
476  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
477 
478  double Mass, Sigma, Dln_Sigma;
479  vector<double> mass, sigma, dlnsigma;
480  while (fin >>Mass>>Sigma>>Dln_Sigma) {
481  if (Mass_min<Mass && Mass<Mass_max) {
482  mass.push_back(Mass);
483  sigma.push_back(Sigma);
484  dlnsigma.push_back(Dln_Sigma);
485  }
486  }
487 
488 
489  // ---------- compute the number of haloes ----------
490 
491  int step_z = 100;
492  double delta_z = (z_max-z_min)/step_z;
493  double redshift = z_min;
494  double N_haloes = 0.;
495 
496  for (int Red=0; Red<step_z; Red++) {
497 
498  double Int = 0.;
499  for (unsigned int k=0; k<mass.size()-1; k++)
500  Int += mass_function(mass[k], sigma[k], dlnsigma[k], redshift, model_MF, store_output, output_root, Delta)*dV_dZdOmega(redshift, angle_rad)*(mass[k+1]-mass[k]);
501 
502  N_haloes += Int*delta_z;
503 
504  redshift += delta_z;
505  }
506 
507  return N_haloes;
508 }
509 
510 
511 // =====================================================================================
512 
513 
514 double cbl::cosmology::Cosmology::n_haloes (const double Mass_min, const double Mass_max, const double Volume, const double redshift, const std::string model_MF, const std::string method_SS, const int nbin_mass, const bool store_output, const std::string output_root, const double Delta, const std::string interpType, const int norm, const double k_min, const double k_max, const double prec, const std::string input_file, const bool is_parameter_file, const bool default_delta, const double delta_t)
515 {
516  glob::STR_MF mass_function_par;
517 
518  mass_function_par.redshift = redshift;
519  mass_function_par.model_MF = model_MF ;
520  mass_function_par.method_SS = method_SS;
521  mass_function_par.store_output = store_output;
522  mass_function_par.output_root = output_root;
523  mass_function_par.Delta = Delta;
524  mass_function_par.interpType = interpType;
525  mass_function_par.norm = norm;
526  mass_function_par.k_min = k_min;
527  mass_function_par.k_max = k_max;
528  mass_function_par.prec = prec;
529  mass_function_par.input_file = input_file;
530  mass_function_par.is_parameter_file = is_parameter_file;
531  mass_function_par.default_delta = default_delta;
532  mass_function_par.delta_t = delta_t;
533 
534  auto pp = make_shared<glob::STR_MF>(mass_function_par);
535 
536  if (nbin_mass==0) {
537  function<double(double)> func = bind(&Cosmology::m_mass_function, this, std::placeholders::_1, pp);
538  return Volume*wrapper::gsl::GSL_integrate_qag(func, Mass_min, Mass_max, prec);
539  }
540 
541  else if (nbin_mass>0)
542  {
543  vector<double> mass = logarithmic_bin_vector(nbin_mass, Mass_min, Mass_max);
544  vector<double> dndm(nbin_mass);
545 
546  for (int i=0; i<nbin_mass; i++)
547  dndm[i] = m_mass_function(mass[i], pp);
548 
549  glob::FuncGrid mass_function_interpolator(mass, dndm, interpType);
550 
551  return Volume*mass_function_interpolator.integrate_qag(Mass_min, Mass_max, prec);
552  }
553 
554  else
555  ErrorCBL("nbinmass must be >=0", "n_haloes", "MassFunction.cpp");
556 
557  return 0;
558 }
559 
560 
561 // =====================================================================================
562 
563 
564 double cbl::cosmology::Cosmology::MhaloMin (const int n_halo, const double Area, const bool angle_rad, const double z_min, const double z_max, const double Mmax, const double lgM1_guess, const double lgM2_guess, const std::string model_MF, const std::string method_SS, const bool store_output, const std::string output_root, const double Delta, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file) const
565 {
566  cbl::classfunc::func_MhaloMin func(m_Omega_matter, m_Omega_baryon, m_Omega_neutrinos, m_massless_neutrinos, m_massive_neutrinos, m_Omega_DE, m_Omega_radiation, m_hh, m_scalar_amp, m_scalar_pivot, m_n_spec, m_w0, m_wa, m_fNL, m_type_NG, m_tau, m_model, m_unit, n_halo, Area, angle_rad, z_min, z_max, Mmax, model_MF, method_SS, store_output, output_root, Delta, interpType, k_max, input_file, is_parameter_file);
567 
568  function<double(double) > ff = bind(&cbl::classfunc::func_MhaloMin::operator(), func, std::placeholders::_1);
569  double prec = 0.0001;
570  double lgM = wrapper::gsl::GSL_root_brent (ff, 0., lgM1_guess, lgM2_guess, prec);
571 
572  return pow(10., lgM);
573 }
574 
575 
576 // =====================================================================================
577 
578 
579 double cbl::cosmology::Cosmology::unevolved_mass_function (const double mass_accr) const
580 {
581  const double aa = -0.8;
582  const double mm = mass_accr/aa; // mass_accr = m/M
583  const double yn = log(0.3); // log(0.21) check!!!
584 
585  return aa*log(-mm)+log(exp(2.*par::pi*pow(mm,3.)))+yn;
586 }
587 
588 
589 // =====================================================================================
590 
591 
592 double cbl::cosmology::Cosmology::n_haloes_selection_function (const double Mass_min, const double Mass_max, const double z_min, const double z_max, const bool angle_rad, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector<int> cols, const bool store_output, const std::string output_root, const double Delta, const bool isDelta_critical, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file)
593 {
594 
595  // ---------- read the grid file ----------
596 
597  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root,interpType, k_max, input_file, is_parameter_file);
598 
599  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
600 
601  double Mass, Sigma, Dln_Sigma;
602  vector<double> mass, sigma, dlnsigma;
603  while (fin >>Mass>>Sigma>>Dln_Sigma) {
604  if (Mass_min<Mass && Mass<Mass_max) {
605  mass.push_back(Mass);
606  sigma.push_back(Sigma);
607  dlnsigma.push_back(Dln_Sigma);
608  }
609  }
610  fin.clear(); fin.close();
611 
612  glob::FuncGrid interp_sigma(mass, sigma, "Spline");
613  glob::FuncGrid interp_dlnsigma(mass, dlnsigma, "Spline");
614 
615 
616  // ---------- read the selection function ----------
617 
618  vector<double> mass_SF, redshift_SF;
619  vector<vector<double>> selection_function;
620 
621  read_matrix(selection_function_file, mass_SF, redshift_SF, selection_function, cols);
622  function<double(const double, const double)> interp_SF = bind(interpolated_2D, std::placeholders::_1, std::placeholders::_2, mass_SF, redshift_SF, selection_function, "Linear");
623 
624 
625  // ---------- compute the number of haloes ----------
626 
627  vector<vector<double>> limits = {{Mass_min, Mass_max}, {z_min, z_max}};
628 
629  int step_z = 100;
630  double delta_z = (z_max-z_min)/step_z;
631  double redshift = z_min;
632  double N_haloes = 0.;
633  //double norm_SF = 0.;
634 
635  for (int Red=0; Red<step_z; Red++) {
636  double Int = 0.;
637  for (unsigned int k=0; k<mass.size()-1; k++) {
638  double DD = (isDelta_critical) ? Delta/OmegaM(redshift) : Delta;
639  double SF = interp_SF(mass[k], redshift);
640  Int += SF*mass_function(mass[k], sigma[k], dlnsigma[k], redshift, model_MF, store_output, output_root, DD)*dV_dZdOmega(redshift, angle_rad)*(mass[k+1]-mass[k]);
641  }
642 
643  N_haloes += Int*delta_z;
644 
645  redshift += delta_z;
646  }
647 
648  return N_haloes;
649 }
650 
651 
652 // =====================================================================================
653 
654 
655 std::vector<double> cbl::cosmology::Cosmology::mass_function (const std::vector<double> mass, const double z_min, const double z_max, const std::string model_MF, const std::string method_SS, const bool store_output, const std::string output_root, const double Delta, const bool isDelta_critical, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file)
656 {
657  vector<double> MF(mass.size());
658 
659  // ---------- read the grid file ----------
660 
661  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root,interpType, k_max, input_file, is_parameter_file);
662 
663  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
664 
665  double Mass, Sigma, Dln_Sigma;
666  vector<double> _mass, _sigma, _dlnsigma;
667  while (fin >>Mass>>Sigma>>Dln_Sigma) {
668  _mass.push_back(Mass);
669  _sigma.push_back(Sigma);
670  _dlnsigma.push_back(Dln_Sigma);
671  }
672  fin.clear(); fin.close();
673 
674  glob::FuncGrid interp_sigma(_mass, _sigma, "Spline");
675  glob::FuncGrid interp_dlnsigma(_mass, _dlnsigma, "Spline");
676 
677  // ---------- compute the number of haloes ----------
678 
679  //int step_z = 100;
680  //double delta_z = (z_max-z_min)/step_z;
681 
682  double VV = Volume(z_min, z_max, 1.);
683 
684  for (unsigned int k=0; k<mass.size(); k++) {
685  double sigma = interp_sigma(mass[k]);
686  double dlnsigma = interp_dlnsigma(mass[k]);
687 
688  auto func = [&] (const double redshift)
689  {
690  double DD = (isDelta_critical) ? Delta/OmegaM(redshift) : Delta;
691  return mass_function(mass[k], sigma, dlnsigma, redshift, model_MF, store_output, output_root, DD)*dV_dZdOmega(redshift, false);
692  };
693  MF[k] = wrapper::gsl::GSL_integrate_qag(func, z_min, z_max)/VV;
694  }
695 
696  return MF;
697 }
698 
699 
700 // =====================================================================================
701 
702 
703 std::vector<double> cbl::cosmology::Cosmology::mass_function_selection_function_vector (const std::vector<double> mass, const double z_min, const double z_max, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector<int> cols, const bool store_output, const std::string output_root, const double Delta, const bool isDelta_critical, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file)
704 {
705  vector<double> MF(mass.size());
706 
707  // ---------- read the grid file ----------
708 
709  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root,interpType, k_max, input_file, is_parameter_file);
710 
711  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
712 
713  double Mass, Sigma, Dln_Sigma;
714  vector<double> _mass, _sigma, _dlnsigma;
715  while (fin >>Mass>>Sigma>>Dln_Sigma) {
716  _mass.push_back(Mass);
717  _sigma.push_back(Sigma);
718  _dlnsigma.push_back(Dln_Sigma);
719  }
720  fin.clear(); fin.close();
721 
722  glob::FuncGrid interp_sigma(_mass, _sigma, "Spline");
723  glob::FuncGrid interp_dlnsigma(_mass, _dlnsigma, "Spline");
724 
725  // ---------- read the selection function ----------
726 
727  vector<double> mass_SF, redshift_SF;
728  vector<vector<double>> selection_function;
729 
730  read_matrix(selection_function_file, mass_SF, redshift_SF, selection_function, cols);
731  function<double(const double, const double)> interp_SF = bind(interpolated_2D, std::placeholders::_1, std::placeholders::_2, mass_SF, redshift_SF, selection_function, "Linear");
732 
733  // ---------- compute the number of haloes ----------
734 
735  //int step_z = 100;
736  //double delta_z = (z_max-z_min)/step_z;
737 
738  double VV = Volume(z_min, z_max, 1.);
739 
740  for (unsigned int k=0; k<mass.size(); k++) {
741  double sigma = interp_sigma(mass[k]);
742  double dlnsigma = interp_dlnsigma(mass[k]);
743 
744  auto func = [&] (const double redshift)
745  {
746  double SF = interp_SF(mass[k], redshift);
747  double DD = (isDelta_critical) ? Delta/OmegaM(redshift) : Delta;
748  return SF*mass_function(mass[k], sigma, dlnsigma, redshift, model_MF, store_output, output_root, DD)*dV_dZdOmega(redshift, false);
749  };
750 
751  MF[k] = wrapper::gsl::GSL_integrate_qag(func, z_min, z_max)/VV;
752  }
753 
754  return MF;
755 }
756 
757 
758 // =====================================================================================
759 
760 
761 std::vector<double> cbl::cosmology::Cosmology::redshift_distribution_haloes (const double z_min, const double z_max, const int step_z, const double Area_degrees, const double Mass_min, const double Mass_max, const std::string model_MF, const std::string method_SS, const bool store_output, const std::string output_root, const double Delta, const bool isDelta_critical, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file)
762 {
763 
764  // ---------- read the grid file ----------
765 
766  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root,interpType, k_max, input_file, is_parameter_file);
767 
768  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
769 
770  double Mass, Sigma, Dln_Sigma;
771  vector<double> mass, sigma, dlnsigma;
772  while (fin >>Mass>>Sigma>>Dln_Sigma) {
773  if (Mass_min<Mass && Mass<Mass_max) {
774  mass.push_back(Mass);
775  sigma.push_back(Sigma);
776  dlnsigma.push_back(Dln_Sigma);
777  }
778  }
779  fin.clear(); fin.close();
780 
781  glob::FuncGrid interp_sigma(mass, sigma, "Spline");
782  glob::FuncGrid interp_dlnsigma(mass, dlnsigma, "Spline");
783 
784  // ---------- compute the number of haloes ----------
785 
786  vector<double> redshift_distribution(step_z, 0);
787  double delta_z = (z_max-z_min)/step_z;
788  double redshift = z_min;
789 
790  for (int Red=0; Red<step_z; Red++) {
791 
792  vector<vector<double>> integration_limits(2);
793  integration_limits[0] = {redshift, redshift+delta_z};
794  integration_limits[1] = {Mass_min, Mass_max};
795 
796  auto integrand = [&] (const vector<double> parameters)
797  {
798  double redshift = parameters[0];
799  double mass = parameters[1];
800 
801  double DD = (isDelta_critical) ? Delta/OmegaM(redshift) : Delta;
802  double sigma = interp_sigma(mass);
803  double dlnsigma = interp_dlnsigma(mass);
804  return Area_degrees*dV_dZdOmega(redshift, false)*mass_function(mass, sigma, dlnsigma, redshift, model_MF, store_output, output_root, DD);
805  };
806 
807  cbl::wrapper::cuba::CUBAwrapper CubaIntegrand(integrand, 2);
808 
809  redshift_distribution[Red] = CubaIntegrand.IntegrateVegas(integration_limits);
810  redshift += delta_z;
811  }
812 
813  return redshift_distribution;
814 }
815 
816 
817 // =====================================================================================
818 
819 
820 std::vector<double> cbl::cosmology::Cosmology::redshift_distribution_haloes_selection_function (const std::vector<double> redshift, const double Area_degrees, const double Mass_min, const double Mass_max, const string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector<int> cols, const bool store_output, const std::string output_root, const double Delta, const bool isDelta_critical, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file)
821 {
822 
823  // ---------- read the grid file ----------
824 
825  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root,interpType, k_max, input_file, is_parameter_file);
826 
827  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
828 
829  double Mass, Sigma, Dln_Sigma;
830  vector<double> mass, sigma, dlnsigma;
831  while (fin >>Mass>>Sigma>>Dln_Sigma) {
832  if (Mass_min<Mass && Mass<Mass_max) {
833  mass.push_back(Mass);
834  sigma.push_back(Sigma);
835  dlnsigma.push_back(Dln_Sigma);
836  }
837  }
838  fin.clear(); fin.close();
839 
840  glob::FuncGrid interp_sigma(mass, sigma, "Spline");
841  glob::FuncGrid interp_dlnsigma(mass, dlnsigma, "Spline");
842 
843 
844  // ---------- read the selection function ----------
845 
846  vector<double> mass_SF, redshift_SF;
847  vector<vector<double>> selection_function;
848 
849  read_matrix(selection_function_file, mass_SF, redshift_SF, selection_function, cols);
850  //function<double(const double, const double)> interp_SF = bind(interpolated_2D, std::placeholders::_1, std::placeholders::_2, mass_SF, redshift_SF, selection_function, "Linear");
851 
852  glob::FuncGrid2D interp_SF(mass_SF, redshift_SF, selection_function, "Linear");
853 
854 
855  // ---------- compute the number of haloes ----------
856 
857  vector<double> redshift_distribution(redshift.size(), 0.);
858 
859  for (size_t zz=0; zz<redshift.size()-1; ++zz) {
860 
861  vector<vector<double>> integration_limits(2);
862  integration_limits[0] = {redshift[zz], redshift[zz+1]};
863  integration_limits[1] = {Mass_min, Mass_max};
864 
865  auto integrand = [&] (const vector<double> parameters)
866  {
867  const double redshift = parameters[0];
868  const double mass = parameters[1];
869 
870  const double DD = (isDelta_critical) ? Delta/OmegaM(redshift) : Delta;
871  const double sigma = interp_sigma(mass);
872  const double dlnsigma = interp_dlnsigma(mass);
873  const double SF = interp_SF(mass, redshift);
874 
875  return SF*Area_degrees*dV_dZdOmega(redshift, false)*mass_function(mass, sigma, dlnsigma, redshift, model_MF, store_output, output_root, DD);
876  };
877 
878  wrapper::cuba::CUBAwrapper CubaIntegrand(integrand, 2);
879  const double distr = CubaIntegrand.IntegrateVegas(integration_limits);
880  redshift_distribution[zz] = (distr!=distr) ? 0. : distr;
881  }
882 
883  return redshift_distribution;
884 }
885 
886 
887 // =====================================================================================
888 
889 
890 double cbl::cosmology::Cosmology::mean_redshift_haloes_selection_function (const double z_min, const double z_max, const double Mass_min, const double Mass_max, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector<int> cols, const bool store_output, const std::string output_root, const double Delta, const bool isDelta_critical, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file)
891 {
892 
893  // ---------- read the grid file ----------
894 
895  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
896 
897  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
898 
899  double Mass, Sigma, Dln_Sigma;
900  vector<double> mass, sigma, dlnsigma;
901  while (fin >>Mass>>Sigma>>Dln_Sigma) {
902  if (Mass_min<Mass && Mass<Mass_max) {
903  mass.push_back(Mass);
904  sigma.push_back(Sigma);
905  dlnsigma.push_back(Dln_Sigma);
906  }
907  }
908  fin.clear(); fin.close();
909 
910  glob::FuncGrid interp_sigma(mass, sigma, "Spline");
911  glob::FuncGrid interp_dlnsigma(mass, dlnsigma, "Spline");
912 
913  // ---------- read the selection function ----------
914 
915  vector<double> mass_SF, redshift_SF;
916  vector<vector<double>> selection_function;
917 
918  read_matrix(selection_function_file, mass_SF, redshift_SF, selection_function, cols);
919 
920  glob::FuncGrid2D interp_SF( mass_SF, redshift_SF, selection_function, "Linear");
921 
922  // ---------- compute the mean redshift ----------
923 
924  auto integrand_num = [&] (const vector<double> parameters) {
925  double redshift = parameters[0];
926  double mass = parameters[1];
927  double sigma = interp_sigma(mass);
928  double dlnsigma = interp_dlnsigma(mass);
929  double SF = interp_SF(mass, redshift);
930 
931  double DD = (isDelta_critical) ? Delta/OmegaM(redshift) : Delta;
932  double MF = mass_function(mass, sigma, dlnsigma, redshift, model_MF, store_output, output_root, DD);
933  return redshift*MF*SF*dV_dZdOmega(redshift, false);
934  };
935 
936  auto integrand_denom = [&] (const vector<double> parameters) {
937 
938  double redshift = parameters[0];
939  double mass = parameters[1];
940  double sigma = interp_sigma(mass);
941  double dlnsigma = interp_dlnsigma(mass);
942  double SF = interp_SF(mass, redshift);
943 
944  double DD = (isDelta_critical) ? Delta/OmegaM(redshift) : Delta;
945  double MF = mass_function(mass, sigma, dlnsigma, redshift, model_MF, store_output, output_root, DD);
946  return MF*SF*dV_dZdOmega(redshift, false);
947  };
948 
949  cbl::wrapper::cuba::CUBAwrapper IntegrandNum(integrand_num, 2);
950  cbl::wrapper::cuba::CUBAwrapper IntegrandDenom(integrand_denom, 2);
951 
952  vector<vector<double>> limits = {{z_min, z_max}, {Mass_min, Mass_max}};
953 
954  return IntegrandNum.IntegrateVegas(limits)/IntegrandDenom.IntegrateVegas(limits);
955 
956 }
957 
958 
959 // =====================================================================================
960 
961 
962 double cbl::cosmology::Cosmology::converted_mass (const double mass, const cosmology::Cosmology cosmology, const double redshift, const double redshift_source) const
963 {
964  const double Dds1 = cosmology.D_A(redshift, redshift_source);
965  const double Dds2 = this->D_A(redshift, redshift_source);
966 
967  const double Ds1 = cosmology.D_A(redshift);
968  const double Ds2 = this->D_A(redshift);
969 
970  const double H1 = cosmology.HH(redshift);
971  const double H2 = this->HH(redshift);
972 
973  return mass*pow(Dds1/Ds1, 1.5)*H1*pow(Dds2/Ds2, -1.5)/H2;
974 }
Cosmology.h
The class Cosmology.
cbl::Path
The class Path.
Definition: Path.h:145
cbl::Path::DirCosmo
std::string DirCosmo()
get the directory where the CosmoBolognaLbi are stored
Definition: Path.h:178
cbl::cosmology::Cosmology
The class Cosmology.
Definition: Cosmology.h:277
cbl::cosmology::Cosmology::unevolved_mass_function
double unevolved_mass_function(const double mass_accr) const
the unevolved mass function
Definition: MassFunction.cpp:579
cbl::cosmology::Cosmology::redshift_distribution_haloes_selection_function
std::vector< double > redshift_distribution_haloes_selection_function(const std::vector< double > redshift, const double Area_degrees, const double Mass_min, const double Mass_max, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector< int > column={}, const bool store_output=true, const std::string output_root="test", const double Delta=200, const bool isDelta_critical=false, const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true)
redshift distribution of dark matter haloes, given a selection function
Definition: MassFunction.cpp:820
cbl::cosmology::Cosmology::converted_mass
double converted_mass(const double mass, const cosmology::Cosmology cosmology, const double redshift, const double redshift_source=-1.) const
convert a cluster mass estimated in a different cosmology
Definition: MassFunction.cpp:962
cbl::cosmology::Cosmology::MhaloMin
double MhaloMin(const int n_halo, const double Area, const bool angle_rad, const double z_min, const double z_max, const double Mmax, const double lgM1_guess, const double lgM2_guess, const std::string model_MF, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const double Delta=200, const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true) const
minimum halo mass, given the number of haloes in a given region of sky
Definition: MassFunction.cpp:564
cbl::cosmology::Cosmology::n_haloes
double n_haloes(const double Mass_min, const double Mass_max, const double z_min, const double z_max, const bool angle_rad, const std::string model_MF, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const double Delta=200, const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true)
number of dark matter haloes per steradian or square degree, for a given redshift range
Definition: MassFunction.cpp:469
cbl::cosmology::Cosmology::HH
double HH(const double redshift=0.) const
the Hubble function
Definition: Cosmology.cpp:570
cbl::cosmology::Cosmology::m_MF_generator
double m_MF_generator(const double Mass, const double Sigma, const double Dln_Sigma, const double redshift, const std::string model_MF, const double Delta=200., const bool default_delta=true, const double delta_t=1.686)
auxiliary function to compute the mass function
Definition: MassFunction.cpp:200
cbl::cosmology::Cosmology::D_A
double D_A(const double redshift) const
the angular diameter distance at a given redshift
Definition: Cosmology.cpp:848
cbl::cosmology::Cosmology::mass_function_fast
double mass_function_fast(const double Mass, const double redshift, const std::string model_MF, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const double Delta=200., const std::string interpType="Linear", const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string input_file=par::defaultString, const bool is_parameter_file=true)
the mass function of dark matter haloes (filaments and sheets) computed quickly using a grid
Definition: MassFunction.cpp:130
cbl::cosmology::Cosmology::mass_function_selection_function_vector
std::vector< double > mass_function_selection_function_vector(const std::vector< double > mass, const double z_min, const double z_max, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector< int > column={}, const bool store_output=true, const std::string output_root="test", const double Delta=200, const bool isDelta_critical=false, const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true)
mass function given a selection function
Definition: MassFunction.cpp:703
cbl::cosmology::Cosmology::mass_function_fR
double mass_function_fR(const double Mass, const double redshift, const std::string model_MF, const double f_R0=0., const bool store_output=true, const double Delta=200., const std::string interpType="Linear", const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string input_file=par::defaultString, const bool is_parameter_file=true, const bool default_delta=true, const double delta_t=1.686)
the mass function of dark matter haloes in f(R) cosmologies (see Hu & Sawicki 2007) computed with the...
Definition: MassFunction.cpp:67
cbl::cosmology::Cosmology::m_mass_function
double m_mass_function(const double Mass, std::shared_ptr< void > mass_function_params)
auxiliary function to compute the mass function of dark matter haloes (filaments and sheets)
Definition: MassFunction.cpp:105
cbl::cosmology::Cosmology::mass_function
double mass_function(const double Mass, const double redshift, const std::string model_MF, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const double Delta=200., const std::string interpType="Linear", const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string input_file=par::defaultString, const bool is_parameter_file=true, const bool default_delta=true, const double delta_t=1.686)
the mass function of dark matter haloes (filaments and sheets)
Definition: MassFunction.cpp:45
cbl::cosmology::Cosmology::mean_redshift_haloes_selection_function
double mean_redshift_haloes_selection_function(const double z_min, const double z_max, const double Mass_min, const double Mass_max, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector< int > column={}, const bool store_output=true, const std::string output_root="test", const double Delta=200, const bool isDelta_critical=false, const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true)
the mean redshift of a dark matter haloe sample, given a selection function
Definition: MassFunction.cpp:890
cbl::cosmology::Cosmology::redshift_distribution_haloes
std::vector< double > redshift_distribution_haloes(const double z_min, const double z_max, const int step_z, const double Area_degrees, const double Mass_min, const double Mass_max, const std::string model_MF, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const double Delta=200, const bool isDelta_critical=false, const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true)
redshift distribution of dark matter haloes
Definition: MassFunction.cpp:761
cbl::cosmology::Cosmology::n_haloes_selection_function
double n_haloes_selection_function(const double Mass_min, const double Mass_max, const double z_min, const double z_max, const bool angle_rad, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector< int > column={}, const bool store_output=true, const std::string output_root="test", const double Delta=200, const bool isDelta_critical=false, const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true)
number of dark matter haloes per steradian or square degree, for a given redshift range and with sele...
Definition: MassFunction.cpp:592
cbl::glob::FuncGrid2D
The class FuncGrid2D.
Definition: FuncGrid.h:302
cbl::glob::FuncGrid
The class FuncGrid.
Definition: FuncGrid.h:55
cbl::glob::FuncGrid::integrate_qag
double integrate_qag(const double a, const double b, const double rel_err=1.e-2, const double abs_err=1.e-6, const int limit_size=1000, const int rule=6)
compute the definite integral with GSL qag method
Definition: FuncGrid.cpp:187
cbl::wrapper::cuba::CUBAwrapper
The class CUBAwrapper.
Definition: CUBAwrapper.h:187
cbl::wrapper::cuba::CUBAwrapper::IntegrateVegas
double IntegrateVegas(std::vector< std::vector< double >> integration_limits, const bool parallelize=true)
integrate using the Vegas routine
Definition: CUBAwrapper.cpp:103
cbl::par::fDP8
static const char fDP8[]
conversion symbol for: double -> std::string
Definition: Constants.h:151
cbl::par::fDP3
static const char fDP3[]
conversion symbol for: double -> std::string
Definition: Constants.h:136
cbl::par::fINT
static const char fINT[]
conversion symbol for: int -> std::string
Definition: Constants.h:121
cbl::par::pi
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
cbl::par::alpha
static const double alpha
: the fine-structure constant
Definition: Constants.h:233
cbl::par::cc
static const double cc
: the speed of light in vacuum (the value is exact) [km/sec]
Definition: Constants.h:221
cbl::modelling::numbercounts::mass_function
double mass_function(const double mass, cosmology::Cosmology cosmology, const double redshift, const std::string model_MF, const bool store_output, const double Delta, const bool isDelta_critical, const cbl::glob::FuncGrid interp_Pk, const double kmax)
compute the mass function
Definition: ModelFunction_NumberCounts.cpp:145
cbl::wrapper::gsl::GSL_root_brent
double GSL_root_brent(gsl_function Func, const double low_guess, const double up_guess, const double rel_err=1.e-3, const double abs_err=0)
function to find roots using GSL qag method
Definition: GSLwrapper.cpp:430
cbl::wrapper::gsl::GSL_integrate_qag
double GSL_integrate_qag(gsl_function Func, const double a, const double b, const double rel_err=1.e-3, const double abs_err=0, const int limit_size=1000, const int rule=6)
integral, computed using the GSL qag method
Definition: GSLwrapper.cpp:180
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
cbl::conv
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
cbl::read_matrix
void read_matrix(const std::string file_matrix, std::vector< double > &xx, std::vector< double > &yy, std::vector< std::vector< double >> &matrix, const std::vector< int > col={})
read a matrix from file
Definition: Func.cpp:612
cbl::Mass
T Mass(const T RR, const T Rho)
the mass of a sphere of a given radius and density
Definition: Func.h:1243
cbl::logarithmic_bin_vector
std::vector< T > logarithmic_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with logarithmically spaced values
Definition: Kernel.h:1621
cbl::checkIO
void checkIO(const std::ifstream &fin, const std::string file="NULL")
check if an input file can be opened
Definition: Kernel.cpp:160
cbl::interpolated_2D
double interpolated_2D(const double _x1, const double _x2, const std::vector< double > x1, const std::vector< double > x2, const std::vector< std::vector< double >> yy, const std::string type)
2D interpolation
Definition: Func.cpp:502
cbl::ErrorCBL
int ErrorCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_)
throw an exception: it is used for handling exceptions inside the CosmoBolognaLib
Definition: Kernel.h:780
cbl::interpolated
double interpolated(const double _xx, const std::vector< double > xx, const std::vector< double > yy, const std::string type)
1D interpolation
Definition: Func.cpp:445
cbl::Sigma
double Sigma(const std::vector< double > vect)
the standard deviation of a std::vector
Definition: Func.cpp:925