CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Modelling_TwoPointCorrelation_wedges.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.h"
38 #include "ModelFunction_TwoPointCorrelation_wedges.h"
39 #include "Modelling_TwoPointCorrelation_wedges.h"
40 
41 using namespace std;
42 
43 using namespace cbl;
44 
45 
46 // ============================================================================================
47 
48 
49 cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::Modelling_TwoPointCorrelation_wedges (const std::shared_ptr<cbl::measure::twopt::TwoPointCorrelation> twop, const int nWedges, const std::vector<std::vector<double>> mu_integral_limits)
50  : Modelling_TwoPointCorrelation1D_monopole(twop), m_nWedges(nWedges), m_mu_integral_limits(mu_integral_limits)
51 {
52  m_ModelIsSet = false;
53 
54  m_wedges_order.erase(m_wedges_order.begin(), m_wedges_order.end());
55 
56  for (int j=0; j<m_nWedges; j++)
57  for (int i=0; i<m_data->ndata()/m_nWedges; i++)
58  m_wedges_order.push_back(j);
59 
60  m_use_wedge.resize(m_nWedges, true);
61 }
62 
63 
64 // ============================================================================================
65 
66 
67 cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::Modelling_TwoPointCorrelation_wedges (const std::shared_ptr<data::Data> twop_dataset, const int nWedges, const std::vector<std::vector<double>> mu_integral_limits)
68  : Modelling_TwoPointCorrelation1D_monopole(twop_dataset), m_nWedges(nWedges), m_mu_integral_limits(mu_integral_limits)
69 {
70  m_ModelIsSet = false;
71 
72  m_wedges_order.erase(m_wedges_order.begin(), m_wedges_order.end());
73  m_use_wedge.resize(m_nWedges, false);
74 
75  for (int j=0; j<m_nWedges; j++) {
76  m_use_wedge[j] = true;
77  for (int i=0; i<m_data->ndata()/m_nWedges; i++)
78  m_wedges_order.push_back(j);
79  }
80 
81 }
82 
83 
84 // ============================================================================================
85 
86 
87 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_fit_range (const double xmin, const double xmax, const int nWedges)
88 {
89  vector<vector<double>> fit_range(m_nWedges, vector<double>(2, -1.));
90 
91  int mp = (0<nWedges && nWedges<m_nWedges) ? nWedges : m_nWedges;
92 
93  for (int i=0; i<mp; i++) {
94  fit_range[i][0] = xmin;
95  fit_range[i][1] = xmax;
96  }
97 
98  set_fit_range(fit_range);
99 }
100 
101 
102 // ============================================================================================
103 
104 
105 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_fit_range (const std::vector<std::vector<double>> fit_range)
106 {
107  if ((int)fit_range.size()!=m_nWedges)
108  ErrorCBL("wrong number of wedges provided, "+conv(fit_range.size(), cbl::par::fINT)+" instead of "+conv(m_nWedges, cbl::par::fINT)+"!", "set_fit_range", "Modelling_TwoPointCorrelation_wedges.cpp");
109 
110  m_use_wedge = {false, false};
111 
112  m_wedges_order.erase(m_wedges_order.begin(), m_wedges_order.end());
113 
114  const int size = m_data->ndata()/m_nWedges;
115  vector<bool> mask(m_data->ndata(), false);
116  vector<double> xx;
117 
118  for (int j=0; j<m_nWedges; j++) {
119  for (int i=0; i<size; i++) {
120  if (fit_range[j][0]<m_data->xx(i+j*size) && m_data->xx(i+j*size)<fit_range[j][1]) {
121  m_wedges_order.push_back(j);
122  xx.push_back(m_data->xx(i+j*size));
123  m_use_wedge[j] = true;
124  mask[i+j*size] = true;
125  }
126  }
127  }
128  m_data_fit = m_data->cut(mask);
129 
130  m_fit_range = true;
131 
132  if (m_ModelIsSet)
133  m_data_model->dataset_order = m_wedges_order;
134 
135 }
136 
137 
138 // ============================================================================================
139 
140 
141 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_fiducial_PkDM ()
142 {
143  m_data_model->nmultipoles = 3;
144  m_data_model->nWedges = 2;
145 
146  m_data_model->kk = logarithmic_bin_vector(m_data_model->step, max(m_data_model->k_min, 1.e-4), min(m_data_model->k_max, 500.));
147  vector<double> Pk = m_data_model->cosmology->Pk_matter(m_data_model->kk, m_data_model->method_Pk, false, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec, m_data_model->file_par);
148 
149  if (m_data_model->Pk_mu_model=="dispersion_dewiggled") {
150  vector<double> PkNW = m_data_model->cosmology->Pk_matter(m_data_model->kk, "EisensteinHu", false, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec, m_data_model->file_par);
151  m_data_model->func_Pk = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk, "Spline"));
152  m_data_model->func_Pk_NW = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, PkNW, "Spline"));
153 
154  m_data_model->funcs_pk.push_back(m_data_model->func_Pk);
155  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_NW);
156  }
157 
158  else if (m_data_model->Pk_mu_model=="dispersion_modecoupling") {
159  vector<double> kk_1loop, Pk_1loop;
160  for (size_t i=0; i<(size_t)m_data_model->step; i++) {
161  if (m_data_model->kk[i]<par::pi) {
162  kk_1loop.push_back(m_data_model->kk[i]);
163  Pk_1loop.push_back(m_data_model->cosmology->Pk_1loop(m_data_model->kk[i], m_data_model->func_Pk, 0, m_data_model->k_min, 5., m_data_model->prec));
164  }
165  }
166  m_data_model->func_Pk = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk, "Spline"));
167  m_data_model->func_Pk1loop = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(kk_1loop, Pk_1loop, "Spline"));
168 
169  m_data_model->funcs_pk.push_back(m_data_model->func_Pk);
170  m_data_model->funcs_pk.push_back(m_data_model->func_Pk1loop);
171  }
172 
173  else if (m_data_model->Pk_mu_model=="dispersion_Gauss" || m_data_model->Pk_mu_model=="dispersion_Lorentz") {
174  m_data_model->func_Pk = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk, "Spline"));
175 
176  m_data_model->funcs_pk.push_back(m_data_model->func_Pk);
177  }
178 
179  else if (m_data_model->Pk_mu_model=="Scoccimarro_Pezzotta_Gauss" || m_data_model->Pk_mu_model=="Scoccimarro_Pezzotta_Lorentz" || m_data_model->Pk_mu_model=="Scoccimarro_Bel_Gauss" || m_data_model->Pk_mu_model=="Scoccimarro_Bel_Lorentz") {
180  vector<double> Pknonlin = m_data_model->cosmology->Pk_matter(m_data_model->kk, m_data_model->method_Pk, true, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec, m_data_model->file_par);
181  m_data_model->func_Pk = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk, "Spline"));
182  m_data_model->func_Pk_nonlin = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pknonlin, "Spline"));
183 
184  m_data_model->funcs_pk.push_back(m_data_model->func_Pk);
185  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_nonlin);
186  }
187 
188  else if (m_data_model->Pk_mu_model=="Scoccimarro_Gauss" || m_data_model->Pk_mu_model=="Scoccimarro_Lorentz") {
189  vector<vector<double>> Pk_terms = m_data_model->cosmology->Pk_TNS_dd_dt_tt(m_data_model->kk, m_data_model->method_Pk, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec);
190 
191  m_data_model->func_Pk_DeltaDelta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_terms[0], "Spline"));
192  m_data_model->func_Pk_DeltaTheta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_terms[1], "Spline"));
193  m_data_model->func_Pk_ThetaTheta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_terms[2], "Spline"));
194 
195  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_DeltaDelta);
196  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_DeltaTheta);
197  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_ThetaTheta);
198  }
199 
200  else if (m_data_model->Pk_mu_model=="TNS_Gauss" || m_data_model->Pk_mu_model=="TNS_Lorentz") {
201  vector<vector<double>> Pk_terms = m_data_model->cosmology->Pk_TNS_dd_dt_tt(m_data_model->kk, m_data_model->method_Pk, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec);
202  vector<vector<double>> Pk_AB = m_data_model->cosmology->Pk_TNS_AB_terms_1loop(m_data_model->kk, m_data_model->method_Pk, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec);
203 
204  m_data_model->func_Pk_DeltaDelta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_terms[0], "Spline"));
205  m_data_model->func_Pk_DeltaTheta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_terms[1], "Spline"));
206  m_data_model->func_Pk_ThetaTheta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_terms[2], "Spline"));
207 
208  m_data_model->func_Pk_A11 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[0], "Spline"));
209  m_data_model->func_Pk_A12 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[1], "Spline"));
210  m_data_model->func_Pk_A22 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[2], "Spline"));
211  m_data_model->func_Pk_A23 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[3], "Spline"));
212  m_data_model->func_Pk_A33 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[4], "Spline"));
213  m_data_model->func_Pk_B12 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[5], "Spline"));
214  m_data_model->func_Pk_B13 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[6], "Spline"));
215  m_data_model->func_Pk_B14 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[7], "Spline"));
216  m_data_model->func_Pk_B22 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[8], "Spline"));
217  m_data_model->func_Pk_B23 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[9], "Spline"));
218  m_data_model->func_Pk_B24 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[10], "Spline"));
219  m_data_model->func_Pk_B33 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[11], "Spline"));
220  m_data_model->func_Pk_B34 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[12], "Spline"));
221  m_data_model->func_Pk_B44 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[13], "Spline"));
222 
223  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_DeltaDelta);
224  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_DeltaTheta);
225  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_ThetaTheta);
226  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A11);
227  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A12);
228  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A22);
229  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A23);
230  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A33);
231  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B12);
232  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B13);
233  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B14);
234  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B22);
235  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B23);
236  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B24);
237  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B33);
238  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B34);
239  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B44);
240  }
241 
242  else if (m_data_model->Pk_mu_model=="eTNS_Gauss" || m_data_model->Pk_mu_model=="eTNS_Lorentz") {
243  vector<vector<double>> Pk_AB = m_data_model->cosmology->Pk_TNS_AB_terms_1loop(m_data_model->kk, m_data_model->method_Pk, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec);
244  vector<vector<double>> Pk_eTNS_terms = m_data_model->cosmology->Pk_eTNS_terms_1loop(m_data_model->kk, m_data_model->method_Pk, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec);
245 
246  m_data_model->func_Pk_DeltaDelta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[0], "Spline"));
247  m_data_model->func_Pk_DeltaTheta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[1], "Spline"));
248  m_data_model->func_Pk_ThetaTheta = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[2], "Spline"));
249 
250  m_data_model->func_Pk_A11 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[0], "Spline"));
251  m_data_model->func_Pk_A12 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[1], "Spline"));
252  m_data_model->func_Pk_A22 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[2], "Spline"));
253  m_data_model->func_Pk_A23 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[3], "Spline"));
254  m_data_model->func_Pk_A33 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[4], "Spline"));
255  m_data_model->func_Pk_B12 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[5], "Spline"));
256  m_data_model->func_Pk_B13 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[6], "Spline"));
257  m_data_model->func_Pk_B14 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[7], "Spline"));
258  m_data_model->func_Pk_B22 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[8], "Spline"));
259  m_data_model->func_Pk_B23 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[9], "Spline"));
260  m_data_model->func_Pk_B24 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[10], "Spline"));
261  m_data_model->func_Pk_B33 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[11], "Spline"));
262  m_data_model->func_Pk_B34 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[12], "Spline"));
263  m_data_model->func_Pk_B44 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_AB[13], "Spline"));
264 
265  m_data_model->func_Pk_b2d = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[3], "Spline"));
266  m_data_model->func_Pk_b2v = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[4], "Spline"));
267  m_data_model->func_Pk_b22 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[5], "Spline"));
268  m_data_model->func_Pk_bs2d = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[6], "Spline"));
269  m_data_model->func_Pk_bs2v = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[7], "Spline"));
270  m_data_model->func_Pk_b2s2 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[8], "Spline"));
271  m_data_model->func_Pk_bs22 = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[9], "Spline"));
272  m_data_model->func_sigma32Pklin = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk_eTNS_terms[10], "Spline"));
273 
274  //-------------------------
275 
276  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_DeltaDelta);
277  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_DeltaTheta);
278  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_ThetaTheta);
279  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A11);
280  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A12);
281  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A22);
282  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A23);
283  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_A33);
284  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B12);
285  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B13);
286  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B14);
287  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B22);
288  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B23);
289  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B24);
290  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B33);
291  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B34);
292  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_B44);
293 
294  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_b2d);
295  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_b2v);
296  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_b22);
297  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_bs2d);
298  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_bs2v);
299  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_b2s2);
300  m_data_model->funcs_pk.push_back(m_data_model->func_Pk_bs22);
301  m_data_model->funcs_pk.push_back(m_data_model->func_sigma32Pklin);
302  }
303 
304  else ErrorCBL("the chosen model ("+m_data_model->Pk_mu_model+") is not currently implemented!", "set_fiducial_PkDM", "Modelling_TwoPointCorrelation_wedges.cpp");
305 }
306 
307 
308 // ============================================================================================
309 
310 
311 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_fiducial_xiDM ()
312 {
313  cout << endl; coutCBL << "Setting up the fiducial dark matter two-point correlation function model" << endl;
314 
315  m_data_model->nmultipoles = 3;
316  m_data_model->nWedges = 2;
317 
318  const vector<double> rad = linear_bin_vector(m_data_model->step, m_data_model->r_min, m_data_model->r_max);
319 
320  m_data_model->rr = rad;
321  m_data_model->kk = logarithmic_bin_vector(m_data_model->step, max(m_data_model->k_min, 1.e-4), min(m_data_model->k_max, 500.));
322 
323  vector<double> Pk = m_data_model->cosmology->Pk_matter(m_data_model->kk, m_data_model->method_Pk, false, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec, m_data_model->file_par);
324  vector<double> PkNW = m_data_model->cosmology->Pk_matter(m_data_model->kk, "EisensteinHu", false, m_data_model->redshift, m_data_model->store_output, m_data_model->output_root, m_data_model->norm, m_data_model->k_min, m_data_model->k_max, m_data_model->prec, m_data_model->file_par);
325 
326 
327  m_data_model->func_Pk = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, Pk, "Spline"));
328  m_data_model->func_Pk_NW = make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(m_data_model->kk, PkNW, "Spline"));
329 
330  std::vector<double> template_parameter = {m_data_model->sigmaNL_perp, m_data_model->sigmaNL_par, m_data_model->linear_growth_rate_z, m_data_model->bias, 0.};
331  cbl::Print(template_parameter);
332 
333  vector<vector<double>> xil = Xi_l(rad, m_data_model->nmultipoles, "dispersion_dewiggled", {m_data_model->sigmaNL_perp, m_data_model->sigmaNL_par, m_data_model->linear_growth_rate_z, m_data_model->bias, 0.}, {m_data_model->func_Pk, m_data_model->func_Pk_NW}, m_data_model->prec, 1, 1);
334 
335  m_data_model->func_multipoles.erase(m_data_model->func_multipoles.begin(), m_data_model->func_multipoles.end());
336 
337  for (int i=0; i< m_data_model->nmultipoles; i++)
338  m_data_model->func_multipoles.push_back(make_shared<cbl::glob::FuncGrid>(cbl::glob::FuncGrid(rad, xil[i], "Spline")));
339 
340 }
341 
342 
343 // ============================================================================================
344 
345 
346 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::write_model (const std::string output_dir, const std::string output_file, const std::vector<double> xx, const std::vector<double> parameter)
347 {
348  if (m_likelihood==NULL) ErrorCBL("this function requires the likelihood to be defined (with the function set_likelihood)!", "write_model", "Modelling_TwoPointCorrelation_wedges.cpp");
349 
350  vector<int> dataset_order_original = m_data_model->dataset_order;
351 
352  vector<int> new_dataset_order;
353  vector<double> new_xx;
354 
355  if (xx.size()==0)
356  new_xx = m_data_fit->xx();
357  else
358  for (int n=0; n<m_data_model->nWedges; n++)
359  for (size_t i=0; i<xx.size(); i++) {
360  new_xx.push_back(xx[i]);
361  new_dataset_order.push_back(n);
362  }
363 
364  m_data_model->dataset_order = new_dataset_order;
365  m_likelihood->write_model(output_dir, output_file, parameter, new_xx);
366  m_data_model->dataset_order = dataset_order_original;
367 }
368 
369 
370 // ============================================================================================
371 
372 
373 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::write_model_at_bestfit (const std::string output_dir, const std::string output_file, const std::vector<double> xx)
374 {
375  if (m_posterior==NULL)
376  ErrorCBL("no posterior found: run maximize_posterior first!", "write_model_at_bestfit", "Modelling_TwoPointCorrelation_wedges.cpp");
377 
378  vector<int> dataset_order_original = m_data_model->dataset_order;
379 
380  vector<int> new_dataset_order;
381  vector<double> new_xx;
382 
383  if (xx.size()==0)
384  new_xx = m_data_fit->xx();
385  else
386  for (int n=0; n<m_data_model->nWedges; n++)
387  for (size_t i=0; i<xx.size(); i++) {
388  new_xx.push_back(xx[i]);
389  new_dataset_order.push_back(n);
390  }
391 
392  m_data_model->dataset_order = new_dataset_order;
393  m_posterior->write_model_at_bestfit(output_dir, output_file, new_xx);
394 
395  m_data_model->dataset_order = dataset_order_original;
396 }
397 
398 
399 // ============================================================================================
400 
401 
402 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::write_model_from_chains (const std::string output_dir, const std::string output_file, const std::vector<double> xx, const int start, const int thin)
403 {
404  if (m_posterior==NULL)
405  ErrorCBL("no posterior found: run sample_posterior first!", "write_model_from_chains", "Modelling_TwoPointCorrelation_wedges.cpp");
406 
407  vector<int> dataset_order_original = m_data_model->dataset_order;
408 
409  vector<int> new_dataset_order;
410  vector<double> new_xx;
411 
412  if (xx.size()==0)
413  new_xx = m_data_fit->xx();
414  else
415  for (int n=0; n<m_data_model->nWedges; n++)
416  for (size_t i=0; i<xx.size(); i++) {
417  new_xx.push_back(xx[i]);
418  new_dataset_order.push_back(n);
419  }
420 
421  m_data_model->dataset_order = new_dataset_order;
422  m_posterior->write_model_from_chain(output_dir, output_file, new_xx, {}, start, thin);
423 
424  m_data_model->dataset_order = dataset_order_original;
425 }
426 
427 
428 // ============================================================================================
429 
430 
431 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_fullShape_DeWiggled (const statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const statistics::PriorDistribution SigmaNL_perpendicular_prior, const statistics::PriorDistribution SigmaNL_parallel_prior, const statistics::PriorDistribution fsigma8_prior, const statistics::PriorDistribution bsigma8_prior, const statistics::PriorDistribution SigmaS_prior, const bool compute_PkDM)
432 {
433  m_data_model->Pk_mu_model = "dispersion_dewiggled";
434 
435  // compute the fiducial dark matter power spectrum terms used to construct the model
436  if (compute_PkDM) set_fiducial_PkDM();
437 
438  // number of wedges to be used
439  m_data_model->nWedges = m_nWedges;
440 
441  // integral limits used to measure the wedges
442  m_data_model->mu_integral_limits = m_mu_integral_limits;
443 
444  // scales to be used for each wedges
445  m_data_model->dataset_order = m_wedges_order;
446 
447  m_data_model->nmultipoles = 3;
448 
449  // set the model parameters
450  const int nparameters = 7;
451 
452  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
453 
454  vector<string> parameterName(nparameters);
455  parameterName[0] = "alpha_perpendicular";
456  parameterName[1] = "alpha_parallel";
457  parameterName[2] = "SigmaNL_perpendicular";
458  parameterName[3] = "SigmaNL_parallel";
459  parameterName[4] = "f*sigma8";
460  parameterName[5] = "b*sigma8";
461  parameterName[6] = "Sigma_S";
462 
463  vector<statistics::PriorDistribution> priors = {alpha_perpendicular_prior, alpha_parallel_prior, SigmaNL_perpendicular_prior, SigmaNL_parallel_prior, fsigma8_prior, bsigma8_prior, SigmaS_prior};
464 
465  //set the priors
466  m_set_prior(priors);
467 
468  // construct the model
469  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges, nparameters, parameterType, parameterName, m_data_model));
470  m_ModelIsSet = true;
471 }
472 
473 // ============================================================================================
474 
475 
476 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_fullShape_ModeCoupling (const statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const statistics::PriorDistribution fsigma8_prior, const statistics::PriorDistribution bsigma8_prior, const statistics::PriorDistribution SigmaV_prior, const statistics::PriorDistribution AMC_prior, const bool compute_PkDM)
477 {
478  m_data_model->Pk_mu_model = "dispersion_modecoupling";
479 
480  // compute the fiducial dark matter power spectrum terms used to construct the model
481  if (compute_PkDM) set_fiducial_PkDM();
482 
483  // scales to be used for each wedges
484  m_data_model->dataset_order = m_wedges_order;
485 
486  m_data_model->nmultipoles = 3;
487 
488  // set the model parameters
489  const int nparameters = 6;
490 
491  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
492 
493  vector<string> parameterName(nparameters);
494  parameterName[0] = "alpha_perpendicular";
495  parameterName[1] = "alpha_parallel";
496  parameterName[2] = "f*sigma8";
497  parameterName[3] = "b*sigma8";
498  parameterName[4] = "sigma_v";
499  parameterName[5] = "AMC";
500 
501  vector<statistics::PriorDistribution> priors = {alpha_perpendicular_prior, alpha_parallel_prior, fsigma8_prior, bsigma8_prior, SigmaV_prior, AMC_prior};
502 
503  //set the priors
504  m_set_prior(priors);
505 
506  // construct the model
507  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges, nparameters, parameterType, parameterName, m_data_model));
508  m_ModelIsSet = true;
509 }
510 
511 
512 // ============================================================================================
513 
514 
515 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_BAO (const statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const statistics::PriorDistribution Bperp_prior, const statistics::PriorDistribution Bpar_prior, const statistics::PriorDistribution Aperp0_prior, const statistics::PriorDistribution Apar0_prior, const statistics::PriorDistribution Aperp1_prior, const statistics::PriorDistribution Apar1_prior, const statistics::PriorDistribution Aperp2_prior, const statistics::PriorDistribution Apar2_prior, const bool compute_XiTemplate, const bool isRealSpace)
516 {
517  // compute the fiducial dark matter two-point correlation function
518 
519  if (m_data_model->nWedges>2)
520  ErrorCBL("BAO modelling can be done only with two wedges!", "set_model_BAO", "Modelling_TwoPointCorrelation_wedges");
521 
522  if (isRealSpace) {
523  double lgf = m_data_model->linear_growth_rate_z;
524  m_data_model->linear_growth_rate_z = 0.;
525 
526  set_fiducial_xiDM();
527 
528  m_data_model->linear_growth_rate_z = lgf;
529  }
530  else if (compute_XiTemplate)
531  set_fiducial_xiDM();
532 
533  m_data_model->dataset_order = m_wedges_order;
534 
535  // set the model parameters
536  const int nparameters = 10;
537 
538  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
539 
540  vector<string> parameterName(nparameters);
541  parameterName[0] = "alpha_perpendicular";
542  parameterName[1] = "alpha_parallel";
543  parameterName[2] = "Bperp";
544  parameterName[3] = "Bpar";
545  parameterName[4] = "Aperp0";
546  parameterName[5] = "Apar0";
547  parameterName[6] = "Aperp1";
548  parameterName[7] = "Apar1";
549  parameterName[8] = "Aperp2";
550  parameterName[9] = "Apar2";
551 
552  vector<statistics::PriorDistribution> priors = {alpha_perpendicular_prior, alpha_parallel_prior, Bperp_prior, Bpar_prior, Aperp0_prior, Apar0_prior, Aperp1_prior, Apar1_prior, Aperp2_prior, Apar2_prior};
553 
554  //set the priors
555  m_set_prior(priors);
556 
557  // construct the model
558  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges_BAO, nparameters, parameterType, parameterName, m_data_model));
559  m_ModelIsSet = true;
560 }
561 
562 
563 // ============================================================================================
564 
565 
566 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_dispersion (const statistics::PriorDistribution fsigma8_prior, const statistics::PriorDistribution bsigma8_prior, const statistics::PriorDistribution sigmav_prior, const statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const bool DFoG, const bool compute_PkDM)
567 {
568  if (DFoG) m_data_model->Pk_mu_model = "dispersion_Gauss";
569  else m_data_model->Pk_mu_model = "dispersion_Lorentz";
570 
571  // compute the fiducial dark matter two-point correlation function
572  if (compute_PkDM) set_fiducial_PkDM();
573 
574  // number of wedges to be used
575  m_data_model->nWedges = m_nWedges;
576 
577  // integral limits used to measure the wedges
578  m_data_model->mu_integral_limits = m_mu_integral_limits;
579 
580  // scales to be used for each wedges
581  m_data_model->dataset_order = m_wedges_order;
582 
583  m_data_model->nmultipoles = 3;
584 
585  // set the model parameters
586  const int nparameters = 5;
587 
588  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
589 
590  vector<string> parameterName(nparameters);
591  parameterName[0] = "f*sigma8";
592  parameterName[1] = "b*sigma8";
593  parameterName[2] = "sigmav";
594  parameterName[3] = "alpha_perpendicular";
595  parameterName[4] = "alpha_parallel";
596 
597  vector<statistics::PriorDistribution> priors = {fsigma8_prior, bsigma8_prior, sigmav_prior, alpha_perpendicular_prior, alpha_parallel_prior};
598 
599  //set the priors
600  m_set_prior(priors);
601 
602  // construct the model
603  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges, nparameters, parameterType, parameterName, m_data_model));
604  m_ModelIsSet = true;
605 }
606 
607 // ============================================================================================
608 
609 
610 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_Scoccimarro (const statistics::PriorDistribution fsigma8_prior, const statistics::PriorDistribution bsigma8_prior, const statistics::PriorDistribution sigmav_prior, const statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const bool DFoG, const bool compute_PkDM)
611 {
612  if (DFoG) m_data_model->Pk_mu_model = "Scoccimarro_Gauss";
613  else m_data_model->Pk_mu_model = "Scoccimarro_Lorentz";
614 
615  // compute the fiducial dark matter two-point correlation function
616  if (compute_PkDM) set_fiducial_PkDM();
617 
618  // number of wedges to be used
619  m_data_model->nWedges = m_nWedges;
620 
621  // integral limits used to measure the wedges
622  m_data_model->mu_integral_limits = m_mu_integral_limits;
623 
624  // scales to be used for each wedges
625  m_data_model->dataset_order = m_wedges_order;
626 
627  m_data_model->nmultipoles = 3;
628 
629  // set the model parameters
630  const int nparameters = 5;
631 
632  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
633 
634  vector<string> parameterName(nparameters);
635  parameterName[0] = "f*sigma8";
636  parameterName[1] = "b*sigma8";
637  parameterName[2] = "sigmav";
638  parameterName[3] = "alpha_perpendicular";
639  parameterName[4] = "alpha_parallel";
640 
641  vector<statistics::PriorDistribution> priors = {fsigma8_prior, bsigma8_prior, sigmav_prior, alpha_perpendicular_prior, alpha_parallel_prior};
642 
643  //set the priors
644  m_set_prior(priors);
645 
646  // construct the model
647  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges, nparameters, parameterType, parameterName, m_data_model));
648  m_ModelIsSet = true;
649 }
650 
651 
652 // ============================================================================================
653 
654 
655 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_Scoccimarro_fitPezzotta (const statistics::PriorDistribution fsigma8_prior, const statistics::PriorDistribution bsigma8_prior, const statistics::PriorDistribution sigmav_prior, const statistics::PriorDistribution kd_prior, const statistics::PriorDistribution kt_prior, const statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const bool DFoG, const bool compute_PkDM)
656 {
657  if (DFoG) m_data_model->Pk_mu_model = "Scoccimarro_Pezzotta_Gauss";
658  else m_data_model->Pk_mu_model = "Scoccimarro_Pezzotta_Lorentz";
659 
660  // compute the fiducial dark matter two-point correlation function
661  if (compute_PkDM) set_fiducial_PkDM();
662 
663  // number of wedges to be used
664  m_data_model->nWedges = m_nWedges;
665 
666  // integral limits used to measure the wedges
667  m_data_model->mu_integral_limits = m_mu_integral_limits;
668 
669  // scales to be used for each wedges
670  m_data_model->dataset_order = m_wedges_order;
671 
672  m_data_model->nmultipoles = 3;
673 
674  // set the model parameters
675  const int nparameters = 7;
676 
677  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
678 
679  vector<string> parameterName(nparameters);
680  parameterName[0] = "f*sigma8";
681  parameterName[1] = "b*sigma8";
682  parameterName[2] = "sigmav";
683  parameterName[3] = "kd";
684  parameterName[4] = "kt";
685  parameterName[5] = "alpha_perpendicular";
686  parameterName[6] = "alpha_parallel";
687 
688  vector<statistics::PriorDistribution> priors = {fsigma8_prior, bsigma8_prior, sigmav_prior, kd_prior, kt_prior, alpha_perpendicular_prior, alpha_parallel_prior};
689 
690  //set the priors
691  m_set_prior(priors);
692 
693  // construct the model
694  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges, nparameters, parameterType, parameterName, m_data_model));
695  m_ModelIsSet = true;
696 }
697 
698 
699 // ============================================================================================
700 
701 
702 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_Scoccimarro_fitBel (const statistics::PriorDistribution fsigma8_prior, const statistics::PriorDistribution bsigma8_prior, const statistics::PriorDistribution sigmav_prior, const statistics::PriorDistribution kd_prior, const statistics::PriorDistribution bb_prior, const statistics::PriorDistribution a1_prior, const statistics::PriorDistribution a2_prior, const statistics::PriorDistribution a3_prior, const statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const bool DFoG, const bool compute_PkDM)
703 {
704  if (DFoG) m_data_model->Pk_mu_model = "Scoccimarro_Bel_Gauss";
705  else m_data_model->Pk_mu_model = "Scoccimarro_Bel_Lorentz";
706 
707  // compute the fiducial dark matter two-point correlation function
708  if (compute_PkDM) set_fiducial_PkDM();
709 
710  // number of wedges to be used
711  m_data_model->nWedges = m_nWedges;
712 
713  // integral limits used to measure the wedges
714  m_data_model->mu_integral_limits = m_mu_integral_limits;
715 
716  // scales to be used for each wedges
717  m_data_model->dataset_order = m_wedges_order;
718 
719  m_data_model->nmultipoles = 3;
720 
721  // set the model parameters
722  const int nparameters = 10;
723 
724  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
725 
726  vector<string> parameterName(nparameters);
727  parameterName[0] = "f*sigma8";
728  parameterName[1] = "b*sigma8";
729  parameterName[2] = "sigmav";
730  parameterName[3] = "kd";
731  parameterName[4] = "bb";
732  parameterName[5] = "a1";
733  parameterName[6] = "a2";
734  parameterName[7] = "a3";
735  parameterName[8] = "alpha_perpendicular";
736  parameterName[9] = "alpha_parallel";
737 
738  vector<statistics::PriorDistribution> priors = {fsigma8_prior, bsigma8_prior, sigmav_prior, kd_prior, bb_prior, a1_prior, a2_prior, a3_prior, alpha_perpendicular_prior, alpha_parallel_prior};
739 
740  //set the priors
741  m_set_prior(priors);
742 
743  // construct the model
744  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges, nparameters, parameterType, parameterName, m_data_model));
745  m_ModelIsSet = true;
746 }
747 
748 
749 // ============================================================================================
750 
751 
752 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_TNS (const statistics::PriorDistribution fsigma8_prior, const statistics::PriorDistribution bsigma8_prior, const statistics::PriorDistribution sigmav_prior, statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const bool DFoG, const bool compute_PkDM)
753 {
754  if (DFoG) m_data_model->Pk_mu_model = "TNS_Gauss";
755  else m_data_model->Pk_mu_model = "TNS_Lorentz";
756 
757  // compute the fiducial dark matter two-point correlation function
758  if (compute_PkDM) set_fiducial_PkDM();
759 
760  // number of wedges to be used
761  m_data_model->nWedges = m_nWedges;
762 
763  // integral limits used to measure the wedges
764  m_data_model->mu_integral_limits = m_mu_integral_limits;
765 
766  // scales to be used for each wedges
767  m_data_model->dataset_order = m_wedges_order;
768 
769  m_data_model->nmultipoles = 3;
770 
771  // set the model parameters
772  const int nparameters = 5;
773 
774  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
775 
776  vector<string> parameterName(nparameters);
777  parameterName[0] = "f*sigma8";
778  parameterName[1] = "b*sigma8";
779  parameterName[2] = "sigmav";
780  parameterName[3] = "alpha_perpendicular";
781  parameterName[4] = "alpha_parallel";
782 
783  vector<statistics::PriorDistribution> priors = {fsigma8_prior, bsigma8_prior, sigmav_prior, alpha_perpendicular_prior, alpha_parallel_prior};
784 
785  //set the priors
786  m_set_prior(priors);
787 
788  // construct the model
789  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges, nparameters, parameterType, parameterName, m_data_model));
790  m_ModelIsSet = true;
791 }
792 
793 
794 // ============================================================================================
795 
796 
797 void cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_eTNS (const statistics::PriorDistribution fsigma8_prior, const statistics::PriorDistribution b1sigma8_prior, const statistics::PriorDistribution b2sigma8_prior, const statistics::PriorDistribution sigmav_prior, statistics::PriorDistribution alpha_perpendicular_prior, const statistics::PriorDistribution alpha_parallel_prior, const statistics::PriorDistribution N_prior, const bool DFoG, const bool compute_PkDM)
798 {
799  if (DFoG) m_data_model->Pk_mu_model = "eTNS_Gauss";
800  else m_data_model->Pk_mu_model = "eTNS_Lorentz";
801 
802  // compute the fiducial dark matter two-point correlation function
803  if (compute_PkDM) set_fiducial_PkDM();
804 
805  // number of wedges to be used
806  m_data_model->nWedges = m_nWedges;
807 
808  // integral limits used to measure the wedges
809  m_data_model->mu_integral_limits = m_mu_integral_limits;
810 
811  // scales to be used for each wedges
812  m_data_model->dataset_order = m_wedges_order;
813 
814  m_data_model->nmultipoles = 3;
815 
816  // set the model parameters
817  const int nparameters = 7;
818 
819  vector<statistics::ParameterType> parameterType(nparameters, statistics::ParameterType::_Base_);
820 
821  vector<string> parameterName(nparameters);
822  parameterName[0] = "f*sigma8";
823  parameterName[1] = "b1*sigma8";
824  parameterName[2] = "b2*sigma8";
825  parameterName[3] = "sigmav";
826  parameterName[4] = "alpha_perpendicular";
827  parameterName[5] = "alpha_parallel";
828  parameterName[6] = "N";
829 
830  vector<statistics::PriorDistribution> priors = {fsigma8_prior, b1sigma8_prior, b2sigma8_prior, sigmav_prior, alpha_perpendicular_prior, alpha_parallel_prior, N_prior};
831 
832  //set the priors
833  m_set_prior(priors);
834 
835  // construct the model
836  m_model = make_shared<statistics::Model1D>(statistics::Model1D(&xiWedges, nparameters, parameterType, parameterName, m_data_model));
837  m_ModelIsSet = true;
838 }
839 
840 
Data1D.h
The class Data1D.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
ModelFunction_TwoPointCorrelation_wedges.h
Functions to model the wedges of the two-point correlation function.
Modelling_TwoPointCorrelation_wedges.h
The class Modelling_TwoPointCorrelatoin_wedges.
cbl::glob::FuncGrid
The class FuncGrid.
Definition: FuncGrid.h:55
cbl::modelling::Modelling::m_data
std::shared_ptr< data::Data > m_data
input data to be modelled
Definition: Modelling.h:69
cbl::modelling::twopt::Modelling_TwoPointCorrelation1D_monopole
The class Modelling_TwoPointCorrelation1D_monopole.
Definition: Modelling_TwoPointCorrelation1D_monopole.h:63
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_BAO
void set_model_BAO(const statistics::PriorDistribution alpha_perpendicular_prior={}, const statistics::PriorDistribution alpha_parallel_prior={}, const statistics::PriorDistribution Bperp_prior={}, const statistics::PriorDistribution Bpar_prior={}, const statistics::PriorDistribution Aperp0_prior={}, const statistics::PriorDistribution Apar0_prior={}, const statistics::PriorDistribution Aperp1_prior={}, const statistics::PriorDistribution Apar1_prior={}, const statistics::PriorDistribution Aperp2_prior={}, const statistics::PriorDistribution Apar2_prior={}, const bool compute_XiDM=true, const bool isRealSpace=false)
set the model to fit the wedges of the two-point correlation function, used for anisotropic BAO measu...
Definition: Modelling_TwoPointCorrelation_wedges.cpp:515
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::m_ModelIsSet
bool m_ModelIsSet
bolean to check if the model has been set
Definition: Modelling_TwoPointCorrelation_wedges.h:80
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::m_use_wedge
std::vector< bool > m_use_wedge
vector of booleans indicating the wedges to be modelled (m_use_wedge[i]=true -> the i-th wedge is mod...
Definition: Modelling_TwoPointCorrelation_wedges.h:74
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::write_model
virtual void write_model(const std::string output_dir, const std::string output_file, const std::vector< double > xx, const std::vector< double > parameter) override
write the model at xx, for the given parameters
Definition: Modelling_TwoPointCorrelation_wedges.cpp:346
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_Scoccimarro
void set_model_Scoccimarro(const statistics::PriorDistribution fsigma8_prior={}, const statistics::PriorDistribution bsigma8_prior={}, const statistics::PriorDistribution sigmav_prior={}, const statistics::PriorDistribution alpha_perpendicular_prior={cbl::glob::DistributionType::_Constant_, 1.}, const statistics::PriorDistribution alpha_parallel_prior={cbl::glob::DistributionType::_Constant_, 1.}, const bool DFoG=true, const bool compute_PkDM=true)
set the Scoccimarro model to fit the wedges of the two-point correlation function
Definition: Modelling_TwoPointCorrelation_wedges.cpp:610
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::m_nWedges
int m_nWedges
number of measured wedges in the input dataset
Definition: Modelling_TwoPointCorrelation_wedges.h:68
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_fiducial_PkDM
void set_fiducial_PkDM()
set the fiducial model for the dark matter power spectrum
Definition: Modelling_TwoPointCorrelation_wedges.cpp:141
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_Scoccimarro_fitPezzotta
void set_model_Scoccimarro_fitPezzotta(const statistics::PriorDistribution fsigma8_prior={}, const statistics::PriorDistribution bsigma8_prior={}, const statistics::PriorDistribution sigmav_prior={}, const statistics::PriorDistribution kd_prior={}, const statistics::PriorDistribution kt_prior={}, const statistics::PriorDistribution alpha_perpendicular_prior={cbl::glob::DistributionType::_Constant_, 1.}, const statistics::PriorDistribution alpha_parallel_prior={cbl::glob::DistributionType::_Constant_, 1.}, const bool DFoG=true, const bool compute_PkDM=true)
set the Scoccimarro model to fit the wedges of the two-point correlation function
Definition: Modelling_TwoPointCorrelation_wedges.cpp:655
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::m_wedges_order
std::vector< int > m_wedges_order
vector containing the ordering of the data vector, which spcifies which wedge correponds to each data...
Definition: Modelling_TwoPointCorrelation_wedges.h:77
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_dispersion
void set_model_dispersion(const statistics::PriorDistribution fsigma8_prior={}, const statistics::PriorDistribution bsigma8_prior={}, const statistics::PriorDistribution sigmav_prior={}, const statistics::PriorDistribution alpha_perpendicular_prior={cbl::glob::DistributionType::_Constant_, 1.}, const statistics::PriorDistribution alpha_parallel_prior={cbl::glob::DistributionType::_Constant_, 1.}, const bool DFoG=true, const bool compute_PkDM=true)
set the dispersion model to fit the wedges of the two-point correlation function
Definition: Modelling_TwoPointCorrelation_wedges.cpp:566
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_fullShape_ModeCoupling
void set_model_fullShape_ModeCoupling(const statistics::PriorDistribution alpha_perpendicular_prior={}, const statistics::PriorDistribution alpha_parallel_prior={}, const statistics::PriorDistribution fsigma8_prior={}, const statistics::PriorDistribution bsigma8_prior={}, const statistics::PriorDistribution SigmaV_prior={}, const statistics::PriorDistribution AMC_prior={}, const bool compute_PkDM=true)
set the mode-coupling model to fit the full shape of the wedges of the two-point correlation function
Definition: Modelling_TwoPointCorrelation_wedges.cpp:476
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_eTNS
void set_model_eTNS(const statistics::PriorDistribution fsigma8_prior={}, const statistics::PriorDistribution b1sigma8_prior={}, const statistics::PriorDistribution b2sigma8_prior={}, const statistics::PriorDistribution sigmav_prior={}, const statistics::PriorDistribution alpha_perpendicular_prior={cbl::glob::DistributionType::_Constant_, 1.}, const statistics::PriorDistribution alpha_parallel_prior={cbl::glob::DistributionType::_Constant_, 1.}, const statistics::PriorDistribution N_prior={cbl::glob::DistributionType::_Constant_, 0.}, const bool DFoG=true, const bool compute_PkDM=true)
set the eTNS model, i.e extended-TNS (Taruya, Nishimichi and Saito) model to fit the wedges of the tw...
Definition: Modelling_TwoPointCorrelation_wedges.cpp:797
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::write_model_from_chains
void write_model_from_chains(const std::string output_dir, const std::string output_file, const std::vector< double > xx, const int start=0, const int thin=1) override
write the model at xx computing 16th, 50th and 84th percentiles from the chains
Definition: Modelling_TwoPointCorrelation_wedges.cpp:402
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_TNS
void set_model_TNS(const statistics::PriorDistribution fsigma8_prior={}, const statistics::PriorDistribution bsigma8_prior={}, const statistics::PriorDistribution sigmav_prior={}, const statistics::PriorDistribution alpha_perpendicular_prior={cbl::glob::DistributionType::_Constant_, 1.}, const statistics::PriorDistribution alpha_parallel_prior={cbl::glob::DistributionType::_Constant_, 1.}, const bool DFoG=true, const bool compute_PkDM=true)
set the TNS (Taruya, Nishimichi and Saito) model to fit the wedges of the two-point correlation funct...
Definition: Modelling_TwoPointCorrelation_wedges.cpp:752
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_fit_range
void set_fit_range(const double xmin, const double xmax, const int nWedges=-1)
set the scale range used for the fit
Definition: Modelling_TwoPointCorrelation_wedges.cpp:87
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_fullShape_DeWiggled
void set_model_fullShape_DeWiggled(const statistics::PriorDistribution alpha_perpendicular_prior={}, const statistics::PriorDistribution alpha_parallel_prior={}, const statistics::PriorDistribution SigmaNL_perpendicular_prior={}, const statistics::PriorDistribution SigmaNL_parallel_prior={}, const statistics::PriorDistribution fsigma8_prior={}, const statistics::PriorDistribution bsigma8_prior={}, const statistics::PriorDistribution SigmaS_prior={}, const bool compute_PkDM=true)
set the de-wiggled model to fit the full shape of the wedges of the two-point correlation function
Definition: Modelling_TwoPointCorrelation_wedges.cpp:431
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::write_model_at_bestfit
void write_model_at_bestfit(const std::string output_dir, const std::string output_file, const std::vector< double > xx) override
Definition: Modelling_TwoPointCorrelation_wedges.cpp:373
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::Modelling_TwoPointCorrelation_wedges
Modelling_TwoPointCorrelation_wedges()=default
default constuctor
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_model_Scoccimarro_fitBel
void set_model_Scoccimarro_fitBel(const statistics::PriorDistribution fsigma8_prior={}, const statistics::PriorDistribution bsigma8_prior={}, const statistics::PriorDistribution sigmav_prior={}, const statistics::PriorDistribution kd_prior={}, const statistics::PriorDistribution bb_prior={}, const statistics::PriorDistribution a1_prior={}, const statistics::PriorDistribution a2_prior={}, const statistics::PriorDistribution a3_prior={}, const statistics::PriorDistribution alpha_perpendicular_prior={cbl::glob::DistributionType::_Constant_, 1.}, const statistics::PriorDistribution alpha_parallel_prior={cbl::glob::DistributionType::_Constant_, 1.}, const bool DFoG=true, const bool compute_PkDM=true)
set the Scoccimarro model to fit the wedges of the two-point correlation function
Definition: Modelling_TwoPointCorrelation_wedges.cpp:702
cbl::modelling::twopt::Modelling_TwoPointCorrelation_wedges::set_fiducial_xiDM
void set_fiducial_xiDM()
set the fiducial model for the dark matter two-point correlation function and associated quantities
Definition: Modelling_TwoPointCorrelation_wedges.cpp:311
cbl::statistics::Model1D
The class Model1D.
Definition: Model1D.h:60
cbl::statistics::PriorDistribution
The class PriorDistribution.
Definition: PriorDistribution.h:53
cbl::par::fINT
static const char fINT[]
conversion symbol for: int -> std::string
Definition: Constants.h:121
cbl::par::pi
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
cbl::par::mp
static const double mp
: the mass of the proton [kg]
Definition: Constants.h:266
cbl::modelling::twopt::xiWedges
std::vector< double > xiWedges(const std::vector< double > rad, const std::shared_ptr< void > inputs, std::vector< double > &parameter)
the model wedges of the two-point correlation function
Definition: ModelFunction_TwoPointCorrelation_wedges.cpp:94
cbl::modelling::twopt::Xi_l
std::vector< std::vector< double > > Xi_l(const std::vector< double > rr, const int nmultipoles, const std::string model, const std::vector< double > parameter, const std::vector< std::shared_ptr< glob::FuncGrid >> pk_interp, const double prec=1.e-5, const double alpha_perp=1., const double alpha_par=1.)
the multipole of order l of the two-point correlation function
Definition: ModelFunction_TwoPointCorrelation.cpp:453
cbl::modelling::twopt::xiWedges_BAO
std::vector< double > xiWedges_BAO(const std::vector< double > rad, const std::shared_ptr< void > inputs, std::vector< double > &parameter)
return the wedges of the two-point correlation function, intended for anisotropic BAO measurements
Definition: ModelFunction_TwoPointCorrelation_wedges.cpp:203
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
cbl::conv
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
cbl::Print
void Print(const T value, const int prec, const int ww, const std::string header="", const std::string end="\n", const bool use_coutCBL=true, std::ostream &stream=std::cout, const std::string colour=cbl::par::col_default)
function to print values with a proper homegenised format
Definition: Kernel.h:1142
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::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