CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Bias.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 
34 #include "Cosmology.h"
35 
36 using namespace std;
37 
38 using namespace cbl;
39 
40 
41 // =====================================================================================
42 
43 
44 double cbl::cosmology::Cosmology::bias_halo (const double Mass, const double redshift, const std::string author, const std::string method_SS, const bool store_output, const std::string output_root, const std::string interpType, const double Delta, const double kk, const int norm, const double k_min, const double k_max, const double prec, const std::string input_file, const bool is_parameter_file)
45 {
46  const double SSS = sigma2M(Mass, method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
47  const double Sigma = sqrt(SSS);
48 
49  double bias = m_bias_halo_generator(Sigma, redshift, author, Delta);
50 
51  if (m_fNL!=0)
52  bias += bias_correction(kk, Mass, method_SS, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)*SSS*pow(bias-1, 2);
53 
54  return bias;
55 }
56 
57 
58 // =====================================================================================
59 
60 
61 double cbl::cosmology::Cosmology::bias_halo (const double Mass, const double Sigma, const double redshift, const std::string model_bias, const bool store_output, const std::string output_root, const std::string interpType, const double Delta, const double kk, 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)
62 {
63  double bias = m_bias_halo_generator(Sigma, redshift, model_bias, Delta);
64 
65  if (m_fNL!=0)
66  bias += bias_correction(kk, Mass, method_SS, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)*Sigma*pow(bias-1, 2); // check!!!
67 
68  return bias;
69 }
70 
71 
72 // =====================================================================================
73 
74 
75 double cbl::cosmology::Cosmology::bias_halo (const double Mass, const double Sigma, const double redshift, const double DN, const std::string model_bias, const bool store_output, const std::string output_root, const std::string interpType, const double Delta, const double kk, 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)
76 {
77  double bias = m_bias_halo_generator(Sigma, redshift, DN, model_bias, Delta);
78 
79  if (m_fNL!=0)
80  bias += bias_correction(kk, Mass, method_SS, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)*Sigma*pow(bias-1, 2); // check!!!
81 
82  return bias;
83 }
84 
85 
86 // =====================================================================================
87 
88 
89 double cbl::cosmology::Cosmology::m_bias_halo_generator (const double Sigma, const double redshift, const std::string author, const double Delta) const
90 {
91  const double D_N = DN(redshift);
92 
93  return m_bias_halo_generator(Sigma, redshift, D_N, author, Delta);
94 }
95 
96 
97 // =====================================================================================
98 
99 
100 double cbl::cosmology::Cosmology::m_bias_halo_generator (const double Sigma, const double redshift, const double D_N, const std::string author, const double Delta) const
101 {
102  const double deltacz = deltac(redshift);
103  const double sigmaz = Sigma*D_N;
104 
105  double bias = -1000.;
106 
107  if (author=="ST99") {
108  double aa = 0.707;
109  double pp = 0.3;
110  double ni = pow(deltacz/sigmaz, 2);
111  bias = 1.+(aa*ni-1.)/deltacz+(2.*pp/deltacz)/(1.+pow(aa*ni,pp));
112  }
113 
114  else if (author=="SMT01") {
115  double aa = 0.707;
116  double bb = 0.5;
117  double cc = 0.6;
118  double ni = deltacz/sigmaz;
119  bias = 1.+1./(sqrt(aa)*deltacz)*(sqrt(aa)*aa*pow(ni,2.)+sqrt(aa)*bb*pow(aa*pow(ni,2.),1.-cc)-pow(aa*pow(ni,2.),cc)/(pow(aa*pow(ni,2.),cc)+bb*(1.-cc)*(1.-cc*0.5)));
120  }
121 
122  else if (author=="SMT01_WL04") {
123  double aa = 0.707;
124  double bb = 0.5;
125  double cc = 0.6;
126  double ni = deltacz/sigmaz;
127  double niI = sqrt(aa)*ni;
128  bias = 1.+1./deltacz*(pow(niI,2.)+bb*pow(niI,2.*(1.-cc))-pow(niI,2.*cc)/sqrt(aa)/(pow(niI,2.*cc)+bb*(1.-cc)*(1.-cc*0.5)));
129  }
130 
131  else if (author=="Tinker") { // Tinker et al. (2010)
132  double yy = log10(Delta);
133  double AA = 1.+0.24*yy*exp(-pow(4./yy,4));
134  double aa = 0.44*yy-0.88;
135  double BB = 0.183;
136  double bb = 1.5;
137  double CC = 0.019+0.107*yy+0.19*exp(-pow(4./yy,4));
138  double ccc = 2.4;
139  double ni = 1.686/sigmaz;
140  bias = 1.-AA*pow(ni,aa)/(pow(ni,aa)+pow(1.686,aa))+BB*pow(ni,bb)+CC*pow(ni,ccc);
141  }
142 
143  else
144  ErrorCBL("author = " + author + "!", "m_bias_halo_generator", "Bias.cpp");
145 
146  return bias;
147 }
148 
149 
150 // =====================================================================================
151 
152 
153 double cbl::cosmology::Cosmology::bias_eff (const double Mass_min, const double Mass_max, const double redshift, const std::string model_bias, const std::string model_MF, const std::string method_SS, const bool store_output, const std::string output_root, const double Delta, const double kk, 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)
154 {
155  // ---------- create/read the grid file with sigma(M) and its derivative ----------
156 
157  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
158 
159  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
160 
161  double Mass, Sigma, Dln_Sigma;
162  vector<double> mass, sigma, dlnsigma;
163 
164  while (fin >>Mass>>Sigma>>Dln_Sigma) {
165  if (Mass_min<Mass && Mass<Mass_max) {
166  mass.push_back(Mass);
167  sigma.push_back(Sigma);
168  dlnsigma.push_back(Dln_Sigma);
169  }
170  }
171 
172  if (mass.size()==0)
173  ErrorCBL("mass.size()=0!", "bias_eff", "Bias.cpp");
174 
175 
176  // ---------- compute the effective bias ----------
177 
178  double Bias_eff = 0., Norm = 0.;
179 
180  for (size_t mm=0; mm<mass.size()-1; mm++) {
181 
182  const double MF = mass_function(mass[mm], sigma[mm], dlnsigma[mm], redshift, model_MF, store_output, output_root, Delta, interpType, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
183 
184  Bias_eff += bias_halo(mass[mm], sigma[mm], redshift, model_bias, store_output, output_root, interpType, Delta, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file)*MF*(mass[mm+1]-mass[mm]);
185 
186  Norm += MF*(mass[mm+1]-mass[mm]);
187  }
188 
189  return Bias_eff/Norm;
190 }
191 
192 
193 // =====================================================================================
194 
195 
196 double cbl::cosmology::Cosmology::bias_eff (const std::vector<double> MM, const std::vector<double> MF, const double redshift, const std::string model_bias, const std::string method_SS, const bool store_output, const std::string output_root, const double Delta, const double kk, 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)
197 {
198  // ---------- create/read the grid file with sigma(M) and its derivative ----------
199 
200  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
201 
202  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
203 
204  double Mass, Sigma, Dln_Sigma;
205  vector<double> mass, sigma;
206 
207  while (fin >>Mass>>Sigma>>Dln_Sigma) {
208  if (Min(MM)<Mass && Mass<Max(MM)) {
209  mass.push_back(Mass);
210  sigma.push_back(Sigma);
211  }
212  }
213 
214  if (mass.size()==0)
215  ErrorCBL("mass.size()=0, Min(MM) = " + conv(Min(MM),par::fDP3) + ", Max(MM) = " + conv(Max(MM),par::fDP3) + ", file_grid = " + file_grid, "bias_eff", "Bias.cpp");
216 
217 
218  // ---------- compute the effective bias ----------
219 
220  double Bias_eff = 0., Norm = 0.;
221  double mf, sig, err = -1;
222 
223  for (size_t k=0; k<MM.size()-1; k++) {
224  mf = MF[k];
225  sig = interpolated(MM[k], mass, sigma, "Linear");
226 
227  if (err/sig>0.1)
228  ErrorCBL("err/sig = " + conv(err/sig, par::fDP3) + "!", "bias_eff", "Bias.cpp");
229 
230  Bias_eff += bias_halo(MM[k], sig, redshift, model_bias, store_output, output_root, interpType, Delta, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file)*mf*(MM[k+1]-MM[k]);
231  Norm += mf*(MM[k+1]-MM[k]);
232  }
233 
234  return Bias_eff/Norm;
235 }
236 
237 
238 // =====================================================================================
239 
240 
241 vector<double> cbl::cosmology::Cosmology::bias_eff_mass_grid (const std::vector<double> MM, const std::vector<double> redshift, const std::string model_bias, const std::string method_SS, const std::string meanType, const bool store_output, const std::string output_root, const double Delta_crit, const double kk, 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)
242 {
243  // ---------- create/read the grid file with sigma(M) and its derivative ----------
244 
245  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
246 
247  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
248 
249  double Mass, Sigma, Dln_Sigma;
250  vector<double> mass, sigma;
251 
252  while (fin >>Mass>>Sigma>>Dln_Sigma) {
253  if (Min(MM)<Mass && Mass<Max(MM)) {
254  mass.push_back(Mass);
255  sigma.push_back(Sigma);
256  }
257  }
258 
259  if (mass.size()==0)
260  ErrorCBL("mass.size()=0, Min(MM) = " + conv(Min(MM),par::fDP3) + ", Max(MM) = " + conv(Max(MM),par::fDP3) + ", file_grid = " + file_grid, "bias_eff_mass_grid", "Bias.cpp");
261 
262 
263  // ---------- compute the effective bias ----------
264 
265  if (meanType!="mean_bias" && meanType!="pair_mean_bias")
266  ErrorCBL("the chosen meanType is not allowed!", "bias_eff_mass_grid", "Bias.cpp");
267 
268  if (meanType=="mean_bias") {
269  vector<double> bias(MM.size());
270 
271  for (size_t k=0; k<MM.size(); k++) {
272  const double zz = (redshift.size()>1) ? redshift[k] : redshift[0];
273  bias[k] = bias_halo(MM[k], interpolated(MM[k], mass, sigma, "Linear"), zz, model_bias, store_output, output_root, interpType, Delta_crit/OmegaM(zz), kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
274  }
275 
276  return {Average(bias), cbl::Sigma(bias)/sqrt(MM.size())};
277  }
278 
279  else {
280  vector<double> bias2(MM.size());
281 
282  for (size_t k=0; k<MM.size(); ++k) {
283  const double z1 = (redshift.size()>1) ? redshift[k] : redshift[0];
284  for (size_t l=k+1; l<MM.size(); ++l) {
285  const double z2 = (redshift.size()>1) ? redshift[l] : redshift[0];
286  bias2[k] = bias_halo(MM[k], interpolated(MM[k], mass, sigma, "Linear"), z1, model_bias, store_output, output_root, interpType, Delta_crit/OmegaM(z1), kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file)*bias_halo(MM[l], interpolated(MM[l], mass, sigma, "Linear"), z2, model_bias, store_output, output_root, interpType, Delta_crit/OmegaM(z2), kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
287  }
288  }
289 
290  return {sqrt(Average(bias2)), sqrt(cbl::Sigma(bias2)/sqrt(MM.size()))};
291  }
292 }
293 
294 
295 // =====================================================================================
296 
297 
298 vector<double> cbl::cosmology::Cosmology::bias_eff_mass (const std::vector<double> MM, const std::vector<double> redshift, const std::string model_bias, const std::string method_SS, const std::string meanType, const bool store_output, const std::string output_root, const double Delta, const double kk, 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)
299 {
300  if (meanType!="mean_bias" && meanType!="pair_mean_bias")
301  ErrorCBL("the chosen meanType is not allowed!", "bias_eff_mass", "Bias.cpp");
302 
303  if (meanType=="mean_bias") {
304 
305  vector<double> bias(MM.size());
306 
307 #pragma omp parallel num_threads(omp_get_max_threads())
308  {
309 
310 #pragma omp for schedule(static, 2)
311  for (size_t k=0; k<MM.size(); k++) {
312  const double sigma = sqrt(sigma2M(MM[k], method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file, true));
313  bias[k] = bias_halo(MM[k], sigma, (redshift.size()>1) ? redshift[k] : redshift[0], model_bias, store_output, output_root, interpType, Delta, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
314  }
315 
316  }
317 
318  return {Average(bias), cbl::Sigma(bias)/sqrt(MM.size())};
319  }
320 
321  else {
322  vector<double> bias2(MM.size());
323 
324 #pragma omp parallel num_threads(omp_get_max_threads())
325  {
326 
327 #pragma omp for schedule(static, 2)
328  for (size_t k=0; k<MM.size(); ++k) {
329  const double z1 = (redshift.size()>1) ? redshift[k] : redshift[0];
330  const double sigma1 = sqrt(sigma2M(MM[k], method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file, true));
331  for (size_t l=k+1; l<MM.size(); ++l) {
332  const double z2 = (redshift.size()>1) ? redshift[l] : redshift[0];
333  const double sigma2 = sqrt(sigma2M(MM[l], method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file, true));
334  bias2[k] = bias_halo(MM[k], sigma1, z1, model_bias, store_output, output_root, interpType, Delta, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file)*bias_halo(MM[l], sigma2, z2, model_bias, store_output, output_root, interpType, Delta, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
335  }
336  }
337 
338  }
339 
340  return {sqrt(Average(bias2)), sqrt(cbl::Sigma(bias2)/sqrt(MM.size()))};
341  }
342 
343 }
344 
345 
346 // =====================================================================================
347 
348 
349 vector<double> cbl::cosmology::Cosmology::bias_eff_mass (const std::vector<double> mass, const std::vector<double> mass_grid, const std::vector<double> redshift, const std::string model_bias, const std::string method_SS, const std::string meanType, const bool store_output, const std::string output_root, const double Delta, const double kk, 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)
350 {
351  if (meanType!="mean_bias" && meanType!="pair_mean_bias")
352  ErrorCBL("the chosen meanType is not allowed!", "bias_eff_mass", "Bias.cpp");
353 
354  vector<double> Sigma;
355 
356  for (size_t k=0; k<mass_grid.size(); k++)
357  Sigma.emplace_back(sqrt(sigma2M(mass_grid[k], method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file, true)));
358  glob::FuncGrid sigma_interp(mass_grid, Sigma, "Spline");
359 
360  if (meanType=="mean_bias") {
361  vector<double> bias(mass.size());
362  for (size_t k=0; k<mass.size(); k++)
363  bias[k] = bias_halo(mass[k], sigma_interp(mass[k]), (redshift.size()>1) ? redshift[k] : redshift[0], model_bias, store_output, output_root, interpType, Delta, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
364  return {Average(bias), cbl::Sigma(bias)/sqrt(mass.size())};
365  }
366 
367  else {
368  vector<double> bias2(mass.size());
369  for (size_t k=0; k<mass.size(); k++) {
370  const double z1 = (redshift.size()>1) ? redshift[k] : redshift[0];
371  for (size_t l=k+1; l<mass.size(); ++l) {
372  const double z2 = (redshift.size()>1) ? redshift[l] : redshift[0];
373  bias2[k] = bias_halo(mass[k], sigma_interp(mass[k]), z1, model_bias, store_output, output_root, interpType, Delta, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file)*bias_halo(mass[l], sigma_interp(mass[l]), z2, model_bias, store_output, output_root, interpType, Delta, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
374  }
375  }
376 
377  return {sqrt(Average(bias2)), sqrt(cbl::Sigma(bias2)/sqrt(mass.size()))};
378  }
379 }
380 
381 
382 // =====================================================================================
383 
384 
385 vector<double> cbl::cosmology::Cosmology::bias_eff_selection_function (const glob::FuncGrid interp_sigma, const glob::FuncGrid interp_DlnSigma, const glob::FuncGrid interp_SF, const double Mass_min, const double Mass_max, const std::vector<double> redshift, const std::string model_bias, const std::string model_MF, const std::string method_SS, const double alpha, const bool store_output, const std::string output_root, const double Delta_crit, const double kk, 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)
386 {
387 
388  vector<double> Bias_eff(redshift.size(), 0.);
389 
390  vector<double> mass_vec = logarithmic_bin_vector(50, Mass_min, Mass_max);
391 
392  for (size_t i=0; i<redshift.size(); ++i) {
393 
394  const double DD = Delta_crit/OmegaM(redshift[i]);
395  vector<double> sigma = interp_sigma.eval_func(mass_vec);
396  vector<double> dlnsigma = interp_DlnSigma.eval_func(mass_vec);
397  vector<double> mass_func = mass_function(mass_vec, sigma, dlnsigma, redshift[i], model_MF, store_output, output_root, DD, interpType, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
398  vector<double> bias_func = bias_halo(mass_vec, sigma, redshift[i], model_bias, store_output, output_root, interpType, DD, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
399 
400  cbl::glob::FuncGrid interp_MF(mass_vec, mass_func, "Spline");
401  cbl::glob::FuncGrid interp_BH(mass_vec, bias_func, "Spline");
402 
403  auto integrand_num = [&] (const double lg_mass)
404  {
405  const double mass = exp(lg_mass);
406 
407  const double SF = interp_SF(mass/alpha);
408 
409  const double BH = interp_BH(mass);
410 
411  const double MF = interp_MF(mass);
412 
413  return SF*BH*MF*mass;
414  };
415 
416  auto integrand_denom = [&] (const double lg_mass)
417  {
418  const double mass = exp(lg_mass);
419 
420  const double SF = interp_SF(mass/alpha);
421 
422  const double MF = interp_MF(mass);
423 
424  return SF*MF*mass;
425  };
426 
427  Bias_eff[i] = wrapper::gsl::GSL_integrate_qag(integrand_num, log(Mass_min), log(Mass_max))/wrapper::gsl::GSL_integrate_qag(integrand_denom, log(Mass_min), log(Mass_max));
428  }
429 
430 
431  return Bias_eff;
432 }
433 
434 
435 // =====================================================================================
436 
437 
438 vector<double> cbl::cosmology::Cosmology::bias_eff_selection_function (const glob::FuncGrid interp_sigma, const glob::FuncGrid interp_DlnSigma, const glob::FuncGrid2D interp_SF, const double Mass_min, const double Mass_max, const std::vector<double> redshift, const std::string model_bias, const std::string model_MF, const std::string method_SS, const double alpha, const bool store_output, const std::string output_root, const double Delta_crit, const double kk, 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)
439 {
440  vector<double> Bias_eff(redshift.size(), 0.);
441 
442  for (size_t i=0; i<redshift.size(); ++i) {
443 
444  const double DD = Delta_crit/OmegaM(redshift[i]);
445 
446  auto integrand_num = [&] (const double mass)
447  {
448  const double sigma = interp_sigma(mass);
449  const double dlnsigma = interp_DlnSigma(mass);
450 
451  const double SF = interp_SF(mass/alpha, redshift[i]);
452 
453  const double BH = bias_halo(mass, sigma, redshift[i], model_bias, store_output, output_root, interpType, DD, kk, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
454 
455  const double MF = mass_function(mass, sigma, dlnsigma, redshift[i], model_MF, store_output, output_root, DD, interpType, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
456 
457  return SF*BH*MF;
458  };
459 
460  auto integrand_denom = [&] (const double mass)
461  {
462  const double sigma = interp_sigma(mass);
463  const double dlnsigma = interp_DlnSigma(mass);
464 
465  const double SF = interp_SF(mass/alpha, redshift[i]);
466 
467  const double MF = mass_function(mass, sigma, dlnsigma, redshift[i], model_MF, store_output, output_root, DD, interpType, norm, k_min, k_max, prec, method_SS, input_file, is_parameter_file);
468 
469  return SF*MF;
470  };
471 
472  Bias_eff[i] = wrapper::gsl::GSL_integrate_qag(integrand_num, Mass_min, Mass_max)/wrapper::gsl::GSL_integrate_qag(integrand_denom, Mass_min, Mass_max);
473  }
474 
475 
476  return Bias_eff;
477 }
478 
479 
480 // =====================================================================================
481 
482 
483 vector<double> cbl::cosmology::Cosmology::bias_eff_selection_function (const double Mass_min, const double Mass_max, const std::vector<double> redshift, const std::string model_bias, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector<int> column, const double alpha, const bool store_output, const std::string output_root, const double Delta_crit, const double kk, 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)
484 {
485  // ---------- create/read the grid file with sigmaM, dlnsigmaM ----------
486 
487  const string file_grid = create_grid_sigmaM(method_SS, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
488 
489  ifstream fin(file_grid.c_str()); checkIO(fin, file_grid);
490 
491  double Mass, Sigma, Dln_Sigma;
492  vector<double> mass, sigma, dlnsigma;
493 
494  while (fin >>Mass>>Sigma>>Dln_Sigma) {
495  if (Mass_min<Mass && Mass<Mass_max) {
496  mass.push_back(Mass);
497  sigma.push_back(Sigma);
498  dlnsigma.push_back(Dln_Sigma);
499  }
500  }
501 
502  if (mass.size()==0)
503  ErrorCBL("mass.size()=0, Mass_min = " + conv(Mass_min,par::fDP3) + ", Mass_max = " + conv(Mass_max,par::fDP3) + ", file_grid = " + file_grid, "bias_eff_selection_function", "Bias.cpp");
504 
505  const glob::FuncGrid interp_sigma(mass, sigma, "Spline");
506  const glob::FuncGrid interp_DlnSigma(mass, dlnsigma, "Spline");
507 
508 
509  // ---------- read the selection function ----------
510 
511  vector<double> mass_SF, redshift_SF;
512  vector<vector<double>> selection_function;
513 
514  read_matrix(selection_function_file, mass_SF, redshift_SF, selection_function, column);
515 
516  const glob::FuncGrid2D interp_SF(mass_SF, redshift_SF, selection_function, "Linear");
517 
518 
519  // ---------- compute the effective bias the given redshifts ----------
520 
521  return bias_eff_selection_function(interp_sigma, interp_DlnSigma, interp_SF, Mass_min, Mass_max, redshift, model_bias, model_MF, method_SS, alpha, store_output, output_root, Delta_crit, kk, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
522 
523 }
524 
525 
526 // =====================================================================================
527 
528 
529 void cbl::cosmology::Cosmology::generate_bias_eff_grid_one_cosmopar (std::vector<double> &parameter, std::vector<double> &bias_eff, const std::string dir_output, const std::string file_bias_eff_grid, const cbl::cosmology::CosmologicalParameter cosmoPar, const double min_par, const double max_par, const int nbin_par, const std::vector<double> mass, const std::vector<double> mass_grid, const std::vector<double> redshift, const std::string model_bias, const std::string method_SS, const std::string meanType, const bool store_output, const std::string output_root, const double Delta, const double kk, 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 cbl::cosmology::Cosmology cosmology_mass, const std::vector<double> redshift_source)
530 {
531  const double defaultValue = value(cosmoPar);
532 
533  const string file = dir_output+file_bias_eff_grid;
534 
535  ifstream fin(file.c_str());
536 
537  if (!fin) {
538 
539  vector<double> pp = linear_bin_vector(nbin_par, min_par, max_par);
540 
541  ofstream fout(file.c_str()); checkIO(fout, file);
542 
543  for (int i=0; i<nbin_par; i++) {
544  set_parameter(cosmoPar, pp[i]);
545 
546  (void)cosmology_mass;
547  (void)redshift_source;
548  /*
549  // convert the masses in the new cosmology, if they were computed in a different cosmology
550  vector<double> _mass(mass.size());
551  for (size_t mm=0; mm<mass.size(); mm++)
552  _mass[mm] = converted_mass(mass[mm], cosmology_mass, redshift[mm], (redshift_source.size()==redshift.size()) ? redshift_source[mm] : 0.);
553  */
554 
555  fout << pp[i] << " " << bias_eff_mass(/*_*/mass, mass_grid, redshift, model_bias, method_SS, meanType, store_output, output_root, Delta, kk, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)[0] << endl;
556 
557  }
558  fout.clear(); fout.close();
559  }
560 
561  fin.clear(); fin.close();
562 
563  fin.open(file.c_str()); checkIO(fin, file);
564 
565  parameter.erase(parameter.begin(), parameter.end());
566  bias_eff.erase(bias_eff.begin(), bias_eff.end());
567 
568  string line;
569  while (getline(fin, line)) {
570  stringstream SS(line); double _p, _b;
571  SS >> _p >> _b;
572  parameter.push_back(_p);
573  bias_eff.push_back(_b);
574  }
575  fin.clear(); fin.close();
576 
577  if (parameter.size()<2) ErrorCBL("parameter.size()<2; check the grid file: "+file+"!", "generate_bias_eff_grid_one_cosmopar", "Bias.cpp");
578 
579  set_parameter(cosmoPar, defaultValue);
580 }
581 
582 
583 // =====================================================================================
584 
585 
586 void cbl::cosmology::Cosmology::generate_bias_eff_grid_one_cosmopar (std::vector<double> &parameter, std::vector<double> &bias_eff, const std::string dir_output, const std::string file_bias_eff_grid, const cbl::cosmology::CosmologicalParameter cosmoPar, const double min_par, const double max_par, const int nbin_par, const double redshift, const double Mass_min, const double Mass_max, const std::string model_bias, const std::string model_MF, const std::string method_SS, const std::string selection_function_file, const std::vector<int> column, const double alpha, const bool store_output, const std::string output_root, const double Delta_crit, const double kk, 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)
587 {
588  const double defaultValue = value(cosmoPar);
589 
590  const string file = dir_output+file_bias_eff_grid;
591 
592  ifstream fin(file.c_str());
593 
594  if (!fin) {
595  const vector<double> pp = linear_bin_vector(nbin_par, min_par, max_par);
596 
597  ofstream fout(file.c_str());
598 
599  if ((int)pp.size()<nbin_par) ErrorCBL(conv(pp.size(), par::fINT)+" < nbin_par!", "generate_bias_eff_grid_one_cosmopar", "Bias.cpp");
600 
601  for (int i=0; i<nbin_par; i++) {
602  set_parameter(cosmoPar, pp[i]);
603 
604  fout << pp[i] << " " << bias_eff_selection_function(Mass_min, Mass_max, {redshift}, model_bias, model_MF, method_SS, selection_function_file, column, alpha, store_output, output_root, Delta_crit, kk, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)[0] << endl;
605 
606  }
607 
608  fout.clear(); fout.close();
609  }
610 
611  fin.clear(); fin.close();
612  fin.open(file.c_str());
613 
614  parameter.erase(parameter.begin(), parameter.end());
615  bias_eff.erase(bias_eff.begin(), bias_eff.end());
616 
617  string line;
618  while (getline(fin, line)) {
619  stringstream SS(line); double _p, _b;
620  SS >> _p >> _b;
621 
622  parameter.push_back(_p);
623  bias_eff.push_back(_b);
624  }
625  fin.clear(); fin.close();
626 
627  if (parameter.size()<2) ErrorCBL("parameter.size()<2; check the grid file: "+file+"!", "generate_bias_eff_grid_one_cosmopar", "Bias.cpp");
628 
629  set_parameter(cosmoPar, defaultValue);
630 }
631 
632 
633 // =====================================================================================
634 
635 
636 void cbl::cosmology::Cosmology::generate_bias_eff_grid_two_cosmopars (vector<double> &parameter1, vector<double> &parameter2, vector<vector<double>> &bias_eff, const std::string dir_output, const std::string file_bias_eff_grid, const cbl::cosmology::CosmologicalParameter cosmoPar1, const double min_par1, const double max_par1, const int nbin_par1, const cbl::cosmology::CosmologicalParameter cosmoPar2, const double min_par2, const double max_par2, const int nbin_par2, const std::vector<double> mass, const std::vector<double> mass_grid, const std::vector<double> redshift, const std::string model_bias, const std::string method_SS, const std::string meanType, const bool store_output, const std::string output_root, const double Delta, const double kk, 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 cbl::cosmology::Cosmology cosmology_mass, const std::vector<double> redshift_source)
637 {
638  double defaultValue1 = value(cosmoPar1);
639  double defaultValue2 = value(cosmoPar2);
640 
641  string file = dir_output+file_bias_eff_grid;
642 
643  ifstream fin(file.c_str());
644 
645  if (!fin) {
646  vector<double> pp1 = linear_bin_vector(nbin_par1, min_par1, max_par1);
647  vector<double> pp2 = linear_bin_vector(nbin_par2, min_par2, max_par2);
648 
649  ofstream fout(file.c_str()); checkIO(fout, file);
650 
651  for (int i=0; i<nbin_par1; i++) {
652  set_parameter(cosmoPar1, pp1[i]);
653  for (int j=0; j<nbin_par2; j++) {
654  set_parameter(cosmoPar2, pp2[j]);
655 
656  // convert the masses in the new cosmology, if they were computed in a different cosmology
657  vector<double> _mass(mass.size());
658  for (size_t mm=0; mm<mass.size(); mm++)
659  _mass[mm] = converted_mass(mass[mm], cosmology_mass, redshift[mm], (redshift_source.size()==redshift.size()) ? redshift_source[mm] : 0.);
660 
661  const double bias = bias_eff_mass(/*_*/mass, mass_grid, redshift, model_bias, method_SS, meanType, store_output, output_root, Delta, kk, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)[0];
662 
663  fout << pp1[i] << " " << pp2[j] << " " << bias << endl;
664  coutCBL << "parameter1 = " << pp1[i] << ", parameter2 = " << pp2[j] << ", bias = " << bias << endl;
665  }
666  fout << endl;
667  }
668  fout.clear(); fout.close();
669  }
670 
671  read_matrix(file, parameter1, parameter2, bias_eff);
672 
673  set_parameter(cosmoPar1, defaultValue1);
674  set_parameter(cosmoPar2, defaultValue2);
675 }
676 
677 
Cosmology.h
The class Cosmology.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
cbl::cosmology::Cosmology
The class Cosmology.
Definition: Cosmology.h:277
cbl::cosmology::Cosmology::generate_bias_eff_grid_one_cosmopar
void generate_bias_eff_grid_one_cosmopar(std::vector< double > &parameter, std::vector< double > &bias_eff, const std::string dir_output, const std::string file_bias_eff_grid, const cbl::cosmology::CosmologicalParameter cosmoPar, const double min_par, const double max_par, const int nbin_par, const std::vector< double > mass, const std::vector< double > mass_grid, const std::vector< double > redshift, const std::string model_bias, const std::string method_SS, const std::string meanType="mean_bias", const bool store_output=true, const std::string output_root="test", const double Delta_crit=200., const double kk=-1., 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 cbl::cosmology::Cosmology cosmology_mass={}, const std::vector< double > redshift_source={})
compute the effective bias of dark matter haloes, by averaging the bias of a set of haloes,...
Definition: Bias.cpp:529
cbl::cosmology::Cosmology::bias_eff_selection_function
std::vector< double > bias_eff_selection_function(const glob::FuncGrid interp_sigma, const glob::FuncGrid interp_DnSigma, const glob::FuncGrid interp_SF, const double Mass_min, const double Mass_max, const std::vector< double > redshift, const std::string model_bias, const std::string model_MF, const std::string method_SS, const double alpha=1., const bool store_output=true, const std::string output_root="test", const double Delta_crit=200., const double kk=-1., 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)
effective bias of dark matter haloes, computed using a given selection function; σ(mass) and dlnσ/dM ...
Definition: Bias.cpp:385
cbl::cosmology::Cosmology::generate_bias_eff_grid_two_cosmopars
void generate_bias_eff_grid_two_cosmopars(std::vector< double > &parameter1, std::vector< double > &parameter2, std::vector< std::vector< double >> &bias_eff, const std::string dir_output, const std::string file_bias_eff_grid, const cbl::cosmology::CosmologicalParameter cosmoPar1, const double min_par1, const double max_par1, const int nbin_par1, const cbl::cosmology::CosmologicalParameter cosmoPar2, const double min_par2, const double max_par2, const int nbin_par2, const std::vector< double > mass, const std::vector< double > mass_grid, const std::vector< double > redshift, const std::string model_bias, const std::string method_SS, const std::string meanType="mean_bias", const bool store_output=true, const std::string output_root="test", const double Delta=200., const double kk=-1., 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 cbl::cosmology::Cosmology cosmology_mass={}, const std::vector< double > redshift_source={})
effective bias of dark matter haloes, computed by averaging the bias of a set of haloes,...
Definition: Bias.cpp:636
cbl::cosmology::Cosmology::bias_eff_mass_grid
std::vector< double > bias_eff_mass_grid(const std::vector< double > MM, const std::vector< double > redshift, const std::string model_bias, const std::string method_SS, const std::string meanType="mean_bias", const bool store_output=true, const std::string output_root="test", const double Delta_crit=200., const double kk=-1., 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)
effective bias of dark matter haloes, computed by averaging the bias of a set of haloes,...
Definition: Bias.cpp:241
cbl::cosmology::Cosmology::m_bias_halo_generator
double m_bias_halo_generator(const double Sigma, const double redshift, const std::string author, const double Delta=200.) const
auxiliary function to compute the halo bias
Definition: Bias.cpp:89
cbl::cosmology::Cosmology::bias_halo
double bias_halo(const double Mass, const double redshift, const std::string author, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const std::string interpType="Linear", const double Delta=200., const double kk=-1., 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)
bias of dark matter haloes
Definition: Bias.cpp:44
cbl::cosmology::Cosmology::bias_eff_mass
std::vector< double > bias_eff_mass(const std::vector< double > MM, const std::vector< double > redshift, const std::string model_bias, const std::string method_SS, const std::string meanType="mean_bias", const bool store_output=true, const std::string output_root="test", const double Delta=200., const double kk=-1., 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)
effective bias of dark matter haloes, computed by averaging the bias of a set of haloes
Definition: Bias.cpp:298
cbl::cosmology::Cosmology::bias_eff
double bias_eff(const double Mass_min, const double Mass_max, const double redshift, const std::string model_bias, 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 double kk=-1., 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 effective bias of dark matter haloes, with masses in a given range and at a given mean redshift
Definition: Bias.cpp:153
cbl::glob::FuncGrid2D
The class FuncGrid2D.
Definition: FuncGrid.h:302
cbl::glob::FuncGrid
The class FuncGrid.
Definition: FuncGrid.h:55
cbl::glob::FuncGrid::eval_func
std::vector< double > eval_func(const std::vector< double > xx) const
evaluate the function at the xx points
Definition: FuncGrid.cpp:153
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::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::cosmology::CosmologicalParameter
CosmologicalParameter
the cosmological parameters
Definition: Cosmology.h:134
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::modelling::twopt::bias
double bias(const double Mmin, const double sigmalgM, const double M0, const double M1, const double alpha, const std::shared_ptr< void > inputs)
the mean galaxy bias
Definition: ModelFunction_TwoPointCorrelation1D_monopole.cpp:1026
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::Min
T Min(const std::vector< T > vect)
minimum element of a std::vector
Definition: Kernel.h:1324
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::Average
double Average(const std::vector< double > vect)
the average of a std::vector
Definition: Func.cpp:870
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::linear_bin_vector
std::vector< T > linear_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with linearly spaced values
Definition: Kernel.h:1604
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::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::Max
T Max(const std::vector< T > vect)
maximum element of a std::vector
Definition: Kernel.h:1336
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