CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
PkXi.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 
36 #include "FuncGrid_Bspline.h"
37 #include "Cosmology.h"
38 #include <regex>
39 
40 using namespace std;
41 
42 using namespace cbl;
43 using namespace cosmology;
44 
45 
46 // =====================================================================================
47 
48 
49 double cbl::cosmology::Cosmology::As (const double sigma8) const //check
50 {
51  return pow(sigma8/1.79e4*pow(m_Omega_baryon*m_hh*m_hh/0.024, 1./3.)*pow(m_Omega_matter*m_hh*m_hh/0.14, -0.563)*pow(7.808*m_hh, 0.5*(1.-m_n_spec))*pow(m_hh/0.72, -0.693)*0.76/gg(0.), 2);
52 }
53 
54 
55 // =====================================================================================
56 
57 
58 double cbl::cosmology::Cosmology::sigma8_interpolated (const double redshift) const
59 {
60  const double wm = m_Omega_matter*m_hh*m_hh;
61  const double wb = m_Omega_baryon*m_hh*m_hh;
62  const double sigma8 = 0.1058*pow(m_scalar_amp/2.196e-9, 0.5)*pow(wm/0.1426, 0.520)*pow(wb/0.02205, -0.294)*pow(m_hh/0.673, 0.683)*pow((m_massless_neutrinos+m_massive_neutrinos)/3.046, -0.24)*exp(0.3727*(m_n_spec-0.96))*pow(1-m_Omega_k, 0.175)*DN(redshift, 9.);
63 
64  return ((m_Omega_neutrinos>0) ? 0.995*sigma8 : sigma8);
65 }
66 
67 
68 // =====================================================================================
69 
70 
71 std::string cbl::cosmology::Cosmology::Pk_output_file (const string code, const bool NL, const double redshift, const bool run, const bool store_output, const string output_root, const double k_max, const string file_par)
72 {
73  cbl::Path path;
74  string dir_loc = path.fullpath(path.DirLoc());
75  string dir_cosmo = path.DirCosmo();
76 
77  string dir_grid;
78  if (NL==0) dir_grid = "output_linear/";
79  else if (NL==1) dir_grid = "output_nonlinear/";
80  else ErrorCBL("", "Pk_output_file", "PkXi.cpp");
81 
82  dir_grid += "h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6)+"/";
83 
84  //dir_grid += (output_root=="test") ? "/" : "_"+output_root;
85 
86  string dir = dir_cosmo+"/External/"+code+"/";
87  string dirC = dir_cosmo+"/External/CAMB/";
88  if (chdir(dirC.c_str())) {}
89 
90  string dir_output = dir+dir_grid+"pk.dat";
91 
92  if (run) {
93  vector<double> lgkk, lgPk;
94  Table_PkCodes(code, NL, lgkk, lgPk, redshift, store_output, output_root, k_max, file_par);
95  }
96 
97  if (chdir(dir_loc.c_str())) {}
98 
99  return dir_output;
100 }
101 
102 
103 // =====================================================================================
104 
105 
106 void cbl::cosmology::Cosmology::run_CAMB (const bool NL, const double redshift, const std::string output_root, const std::string output_dir, const double k_max, const std::string file_par) const
107 {
108  cbl::Path path;
109  string dir_CAMB = path.DirCosmo()+"/External/CAMB/fortran/";
110 
111  string File_par = file_par;
112 
113  bool delete_output = (output_dir==par::defaultString) ? true : false;
114 
115  if (output_dir!=par::defaultString) if (system(("mkdir -p "+output_dir).c_str())) {}
116 
117  string OutputRoot = (output_dir==par::defaultString) ? dir_CAMB+"../inifiles/"+output_root : path.DirLoc()+output_dir+"/"+output_root;
118 
119  OutputRoot = (omp_get_max_threads()>1) ? OutputRoot+"_t"+conv(omp_get_thread_num(), par::fINT) : OutputRoot;
120 
121  if (File_par==par::defaultString) {
122 
123  // --------------------------------------------------------------------------
124  // --------- set the cosmological parameters in the file params.ini ---------
125  // --------------------------------------------------------------------------
126 
127  string file_par_default = dir_CAMB+"../inifiles/params_cut.ini";
128  File_par = OutputRoot+"_params.ini";
129 
130  if (system(("cp "+file_par_default+" "+File_par).c_str())) {}
131  ofstream fout(File_par.c_str(), std::ios_base::app | std::ios_base::out);
132  double HH0 = m_hh*100.;
133 
134  fout << "output_root = " << OutputRoot << endl;
135  fout << "do_nonlinear = " << conv(NL, par::fINT) << endl;
136  fout << "hubble = " << conv(HH0, par::fDP6) << endl;
137  fout << "ombh2 =" << conv(m_Omega_baryon*m_hh*m_hh, par::fDP6) << endl;
138  fout << "omch2 = " << conv(m_Omega_CDM*m_hh*m_hh, par::fDP6) << endl;
139  fout << "omk = " << conv(m_Omega_k, par::fDP6) << endl;
140  fout << "omnuh2 = " << conv(m_Omega_neutrinos*m_hh*m_hh, par::fDP6) << endl;
141  fout << "transfer_redshift(1) = " << conv(redshift, par::fDP6) << endl;
142  fout << "massless_neutrinos = " << conv(m_massless_neutrinos, par::fDP6) << endl;
143  fout << "massive_neutrinos = " << conv(m_massive_neutrinos, par::fINT) << endl;
144  fout << "scalar_spectral_index(1) = " << conv(m_n_spec, par::fDP6) << endl;
145  fout << "w = " << conv(m_w0, par::fDP6) << endl;
146  fout << "wa = " << conv(m_wa, par::fDP6) << endl;
147  if (m_scalar_amp>0) {
148  fout << "scalar_amp(1) = " << conv(m_scalar_amp, par::ee3) << endl;
149  fout << "pivot_scalar = " << conv(m_scalar_pivot, par::fDP6) << endl;
150  }
151  fout << "transfer_kmax = "+conv(k_max, par::fDP6) << endl;
152  fout << "re_optical_depth = "+conv(m_tau, par::fDP6) << endl;
153  fout << "feedback_level = 1" << endl;
154  fout << "print_sigma8 = T" << endl;
155  fout << endl;
156 
157  fout.clear(); fout.close();
158 
159  }
160 
161  // --------------------------------------------------------------------------
162 
163  if (system((dir_CAMB+"camb "+File_par).c_str())) {}
164 
165  if (delete_output) {
166  string RM = "rm -f "+OutputRoot+"*";
167  if (system(RM.c_str())) {}
168  }
169 
170 }
171 
172 
173 // =====================================================================================
174 
175 
176 void cbl::cosmology::Cosmology::run_CAMB (std::vector<double> &kk, std::vector<double> &Pk, const bool NL, const double redshift, const std::string output_root, const std::string output_dir, const double k_max, const std::string file_par) const
177 {
178  cbl::Path path;
179  string dir_CAMB = path.DirCosmo()+"/External/CAMB/fortran/";
180  string File_par = file_par;
181  bool delete_output = (output_dir==par::defaultString) ? true : false;
182 
183  if (output_dir!=par::defaultString) if (system(("mkdir -p "+output_dir).c_str())) {}
184 
185  string OutputRoot = (output_dir==par::defaultString) ? dir_CAMB+"../inifiles/"+output_root : path.DirLoc()+output_dir+"/"+output_root;
186 
187  OutputRoot = (omp_get_max_threads() > 1) ? OutputRoot+"_t"+conv(omp_get_thread_num(), par::fINT) : OutputRoot;
188 
189  if (File_par==par::defaultString) {
190 
191  // --------------------------------------------------------------------------
192  // --------- set the cosmological parameters in the file params.ini ---------
193  // --------------------------------------------------------------------------
194 
195  string file_par_default = dir_CAMB+"../inifiles/params_cut.ini";
196  File_par = OutputRoot+"_params.ini";
197 
198  if (system(("cp "+file_par_default+" "+File_par).c_str())) {}
199  ofstream fout(File_par.c_str(), std::ios_base::app | std::ios_base::out);
200  double HH0 = m_hh*100.;
201 
202  fout << "output_root = " << OutputRoot << endl;
203  fout << "do_nonlinear = " << conv(NL, par::fINT) << endl;
204  fout << "hubble = " << conv(HH0, par::fDP6) << endl;
205  fout << "ombh2 =" << conv(m_Omega_baryon*m_hh*m_hh, par::fDP6) << endl;
206  fout << "omch2 = " << conv(m_Omega_CDM*m_hh*m_hh, par::fDP6) << endl;
207  fout << "omk = " << conv(m_Omega_k, par::fDP6) << endl;
208  fout << "omnuh2 = " << conv(m_Omega_neutrinos*m_hh*m_hh, par::fDP6) << endl;
209  fout << "transfer_redshift(1) = " << conv(redshift, par::fDP6) << endl;
210  fout << "massless_neutrinos = " << conv(m_massless_neutrinos, par::fDP6) << endl;
211  fout << "massive_neutrinos = " << conv(m_massive_neutrinos, par::fINT) << endl;
212  fout << "scalar_spectral_index(1) = " << conv(m_n_spec, par::fDP6) << endl;
213  fout << "w = " << conv(m_w0, par::fDP6) << endl;
214  fout << "wa = " << conv(m_wa, par::fDP6) << endl;
215  if (m_scalar_amp>0) {
216  fout << "scalar_amp(1) = " << conv(m_scalar_amp, par::ee3) << endl;
217  fout << "pivot_scalar = " << conv(m_scalar_pivot, par::fDP6) << endl;
218  }
219  fout << "transfer_kmax = "+conv(k_max, par::fDP6) << endl;
220  fout << "re_optical_depth = "+conv(m_tau, par::fDP6) << endl;
221  fout << "feedback_level = -1" << endl;
222  fout << "print_sigma8 = F" << endl;
223  fout << endl;
224 
225  fout.clear(); fout.close();
226  }
227 
228 
229  // --------------------------------------------------------------------------
230 
231  if (system((dir_CAMB+"camb "+File_par).c_str())) {}
232 
233  // read the power spectrum
234  string file_inCAMB = OutputRoot+"_matterpower.dat";
235  ifstream fin(file_inCAMB.c_str()); checkIO(fin, file_inCAMB);
236 
237  kk.erase(kk.begin(), kk.end());
238  Pk.erase(Pk.begin(), Pk.end());
239 
240  string line;
241  getline(fin, line);
242 
243  double KK, PK;
244  while (fin >> KK >> PK)
245  if (KK>0 && PK>0) {
246  kk.push_back(KK);
247  Pk.push_back(PK);
248  }
249  fin.clear(); fin.close();
250 
251  if (delete_output) if (system(("rm -f "+OutputRoot+"*").c_str())) {}
252 }
253 
254 
255 // =====================================================================================
256 
257 
258 void cbl::cosmology::Cosmology::Table_PkCodes (const std::string code, const bool NL, std::vector<double> &lgkk, std::vector<double> &lgPk, const double redshift, const bool store_output, const std::string output_root, const double k_max, const std::string file_par) const
259 {
260  if (code=="MPTbreeze-v1") {
261  if (m_sigma8<0)
262  ErrorCBL("sigma8<0! The function set_sigma8() can be used to set the value of sigma8!", "Table_PkCodes", "PkXi.cpp");
263  if (NL)
264  WarningMsgCBL("NL is ignored by MPTbreeze-v1, that provides in output the non-linear power spectrum", "Table_PkCodes", "PkXi.cpp");
265  }
266 
267  if (file_par==par::defaultString) {
268 
269  if (code=="CAMB" || code=="MGCAMB" || code=="MPTbreeze-v1")
270  m_Table_Pk_CAMB_MPTbreeze(code, NL, lgkk, lgPk, redshift, store_output, output_root, k_max);
271 
272  else if (code=="CLASS")
273  m_Table_Pk_CLASS(NL, lgkk, lgPk, redshift, store_output, output_root, k_max);
274 
275  else
276  ErrorCBL("the choosen code is not allowed!", "Table_PkCodes", "PkXi.cpp");
277  }
278 
279  else {
280  WarningMsgCBL("the input k_max parameters will not be used", "Table_PkCodes", "PkXi.cpp");
281  m_Table_Pk_parameterFile(code, file_par, NL, lgkk, lgPk, redshift, output_root);
282  }
283 
284 }
285 
286 
287 // =====================================================================================
288 
289 
290 void cbl::cosmology::Cosmology::Table_PkCodes (const std::string code, const bool NL, std::vector<std::vector<double>> &lgkk, std::vector<std::vector<double>> &lgPk, const std::vector<double> redshift, const bool store_output, const std::string output_root, const double k_max, const std::string file_par) const
291 {
292  if (code=="MPTbreeze-v1") {
293  if (m_sigma8<0)
294  ErrorCBL("sigma8<0! The function set_sigma8() can be used to set the value of sigma8!", "Table_PkCodes", "PkXi.cpp");
295  if (NL)
296  WarningMsgCBL("NL is ignored by MPTbreeze-v1, that provides in output the non-linear power spectrum", "Table_PkCodes", "PkXi.cpp");
297  }
298 
299  lgkk.resize(redshift.size(), vector<double>(1, 0.));
300  lgPk.resize(redshift.size(), vector<double>(1, 0.));
301 
302  if (file_par==par::defaultString) {
303 
304  if (code=="CAMB" || code=="MGCAMB" || code=="MPTbreeze-v1")
305  m_Table_Pk_CAMB_MPTbreeze(code, NL, lgkk, lgPk, redshift, store_output, output_root, k_max);
306 
307  else if (code=="CLASS")
308  m_Table_Pk_CLASS(NL, lgkk, lgPk, redshift, store_output, output_root, k_max);
309 
310  else
311  ErrorCBL("the choosen code is not allowed!", "Table_PkCodes", "PkXi.cpp");
312 
313  }
314 
315  else {
316  WarningMsgCBL("the input k_max parameters will not be used", "Table_PkCodes", "PkXi.cpp");
317  for (size_t zz=0; redshift.size(); zz++)
318  m_Table_Pk_parameterFile(code, file_par, NL, lgkk[zz], lgPk[zz], redshift[zz], output_root);
319  }
320 
321 }
322 
323 
324 // =====================================================================================
325 
326 
327 void cbl::cosmology::Cosmology::m_Table_Pk_CAMB_MPTbreeze (const std::string code, const bool NL, std::vector<double> &lgkk, std::vector<double> &lgPk, const double redshift, const bool store_output, const std::string output_root, const double k_max) const
328 {
329  lgkk.erase(lgkk.begin(), lgkk.end());
330  lgPk.erase(lgPk.begin(), lgPk.end());
331 
332  string dir_grid;
333  if (code=="CAMB" || code=="MGCAMB")
334  dir_grid = (NL) ? "output_nonlinear/" : "output_linear/";
335  else
336  dir_grid = "output/";
337 
338  string filename = "h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
339 
340  filename += (output_root=="test") ? "/" : "_"+output_root+"/";
341 
342  string file_par;
343 
344  cbl::Path path;
345  const string dirBZ = path.DirCosmo()+"/External/"+code+"/";
346  const string dirBZ_output = dirBZ+dir_grid+filename;
347  string fileBZ_in = dirBZ_output+"Pk.dat";
348  ifstream fin_BZ;
349  fin_BZ.open(fileBZ_in.c_str());
350 
351  const string nn = output_root+"_t"+conv(omp_get_thread_num(), par::fINT);
352  const string tot_output_root = dirBZ+nn;
353  string dir_output_root = dirBZ+nn;
354  dir_output_root = regex_replace(dir_output_root, std::regex("/"), "\\/");
355 
356 
357  // ------------------------------------------------------------------------------------------------------------
358  // --------- set the CAMB parameter file if it does not exist, or if the MPTbreeze one does not exist ---------
359  // ------------------------------------------------------------------------------------------------------------
360 
361  if (!fin_BZ && (code=="CAMB" || code=="MPTbreeze-v1")) {
362 
363  const string dirCAMB = path.DirCosmo()+"/External/CAMB/fortran/";
364 
365  file_par = dirCAMB+"../inifiles/params_"+nn+".ini";
366  if (system(("cp "+dirCAMB+"../inifiles/params.ini "+file_par).c_str())) {}
367 
368  string sed;
369 
370  sed = "sed '/test/s//"+dir_output_root+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par;
371  if (system(sed.c_str())) {}
372  sed = "sed '/do_nonlinear = 0/s//do_nonlinear = "+conv(NL, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
373  sed = "sed '/hubble = 70/s//hubble = "+conv(m_hh*100., par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
374  sed = "sed '/ombh2 = 0.0226/s//ombh2 = "+conv(m_Omega_baryon*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
375  sed = "sed '/omch2 = 0.112/s//omch2 = "+conv(m_Omega_CDM*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
376  sed = "sed '/omk = 0/s//omk = "+conv(m_Omega_k, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
377  sed = "sed '/omnuh2 = 0.00064/s//omnuh2 = "+conv(m_Omega_neutrinos*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
378  sed = "sed '/transfer_redshift(1) = 0/s//transfer_redshift(1) = "+conv(redshift, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
379  sed = "sed '/massless_neutrinos = 2.046/s//massless_neutrinos = "+conv(m_massless_neutrinos, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
380  sed = "sed '/massive_neutrinos = 1/s//massive_neutrinos = "+conv(m_massive_neutrinos, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
381  sed = "sed '/scalar_spectral_index(1) = 0.96/s//scalar_spectral_index(1) = "+conv(m_n_spec, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
382  sed = "sed '/w = -1/s//w = "+conv(m_w0, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
383  sed = "sed '/wa = 0/s//wa = "+conv(m_wa, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
384  if (m_scalar_amp>0) {
385  sed = "sed '/scalar_amp(1) = 2.1e-9/s//scalar_amp(1) = "+conv(m_scalar_amp, par::ee3)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
386  sed = "sed '/pivot_scalar = 0.05/s//pivot_scalar = "+conv(m_scalar_pivot, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
387  }
388  sed = "sed '/transfer_kmax = 2/s//transfer_kmax = "+conv(k_max, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
389  sed = "sed '/re_optical_depth = 0.09/s//re_optical_depth = "+conv(m_tau, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
390 
391  // ----------------------------
392  // --------- run CAMB ---------
393  // ----------------------------
394  if (system((dirCAMB+"camb "+file_par).c_str())) {}
395  }
396 
397  else if (!fin_BZ && code=="MGCAMB") {
398 
399  coutCBL << "Set the parameters for modified gravity in the file External/MGCAMB/params_MG.ini" << endl;
400 
401  file_par = dirBZ+"params_"+nn+".ini";
402  if (system(("cp "+dirBZ+"/params.ini "+file_par).c_str())) {}
403 
404  string sed;
405 
406  sed = "sed '/test/s//"+dir_output_root+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
407  sed = "sed '/do_nonlinear = 0/s//do_nonlinear = "+conv(NL, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
408  sed = "sed '/hubble = 70/s//hubble = "+conv(m_hh*100., par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
409  sed = "sed '/ombh2 = 0.0226/s//ombh2 = "+conv(m_Omega_baryon*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
410  sed = "sed '/omch2 = 0.112/s//omch2 = "+conv(m_Omega_CDM*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
411  sed = "sed '/omk = 0/s//omk = "+conv(m_Omega_k, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
412  sed = "sed '/omnuh2 = 0.00064/s//omnuh2 = "+conv(m_Omega_neutrinos*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
413  sed = "sed '/transfer_redshift(1) = 0/s//transfer_redshift(1) = "+conv(redshift, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
414  sed = "sed '/massless_neutrinos = 2.046/s//massless_neutrinos = "+conv(m_massless_neutrinos, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
415  sed = "sed '/massive_neutrinos = 1/s//massive_neutrinos = "+conv(m_massive_neutrinos, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
416  sed = "sed '/scalar_spectral_index(1) = 0.96/s//scalar_spectral_index(1) = "+conv(m_n_spec, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
417  sed = "sed '/w = -1/s//w = "+conv(m_w0, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
418  sed = "sed '/wa = 0/s//wa = "+conv(m_wa, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
419  if (m_scalar_amp>0) {
420  sed = "sed '/scalar_amp(1) = 2.1e-9/s//scalar_amp(1) = "+conv(m_scalar_amp, par::ee3)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
421  sed = "sed '/pivot_scalar = 0.05/s//pivot_scalar = "+conv(m_scalar_pivot, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
422  }
423  sed = "sed '/transfer_kmax = 2/s//transfer_kmax = "+conv(k_max, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
424  sed = "sed '/re_optical_depth = 0.09/s//re_optical_depth = "+conv(m_tau, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
425 
426  // ----------------------------
427  // -------- run MGCAMB --------
428  // ----------------------------
429  if (system((dirBZ+"camb "+file_par).c_str())) {}
430 
431  }
432 
433  // -----------------------------------------------------------------------------
434  // --------- check if the MPTbreeze table exists, if not run MPTbreeze ---------
435  // -----------------------------------------------------------------------------
436 
437  if (!fin_BZ && code=="MPTbreeze-v1") {
438 
439  // ---------------------------------
440  // --------- run MPTbreeze ---------
441  // ---------------------------------
442 
443  if (system((dirBZ+"mptbreeze -noverbose -camb "+file_par+" -fileTF "+tot_output_root+"_transfer_out.dat -sigma8 "+conv(m_sigma8, par::fDP6)+" -omegam "+conv(m_Omega_matter, par::fDP6)+" -ns "+conv(m_n_spec, par::fDP6)+" -w "+conv(m_w0, par::fDP6)+" -redshift "+conv(redshift, par::fDP6)+" -filePk "+tot_output_root+"_matterpower.dat").c_str())) {}
444 
445  }
446 
447  bool remove = false;
448 
449  if (!fin_BZ) {
450  if (store_output) {
451  if (system(("mkdir -p "+dirBZ_output).c_str())) {}
452  if (system(("mv "+tot_output_root+"_matterpower.dat "+dirBZ_output+"Pk.dat").c_str())) {}
453  if (system(("mv "+tot_output_root+"_transfer_out.dat "+dirBZ_output+"transfer_out.dat").c_str())) {}
454  }
455  else {
456  fileBZ_in = tot_output_root+"_matterpower.dat";
457  remove = true;
458  }
459  }
460 
461  if (system(("rm -f "+file_par+" "+tot_output_root+"_params.ini").c_str())) {}
462  fin_BZ.clear(); fin_BZ.close();
463 
464 
465  // ---------------------------------------
466  // --------- get the output P(k) ---------
467  // ---------------------------------------
468 
469  ifstream fin(fileBZ_in.c_str()); checkIO(fin, fileBZ_in);
470 
471  string line;
472 
473  while (getline(fin, line)) {
474  if (line.find("#")!=0) {
475  stringstream ss(line);
476  vector<double> num;
477  double aa;
478  while (ss>>aa) num.push_back(aa);
479 
480  if (num[0]>0 && num[1]>0) {
481 
482  lgkk.push_back(log10(num[0]));
483  if (code=="MPTbreeze-v1")
484  lgPk.push_back(log10(pow(2.*par::pi, 3)*(num[1]+num[2]+num[3])));
485  else
486  lgPk.push_back(log10(num[1]));
487 
488  }
489  }
490  }
491 
492  if (remove) if (system(("rm "+fileBZ_in+" "+tot_output_root+"_transfer_out.dat").c_str())) {}
493  fin.clear(); fin.close();
494 
495  if (lgkk.size()==0 || lgPk.size()==0)
496  ErrorCBL("lgkk.size()="+conv(lgkk.size(), par::fINT)+", lgPk.size()="+conv(lgPk.size(), par::fINT), "m_Table_Pk_CAMB_MPTbreeze", "PkXi.cpp");
497 
498 }
499 
500 
501 // =====================================================================================
502 
503 
504 void cbl::cosmology::Cosmology::m_Table_Pk_CAMB_MPTbreeze (const std::string code, const bool NL, std::vector<std::vector<double>> &lgkk, std::vector<std::vector<double>> &lgPk, const std::vector<double> redshift, const bool store_output, const std::string output_root, const double k_max) const
505 {
506  string dir_grid;
507  if (code=="CAMB" || code=="MGCAMB")
508  dir_grid = (NL) ? "output_nonlinear/" : "output_linear/";
509  else
510  dir_grid = "output/";
511 
512  string add_string;
513  add_string = (output_root=="test") ? "/" : "_"+output_root+"/";
514 
515  vector<double> zz = redshift;
516  std::sort(zz.begin(), zz.end(), greater<double>());
517  vector<string> filename_z(zz.size());
518  vector<string> filename_z_original(redshift.size());
519 
520  for (size_t ii=0; ii<zz.size(); ii++) {
521  filename_z[ii] = "h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(zz[ii], par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
522  filename_z[ii]+=add_string;
523  }
524 
525  for (size_t ii=0; ii<redshift.size(); ii++) {
526  filename_z_original[ii] = "h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift[ii], par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
527  filename_z_original[ii]+=add_string;
528  }
529 
530  string file_par;
531 
532  cbl::Path path;
533  const string dirBZ = path.DirCosmo()+"/External/"+code+"/";
534  const string dirBZ_output = dirBZ+dir_grid;
535 
536  // check if the power spectrum has been already stored for that cosmology and redshift
537  for (size_t ii=zz.size(); ii-->0;) {
538  string fileBZ_in = dirBZ+dir_grid+filename_z[ii]+"Pk.dat";
539  ifstream fin_BZ;
540  fin_BZ.open(fileBZ_in.c_str());
541  if (fin_BZ) { // remove the redshift from the vector if already computed
542  zz.erase(zz.begin()+ii);
543  filename_z.erase(filename_z.begin()+ii);
544  }
545  fin_BZ.clear(); fin_BZ.close();
546  }
547 
548  const string nn = output_root+"_t"+conv(omp_get_thread_num(), par::fINT);
549  const string tot_output_root = dirBZ+nn;
550  string dir_output_root = dirBZ+nn;
551  dir_output_root = regex_replace(dir_output_root, std::regex("/"), "\\/");
552 
553  // ------------------------------------------------------------------------------------------------------------
554  // --------- set the CAMB parameter file if it does not exist, or if the MPTbreeze one does not exist ---------
555  // ------------------------------------------------------------------------------------------------------------
556 
557  if (zz.size()>0 && (code=="CAMB" || code=="MPTbreeze-v1")) {
558 
559  cbl::Path path;
560  const string dirCAMB = path.DirCosmo()+"/External/CAMB/fortran/";
561 
562  file_par = dirCAMB+"../inifiles/params_"+nn+".ini";
563  if (system(("cp "+dirCAMB+"../inifiles/params.ini "+file_par).c_str())) {}
564 
565  string sed;
566 
567  sed = "sed '/test/s//"+dir_output_root+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par;
568  if (system(sed.c_str())) {}
569  sed = "sed '/do_nonlinear = 0/s//do_nonlinear = "+conv(NL, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
570  sed = "sed '/hubble = 70/s//hubble = "+conv(m_hh*100., par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
571  sed = "sed '/ombh2 = 0.0226/s//ombh2 = "+conv(m_Omega_baryon*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
572  sed = "sed '/omch2 = 0.112/s//omch2 = "+conv(m_Omega_CDM*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
573  sed = "sed '/omk = 0/s//omk = "+conv(m_Omega_k, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
574  sed = "sed '/omnuh2 = 0.00064/s//omnuh2 = "+conv(m_Omega_neutrinos*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
575  sed = "sed '/massless_neutrinos = 2.046/s//massless_neutrinos = "+conv(m_massless_neutrinos, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
576  sed = "sed '/massive_neutrinos = 1/s//massive_neutrinos = "+conv(m_massive_neutrinos, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
577  sed = "sed '/scalar_spectral_index(1) = 0.96/s//scalar_spectral_index(1) = "+conv(m_n_spec, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
578 
579  string transfer_z = "//";
580  string transfer_file = "//";
581  string transfer_matter = "//";
582  for (size_t ii=0; ii<zz.size(); ii++) {
583  transfer_z+="transfer_redshift("+std::to_string(ii+1)+") = "+conv(zz[ii], par::fDP6)+" \\\n";
584  transfer_file+="transfer_filename("+std::to_string(ii+1)+") = transfer_out_"+std::to_string(ii+1)+".dat \\\n";
585  transfer_matter+="transfer_matterpower("+std::to_string(ii+1)+") = matterpower_"+std::to_string(ii+1)+".dat \\\n";
586  }
587 
588  sed = "sed '/transfer_num_redshifts = 1/s//transfer_num_redshifts = "+conv(zz.size(), par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
589  sed = "sed '/transfer_redshift(1) = 0/s"+transfer_z+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
590  sed = "sed '/transfer_filename(1) = transfer_out.dat/s"+transfer_file+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
591  sed = "sed '/transfer_matterpower(1) = matterpower.dat/s"+transfer_matter+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
592 
593  sed = "sed '/w = -1/s//w = "+conv(m_w0, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
594  sed = "sed '/wa = 0/s//wa = "+conv(m_wa, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
595  if (m_scalar_amp>0) {
596  sed = "sed '/scalar_amp(1) = 2.1e-9/s//scalar_amp(1) = "+conv(m_scalar_amp, par::ee3)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
597  sed = "sed '/pivot_scalar = 0.05/s//pivot_scalar = "+conv(m_scalar_pivot, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
598  }
599  sed = "sed '/transfer_kmax = 2/s//transfer_kmax = "+conv(k_max, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
600  sed = "sed '/re_optical_depth = 0.09/s//re_optical_depth = "+conv(m_tau, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
601 
602  // ----------------------------
603  // --------- run CAMB ---------
604  // ----------------------------
605 
606  if (system((dirCAMB+"camb "+file_par).c_str())) {}
607  }
608 
609  else if (zz.size()>0 && code=="MGCAMB") {
610 
611  coutCBL << "Set the parameters for modified gravity in the file External/MGCAMB/params_MG.ini" << endl;
612 
613  file_par = dirBZ+"params_"+nn+".ini";
614  if (system(("cp "+dirBZ+"/params.ini "+file_par).c_str())) {}
615 
616  string sed;
617 
618  sed = "sed '/test/s//"+dir_output_root+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
619  sed = "sed '/do_nonlinear = 0/s//do_nonlinear = "+conv(NL, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
620  sed = "sed '/hubble = 70/s//hubble = "+conv(m_hh*100., par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
621  sed = "sed '/ombh2 = 0.0226/s//ombh2 = "+conv(m_Omega_baryon*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
622  sed = "sed '/omch2 = 0.112/s//omch2 = "+conv(m_Omega_CDM*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
623  sed = "sed '/omk = 0/s//omk = "+conv(m_Omega_k, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
624  sed = "sed '/omnuh2 = 0.00064/s//omnuh2 = "+conv(m_Omega_neutrinos*m_hh*m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
625 
626  string transfer_z = "//";
627  string transfer_file = "//";
628  string transfer_matter = "//";
629  for (size_t ii=0; ii<zz.size(); ii++) {
630  transfer_z+="transfer_redshift("+std::to_string(ii+1)+") = "+conv(zz[ii], par::fDP6)+" \\\n";
631  transfer_file+="transfer_filename("+std::to_string(ii+1)+") = transfer_out_"+std::to_string(ii+1)+".dat \\\n";
632  transfer_matter+="transfer_matterpower("+std::to_string(ii+1)+") = matterpower_"+std::to_string(ii+1)+".dat \\\n";
633  }
634  sed = "sed '/transfer_num_redshifts = 1/s//transfer_num_redshifts = "+conv(zz.size(), par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
635  sed = "sed '/transfer_redshift(1) = 0/s"+transfer_z+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
636  sed = "sed '/transfer_filename(1) = transfer_out.dat/s"+transfer_file+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
637  sed = "sed '/transfer_matterpower(1) = matterpower.dat/s"+transfer_matter+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
638 
639  sed = "sed '/massless_neutrinos = 2.046/s//massless_neutrinos = "+conv(m_massless_neutrinos, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
640  sed = "sed '/massive_neutrinos = 1/s//massive_neutrinos = "+conv(m_massive_neutrinos, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
641  sed = "sed '/scalar_spectral_index(1) = 0.96/s//scalar_spectral_index(1) = "+conv(m_n_spec, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
642  sed = "sed '/w = -1/s//w = "+conv(m_w0, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
643  sed = "sed '/wa = 0/s//wa = "+conv(m_wa, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
644  if (m_scalar_amp>0) {
645  sed = "sed '/scalar_amp(1) = 2.1e-9/s//scalar_amp(1) = "+conv(m_scalar_amp, par::ee3)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
646  sed = "sed '/pivot_scalar = 0.05/s//pivot_scalar = "+conv(m_scalar_pivot, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
647  }
648  sed = "sed '/transfer_kmax = 2/s//transfer_kmax = "+conv(k_max, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
649  sed = "sed '/re_optical_depth = 0.09/s//re_optical_depth = "+conv(m_tau, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
650 
651  // ----------------------------
652  // -------- run MGCAMB --------
653  // ----------------------------
654 
655  if (system((dirBZ+"camb "+file_par).c_str())) {}
656 
657  }
658 
659  // -----------------------------------------------------------------------------
660  // --------- check if the MPTbreeze table exists, if not run MPTbreeze ---------
661  // -----------------------------------------------------------------------------
662 
663  if (zz.size()>0 && code=="MPTbreeze-v1") {
664 
665  // ---------------------------------
666  // --------- run MPTbreeze ---------
667  // ---------------------------------
668 
669  for (size_t ii=0; ii<zz.size(); ii++) {
670  if (system((dirBZ+"mptbreeze -noverbose -camb "+file_par+" -fileTF "+tot_output_root+"_transfer_out_"+std::to_string(ii+1)+".dat -sigma8 "+conv(m_sigma8, par::fDP6)+" -omegam "+conv(m_Omega_matter, par::fDP6)+" -ns "+conv(m_n_spec, par::fDP6)+" -w "+conv(m_w0, par::fDP6)+" -redshift "+conv(zz[ii], par::fDP6)+" -filePk "+tot_output_root+"_matterpower_"+std::to_string(ii+1)+".dat").c_str())) {}
671  }
672 
673  }
674 
675  if (store_output)
676  for (size_t ii=0; ii<zz.size(); ii++) {
677  if (system(("mkdir -p "+dirBZ_output+filename_z[ii]).c_str())) {}
678  if (system(("mv "+tot_output_root+"_matterpower_"+std::to_string(ii+1)+".dat "+dirBZ_output+filename_z[ii]+"Pk.dat").c_str())) {}
679  if (system(("mv "+tot_output_root+"_transfer_out_"+std::to_string(ii+1)+".dat "+dirBZ_output+filename_z[ii]+"transfer_out.dat").c_str())) {}
680  }
681 
682  if (system(("rm -f "+file_par+" "+tot_output_root+"_params.ini").c_str())) {}
683 
684  // ---------------------------------------
685  // --------- get the output P(k) ---------
686  // ---------------------------------------
687 
688  int counter = 0;
689 
690  for (size_t ii=redshift.size(); ii-->0;) {
691 
692  lgkk[ii].erase(lgkk[ii].begin(), lgkk[ii].end());
693  lgPk[ii].erase(lgPk[ii].begin(), lgPk[ii].end());
694 
695  string fileBZ_in = dirBZ+dir_grid+filename_z_original[ii]+"Pk.dat";
696  ifstream fin_check(fileBZ_in.c_str());
697 
698  if (fin_check and !store_output) counter++;
699  if (!fin_check) fileBZ_in = tot_output_root+"_matterpower_"+std::to_string(redshift.size()-ii-counter)+".dat";
700 
701  ifstream fin(fileBZ_in.c_str());
702  string line;
703 
704  while (getline(fin, line)) {
705 
706  if (line.find("#")!=0) {
707  stringstream ss(line);
708  vector<double> num;
709  double aa;
710  while (ss>>aa) num.push_back(aa);
711 
712  if (num[0]>0 && num[1]>0) {
713 
714  lgkk[ii].push_back(log10(num[0]));
715  if (code=="MPTbreeze-v1")
716  lgPk[ii].push_back(log10(pow(2.*par::pi, 3)*(num[1]+num[2]+num[3])));
717  else
718  lgPk[ii].push_back(log10(num[1]));
719 
720  }
721  }
722  }
723 
724  if (!fin_check) if (system(("rm "+fileBZ_in+" "+tot_output_root+"_transfer_out_"+std::to_string(redshift.size()-ii-counter)+".dat").c_str())) {}
725 
726  fin_check.clear(); fin_check.close();
727  fin.clear(); fin.close();
728 
729  if (lgkk[ii].size()==0 || lgPk[ii].size()==0)
730  ErrorCBL("lgkk.size()="+conv(lgkk[ii].size(), par::fINT)+", lgPk.size()="+conv(lgPk[ii].size(), par::fINT), "m_Table_Pk_CAMB_MPTbreeze", "PkXi.cpp");
731 
732  }
733 
734 }
735 
736 
737 // =====================================================================================
738 
739 
740 void cbl::cosmology::Cosmology::m_Table_Pk_CLASS (const bool NL, std::vector<double> &lgkk, std::vector<double> &lgPk, const double redshift, const bool store_output, const std::string output_root, const double k_max) const
741 {
742  lgkk.erase(lgkk.begin(), lgkk.end());
743  lgPk.erase(lgPk.begin(), lgPk.end());
744 
745  cbl::Path path;
746  const string dirC = path.DirCosmo()+"/External/CLASS/";
747 
748  string dir_grid = (NL) ? dirC+"output_nonlinear/" : dirC+"output_linear/";
749 
750  dir_grid += "h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
751 
752  dir_grid += (output_root=="test") ? "/" : "_"+output_root+"/";
753 
754  const string nn = output_root+"_t"+conv(omp_get_thread_num(), par::fINT);
755  const string dir_output = dirC+nn+"_00";
756 
757  string file_in = dir_grid+"Pk.dat";
758  ifstream fin;
759  fin.open(file_in.c_str());
760 
761  string file_par;
762 
763  if (!fin) {
764 
765  // --------------------------------------------------------------------------
766  // --------- set the cosmological parameters in the file params.ini ---------
767  // --------------------------------------------------------------------------
768 
769  string dir_output_root = dirC+nn;
770  dir_output_root = regex_replace(dir_output_root, std::regex("/"), "\\/");
771 
772  file_par = dirC+"params_"+nn+".ini";
773  if (system(("cp "+dirC+"params.ini "+file_par).c_str())) {}
774 
775  string sed;
776 
777  sed = "sed '/output\\/test_/s//"+dir_output_root+"_/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
778 
779  if (NL) {
780  sed = "sed '/non linear =/s//non linear = halofit/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+"";
781  if (system(sed.c_str())) {}
782  }
783 
784  sed = "sed '/h = 0.67556/s//h = "+conv(m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
785  sed = "sed '/Omega_b = 0.022032/s//Omega_b = "+conv(m_Omega_baryon, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
786  sed = "sed '/Omega_cdm = 0.12038/s//Omega_cdm = "+conv(m_Omega_CDM, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
787  sed = "sed '/N_ncdm = 0/s//N_ncdm = "+conv(m_massive_neutrinos, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
788  sed = "sed '/Omega_ncdm = 0/s//Omega_ncdm = "+conv(m_Omega_neutrinos, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
789 
790  if (m_Omega_neutrinos>0) { // check!!!
791  sed = "sed '/m_ncdm = 0.04, 0.04, 0.04/s//#m_ncdm/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+"";
792  if (system(sed.c_str())) {}
793  }
794 
795  sed = "sed '/Omega_Lambda = 0.7/s//Omega_Lambda = "+conv(m_Omega_DE, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
796  sed = "sed '/Omega_k = 0./s//Omega_k = "+conv(m_Omega_k, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
797  sed = "sed '/z_pk = 0/s//z_pk = "+conv(redshift, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
798  sed = "sed '/N_ur = 3.046/s//N_ur = "+conv(m_massless_neutrinos, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
799  sed = "sed '/N_ncdm = 0/s//N_ncdm = "+conv(m_massive_neutrinos, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
800  sed = "sed '/n_s = 0.9619/s//n_s = "+conv(m_n_spec, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
801 
802  double w00 = max(-0.999,m_w0); // check!!!
803  sed = "sed '/w0_fld = -0.9/s//w0_fld = "+conv(w00, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
804 
805  sed = "sed '/wa_fld = 0./s//wa_fld = "+conv(m_wa, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
806  if (m_scalar_amp>0) {
807  sed = "sed '/A_s = 2.215e-9/s//A_s = "+conv(m_scalar_amp, par::ee3)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
808  sed = "sed '/k_pivot = 0.05/s//k_pivot = "+conv(m_scalar_pivot, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
809  }
810  sed = "sed '/P_k_max_h\\/Mpc = 1/s//P_k_max_h\\/Mpc = "+conv(k_max, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
811  sed = "sed '/tau_reio = 0.0925/s//tau_reio = "+conv(m_tau, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
812 
813  // -----------------------------
814  // --------- run CLASS ---------
815  // -----------------------------
816 
817  if (system((dirC+"/class "+file_par).c_str())) {}
818  if (system(("rm -f "+file_par).c_str())) {}
819  }
820 
821  fin.clear(); fin.close();
822 
823  bool remove = false;
824 
825  if (store_output) {
826  if (system(("mkdir -p "+dir_grid).c_str())) {}
827  if (system(("mv "+dir_output+"_pk"+((NL) ? "_nl.dat" : ".dat")+" "+dir_grid+"Pk.dat").c_str())) {}
828  }
829 
830  else {
831  file_in = dir_output+"_pk"+((NL) ? "_nl.dat" : ".dat");
832  remove = true;
833  }
834 
835  // ---------------------------------------
836  // --------- get the output P(k) ---------
837  // ---------------------------------------
838 
839  fin.open(file_in.c_str()); checkIO(fin, file_in);
840  string line;
841 
842  while (getline(fin, line)) {
843  stringstream ss(line);
844  vector<double> num;
845  double aa;
846  while (ss>>aa) num.push_back(aa);
847 
848  if (num[0]>0 && num[1]>0) {
849  lgkk.push_back(log10(num[0]));
850  lgPk.push_back(log10(num[1]));
851  }
852  }
853 
854  if (remove) if (system(("rm -f "+dir_output+"_pk*.dat").c_str())) {}
855  fin.clear(); fin.close();
856 
857  if (lgkk.size()==0 || lgPk.size()==0)
858  ErrorCBL("lgkk.size()="+conv(lgkk.size(), par::fINT)+", lgPk.size()="+conv(lgPk.size(), par::fINT), "m_Table_Pk_CLASS", "PkXi.cpp");
859 
860 }
861 
862 
863 // =====================================================================================
864 
865 
866 void cbl::cosmology::Cosmology::m_Table_Pk_CLASS (const bool NL, std::vector<std::vector<double>> &lgkk, std::vector<std::vector<double>> &lgPk, const std::vector<double> redshift, const bool store_output, const std::string output_root, const double k_max) const
867 {
868  cbl::Path path;
869  const string dirC = path.DirCosmo()+"/External/CLASS/";
870  string dir_grid = (NL) ? dirC+"output_nonlinear/" : dirC+"output_linear/";
871 
872  vector<string> filename_z_original(redshift.size());
873  vector<string> dir_output_z_original(redshift.size());
874 
875 
876  for (size_t ii=0; ii<redshift.size(); ii++) {
877  filename_z_original[ii] = dir_grid+"h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift[ii], par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
878  filename_z_original[ii] += (output_root=="test") ? "/" : "_"+output_root+"/";
879  dir_output_z_original[ii] = filename_z_original[ii];
880  }
881 
882  vector<double> zz = redshift;
883 
884  for (size_t ii=zz.size(); ii-->0;) {
885  string file_in = dir_output_z_original[ii]+"Pk.dat";
886  ifstream fin;
887  fin.open(file_in.c_str());
888  if (fin) zz.erase(zz.begin()+ii); // remove the redshift from the vector if already computed
889  fin.clear(); fin.close();
890  }
891 
892  vector<string> filename_z(zz.size());
893  vector<string> dir_output_z(zz.size());
894 
895  for (size_t ii=0; ii<zz.size(); ii++) {
896  filename_z[ii] = dir_grid+"h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(zz[ii], par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
897  filename_z[ii] += (output_root=="test") ? "/" : "_"+output_root+"/";
898  dir_output_z[ii] = filename_z[ii];
899  }
900 
901  const string nn = output_root+"_t"+conv(omp_get_thread_num(), par::fINT);
902 
903  if (zz.size()>0) {
904 
905  // --------------------------------------------------------------------------
906  // --------- set the cosmological parameters in the file params.ini ---------
907  // --------------------------------------------------------------------------
908 
909  string dir_output_root = dirC+nn;
910  dir_output_root = regex_replace(dir_output_root, std::regex("/"), "\\/");
911 
912  string file_par = dirC+"params_"+output_root+"_t"+conv(omp_get_thread_num(), par::fINT)+".ini";
913  if (system(("cp "+dirC+"params.ini "+file_par).c_str())) {}
914 
915  string sed;
916 
917  sed = "sed '/output\\/test_/s//"+dir_output_root+"_/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
918 
919  if (NL) {
920  sed = "sed '/non linear =/s//non linear = halofit/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+"";
921  if (system(sed.c_str())) {}
922  }
923 
924  sed = "sed '/h = 0.67556/s//h = "+conv(m_hh, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
925  sed = "sed '/Omega_b = 0.022032/s//Omega_b = "+conv(m_Omega_baryon, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
926  sed = "sed '/Omega_cdm = 0.12038/s//Omega_cdm = "+conv(m_Omega_CDM, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
927  sed = "sed '/N_ncdm = 0/s//N_ncdm = "+conv(m_massive_neutrinos, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
928  sed = "sed '/Omega_ncdm = 0/s//Omega_ncdm = "+conv(m_Omega_neutrinos, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
929 
930  if (m_Omega_neutrinos>0) { // check!!!
931  sed = "sed '/m_ncdm = 0.04, 0.04, 0.04/s//#m_ncdm/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+"";
932  if (system(sed.c_str())) {}
933  }
934 
935  sed = "sed '/Omega_Lambda = 0.7/s//Omega_Lambda = "+conv(m_Omega_DE, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
936  sed = "sed '/Omega_k = 0./s//Omega_k = "+conv(m_Omega_k, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par+""; if (system(sed.c_str())) {}
937 
938  string z_pk = "z_pk = ";
939  for (size_t ii=0; ii<zz.size(); ii++) {
940  z_pk+=conv(zz[ii], par::fDP6);
941  if (ii<zz.size()-1) z_pk += ",";
942  }
943 
944  sed = "sed '/z_pk = 0/s//"+z_pk+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
945 
946  sed = "sed '/N_ur = 3.046/s//N_ur = "+conv(m_massless_neutrinos, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
947  sed = "sed '/N_ncdm = 0/s//N_ncdm = "+conv(m_massive_neutrinos, par::fINT)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
948  sed = "sed '/n_s = 0.9619/s//n_s = "+conv(m_n_spec, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
949 
950  double w00 = max(-0.999,m_w0); // check!!!
951  sed = "sed '/w0_fld = -0.9/s//w0_fld = "+conv(w00, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
952 
953  sed = "sed '/wa_fld = 0./s//wa_fld = "+conv(m_wa, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
954  if (m_scalar_amp>0) {
955  sed = "sed '/A_s = 2.215e-9/s//A_s = "+conv(m_scalar_amp, par::ee3)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
956  sed = "sed '/k_pivot = 0.05/s//k_pivot = "+conv(m_scalar_pivot, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
957  }
958  sed = "sed '/P_k_max_h\\/Mpc = 1/s//P_k_max_h\\/Mpc = "+conv(k_max, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
959  sed = "sed '/tau_reio = 0.0925/s//tau_reio = "+conv(m_tau, par::fDP6)+"/g' "+file_par+" > temp_"+nn+"; mv temp_"+nn+" "+file_par; if (system(sed.c_str())) {}
960 
961 
962  // -----------------------------
963  // --------- run CLASS ---------
964  // -----------------------------
965 
966  if (system((dirC+"/class "+file_par).c_str())) {}
967  if (system(("rm -f "+file_par).c_str())) {}
968 
969  if (store_output) {
970  for (size_t ii=0; ii<zz.size(); ii++) {
971  if (system(("mkdir -p "+dir_output_z[ii]).c_str())) {}
972 
973  if (zz.size()==1) {if (system(("mv "+dirC+nn+"_00_pk"+((NL) ? "_00_nl.dat" : ".dat")+" "+dir_output_z[ii]+"Pk.dat").c_str())) {}}
974  else {if (system(("mv "+dirC+nn+"_00_z"+std::to_string(ii+1)+((NL) ? "_pk_nl.dat" : "_pk.dat")+" "+dir_output_z[ii]+"Pk.dat").c_str())) {}}
975  }
976  }
977 
978  }
979 
980  // ---------------------------------------
981  // --------- get the output P(k) ---------
982  // ---------------------------------------
983 
984  int counter = 0;
985 
986  for (size_t ii=0; ii<redshift.size(); ii++) {
987 
988  lgkk[ii].erase(lgkk[ii].begin(), lgkk[ii].end());
989  lgPk[ii].erase(lgPk[ii].begin(), lgPk[ii].end());
990 
991  string file_in = dir_output_z_original[ii]+"Pk.dat";
992  ifstream fin_check(file_in.c_str());
993 
994  if (fin_check and !store_output) counter++;
995  if (!fin_check) {
996  if (zz.size()==1) file_in = dirC+nn+"_00_pk"+((NL) ? "_00_nl.dat" : ".dat");
997  else file_in = dirC+nn+"_00_z"+std::to_string(ii+1-counter)+((NL) ? "_pk_nl.dat" : "_pk.dat");
998  }
999 
1000  ifstream fin(file_in.c_str());
1001  string line;
1002 
1003  while (getline(fin, line)) {
1004  stringstream ss(line);
1005  vector<double> num;
1006  double aa;
1007  while (ss>>aa) num.push_back(aa);
1008 
1009  if (num[0]>0 && num[1]>0) {
1010  lgkk[ii].push_back(log10(num[0]));
1011  lgPk[ii].push_back(log10(num[1]));
1012  }
1013  }
1014 
1015  fin_check.clear(); fin_check.close();
1016  fin.clear(); fin.close();
1017 
1018  if (lgkk[ii].size()==0 || lgPk[ii].size()==0)
1019  ErrorCBL("lgkk.size()="+conv(lgkk[ii].size(), par::fINT)+", lgPk.size()="+conv(lgPk[ii].size(), par::fINT), "m_Table_Pk_CLASS", "PkXi.cpp");
1020  }
1021 
1022  if (system(("rm -f "+dirC+nn+"_00_*.dat").c_str())) {} // to remove also *_cb.dat files
1023 
1024 }
1025 
1026 
1027 // =====================================================================================
1028 
1029 
1030 void cbl::cosmology::Cosmology::remove_output_Pk_tables (const string code, const bool NL, const double redshift, const string output_root) const
1031 {
1032  string dir_grid;
1033  if (code=="CAMB" || code=="MGCAMB" || code=="CLASS")
1034  dir_grid = (NL) ? "output_nonlinear/" : "output_linear/";
1035  else
1036  dir_grid = "output/";
1037 
1038  dir_grid += "h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
1039 
1040  if (output_root!="test") dir_grid += "_"+output_root;
1041 
1042  string dir_output;
1043 
1044  cbl::Path path;
1045  if (code=="CAMB") dir_output = path.DirCosmo()+"/External/CAMB/"+dir_grid;
1046  if (code=="MGCAMB") dir_output = path.DirCosmo()+"/External/MGCAMB/"+dir_grid;
1047  else if (code=="CLASS") dir_output = path.DirCosmo()+"/External/CLASS/"+dir_grid;
1048  else dir_output = path.DirCosmo()+"/External/MPTbreeze-v1/"+dir_grid;
1049 
1050  if (system(("rm -rf "+dir_output+" > /dev/null 2>&1").c_str())) {}
1051 
1052 }
1053 
1054 
1055 // =====================================================================================
1056 
1057 
1058 void cbl::cosmology::Cosmology::m_Table_Pk_parameterFile (const std::string code, const std::string file_par, const bool NL, std::vector<double> &lgkk, std::vector<double> &lgPk, const double redshift, const std::string output_root) const
1059 {
1060  WarningMsgCBL("Check the consistency between the output_root/root parameter given in input and the one set in the parameter file", "m_Table_Pk_parameterFile", "PkXi.cpp");
1061 
1062  lgkk.erase(lgkk.begin(), lgkk.end());
1063  lgPk.erase(lgPk.begin(), lgPk.end());
1064 
1065  cbl::Path path;
1066  const string dir_loc = path.fullpath(path.DirLoc());
1067  const string dirB = (code=="CAMB" || code=="MPTbreeze-v1") ? path.DirCosmo()+"/External/CAMB/fortran/"
1068  : path.DirCosmo()+"/External/"+code+"/";
1069 
1070  const string dir_output = (code=="CAMB" || code=="MPTbreeze-v1") ? dirB+"../output_ParameterFiles/"+file_par+"/"
1071  : dirB+"output_ParameterFiles/"+file_par+"/";
1072 
1073  if (system(("mkdir -p "+dir_output).c_str())) {}
1074 
1075  const string file_out = dir_output+"Pk.dat";
1076 
1077  const string dirPar = (code=="CAMB" || code=="MPTbreeze-v1") ? "../inifiles/" : "./";
1078  const string File_par = dirB+dirPar+file_par;
1079 
1080  ifstream check;
1081  check.open(File_par.c_str());
1082  if (!check) ErrorCBL("Parameter file not found! Make sure that this file is located in CosmoBolognaLib/External/CAMB/inifiles/ for CAMB and MPTbreeze-v1, and in CosmoBolognaLib/External/_NameBoltzmannSolver_/ for the remaining cases", "m_Table_Pk_parameterFile", "PkXi.cpp");
1083  check.clear(); check.close();
1084 
1085  ifstream fin;
1086  fin.open(file_out.c_str());
1087 
1088  if (!fin) {
1089 
1090  if (chdir(dirB.c_str())) {}
1091 
1092  // --------------------------------------------
1093  // --------- run the Boltzmann solver ---------
1094  // --------------------------------------------
1095 
1096  if (code=="CAMB" || code=="MGCAMB" || code=="MPTbreeze-v1") {
1097  if (system(("./camb "+File_par).c_str())) {}
1098  if (system(("mv "+output_root+"_matterpower*dat "+dir_output+"Pk.dat").c_str())) {}
1099 
1100  if (code=="MPTbreeze-v1") {
1101  WarningMsgCBL("Check the consistency of the input redshift with the one set in the parameter file", "m_Table_Pk_parameterFile", "PkXi.cpp");
1102  if (chdir((path.DirCosmo()+"/External/"+code).c_str())) {}
1103 
1104  if (system(("./mptbreeze -noverbose -camb ../../CAMB/"+File_par+" -fileTF ../../CAMB/"+output_root+"_transfer_out.dat -sigma8 "+conv(m_sigma8, par::fDP6)+" -omegam "+conv(m_Omega_matter, par::fDP6)+" -ns "+conv(m_n_spec, par::fDP6)+" -w "+conv(m_w0, par::fDP6)+" -redshift "+conv(redshift, par::fDP6)+" -filePk ../../CAMB/"+output_root+"_matterpower.dat").c_str())) {}
1105  if (system(("mv ../../CAMB/"+output_root+"_matterpower*dat "+dir_output+"Pk.dat").c_str())) {}
1106  }
1107 
1108  }
1109 
1110  else if (code=="CLASS") {
1111  WarningMsgCBL("Check the consistency of the input NL value with the one set in the parameter file", "m_Table_Pk_parameterFile", "PkXi.cpp");
1112  if (system(("./class "+File_par).c_str())) {}
1113  if (system(("mv "+output_root+"pk"+((NL) ? ".dat" : "_nl.dat")+" "+dir_output+"Pk.dat").c_str())) {}
1114  }
1115 
1116  else ErrorCBL("the chosen code is not allowed!", "m_Table_parameterFile", "PkXi.cpp");
1117 
1118  if (system(("rm -f "+output_root+"*_params.ini "+output_root+"*_transfer*" ).c_str())) {}
1119  if (chdir(dir_loc.c_str())) {}
1120 
1121  }
1122 
1123  fin.clear(); fin.close();
1124 
1125  // ---------------------------------------
1126  // --------- get the output P(k) ---------
1127  // ---------------------------------------
1128 
1129  fin.open(file_out.c_str()); checkIO(fin, file_out);
1130  string line;
1131 
1132  while (getline(fin, line)) {
1133  if (line.find("#")!=0) {
1134  stringstream ss(line);
1135  vector<double> num;
1136  double aa;
1137  while (ss>>aa) num.push_back(aa);
1138 
1139  if (num[0]>0 && num[1]>0) {
1140 
1141  lgkk.push_back(log10(num[0]));
1142  if (code=="MPTbreeze-v1")
1143  lgPk.push_back(log10(pow(2.*par::pi, 3)*(num[1]+num[2]+num[3])));
1144  else
1145  lgPk.push_back(log10(num[1]));
1146 
1147  }
1148  }
1149  }
1150 
1151  fin.clear(); fin.close();
1152 
1153  if (lgkk.size()==0 || lgPk.size()==0)
1154  ErrorCBL("lgkk.size()="+conv(lgkk.size(), par::fINT)+", lgPk.size()="+conv(lgPk.size(), par::fINT), "Table_PkCodes", "PkXi.cpp");
1155 
1156  if (chdir(dir_loc.c_str())) {}
1157 }
1158 
1159 
1160 
1161 // =====================================================================================
1162 
1163 
1164 void cbl::cosmology::Cosmology::Table_XiCodes (const std::string code, const bool NL, std::vector<double> &rr, std::vector<double> &xi, const double redshift, const bool store_output, const std::string output_root, const double k_max, const std::string file_par) const
1165 {
1166  vector<double> _lgkk, _lgPk;
1167  Table_PkCodes(code, NL, _lgkk, _lgPk, redshift, store_output, output_root, k_max, file_par);
1168 
1169  if (_lgkk.size()==0 || _lgPk.size()==0)
1170  ErrorCBL("_lgkk.size()="+conv(_lgkk.size(), par::fINT)+", _lgPk.size()="+conv(_lgPk.size(), par::fINT), "Table_XiCodes", "PkXi.cpp");
1171 
1172  // if the size of the input dataset is less than 1000, interpolate the values in a denser grid
1173  //WarningMsgCBL("the number of k-P(k) grid points will be increased to 1000, by interpolation", "Table_XiCodes", "PkXi.cpp");
1174  vector<double> lgkk, lgPk;
1175  if (_lgkk.size()<1000) {
1176  lgkk = linear_bin_vector(1000, min(1.e-4, Min(_lgkk)), max(1., Max(_lgkk)));
1177  for (size_t i=0; i<lgkk.size(); ++i)
1178  lgPk.emplace_back(interpolated(lgkk[i], _lgkk, _lgPk, "Linear"));
1179  }
1180  else {
1181  lgkk = _lgkk;
1182  lgPk = _lgPk;
1183  }
1184  cbl::Path path;
1185  const string dir_loc = path.DirLoc();
1186  const string dir_cosmo = path.DirCosmo();
1187 
1188  string dir_grid;
1189  if (code=="MPTbreeze-v1") dir_grid = "output/";
1190  else if (NL==false) dir_grid = "output_linear/";
1191  else dir_grid = "output_nonlinear/";
1192 
1193  dir_grid += "h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
1194 
1195  if (output_root!="test") dir_grid += "_"+output_root;
1196 
1197  string dir = dir_cosmo+"/External/"+code+"/";
1198  string dirFFT = dir_cosmo+"/External/fftlog-f90-master/";
1199  string dir_output = dir+dir_grid+"/";
1200 
1201  string file_in = (code=="MPTbreeze-v1") ? dir_output+"Pk2.dat" : dir_output+"Pk.dat";
1202  string file_out = dir_output+"Xi.dat";
1203 
1204  const string mkdir = "mkdir -p "+dir_output; if (system(mkdir.c_str())) {}
1205 
1206  ifstream fin;
1207  fin.open(file_out.c_str());
1208 
1209  if (!fin) {
1210  vector<double> kk = lgkk, Pk = lgPk;
1211  for_each( kk.begin(), kk.end(), [] (double &vv) { vv = pow(10., vv); } );
1212  for_each( Pk.begin(), Pk.end(), [] (double &vv) { vv = pow(10., vv); } );
1213 
1214  wrapper::fftlog::transform_FFTlog(rr, xi, 1, kk, Pk);
1215 
1216  ofstream fout(file_out.c_str()); checkIO(fout, file_out);
1217 
1218  for (size_t i=0; i<rr.size(); ++i)
1219  fout << rr[i] << " " << xi[i] << endl;
1220 
1221  fout.clear(); fout.close();
1222  }
1223 
1224  if (chdir(dir_loc.c_str())) {}
1225 
1226  fin.clear(); fin.close();
1227 
1228  rr.erase(rr.begin(), rr.end());
1229  xi.erase(xi.begin(), xi.end());
1230 
1231  fin.open(file_out.c_str());
1232 
1233  double RR, XI;
1234  while (fin >> RR >> XI) {
1235  rr.push_back(RR);
1236  xi.push_back(XI);
1237  }
1238 
1239  fin.clear(); fin.close();
1240 
1241  if (!store_output) {
1242  string RM = "rm -rf "+dir_grid;
1243  if (system(RM.c_str())) {}
1244  }
1245 }
1246 
1247 
1248 // =====================================================================================
1249 
1250 
1251 void cbl::cosmology::Cosmology::Pk_0 (const std::string method_Pk, const double redshift, const bool store_output, const std::string output_root, const double k_min, const double k_max, const double prec, const std::string file_par)
1252 {
1253  if (m_sigma8<0) ErrorCBL("sigma8<0!", "Pk_0", "PkXi.cpp");
1254 
1255  double RR = 8.; // sigma_8 = sigma(8Mpc/h)
1256  double RHO = rho_m(0., true);
1257  double MM = Mass(RR, RHO);
1258 
1259  bool NL = false;
1260  double Int = -1.;
1261  double error = -1.;
1262 
1263  if (method_Pk=="EisensteinHu") {
1264 
1265  EisensteinHu eh;
1266 
1267  eh.TFmdm_set_cosm(m_Omega_matter, m_Omega_baryon, m_Omega_neutrinos, m_massive_neutrinos, m_Omega_DE, m_hh, redshift, m_scalar_amp, m_scalar_pivot, m_n_spec);
1268 
1269  auto func = [&] (const double kk)
1270  {
1271  return pow(TopHat_WF(kk*RR)*kk, 2)*eh.Pk(kk);
1272  };
1273 
1274  Int = wrapper::gsl::GSL_integrate_qag(func, k_min, k_max, prec);
1275  }
1276 
1277  else if (method_Pk == "CAMB" || method_Pk=="CAMB_wrapper" || method_Pk=="MGCAMB" || method_Pk=="MPTbreeze-v1" || method_Pk=="CLASS") {
1278 
1279  vector<double> lgkk, lgPk;
1280 
1281  if (method_Pk == "CAMB_wrapper") {
1282 
1283  const int npoints = 500;
1284  vector<double> Pk = wrapper::camb::Pk_CAMB (false, redshift, k_min, k_max, npoints, m_Omega_baryon*m_hh*m_hh, m_Omega_CDM*m_hh*m_hh, m_Omega_neutrinos*m_hh*m_hh, m_massless_neutrinos, m_massive_neutrinos, m_Omega_k, m_hh*100., m_n_spec, m_scalar_amp, m_scalar_pivot, m_w0, m_wa, m_tau);
1285  vector<double> kk = logarithmic_bin_vector(npoints, k_min, k_max);
1286 
1287  for (size_t i=0; i<Pk.size(); i++) {
1288  lgkk.emplace_back(log10(kk[i]/m_hh));
1289  lgPk.emplace_back(log10(Pk[i]));
1290  }
1291 
1292  }
1293 
1294  else
1295  Table_PkCodes(method_Pk, NL, lgkk, lgPk, redshift, store_output, output_root, k_max, file_par);
1296 
1297  int limit_size = 1000;
1298  gsl_integration_workspace *ww = gsl_integration_workspace_alloc(limit_size);
1299  gsl_function Func;
1300 
1301  glob::STR_SSM str;
1302  str.unit = true;
1303  str.hh = m_hh;
1304  str.n_spec = m_n_spec;
1305  str.mass = MM;
1306  str.rho = RHO;
1307  str.lgkk = lgkk;
1308  str.lgPk = lgPk;
1309 
1310  Func.function = &glob::func_SSM_GSL;
1311  Func.params = &str;
1312  gsl_integration_qag(&Func, k_min, k_max, 0., prec, limit_size, 6, ww, &Int, &error);
1313  gsl_integration_workspace_free (ww);
1314  }
1315 
1316  else ErrorCBL("method_Pk is wrong!", "Pk_0", "PkXi.cpp");
1317 
1318  const double D_N = DN(redshift);
1319 
1320  if (method_Pk=="EisensteinHu") m_Pk0_EH = 2.*pow(par::pi*m_sigma8,2)/Int*D_N*D_N;
1321  if (method_Pk=="CAMB" || method_Pk=="CAMB_wrapper" || method_Pk=="MGCAMB") m_Pk0_CAMB = 2.*pow(par::pi*m_sigma8,2)/Int*D_N*D_N;
1322  if (method_Pk=="MPTbreeze-v1") m_Pk0_MPTbreeze = 2.*pow(par::pi*m_sigma8,2)/Int*D_N*D_N;
1323  if (method_Pk=="CLASS") m_Pk0_CLASS = 2.*pow(par::pi*m_sigma8,2)/Int*D_N*D_N;
1324 
1325 }
1326 
1327 
1328 // =====================================================================================
1329 
1330 
1331 std::vector<double> cbl::cosmology::Cosmology::Pk_matter (const std::vector<double> kk, const std::string method_Pk, const bool NL, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec, const std::string file_par, const bool unit1)
1332 {
1333  double fact1 = (m_unit || unit1) ? 1. : 1./m_hh;
1334  double fact2 = pow(fact1, 3);
1335 
1336  vector<double> newk = kk;
1337  if (fact1!=1)
1338  for(size_t i=0; i<newk.size(); i++)
1339  newk[i] *= fact1;
1340 
1341  // define the normalization
1342  int Norm = norm;
1343  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
1344 
1345  if (Norm==1) Pk_0(method_Pk, redshift, store_output, output_root, k_min, k_max, prec, file_par);
1346  else { m_Pk0_EH = 1.; m_Pk0_CAMB = 1.; m_Pk0_MPTbreeze = 1.; m_Pk0_CLASS = 1.; }
1347 
1348  vector<double> Pk(kk.size());
1349 
1350  if (method_Pk=="EisensteinHu") { // NL is not used!!!
1351 
1352  EisensteinHu eh;
1353 
1354  eh.TFmdm_set_cosm(m_Omega_matter, m_Omega_baryon, m_Omega_neutrinos, m_massive_neutrinos, m_Omega_DE, m_hh, redshift, m_scalar_amp, m_scalar_pivot, m_n_spec);
1355 
1356  for (size_t i=0; i<kk.size(); i++) {
1357  Pk[i] = m_Pk0_EH*eh.Pk(newk[i])*fact2;
1358  }
1359  }
1360 
1361  else if (method_Pk=="CAMB" || method_Pk=="MGCAMB" || method_Pk=="MPTbreeze-v1" || method_Pk=="CLASS" || (method_Pk=="CAMB_wrapper" && kk.size()==1)) {
1362 
1363  vector<double> _kk, _pk;
1364  Table_PkCodes(method_Pk, NL, _kk, _pk, redshift, store_output, output_root, k_max, file_par);
1365 
1366  for(size_t i=0; i<_kk.size(); i++) {
1367  _kk[i] = pow(10., _kk[i]);
1368  _pk[i] = pow(10., _pk[i]);
1369  }
1370 
1371  glob::FuncGrid interp_Pk(_kk, _pk, "Spline");
1372 
1373  Pk = interp_Pk.eval_func(newk);
1374 
1375  double Pk0 = 1.;
1376  if (method_Pk=="CAMB") Pk0 = m_Pk0_CAMB;
1377  else if (method_Pk=="MGCAMB") Pk0 = m_Pk0_CAMB;
1378  else if (method_Pk=="MPTbreeze-v1") Pk0 = m_Pk0_MPTbreeze;
1379  else if (method_Pk=="CLASS") Pk0 = m_Pk0_CLASS;
1380 
1381  for (size_t i=0; i<kk.size(); i++)
1382  Pk[i] *= Pk0*fact2;
1383  }
1384 
1385  else if (method_Pk=="CAMB_wrapper" && kk.size()>1) {
1386 
1387  double fact3 = (m_unit || unit1) ? m_hh : 1.;
1388 
1389  Pk = wrapper::camb::Pk_CAMB (NL, redshift, Min(kk)*fact3, Max(kk)*fact3, (int)(kk.size()), m_Omega_baryon*m_hh*m_hh, m_Omega_CDM*m_hh*m_hh, m_Omega_neutrinos*m_hh*m_hh, m_massless_neutrinos, m_massive_neutrinos, m_Omega_k, m_hh*100., m_n_spec, m_scalar_amp, m_scalar_pivot, m_w0, m_wa, m_tau);
1390 
1391  for (size_t i=0; i<kk.size(); i++)
1392  Pk[i] *= m_Pk0_CAMB*fact2;
1393 
1394  }
1395 
1396  else { ErrorCBL("method_Pk is wrong!", "Pk", "PkXi.cpp"); vector<double> vv; return vv; }
1397 
1398  return Pk;
1399 }
1400 
1401 
1402 // =====================================================================================
1403 
1404 
1405 std::vector<std::vector<double>> cbl::cosmology::Cosmology::Pk_matter (const std::vector<double> kk, const std::string method_Pk, const bool NL, const std::vector<double> redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec, const std::string file_par, const bool unit1)
1406 {
1407 
1408  double fact1 = (m_unit || unit1) ? 1. : 1./m_hh;
1409  double fact2 = pow(fact1, 3);
1410 
1411  vector<double> newk = kk;
1412  if (fact1!=1)
1413  for(size_t i=0; i<newk.size(); i++)
1414  newk[i] *= fact1;
1415 
1416  // define the normalization
1417  int Norm = norm;
1418  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
1419 
1420  vector<vector<double>> Pk(redshift.size(), vector<double>(kk.size()));
1421  double Pk0 = 1.;
1422 
1423  if (method_Pk=="EisensteinHu") { // NL is not used!!!
1424 
1425  for (size_t zz=0; zz<redshift.size(); zz++) {
1426 
1427  if (Norm==1) Pk_0(method_Pk, redshift[zz], store_output, output_root, k_min, k_max, prec, file_par);
1428  else { m_Pk0_EH = 1.;}
1429 
1430  EisensteinHu eh;
1431 
1432  eh.TFmdm_set_cosm(m_Omega_matter, m_Omega_baryon, m_Omega_neutrinos, m_massive_neutrinos, m_Omega_DE, m_hh, redshift[zz], m_scalar_amp, m_scalar_pivot, m_n_spec);
1433 
1434  for (size_t i=0; i<kk.size(); i++) {
1435  Pk[zz][i] = m_Pk0_EH*eh.Pk(newk[i])*fact2;
1436  }
1437  }
1438  }
1439 
1440  else if (method_Pk=="CAMB" || method_Pk=="MGCAMB" || method_Pk=="MPTbreeze-v1" || method_Pk=="CLASS") {
1441 
1442  vector<vector<double>> _kk, _pk;
1443  Table_PkCodes(method_Pk, NL, _kk, _pk, redshift, store_output, output_root, k_max, file_par);
1444 
1445  for (size_t zz=0; zz<redshift.size(); zz++) {
1446 
1447  for(size_t i=0; i<_kk[zz].size(); i++) {
1448  _kk[zz][i] = pow(10., _kk[zz][i]);
1449  _pk[zz][i] = pow(10., _pk[zz][i]);
1450  }
1451 
1452  glob::FuncGrid interp_Pk(_kk[zz], _pk[zz], "Spline");
1453 
1454  if (Norm==1) {
1455 
1456  double sigma8;
1457 
1458  if (NL==true)
1459  sigma8 = sigma8_Pk("CAMB", 0., store_output, "test", false, k_min, k_max, prec);
1460 
1461  else {
1462  const double RR = 8.;
1463  auto func_sigma = [&] (double _k){
1464  return pow(TopHat_WF(_k*RR)*_k, 2)*interp_Pk(_k);
1465  };
1466  sigma8 = sqrt(1./(2.*pow(par::pi, 2))*wrapper::gsl::GSL_integrate_qag (func_sigma, k_min, k_max, 1.e-5))/DN(redshift[zz], 0.);
1467  }
1468 
1469  Pk0 = pow(m_sigma8/sigma8,2);
1470  }
1471 
1472  Pk[zz] = interp_Pk.eval_func(newk);
1473 
1474  for (size_t i=0; i<kk.size(); i++)
1475  Pk[zz][i] *= Pk0*fact2;
1476 
1477  }
1478  }
1479 
1480  else { ErrorCBL("method_Pk is wrong!", "Pk", "PkXi.cpp"); vector<vector<double>> vv; return vv; }
1481 
1482  return Pk;
1483 }
1484 
1485 
1486 // =====================================================================================
1487 
1489 
1490 double cbl::glob::func_xi_EH_GSL (double kk, void *params)
1491 {
1492  struct glob::STR_xi_EH *pp = (struct glob::STR_xi_EH *) params;
1493 
1494  EisensteinHu eh;
1495 
1496  eh.TFmdm_set_cosm(pp->Omega_matter, pp->Omega_baryon, pp->Omega_neutrinos, pp->massive_neutrinos, pp->Omega_DE, pp->hh, pp->redshift, pp->scalar_amp, pp->scalar_pivot, pp->n_spec);
1497 
1498  double Int = eh.Pk(kk)*sin(kk*pp->rr)*kk/pp->rr;
1499 
1500  return Int*exp(-kk*kk*pp->aa*pp->aa); // eq. 24 of Anderson et al. 2012
1501 }
1502 
1504 
1505 
1506 // =====================================================================================
1507 
1508 
1509 double cbl::cosmology::Cosmology::xi_matter (const double rr, const std::string method_Pk, const bool NL, const double redshift, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double aa, const bool GSL, const double prec, const std::string file_par)
1510 {
1511  bool gsl = GSL;
1512  if (gsl==false && method_Pk=="EisensteinHu") {
1513  //WarningMsgCBL("EisensteinHu method only works with GSL integration", "xi_matter", "PkXi.cpp");
1514  gsl = true;
1515  }
1516 
1517  int Norm = norm;
1518  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
1519 
1520  if (method_Pk=="MPTbreeze-v1") Norm = 0; // check!!!
1521 
1522  double Int = -1.;
1523  double fact = (gsl) ? 1./(2.*pow(par::pi,2)) : 1.;
1524 
1525  if (gsl) {
1526  int limit_size = 1000;
1527  gsl_integration_workspace *ww = gsl_integration_workspace_alloc(limit_size);
1528  gsl_function Func;
1529 
1530  double error = -1.;
1531 
1532  if (method_Pk=="EisensteinHu") {
1533 
1534  if (m_sigma8<0) ErrorCBL("sigma8<0!", "xi_matter", "PkXi.cpp");
1535  if (NL==1) WarningMsgCBL("the correlation function by Eisenstein&Hu is linear (see xi_matter of PkXi.cpp)!", "xi_matter", "PkXi.cpp");
1536 
1537  glob::STR_xi_EH str;
1538  str.Omega_matter = m_Omega_matter;
1539  str.Omega_baryon = m_Omega_baryon;
1540  str.Omega_neutrinos = m_Omega_neutrinos;
1541  str.massless_neutrinos = m_massless_neutrinos;
1542  str.massive_neutrinos = m_massive_neutrinos;
1543  str.Omega_DE = m_Omega_DE;
1544  str.Omega_radiation = m_Omega_radiation;
1545  str.hh = m_hh;
1546  str.scalar_amp = m_scalar_amp;
1547  str.scalar_pivot = m_scalar_pivot;
1548  str.n_spec = m_n_spec;
1549  str.w0 = m_w0;
1550  str.wa = m_wa;
1551  str.fNL = m_fNL;
1552  str.type_NG = m_type_NG;
1553  str.tau = m_tau;
1554  str.model = m_model;
1555  str.unit = m_unit;
1556  str.rr = rr;
1557  str.aa = aa;
1558  str.redshift = redshift;
1559  str.method_Pk = method_Pk;
1560 
1561  Func.function = &glob::func_xi_EH_GSL;
1562  Func.params = &str;
1563  gsl_integration_qag(&Func, k_min, k_max, 0., prec, limit_size, 6, ww, &Int, &error);
1564  }
1565 
1566  else if (method_Pk=="CAMB" || method_Pk=="MGCAMB" || method_Pk=="MPTbreeze-v1" || method_Pk=="CLASS") {
1567  vector<double> lgkk, lgPk;
1568  Table_PkCodes(method_Pk, NL, lgkk, lgPk, redshift, store_output, output_root, k_max, file_par);
1569 
1570  glob::STR_xi str;
1571  str.rr = rr;
1572  str.aa = aa;
1573  str.lgkk = lgkk;
1574  str.lgPk = lgPk;
1575 
1576  Func.function = &glob::func_xi_GSL;
1577  Func.params = &str;
1578  gsl_integration_qag(&Func, k_min, k_max, 0., prec, limit_size, 5, ww, &Int, &error);
1579  }
1580 
1581  else ErrorCBL("method_Pk is wrong!", "xi_matter", "PkXi.cpp");
1582 
1583  gsl_integration_workspace_free(ww);
1584  }
1585 
1586 
1587  else { // using FFTLOG
1588  if (method_Pk=="CAMB" || method_Pk=="MGCAMB" || method_Pk=="MPTbreeze-v1" || method_Pk=="CLASS") {
1589  vector<double> r, xi;
1590  Table_XiCodes(method_Pk, NL, r, xi, redshift, store_output, output_root, k_max, file_par);
1591  Int = interpolated(rr, r, xi, "Spline");
1592  }
1593 
1594  else ErrorCBL("method_Pk is wrong!", "xi_matter", "PkXi.cpp");
1595 
1596  }
1597 
1598 
1599  if (Norm==1) Pk_0(method_Pk, redshift, store_output, output_root, k_min, k_max, prec, file_par);
1600 
1601  double PP0 = -1.;
1602  if (method_Pk=="EisensteinHu") PP0 = m_Pk0_EH;
1603  if (method_Pk=="CAMB" || method_Pk=="MGCAMB") PP0 = m_Pk0_CAMB;
1604  if (method_Pk=="MPTbreeze-v1") PP0 = m_Pk0_MPTbreeze;
1605  if (method_Pk=="CLASS") PP0 = m_Pk0_CLASS;
1606 
1607  return PP0*fact*Int;
1608 }
1609 
1610 
1611 // =====================================================================================
1612 
1613 
1614 double cbl::cosmology::Cosmology::wp_DM (const double rp, const std::string method_Pk, const bool NL, const double redshift, const double pimax, const bool store_output, const std::string output_root, const int norm, const double r_min, const double r_max, const double k_min, const double k_max, const double aa, const bool GSL, const double prec, const std::string file_par)
1615 {
1616  int Norm = norm;
1617  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
1618 
1619  if (method_Pk=="MPTbreeze-v1") Norm = 0; // check!!!
1620  if (method_Pk=="MPTbreeze-v1" && NL==0) ErrorCBL("MPTbreeze is non-linear!", "wp_DM", "PkXi.cpp");
1621 
1622  // check if the table with lg(r)-lg(xi) already exists
1623 
1624  string mDir = (GSL==0) ? "fftlog" : "GSL";
1625 
1626  cbl::Path path;
1627  string dir_grid = path.DirCosmo()+"/Cosmology/Tables/"+mDir+"/"+method_Pk+"/h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6)+"/";
1628 
1629  //dir_grid += (output_root=="test") ? "/" : "_"+output_root;
1630 
1631  string file_table = (NL) ? dir_grid+"xiDM_NL.dat" : dir_grid+"xiDM_Lin.dat";
1632  ifstream fin;
1633  fin.open(file_table.c_str());
1634 
1635  double RR, XI;
1636  vector<double> rr, Xi;
1637 
1638  if (fin) // read the table
1639  while (fin >> RR >> XI) {
1640  rr.push_back(RR);
1641  Xi.push_back(XI);
1642  }
1643 
1644  else { // create the table
1645  coutCBL <<"I'm writing the file: "<<file_table<<endl;
1646 
1647  string MK = "mkdir -p "+dir_grid; if (system(MK.c_str())) {}
1648  ofstream fout(file_table.c_str()); checkIO(fout, file_table);
1649 
1650  int step = 200;
1651  vector<double> rad = logarithmic_bin_vector(step, r_min, r_max);
1652 
1653  int index = 0;
1654  while (index<int(rad.size())) {
1655  RR = rad[index];
1656  XI = xi_matter(rad[index], method_Pk, NL, redshift, store_output, output_root, 0, k_min, k_max, aa, GSL, prec, file_par);
1657  fout << RR << " " << XI << endl;
1658  coutCBL << "xi(" << RR << ") = " << XI << endl;
1659  rr.push_back(RR);
1660  Xi.push_back(XI);
1661  index ++;
1662  }
1663  fout.clear(); fout.close(); coutCBL <<"I wrote the file: "<<file_table<<endl;
1664  }
1665 
1666  fin.clear(); fin.close();
1667 
1668  double rmax_integral = sqrt(rp*rp+pimax*pimax);
1669  double Int = wp(rp, rr, Xi, rmax_integral);
1670 
1671  if (Norm==1) Pk_0(method_Pk, redshift, store_output, output_root, k_min, k_max, prec, file_par);
1672 
1673  double PP0 = -1.;
1674  if (method_Pk=="EisensteinHu") PP0 = m_Pk0_EH;
1675  if (method_Pk=="CAMB" || method_Pk=="MGCAMB") PP0 = m_Pk0_CAMB;
1676  if (method_Pk=="MPTbreeze-v1") PP0 = m_Pk0_MPTbreeze;
1677  if (method_Pk=="CLASS") PP0 = m_Pk0_CLASS;
1678 
1679  return PP0*Int;
1680 }
1681 
1682 
1683 // =====================================================================================
1684 
1685 
1686 double cbl::cosmology::Cosmology::sigmaR_DM (const double RR, const int corrType, const std::string method_Pk, const double redshift, const double pimax, const bool store_output, const std::string output_root, const bool NL, const int norm, const double r_min, const double r_max, const double k_min, const double k_max, const double aa, const bool GSL, const double prec, const std::string file_par)
1687 {
1688  // check if the table with lg(r)-lg(xi) already exists
1689 
1690  string mDir = "GSL";
1691 
1692  cbl::Path path;
1693  string dir_grid = path.DirCosmo()+"/Cosmology/Tables/"+mDir+"/"+method_Pk+"/h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6)+"/";
1694 
1695  //dir_grid += (output_root=="test") ? "/" : "_"+output_root;
1696 
1697  string file_table = (corrType==1) ? dir_grid+"xi_matter.dat" : dir_grid+"wp_DM.dat";
1698  ifstream fin;
1699  fin.open(file_table.c_str());
1700 
1701  double RRR, XI;
1702  vector<double> rr, Xi;
1703 
1704  if (fin) { // read the table
1705  while (fin>>RRR>>XI) {
1706  rr.push_back(RRR);
1707  Xi.push_back(XI);
1708  }
1709  }
1710 
1711  else { // create the table
1712  coutCBL <<"I'm writing the file: "<<file_table<<endl;
1713 
1714  string MK = "mkdir -p "+dir_grid; if (system(MK.c_str())) {}
1715  ofstream fout(file_table.c_str()); checkIO(fout, file_table);
1716 
1717  int step = 200;
1718  vector<double> rad = logarithmic_bin_vector(step, r_min, r_max);
1719 
1720  int index = 0;
1721  while (index<int(rad.size())) {
1722  RRR = rad[index];
1723  XI = (corrType==1) ? xi_matter(rad[index], method_Pk, NL, redshift, store_output, output_root, norm, k_min, k_max, aa, GSL, prec, file_par) : wp_DM(rad[index], method_Pk, NL, redshift, pimax, store_output, output_root, norm, r_min, r_max, k_min, k_max, aa, GSL, prec, file_par);
1724  fout <<RRR<<" "<<XI<<endl;
1725  coutCBL <<"xi("<<RRR<<") = "<<XI<<endl;
1726  rr.push_back(RRR);
1727  Xi.push_back(XI);
1728  index ++;
1729  }
1730  fout.clear(); fout.close(); coutCBL <<"I wrote the file: "<<file_table<<endl;
1731  }
1732 
1733  fin.clear(); fin.close();
1734 
1735  return sigmaR(RR, corrType, rr, Xi);
1736 }
1737 
1738 
1739 // =====================================================================================
1740 
1742 
1743 double cbl::glob::func_sigma2M_EH_GSL (double kk, void *params)
1744 {
1745  struct glob::STR_sigma2M_EH *pp = (struct glob::STR_sigma2M_EH *) params;
1746 
1747  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);
1748 
1749  double RHO = cosm.rho_m(0., true);
1750  double rr = Radius(pp->mass,RHO);
1751 
1752  EisensteinHu eh;
1753 
1754  eh.TFmdm_set_cosm(pp->Omega_matter, pp->Omega_baryon, pp->Omega_neutrinos, pp->massive_neutrinos, pp->Omega_DE, pp->hh, pp->redshift, pp->scalar_amp, pp->scalar_pivot, pp->n_spec);
1755 
1756  return eh.Pk(kk)*pow(TopHat_WF(kk*rr)*kk, 2);
1757 }
1758 
1760 
1761 // =====================================================================================
1762 
1763 
1764 double cbl::cosmology::Cosmology::sigma8_Pk (const std::string method_Pk, const double redshift, const bool store_output, const std::string output_root, const bool NL, const double k_min, const double k_max, const double prec, const std::string file_par) const
1765 {
1766  if (NL) WarningMsgCBL("sigma8 is defined for the linear P(k)!", "sigma8_Pk", "PkXi.cpp");
1767 
1768  if (m_sigma8>0) return m_sigma8*DN(redshift);
1769 
1770  else {
1771 
1772  const double RR = 8.; // sigma_8 = sigma(8Mpc/h)
1773  const double RHO = rho_m(0., true);
1774  const double MM = Mass(RR, RHO);
1775  double Int = -1., error = -1.;
1776 
1777  const int limit_size = 1000;
1778  gsl_integration_workspace *ww = gsl_integration_workspace_alloc(limit_size);
1779  gsl_function Func;
1780 
1781  if (method_Pk=="EisensteinHu") {
1782 
1783  EisensteinHu eh;
1784 
1785  eh.TFmdm_set_cosm(m_Omega_matter, m_Omega_baryon, m_Omega_neutrinos, m_massive_neutrinos, m_Omega_DE, m_hh, redshift, m_scalar_amp, m_scalar_pivot, m_n_spec);
1786 
1787  auto func = [&] (double kk)
1788  {
1789  return pow(TopHat_WF(kk*RR)*kk, 2)*eh.Pk(kk);
1790  };
1791 
1792  Int = wrapper::gsl::GSL_integrate_qag(func, k_min, k_max, prec, limit_size);
1793 
1794  }
1795 
1796  else if (method_Pk=="CAMB" || method_Pk=="MPTbreeze-v1" || method_Pk=="CLASS") {
1797  vector<double> lgkk, lgPk;
1798  Table_PkCodes(method_Pk, NL, lgkk, lgPk, redshift, store_output, output_root, k_max, file_par);
1799 
1800  glob::STR_SSM str;
1801  str.unit = true;
1802  str.hh = m_hh;
1803  str.n_spec = m_n_spec;
1804  str.mass = MM;
1805  str.rho = RHO;
1806  str.lgkk = lgkk;
1807  str.lgPk = lgPk;
1808 
1809  Func.function = &glob::func_SSM_GSL;
1810  Func.params = &str;
1811  gsl_integration_qag(&Func, k_min, k_max, 0., prec, limit_size, 6, ww, &Int, &error);
1812  }
1813 
1814  else ErrorCBL("method_Pk is wrong!", "sigma8_Pk", "PkXi.cpp");
1815 
1816  gsl_integration_workspace_free(ww);
1817 
1818  return sqrt(1./(2.*par::pi*par::pi)*Int);
1819  }
1820 }
1821 
1822 
1823 // =====================================================================================
1824 
1825 
1826 double cbl::cosmology::Cosmology::Sn_PT (const int nn, const double RR, const std::string method_SS, const bool store_output, const std::string output_root, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file) const
1827 {
1828  if (3>nn || nn>5) ErrorCBL("nn = " + conv(nn, par::fINT), "Sn_PT", "PkXi.cpp");
1829 
1830  double redshift = 0.; // (the hierarchical moments predicted by the PT do not depend on the redshift)
1831 
1832  double gamma1 = 1., gamma2 = -1., gamma3 = -1., d2S = -1., d3S = -1., Sn = 1.;
1833 
1834  double RHO = rho_m(0., true);
1835  double MASS = Mass(RR,RHO);
1836  double SSS = sigma2M(MASS, method_SS, redshift, store_output, output_root, interpType, k_max, input_file, is_parameter_file);
1837 
1838  gamma1 = RR/SSS*dnsigma2R(1, RR, method_SS, redshift, store_output, output_root, interpType, k_max, input_file, is_parameter_file);
1839 
1840  if (nn>3) {
1841  d2S = dnsigma2R(2, RR, method_SS, redshift, store_output, output_root, interpType, k_max, input_file, is_parameter_file);
1842  gamma2 = gamma1+pow(RR, 2)/SSS*d2S;
1843  }
1844 
1845  if (nn>4) {
1846  d3S = dnsigma2R(3, RR, method_SS, redshift, store_output, output_root, interpType, k_max, input_file, is_parameter_file);
1847  gamma3 = gamma2+pow(RR, 2)/SSS*(2.*d2S+RR*d3S);
1848  }
1849 
1850  if (nn==3) Sn = 34./7.+gamma1;
1851  if (nn==4) Sn = 60712./1323.+62./3.*gamma1+7./3.*pow(gamma1, 2)+2./3.*gamma2;
1852  if (nn==5) Sn = 200575880./305613.+1847200./3969.*gamma1+6940./63.*pow(gamma1, 2)+235./27.*pow(gamma1, 3)+1490./63.*gamma2+50./9.*gamma1*gamma2+10./27.*gamma3;
1853 
1854  return Sn;
1855 }
1856 
1857 
1858 // =====================================================================================
1859 
1860 
1861 double cbl::cosmology::Cosmology::Sigman_PT (const int nn, const double RR, const std::string method_SS, const bool store_output, const std::string output_root, const std::string interpType, const double k_max, const std::string input_file, const bool is_parameter_file) const
1862 {
1863  if (3>nn || nn>5) ErrorCBL("nn = " + conv(nn, par::fINT), "Sigman_PT", "PkXi.cpp");
1864 
1865  double redshift = 0.; // (the hierarchical moments predicted by the PT do not depend on the redshift)
1866 
1867  double RHO = rho_m(redshift, true);
1868  double MASS = Mass(RR, RHO);
1869  double SSS = sigma2M(MASS, method_SS, redshift, store_output, output_root, interpType, k_max, input_file, is_parameter_file);
1870 
1871  double gamma1 = RR/SSS*dnsigma2R(1, RR, method_SS, redshift, store_output, output_root, interpType, k_max, input_file, is_parameter_file);
1872 
1873  double Sn = -1.;
1874 
1875  if (nn==3) Sn = 36./7.+3./2.*(gamma1+1.);
1876  if (nn==4) Sn = 2540./49.+33.*(gamma1+1.)+21./4.*pow(gamma1+1.,2);
1877  if (nn==5) Sn = 793.+794.*(gamma1+1.)+265.*pow(gamma1+1.,2)+29.4*pow(gamma1+1.,3);
1878 
1879  return Sn;
1880 }
1881 
1882 
1883 // =====================================================================================
1884 
1885 
1886 double cbl::cosmology::Cosmology::k_star (const std::string method_Pk, const double redshift, const bool store_output, const std::string output_root, const double k_max, const std::string file_par) const
1887 {
1888  if (method_Pk=="EisensteinHu") ErrorCBL("", "k_star", "PkXi.cpp", glob::ExitCode::_workInProgress_);
1889 
1890  vector<double> lgkk, lgPk;
1891  bool do_nonlinear = 0;
1892  Table_PkCodes(method_Pk, do_nonlinear, lgkk, lgPk, redshift, store_output, output_root, k_max, file_par);
1893 
1894  classfunc::func_kstar func (m_hh, m_unit, lgkk, lgPk);
1895 
1896  function<double(double)> ff = bind(&classfunc::func_kstar::operator(), func, std::placeholders::_1);
1897  double Int1 = wrapper::gsl::GSL_integrate_qag(ff, 0., 1., 1.e-4);
1898  double Int2 = wrapper::gsl::GSL_integrate_qag(ff, 1., 1.e30, 1.e-4);
1899 
1900  double Int = Int1+Int2;
1901 
1902  return pow(1./(3.*par::pi*par::pi)*Int,-0.5);
1903 }
1904 
1905 
1906 // =====================================================================================
1907 
1908 
1909 void cbl::cosmology::Cosmology::get_xi (std::vector<double> &rr, std::vector<double> &Xi, const std::string method_Pk, const double redshift, const bool store_output, const std::string output_root, const bool xiType, const double k_star, const bool xiNL, const int norm, const double r_min, const double r_max, const double k_min, const double k_max, const double aa, const bool GSL, const double prec, const std::string file_par)
1910 {
1911  int Norm = norm;
1912  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
1913 
1914  if (method_Pk=="MPTbreeze-v1") Norm = 0; // check!!!
1915 
1916  bool XiNL = xiNL;
1917  if (XiNL && method_Pk=="EisensteinHu")
1918  { WarningMsgCBL("The P(k) of EisensteinHu is linear! --> XiNL = 0", "get_xi", "PkXi.cpp"); XiNL = 0; }
1919 
1920 
1921  // ----- compute the real space DM xi(r) -----
1922 
1923  // check if the table with lg(r)-lg(xi) already exists
1924 
1925  string mDir = "GSL";
1926  string nDir = (xiType==0) ? method_Pk : "CWmodel";
1927 
1928  cbl::Path path;
1929  string dir_cosmo = path.DirCosmo();
1930 
1931  string dir_grid = dir_cosmo+"/Cosmology/Tables/"+mDir+"/"+nDir+"/h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6)+"/";
1932 
1933  //dir_grid += (output_root=="test") ? "/" : "_"+output_root;
1934 
1935  string file_table = (XiNL) ? dir_grid+"xi_matter.dat": dir_grid+"xi_matter_lin.dat";
1936 
1937  //coutCBL <<endl<<"file with tabulated values of xi(r): "<<file_table<<endl<<endl;
1938 
1939  ifstream fin (file_table.c_str());
1940 
1941  double RR, XI;
1942 
1943  if (fin) { // read the table
1944  while (fin >>RR>>XI) {
1945  rr.push_back(RR);
1946  Xi.push_back(XI);
1947  }
1948  fin.clear(); fin.close();
1949  }
1950 
1951  else { // create the table
1952 
1953  string MK = "mkdir -p "+dir_grid; if (system(MK.c_str())) {}
1954  ofstream fout(file_table.c_str()); checkIO(fout, file_table);
1955 
1956  int step = 1000;
1957  vector<double> rad = linear_bin_vector(step, r_min, r_max);
1958 
1959  int index = 0;
1960  while (index<int(rad.size())) {
1961  RR = rad[index];
1962  XI = (xiType==0) ? xi_matter(rad[index], method_Pk, XiNL, redshift, store_output, output_root, Norm, k_min, k_max, aa, GSL, prec, file_par) :
1963  xi_star(rad[index], redshift, store_output, output_root, k_star, k_max, k_max, prec, file_par);
1964 
1965  fout <<RR<<" "<<XI<<endl;
1966  coutCBL <<"xi("<<RR<<") = "<<XI<<endl;
1967 
1968  rr.push_back(RR);
1969  Xi.push_back(XI);
1970  index ++;
1971  }
1972  fout.clear(); fout.close(); coutCBL <<"I wrote the file: "<<file_table<<endl;
1973  }
1974 
1975 }
1976 
1977 
1978 // =====================================================================================
1979 
1980 
1981 void cbl::cosmology::Cosmology::get_barred_xi (std::vector<double> rr, std::vector<double> Xi, std::vector<double> &Xi_, std::vector<double> &Xi__, const std::string method_Pk, const double redshift, const bool xiType, const double k_star, const bool xiNL, const int norm, const double r_min, const double r_max, const double k_min, const double k_max, const double aa, const double prec, const std::string file_par) const
1982 {
1983  (void)k_star; (void)k_min; (void)k_max; (void)aa; (void)prec; (void)file_par;
1984 
1985  int Norm = norm;
1986  if (Norm==-1) Norm = (m_sigma8>0) ? 1 : 0;
1987 
1988  if (method_Pk=="MPTbreeze-v1") Norm = 0; // check!!!
1989 
1990  bool XiNL = xiNL;
1991  if (XiNL && method_Pk=="EisensteinHu")
1992  { WarningMsgCBL("The P(k) of EisensteinHu is linear! --> XiNL = 0", "get_barred_xi", "PkXi.cpp"); XiNL = 0;}
1993 
1994 
1995  // ----- compute the barred functions: xi_(r) and xi__(r) -----
1996 
1997  // check if the table with lg(r)-lg(xi) already exists
1998 
1999  string mDir = "GSL";
2000  string nDir = (xiType==0) ? method_Pk : "CWmodel";
2001 
2002  cbl::Path path;
2003  string dir_cosmo = path.DirCosmo();
2004 
2005  string dir_grid = dir_cosmo+"/Cosmology/Tables/"+mDir+"/"+nDir+"/h"+conv(m_hh, par::fDP6)+"_OmB"+conv(m_Omega_baryon, par::fDP6)+"_OmCDM"+conv(m_Omega_CDM, par::fDP6)+"_OmL"+conv(m_Omega_DE, par::fDP6)+"_OmN"+conv(m_Omega_neutrinos, par::fDP6)+"_Z"+conv(redshift, par::fDP6)+"_scalar_amp"+conv(m_scalar_amp, par::ee3)+"_scalar_pivot"+conv(m_scalar_pivot, par::fDP6)+"_n"+conv(m_n_spec, par::fDP6)+"_w0"+conv(m_w0, par::fDP6)+"_wa"+conv(m_wa, par::fDP6);
2006 
2007  string file_table = (XiNL) ? dir_grid+"xi_matter.dat": dir_grid+"xi_matter_lin.dat";
2008 
2009  ifstream fin(file_table.c_str());
2010 
2011  string file_tableb = (XiNL) ? dir_grid+"xibarred_DM.dat": dir_grid+"xibarred_DM_lin.dat";
2012 
2013  ifstream finb (file_tableb.c_str());
2014 
2015  double RR, XI_, XI__;
2016 
2017  bool DO = (rr.size()==0) ? 1 : 0;
2018 
2019  if (finb) { // read the table
2020  while (finb >>RR>>XI_>>XI__) {
2021  if (DO) rr.push_back(RR);
2022  Xi_.push_back(XI_);
2023  Xi__.push_back(XI__);
2024  }
2025  finb.clear(); finb.close();
2026  }
2027 
2028  else { // create the table
2029  coutCBL <<"I'm writing the file: "<<file_tableb<<endl;
2030 
2031  string MK = "mkdir -p "+dir_grid; if (system(MK.c_str())) {}
2032  ofstream foutb(file_tableb.c_str()); checkIO(foutb, file_tableb);
2033 
2034  vector<double> rad;
2035 
2036  if (DO) {
2037  int step = 1000;
2038  rad = linear_bin_vector(step, r_min, r_max);
2039  }
2040  else rad = rr;
2041 
2042  int index = 0;
2043  while (index<int(rad.size())) {
2044  XI_ = barred_xi_direct(rad[index], rr, Xi, 0, -1, -1);
2045  XI__ = barred_xi__direct(rad[index], rr, Xi, 0, -1, -1);
2046 
2047  foutb <<rad[index]<<" "<<XI_<<" "<<XI__<<endl;
2048  coutCBL <<"r = "<<rad[index]<<" --> xi_ = "<<XI_<<", xi__ = "<<XI__<<endl;
2049 
2050  Xi_.push_back(XI_);
2051  Xi__.push_back(XI__);
2052 
2053  index ++;
2054  }
2055  foutb.clear(); foutb.close(); coutCBL <<"I wrote the file: "<<file_tableb<<endl;
2056  }
2057 }
2058 
2059 
2060 // =====================================================================================
2061 
2062 
2063 vector<double> cbl::cosmology::Cosmology::Pk_matter_NoWiggles_gaussian (const vector<double> kk, const vector<double> PkLin, const vector<double> PkApprox, const double lambda, const string kind)
2064 {
2065  vector<double> PkNW(kk.size());
2066  vector<double> OF(kk.size());
2067 
2068  for (size_t i=0; i<kk.size(); i++)
2069  OF[i] = PkLin[i]/PkApprox[i];
2070 
2071  // Get Pk normalization
2072  glob::FuncGrid interp_OF(kk, OF, "Spline");
2073 
2074  // Smooth the oscillatory part
2075  if (kind=="gaussian_3d")
2076  {
2077  double norm = sqrt(2/par::pi)/lambda*log(10.);
2078  for (size_t i=0; i<kk.size(); i++)
2079  {
2080  auto integrand = [&] (const double log_q) {
2081  double qq = pow(10, log_q);
2082  double x = qq*kk[i]/(lambda*lambda);
2083  double fact = -(qq*qq+kk[i]*kk[i])/(2*lambda*lambda)+gsl_sf_lnsinh(x)-log(kk[i]*qq)+3*log(qq)+log(interp_OF(qq));
2084  return gsl_sf_exp(fact);
2085  };
2086 
2087  PkNW[i] = wrapper::gsl::GSL_integrate_cquad(integrand, -5, 3)*PkApprox[i]*norm;
2088  }
2089  }
2090  else if (kind=="gaussian_1d")
2091  {
2092  double norm = 1./(sqrt(2*par::pi)*lambda);
2093  for (size_t i=0; i<kk.size(); i++)
2094  {
2095  double log_k = log10(kk[i]);
2096 
2097  double log_qmin = log_k-4*lambda;
2098  double log_qmax = log_k+4*lambda;
2099 
2100  auto integrand = [&] (const double log_q) {
2101  return interp_OF(pow(10., log_q))*exp(-pow(log_k-log_q, 2)/(2*lambda*lambda));
2102  };
2103 
2104  PkNW[i] = wrapper::gsl::GSL_integrate_cquad(integrand, log_qmin, log_qmax)*norm*PkApprox[i];
2105  }
2106  }
2107  else
2108  ErrorCBL("wrong name for gaussian wiggles smoothing!", "Pk_matter_NoWiggles_gaussian", "PkXi.cpp", glob::ExitCode::_error_);
2109 
2110  return PkNW;
2111 }
2112 
2113 
2114 // =====================================================================================
2115 
2116 
2117 vector<double> cbl::cosmology::Cosmology::Pk_matter_NoWiggles_bspline (const vector<double> kk, const vector<double> PkLin, const vector<double> PkApprox, const int order, const int nknots)
2118 {
2119  vector<double> log_kk(kk.size());
2120  vector<double> PkNW(kk.size());
2121  vector<double> OF(kk.size());
2122 
2123  for (size_t i=0; i<kk.size(); i++) {
2124  log_kk[i] = log10(kk[i]);
2125  OF[i] = PkLin[i]/PkApprox[i];
2126  }
2127 
2128  // Create b-spline
2129  glob::FuncGrid_Bspline bspline(log_kk, OF, nknots, order);
2130  vector<double> interp = bspline.eval_func(log_kk);
2131  for (size_t i=0; i<kk.size(); i++)
2132  PkNW[i] = interp[i]*PkApprox[i];
2133 
2134  return PkNW;
2135 }
2136 
2137 
2138 // =====================================================================================
2139 
2140 
2141 vector<double> cbl::cosmology::Cosmology::Pk_matter_NoWiggles (const string method, const vector<double> kk, const double redshift, const string linear_method, const int order, const int nknots, const double lambda, const bool store_output, const std::string output_root, const bool norm, const double prec)
2142 {
2143  vector<double> PkNW;
2144  if (method == "EisensteinHu") {
2145  PkNW = Pk_matter(kk, "EisensteinHu", false, redshift, store_output, output_root, norm, 1.e-4, 100., prec);
2146  }
2147  else if (method == "bspline" or method=="gaussian_1d" or method=="gaussian_3d"){
2148  vector<double> PkLin = Pk_matter(kk, linear_method, false, redshift, store_output, output_root, norm, 1.e-4, 100., prec);
2149  vector<double> PkApprox = Pk_matter(kk, "EisensteinHu", false, redshift, store_output, output_root, norm, 1.e-4, 100., prec);
2150 
2151  if (method == "bspline") {
2152  PkNW = Pk_matter_NoWiggles_bspline(kk, PkLin, PkApprox, order, nknots);
2153  }
2154  else if (method == "gaussian_3d" or method == "gaussian_1d") {
2155  PkNW = Pk_matter_NoWiggles_gaussian(kk, PkLin, PkApprox, lambda, method);
2156  }
2157 
2158  // Get Pk normalization
2159  glob::FuncGrid interp_camb(kk, PkLin, "Spline");
2160  double integral = interp_camb.integrate_qag(1.e-5, 1.e3);
2161 
2162  glob::FuncGrid interp_nw(kk, PkNW, "Spline");
2163  double interp_integral = interp_nw.integrate_qag(1.e-5, 1.e3);
2164 
2165  for (size_t i=0; i<kk.size(); i++)
2166  PkNW[i] = PkNW[i]*integral/interp_integral;
2167  }
2168  else
2169  ErrorCBL("wrong name for pk no wiggles!", "Pk_NoWiggles", "PkXi.cpp", glob::ExitCode::_error_);
2170 
2171  return PkNW;
2172 }
2173 
2174 
2175 // =====================================================================================
2176 
2177 
2178 vector<double> cbl::cosmology::Cosmology::Pk_matter_Linear (const string method, const vector<double> kk, const double redshift, const bool store_output, const std::string output_root, const bool norm, const double prec)
2179 {
2180  vector<double> pk;
2181  if (method=="CAMB" or method=="CLASS")
2182  pk = Pk_matter(kk, method, false, redshift, store_output, output_root, norm, 1.e-4, 100., prec);
2183  else
2184  ErrorCBL("wrong name for linear!", "Pk_matter_Linear", "PkXi.cpp", glob::ExitCode::_error_);
2185 
2186  return pk;
2187 }
2188 
2189 
2190 // =====================================================================================
2191 
2192 
2193 vector<double> cbl::cosmology::Cosmology::Pk_matter_DeWiggled (const string linear_method, const string nowiggles_method, const vector<double> kk, const double redshift, const double sigma_NL, const int order, const int nknots, const double lambda, const bool store_output, const std::string output_root, const bool norm, const double prec)
2194 {
2195  vector<double> PkLin = Pk_matter_Linear(linear_method, kk, redshift, store_output, output_root, norm, prec);
2196  vector<double> PkNW = Pk_matter_NoWiggles(nowiggles_method, kk, redshift, linear_method, order, nknots, lambda, store_output, output_root, norm, prec);
2197  vector<double> PkDW(kk.size(), 0);
2198 
2199  for (size_t i=0; i<kk.size(); i++)
2200  PkDW[i] = PkNW[i]+exp(0.5*pow(kk[i]*sigma_NL, 2))*(PkLin[i]-PkNW[i]);
2201 
2202  return PkDW;
2203 }
2204 
2205 
2206 // =====================================================================================
2207 
2208 
2209 double cbl::cosmology::Cosmology::xi_matter_DeWiggle (const double rr, const double redshift, const double sigma_NL, const bool store_output, const std::string output_root, const bool norm, const double k_min, const double k_max, const double aa, const double prec)
2210 {
2211  bool NL = false;
2212 
2213  string author1 = "CAMB";
2214  string author2 = "EisensteinHu";
2215 
2216  vector<double> kk, PkCamb, PkM, PkEH;
2217  Table_PkCodes(author1, NL, kk, PkCamb, redshift, store_output, output_root, k_max);
2218 
2219  for (size_t i = 0; i<kk.size(); i++)
2220  kk[i] = pow(10,kk[i]);
2221 
2222  PkCamb = Pk_matter(kk, author1, NL, redshift, store_output, output_root, norm, k_min, k_max, prec);
2223  PkEH = Pk_matter(kk, author2, NL, redshift, store_output, output_root, norm, k_min, k_max, prec);
2224 
2225  for (size_t i = 0; i<kk.size(); i++)
2226  PkM.push_back(PkEH[i]*(1+(PkCamb[i]/PkEH[i]-1)*exp(-0.5*pow(kk[i]*sigma_NL, 2))));
2227 
2228  return xi_from_Pk(rr, kk, PkM, k_min, k_max, aa, prec);
2229 }
2230 
2231 
2232 // =====================================================================================
2233 
2234 
2235 std::vector<std::vector<double> > cbl::cosmology::Cosmology::XiMonopole_covariance (const int nbins, const double rMin, const double rMax, const double nn, const double Volume, const std::vector<double> kk, const std::vector<double> Pk0, const int IntegrationMethod)
2236 {
2237  int nbins_k = kk.size();
2238  vector<double> r = linear_bin_vector(nbins,rMin,rMax);
2239  vector<vector<double>> covariance(nbins,vector<double>(nbins,0));
2240  double dr=r[1]-r[0];
2241 
2242  vector<vector<double>> sigma2 = cbl::sigma2_k(nn, Volume,kk,{Pk0},{0});
2243  vector<vector<double>> jr(nbins,vector<double>(nbins_k,0));
2244 
2245  for (size_t j=0; j<r.size(); j++) {
2246  for (size_t i=0; i<kk.size(); i++) {
2247  jr[j][i] = jl_distance_average(kk[i], 0, r[j]-dr, r[j]+dr);
2248  }
2249  }
2250 
2251  if (IntegrationMethod==0) //Perform trapezoid integration
2252  {
2253  vector<double> integrand(kk.size(),0);
2254 
2255  for (int i=0; i<nbins; i++) {
2256  for (int j=i; j<nbins; j++) {
2257  for (int ii=0; ii<nbins_k; ii++) {
2258  integrand[ii] =kk[ii]*kk[ii]*sigma2[0][ii]*jr[i][ii]*jr[j][ii];
2259  }
2260 
2261  double Int = trapezoid_integration(kk,integrand);
2262 
2263  Int = Int/(2.*par::pi*par::pi);
2264  covariance[i][j] = Int;
2265  covariance[j][i] = Int;
2266  }
2267  }
2268  }
2269  else if (IntegrationMethod==1) //Perform integration with GSL
2270  {
2271  glob::STR_covariance_XiMultipoles_integrand params;
2272  int limit_size = 1000;
2273 
2274  gsl_function Func;
2275  Func.function = &covariance_XiMultipoles_integrand;
2276  Func.params = &params;
2277 
2278  double k_min=1.e-4;
2279  double k_max = 1.e0;
2280  double prec = 1.e-2;
2281 
2282  glob::FuncGrid s2(kk,sigma2[0], "Spline");
2283  params.s2 = &s2;
2284 
2285  for (int i=0; i<nbins; i++) {
2286  glob::FuncGrid jl1r1(kk,jr[i], "Spline");
2287  for (int j=i; j<nbins; j++) {
2288  glob::FuncGrid jl2r2(kk,jr[j], "Spline");
2289  params.jl1r1 = &jl1r1;
2290  params.jl2r2 = &jl2r2;
2291 
2292  double Int = wrapper::gsl::GSL_integrate_qag(Func,k_min,k_max, prec,limit_size,6);
2293  Int = Int/(2.*par::pi*par::pi);
2294  covariance[i][j] = Int;
2295  covariance[j][i] = Int;
2296  jl2r2.free();
2297  }
2298  jl1r1.free();
2299  }
2300  s2.free();
2301  }
2302 
2303  return covariance;
2304 
2305 }
2306 
2307 
2308 // =====================================================================================
2309 
2310 
2311 std::vector<std::vector<double> > cbl::cosmology::Cosmology::XiMultipoles_covariance (const int nbins, const double rMin, const double rMax, const double nn, const double Volume, const std::vector<double> kk, const std::vector<double> Pk0, const std::vector<double> Pk2, const std::vector<double> Pk4, const int IntegrationMethod)
2312 {
2313  int n_leg = 3;
2314  int nbins_k = kk.size();
2315  vector<double> r = linear_bin_vector(nbins, rMin, rMax);
2316  vector<vector<double>> covariance(n_leg*nbins, vector<double>(n_leg*nbins, 0));
2317 
2318  vector<vector<double>> sigma2 = sigma2_k(nn, Volume, kk, {Pk0, Pk2, Pk4}, {0,2,4});
2319  double dr = r[1]-r[0];
2320 
2321  vector<vector<vector<double> >> jr(n_leg, vector<vector<double>>(nbins, vector<double>(nbins_k, 0)));
2322 
2323  for (int ll=0; ll<n_leg; ll++)
2324  for (size_t j=0; j<r.size(); j++)
2325  for (size_t i=0; i<kk.size(); i++)
2326  jr[ll][j][i] = jl_distance_average(kk[i], ll*2, r[j]-dr, r[j]+dr);
2327 
2328  if (IntegrationMethod==0) { // perform trapezoid integration
2329 
2330  vector<double> integrand(kk.size(), 0);
2331 
2332  for (int l1=0; l1<n_leg; l1++) {
2333  for (int l2=l1; l2<n_leg; l2++) {
2334  int index = l2+n_leg*l1;
2335  for (int i=0; i<nbins; i++) {
2336  for (int j=i; j<nbins; j++) {
2337  for (int ii=0; ii<nbins_k; ii++) {
2338  integrand[ii] = kk[ii]*kk[ii]*sigma2[index][ii]*jr[l1][i][ii]*jr[l2][j][ii];
2339  }
2340 
2341  double Int = trapezoid_integration(kk, integrand);
2342 
2343  Int = Int/(2.*par::pi*par::pi);
2344  covariance[i+nbins*l1][j+nbins*l2] = Int;
2345  covariance[j+nbins*l1][i+nbins*l2] = Int;
2346  covariance[i+nbins*l2][j+nbins*l1] = Int;
2347  covariance[j+nbins*l2][i+nbins*l1] = Int;
2348 
2349  }
2350  }
2351  }
2352  }
2353  }
2354 
2355  else if (IntegrationMethod==1) // perform integration with GSL
2356  {
2357  glob::STR_covariance_XiMultipoles_integrand params;
2358  int limit_size = 1000;
2359 
2360  gsl_function Func;
2361  Func.function = &covariance_XiMultipoles_integrand;
2362  Func.params = &params;
2363 
2364  double k_min = 1.e-4;
2365  double k_max = 1.e0;
2366  double prec = 1.e-2;
2367 
2368  for (int l1=0; l1<n_leg; l1++) {
2369  for (int l2=l1; l2<n_leg; l2++) {
2370  int index = l2+n_leg*l1;
2371  glob::FuncGrid s2(kk, sigma2[index], "Spline");
2372  params.s2 = &s2;
2373 
2374  for (int i=0; i<nbins; i++) {
2375  glob::FuncGrid jl1r1(kk, jr[l1][i], "Spline");
2376  for (int j=i; j<nbins; j++) {
2377  glob::FuncGrid jl2r2(kk, jr[l2][j], "Spline");
2378  params.jl1r1 = &jl1r1;
2379  params.jl2r2 = &jl2r2;
2380 
2381  double Int = wrapper::gsl::GSL_integrate_qag(Func, k_min, k_max, prec, limit_size, 6);
2382  Int = Int/(2.*par::pi*par::pi);
2383  covariance[i+nbins*l1][j+nbins*l2] = Int;
2384  covariance[j+nbins*l1][i+nbins*l2] = Int;
2385  covariance[i+nbins*l2][j+nbins*l1] = Int;
2386  covariance[j+nbins*l2][i+nbins*l1] = Int;
2387  jl2r2.free();
2388  }
2389  jl1r1.free();
2390  }
2391  s2.free();
2392  }
2393  }
2394  }
2395 
2396  return covariance;
2397 }
2398 
2399 
2400 // =====================================================================================
2401 
2402 
2403 std::vector<std::vector<double> > cbl::cosmology::Cosmology::XiMultipoles (const int nbins, const double rMin, const double rMax, const std::vector<double> kk, const std::vector<double> Pk0, const std::vector<double> Pk2, const std::vector<double> Pk4, const int IntegrationMethod)
2404 {
2405  int nbins_k = kk.size();
2406  vector<double> r = linear_bin_vector(nbins,rMin,rMax);
2407  double f0 = 1./(2.*par::pi*par::pi), f2 = -1./(2.*par::pi*par::pi), f4 =1./(2.*par::pi*par::pi);
2408 
2409  vector<double> xi0(nbins,0), xi2(nbins,0), xi4(nbins,0);
2410 
2411  vector<vector<double> > j0_vec(nbins,vector<double>(nbins_k,0)), j2_vec(nbins,vector<double>(nbins_k,0)), j4_vec(nbins,vector<double>(nbins_k,0));
2412 
2413  for (size_t i=0; i<r.size(); i++) {
2414  for (size_t j=0; j<kk.size(); j++) {
2415  double xx = kk[j]*r[i];
2416  j0_vec[i][j] = j0(xx);
2417  j2_vec[i][j] = j2(xx);
2418  j4_vec[i][j] = j4(xx);
2419  }
2420  }
2421 
2422  if (IntegrationMethod==0) { // perform trapezoid integration
2423 
2424  for (int i=0; i<nbins; i++) {
2425 
2426  vector<double> i0(kk.size(),0), i2(kk.size(),0), i4(kk.size(),0);
2427 
2428  for (int j=0; j<nbins_k; j++) {
2429  double xx = kk[j]*r[i];
2430  i0[j] =kk[j]*kk[j]*Pk0[j]*j0(xx);
2431  i2[j] =kk[j]*kk[j]*Pk2[j]*j2(xx);
2432  i4[j] =kk[j]*kk[j]*Pk4[j]*j4(xx);
2433  }
2434 
2435  xi0[i] = trapezoid_integration(kk,i0)*f0;
2436  xi2[i] = trapezoid_integration(kk,i2)*f2;
2437  xi4[i] = trapezoid_integration(kk,i4)*f4;
2438  }
2439 
2440  }
2441 
2442  else if (IntegrationMethod==1) { // perform integration with GSL
2443 
2444  glob::STR_XiMultipoles_integrand params;
2445  int limit_size = 1000;
2446 
2447  gsl_function Func;
2448  Func.function = &XiMultipoles_integrand;
2449  Func.params = &params;
2450 
2451  double k_min = 1.e-4;
2452  double k_max = 1.e0;
2453  double prec = 1.e-3;
2454 
2455  glob::FuncGrid Pk0_interp(kk, Pk0, "Spline");
2456  glob::FuncGrid Pk2_interp(kk, Pk2, "Spline");
2457  glob::FuncGrid Pk4_interp(kk, Pk4, "Spline");
2458 
2459  for (int i=0; i<nbins; i++) {
2460  params.r = r[i];
2461 
2462  params.Pkl = &Pk0_interp;
2463  params.l = 0;
2464  params.k_cut = 0.7;
2465  params.cut_pow = 2.;
2466 
2467  xi0[i] = wrapper::gsl::GSL_integrate_qag(Func,k_min,k_max, prec,limit_size,6)*f0;
2468 
2469  params.Pkl = &Pk2_interp;
2470  params.l = 2;
2471  params.k_cut = 0.58;
2472  params.cut_pow = 4.;
2473  xi2[i] = wrapper::gsl::GSL_integrate_qag(Func, k_min, k_max, prec, limit_size, 6)*f2;
2474 
2475  params.Pkl = &Pk4_interp;
2476  params.l = 4;
2477  params.k_cut = 0.6;
2478  params.cut_pow = 2.;
2479  xi4[i] = wrapper::gsl::GSL_integrate_qag(Func, k_min, k_max, prec, limit_size, 6)*f4;
2480  }
2481 
2482  Pk0_interp.free();
2483  Pk2_interp.free();
2484  Pk4_interp.free();
2485 
2486  }
2487 
2488  return {xi0, xi2, xi4};
2489 }
2490 
2491 
2492 // =====================================================================================
2493 
2494 
2495 double cbl::cosmology::Cosmology::wtheta_DM (const double theta, const std::vector<double> zz, const std::vector<double> phiz, const std::string interpolationType, const CoordinateUnits coordUnits, const bool GSL, const std::string method_Pk, const bool NL, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec, const std::string file_par)
2496 {
2497  if (NL)
2498  ErrorCBL("non linearities in angular correlation function not yet implemented!", "wtheta_DM", "PkXi.cpp", glob::ExitCode::_workInProgress_);
2499 
2500  double theta_rad = converted_angle (theta, coordUnits, CoordinateUnits::_radians_);
2501 
2502  vector<double> kk = logarithmic_bin_vector(200, k_min, k_max);
2503  vector<double> Pk;
2504 
2505  vector<double> dc;
2506  vector<double> _zz = linear_bin_vector(1000, 0., 2*Max(zz));
2507 
2508  for(size_t i=0; i<_zz.size(); i++)
2509  dc.push_back(this->D_C(_zz[i]));
2510  glob::FuncGrid DC_interp(_zz, dc, interpolationType);
2511 
2512  glob::Distribution phi(glob::DistributionType::_Interpolated_, zz, phiz, 0, interpolationType);
2513 
2514  double zmin, zmax, zmean;
2515  zmin = Min(zz);
2516  zmax = Max(zz);
2517  zmean = phi.mean();
2518 
2519  Pk = this->Pk_matter(kk, method_Pk, NL, zmean, store_output, output_root, norm, k_min, k_max, prec, file_par);
2520 
2521  vector<double> r, xi;
2522  wrapper::fftlog::transform_FFTlog(r, xi, 1, kk, Pk);
2523  glob::FuncGrid xi_interp(r, xi, interpolationType);
2524 
2525  if (GSL) {
2526  auto integrand = [&] (double z1)
2527  {
2528  double r1 = DC_interp(z1);
2529 
2530  auto integrand_z2 = [&] (double z2) {
2531  double r2 = DC_interp(z2);
2532  double ss = sqrt(pow(r1,2)+pow(r2,2)-2*r1*r2*cos(theta_rad));
2533  return xi_interp(ss)*phi(z2);
2534  };
2535 
2536  return wrapper::gsl::GSL_integrate_qag(integrand_z2, zmin, zmax)*phi(z1);
2537  };
2538  return wrapper::gsl::GSL_integrate_qag(integrand, zmin, zmax);
2539  }
2540  else{
2541  auto integrand = [&] (vector<double> zz)
2542  {
2543  double r1 = DC_interp(zz[0]);
2544  double r2 = DC_interp(zz[1]);
2545  double ss = sqrt(pow(r1,2)+pow(r2,2)-2*r1*r2*cos(theta_rad));
2546  return xi_interp(ss)*phi(zz[0])*phi(zz[1]);
2547  };
2548 
2549  wrapper::cuba::CUBAwrapper integrator(integrand, 2);
2550 
2551  return integrator.IntegrateCuhre( {{zmin, zmax}, {zmin, zmax}});
2552  }
2553 }
2554 
2555 
2556 // =====================================================================================
2557 
2558 
2559 double cbl::cosmology::Cosmology::wtheta_DM (const double theta, const std::vector<double> kk, const std::vector<double> Pk, const std::vector<double> zz, const std::vector<double> nz, const std::vector<double> phiz, const std::string interpolationType, const CoordinateUnits coordUnits, const bool GSL, const double redshift_Pk)
2560 {
2561  const double theta_rad = converted_angle (theta, coordUnits, CoordinateUnits::_radians_);
2562 
2563  const double zmin = Min(zz);
2564  const double zmax = Max(zz);
2565 
2566 
2567  // set the distribution function
2568 
2569  vector<double> nphi(zz.size(), 0);
2570 
2571  for(size_t i=0; i<nphi.size(); i++)
2572  nphi[i] = (phiz.size()!=0) ? phiz[i]*nz[i] : nz[i];
2573 
2574  glob::Distribution phi(glob::DistributionType::_Interpolated_, zz, nphi, 0, interpolationType);
2575  //const double zmean = phi.mean();
2576 
2577 
2578  // set the distirbution integrand and normalization
2579 
2580  auto normalization_integrand = [&] (const double redshift)
2581  {
2582  return phi(redshift)*this->dV_dZdOmega(redshift, 1);
2583  };
2584 
2585  double normalization = wrapper::gsl::GSL_integrate_qag(normalization_integrand, zmin, zmax);
2586 
2587 
2588  // compute the 2PCF model
2589 
2590  vector<double> r, xi;
2591  wrapper::fftlog::transform_FFTlog(r, xi, 1, kk, Pk);
2592  glob::FuncGrid xi_interp(r, xi, interpolationType);
2593  double integral;
2594 
2595  if (GSL) {
2596  auto integrand = [&] (double z1)
2597  {
2598  double r1 = this->D_C(z1);
2599 
2600  auto integrand_z2 = [&] (double z2) {
2601  double DD = this->DN((zz[0]+zz[1])*0.5, redshift_Pk);
2602  double r2 = this->D_C(z2);
2603  double ss = sqrt(pow(r1,2)+pow(r2,2)-2*r1*r2*cos(theta_rad));
2604  return DD*DD*xi_interp(ss)*normalization_integrand(z2);
2605  };
2606 
2607  return wrapper::gsl::GSL_integrate_qag(integrand_z2, zmin, zmax)*normalization_integrand(z1);
2608  };
2609 
2610  integral = wrapper::gsl::GSL_integrate_qag(integrand, zmin, zmax);
2611 
2612  }
2613 
2614  else {
2615  auto integrand = [&] (vector<double> zz)
2616  {
2617  double DD = this->DN((zz[0]+zz[1])*0.5, redshift_Pk);
2618  double r1 = this->D_C(zz[0]);
2619  double r2 = this->D_C(zz[1]);
2620  double ss = sqrt(pow(r1,2)+pow(r2,2)-2*r1*r2*cos(theta_rad));
2621  return DD*DD*xi_interp(ss)*normalization_integrand(zz[0])*normalization_integrand(zz[1]);
2622  };
2623 
2624  wrapper::cuba::CUBAwrapper integrator(integrand, 2);
2625 
2626  integral = integrator.IntegrateCuhre( {{zmin, zmax}, {zmin, zmax}});
2627  }
2628 
2629  return integral/normalization;
2630 }
2631 
2632 
2633 // =====================================================================================
2634 
2635 
2636 std::vector<double> cbl::cosmology::Cosmology::C_l_DM (const int lmax, const std::vector<double> zz, const std::vector<double> phiz, const std::string interpolationMethod, const std::string method_Pk, const bool store_output, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec, const std::string file_par)
2637 {
2638  const double zmin = Min(zz);
2639  const double zmax = Max(zz);
2640 
2641  vector<double> kk = logarithmic_bin_vector(200, k_min, k_max);
2642 
2643  vector<double> Pk = this->Pk_matter(kk, method_Pk, false, 0., store_output, output_root, norm, k_min, k_max, prec, file_par);
2644  glob::FuncGrid Pk_interp(kk, Pk, interpolationMethod);
2645 
2646  auto integrand_sbao = [&] (const double kk)
2647  {
2648  return Pk_interp(kk);
2649  };
2650 
2651  double sbao = sqrt(4.*par::pi*wrapper::gsl::GSL_integrate_qag(integrand_sbao, 1.e-4, 1)/3./pow(2*par::pi,3));
2652 
2653  glob::FuncGrid phi(zz, phiz, interpolationMethod);
2654 
2655  vector<double> C_l;
2656 
2657  ErrorCBL("Check the normalization of D(z)", "C_l_DM", "PkXi.cpp", glob::ExitCode::_workInProgress_);
2658 
2659  for (int l=0; l<lmax+1; l++) {
2660  double integral;
2661 
2662  if (l<60) {
2663  auto integrand = [&] ( const double kk)
2664  {
2665  auto integrand_z = [&] (const double zz)
2666  {
2667  return DN(zz)*jl(kk*D_C(zz), l);
2668  };
2669 
2670  double integral_z = wrapper::gsl::GSL_integrate_qag(integrand_z, zmin, zmax);
2671  return kk*kk*Pk_interp(kk)*pow(integral_z, 2)*exp(-kk*kk*sbao*sbao);
2672  };
2673 
2674  integral = 2*wrapper::gsl::GSL_integrate_qag(integrand, 1.e-4, 10)/par::pi;
2675  }
2676  else {
2677 
2678  auto integrand = [&] (const double zz)
2679  {
2680  double dc = D_C(zz);
2681  double kk = (l+0.5)/dc;
2682  return pow(DN(zz)*phi(zz), 2)*Pk_interp(kk)*HH(zz)/(par::cc*dc*dc);
2683  };
2684 
2685  integral = wrapper::gsl::GSL_integrate_qag(integrand, zmin, zmax);
2686  }
2687  C_l.push_back(integral);
2688  }
2689 
2690  return C_l;
2691 }
2692 
2693 /*
2694  std::vector<double> cbl::cosmology::Cosmology::C_l_DM (const int lmax, const std::vector<double> zz, const std::vector<double> phiz, const std::string interpolationMethod, const std::string method_Pk, const std::string output_root, const int norm, const double k_min, const double k_max, const double prec, const std::string file_par)
2695  {
2696  const double zmin = Min(zz);
2697  const double zmax = Max(zz);
2698 
2699  vector<double> kk = logarithmic_bin_vector(200, k_min, k_max);
2700  vector<double> Pk;
2701  for(size_t i=0; i<kk.size(); i++)
2702  Pk.push_back( this->Pk_matter(kk[i], method_Pk, false, 0., output_root, norm, k_min, k_max, prec, file_par));
2703 
2704  glob::FuncGrid Pk_interp(kk, Pk, interpolationMethod);
2705  glob::FuncGrid phi(zz, phiz, interpolationMethod);
2706 
2707  vector<double> C_l;
2708 
2709  for (int l=0; l<lmax+1; l++) {
2710  auto integrand = [&] ( const double redshift)
2711  {
2712  double dc = D_C(redshift);
2713  double _kk = double(l)/dc;
2714  return pow(phi(redshift), 2)*HH(redshift)/pow(dc, 2)*pow(DN(redshift),2)*Pk_interp(_kk);
2715  };
2716  C_l.push_back(wrapper::gsl::GSL_integrate_qag(integrand, zmin, zmax)/par::cc);
2717  }
2718 
2719  return C_l;
2720  }
2721 */
2722 
2723 
2724 
Cosmology.h
The class Cosmology.
FuncGrid_Bspline.h
Class used to handle functions interpolated using a basis spline http://mathworld....
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::Path::fullpath
std::string fullpath(std::string path, const bool isDir=true)
substitute ~ with the full path
Definition: Path.cpp:45
cbl::Path::DirLoc
std::string DirLoc()
get the local directory
Definition: Path.h:185
cbl::cosmology::Cosmology
The class Cosmology.
Definition: Cosmology.h:277
cbl::cosmology::Cosmology::Table_PkCodes
void Table_PkCodes(const std::string code, const bool NL, std::vector< double > &lgkk, std::vector< double > &lgPk, const double redshift, const bool store_output=true, const std::string output_root="test", const double k_max=100., const std::string file_par=par::defaultString) const
write or read the table where the dark matter power spectrum is stored
Definition: PkXi.cpp:258
cbl::cosmology::Cosmology::Pk_matter_DeWiggled
std::vector< double > Pk_matter_DeWiggled(const std::string linear_method, const std::string nowiggles_method, const std::vector< double > kk, const double redshift, const double sigma_NL, const int order=4, const int nknots=10, const double lambda=0.25, const bool store_output=true, const std::string output_root="test", const bool norm=1, const double prec=1.e-4)
the dark matter power spectrum, de-wiggled (see e.g. Anderson et al 2014)
Definition: PkXi.cpp:2193
cbl::cosmology::Cosmology::Pk_matter_NoWiggles
std::vector< double > Pk_matter_NoWiggles(const std::string method, const std::vector< double > kk, const double redshift, const std::string linear_method="CAMB", const int order=4, const int nknots=10, const double lambda=0.25, const bool store_output=true, const std::string output_root="test", const bool norm=1, const double prec=1.e-4)
the dark matter power spectrum without BAO wiggles
Definition: PkXi.cpp:2141
cbl::cosmology::Cosmology::get_barred_xi
void get_barred_xi(std::vector< double > rr, std::vector< double > Xi, std::vector< double > &Xi_, std::vector< double > &Xi__, const std::string method_Pk, const double redshift, const bool xiType=0, const double k_star=-1., const bool xiNL=0, const int norm=-1, const double r_min=0.1, const double r_max=150., const double k_min=0.001, const double k_max=100., const double aa=0., const double prec=1.e-2, const std::string file_par=par::defaultString) const
get the barred dark matter correlation functions
Definition: PkXi.cpp:1981
cbl::cosmology::Cosmology::wtheta_DM
double wtheta_DM(const double theta, const std::vector< double > zz, const std::vector< double > phiz, const std::string interpolationMethod, const CoordinateUnits coordUnits=CoordinateUnits::_degrees_, const bool GSL=false, const std::string method_Pk="CAMB", const bool NL=false, const bool store_output=true, const std::string output_root="test", const int norm=-1, const double k_min=1.e-4, const double k_max=100, const double prec=1.e-2, const std::string file_par=par::defaultString)
the dark matter angular two-point correlation function
Definition: PkXi.cpp:2495
cbl::cosmology::Cosmology::get_xi
void get_xi(std::vector< double > &rr, std::vector< double > &Xi, const std::string method_Pk, const double redshift, const bool store_output=true, const std::string output_root="test", const bool xiType=0, const double k_star=-1., const bool xiNL=0, const int norm=-1, const double r_min=0.1, const double r_max=150., const double k_min=0.001, const double k_max=100., const double aa=0., const bool GSL=false, const double prec=1.e-2, const std::string file_par=par::defaultString)
get the dark matter two-point correlation function
Definition: PkXi.cpp:1909
cbl::cosmology::Cosmology::XiMonopole_covariance
std::vector< std::vector< double > > XiMonopole_covariance(const int nbins, const double rMin, const double rMax, const double nn, const double Volume, const std::vector< double > kk, const std::vector< double > Pk0, const int IntegrationMethod=1)
the covariance matrix of the first three non-null multipoles of the two-point correlation function
Definition: PkXi.cpp:2235
cbl::cosmology::Cosmology::C_l_DM
std::vector< double > C_l_DM(const int lmax, const std::vector< double > zz, const std::vector< double > phiz, const std::string interpolationMethod, const std::string method_Pk="CAMB", const bool store_output=true, const std::string output_root="test", const int norm=-1, const double k_min=1.e-4, const double k_max=100, const double prec=1.e-2, const std::string file_par=par::defaultString)
the dark matter angular linear power spectrum .
Definition: PkXi.cpp:2636
cbl::cosmology::Cosmology::wp_DM
double wp_DM(const double rp, const std::string method_Pk, const bool NL, const double redshift, const double pimax, const bool store_output=true, const std::string output_root="test", const int norm=-1, const double r_min=1.e-3, const double r_max=350., const double k_min=0.001, const double k_max=100., const double aa=0., const bool GSL=false, const double prec=1.e-2, const std::string file_par=cbl::par::defaultString)
the dark matter projected correlation function
Definition: PkXi.cpp:1614
cbl::cosmology::Cosmology::remove_output_Pk_tables
void remove_output_Pk_tables(const std::string code, const bool NL, const double redshift, const std::string output_root="test") const
remove the output generated by the methods CAMB, MPTbreeze or CLASS
Definition: PkXi.cpp:1030
cbl::cosmology::Cosmology::Table_XiCodes
void Table_XiCodes(const std::string code, const bool NL, std::vector< double > &rr, std::vector< double > &xi, const double redshift, const bool store_output, const std::string output_root, const double k_max, std::string file_par) const
write or read the table where the dark matter two-point correlation function is stored
Definition: PkXi.cpp:1164
cbl::cosmology::Cosmology::As
double As(const double sigma8) const
amplitude of the curvature perturbations
Definition: PkXi.cpp:49
cbl::cosmology::Cosmology::Pk_matter_NoWiggles_gaussian
std::vector< double > Pk_matter_NoWiggles_gaussian(const std::vector< double > kk, const std::vector< double > PkLin, const std::vector< double > PkApprox, const double lambda, const std::string method)
the dark matter power spectrum without BAO wiggles.
Definition: PkXi.cpp:2063
cbl::cosmology::Cosmology::Pk_matter
std::vector< double > Pk_matter(const std::vector< double > kk, const std::string method_Pk, const bool NL, const double redshift, 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, const bool unit1=false)
the dark matter power spectrum
Definition: PkXi.cpp:1331
cbl::cosmology::Cosmology::xi_matter
double xi_matter(const double rr, const std::string method_Pk, const bool NL, const double redshift, 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 aa=0., const bool GSL=false, const double prec=1.e-2, const std::string file_par=par::defaultString)
the dark matter two-point correlation function
Definition: PkXi.cpp:1509
cbl::cosmology::Cosmology::xi_matter_DeWiggle
double xi_matter_DeWiggle(const double rr, const double redshift, const double sigma_NL, const bool store_output=true, const std::string output_root="test", const bool norm=1, const double k_min=0.001, const double k_max=100., const double aa=1., const double prec=1.e-2)
the dark matter two-point correlation function, de-wiggled (see e.g. Anderson et al 2014)
Definition: PkXi.cpp:2209
cbl::cosmology::Cosmology::m_Table_Pk_parameterFile
void m_Table_Pk_parameterFile(const std::string code, const std::string file_par, const bool NL, std::vector< double > &lgkk, std::vector< double > &lgPk, const double redshift, const std::string output_root="test") const
write and read the table where the dark matter power spectrum is stored; it is used when a parameter ...
Definition: PkXi.cpp:1058
cbl::cosmology::Cosmology::Pk_0
void Pk_0(const std::string method_Pk, const double redshift, const bool store_output=true, const std::string output_root="test", const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string file_par=par::defaultString)
normalisation of the power spectrum
Definition: PkXi.cpp:1251
cbl::cosmology::Cosmology::XiMultipoles
std::vector< std::vector< double > > XiMultipoles(const int nbins, const double rMin, const double rMax, const std::vector< double > kk, const std::vector< double > Pk0, const std::vector< double > Pk2, const std::vector< double > Pk4, const int IntegrationMethod=1)
the first three non-null multipoles of the two-point correlation function
Definition: PkXi.cpp:2403
cbl::cosmology::Cosmology::sigmaR_DM
double sigmaR_DM(const double RR, const int corrType, const std::string method_Pk, const double redshift, const double pimax=40, const bool store_output=true, const std::string output_root="test", const bool NL=1, const int norm=-1, const double r_min=1.e-3, const double r_max=350., const double k_min=0.001, const double k_max=100., const double aa=0., const bool GSL=false, const double prec=1.e-2, const std::string file_par=par::defaultString)
the dark matter rms mass fluctuation
Definition: PkXi.cpp:1686
cbl::cosmology::Cosmology::Pk_matter_Linear
std::vector< double > Pk_matter_Linear(const std::string method, const std::vector< double > kk, const double redshift, const bool store_output=true, const std::string output_root="test", const bool norm=1, const double prec=1.e-4)
the dark matter linear power spectrum.
Definition: PkXi.cpp:2178
cbl::cosmology::Cosmology::k_star
double k_star(const std::string method_Pk, const double redshift, const bool store_output=true, const std::string output_root="test", const double k_max=100., const std::string file_par=par::defaultString) const
the k* parameter
Definition: PkXi.cpp:1886
cbl::cosmology::Cosmology::Sn_PT
double Sn_PT(const int nn, const double RR, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true) const
the hierarchical moments Sn
Definition: PkXi.cpp:1826
cbl::cosmology::Cosmology::Pk_output_file
std::string Pk_output_file(const std::string code, const bool NL, const double redshift, const bool run=0, const bool store_output=true, const std::string output_root="test", const double k_max=100., const std::string file_par=par::defaultString)
return the path to the power spectrum output
Definition: PkXi.cpp:71
cbl::cosmology::Cosmology::Sigman_PT
double Sigman_PT(const int nn, const double RR, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const std::string interpType="Linear", const double k_max=100., const std::string input_file=par::defaultString, const bool is_parameter_file=true) const
the deprojected hierarchical moments Σn
Definition: PkXi.cpp:1861
cbl::cosmology::Cosmology::m_Table_Pk_CLASS
void m_Table_Pk_CLASS(const bool NL, std::vector< double > &lgkk, std::vector< double > &lgPk, const double redshift, const bool store_output=true, const std::string output_root="test", const double k_max=100.) const
write and read the table where the dark matter power spectrum computed with CLASS is stored
Definition: PkXi.cpp:740
cbl::cosmology::Cosmology::run_CAMB
void run_CAMB(const bool NL, const double redshift, const std::string output_root=par::defaultString, const std::string output_dir=par::defaultString, const double k_max=100., const std::string file_par=par::defaultString) const
run CAMB [http://camb.info/]
Definition: PkXi.cpp:106
cbl::cosmology::Cosmology::m_Table_Pk_CAMB_MPTbreeze
void m_Table_Pk_CAMB_MPTbreeze(const std::string code, const bool NL, std::vector< double > &lgkk, std::vector< double > &lgPk, const double redshift, const bool store_output=true, const std::string output_root="test", const double k_max=100.) const
write and read the table where the dark matter power spectrum, computed with either CAMB or MPTbreeze...
Definition: PkXi.cpp:327
cbl::cosmology::Cosmology::sigma8_Pk
double sigma8_Pk(const std::string method_Pk, const double redshift, const bool store_output=true, const std::string output_root="test", const bool NL=0, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string file_par=par::defaultString) const
the dark matter rms mass fluctuation within 8 Mpc/h
Definition: PkXi.cpp:1764
cbl::cosmology::Cosmology::Pk_matter_NoWiggles_bspline
std::vector< double > Pk_matter_NoWiggles_bspline(const std::vector< double > kk, const std::vector< double > PkLin, const std::vector< double > PkApprox, const int order, const int nknots)
the dark matter power spectrum without BAO wiggles.
Definition: PkXi.cpp:2117
cbl::cosmology::Cosmology::XiMultipoles_covariance
std::vector< std::vector< double > > XiMultipoles_covariance(const int nbins, const double rMin, const double rMax, const double nn, const double Volume, const std::vector< double > kk, const std::vector< double > Pk0, const std::vector< double > Pk2, const std::vector< double > Pk4, const int IntegrationMethod=1)
the covariance matrix of the first three non-null multipole moments of the two-point correlation func...
Definition: PkXi.cpp:2311
cbl::cosmology::Cosmology::sigma8_interpolated
double sigma8_interpolated(const double redshift) const
σ8
Definition: PkXi.cpp:58
cbl::cosmology::EisensteinHu
The class EisensteinHu.
Definition: EisensteinHu.h:54
cbl::cosmology::EisensteinHu::Pk
double Pk(double kk)
get the power spectrum in unit of h^{-3}Mpc^{-3}
Definition: EisensteinHu.h:407
cbl::cosmology::EisensteinHu::TFmdm_set_cosm
int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm, int degen_hdm, double omega_lambda, double hubble, double redshift, double As=2.56e-9, double k_pivot=0.05, double n_spec=0.96)
set the cosmological parameters
Definition: EisensteinHu.h:231
cbl::glob::Distribution
The class CombinedDistribution.
Definition: Distribution.h:124
cbl::glob::Distribution::mean
double mean()
return the mean value of the distribution
Definition: Distribution.cpp:408
cbl::glob::FuncGrid
The class FuncGrid.
Definition: FuncGrid.h:55
cbl::glob::FuncGrid::free
void free()
free the GSL objects
Definition: FuncGrid.cpp:110
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::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::IntegrateCuhre
double IntegrateCuhre(std::vector< std::vector< double >> integration_limits, const bool parallelize=true)
integrate using the Cuhre routine
Definition: CUBAwrapper.cpp:201
cbl::par::fDP6
static const char fDP6[]
conversion symbol for: double -> std::string
Definition: Constants.h:145
cbl::par::ee3
static const char ee3[]
conversion symbol for: double -> std::string
Definition: Constants.h:160
cbl::par::fINT
static const char fINT[]
conversion symbol for: int -> std::string
Definition: Constants.h:121
cbl::par::defaultString
static const std::string defaultString
default std::string value
Definition: Constants.h:336
cbl::par::pi
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
cbl::par::cc
static const double cc
: the speed of light in vacuum (the value is exact) [km/sec]
Definition: Constants.h:221
cbl::wrapper::camb::Pk_CAMB
std::vector< double > Pk_CAMB(const bool nonlinear, const double redshift, const double kmin, const double kmax, const int npoints, const double ombh2, const double omch2, const double omnuh2, const double massless_nu, const int massive_nu, const double omk, const double H0, const double ns, const double As, const double pivot_scalar=0.05, const double w=-1., const double wa=0., const double tau=2.1e-9, const bool accurate_massive_nu=false, const int neutrino_hierarchy=3, const int dark_energy_model=1, const double cs2_lam=1., const double T_cmb=2.7255, const double helium_fraction=0.24)
Get the matter power spectrum from CAMB.
Definition: CAMB.cpp:42
cbl::wrapper::fftlog::transform_FFTlog
std::vector< double > transform_FFTlog(const std::vector< double > yy, const int dir, const std::vector< double > xx, const std::vector< double > fx, const double mu=0, const double q=0, const double kr=1, const int kropt=0)
wrapper of the FFTlog to compute the discrete Fourier sine or cosine transform of a logarithmically...
Definition: FFTlog.cpp:43
cbl::wrapper::gsl::GSL_integrate_cquad
double GSL_integrate_cquad(gsl_function func, const double a, const double b, const double rel_err=1.e-3, const double abs_err=0, const int nevals=100)
integral, using the gsl cquad method
Definition: GSLwrapper.cpp:159
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::sigmaR
double sigmaR(const double RR, const int corrType, const std::vector< double > rr, const std::vector< double > corr)
the rms mass fluctuation within a given radius
Definition: FuncXi.cpp:227
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::barred_xi_direct
double barred_xi_direct(const double RR, const std::vector< double > rr, const std::vector< double > xi, const double rAPP=0., const double r0=-1., const double gamma=1.)
the barred correlation function
Definition: FuncXi.cpp:331
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::j2
double j2(const double xx)
the l=2 spherical Bessel function
Definition: Func.cpp:2048
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::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::converted_angle
double converted_angle(const double angle, const CoordinateUnits inputUnits=CoordinateUnits::_radians_, const CoordinateUnits outputUnits=CoordinateUnits::_degrees_)
conversion to angle units
Definition: Func.cpp:192
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::XiMultipoles_integrand
double XiMultipoles_integrand(const double kk, void *parameters)
integrand to obtain covariance for the 2PCF multipoles
Definition: FuncMultipoles.cpp:420
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::sigma2_k
std::vector< std::vector< double > > sigma2_k(const double nObjects, const double Volume, const std::vector< double > kk, const std::vector< std::vector< double >> Pk_multipoles, const std::vector< int > orders)
multipole expansion of the per-mode covariance sigma2_k (see Grieb et al. 2016, eq....
Definition: FuncMultipoles.cpp:944
cbl::Max
T Max(const std::vector< T > vect)
maximum element of a std::vector
Definition: Kernel.h:1336
cbl::barred_xi__direct
double barred_xi__direct(const double RR, const std::vector< double > rr, const std::vector< double > xi, const double rAPP=0., const double r0=-1., const double gamma=1.)
the double barred correlation function
Definition: FuncXi.cpp:361
cbl::j0
double j0(const double xx)
the l=0 spherical Bessel function
Definition: Func.cpp:2039
cbl::trapezoid_integration
double trapezoid_integration(const std::vector< double > xx, const std::vector< double > yy)
integral, computed with the trapezoid rule, using ordered data
Definition: Func.cpp:2203
cbl::xi_from_Pk
double xi_from_Pk(const double rr, const std::vector< double > lgkk, const std::vector< double > lgPk, const double k_min=0., const double k_max=100., const double aa=0., const double prec=1.e-2)
the two-point correlation function computed from the Fourier transform of the power spectrum
Definition: FuncXi.cpp:80
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::j4
double j4(const double xx)
the l=4 spherical Bessel function
Definition: Func.cpp:2057
cbl::jl_distance_average
double jl_distance_average(const double kk, const int order, const double r_down, const double r_up)
the distance average for the order l-th spherical Bessel function
Definition: Func.cpp:2109
cbl::CoordinateUnits
CoordinateUnits
the coordinate units
Definition: Kernel.h:562
cbl::WarningMsgCBL
void WarningMsgCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL)
internal CBL warning message
Definition: Kernel.h:747
cbl::covariance_XiMultipoles_integrand
double covariance_XiMultipoles_integrand(const double kk, void *parameters)
integrand to obtain the 2PCF multipoles
Definition: FuncMultipoles.cpp:466
cbl::jl
double jl(const double xx, const int order)
the order l spherical Bessel function
Definition: Func.cpp:2066
cbl::wp
double wp(const double rp, const std::vector< double > rr, const std::vector< double > xi, const double r_max=100.)
the projected two-point correlation function
Definition: FuncXi.cpp:192