CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Modelling_TwoPointCorrelation1D.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2016 by Federico Marulli and Alfonso Veropalumbo *
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 
37 #include "Data1D_extra.h"
39 
40 using namespace std;
41 
42 using namespace cbl;
43 
44 
45 // ============================================================================================
46 
47 
48 cbl::modelling::twopt::Modelling_TwoPointCorrelation1D::Modelling_TwoPointCorrelation1D (const std::shared_ptr<cbl::measure::twopt::TwoPointCorrelation> twop)
49 {
50  m_data = twop->dataset();
51  m_twoPType = twop->twoPType();
52 }
53 
54 
55 // ============================================================================================
56 
57 
59 {
60  m_data = dataset;
61  m_twoPType = twoPType;
62 }
63 
64 // ============================================================================================
65 
66 
67 void cbl::modelling::twopt::Modelling_TwoPointCorrelation1D::set_data_model (const cbl::cosmology::Cosmology cosmology, const double redshift, const std::string method_Pk, const double sigmaNL_perp, const double sigmaNL_par, const bool NL, const double bias, const double pimax, const double r_min, const double r_max, const double k_min, const double k_max, const int step, const std::string output_dir, const std::string output_root, const int norm, const double aa, const bool GSL, const double prec, const std::string file_par, const double Delta, const bool isDelta_critical, const std::vector<double> cluster_redshift, const std::vector<double> cluster_mass_proxy, const std::vector<double> cluster_mass_proxy_error, const std::string model_bias, const std::string meanType, const int seed, const cbl::cosmology::Cosmology cosmology_mass, const std::vector<double> redshift_source)
68 {
69  if (file_par!=par::defaultString)
70  WarningMsgCBL("check the consistency between the parameters of the object cosmology, provided in input, and the ones in the parameter file", "set_data_model", "Modelling_TwoPointCorrelation1D.cpp");
71 
72  m_data_model = make_shared<STR_data_model>(STR_data_model());
73 
74  m_data_model->cosmology = make_shared<cosmology::Cosmology>(cosmology);
75  m_data_model->redshift = redshift;
76  m_data_model->method_Pk = method_Pk;
77  m_data_model->sigmaNL_perp = sigmaNL_perp;
78  m_data_model->sigmaNL_par = sigmaNL_par;
79  m_data_model->sigmaNL = sqrt(0.5*(sigmaNL_par*sigmaNL_par+sigmaNL_perp*sigmaNL_perp));
80  m_data_model->NL = NL;
81  m_data_model->pi_max = pimax;
82  m_data_model->r_min = r_min;
83  m_data_model->r_max = r_max;
84  m_data_model->k_min = k_min;
85  m_data_model->k_max = k_max;
86  m_data_model->step = step;
87  m_data_model->output_dir = output_dir;
88  m_data_model->output_root = output_root;
89  m_data_model->norm = norm;
90  m_data_model->aa = aa;
91  m_data_model->GSL = GSL;
92  m_data_model->prec = prec;
93  m_data_model->file_par = file_par;
94  m_data_model->bias = bias;
95  m_data_model->Delta_input = Delta;
96  m_data_model->isDelta_critical = isDelta_critical;
97  m_data_model->Delta = (isDelta_critical) ? Delta/cosmology.OmegaM(redshift) : Delta;
98  m_data_model->model_bias = model_bias;
99  m_data_model->meanType = meanType;
100  m_data_model->gau_ran = make_shared<random::NormalRandomNumbers>(random::NormalRandomNumbers(0., 1., seed));
101 
102  if (cosmology.sigma8()>0)
103  m_data_model->sigma8_z = cosmology.sigma8(redshift);
104 
105  else {
106  coutCBL << "sigma8 is not set, it will be computed from the power spectrum with " << method_Pk << endl;
107 
108  m_data_model->sigma8_z = cosmology.sigma8_Pk(method_Pk, redshift, true, output_root, NL, k_min, k_max, prec, file_par);
109 
110  coutCBL << "--> sigma8(z=" << redshift << ") = " << m_data_model->sigma8_z << endl << endl;
111  }
112  m_data_model->linear_growth_rate_z = cosmology.linear_growth_rate(redshift, 1.);
113 
114  m_data_model->DVfid = cosmology.D_V(redshift);
115  m_data_model->DAfid = cosmology.D_A(redshift);
116  m_data_model->HHfid = cosmology.HH(redshift);
117 
118  if (cluster_mass_proxy.size()>0) {
119  vector<double> cl_zz(cluster_mass_proxy.size(), redshift), cl_sigma(cluster_mass_proxy.size());
120  if (cluster_redshift.size()==cluster_mass_proxy.size())
121  cl_zz = cluster_redshift;
122  for (size_t i=0; i<cl_sigma.size(); ++i)
123  cl_sigma[i] = sqrt(cosmology.sigma2M(cluster_mass_proxy[i], method_Pk, 0., true, output_root, "Spline", k_max));
124 
125  m_data_model->cluster_mass_proxy = make_shared<data::Data1D_extra>(data::Data1D_extra(cl_zz, cluster_mass_proxy, cluster_mass_proxy_error, {cl_sigma}));
126 
127  m_data_model->cosmology_mass = cosmology_mass;
128  m_data_model->redshift_source = redshift_source;
129  }
130 }
131 
132 
133 // ============================================================================================
134 
135 
136 void cbl::modelling::twopt::Modelling_TwoPointCorrelation1D::set_data_model (const cbl::cosmology::Cosmology cosmology, const double redshift, const std::vector<double> cluster_redshift, const std::vector<double> cluster_mass_proxy, const double redshift_pivot, const double proxy_pivot, double mass_pivot, const double log_base, const std::string method_Pk, const bool NL, const double r_min, const double r_max, const double k_min, const double k_max, const int step, const std::string output_dir, const std::string output_root, const int norm, const double prec, const std::string file_par, const double Delta, const bool isDelta_critical, const std::string model_bias)
137 {
138  if (file_par!=par::defaultString)
139  WarningMsgCBL("check the consistency between the parameters of the object cosmology, provided in input, and the ones in the parameter file", "set_data_model", "Modelling_TwoPointCorrelation1D.cpp");
140 
141  m_data_model = make_shared<STR_data_model>(STR_data_model());
142 
143  m_data_model->cosmology = make_shared<cosmology::Cosmology>(cosmology);
144  m_data_model->redshift = redshift;
145  m_data_model->method_Pk = method_Pk;
146  m_data_model->NL = NL;
147  m_data_model->r_min = r_min;
148  m_data_model->r_max = r_max;
149  m_data_model->k_min = k_min;
150  m_data_model->k_max = k_max;
151  m_data_model->step = step;
152  m_data_model->output_dir = output_dir;
153  m_data_model->output_root = output_root;
154  m_data_model->norm = norm;
155  m_data_model->prec = prec;
156  m_data_model->file_par = file_par;
157  m_data_model->Delta_input = Delta;
158  m_data_model->isDelta_critical = isDelta_critical;
159  m_data_model->Delta = (isDelta_critical) ? Delta/cosmology.OmegaM(redshift) : Delta;
160  m_data_model->model_bias = model_bias;
161 
162  if (cosmology.sigma8()>0)
163  m_data_model->sigma8_z = cosmology.sigma8(redshift);
164  else {
165  coutCBL << "sigma8 is not set, it will be computed from the power spectrum with " << method_Pk << endl;
166  m_data_model->sigma8_z = cosmology.sigma8_Pk(method_Pk, redshift, false, output_root, NL, k_min, k_max, prec, file_par);
167  coutCBL << "--> sigma8(z=" << redshift << ") = " << m_data_model->sigma8_z << endl << endl;
168  }
169 
170  if (cluster_mass_proxy.size() == 0 || cluster_redshift.size() == 0 || cluster_mass_proxy.size() != cluster_redshift.size())
171  ErrorCBL("The mass proxy and redshift vectors must contain at least one object each, and must have the same size.", "set_data_model", "Modelling_TwoPointCorrelation1D.cpp");
172 
173  // Build a dummy dataset for the scaling relation Modelling object.
174  // Only the x values of the dataset (the mass proxy values) are used in the computation
175  std::shared_ptr<cbl::data::Data> dataset = std::make_shared<cbl::data::Data1D>(cbl::data::Data1D(cluster_mass_proxy, cluster_mass_proxy, cluster_mass_proxy));
176 
177  // Build the scaling relation Modelling object
179  m_data_model->scaling_relation = make_shared<modelling::massobsrel::Modelling_MassObservableRelation>(scaling_relation);
180 
181  (m_data_model->scaling_relation)->set_data_model(cosmology, cluster_redshift, redshift_pivot, proxy_pivot, log_base);
182  (m_data_model->scaling_relation)->set_mass_pivot(mass_pivot);
183 }
184 
185 
186 // ============================================================================================
187 
188 
189 void cbl::modelling::twopt::Modelling_TwoPointCorrelation1D::set_data_model (const double z_abs_err, const double proxy_rel_err, const cbl::cosmology::Cosmology cosmology, const double redshift, const std::vector<double> cluster_redshift, const std::vector<double> cluster_mass_proxy, const double redshift_pivot, const double proxy_pivot, double mass_pivot, const double log_base, const std::string method_Pk, const bool NL, const double r_min, const double r_max, const double k_min, const double k_max, const int step, const std::string output_dir, const std::string output_root, const int norm, const double prec, const std::string file_par, const double Delta, const bool isDelta_critical, const std::string model_bias)
190 {
191  set_data_model(cosmology, redshift, cluster_redshift, cluster_mass_proxy, redshift_pivot, proxy_pivot, mass_pivot, log_base, method_Pk, NL, r_min, r_max, k_min, k_max, step, output_dir, output_root, norm, prec, file_par, Delta, isDelta_critical, model_bias);
192 
193  if (z_abs_err < 0 || proxy_rel_err < 0)
194  ErrorCBL("The errors on z and proxy must be > 0.", "set_data_model", "Modelling_TwoPointCorrelation1D.cpp");
195 
196  m_data_model->z_abs_err = z_abs_err;
197  m_data_model->proxy_rel_err = proxy_rel_err;
198 }
199 
200 
201 // ============================================================================================
202 
203 
204 void cbl::modelling::twopt::Modelling_TwoPointCorrelation1D::set_data_HOD (const cosmology::Cosmology cosmology, const double redshift, const string model_MF, const string model_bias, const double Mh_min, const double Mh_max, const double pi_max, const double r_max_int, const double r_min, const double r_max, const double k_min, const double k_max, const int step, const double m_min, const double m_max, const int m_step, const string method_Pk, const bool NL, const string output_root, const double Delta, const double kk, const string interpType, const int norm, const double prec, const string input_file, const bool is_parameter_file, const string model_cM, const string profile, const string halo_def)
205 {
206  m_data_HOD = make_shared<STR_data_HOD>(STR_data_HOD());
207 
208  m_data_HOD->cosmology = make_shared<cosmology::Cosmology>(cosmology);
209  m_data_HOD->redshift = redshift;
210  m_data_HOD->model_MF = model_MF;
211  m_data_HOD->model_bias = model_bias;
212  m_data_HOD->Mh_min = Mh_min;
213  m_data_HOD->Mh_max = Mh_max;
214  m_data_HOD->pi_max = pi_max;
215  m_data_HOD->r_max_int = r_max_int;
216  m_data_HOD->r_min = r_min;
217  m_data_HOD->r_max = r_max;
218 
219  // creation of the mass vector for interpolating pk
220  m_data_HOD->k_min = k_min;
221  m_data_HOD->k_max = k_max;
222  m_data_HOD->step = step;
223  m_data_HOD->kkvec = logarithmic_bin_vector(step, k_min, k_max);
224 
225  // creation of the mass vector for interpolating mass functiom and bias
226  m_data_HOD->m_min = m_min;
227  m_data_HOD->m_max = m_max;
228  m_data_HOD->m_step = m_step;
229  m_data_HOD->massvec = logarithmic_bin_vector(m_step, m_min, m_max);
230 
231  // creation of vectors containing mass function and bias values for interpolation
232  vector<double> mass_function_vec;
233  vector<double> bias_vec;
234  for (int i=0; i<m_step; i++) {
235  mass_function_vec.emplace_back(m_data_HOD->cosmology->mass_function(m_data_HOD->massvec[i], redshift, model_MF, "CAMB"));
236  bias_vec.emplace_back(m_data_HOD->cosmology->bias_halo(m_data_HOD->massvec[i], redshift, model_MF, "CAMB"));
237  }
238  m_data_HOD->mass_function_vec = mass_function_vec;
239  m_data_HOD->bias_vec = bias_vec;
240 
241  // creation of vector containing pk values for interpolation
242  vector<double> pk_vec = m_data_HOD->cosmology->Pk_matter(m_data_HOD->kkvec, method_Pk, NL, redshift, false, output_root, norm, k_min, k_max, prec, input_file);
243 
244  m_data_HOD->pk_vec = pk_vec;
245 
246  // creation of FuncGrid objects for interpolation
247  glob::FuncGrid MF_grid(m_data_HOD->massvec, mass_function_vec, "Spline");
248  m_data_HOD->interpMF = MF_grid;
249  glob::FuncGrid BI_grid(m_data_HOD->massvec, bias_vec, "Spline");
250  m_data_HOD->interpBias = BI_grid;
251  glob::FuncGrid PK_grid(m_data_HOD->kkvec, pk_vec, "Spline");
252  m_data_HOD->interpPk = PK_grid;
253 
254  m_data_HOD->method_Pk = method_Pk;
255  m_data_HOD->NL = NL;
256  m_data_HOD->output_root = output_root;
257  m_data_HOD->Delta = Delta;
258  m_data_HOD->kk = kk;
259  m_data_HOD->interpType = interpType;
260  m_data_HOD->norm = norm;
261  m_data_HOD->prec = prec;
262  m_data_HOD->input_file = input_file;
263  m_data_HOD->is_parameter_file = is_parameter_file;
264  m_data_HOD->model_cM = model_cM;
265  m_data_HOD->profile = profile;
266  m_data_HOD->halo_def = halo_def;
267 }
268 
269 
270 // ============================================================================================
271 
272 
273 void cbl::modelling::twopt::Modelling_TwoPointCorrelation1D::set_data_model_cluster_selection_function (const cosmology::Cosmology cosmology, const cosmology::Cosmology test_cosmology, const double mean_redshift, const string model_MF, const string model_bias, const string selection_function_file, const std::vector<int> selection_function_column, const double z_min, const double z_max, const double Mass_min, const double Mass_max, const string file_par, const double Delta, const bool isDelta_critical, const string method_Pk, const string output_dir, const double k_min, const double k_max, const double prec, const int step, const int mass_step)
274 {
275  if (file_par!=par::defaultString)
276  WarningMsgCBL("check the consistency between the parameters of the object cosmology, provided in input, and the ones in the parameter file", "set_data_model_cluster_selection_function", "Modelling_TwoPointCorrelation1D.cpp");
277 
278  m_data_model = make_shared<STR_data_model>(STR_data_model());
279 
280  m_data_model->cosmology = make_shared<cosmology::Cosmology>(cosmology);
281  m_data_model->test_cosmology = make_shared<cosmology::Cosmology>(test_cosmology);
282  m_data_model->redshift = mean_redshift;
283  m_data_model->method_Pk = method_Pk;
284  m_data_model->k_min = k_min;
285  m_data_model->k_max = k_max;
286  m_data_model->step = step;
287  m_data_model->prec = prec;
288  m_data_model->output_root = "test";
289  m_data_model->Delta = Delta;
290  m_data_model->isDelta_critical = isDelta_critical;
291  m_data_model->model_MF = model_MF;
292  m_data_model->model_bias = model_bias;
293  m_data_model->kk = logarithmic_bin_vector(step, k_min, k_max);
294  m_data_model->mass = logarithmic_bin_vector(mass_step, Mass_min, Mass_max);
295  m_data_model->file_par = par::defaultString;
296  m_data_model->output_dir = output_dir;
297 
298  vector<double> mass_grid, redshift_grid;
299  vector<vector<double>> selection_function;
300  read_matrix(selection_function_file, mass_grid, redshift_grid, selection_function, selection_function_column);
301  m_data_model->interp_SelectionFunction = make_shared<glob::FuncGrid2D>(glob::FuncGrid2D(mass_grid, redshift_grid, selection_function, "Linear"));
302 
303  vector<double> sf_monoz(mass_grid.size());
304  for (size_t i=0; i<mass_grid.size(); i++)
305  sf_monoz[i] = m_data_model->interp_SelectionFunction->operator()(mass_grid[i], mean_redshift);
306  m_data_model->interp_SelectionFunction_cut = make_shared<glob::FuncGrid>(glob::FuncGrid(mass_grid, sf_monoz, "Spline"));
307 
308  m_data_model->z_min = (z_min>par::defaultDouble) ? z_min : Min(redshift_grid);
309  m_data_model->z_max = (z_max>par::defaultDouble) ? z_max : Max(redshift_grid);
310  m_data_model->Mass_min = (Mass_min>par::defaultDouble) ? Mass_min : Min(mass_grid);
311  m_data_model->Mass_max = (Mass_max>par::defaultDouble) ? Mass_max : Max(mass_grid);
312 
313  m_data_model->DVfid = cosmology.D_V(mean_redshift);
314 }
The class Data1D_extra.
#define coutCBL
CBL print message.
Definition: Kernel.h:734
The class Modelling_TwoPointCorrelation1D.
The class Cosmology.
Definition: Cosmology.h:277
double sigma8() const
get the private member Cosmology::m_sigma8
Definition: Cosmology.h:1212
double D_V(const double redshift) const
the average distance at a given redshift, used to rescale the correlation function
Definition: Cosmology.cpp:884
double HH(const double redshift=0.) const
the Hubble function
Definition: Cosmology.cpp:570
double D_A(const double redshift) const
the angular diameter distance at a given redshift
Definition: Cosmology.cpp:848
double sigma2M(const double mass, const std::string method_Pk, const double redshift, 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 bool unit1=false) const
the mass variance,
Definition: Sigma.cpp:230
double linear_growth_rate(const double redshift, const double prec=1.e-4) const
the linear growth rate at a given redshift,
Definition: Cosmology.cpp:662
double OmegaM(const double redshift=0.) const
the matter density at a given redshift
Definition: Cosmology.cpp:579
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
The class Data1D_extra.
Definition: Data1D_extra.h:52
The class Data1D.
Definition: Data1D.h:51
The class FuncGrid2D.
Definition: FuncGrid.h:302
The class FuncGrid.
Definition: FuncGrid.h:55
void set_data_model(const cbl::cosmology::Cosmology cosmology, const double redshift, const std::string method_Pk="CAMB", const double sigmaNL_perp=0., const double sigmaNL_par=0., const bool NL=true, const double bias=1., const double pimax=40., const double r_min=1., const double r_max=350., const double k_min=1.e-4, const double k_max=100., const int step=500, const std::string output_dir=par::defaultString, const std::string output_root="test", const int norm=-1, const double aa=0., const bool GSL=true, const double prec=1.e-3, const std::string file_par=par::defaultString, const double Delta=200., const bool isDelta_critical=true, const std::vector< double > cluster_redshift={}, const std::vector< double > cluster_mass_proxy={}, const std::vector< double > cluster_mass_proxy_error={}, const std::string model_bias="Tinker", const std::string meanType="mean_bias", const int seed=666, const cbl::cosmology::Cosmology cosmology_mass={}, const std::vector< double > redshift_source={})
set the data used to construct generic models of the two-point correlation function
void set_data_HOD(const cbl::cosmology::Cosmology cosmology={}, const double redshift=0., const std::string model_MF="Tinker", const std::string model_bias="Tinker", const double Mh_min=0., const double Mh_max=1.e16, const double pi_max=100., const double r_max_int=100., const double r_min=1.e-3, const double r_max=350., const double k_min=1.e-4, const double k_max=100., const int step=200, const double m_min=1.e7, const double m_max=1.e17, const int m_step=100, const std::string method_Pk="CAMB", const bool NL=true, const std::string output_root="test", const double Delta=200., const double kk=0., const std::string interpType="Linear", const int norm=-1, const double prec=1.e-2, const std::string input_file=par::defaultString, const bool is_parameter_file=true, const std::string model_cM="Duffy", const std::string profile="NFW", const std::string halo_def="vir")
set the data used to construct the HOD model
void set_data_model_cluster_selection_function(const cbl::cosmology::Cosmology cosmology, const cbl::cosmology::Cosmology test_cosmology, const double mean_redshift, const std::string model_MF, const std::string model_bias, const std::string selection_function_file, const std::vector< int > selection_function_column={}, const double z_min=par::defaultDouble, const double z_max=par::defaultDouble, const double Mass_min=par::defaultDouble, const double Mass_max=par::defaultDouble, const std::string file_par=par::defaultString, const double Delta=200, const bool isDelta_critical=false, const std::string method_Pk="CAMB", const std::string output_dir=par::defaultString, const double k_min=1.e-4, const double k_max=100, const double prec=1.e-2, const int step=200, const int mass_step=50)
set the data used to construct the model for clustering analyses using the selection function
Modelling_TwoPointCorrelation1D()=default
default constuctor
The class NormalRandomNumbers.
static const std::string defaultString
default std::string value
Definition: Constants.h:336
static const double defaultDouble
default double value
Definition: Constants.h:348
TwoPType
the two-point correlation function type
double scaling_relation(const double alpha, const double beta, const double gamma, const double div_proxy, const double f_z)
compute the mass - mass proxy scaling relation
double bias(const double Mmin, const double sigmalgM, const double M0, const double M1, const double alpha, const std::shared_ptr< void > inputs)
the mean galaxy bias
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
T Min(const std::vector< T > vect)
minimum element of a std::vector
Definition: Kernel.h:1324
void read_matrix(const std::string file_matrix, std::vector< double > &xx, std::vector< double > &yy, std::vector< std::vector< double >> &matrix, const std::vector< int > col={})
read a matrix from file
Definition: Func.cpp:612
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
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
T Max(const std::vector< T > vect)
maximum element of a std::vector
Definition: Kernel.h:1336
void WarningMsgCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL)
internal CBL warning message
Definition: Kernel.h:747