CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
NG.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2010 by Federico Marulli and Cosimo Fedeli *
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 
36 #include "Cosmology.h"
37 
38 using namespace std;
39 
40 using namespace cbl;
41 using namespace cosmology;
42 
43 
44 // =====================================================================================
45 
46 
47 double cbl::cosmology::Cosmology::Am (const string method_Pk, const bool store_output, const string output_root, const int norm, const double k_min, const double k_max, const double prec, const string file_par)
48 {
49  double kk = 1.e-4;
50  bool NL = false;
51  double redshift = 0.;
52  return Pk_matter({kk}, method_Pk, NL, redshift, store_output, output_root, norm, k_min, k_max, prec, file_par)[0] / pow(kk, m_n_spec);
53 }
54 
55 
56 // =====================================================================================
57 
58 
59 double cbl::cosmology::Cosmology::potential_spectral_amplitude (const string method_Pk, const bool store_output, const string output_root, const int norm, const double k_min, const double k_max, const double prec, const string file_par)
60 {
61  return 2.78548e-14 * gsl_pow_2(m_Omega_matter) * Am(method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par);
62 }
63 
64 
65 // =====================================================================================
66 
67 
68 double cbl::cosmology::Cosmology::bispectrum (const vector<double> kk, const string method_Pk, const bool store_output, const string output_root, const int norm, const double k_min, const double k_max, const double prec, const string file_par)
69 {
70  double bs = 0.0;
71 
72  double bm = potential_spectral_amplitude(method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par);
73 
74  switch (m_type_NG)
75  {
76  case 1: // Local shape
77  bs = (pow(kk[0]*kk[1],m_n_spec-4.0)+pow(kk[0]*kk[2],m_n_spec-4.0)+pow(kk[1]*kk[2],m_n_spec-4.0));
78  bs *= 2.0*gsl_pow_2(bm);
79  break;
80 
81  case 2: // Equilateral shape
82  bs = pow(kk[0],(m_n_spec-4.0)/3.0)*pow(kk[1],2.0*(m_n_spec-4.0)/3.0)*pow(kk[2],m_n_spec-4.0)+pow(kk[2],(m_n_spec-4.0)/3.0)*pow(kk[0],2.0*(m_n_spec-4.0)/3.0)*pow(kk[1],m_n_spec-4.0)+
83  pow(kk[1],(m_n_spec-4.0)/3.0)*pow(kk[2],2.0*(m_n_spec-4.0)/3.0)*pow(kk[0],m_n_spec-4.0)+pow(kk[1],(m_n_spec-4.0)/3.0)*pow(kk[0],2.0*(m_n_spec-4.0)/3.0)*pow(kk[2],m_n_spec-4.0)+
84  pow(kk[2],(m_n_spec-4.0)/3.0)*pow(kk[1],2.0*(m_n_spec-4.0)/3.0)*pow(kk[0],m_n_spec-4.0)+pow(kk[0],(m_n_spec-4.0)/3.0)*pow(kk[2],2.0*(m_n_spec-4.0)/3.0)*pow(kk[1],m_n_spec-4.0)-
85  pow(kk[0]*kk[1],m_n_spec-4.0)-pow(kk[0]*kk[2],m_n_spec-4.0)-pow(kk[1]*kk[2],m_n_spec-4.0)-2.0*pow(kk[0]*kk[1]*kk[2],2.0*(m_n_spec-4.0)/3.0);
86  bs *= 6.0*gsl_pow_2(bm);
87  break;
88 
89  case 3: // Enfolded shape
90  bs = pow(kk[0]*kk[1],m_n_spec-4.0)+pow(kk[0]*kk[2],m_n_spec-4.0)+pow(kk[1]*kk[2],m_n_spec-4.0)+3.0*pow(kk[0]*kk[1]*kk[2],2.0*(m_n_spec-4.0)/3.0)-
91  pow(kk[0],(m_n_spec-4.0)/3.0)*pow(kk[1],2.0*(m_n_spec-4.0)/3.0)*pow(kk[2],m_n_spec-4.0)-pow(kk[2],(m_n_spec-4.0)/3.0)*pow(kk[0],2.0*(m_n_spec-4.0)/3.0)*pow(kk[1],m_n_spec-4.0)-
92  pow(kk[1],(m_n_spec-4.0)/3.0)*pow(kk[2],2.0*(m_n_spec-4.0)/3.0)*pow(kk[0],m_n_spec-4.0)-pow(kk[1],(m_n_spec-4.0)/3.0)*pow(kk[0],2.0*(m_n_spec-4.0)/3.0)*pow(kk[2],m_n_spec-4.0)-
93  pow(kk[2],(m_n_spec-4.0)/3.0)*pow(kk[1],2.0*(m_n_spec-4.0)/3.0)*pow(kk[0],m_n_spec-4.0)-pow(kk[0],(m_n_spec-4.0)/3.0)*pow(kk[2],2.0*(m_n_spec-4.0)/3.0)*pow(kk[1],m_n_spec-4.0);
94  bs *= 6.0*gsl_pow_2(bm);
95  break;
96 
97  case 4: // Orthogonal shape
98  bs = 3.0*pow(kk[0],(m_n_spec-4.0)/3.0)*pow(kk[1],2.0*(m_n_spec-4.0)/3.0)*pow(kk[2],m_n_spec-4.0)+3.0*pow(kk[2],(m_n_spec-4.0)/3.0)*pow(kk[0],2.0*(m_n_spec-4.0)/3.0)*pow(kk[1],m_n_spec-4.0)+
99  3.0*pow(kk[1],(m_n_spec-4.0)/3.0)*pow(kk[2],2.0*(m_n_spec-4.0)/3.0)*pow(kk[0],m_n_spec-4.0)+3.0*pow(kk[1],(m_n_spec-4.0)/3.0)*pow(kk[0],2.0*(m_n_spec-4.0)/3.0)*pow(kk[2],m_n_spec-4.0)+
100  3.0*pow(kk[2],(m_n_spec-4.0)/3.0)*pow(kk[1],2.0*(m_n_spec-4.0)/3.0)*pow(kk[0],m_n_spec-4.0)+3.0*pow(kk[0],(m_n_spec-4.0)/3.0)*pow(kk[2],2.0*(m_n_spec-4.0)/3.0)*pow(kk[1],m_n_spec-4.0)-
101  3.0*pow(kk[0]*kk[1],m_n_spec-4.0)-3.0*pow(kk[1]*kk[2],m_n_spec-4.0)-3.0*pow(kk[0]*kk[2],m_n_spec-4.0)-8.0*pow(kk[0]*kk[1]*kk[2],2.0*(m_n_spec-4)/3.0);
102  bs *= 6.0*gsl_pow_2(bm);
103  break;
104 
105  default:
106  ErrorCBL("", "bispectrum", "NG.cpp");
107  }
108 
109  return bs;
110 }
111 
112 
113 // =====================================================================================
114 
115 
116 double cbl::cosmology::Cosmology::mrk (const double kk, const double mass, const string method_Pk, const bool store_output, const string output_root, const int norm, const double k_min, const double k_max, const double prec, const string file_par)
117 {
118  double xx = kk * Radius(mass, m_RhoZero);
119 
120  double AA = Am(method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par);
121 
122  double TT = sqrt(Pk_matter({kk}, method_Pk, false, 0., store_output, output_root, norm, k_min, k_max, prec, file_par)[0] / AA / pow(kk,m_n_spec));
123 
124  return 5.99170e6 * gsl_pow_2(kk) * TopHat_WF(xx) * TT / m_Omega_matter;
125 }
126 
127 
128 // =====================================================================================
129 
131 
132 double cbl::cosmology::Cosmology::bias_kernel (double xx, void *params)
133 {
134  cbl::glob::GSL_f_pars *pp = (cbl::glob::GSL_f_pars *)params;
135 
136  int ni = 16;
137  double *XX = new double[ni];
138  double *Weight = new double[ni];
139  gauleg (0., 1., XX, Weight, ni);
140 
141  vector<double> km(3);
142  km[0] = xx;
143  km[2] = pp->kt;
144 
145  double mass = pp->mass;
146  string method_Pk = pp->method_Pk;
147  bool store_output = pp->store_output;
148  string output_root = pp->output_root;
149  int norm = pp->norm;
150  double k_min = pp->k_min;
151  double k_max = pp->k_max;
152  double prec = pp->prec;
153  string file_par = pp->file_par;
154 
155 
156  const double xi1 = -1.0, xi2 = 1.0;
157  double yi = 0.0;
158 
159  for (int i=0; i<ni; i++) {
160  double xi = xi1+(xi2-xi1)*XX[i];
161  km[1] = sqrt(gsl_pow_2(xx)+gsl_pow_2(pp->kt)+2.0*xx*pp->kt*xi);
162  yi += bispectrum(km, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par)*mrk(km[1], mass, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par)*Weight[i];
163  }
164 
165  yi *= (xi2-xi1)*gsl_pow_2(xx)*mrk(xx, mass, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par);
166 
167  return yi;
168 }
169 
170 
171 // =====================================================================================
172 
173 
174 double cbl::cosmology::Cosmology::frk_test (const double kk, const double mass, const string method_Pk, const bool store_output, const string output_root, const string interpType, const int norm, const double k_min, const double k_max, const double prec, const string input_file, const bool is_parameter_file)
175 {
176  cbl::glob::GSL_f_pars pp;
177  struct cbl::glob::GSL_f_pars *ppp = &pp;
178 
179  ppp->kt = kk;
180  ppp->mass = mass;
181  ppp->method_Pk = method_Pk;
182  ppp->output_root = output_root;
183  ppp->norm = norm;
184  ppp->k_min = k_min;
185  ppp->k_max = k_max;
186  ppp->prec = prec;
187  ppp->file_par = input_file; // check!
188  ppp->pt_Cosmology = this;
189 
190  gsl_function Func;
191 
192  Func.params = ppp;
193 
194  Func.function = &glob::GSL_bias_kernel_wrapper;
195 
196  double ibs = -1., err = -1.;
197 
198  gsl_integration_workspace *ww = gsl_integration_workspace_alloc(1000);
199  gsl_integration_qagiu (&Func, 0., 1.e-8, 1e-3, 1000, ww, &ibs, &err);
200 
201  gsl_integration_workspace_free(ww);
202 
203  double bm = potential_spectral_amplitude(method_Pk, store_output, output_root, norm, k_min, k_max, prec, input_file);
204 
205  double var = sigma2M(mass, method_Pk, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
206 
207  ibs /= 8.0*gsl_pow_2(par::pi)*var;
208  ibs /= bm*pow(kk, m_n_spec-4.0);
209 
210  return ibs;
211 }
212 
214 
215 // =====================================================================================
216 
218 
219 double cbl::glob::bias_kernel2 (const double xx, void *params)
220 {
221  struct cbl::glob::STR_NG *pp = (struct cbl::glob::STR_NG *) params;
222 
223  Cosmology cosm (pp->Omega_matter, pp->Omega_baryon, pp->Omega_neutrinos, pp->massless_neutrinos, pp->massive_neutrinos, pp->Omega_DE, pp->Omega_radiation, pp->hh, pp->scalar_amp, pp->scalar_pivot, pp->n_spec, pp->w0, pp->wa, pp->fNL, pp->type_NG, pp->tau, pp->model, pp->unit);
224 
225  int ni = 16;
226  double *XX = new double[ni];
227  double *Weight = new double[ni];
228  gauleg (0., 1., XX, Weight, ni);
229 
230  vector<double> km(3);
231  km[0] = xx;
232  km[2] = pp->kt;
233 
234  double mass = pp->mass;
235  string method_Pk = pp->method_Pk;
236  bool store_output = pp->store_output;
237  string output_root = pp->output_root;
238  int norm = pp->norm;
239  double k_min = pp->k_min;
240  double k_max = pp->k_max;
241  double prec = pp->prec;
242  string file_par = pp->file_par;
243 
244 
245  const double xi1 = -1.0, xi2 = 1.0;
246  double yi = 0.0;
247 
248  for (int i=0; i<ni; i++)
249  {
250  double xi = xi1+(xi2-xi1)*XX[i];
251  km[1] = sqrt(gsl_pow_2(xx)+gsl_pow_2(pp->kt)+2.0*xx*pp->kt*xi);
252  yi += cosm.bispectrum(km, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par)*cosm.mrk(km[1], mass, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par)*Weight[i];
253  }
254 
255  yi *= (xi2-xi1)*gsl_pow_2(xx)*cosm.mrk(xx, mass, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par);
256 
257  return yi;
258 }
259 
261 
262 
263 // =====================================================================================
264 
265 
266 double cbl::cosmology::Cosmology::frk (const double kk, const double mass, const string method_Pk, const bool store_output, const string output_root, const string interpType, const int norm, const double k_min, const double k_max, const double prec, const string input_file, const bool is_parameter_file)
267 {
268  cbl::Path path;
269  string dir_grid = path.DirCosmo()+"/Cosmology/Tables/grid_NG/bias_kernel/unit"+conv(m_unit,par::fINT)+"/";
270  string MK = "mkdir -p "+dir_grid; if (system (MK.c_str())) {};
271 
272  string Norm = (m_sigma8>0) ? "_sigma8"+conv(m_sigma8,par::fDP3) : "_scalar_amp"+conv(m_scalar_amp,par::ee3);
273  string file_grid = dir_grid+"grid"+Norm+"_h"+conv(m_hh,par::fDP3)+"_OmB"+conv(m_Omega_baryon,par::fDP3)+"_OmCDM"+conv(m_Omega_CDM,par::fDP3)+"_OmL"+conv(m_Omega_DE,par::fDP3)+"_OmN"+conv(m_Omega_neutrinos,par::fDP3)+"_typeNG"+conv(m_type_NG,par::fINT)+"_k"+conv(kk,par::fDP5)+".dat";
274 
275  int bin = 100;
276  double x_min = 1.e-3;
277  double x_max = 1.e3;
278  vector<double> xx, yy;
279 
280  cbl::glob::STR_NG str;
281  str.Omega_matter = m_Omega_matter;
282  str.Omega_baryon = m_Omega_baryon;
283  str.Omega_neutrinos = m_Omega_neutrinos;
284  str.massless_neutrinos = m_massless_neutrinos;
285  str.massive_neutrinos = m_massive_neutrinos;
286  str.Omega_DE = m_Omega_DE;
287  str.Omega_radiation = m_Omega_radiation;
288  str.hh = m_hh;
289  str.scalar_amp = m_scalar_amp;
290  str.scalar_pivot = m_scalar_pivot;
291  str.n_spec = m_n_spec;
292  str.w0 = m_w0;
293  str.wa = m_wa;
294  str.fNL = m_fNL;
295  str.type_NG = m_type_NG;
296  str.tau = m_tau;
297  str.model = m_model;
298  str.unit = m_unit;
299  str.kt = kk;
300  str.mass = mass;
301  str.method_Pk = method_Pk;
302  str.store_output = store_output;
303  str.output_root = output_root;
304  str.norm = norm;
305  str.k_min = k_min;
306  str.k_max = k_max;
307  str.prec = prec;
308  str.file_par = input_file; // check!
309 
310  bin_function(file_grid, glob::bias_kernel2, &str, bin, x_min, x_max, "loglin", xx, yy);
311 
312 
313  cbl::glob::STR_grid str_grid;
314  str_grid._xx = xx;
315  str_grid._yy = yy;
316 
317  gsl_function Func;
318 
319  Func.function = func_grid_loglin;
320  Func.params = &str_grid;
321 
322  double ibs = -1., err = -1.;
323 
324  gsl_integration_workspace *ww = gsl_integration_workspace_alloc(1000);
325  //gsl_integration_qagiu (&Func, 0., 1.e-8, 1e-3, 1000, ww, &ibs, &err);
326  gsl_integration_qag(&Func, 1.e-3, 1.e3, 0., 1.e-3, 1000, 6, ww, &ibs, &err);
327 
328  gsl_integration_workspace_free(ww);
329 
330  double bm = potential_spectral_amplitude(method_Pk, store_output, output_root, norm, k_min, k_max, prec, input_file);
331 
332  double var = sigma2M(mass, method_Pk, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
333 
334  ibs /= 8.0*gsl_pow_2(par::pi)*var;
335  ibs /= bm*pow(kk,m_n_spec-4.0);
336 
337  return ibs;
338 }
339 
340 
341 // =====================================================================================
342 
344 
345 double cbl::glob::GSL_bias_kernel_wrapper (const double xx, void *params)
346 {
347  cbl::glob::GSL_f_pars *pp = (cbl::glob::GSL_f_pars *)params;
348  return pp->pt_Cosmology->bias_kernel(xx,params);
349 }
350 
352 
353 // =====================================================================================
354 
355 
356 double cbl::cosmology::Cosmology::bias_correction (const double kk, const double mass, const string method_Pk, const bool store_output, const string output_root, const string interpType, const int norm, const double k_min, const double k_max, const double prec, const string input_file, const bool is_parameter_file)
357 {
358  return m_fNL * 0.8 * frk(kk, mass, method_Pk, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file) / mrk(kk, mass, method_Pk, store_output, output_root, norm, k_min, k_max, prec, input_file);
359 }
360 
361 
362 // =====================================================================================
363 
365 
366 double cbl::glob::skewness_kernel (double *kk, size_t dim, void *params)
367 {
368  (void)dim;
369 
370  struct cbl::glob::STR_NG *pp = (struct cbl::glob::STR_NG *) params;
371 
372  Cosmology cosm (pp->Omega_matter, pp->Omega_baryon, pp->Omega_neutrinos, pp->massless_neutrinos, pp->massive_neutrinos, pp->Omega_DE, pp->Omega_radiation, pp->hh, pp->scalar_amp, pp->scalar_pivot, pp->n_spec, pp->w0, pp->wa, pp->fNL, pp->type_NG, pp->tau, pp->output_root, pp->unit);
373 
374  int ni = 16;
375  double *XX = new double[ni];
376  double *Weight = new double[ni];
377  gauleg (0., 1., XX, Weight, ni);
378 
379  vector<double> km(3);
380  km[0] = kk[0];
381  km[2] = kk[1];
382 
383  double mass = pp->mass;
384  string method_Pk = pp->method_Pk;
385  bool store_output = pp->store_output;
386  string output_root = pp->output_root;
387  int norm = pp->norm;
388  double k_min = pp->k_min;
389  double k_max = pp->k_max;
390  double prec = pp->prec;
391  string file_par = pp->file_par;
392 
393 
394  const double xi1 = -1.0, xi2 = 1.0;
395  double xi, yi = 0.0;
396 
397  for (int i=0; i<ni; i++) {
398  xi = xi1+(xi2-xi1)*XX[i];
399  km[1] = sqrt(gsl_pow_2(kk[0])+gsl_pow_2(kk[1])+2.0*kk[0]*kk[1]*xi);
400  yi += cosm.bispectrum(km, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par)*cosm.mrk(km[1], mass, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par)*Weight[i];
401  }
402 
403  yi *= (xi2-xi1);
404 
405  return cosm.mrk(kk[0], mass, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par)*cosm.mrk(kk[1], mass, method_Pk, store_output, output_root, norm, k_min, k_max, prec, file_par)*yi*gsl_pow_2(kk[0]*kk[1]);
406 }
407 
409 
410 // =====================================================================================
411 
412 
413 
414 double cbl::cosmology::Cosmology::skewness (const double mass, const string method_Pk, const bool store_output, const string output_root, const string interpType, const int norm, const double k_min, const double k_max, const double prec, const string input_file, const bool is_parameter_file)
415 {
416  cbl::Path path;
417  string dir_grid = path.DirCosmo()+"/Cosmology/Tables/grid_NG/skewness_kernel/unit"+conv(m_unit,par::fINT)+"/";
418  string MK = "mkdir -p "+dir_grid; if (system (MK.c_str())) {};
419 
420  string Norm = (m_sigma8>0) ? "_sigma8"+conv(m_sigma8,par::fDP3) : "_scalar_amp"+conv(m_scalar_amp,par::ee3);
421  string file_grid = dir_grid+"grid"+Norm+"_h"+conv(m_hh,par::fDP3)+"_OmB"+conv(m_Omega_baryon,par::fDP3)+"_OmCDM"+conv(m_Omega_CDM,par::fDP3)+"_OmL"+conv(m_Omega_DE,par::fDP3)+"_OmN"+conv(m_Omega_neutrinos,par::fDP3)+"_typeNG"+conv(m_type_NG,par::fINT)+"_lgMass"+conv(log10(mass),par::fDP2)+".dat";
422 
423  // check !!!
424  int bin = 20;
425  double x1_min = 1.e-5;
426  double x1_max = 1.e2;
427  double x2_min = 1.e-5;
428  double x2_max = 1.e2;
429  vector<double> xx1, xx2;
430  vector< vector<double> > yy;
431 
432  cbl::glob::STR_NG str;
433  str.Omega_matter = m_Omega_matter;
434  str.Omega_baryon = m_Omega_baryon;
435  str.Omega_neutrinos = m_Omega_neutrinos;
436  str.massless_neutrinos = m_massless_neutrinos;
437  str.massive_neutrinos = m_massive_neutrinos;
438  str.Omega_DE = m_Omega_DE;
439  str.Omega_radiation = m_Omega_radiation;
440  str.hh = m_hh;
441  str.scalar_amp = m_scalar_amp;
442  str.scalar_pivot = m_scalar_pivot;
443  str.n_spec = m_n_spec;
444  str.w0 = m_w0;
445  str.wa = m_wa;
446  str.fNL = m_fNL;
447  str.type_NG = m_type_NG;
448  str.model = m_model;
449  str.unit = m_unit;
450  str.kt = -1.; // check!!!
451  str.mass = mass;
452  str.method_Pk = method_Pk;
453  str.store_output = store_output;
454  str.output_root = output_root;
455  str.norm = norm;
456  str.k_min = k_min;
457  str.k_max = k_max;
458  str.prec = prec;
459  str.file_par = input_file; // check!!!
460 
461  bin_function_2D(file_grid, glob::skewness_kernel, &str, bin, x1_min, x1_max, x2_min, x2_max, "loglin", xx1, xx2, yy);
462 
463  cbl::glob::STR_grid_2D str_grid_2D;
464  str_grid_2D._xx1 = xx1;
465  str_grid_2D._xx2 = xx2;
466  str_grid_2D._yy = yy;
467 
468 
469  int dim = 2;
470  double *kl = new double[dim];
471  double *kh = new double[dim];
472 
473  for (int i=0; i<dim; i++) {
474  kl[i] = 0.0;
475  kh[i] = 30.0/pow(Radius(mass,m_RhoZero),0.7);
476  }
477 
478 
479  gsl_monte_function Func;
480  Func.f = func_grid_loglin_2D;
481  Func.dim = dim;
482  Func.params = &str_grid_2D;
483 
484 
485  double ibs = -1., err = -1.;
486 
487  const size_t cl = 75000;
488 
489  // random number generator
490  const gsl_rng_type *TT = gsl_rng_default;
491  gsl_rng *rn = gsl_rng_alloc (TT);
492 
493  coutCBL <<"Please wait, I'm computing the Monte Carlo integral..."<<endl;
494 
495  gsl_monte_vegas_state *st = gsl_monte_vegas_alloc (dim);
496  gsl_monte_vegas_integrate(&Func, kl, kh, dim, cl/100, rn, st, &ibs, &err);
497 
498  gsl_monte_vegas_free(st);
499 
500  double var = sigma2M(mass, method_Pk, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
501 
502  return ibs/(8.0*gsl_pow_4(par::pi))/gsl_pow_2(var)*m_fNL;
503 }
504 
505 
506 // =====================================================================================
507 
508 
509 double cbl::cosmology::Cosmology::dskewnessdM (const double mass, const string method_Pk, const bool store_output, const string output_root, const string interpType, const int norm, const double k_min, const double k_max, const double prec, const string input_file, const bool is_parameter_file)
510 {
511  double dlogm = 0.1;
512  double mInf = 6.0, mSup = 16.0;
513 
514  if (log10(mass) < mInf) ErrorCBL("mass should be > 10^6 Msun/h", "dskewnessdM", "NG.cpp");
515 
516  if (log10(mass) > mSup) ErrorCBL("mass should be < 10^16 Msun/h", "dskewnessdM", "NG.cpp");
517 
518 
519  double logm1 = log10(mass)-dlogm;
520  double logm2 = log10(mass)+dlogm;
521 
522  if (logm1 < mInf) logm1 = mInf;
523  if (logm2 > mSup) logm2 = mSup;
524 
525  double sk = skewness(mass, method_Pk, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
526  double logsk1 = 0.0;
527  double logsk2 = 0.0;
528 
529  double M1 = pow(10.0,logm1);
530  double M2 = pow(10.0,logm2);
531 
532  if (m_fNL != 0.0) {
533  switch (m_type_NG) {
534  default:
535  logsk1 = log10(skewness(M1, method_Pk, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)/m_fNL);
536  logsk2 = log10(skewness(M2, method_Pk, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)/m_fNL);
537  break;
538  case 4:
539  logsk1 = log10(-skewness(M1, method_Pk, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)/m_fNL);
540  logsk2 = log10(-skewness(M2, method_Pk, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file)/m_fNL);
541  break;
542  }
543  }
544 
545  return sk*(logsk2-logsk1)/(logm2-logm1)/mass;
546 }
547 
548 
549 // =====================================================================================
550 
551 
552 double cbl::cosmology::Cosmology::MF_correction (const double mass, const double redshift, const string method_Pk, const bool store_output, const string output_root, const string interpType, const int norm, const double k_min, const double k_max, const double prec, const string input_file, const bool is_parameter_file)
553 {
554  double dc = deltac(redshift)*sqrt(0.8);
555  double gf = DN(redshift); // check the normalization of D(z)!!!
556 
557  double SSS = sigma2M(mass, method_Pk, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file);
558  double sm = sqrt(SSS);
559 
560  double DlnSigmaDlnM = dnsigma2M(1, mass, method_Pk, 0., store_output, output_root, interpType, k_max, input_file, is_parameter_file)*(mass/(2.*SSS));
561 
562  double sk = skewness(mass, method_Pk, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
563 
564  double dsk = dskewnessdM(mass, method_Pk, store_output, output_root, interpType, norm, k_min, k_max, prec, input_file, is_parameter_file);
565 
566  double mfc1 = DlnSigmaDlnM*(dc/gf/sm+sk*sm/6.0*(gsl_pow_4(dc/gf/sm)-2.0*gsl_pow_2(dc/gf/sm)-1.0))/mass;
567  double mfc2 = dsk*sm/6.0*(gsl_pow_2(dc/gf/sm)-1.0);
568  double mfc = mfc1+mfc2;
569  double psg = 1.;
570 
571  if (mfc > 0.0) mfc = 0.0;
572  else psg = DlnSigmaDlnM*(dc/gf/sm)/mass;
573 
574  return mfc/psg;
575 }
Cosmology.h
The class Cosmology.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
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::MF_correction
double MF_correction(const double mass, const double redshift, const std::string method_Pk, const bool store_output=true, const std::string output_root="test", 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)
correction to the halo mass in non-Gaussian cosmologies
Definition: NG.cpp:552
cbl::cosmology::Cosmology::dskewnessdM
double dskewnessdM(const double mass, const std::string method_Pk, const bool store_output=true, const std::string output_root="test", 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 derivative of the skewness, ds/dM
Definition: NG.cpp:509
cbl::cosmology::Cosmology::bispectrum
double bispectrum(const std::vector< double > kk, const std::string method_Pk, const bool store_output=true, const std::string output_root="test", const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string file_par=par::defaultString)
the bispectrum
Definition: NG.cpp:68
cbl::cosmology::Cosmology::potential_spectral_amplitude
double potential_spectral_amplitude(const std::string method_Pk, const bool store_output=true, const std::string output_root="test", const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string file_par=par::defaultString)
the potential spectral amplitude
Definition: NG.cpp:59
cbl::cosmology::Cosmology::bias_correction
double bias_correction(const double kk, const double mass, const std::string method_Pk, const bool store_output=true, const std::string output_root="test", 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)
correction to the halo bias in non-Gaussian cosmologies
Definition: NG.cpp:356
cbl::cosmology::Cosmology::frk
double frk(const double kk, const double mass, const std::string method_Pk, const bool store_output=true, const std::string output_root="test", 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)
auxiliary function to estimate cosmological quantities in non-Gaussian cosmologies
Definition: NG.cpp:266
cbl::cosmology::Cosmology::skewness
double skewness(const double mass, const std::string method_Pk, const bool store_output=true, const std::string output_root="test", 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 skewness
Definition: NG.cpp:414
cbl::cosmology::Cosmology::mrk
double mrk(const double kk, const double mass, const std::string method_Pk, const bool store_output=true, const std::string output_root="test", const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string file_par=par::defaultString)
auxiliary function to estimate cosmological quantities in non-Gaussian cosmologies
Definition: NG.cpp:116
cbl::cosmology::Cosmology::Am
double Am(const std::string method_Pk, const bool store_output=true, const std::string output_root="test", const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string file_par=par::defaultString)
the amplitude of the matter power spectrum
Definition: NG.cpp:47
cbl::par::fDP2
static const char fDP2[]
conversion symbol for: double -> std::string
Definition: Constants.h:133
cbl::par::ee3
static const char ee3[]
conversion symbol for: double -> std::string
Definition: Constants.h:160
cbl::par::fDP3
static const char fDP3[]
conversion symbol for: double -> std::string
Definition: Constants.h:136
cbl::par::fDP5
static const char fDP5[]
conversion symbol for: double -> std::string
Definition: Constants.h:142
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
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::Radius
T Radius(const T Mass, const T Rho)
the radius of a sphere of a given mass and density
Definition: Func.h:1220
cbl::TopHat_WF
T TopHat_WF(const T kR)
the top-hat window function
Definition: Func.h:1195
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::bin_function_2D
void bin_function_2D(const std::string file_grid, double func(double *, size_t, void *), void *par, const int bin, const double x1_min, const double x1_max, const double x2_min, const double x2_max, const std::string binning, std::vector< double > &xx1, std::vector< double > &xx2, std::vector< std::vector< double >> &yy)
create a 2D grid given an input function
Definition: Func.cpp:1183
cbl::bin_function
void bin_function(const std::string file_grid, double func(double, void *), void *par, const int bin, const double x_min, const double x_max, const std::string binning, std::vector< double > &xx, std::vector< double > &yy)
create a 1D grid given an input function
Definition: Func.cpp:1118