CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Modelling_DensityProfile.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2021 by Federico Marulli and Giorgio Lesci *
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 "Modelling_DensityProfile.h"
38 #include "Data1D.h"
39 
40 using namespace std;
41 
42 using namespace cbl;
43 
44 
45 // ===========================================================================================
46 
47 cbl::modelling::densityprofile::Modelling_DensityProfile::Modelling_DensityProfile (const std::shared_ptr<cbl::measure::stackprofile::StackedDensityProfile> profile, const std::string profile_author, const bool _2halo, const std::string halo_def, const double Delta)
48 {
49  m_data = profile->dataset();
50  m_profile_author = profile_author;
51  m_2halo = _2halo;
52  m_mass_is_derived = false;
53  m_halo_def = halo_def;
54  m_Delta = Delta;
55 }
56 
57 
58 // ===========================================================================================
59 
60 cbl::modelling::densityprofile::Modelling_DensityProfile::Modelling_DensityProfile (const std::shared_ptr<cbl::data::Data> dataset, const std::string profile_author, const bool _2halo, const std::string halo_def, const double Delta)
61 {
62  m_data = dataset;
63  m_profile_author = profile_author;
64  m_2halo = _2halo;
65  m_mass_is_derived = false;
66  m_halo_def = halo_def;
67  m_Delta = Delta;
68 }
69 
70 
71 // ===========================================================================================
72 
73 void cbl::modelling::densityprofile::Modelling_DensityProfile::set_data_model (const cosmology::Cosmology cosmology, const double redshift, const double contrast, const double logM_base, const double mass_pivot, const std::string bias_author, const std::string method_Pk, std::string interp_type)
74 {
75  m_data_model.cosmology = make_shared<cosmology::Cosmology>(cosmology);
76  m_data_model.cosmology->set_unit(true); // Force cosmological units
77 
78  m_data_model.redshift = redshift;
79  m_data_model.contrast = contrast;
80  m_data_model.logM_base = logM_base;
81  m_data_model.mass_pivot = mass_pivot;
82 
83  m_data_model.bias_author = bias_author;
84  m_data_model.method_Pk = method_Pk;
85  m_data_model.interp_type = interp_type;
86 }
87 
88 
89 // ===========================================================================================
90 
91 void cbl::modelling::densityprofile::Modelling_DensityProfile::set_data_model (const cosmology::Cosmology cosmology, const double redshift, const double mass_proxy, const double redshift_pivot, const double proxy_pivot, const double contrast, const double logM_base, const double mass_pivot, const double Nclusters, const std::string bias_author, const std::string method_Pk, std::string interp_type)
92 {
93  m_data_model.cosmology = make_shared<cosmology::Cosmology>(cosmology);
94  m_data_model.cosmology->set_unit(true); // Force cosmological units
95 
96  m_data_model.redshift = redshift;
97  m_data_model.mass_proxy = mass_proxy;
98  m_data_model.redshift_pivot = redshift_pivot;
99  m_data_model.proxy_pivot = proxy_pivot;
100 
101  m_data_model.contrast = contrast;
102  m_data_model.logM_base = logM_base;
103  m_data_model.mass_pivot = mass_pivot;
104 
105  m_data_model.bias_author = bias_author;
106  m_data_model.method_Pk = method_Pk;
107  m_data_model.interp_type = interp_type;
108 
109  m_mass_is_derived = true;
110 
111  // Build a dummy dataset for the scaling relation Modelling object, useful only to avoid internal errors
112  std::vector<double> dummy_vec = {1.};
113  std::shared_ptr<cbl::data::Data> dataset = std::make_shared<cbl::data::Data1D>(cbl::data::Data1D(dummy_vec, dummy_vec, dummy_vec));
114 
115  // Build the scaling relation object
116  modelling::massobsrel::Modelling_MassObservableRelation scaling_relation (dataset);
117  m_data_model.scaling_relation = make_shared<modelling::massobsrel::Modelling_MassObservableRelation>(scaling_relation);
118 
119  (m_data_model.scaling_relation)->set_data_model(cosmology, {redshift}, redshift_pivot, proxy_pivot, logM_base, {Nclusters});
120 }
121 
122 
123 // ===========================================================================================
124 
125 void cbl::modelling::densityprofile::Modelling_DensityProfile::set_model_DensityProfile_cosmology (const std::vector<cbl::cosmology::CosmologicalParameter> cosmo_param, const std::vector<statistics::PriorDistribution> cosmo_prior, const statistics::PriorDistribution Rt_prior, const statistics::PriorDistribution concentration_prior, const statistics::PriorDistribution logM_prior, const statistics::PriorDistribution f_off_prior, const statistics::PriorDistribution sigma_off_prior, const statistics::PriorDistribution anisotropic_boost_prior, const statistics::PriorDistribution orientation_boost_prior)
126 {
127  m_data_model.Cpar = cosmo_param;
128 
129  const size_t nParams = cosmo_param.size()+7; // The total number of parameters is given by the cosmological ones + 7, since the density profile has 7 parameters (Rt, conc, logM, f_off, sigma_off, AB_fact, OB_fact)
130 
131  vector<statistics::ParameterType> Par_type (nParams, statistics::ParameterType::_Base_);
132  vector<string> Par_string (nParams);
133  std::vector<statistics::PriorDistribution> param_prior (nParams);
134 
135  // Set the names and priors of the cosmological parameters
136  for (size_t i=0; i<cosmo_param.size(); i++){
137  Par_string[i] = CosmologicalParameter_name(cosmo_param[i]);
138  param_prior[i] = cosmo_prior[i];
139  }
140 
141  // Set the names and priors for the density profile parameters
142  Par_string[cosmo_param.size()] = "Rt";
143  param_prior[cosmo_param.size()] = Rt_prior;
144  Par_string[cosmo_param.size()+1] = "concentration";
145  param_prior[cosmo_param.size()+1] = concentration_prior;
146  Par_string[cosmo_param.size()+2] = "logM";
147  param_prior[cosmo_param.size()+2] = logM_prior;
148  Par_string[cosmo_param.size()+3] = "f_off";
149  param_prior[cosmo_param.size()+3] = f_off_prior;
150  Par_string[cosmo_param.size()+4] = "sigma_off";
151  param_prior[cosmo_param.size()+4] = sigma_off_prior;
152 
153  Par_string[cosmo_param.size()+5] = "AB_fact";
154  param_prior[cosmo_param.size()+5] = anisotropic_boost_prior;
155  Par_string[cosmo_param.size()+6] = "OB_fact";
156  param_prior[cosmo_param.size()+6] = orientation_boost_prior;
157 
158  // set the parameter indices, used in the model function
159  m_data_model.i_Rt = cosmo_param.size();
160  m_data_model.i_conc = cosmo_param.size()+1;
161  m_data_model.i_logM = cosmo_param.size()+2;
162  m_data_model.i_foff = cosmo_param.size()+3;
163  m_data_model.i_sigmaoff = cosmo_param.size()+4;
164  m_data_model.i_AB = cosmo_param.size()+5;
165  m_data_model.i_OB = cosmo_param.size()+6;
166 
167  m_data_model.i_Rt_func = 0;
168  m_data_model.i_foff_func = 1;
169  m_data_model.i_sigmaoff_func = 2;
170  m_data_model.i_AB_func = 3;
171  m_data_model.i_OB_func = 4;
172 
173  // set the functions returning the model parameters
174  m_data_model.get_parameter.resize(5);
175  m_data_model.priors_excluded.resize(5);
176 
177  for (size_t i=0; i<m_data_model.get_parameter.size(); i++)
178  m_data_model.get_parameter[i] = [] (std::vector<double> &par, const int idx, std::vector<statistics::PriorDistribution> prior, const int i_prior) {
179  (void)prior; (void)i_prior; return par[idx];
180  };
181 
182  m_data_model.priors_excluded[0] = Rt_prior;
183  m_data_model.priors_excluded[1] = f_off_prior;
184  m_data_model.priors_excluded[2] = sigma_off_prior;
185  m_data_model.priors_excluded[3] = anisotropic_boost_prior;
186  m_data_model.priors_excluded[4] = orientation_boost_prior;
187 
188  // Build the HaloProfile object
189  cosmology::HaloProfile halo_profile (*(m_data_model.cosmology), m_data_model.redshift, 0., 0., m_Delta, m_profile_author, m_halo_def, 0., true, false, 0., 0.);
190  m_data_model.halo_profile = make_shared<cosmology::HaloProfile>(halo_profile);
191 
192  // Set the function returning the concentration
193  m_data_model.conc_func = [] (const double conc, cbl::cosmology::HaloProfile halo_profile) {(void)halo_profile; return conc;};
194 
195  // Set the function returning the 2-halo term
196  if (m_2halo)
197  m_data_model.two_halo_func = [] (const std::vector<double> radius, cbl::cosmology::HaloProfile halo_profile, const std::string bias_author, const std::string method_Pk, const std::string interp_type) {
198  return halo_profile.DeltaSigma_2h(radius, bias_author, method_Pk, interp_type);};
199  else
200  m_data_model.two_halo_func = [] (const std::vector<double> radius, cbl::cosmology::HaloProfile halo_profile, const std::string bias_author, const std::string method_Pk, const std::string interp_type) {
201  (void)halo_profile; (void)bias_author; (void)method_Pk; (void)interp_type;
202  std::vector<double> res (radius.size(), 0.);
203  return res;};
204 
205  // input data used to construct the model
206  auto inputs = make_shared<STR_Profile_data_model>(m_data_model);
207 
208  // set prior
209  m_set_prior(param_prior);
210 
211  m_data_model.Par_type = Par_type;
212  m_data_model.Par_string = Par_string;
213 
214  // construct the model
215  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&model_density, nParams, Par_type, Par_string, inputs));
216 }
217 
218 // ===========================================================================================
219 
220 void cbl::modelling::densityprofile::Modelling_DensityProfile::set_model_DensityProfile_cosmology (const std::vector<cbl::cosmology::CosmologicalParameter> cosmo_param, const std::vector<statistics::PriorDistribution> cosmo_prior, const statistics::PriorDistribution Rt_prior, const std::string cM_author, const statistics::PriorDistribution logM_prior, const statistics::PriorDistribution f_off_prior, const statistics::PriorDistribution sigma_off_prior, const statistics::PriorDistribution anisotropic_boost_prior, const statistics::PriorDistribution orientation_boost_prior)
221 {
222  m_data_model.Cpar = cosmo_param;
223 
224  const size_t nParams = cosmo_param.size()+7; // The total number of parameters is given by the cosmological ones + 7, since the density profile has 6 base parameters (Rt, logM, f_off, sigma_off, AB_fact, OB_fact) and 1 derived parameter (conc)
225  const int n_derivedPars = 1;
226 
227  vector<statistics::ParameterType> Par_type (nParams, statistics::ParameterType::_Base_);
228  Par_type[cosmo_param.size()+1] = statistics::ParameterType::_Derived_;
229 
230  vector<string> Par_string (nParams);
231  std::vector<statistics::PriorDistribution> param_prior (nParams-n_derivedPars);
232 
233  // Set the names and priors of the cosmological parameters
234  for (size_t i=0; i<cosmo_param.size(); i++){
235  Par_string[i] = CosmologicalParameter_name(cosmo_param[i]);
236  param_prior[i] = cosmo_prior[i];
237  }
238 
239  // Set the names and priors for the density profile parameters
240  Par_string[cosmo_param.size()] = "Rt";
241  param_prior[cosmo_param.size()] = Rt_prior;
242 
243  Par_string[cosmo_param.size()+1] = "concentration";
244 
245  Par_string[cosmo_param.size()+2] = "logM";
246  param_prior[cosmo_param.size()+1] = logM_prior;
247  Par_string[cosmo_param.size()+3] = "f_off";
248  param_prior[cosmo_param.size()+2] = f_off_prior;
249  Par_string[cosmo_param.size()+4] = "sigma_off";
250  param_prior[cosmo_param.size()+3] = sigma_off_prior;
251 
252  Par_string[cosmo_param.size()+5] = "AB_fact";
253  param_prior[cosmo_param.size()+4] = anisotropic_boost_prior;
254  Par_string[cosmo_param.size()+6] = "OB_fact";
255  param_prior[cosmo_param.size()+5] = orientation_boost_prior;
256 
257  // set the parameter indices, used in the model function
258  m_data_model.i_Rt = cosmo_param.size();
259  m_data_model.i_conc = cosmo_param.size()+1;
260  m_data_model.i_logM = cosmo_param.size()+2;
261  m_data_model.i_foff = cosmo_param.size()+3;
262  m_data_model.i_sigmaoff = cosmo_param.size()+4;
263  m_data_model.i_AB = cosmo_param.size()+5;
264  m_data_model.i_OB = cosmo_param.size()+6;
265 
266  m_data_model.i_Rt_func = 0;
267  m_data_model.i_foff_func = 1;
268  m_data_model.i_sigmaoff_func = 2;
269  m_data_model.i_AB_func = 3;
270  m_data_model.i_OB_func = 4;
271 
272  // set the functions returning the model parameters
273  m_data_model.get_parameter.resize(5);
274  m_data_model.priors_excluded.resize(5);
275 
276  for (size_t i=0; i<m_data_model.get_parameter.size(); i++)
277  m_data_model.get_parameter[i] = [] (std::vector<double> &par, const int idx, std::vector<statistics::PriorDistribution> prior, const int i_prior) {
278  (void)prior; (void)i_prior; return par[idx];
279  };
280 
281  m_data_model.priors_excluded[0] = Rt_prior;
282  m_data_model.priors_excluded[1] = f_off_prior;
283  m_data_model.priors_excluded[2] = sigma_off_prior;
284  m_data_model.priors_excluded[3] = anisotropic_boost_prior;
285  m_data_model.priors_excluded[4] = orientation_boost_prior;
286 
287  // Build the HaloProfile object
288  cosmology::HaloProfile halo_profile (*(m_data_model.cosmology), m_data_model.redshift, cM_author, 0., m_Delta, m_profile_author, m_halo_def, 0., true, false, 0., 0.);
289  m_data_model.halo_profile = make_shared<cosmology::HaloProfile>(halo_profile);
290 
291  // Set the function returning the concentration
292  m_data_model.conc_func = [] (const double conc, cbl::cosmology::HaloProfile halo_profile) { (void)conc; return halo_profile.concentration(); };
293 
294  // Set the function returning the 2-halo term
295  if (m_2halo)
296  m_data_model.two_halo_func = [] (const std::vector<double> radius, cbl::cosmology::HaloProfile halo_profile, const std::string bias_author, const std::string method_Pk, const std::string interp_type) {
297  return halo_profile.DeltaSigma_2h(radius, bias_author, method_Pk, interp_type);};
298  else
299  m_data_model.two_halo_func = [] (const std::vector<double> radius, cbl::cosmology::HaloProfile halo_profile, const std::string bias_author, const std::string method_Pk, const std::string interp_type) {
300  (void)halo_profile; (void)bias_author; (void)method_Pk; (void)interp_type;
301  std::vector<double> res (radius.size(), 0.);
302  return res;};
303 
304  // input data used to construct the model
305  auto inputs = make_shared<STR_Profile_data_model>(m_data_model);
306 
307  // set prior
308  m_set_prior(param_prior);
309 
310  m_data_model.Par_type = Par_type;
311  m_data_model.Par_string = Par_string;
312 
313  // construct the model
314  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&model_density, nParams, Par_type, Par_string, inputs));
315 }
316 
317 
318 // ===========================================================================================
319 
320 void cbl::modelling::densityprofile::Modelling_DensityProfile::set_model_DensityProfile_cosmology (const std::vector<cbl::cosmology::CosmologicalParameter> cosmo_param, const std::vector<statistics::PriorDistribution> cosmo_prior, const std::string z_evo, const statistics::PriorDistribution Rt_prior, const statistics::PriorDistribution concentration_prior, const statistics::PriorDistribution f_off_prior, const statistics::PriorDistribution sigma_off_prior, const statistics::PriorDistribution anisotropic_boost_prior, const statistics::PriorDistribution orientation_boost_prior, const statistics::PriorDistribution alpha_prior, const statistics::PriorDistribution beta_prior, const statistics::PriorDistribution gamma_prior, const statistics::PriorDistribution scatter0_prior, const statistics::PriorDistribution scatterM_prior, const statistics::PriorDistribution scatterM_exponent_prior, const statistics::PriorDistribution scatterz_prior, const statistics::PriorDistribution scatterz_exponent_prior)
321 {
322  if (m_mass_is_derived == false)
323  ErrorCBL("If the mass is derived from the scaling relation, you must use the correct set_data_model!", "set_model_DensityProfile_cosmology", "Modelling_DensityProfile.cpp");
324 
325  m_data_model.Cpar = cosmo_param;
326 
327  const size_t nParams = cosmo_param.size()+14; // The total number of parameters is given by the cosmological ones, + 14 base parameters.
328 
329  vector<statistics::ParameterType> Par_type (nParams, statistics::ParameterType::_Base_);
330 
331  vector<string> Par_string (nParams);
332  std::vector<statistics::PriorDistribution> param_prior (nParams);
333 
334  // Set the names and priors of the cosmological parameters
335  for (size_t i=0; i<cosmo_param.size(); i++){
336  Par_string[i] = CosmologicalParameter_name(cosmo_param[i]);
337  param_prior[i] = cosmo_prior[i];
338  }
339 
340  // Set the names and priors for the density profile parameters
341  Par_string[cosmo_param.size()] = "Rt";
342  param_prior[cosmo_param.size()] = Rt_prior;
343  Par_string[cosmo_param.size()+1] = "concentration";
344  param_prior[cosmo_param.size()+1] = concentration_prior;
345 
346  Par_string[cosmo_param.size()+2] = "f_off";
347  param_prior[cosmo_param.size()+2] = f_off_prior;
348  Par_string[cosmo_param.size()+3] = "sigma_off";
349  param_prior[cosmo_param.size()+3] = sigma_off_prior;
350 
351  Par_string[cosmo_param.size()+4] = "AB_fact";
352  param_prior[cosmo_param.size()+4] = anisotropic_boost_prior;
353  Par_string[cosmo_param.size()+5] = "OB_fact";
354  param_prior[cosmo_param.size()+5] = orientation_boost_prior;
355 
356  Par_string[cosmo_param.size()+6] = "alpha";
357  param_prior[cosmo_param.size()+6] = alpha_prior;
358  Par_string[cosmo_param.size()+7] = "beta";
359  param_prior[cosmo_param.size()+7] = beta_prior;
360  Par_string[cosmo_param.size()+8] = "gamma";
361  param_prior[cosmo_param.size()+8] = gamma_prior;
362  Par_string[cosmo_param.size()+9] = "scatter0";
363  param_prior[cosmo_param.size()+9] = scatter0_prior;
364  Par_string[cosmo_param.size()+10] = "scatterM";
365  param_prior[cosmo_param.size()+10] = scatterM_prior;
366  Par_string[cosmo_param.size()+11] = "scatterM_exponent";
367  param_prior[cosmo_param.size()+11] = scatterM_exponent_prior;
368  Par_string[cosmo_param.size()+12] = "scatterz";
369  param_prior[cosmo_param.size()+12] = scatterz_prior;
370  Par_string[cosmo_param.size()+13] = "scatterz_exponent";
371  param_prior[cosmo_param.size()+13] = scatterz_exponent_prior;
372 
373  // set the parameter indices, used in the model function
374  m_data_model.i_Rt = cosmo_param.size();
375  m_data_model.i_conc = cosmo_param.size()+1;
376  m_data_model.i_foff = cosmo_param.size()+2;
377  m_data_model.i_sigmaoff = cosmo_param.size()+3;
378  m_data_model.i_AB = cosmo_param.size()+4;
379  m_data_model.i_OB = cosmo_param.size()+5;
380 
381  m_data_model.i_Rt_func = 0;
382  m_data_model.i_foff_func = 1;
383  m_data_model.i_sigmaoff_func = 2;
384  m_data_model.i_AB_func = 3;
385  m_data_model.i_OB_func = 4;
386 
387  // set the functions returning the model parameters
388  m_data_model.get_parameter.resize(5);
389  m_data_model.priors_excluded.resize(5);
390 
391  for (size_t i=0; i<m_data_model.get_parameter.size(); i++)
392  m_data_model.get_parameter[i] = [] (std::vector<double> &par, const int idx, std::vector<statistics::PriorDistribution> prior, const int i_prior) {
393  (void)prior; (void)i_prior; return par[idx];
394  };
395 
396  m_data_model.priors_excluded[0] = Rt_prior;
397  m_data_model.priors_excluded[1] = f_off_prior;
398  m_data_model.priors_excluded[2] = sigma_off_prior;
399  m_data_model.priors_excluded[3] = anisotropic_boost_prior;
400  m_data_model.priors_excluded[4] = orientation_boost_prior;
401 
402  // Build the HaloProfile object
403  cosmology::HaloProfile halo_profile (*(m_data_model.cosmology), m_data_model.redshift, 0., 0., m_Delta, m_profile_author, m_halo_def, 0., true, false, 0., 0.);
404  m_data_model.halo_profile = make_shared<cosmology::HaloProfile>(halo_profile);
405 
406  // Set the scaling relation object
407  (m_data_model.scaling_relation)->set_model_MassObservableRelation_cosmology (z_evo, cosmo_param, cosmo_prior, alpha_prior, beta_prior, gamma_prior, scatter0_prior, scatterM_prior, scatterM_exponent_prior, scatterz_prior, scatterz_exponent_prior);
408 
409  // Set the likelihood for the scaling relation (only to avoid internal errors, of course it is not used)
410  (m_data_model.scaling_relation)->set_likelihood(cbl::statistics::LikelihoodType::_Gaussian_Error_, {});
411 
412  // Set the function returning the concentration
413  m_data_model.conc_scaling_relation_func = [] (const double conc, const double c0=0, const double cM=0, const double cz=0, const double logM=0, const double logz=0) {
414  (void)c0; (void)cM; (void)cz; (void)logM; (void)logz; return conc;
415  };
416 
417  // Set the function returning the 2-halo term
418  if (m_2halo)
419  m_data_model.two_halo_func_fast = [] (const std::vector<double> radius, cbl::cosmology::HaloProfile halo_profile, const double bias, const std::string method_Pk, const std::string interp_type) {
420  return halo_profile.DeltaSigma_2h(radius, bias, method_Pk, interp_type);};
421  else
422  m_data_model.two_halo_func_fast = [] (const std::vector<double> radius, cbl::cosmology::HaloProfile halo_profile, const double bias, const std::string method_Pk, const std::string interp_type) {
423  (void)halo_profile; (void)bias; (void)method_Pk; (void)interp_type;
424  std::vector<double> res (radius.size(), 0.);
425  return res;};
426 
427  // input data used to construct the model
428  auto inputs = make_shared<STR_Profile_data_model>(m_data_model);
429 
430  // set prior
431  m_set_prior(param_prior);
432 
433  m_data_model.Par_type = Par_type;
434  m_data_model.Par_string = Par_string;
435 
436  // construct the model
437  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&model_density_scaling_relation, nParams, Par_type, Par_string, inputs));
438 }
439 
440 
441 // ===========================================================================================
442 
443 void cbl::modelling::densityprofile::Modelling_DensityProfile::set_model_DensityProfile_cosmology (const std::vector<cbl::cosmology::CosmologicalParameter> cosmo_param, const std::vector<statistics::PriorDistribution> cosmo_prior, const std::string z_evo, const statistics::PriorDistribution Rt_prior, const statistics::PriorDistribution c0_prior, const statistics::PriorDistribution cM_prior, const statistics::PriorDistribution cz_prior, const statistics::PriorDistribution f_off0_prior, const statistics::PriorDistribution f_offM_prior, const statistics::PriorDistribution f_offz_prior, const statistics::PriorDistribution sigma_off0_prior, const statistics::PriorDistribution sigma_offM_prior, const statistics::PriorDistribution sigma_offz_prior, const statistics::PriorDistribution anisotropic_boost_prior, const statistics::PriorDistribution orientation_boost_prior, const statistics::PriorDistribution alpha_prior, const statistics::PriorDistribution beta_prior, const statistics::PriorDistribution gamma_prior, const statistics::PriorDistribution scatter0_prior, const statistics::PriorDistribution scatterM_prior, const statistics::PriorDistribution scatterM_exponent_prior, const statistics::PriorDistribution scatterz_prior, const statistics::PriorDistribution scatterz_exponent_prior)
444 {
445  if (m_mass_is_derived == false)
446  ErrorCBL("If the mass is derived from the scaling relation, you must use the correct set_data_model!", "set_model_DensityProfile_cosmology", "Modelling_DensityProfile.cpp");
447 
448  m_data_model.Cpar = cosmo_param;
449 
450  const size_t nParams = cosmo_param.size()+20; // The total number of parameters is given by the cosmological ones, + 20 base parameters.
451 
452  vector<statistics::ParameterType> Par_type (nParams, statistics::ParameterType::_Base_);
453 
454  vector<string> Par_string (nParams);
455  std::vector<statistics::PriorDistribution> param_prior (nParams);
456 
457  // Set the names and priors of the cosmological parameters
458  for (size_t i=0; i<cosmo_param.size(); i++){
459  Par_string[i] = CosmologicalParameter_name(cosmo_param[i]);
460  param_prior[i] = cosmo_prior[i];
461  }
462 
463  // Set the names and priors for the density profile parameters
464  Par_string[cosmo_param.size()] = "Rt";
465  param_prior[cosmo_param.size()] = Rt_prior;
466 
467  Par_string[cosmo_param.size()+1] = "c0";
468  param_prior[cosmo_param.size()+1] = c0_prior;
469  Par_string[cosmo_param.size()+2] = "cM";
470  param_prior[cosmo_param.size()+2] = cM_prior;
471  Par_string[cosmo_param.size()+3] = "cz";
472  param_prior[cosmo_param.size()+3] = cz_prior;
473 
474  Par_string[cosmo_param.size()+4] = "f_off0";
475  param_prior[cosmo_param.size()+4] = f_off0_prior;
476  Par_string[cosmo_param.size()+5] = "f_offM";
477  param_prior[cosmo_param.size()+5] = f_offM_prior;
478  Par_string[cosmo_param.size()+6] = "f_offz";
479  param_prior[cosmo_param.size()+6] = f_offz_prior;
480  Par_string[cosmo_param.size()+7] = "sigma_off0";
481  param_prior[cosmo_param.size()+7] = sigma_off0_prior;
482  Par_string[cosmo_param.size()+8] = "sigma_offM";
483  param_prior[cosmo_param.size()+8] = sigma_offM_prior;
484  Par_string[cosmo_param.size()+9] = "sigma_offz";
485  param_prior[cosmo_param.size()+9] = sigma_offz_prior;
486 
487  Par_string[cosmo_param.size()+10] = "AB_fact";
488  param_prior[cosmo_param.size()+10] = anisotropic_boost_prior;
489  Par_string[cosmo_param.size()+11] = "OB_fact";
490  param_prior[cosmo_param.size()+11] = orientation_boost_prior;
491 
492  Par_string[cosmo_param.size()+12] = "alpha";
493  param_prior[cosmo_param.size()+12] = alpha_prior;
494  Par_string[cosmo_param.size()+13] = "beta";
495  param_prior[cosmo_param.size()+13] = beta_prior;
496  Par_string[cosmo_param.size()+14] = "gamma";
497  param_prior[cosmo_param.size()+14] = gamma_prior;
498  Par_string[cosmo_param.size()+15] = "scatter0";
499  param_prior[cosmo_param.size()+15] = scatter0_prior;
500  Par_string[cosmo_param.size()+16] = "scatterM";
501  param_prior[cosmo_param.size()+16] = scatterM_prior;
502  Par_string[cosmo_param.size()+17] = "scatterM_exponent";
503  param_prior[cosmo_param.size()+17] = scatterM_exponent_prior;
504  Par_string[cosmo_param.size()+18] = "scatterz";
505  param_prior[cosmo_param.size()+18] = scatterz_prior;
506  Par_string[cosmo_param.size()+19] = "scatterz_exponent";
507  param_prior[cosmo_param.size()+19] = scatterz_exponent_prior;
508 
509  // Build the HaloProfile object
510  cosmology::HaloProfile halo_profile (*(m_data_model.cosmology), m_data_model.redshift, 0., 0., m_Delta, m_profile_author, m_halo_def, 0., true, false, 0., 0.);
511  m_data_model.halo_profile = make_shared<cosmology::HaloProfile>(halo_profile);
512 
513  // Set the scaling relation object
514  (m_data_model.scaling_relation)->set_model_MassObservableRelation_cosmology (z_evo, cosmo_param, cosmo_prior, alpha_prior, beta_prior, gamma_prior, scatter0_prior, scatterM_prior, scatterM_exponent_prior, scatterz_prior, scatterz_exponent_prior);
515 
516  // Set the likelihood for the scaling relation (only to avoid internal errors, of course it is not used)
517  (m_data_model.scaling_relation)->set_likelihood(cbl::statistics::LikelihoodType::_Gaussian_Error_, {});
518 
519  // Set the function returning the concentration
520  m_data_model.conc_scaling_relation_func = [] (const double conc, const double c0, const double cM, const double cz, const double logM, const double logz) {
521  (void)conc;
522  const double logc = c0 + cM*logM + cz*logz;
523  return pow(10, logc);
524  };
525 
526  // Set the function returning the 2-halo term
527  if (m_2halo)
528  m_data_model.two_halo_func_fast = [] (const std::vector<double> radius, cbl::cosmology::HaloProfile halo_profile, const double bias, const std::string method_Pk, const std::string interp_type) {
529  return halo_profile.DeltaSigma_2h(radius, bias, method_Pk, interp_type);};
530  else
531  m_data_model.two_halo_func_fast = [] (const std::vector<double> radius, cbl::cosmology::HaloProfile halo_profile, const double bias, const std::string method_Pk, const std::string interp_type) {
532  (void)halo_profile; (void)bias; (void)method_Pk; (void)interp_type;
533  std::vector<double> res (radius.size(), 0.);
534  return res;};
535 
536  // input data used to construct the model
537  auto inputs = make_shared<STR_Profile_data_model>(m_data_model);
538 
539  // set prior
540  m_set_prior(param_prior);
541 
542  m_data_model.Par_type = Par_type;
543  m_data_model.Par_string = Par_string;
544 
545  // construct the model
546  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&model_density_scaling_relation_evolving_concentration_offcentering, nParams, Par_type, Par_string, inputs));
547 }
548 
549 
550 // ===========================================================================================
551 
552 void cbl::modelling::densityprofile::Modelling_DensityProfile::exclude_parameter_from_MCMC (const std::string parameter) {
553 
554  for (size_t i=0; i<m_data_model.Par_type.size(); i++) {
555 
556  switch (m_data_model.Par_type[i])
557  {
558 
559  case (statistics::ParameterType::_Base_):
560 
561  break;
562 
563  default:
564  ErrorCBL("Work in progress. It is not possible to use this function if a derived parameter is set!", "exclude_parameter_from_MCMC", "Modelling_DensityProfile.cpp");
565  break;
566 
567  }
568 
569  }
570 
571  // »»»
572  if (m_evolving_offcentering)
573  ErrorCBL("Work in progress. If the offcentering parameters evolve with redshift and mass proxy, this function can not be used!", "exclude_parameter_from_MCMC", "Modelling_DensityProfile.cpp");
574 
575  if (m_likelihood != NULL)
576  ErrorCBL("This function must be used before setting the likelihood!", "exclude_parameter_from_MCMC", "Modelling_DensityProfile.cpp");
577 
578  // »»»
579  std::vector<statistics::PriorDistribution> priors;
580  for (size_t i=0; i<m_parameter_priors.size(); i++)
581  priors.push_back(*m_parameter_priors[i]);
582 
583  if (parameter == "Rt") {
584 
585  m_data_model.Par_type.erase(m_data_model.Par_type.begin()+m_data_model.i_Rt);
586  m_data_model.Par_string.erase(m_data_model.Par_string.begin()+m_data_model.i_Rt);
587  priors.erase(priors.begin()+m_data_model.i_Rt);
588 
589  m_data_model.get_parameter[m_data_model.i_Rt_func] = [] (std::vector<double> &par, const int idx, std::vector<statistics::PriorDistribution> prior, const int i_prior) {
590  (void)par; (void)idx;
591  srand(time(0));
592  return prior[i_prior].sample(rand());
593  };
594 
595  m_data_model.i_conc = m_data_model.i_conc - 1;
596  m_data_model.i_logM = m_data_model.i_logM - 1;
597  m_data_model.i_foff = m_data_model.i_foff - 1;
598  m_data_model.i_sigmaoff = m_data_model.i_sigmaoff - 1;
599  m_data_model.i_AB = m_data_model.i_AB - 1;
600  m_data_model.i_OB = m_data_model.i_OB - 1;
601 
602  }
603  else if (parameter == "f_off") {
604 
605  m_data_model.Par_type.erase(m_data_model.Par_type.begin()+m_data_model.i_foff);
606  m_data_model.Par_string.erase(m_data_model.Par_string.begin()+m_data_model.i_foff);
607  priors.erase(priors.begin()+m_data_model.i_foff);
608 
609  m_data_model.get_parameter[m_data_model.i_foff_func] = [] (std::vector<double> &par, const int idx, std::vector<statistics::PriorDistribution> prior, const int i_prior) {
610  (void)par; (void)idx;
611  srand(time(0));
612  return prior[i_prior].sample(rand());
613  };
614 
615  m_data_model.i_sigmaoff = m_data_model.i_sigmaoff - 1;
616  m_data_model.i_AB = m_data_model.i_AB - 1;
617  m_data_model.i_OB = m_data_model.i_OB - 1;
618 
619  }
620  else if (parameter == "sigma_off") {
621 
622  m_data_model.Par_type.erase(m_data_model.Par_type.begin()+m_data_model.i_sigmaoff);
623  m_data_model.Par_string.erase(m_data_model.Par_string.begin()+m_data_model.i_sigmaoff);
624  priors.erase(priors.begin()+m_data_model.i_sigmaoff);
625 
626  m_data_model.get_parameter[m_data_model.i_sigmaoff_func] = [] (std::vector<double> &par, const int idx, std::vector<statistics::PriorDistribution> prior, const int i_prior) {
627  (void)par; (void)idx;
628  srand(time(0));
629  return prior[i_prior].sample(rand());
630  };
631 
632  m_data_model.i_AB = m_data_model.i_AB - 1;
633  m_data_model.i_OB = m_data_model.i_OB - 1;
634 
635  }
636  else if (parameter == "AB_fact") {
637 
638  m_data_model.Par_type.erase(m_data_model.Par_type.begin()+m_data_model.i_AB);
639  m_data_model.Par_string.erase(m_data_model.Par_string.begin()+m_data_model.i_AB);
640  priors.erase(priors.begin()+m_data_model.i_AB);
641 
642  m_data_model.get_parameter[m_data_model.i_AB_func] = [] (std::vector<double> &par, const int idx, std::vector<statistics::PriorDistribution> prior, const int i_prior) {
643  (void)par; (void)idx;
644  srand(time(0));
645  return prior[i_prior].sample(rand());
646  };
647 
648  m_data_model.i_OB = m_data_model.i_OB - 1;
649 
650  }
651  else if (parameter == "OB_fact") {
652 
653  m_data_model.Par_type.erase(m_data_model.Par_type.begin()+m_data_model.i_OB);
654  m_data_model.Par_string.erase(m_data_model.Par_string.begin()+m_data_model.i_OB);
655  priors.erase(priors.begin()+m_data_model.i_OB);
656 
657  m_data_model.get_parameter[m_data_model.i_OB_func] = [] (std::vector<double> &par, const int idx, std::vector<statistics::PriorDistribution> prior, const int i_prior) {
658  (void)par; (void)idx;
659  srand(time(0));
660  return prior[i_prior].sample(rand());
661  };
662 
663  }
664 
665  else
666  ErrorCBL("Wrong parameter declaration ("+parameter+").", "exclude_parameter_from_MCMC", "Modelling_DensityProfile.cpp");
667 
668  auto inputs = make_shared<STR_Profile_data_model>(m_data_model);
669  m_set_prior(priors);
670 
671  m_model.reset();
672 
673  if (m_mass_is_derived)
674  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&model_density_scaling_relation, m_data_model.Par_string.size(), m_data_model.Par_type, m_data_model.Par_string, inputs));
675  else
676  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&model_density, m_data_model.Par_string.size(), m_data_model.Par_type, m_data_model.Par_string, inputs));
677 
678  WarningMsgCBL("New set of MCMC parameters:", "exclude_parameter_from_MCMC", "Modelling_DensityProfile.cpp");
679  for (size_t i=0; i<m_data_model.Par_string.size(); i++)
680  coutCBL << m_data_model.Par_string[i] << endl;
681 
682 }
683 
684 // ===========================================================================================
685 
686 std::vector<double> cbl::modelling::densityprofile::model_density (const std::vector<double> radius, const std::shared_ptr<void> inputs, std::vector<double> &parameter)
687 {
688  // structure contaning the required input data
689  shared_ptr<STR_Profile_data_model> pp = static_pointer_cast<STR_Profile_data_model>(inputs);
690 
691  // redefine the cosmology
692  cbl::cosmology::Cosmology cosmo = *pp->cosmology;
693 
694  // redefine the HaloProfile object
695  cbl::cosmology::HaloProfile halo_profile = *pp->halo_profile;
696 
697  // set the cosmological parameters
698  for (size_t i=0; i<pp->Cpar.size(); ++i)
699  cosmo.set_parameter(pp->Cpar[i], parameter[i]);
700 
701  // set the cluster parameters
702  halo_profile.set_cosmology(cosmo);
703 
704  const double Rt = pp->get_parameter[pp->i_Rt_func](parameter, pp->i_Rt, pp->priors_excluded, pp->i_Rt_func);
705  const double f_off = pp->get_parameter[pp->i_foff_func](parameter, pp->i_foff, pp->priors_excluded, pp->i_foff_func);
706  const double sigma_off = pp->get_parameter[pp->i_sigmaoff_func](parameter, pp->i_sigmaoff, pp->priors_excluded, pp->i_sigmaoff_func);
707  const double AB_fact = pp->get_parameter[pp->i_AB_func](parameter, pp->i_AB, pp->priors_excluded, pp->i_AB_func);
708  const double OB_fact = pp->get_parameter[pp->i_OB_func](parameter, pp->i_OB, pp->priors_excluded, pp->i_OB_func);
709 
710  const double mass = pow(pp->logM_base, parameter[pp->i_logM])*pp->mass_pivot;
711 
712  halo_profile.set_trunc_fact(Rt);
713  halo_profile.set_mass(mass * (1+OB_fact));
714  halo_profile.set_f_off(f_off);
715  halo_profile.set_sigma_off(sigma_off);
716 
717  halo_profile.set_concentration(pp->conc_func(parameter[pp->i_conc], halo_profile));
718  parameter[pp->i_conc] = halo_profile.concentration();
719 
720  // Compute DeltaSigma
721  std::vector<double> one_halo = halo_profile.DeltaSigma(radius);
722  std::vector<double> two_halo = pp->two_halo_func(radius, halo_profile, pp->bias_author, pp->method_Pk, pp->interp_type);
723 
724  std::vector<double> total_profile (radius.size());
725  for (size_t j=0; j<radius.size(); j++)
726  total_profile[j] = one_halo[j] + two_halo[j] * (1.+AB_fact);
727 
728  return total_profile;
729 }
730 
731 // ===========================================================================================
732 
733 std::vector<double> cbl::modelling::densityprofile::compute_model_density_scaling_relation (const std::vector<double> radius, cbl::cosmology::Cosmology cosmo, cbl::cosmology::HaloProfile halo_profile, shared_ptr<STR_Profile_data_model> pp, const double conc, const double c0, const double cM, const double cz, const double AB_fact, const double OB_fact, const double alpha, const double beta, const double gamma, const double scatter0, const double scatterM, const double scatterM_exp, const double scatterz, const double scatterz_exp)
734 {
735  auto cosmo_ptr = std::make_shared<cbl::cosmology::Cosmology>(cosmo);
736 
737  // Define the integrand
738  double sqrt_Nclusters = sqrt( (pp->scaling_relation)->data_model().Nclusters[0] );
739 
740  cbl::glob::FuncGrid DeltaSigma_Rj_interp;
741  std::shared_ptr<void> ptr;
742 
743  auto integrand = [&] (const double x)
744  {
745  double mass = pow(pp->logM_base,x)*pp->mass_pivot;
746 
747  // Compute P(M|lambda,z)
748  double log_lambda = log(pp->mass_proxy/pp->proxy_pivot)/log(pp->logM_base);
749  double log_f_z = log( (pp->scaling_relation)->data_model().fz(pp->redshift, pp->redshift_pivot, cosmo_ptr) )/log(pp->logM_base);
750 
751  double mean = alpha + beta*log_lambda + gamma*log_f_z;
752  double scatter_intr = std::abs(scatter0 + scatterM*pow(log_lambda, scatterM_exp) + scatterz*pow(log_f_z, scatterz_exp)) / sqrt_Nclusters;
753  double P_M__lambda_z = cbl::gaussian(x, ptr, {mean,scatter_intr});
754 
755  // Compute the halo excess surface profile
756  double DeltaSigma = DeltaSigma_Rj_interp(mass);
757 
758  return DeltaSigma * P_M__lambda_z;
759  };
760 
761  // Find the minimum and maximum masses, given the parameters of the scaling relation and the intrinsic scatter
762  double log_lambda = log(pp->mass_proxy/pp->proxy_pivot)/log(pp->logM_base);
763  double log_f_z = log( (pp->scaling_relation)->data_model().fz(pp->redshift, pp->redshift_pivot, cosmo_ptr) )/log(pp->logM_base);
764 
765  double logM = alpha + beta*log_lambda + gamma*log_f_z;
766 
767  double scatter = std::abs( scatter0 + scatterM*pow(log_lambda, scatterM_exp) + scatterz*pow(log_f_z, scatterz_exp) ) / sqrt_Nclusters;
768 
769  double min_logM = logM-3.5*scatter;
770  double max_logM = logM+3.5*scatter;
771 
772  // »»»»»»»»»»»»»»»»»»
773  // Compute DeltaSigma
774  // »»»»»»»»»»»»»»»»»»
775 
776  const std::vector<double> Mass_vector = cbl::logarithmic_bin_vector(5, pow(pp->logM_base,min_logM)*pp->mass_pivot * (1-std::abs(OB_fact)), pow(pp->logM_base,max_logM)*pp->mass_pivot * (1+std::abs(OB_fact)));
777 
778  // For the halo bias in the two-halo term, compute:
779  // - Delta
780  // - growth factor
781  // - the interpolated sigmaM
782  double Delta = halo_profile.Delta() / cosmo.OmegaM(pp->redshift);
783  double DN = cosmo.DN(pp->redshift);
784 
785  std::vector<double> sigmaM (Mass_vector.size(), 0.);
786  for (size_t i=0; i<sigmaM.size(); i++)
787  sigmaM[i] = sqrt( cosmo.sigma2M(Mass_vector[i], pp->method_Pk, 0., false, "test", "Linear", 100.) );
788  cbl::glob::FuncGrid sigmaM_interp (Mass_vector, sigmaM, "Spline");
789 
790 
791  // »»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»
792  // Compute DeltaSigma between maximum and minimum masses, and interpolate it for a given radius
793  // »»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»
794 
795  std::vector<std::vector<double>> total_profile_4interp(Mass_vector.size(), std::vector<double>(radius.size()));
796 
797  // 0) Since DeltaSigma_2h / bias_halo is constant, compute it only once
798  const std::vector<double> normalised_2h = pp->two_halo_func_fast(radius, halo_profile, 1., pp->method_Pk, pp->interp_type);
799 
800  // 1) For a given mass, compute DeltaSigma at all radii
801  for (size_t i=0; i<total_profile_4interp.size(); i++) {
802 
803  double mass = Mass_vector[i] * (1+OB_fact);
804  halo_profile.set_mass(mass);
805  halo_profile.set_concentration( pp->conc_scaling_relation_func( conc, c0, cM, cz, log10(mass/pp->mass_pivot), log10((1+pp->redshift)/(1+pp->redshift_pivot)) ) );
806 
807  double bias = cosmo.bias_halo(mass, sigmaM_interp(mass), pp->redshift, DN, pp->bias_author, false, par::defaultString, "Linear", Delta, -1, -1, 0.001, 100., 1.e-2, pp->method_Pk);
808 
809  std::vector<double> one_halo = halo_profile.DeltaSigma(radius);
810  std::vector<double> two_halo = normalised_2h;
811 
812  for (size_t j=0; j<total_profile_4interp[i].size(); j++)
813  total_profile_4interp[i][j] = one_halo[j] + two_halo[j] * bias * (1.+AB_fact);
814 
815  }
816 
817  std::vector<double> total_profile(radius.size());
818 
819  for (size_t j=0; j<radius.size(); j++) {
820 
821  // 2) Given a radius, interpolate DeltaSigma as a function of mass
822  std::vector<double> DeltaSigma_Rj_4interp(Mass_vector.size());
823 
824  for (size_t i=0; i<Mass_vector.size(); i++)
825  DeltaSigma_Rj_4interp[i] = total_profile_4interp[i][j];
826 
827  DeltaSigma_Rj_interp = cbl::glob::FuncGrid (Mass_vector, DeltaSigma_Rj_4interp, "Spline");
828 
829  // Integrate
830  total_profile[j] = wrapper::gsl::GSL_integrate_qag(integrand, min_logM, max_logM);
831 
832  }
833 
834  return total_profile;
835 }
836 
837 // ===========================================================================================
838 
839 std::vector<double> cbl::modelling::densityprofile::model_density_scaling_relation (const std::vector<double> radius, const std::shared_ptr<void> inputs, std::vector<double> &parameter)
840 {
841  // structure contaning the required input data
842  shared_ptr<STR_Profile_data_model> pp = static_pointer_cast<STR_Profile_data_model>(inputs);
843 
844  // redefine the cosmology
845  cbl::cosmology::Cosmology cosmo = *pp->cosmology;
846 
847  // redefine the HaloProfile object
848  cbl::cosmology::HaloProfile halo_profile = *pp->halo_profile;
849 
850  // set the cosmological parameters
851  for (size_t i=0; i<pp->Cpar.size(); ++i)
852  cosmo.set_parameter(pp->Cpar[i], parameter[i]);
853 
854  // set the cluster parameters
855  halo_profile.set_cosmology(cosmo);
856 
857  const double Rt = pp->get_parameter[pp->i_Rt_func](parameter, pp->i_Rt, pp->priors_excluded, pp->i_Rt_func);
858  const double f_off = pp->get_parameter[pp->i_foff_func](parameter, pp->i_foff, pp->priors_excluded, pp->i_foff_func);
859  const double sigma_off = pp->get_parameter[pp->i_sigmaoff_func](parameter, pp->i_sigmaoff, pp->priors_excluded, pp->i_sigmaoff_func);
860  const double AB_fact = pp->get_parameter[pp->i_AB_func](parameter, pp->i_AB, pp->priors_excluded, pp->i_AB_func);
861  const double OB_fact = pp->get_parameter[pp->i_OB_func](parameter, pp->i_OB, pp->priors_excluded, pp->i_OB_func);
862 
863  halo_profile.set_trunc_fact(Rt);
864  halo_profile.set_f_off(f_off);
865  halo_profile.set_sigma_off(sigma_off);
866 
867  const double conc = parameter[pp->i_conc];
868 
869  // set the scaling relation parameters
870  const double alpha = parameter[parameter.size()-8];
871  const double beta = parameter[parameter.size()-7];
872  const double gamma = parameter[parameter.size()-6];
873  const double scatter0 = parameter[parameter.size()-5];
874  const double scatterM = parameter[parameter.size()-4];
875  const double scatterM_exp = parameter[parameter.size()-3];
876  const double scatterz = parameter[parameter.size()-2];
877  const double scatterz_exp = parameter[parameter.size()-1];
878 
879  return cbl::modelling::densityprofile::compute_model_density_scaling_relation(radius, cosmo, halo_profile, pp, conc, 0., 0., 0., AB_fact, OB_fact, alpha, beta, gamma, scatter0, scatterM, scatterM_exp, scatterz, scatterz_exp);
880 }
881 
882 // ===========================================================================================
883 
884 std::vector<double> cbl::modelling::densityprofile::model_density_scaling_relation_evolving_concentration_offcentering (const std::vector<double> radius, const std::shared_ptr<void> inputs, std::vector<double> &parameter)
885 {
886  // structure contaning the required input data
887  shared_ptr<STR_Profile_data_model> pp = static_pointer_cast<STR_Profile_data_model>(inputs);
888 
889  // redefine the cosmology
890  cbl::cosmology::Cosmology cosmo = *pp->cosmology;
891 
892  // redefine the HaloProfile object
893  cbl::cosmology::HaloProfile halo_profile = *pp->halo_profile;
894 
895  // set the cosmological parameters
896  for (size_t i=0; i<pp->Cpar.size(); ++i)
897  cosmo.set_parameter(pp->Cpar[i], parameter[i]);
898 
899  // set the cluster parameters
900  halo_profile.set_cosmology(cosmo);
901 
902  const double Rt = parameter[pp->Cpar.size()];
903  const double c0 = parameter[pp->Cpar.size()+1];
904  const double cM = parameter[pp->Cpar.size()+2];
905  const double cz = parameter[pp->Cpar.size()+3];
906  const double f_off = parameter[pp->Cpar.size()+4] * pow(pp->mass_proxy/pp->proxy_pivot, parameter[pp->Cpar.size()+5]) * pow((1+pp->redshift)/(1+pp->redshift_pivot), parameter[pp->Cpar.size()+6]);
907  const double sigma_off = parameter[pp->Cpar.size()+7] * pow(pp->mass_proxy/pp->proxy_pivot, parameter[pp->Cpar.size()+8]) * pow((1+pp->redshift)/(1+pp->redshift_pivot), parameter[pp->Cpar.size()+9]);
908  const double AB_fact = parameter[pp->Cpar.size()+10];
909  const double OB_fact = parameter[pp->Cpar.size()+11];
910 
911  halo_profile.set_trunc_fact(Rt);
912  halo_profile.set_f_off(f_off);
913  halo_profile.set_sigma_off(sigma_off);
914 
915  // set the scaling relation parameters
916  const double alpha = parameter[parameter.size()-8];
917  const double beta = parameter[parameter.size()-7];
918  const double gamma = parameter[parameter.size()-6];
919  const double scatter0 = parameter[parameter.size()-5];
920  const double scatterM = parameter[parameter.size()-4];
921  const double scatterM_exp = parameter[parameter.size()-3];
922  const double scatterz = parameter[parameter.size()-2];
923  const double scatterz_exp = parameter[parameter.size()-1];
924 
925  return cbl::modelling::densityprofile::compute_model_density_scaling_relation(radius, cosmo, halo_profile, pp, 0., c0, cM, cz, AB_fact, OB_fact, alpha, beta, gamma, scatter0, scatterM, scatterM_exp, scatterz, scatterz_exp);
926 }
927 
Data1D.h
The class Data1D.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
Modelling_DensityProfile.h
The class Modelling_DensityProfile.
cbl::cosmology::Cosmology
The class Cosmology.
Definition: Cosmology.h:277
cbl::cosmology::Cosmology::DN
double DN(const double redshift, const double redshift_norm=0., const double prec=1.e-4) const
the normalised amplitude of the growing mode at a given redshift,
Definition: Cosmology.cpp:691
cbl::cosmology::Cosmology::set_parameter
void set_parameter(const CosmologicalParameter parameter, const double value)
set the value of one cosmological paramter
Definition: Cosmology.cpp:424
cbl::cosmology::Cosmology::bias_halo
double bias_halo(const double Mass, const double redshift, const std::string author, const std::string method_SS, const bool store_output=true, const std::string output_root="test", const std::string interpType="Linear", const double Delta=200., const double kk=-1., const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string input_file=par::defaultString, const bool is_parameter_file=true)
bias of dark matter haloes
Definition: Bias.cpp:44
cbl::cosmology::Cosmology::sigma2M
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
cbl::cosmology::Cosmology::OmegaM
double OmegaM(const double redshift=0.) const
the matter density at a given redshift
Definition: Cosmology.cpp:579
cbl::cosmology::HaloProfile
The class HaloProfile.
Definition: HaloProfile.h:52
cbl::cosmology::HaloProfile::set_mass
void set_mass(const double mass=par::defaultDouble)
set the private member m_mass
Definition: HaloProfile.h:977
cbl::cosmology::HaloProfile::set_concentration
void set_concentration(const double conc=par::defaultDouble)
set the private member m_concentration
Definition: HaloProfile.h:984
cbl::cosmology::HaloProfile::DeltaSigma_2h
std::vector< double > DeltaSigma_2h(const std::vector< double > rad, const std::string bias_author, const std::string method_Pk="EisensteinHu", const std::string interp_type="Linear", const bool NL=false)
the 2-halo contribution to the excess surface density profile of the halo, computed by assuming a hal...
Definition: HaloProfile.cpp:835
cbl::cosmology::HaloProfile::DeltaSigma
std::vector< double > DeltaSigma(const std::vector< double > rad)
the total excess surface density profile of the halo.
Definition: HaloProfile.cpp:757
cbl::cosmology::HaloProfile::concentration
double concentration()
get the concentration
Definition: HaloProfile.h:567
cbl::cosmology::HaloProfile::Delta
double Delta()
get the overdensity
Definition: HaloProfile.h:592
cbl::cosmology::HaloProfile::set_f_off
void set_f_off(const double f_off=par::defaultDouble)
set the private member m_f_off
Definition: HaloProfile.h:991
cbl::cosmology::HaloProfile::set_sigma_off
void set_sigma_off(const double sigma_off=par::defaultDouble)
set the private member m_sigma_off
Definition: HaloProfile.h:998
cbl::cosmology::HaloProfile::set_trunc_fact
void set_trunc_fact(const double trunc_fact=par::defaultDouble)
set the private member m_trunc_fact
Definition: HaloProfile.h:1005
cbl::cosmology::HaloProfile::set_cosmology
void set_cosmology(const cbl::cosmology::Cosmology cosmology)
set the cosmological model
Definition: HaloProfile.cpp:638
cbl::data::Data1D
The class Data1D.
Definition: Data1D.h:51
cbl::glob::FuncGrid
The class FuncGrid.
Definition: FuncGrid.h:55
cbl::modelling::densityprofile::Modelling_DensityProfile::Modelling_DensityProfile
Modelling_DensityProfile()=default
default constuctor _DensityProfile
cbl::modelling::densityprofile::Modelling_DensityProfile::set_model_DensityProfile_cosmology
void set_model_DensityProfile_cosmology(const std::vector< cbl::cosmology::CosmologicalParameter > cosmo_param, const std::vector< statistics::PriorDistribution > cosmo_prior, const statistics::PriorDistribution Rt_prior, const statistics::PriorDistribution concentration_prior, const statistics::PriorDistribution logM_prior, const statistics::PriorDistribution f_off_prior, const statistics::PriorDistribution sigma_off_prior, const statistics::PriorDistribution anisotropic_boost_prior, const statistics::PriorDistribution orientation_boost_prior)
Set the profile and cosmological parameters used to model the cluster density profile.
Definition: Modelling_DensityProfile.cpp:125
cbl::modelling::densityprofile::Modelling_DensityProfile::exclude_parameter_from_MCMC
void exclude_parameter_from_MCMC(const std::string parameter)
Exclude the selected parameter from the MCMC parameter space, and randomly extract its value from the...
Definition: Modelling_DensityProfile.cpp:552
cbl::modelling::densityprofile::Modelling_DensityProfile::set_data_model
void set_data_model(const cbl::cosmology::Cosmology cosmology, const double redshift, const double contrast, const double logM_base, const double mass_pivot, const std::string bias_author="Tinker", const std::string method_Pk="EisensteinHu", std::string interp_type="Linear")
Set the data used to construct generic models of halo profiles.
Definition: Modelling_DensityProfile.cpp:73
cbl::modelling::massobsrel::Modelling_MassObservableRelation
The class Modelling_MassObservableRelation.
Definition: Modelling_MassObservableRelation.h:179
cbl::statistics::Model1D
The class Model1D.
Definition: Model1D.h:60
cbl::statistics::PriorDistribution
The class PriorDistribution.
Definition: PriorDistribution.h:53
cbl::par::defaultString
static const std::string defaultString
default std::string value
Definition: Constants.h:336
cbl::par::alpha
static const double alpha
: the fine-structure constant
Definition: Constants.h:233
cbl::cosmology::CosmologicalParameter_name
std::string CosmologicalParameter_name(const CosmologicalParameter parameter)
name of the cosmological parameter
Definition: Cosmology.cpp:45
cbl::modelling::densityprofile::model_density_scaling_relation
std::vector< double > model_density_scaling_relation(const std::vector< double > radius, const std::shared_ptr< void > inputs, std::vector< double > &parameter)
Compute the excess density profile model in all the radial bins. The mass is derived from the scaling...
Definition: Modelling_DensityProfile.cpp:839
cbl::modelling::densityprofile::model_density
std::vector< double > model_density(const std::vector< double > radius, const std::shared_ptr< void > inputs, std::vector< double > &parameter)
Compute the excess density profile model in all the radial bins.
Definition: Modelling_DensityProfile.cpp:686
cbl::modelling::densityprofile::model_density_scaling_relation_evolving_concentration_offcentering
std::vector< double > model_density_scaling_relation_evolving_concentration_offcentering(const std::vector< double > radius, const std::shared_ptr< void > inputs, std::vector< double > &parameter)
Compute the excess density profile model in all the radial bins. The mass is derived from the scaling...
Definition: Modelling_DensityProfile.cpp:884
cbl::modelling::densityprofile::compute_model_density_scaling_relation
std::vector< double > compute_model_density_scaling_relation(const std::vector< double > radius, cbl::cosmology::Cosmology cosmo, cbl::cosmology::HaloProfile halo_profile, std::shared_ptr< STR_Profile_data_model > pp, const double conc, const double c0, const double cM, const double cz, const double AB_fact, const double OB_fact, const double alpha, const double beta, const double gamma, const double scatter0, const double scatterM, const double scatterM_exp, const double scatterz, const double scatterz_exp)
Compute the excess density profile model in all the radial bins when mass is derived from the scaling...
Definition: Modelling_DensityProfile.cpp:733
cbl::modelling::massobsrel::scaling_relation
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
Definition: Modelling_MassObservableRelation.cpp:143
cbl::modelling::twopt::bias
double bias(const double Mmin, const double sigmalgM, const double M0, const double M1, const double alpha, const std::shared_ptr< void > inputs)
the mean galaxy bias
Definition: ModelFunction_TwoPointCorrelation1D_monopole.cpp:1026
cbl::statistics::LikelihoodType::_Gaussian_Error_
@ _Gaussian_Error_
Gaussian likelihood error.
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::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::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::gaussian
T gaussian(T xx, std::shared_ptr< void > pp, std::vector< double > par)
the Gaussian function
Definition: Func.h:1127
cbl::WarningMsgCBL
void WarningMsgCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL)
internal CBL warning message
Definition: Kernel.h:747