48 const double FF = alpha_par/alpha_perp;
49 const double fact = sqrt(1.+mu*mu*(1./(FF*FF)-1.));
51 return {kk/alpha_perp*fact, mu/FF/fact};
57 double cbl::modelling::twopt::Pkmu_DeWiggled (
const double kk,
const double mu,
const double sigmaNL_perp,
const double sigmaNL_par,
const double linear_growth_rate,
const double bias,
const double SigmaS,
const std::shared_ptr<cbl::glob::FuncGrid> Pk,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_NW)
59 const double beta = linear_growth_rate/
bias;
61 const double KaiserBoost = pow(1.+mu*mu*beta, 2);
62 const double sigmaNL2 = 0.5*((1.-mu*mu)*sigmaNL_perp*sigmaNL_perp+mu*mu*sigmaNL_par*sigmaNL_par);
64 const double Fstreaming = pow(1+kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*SigmaS*SigmaS, -2);
66 const double sigmaNL = sqrt(sigmaNL_perp*sigmaNL_perp+sigmaNL_par*sigmaNL_par);
67 const double _Pk = (sigmaNL<1.e-5) ? Pk->operator()(kk) : (Pk->operator()(kk)-Pk_NW->operator()(kk))*exp(-kk*kk*sigmaNL2)+Pk_NW->operator()(kk);
69 return bias*
bias*KaiserBoost*Fstreaming*_Pk;
76 double cbl::modelling::twopt::Pkmu_ModeCoupling (
const double kk,
const double mu,
const double linear_growth_rate,
const double bias,
const double SigmaV,
const double AMC,
const std::shared_ptr<cbl::glob::FuncGrid> Pk,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_1loop)
78 const double beta = linear_growth_rate/
bias;
80 const double KaiserBoost = pow(1.+mu*mu*beta, 2);
81 const double Fstreaming = pow(1+kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*SigmaV*SigmaV, -2);
82 double Pk_NL =
bias*
bias*(Pk->operator()(kk)*exp(-kk*kk*SigmaV*SigmaV));
87 return KaiserBoost*Fstreaming*Pk_NL;
94 double cbl::modelling::twopt::Pkmu_dispersion (
const double kk,
const double mu,
const std::string DFoG,
const double linear_growth_rate,
const double bias,
const double sigmav,
const std::shared_ptr<cbl::glob::FuncGrid> Pklin)
99 DispFactor = exp(-kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav);
101 else if (DFoG==
"Lorentzian")
102 DispFactor = 1./(1.+(kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav));
105 ErrorCBL(
"the chosen DFoG ("+DFoG+
") is not currently implemented!",
"Pkmu_dispersion",
"ModelFunction_TwoPointCorrelation.cpp");
107 return DispFactor*
bias*
bias*pow(1+linear_growth_rate/
bias*mu*mu, 2)*Pklin->operator()(kk);
113 double cbl::modelling::twopt::Pkmu_Scoccimarro (
const double kk,
const double mu,
const std::string DFoG,
const double linear_growth_rate,
const double bias,
const double sigmav,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_DeltaDelta,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_DeltaTheta,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_ThetaTheta)
115 const double Pk_dd = Pk_DeltaDelta->operator()(kk);
116 const double Pk_dt = Pk_DeltaTheta->operator()(kk);
117 const double Pk_tt = Pk_ThetaTheta->operator()(kk);
121 if (DFoG==
"Gaussian")
122 DispFactor = exp(-kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav);
124 else if (DFoG==
"Lorentzian")
125 DispFactor = 1./(1.+(kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav));
128 ErrorCBL(
"the chosen DFoG ("+DFoG+
") is not currently implemented!",
"Pkmu_Scoccimarro",
"ModelFunction_TwoPointCorrelation.cpp");
130 const double beta = linear_growth_rate/
bias;
132 return DispFactor*
bias*
bias*(Pk_dd+2.*beta*pow(mu, 2)*Pk_dt+pow(beta, 2)*pow(mu, 4)*Pk_tt);
139 double cbl::modelling::twopt::Pkmu_Scoccimarro_fitPezzotta (
const double kk,
const double mu,
const std::string DFoG,
const double linear_growth_rate,
const double bias,
const double sigmav,
const double kd,
const double kt,
const std::shared_ptr<cbl::glob::FuncGrid> Pklin,
const std::shared_ptr<cbl::glob::FuncGrid> Pknonlin)
141 const double beta = linear_growth_rate/
bias;
143 double Pk_dd = Pknonlin->operator()(kk);
144 double Pk_dt = sqrt(Pknonlin->operator()(kk)*Pklin->operator()(kk)*exp(-kk/kd));
145 double Pk_tt = Pklin->operator()(kk)*exp(-kk/kt);
147 if (DFoG==
"Gaussian")
148 DispFactor = exp(-kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav);
150 else if (DFoG==
"Lorentzian")
151 DispFactor = 1./(1.+(kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav));
154 ErrorCBL(
"the chosen DFoG ("+DFoG+
") is not currently implemented!",
"Pkmu_Scoccimarro_fitPezzotta",
"ModelFunction_TwoPointCorrelation.cpp");
156 return DispFactor*
bias*
bias*(Pk_dd+2.*beta*pow(mu, 2)*Pk_dt+pow(beta, 2)*pow(mu, 4)*Pk_tt);
163 double cbl::modelling::twopt::Pkmu_Scoccimarro_fitBel (
const double kk,
const double mu,
const std::string DFoG,
const double linear_growth_rate,
const double bias,
const double sigmav,
const double kd,
const double bb,
const double a1,
const double a2,
const double a3,
const std::shared_ptr<cbl::glob::FuncGrid> Pklin,
const std::shared_ptr<cbl::glob::FuncGrid> Pknonlin)
165 const double beta = linear_growth_rate/
bias;
167 double Pk_dd = Pknonlin->operator()(kk);
168 double Pk_dt = sqrt(Pknonlin->operator()(kk)*Pklin->operator()(kk))*exp(-kk/kd-bb*pow(kk, 6.0));
169 double Pk_tt = Pklin->operator()(kk)*exp(-kk*(a1 + a2*kk + a3*kk*kk));
171 if (DFoG==
"Gaussian")
172 DispFactor = exp(-kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav);
173 else if (DFoG==
"Lorentzian")
174 DispFactor = 1./(1.+(kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav));
176 ErrorCBL(
"the chosen DFoG ("+DFoG+
") is not currently implemented!",
"Pkmu_Scoccimarro_fitBel",
"ModelFunction_TwoPointCorrelation.cpp");
178 return DispFactor*
bias*
bias*(Pk_dd+2.*beta*pow(mu, 2)*Pk_dt+pow(beta, 2)*pow(mu, 4)*Pk_tt);
185 double cbl::modelling::twopt::Pkmu_TNS (
const double kk,
const double mu,
const std::string DFoG,
const double linear_growth_rate,
const double bias,
const double sigmav,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_DeltaDelta,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_DeltaTheta,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_ThetaTheta,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A11,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A12,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A22,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A23,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A33,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B12,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B13,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B14,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B22,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B23,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B24,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B33,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B34,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B44)
187 const double beta = linear_growth_rate/
bias;
189 const double Pk_dd = Pk_DeltaDelta->operator()(kk);
190 const double Pk_dt = Pk_DeltaTheta->operator()(kk);
191 const double Pk_tt = Pk_ThetaTheta->operator()(kk);
193 const double pkA11 = Pk_A11->operator()(kk);
194 const double pkA12 = Pk_A12->operator()(kk);
195 const double pkA22 = Pk_A22->operator()(kk);
196 const double pkA23 = Pk_A23->operator()(kk);
197 const double pkA33 = Pk_A33->operator()(kk);
198 const double pkB12 = Pk_B12->operator()(kk);
199 const double pkB13 = Pk_B13->operator()(kk);
200 const double pkB14 = Pk_B14->operator()(kk);
201 const double pkB22 = Pk_B22->operator()(kk);
202 const double pkB23 = Pk_B23->operator()(kk);
203 const double pkB24 = Pk_B24->operator()(kk);
204 const double pkB33 = Pk_B33->operator()(kk);
205 const double pkB34 = Pk_B34->operator()(kk);
206 const double pkB44 = Pk_B44->operator()(kk);
209 if (DFoG==
"Gaussian")
210 DispFactor = exp(-kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav);
211 else if (DFoG==
"Lorentzian")
212 DispFactor = 1./(1.+(kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav));
214 ErrorCBL(
"the chosen DFoG ("+DFoG+
") is not currently implemented!",
"Pkmu_TNS",
"ModelFunction_TwoPointCorrelation.cpp");
216 const double Pk_Scoccimarro_large_scales =
bias*
bias*(Pk_dd+2*beta*pow(mu, 2)*Pk_dt+pow(beta, 2)*pow(mu, 4)*Pk_tt);
218 const double A2 = pow(mu, 2)*(beta*pkA11+pow(beta, 2)*pkA12);
219 const double A4 = pow(mu, 4)*pow(beta, 2)*(pkA22+beta*pkA23);
220 const double A6 = pow(mu, 6)*pow(beta, 3)*pkA33;
221 const double B2 = pow(mu, 2)*(pow(beta, 2)*pkB12+pow(beta, 3)*pkB13+pow(beta, 4)*pkB14);
222 const double B4 = pow(mu, 4)*(pow(beta, 2)*pkB22+pow(beta, 3)*pkB23+pow(beta, 4)*pkB24);
223 const double B6 = pow(mu, 6)*(pow(beta, 3)*pkB33+pow(beta, 4)*pkB34);
224 const double B8 = pow(mu, 8)*pow(beta, 4)*pkB44;
226 const double Pk_A = pow(
bias, 3)*(A2+A4+A6);
227 const double Pk_B = pow(
bias, 4)*(B2+B4+B6+B8);
229 return DispFactor*(Pk_Scoccimarro_large_scales+Pk_A+Pk_B);
236 double cbl::modelling::twopt::Pkmu_eTNS (
const double kk,
const double mu,
const std::string DFoG,
const double linear_growth_rate,
const double bias,
const double bias2,
const double sigmav,
const double Ncorr,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_DeltaDelta,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_DeltaTheta,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_ThetaTheta,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A11,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A12,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A22,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A23,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_A33,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B12,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B13,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B14,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B22,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B23,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B24,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B33,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B34,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_B44,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_b2d,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_b2v,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_b22,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_bs2d,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_bs2v,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_b2s2,
const std::shared_ptr<cbl::glob::FuncGrid> Pk_bs22,
const std::shared_ptr<cbl::glob::FuncGrid> sigma32Pklin)
238 const double beta = linear_growth_rate/
bias;
240 const double pk_dd = Pk_DeltaDelta->operator()(kk);
241 const double pk_dt = Pk_DeltaTheta->operator()(kk);
242 const double pk_tt = Pk_ThetaTheta->operator()(kk);
244 const double pkA11 = Pk_A11->operator()(kk);
245 const double pkA12 = Pk_A12->operator()(kk);
246 const double pkA22 = Pk_A22->operator()(kk);
247 const double pkA23 = Pk_A23->operator()(kk);
248 const double pkA33 = Pk_A33->operator()(kk);
249 const double pkB12 = Pk_B12->operator()(kk);
250 const double pkB13 = Pk_B13->operator()(kk);
251 const double pkB14 = Pk_B14->operator()(kk);
252 const double pkB22 = Pk_B22->operator()(kk);
253 const double pkB23 = Pk_B23->operator()(kk);
254 const double pkB24 = Pk_B24->operator()(kk);
255 const double pkB33 = Pk_B33->operator()(kk);
256 const double pkB34 = Pk_B34->operator()(kk);
257 const double pkB44 = Pk_B44->operator()(kk);
259 const double pk_b2d = Pk_b2d->operator()(kk);
260 const double pk_b2v = Pk_b2v->operator()(kk);
261 const double pk_b22 = Pk_b22->operator()(kk);
262 const double pk_bs2d = Pk_bs2d->operator()(kk);
263 const double pk_bs2v = Pk_bs2v->operator()(kk);
264 const double pk_b2s2 = Pk_b2s2->operator()(kk);
265 const double pk_bs22 = Pk_bs22->operator()(kk);
266 const double sigma32pklin = sigma32Pklin->operator()(kk);
269 if (DFoG==
"Gaussian")
270 DispFactor = exp(-kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav);
271 else if (DFoG==
"Lorentzian")
272 DispFactor = 1./(1.+(kk*kk*mu*mu*linear_growth_rate*linear_growth_rate*sigmav*sigmav));
274 ErrorCBL(
"the chosen DFoG ("+DFoG+
") is not currently implemented!",
"Pkmu_TNS",
"ModelFunction_TwoPointCorrelation.cpp");
276 const double A2 = pow(mu, 2)*(beta*pkA11+pow(beta, 2)*pkA12);
277 const double A4 = pow(mu, 4)*pow(beta, 2)*(pkA22+beta*pkA23);
278 const double A6 = pow(mu, 6)*pow(beta, 3)*pkA33;
279 const double B2 = pow(mu, 2)*(pow(beta, 2)*pkB12+pow(beta, 3)*pkB13+pow(beta, 4)*pkB14);
280 const double B4 = pow(mu, 4)*(pow(beta, 2)*pkB22+pow(beta, 3)*pkB23+pow(beta, 4)*pkB24);
281 const double B6 = pow(mu, 6)*(pow(beta, 3)*pkB33+pow(beta, 4)*pkB34);
282 const double B8 = pow(mu, 8)*pow(beta, 4)*pkB44;
284 const double Pk_A = pow(
bias, 3)*(A2+A4+A6);
285 const double Pk_B = pow(
bias, 4)*(B2+B4+B6+B8);
288 const double bs2 = (-4./7.)*(
bias-1.);
289 const double b3nl = (32./315.)*(
bias-1.);
291 const double Pgdd =
bias*
bias*pk_dd+2*bias2*
bias*pk_b2d+2*bs2*
bias*pk_bs2d+2*b3nl*
bias*sigma32pklin+bias2*bias2*pk_b22+2*bias2*bs2*pk_b2s2+bs2*bs2*pk_bs22+Ncorr;
292 const double Pgdv =
bias*pk_dt+bias2*pk_b2v+bs2*pk_bs2v+b3nl*sigma32pklin;
293 const double Pk_Scoccimarro_large_scales = Pgdd+2.*linear_growth_rate*pow(mu, 2)*Pgdv+pow(linear_growth_rate, 2)*pow(mu, 4)*pk_tt;
295 return DispFactor*(Pk_Scoccimarro_large_scales+Pk_A+Pk_B);
302 double cbl::modelling::twopt::Pkmu (
const double kk,
const double mu,
const std::string model,
const std::vector<double> parameter,
const std::vector<std::shared_ptr<glob::FuncGrid>> pk_interp,
const double alpha_perp,
const double alpha_par)
304 const double kp =
true_k_mu_AP(kk, mu, alpha_perp, alpha_par)[0];
305 const double mup =
true_k_mu_AP(kk, mu, alpha_perp, alpha_par)[1];
309 if (model==
"dispersion_dewiggled") {
310 if (parameter.size()!=5)
311 ErrorCBL(
"the "+model+
" model has 5 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
312 pkmu =
Pkmu_DeWiggled(kp, mup, parameter[0], parameter[1], parameter[2], parameter[3], parameter[4], pk_interp[0], pk_interp[1]);
315 else if (model==
"dispersion_modecoupling") {
316 if (parameter.size()!=4)
317 ErrorCBL(
"the "+model+
" model has 4 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
318 pkmu =
Pkmu_ModeCoupling(kp, mup, parameter[0], parameter[1], parameter[2], parameter[3], pk_interp[0], pk_interp[1]);
321 else if (model==
"dispersion_Gauss") {
322 if (parameter.size()!=3)
323 ErrorCBL(
"the "+model+
" model has 3 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
324 pkmu =
Pkmu_dispersion(kp, mup,
"Gaussian", parameter[0], parameter[1], parameter[2], pk_interp[0]);
327 else if (model==
"dispersion_Lorentz") {
328 if (parameter.size()!=3)
329 ErrorCBL(
"the "+model+
" model has 3 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
330 pkmu =
Pkmu_dispersion(kp, mup,
"Lorentzian", parameter[0], parameter[1], parameter[2], pk_interp[0]);
333 else if (model==
"Scoccimarro_Gauss") {
334 if (parameter.size()!=3)
335 ErrorCBL(
"the "+model+
" model has 3 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
336 pkmu =
Pkmu_Scoccimarro(kp, mup,
"Gaussian", parameter[0], parameter[1], parameter[2], pk_interp[0], pk_interp[1], pk_interp[2]);
339 else if (model==
"Scoccimarro_Lorentz") {
340 if (parameter.size()!=3)
341 ErrorCBL(
"the "+model+
" model has 3 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
342 pkmu =
Pkmu_Scoccimarro(kp, mup,
"Lorentzian", parameter[0], parameter[1], parameter[2], pk_interp[0], pk_interp[1], pk_interp[2]);
345 else if (model==
"Scoccimarro_Pezzotta_Gauss") {
346 if (parameter.size()!=5)
347 ErrorCBL(
"the "+model+
" model has 5 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
348 pkmu =
Pkmu_Scoccimarro_fitPezzotta(kp, mup,
"Gaussian", parameter[0], parameter[1], parameter[2], parameter[3], parameter[4], pk_interp[0], pk_interp[1]);
351 else if (model==
"Scoccimarro_Pezzotta_Lorentz") {
352 if (parameter.size()!=5)
353 ErrorCBL(
"the "+model+
" model has 5 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
354 pkmu =
Pkmu_Scoccimarro_fitPezzotta(kp, mup,
"Lorentzian", parameter[0], parameter[1], parameter[2], parameter[3], parameter[4], pk_interp[0], pk_interp[1]);
357 else if (model==
"Scoccimarro_Bel_Gauss") {
358 if (parameter.size()!=8)
359 ErrorCBL(
"the "+model+
" model has 8 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
360 pkmu =
Pkmu_Scoccimarro_fitBel(kp, mup,
"Gaussian", parameter[0], parameter[1], parameter[2], parameter[3], parameter[4], parameter[5], parameter[6], parameter[7], pk_interp[0], pk_interp[1]);
363 else if (model==
"Scoccimarro_Bel_Lorentz") {
364 if (parameter.size()!=8)
365 ErrorCBL(
"the "+model+
" model has 8 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
366 pkmu =
Pkmu_Scoccimarro_fitBel(kp, mup,
"Lorentzian", parameter[0], parameter[1], parameter[2], parameter[3], parameter[4], parameter[5], parameter[6], parameter[7], pk_interp[0], pk_interp[1]);
369 else if (model==
"TNS_Gauss") {
370 if (parameter.size()!=3)
371 ErrorCBL(
"the "+model+
" model has 3 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
372 pkmu =
Pkmu_TNS(kp, mup,
"Gaussian", parameter[0], parameter[1], parameter[2], pk_interp[0], pk_interp[1], pk_interp[2], pk_interp[3], pk_interp[4], pk_interp[5], pk_interp[6], pk_interp[7], pk_interp[8], pk_interp[9], pk_interp[10], pk_interp[11], pk_interp[12], pk_interp[13], pk_interp[14], pk_interp[15], pk_interp[16]);
375 else if (model==
"TNS_Lorentz") {
376 if (parameter.size()!=3)
377 ErrorCBL(
"the "+model+
" model has 3 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
378 pkmu =
Pkmu_TNS(kp, mup,
"Lorentzian", parameter[0], parameter[1], parameter[2], pk_interp[0], pk_interp[1], pk_interp[2], pk_interp[3], pk_interp[4], pk_interp[5], pk_interp[6], pk_interp[7], pk_interp[8], pk_interp[9], pk_interp[10], pk_interp[11], pk_interp[12], pk_interp[13], pk_interp[14], pk_interp[15], pk_interp[16]);
381 else if (model==
"eTNS_Gauss") {
382 if (parameter.size()!=5)
383 ErrorCBL(
"the "+model+
" model has 5 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
384 pkmu =
Pkmu_eTNS(kp, mup,
"Gaussian", parameter[0], parameter[1], parameter[2], parameter[3], parameter[4], pk_interp[0], pk_interp[1], pk_interp[2], pk_interp[3], pk_interp[4], pk_interp[5], pk_interp[6], pk_interp[7], pk_interp[8], pk_interp[9], pk_interp[10], pk_interp[11], pk_interp[12], pk_interp[13], pk_interp[14], pk_interp[15], pk_interp[16], pk_interp[17], pk_interp[18], pk_interp[19], pk_interp[20], pk_interp[21], pk_interp[22], pk_interp[23], pk_interp[24]);
387 else if (model==
"eTNS_Lorentz") {
388 if (parameter.size()!=5)
389 ErrorCBL(
"the "+model+
" model has 5 parameters, while in the parameter vector in input has "+
conv(parameter.size(),
par::fINT)+
" parameters!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
390 pkmu =
Pkmu_eTNS(kp, mup,
"Lorentzian", parameter[0], parameter[1], parameter[2], parameter[3], parameter[4], pk_interp[0], pk_interp[1], pk_interp[2], pk_interp[3], pk_interp[4], pk_interp[5], pk_interp[6], pk_interp[7], pk_interp[8], pk_interp[9], pk_interp[10], pk_interp[11], pk_interp[12], pk_interp[13], pk_interp[14], pk_interp[15], pk_interp[16], pk_interp[17], pk_interp[18], pk_interp[19], pk_interp[20], pk_interp[21], pk_interp[22], pk_interp[23], pk_interp[24]);
394 ErrorCBL(
"the chosen model ("+model+
") is not currently implemented!",
"Pkmu",
"ModelFunction_TwoPointCorrelation.cpp");
396 return pkmu/(alpha_perp*alpha_perp*alpha_par);
403 double cbl::modelling::twopt::Pk_l (
const double kk,
const int l,
const std::string model,
const std::vector<double> parameter,
const std::vector<std::shared_ptr<glob::FuncGrid>> pk_interp,
const double prec,
const double alpha_perp,
const double alpha_par)
405 auto integrand = [&] (
const double mu)
416 std::vector<double>
cbl::modelling::twopt::Pk_l (
const std::vector<double> kk,
const int l,
const std::string model,
const std::vector<double> parameter,
const std::vector<std::shared_ptr<glob::FuncGrid>> pk_interp,
const double prec,
const double alpha_perp,
const double alpha_par)
418 vector<double> Pkl(kk.size(), 0);
420 for (
size_t i=0; i<kk.size(); i++) {
422 auto integrand = [&] (
const double mu)
437 cbl::glob::FuncGrid cbl::modelling::twopt::Xil_interp (
const std::vector<double> kk,
const int l,
const std::string model,
const std::vector<double> parameter,
const std::vector<std::shared_ptr<glob::FuncGrid>> pk_interp,
const double prec,
const double alpha_perp,
const double alpha_par)
439 vector<double> Pkl =
Pk_l(kk, l, model, parameter, pk_interp, prec, alpha_perp, alpha_par);
440 vector<double> rr, Xil;
453 std::vector<std::vector<double>>
cbl::modelling::twopt::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,
const double alpha_perp,
const double alpha_par)
455 vector<vector<double>> Xil(3);
457 for (
int i=0; i<nmultipoles; i++) {
458 double sign = (i%2==0) ? 1 : -1;
459 Xil[i] =
Xil_interp(pk_interp[0]->x(), 2*i, model, parameter, pk_interp, prec, alpha_perp, alpha_par).
eval_func(rr);
460 for (
size_t j=0; j<rr.size(); j++)
471 std::vector<double>
cbl::modelling::twopt::Xi_l (
const std::vector<double> rr,
const std::vector<int> dataset_order,
const std::vector<bool> use_pole,
const std::string model,
const std::vector<double> parameter,
const std::vector<std::shared_ptr<glob::FuncGrid>> pk_interp,
const double prec,
const double alpha_perp,
const double alpha_par)
473 vector<cbl::glob::FuncGrid> interp_Xil(3);
474 vector<double> sign={1., -1., 1.};
476 for (
size_t i=0; i<3; i++)
478 interp_Xil[i] =
Xil_interp(pk_interp[0]->x(), 2*i, model, parameter, pk_interp, prec, alpha_perp, alpha_par);
480 vector<double> Xil(rr.size());
482 for (
size_t i=0; i<rr.size(); i++)
483 Xil[i] = sign[dataset_order[i]]*interp_Xil[dataset_order[i]](rr[i]);
492 double cbl::modelling::twopt::Xi_polar(
const double rad_fid,
const double mu_fid,
const double alpha_perpendicular,
const double alpha_parallel,
const std::vector<std::shared_ptr<cbl::glob::FuncGrid>> xi_multipoles)
494 const double apar2 = alpha_parallel*alpha_parallel;
495 const double aperp2 = alpha_perpendicular*alpha_perpendicular;
497 double mu_fid_sq = mu_fid*mu_fid;
498 const double factor = sqrt(apar2*mu_fid_sq+(1.-mu_fid_sq)*aperp2);
499 const double mu_true = mu_fid * alpha_parallel/factor;
500 const double mu_true_sq = mu_true*mu_true;
501 const double st = rad_fid*factor;
503 return xi_multipoles[0]->operator()(st)+
504 xi_multipoles[1]->
operator()(st)*0.5*(3*mu_true_sq-1)+
505 xi_multipoles[2]->operator()(st)*0.125*(35*mu_true_sq*mu_true_sq-30*mu_true_sq+3);
512 std::vector<std::vector<double>>
cbl::modelling::twopt::Xi_rppi (
const std::vector<double> rp,
const std::vector<double>
pi,
const std::string model,
const std::vector<double> parameter,
const std::vector<std::shared_ptr<glob::FuncGrid>> pk_interp,
const double prec,
const double alpha_perp,
const double alpha_par)
515 const size_t nmultipoles = pk_interp.size();
516 vector<double> sign = {1., -1., 1., 1, -1};
519 vector<glob::FuncGrid> interp_Xil(nmultipoles);
521 for (
size_t i=0; i<nmultipoles; i++)
522 interp_Xil[i] =
Xil_interp(pk_interp[0]->x(), 2*i, model, parameter, pk_interp, prec, alpha_perp, alpha_par);
524 vector<vector<double>> xi_rppi(rp.size(), vector<double>(
pi.size(), 0));
526 for (
size_t i =0; i<rp.size(); i++)
527 for (
size_t j =0; j<
pi.size(); j++) {
528 double s = sqrt(rp[i]*rp[i]+
pi[j]*
pi[j]);
530 for (
size_t l=0; l<nmultipoles; l++)
541 std::vector<double>
cbl::modelling::twopt::wp_from_Xi_rppi (
const std::vector<double> rp,
const double pimax,
const std::string model,
const std::vector<double> parameter,
const std::vector<std::shared_ptr<glob::FuncGrid>> pk_interp,
const double prec,
const double alpha_perp,
const double alpha_par) {
546 vector<double>
wp(rp.size());
548 for (
size_t i=0; i<rp.size(); i++) {
560 std::vector<std::vector<double>>
cbl::modelling::twopt::damped_Pk_terms (
const std::vector<double> kk,
const double linear_growth_rate,
const double SigmaS,
const std::shared_ptr<cbl::glob::FuncGrid> PkDM)
563 double sqrt_pi = sqrt(
par::pi);
564 vector<vector<double>> pk(3, vector<double>(kk.size(), 0));
566 for (
size_t i=0; i< kk.size(); i++) {
567 double kSigma = kk[i]*SigmaS;
568 double pk_dm = PkDM->operator()(kk[i]);
569 double erf_kSigma = erf(kSigma);
571 pk[0][i] = pk_dm*sqrt_pi/(2.*kSigma)*erf_kSigma;
572 pk[1][i] = linear_growth_rate*pow(kSigma, -3.)*pk_dm*(0.5*sqrt_pi*erf_kSigma-kSigma*exp(-kSigma*kSigma));
573 pk[2][i] = linear_growth_rate*linear_growth_rate*pow(kSigma, -5.)*pk_dm*(3.*sqrt_pi/8.*erf_kSigma-0.25*kSigma*(2*kSigma*kSigma+3)*exp(-kSigma*kSigma));
583 std::vector<double>
cbl::modelling::twopt::damped_Xi (
const std::vector<double> ss,
const double bias,
const double linear_growth_rate,
const double SigmaS,
const std::vector<double> kk,
const std::shared_ptr<cbl::glob::FuncGrid> PkDM)
587 vector<double> xi(ss.size(), 0);
589 for (
size_t i=0; i<pk_terms.size(); i++) {
591 for (
size_t j=0; j<ss.size(); j++)
592 xi[j] += pow(
bias, 2-i)*xi_term[j];
Global functions to model two-point correlation functions of any type.
std::vector< double > eval_func(const std::vector< double > xx) const
evaluate the function at the xx points
double integrate_qag(const double a, const double b, const double rel_err=1.e-2, const double abs_err=1.e-6, const int limit_size=1000, const int rule=6)
compute the definite integral with GSL qag method
static const char fINT[]
conversion symbol for: int -> std::string
static const double pi
: the ratio of a circle's circumference to its diameter
std::vector< double > damped_Xi(const std::vector< double > ss, const double bias, const double linear_growth_rate, const double SigmaS, const std::vector< double > kk, const std::shared_ptr< cbl::glob::FuncGrid > PkDM)
the damped two-point correlation monopole; from Sereno et al. 2015
std::vector< double > wp_from_Xi_rppi(const std::vector< double > rp, const double pimax, 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 projected two-point correlation function
double Pkmu_Scoccimarro(const double kk, const double mu, const std::string DFoG, const double linear_growth_rate, const double bias, const double sigmav, const std::shared_ptr< cbl::glob::FuncGrid > Pk_DeltaDelta, const std::shared_ptr< cbl::glob::FuncGrid > Pk_DeltaTheta, const std::shared_ptr< cbl::glob::FuncGrid > Pk_ThetaTheta)
the redshift-space galaxy power spectrum, as a function of and , predicted by the Scoccimarro model
double Pkmu_DeWiggled(const double kk, const double mu, const double sigmaNL_perp, const double sigmaNL_par, const double linear_growth_rate, const double bias, const double SigmaS, const std::shared_ptr< cbl::glob::FuncGrid > Pk, const std::shared_ptr< cbl::glob::FuncGrid > Pk_NW)
the redshift-space galaxy power spectrum, as a function of and , predicted by the de-wiggled model
double Pkmu_dispersion(const double kk, const double mu, const std::string DFoG, const double linear_growth_rate, const double bias, const double sigmav, const std::shared_ptr< cbl::glob::FuncGrid > Pklin)
the redshift-space galaxy power spectrum, as a function of and , predicted by the dispersion model
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
double Pkmu_eTNS(const double kk, const double mu, const std::string DFoG, const double linear_growth_rate, const double bias, const double bias2, const double sigmav, const double Ncorr, const std::shared_ptr< cbl::glob::FuncGrid > Pk_DeltaDelta, const std::shared_ptr< cbl::glob::FuncGrid > Pk_DeltaTheta, const std::shared_ptr< cbl::glob::FuncGrid > Pk_ThetaTheta, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A11, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A12, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A22, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A23, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A33, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B12, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B13, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B14, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B22, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B23, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B24, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B33, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B34, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B44, const std::shared_ptr< cbl::glob::FuncGrid > Pk_b2d, const std::shared_ptr< cbl::glob::FuncGrid > Pk_b2v, const std::shared_ptr< cbl::glob::FuncGrid > Pk_b22, const std::shared_ptr< cbl::glob::FuncGrid > Pk_bs2d, const std::shared_ptr< cbl::glob::FuncGrid > Pk_bs2v, const std::shared_ptr< cbl::glob::FuncGrid > Pk_b2s2, const std::shared_ptr< cbl::glob::FuncGrid > Pk_bs22, const std::shared_ptr< cbl::glob::FuncGrid > sigma32Pklin)
the redshift-space galaxy power spectrum, as a function of and , predicted by the extended TNS (Taru...
double Pkmu_Scoccimarro_fitBel(const double kk, const double mu, const std::string DFoG, const double linear_growth_rate, const double bias, const double sigmav, const double kd, const double bb, const double a1, const double a2, const double a3, const std::shared_ptr< cbl::glob::FuncGrid > Pklin, const std::shared_ptr< cbl::glob::FuncGrid > Pknonlin)
the redshift-space galaxy power spectrum, as a function of and , predicted by the Scoccimarro model
std::vector< std::vector< double > > Xi_rppi(const std::vector< double > rp, const std::vector< double > pi, 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 cartesian two-point correlation function
std::vector< std::vector< double > > damped_Pk_terms(const std::vector< double > kk, const double linear_growth_rate, const double SigmaS, const std::shared_ptr< cbl::glob::FuncGrid > PkDM)
the power spectrum terms obtained integrating the redshift space 2D power spectrum
cbl::glob::FuncGrid Xil_interp(const std::vector< double > kk, const int l, 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 interpolating function of multipole expansion of the two-point correlation function at a given or...
double Pkmu_ModeCoupling(const double kk, const double mu, const double linear_growth_rate, const double bias, const double sigmav, const double AMC, const std::shared_ptr< cbl::glob::FuncGrid > PkLin, const std::shared_ptr< cbl::glob::FuncGrid > PkMC)
the redshift-space galaxy power spectrum, as a function of and , predicted by the mode-coupling mode...
double Pkmu(const double kk, const double mu, const std::string model, const std::vector< double > parameter, const std::vector< std::shared_ptr< glob::FuncGrid >> pk_interp, const double alpha_perp=1., const double alpha_par=1.)
the power spectrum as a function of k and
double Xi_polar(const double rad_fid, const double mu_fid, const double alpha_perpendicular, const double alpha_parallel, const std::vector< std::shared_ptr< cbl::glob::FuncGrid >> xi_multipoles)
the polar two-point correlation function
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
double Pkmu_Scoccimarro_fitPezzotta(const double kk, const double mu, const std::string DFoG, const double linear_growth_rate, const double bias, const double sigmav, const double kd, const double kt, const std::shared_ptr< cbl::glob::FuncGrid > Pklin, const std::shared_ptr< cbl::glob::FuncGrid > Pknonlin)
the redshift-space galaxy power spectrum, as a function of and , predicted by the Scoccimarro model
std::vector< double > true_k_mu_AP(const double kk, const double mu, const double alpha_perp, const double alpha_par)
true k and power spectrum coordinates as a function of observed ones
double Pk_l(const double kk, const int l, 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 power spectrum
double Pkmu_TNS(const double kk, const double mu, const std::string DFoG, const double linear_growth_rate, const double bias, const double sigmav, const std::shared_ptr< cbl::glob::FuncGrid > Pk_DeltaDelta, const std::shared_ptr< cbl::glob::FuncGrid > Pk_DeltaTheta, const std::shared_ptr< cbl::glob::FuncGrid > Pk_ThetaTheta, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A11, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A12, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A22, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A23, const std::shared_ptr< cbl::glob::FuncGrid > Pk_A33, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B12, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B13, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B14, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B22, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B23, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B24, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B33, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B34, const std::shared_ptr< cbl::glob::FuncGrid > Pk_B44)
the redshift-space galaxy power spectrum, as a function of and , predicted by the TNS (Taruya,...
std::vector< double > transform_FFTlog(const std::vector< double > yy, const int dir, const std::vector< double > xx, const std::vector< double > fx, const double mu=0, const double q=0, const double kr=1, const int kropt=0)
wrapper of the FFTlog to compute the discrete Fourier sine or cosine transform of a logarithmically...
double GSL_integrate_cquad(gsl_function func, const double a, const double b, const double rel_err=1.e-3, const double abs_err=0, const int nevals=100)
integral, using the gsl cquad method
The global namespace of the CosmoBolognaLib
T Min(const std::vector< T > vect)
minimum element of a std::vector
std::string conv(const T val, const char *fact)
convert a number to a std::string
std::vector< T > linear_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with linearly spaced values
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
T Max(const std::vector< T > vect)
maximum element of a std::vector
double legendre_polynomial(const double mu, const int l)
the order l Legendre polynomial
double wp(const double rp, const std::vector< double > rr, const std::vector< double > xi, const double r_max=100.)
the projected two-point correlation function