Skip to content

2 filters for kSZ #10

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
olakusiak wants to merge 3 commits into
base: master
Choose a base branch
from
Open

2 filters for kSZ #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f7394ba
2 filters for kSZ
olakusiak Nov 25, 2025
e5222e6
forgot this file
olakusiak Nov 25, 2025
3d9f52e
2-filter for ksz corrected
olakusiak May 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions class-sz/include/class_sz.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1594,6 +1594,10 @@ double szcounts_ntot;
double * f_unwise_filter;
int unwise_filter_size;

double * l_unwise_filter2;
double * f_unwise_filter2;
int unwise_filter2_size;

double * M_min_of_z_z;
double * M_min_of_z_M_min;
int M_min_of_z_size;
Expand Down Expand Up @@ -2497,6 +2501,8 @@ double get_tau_profile_at_z_m_l(double z,

double get_ksz_filter_at_l(double l,
struct class_sz_structure * pclass_sz);
double get_ksz_filter_at_l2(double l,
struct class_sz_structure * pclass_sz);

double get_M_min_of_z(double l,
struct class_sz_structure * pclass_sz);
Expand Down
1 change: 1 addition & 0 deletions class-sz/include/class_sz_precisions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ class_sz_string_parameter(cib_Snu_file_nu,"/class_sz_auxiliary_files/includes/fi


class_sz_string_parameter(ksz_filter_file,"/class_sz_auxiliary_files/UNWISE_galaxy_distributions/unwise_filter_functions_l_fl.txt","ksz filter file")
class_sz_string_parameter(ksz_filter_file2,"/class_sz_auxiliary_files/UNWISE_galaxy_distributions/unwise_filter_functions_l_fl2.txt","ksz filter file2")
class_sz_string_parameter(ksz_template_file,"/class_sz_auxiliary_files/cl_ksz_bat.dat","ksz template file")
class_sz_string_parameter(ksz_reio_template_file,"/class_sz_auxiliary_files/FBN_kSZ_PS_patchy.d.txt","ksz template file, reio contribution")

Expand Down
1 change: 1 addition & 0 deletions class-sz/include/class_sz_tools.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -315,6 +315,7 @@ double delta_to_delta_prime_nfw(
int load_unbinned_nl_yy(struct class_sz_structure * pclass_sz);
int load_unbinned_nl_tt(struct class_sz_structure * pclass_sz);
int load_ksz_filter(struct class_sz_structure * pclass_sz);
int load_ksz_filter2(struct class_sz_structure * pclass_sz);
int load_M_min_of_z(struct class_sz_structure * pclass_sz);
double get_T10_alpha_at_z(double z_asked, struct class_sz_structure * pclass_sz);
double get_knl_at_z(double z_asked, struct class_sz_structure * pclass_sz);
Expand Down
127 changes: 86 additions & 41 deletions class-sz/source/class_sz.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1133,7 +1133,7 @@ if (pclass_sz->has_kSZ_kSZ_gal_1h
|| pclass_sz->has_kSZ_kSZ_gal_3h
|| pclass_sz->has_kSZ_kSZ_gal_hf)
load_ksz_filter(pclass_sz);

load_ksz_filter2(pclass_sz);

if (pclass_sz->has_tSZ_gal_1h
|| pclass_sz->has_tSZ_gal_2h
Expand Down Expand Up @@ -2194,6 +2194,7 @@ double cl_tt_lensed_fft[N];
double cl_ksz_fft[N];
double cl_noise_fft[N];
double Fl_bl;
double Fl_bl2;
double cl_ttf[N];
int i;
double * cl_tot;
Expand Down Expand Up @@ -2238,15 +2239,15 @@ cl_noise_fft[i] = pwl_value_1d(pclass_sz->unbinned_nl_tt_size,pclass_sz->unbinne
cl_ksz_fft[i] = cl_ksz_fft[i]/l[i]/(l[i]+1.)*2.*_PI_*1e-12/pba->T_cmb/pba->T_cmb;
cl_noise_fft[i] = cl_noise_fft[i]*1e-12/pba->T_cmb/pba->T_cmb;
Fl_bl = get_ksz_filter_at_l(l[i],pclass_sz);

Fl_bl2 = get_ksz_filter_at_l2(l[i],pclass_sz);

double result_ttf;
if (pclass_sz->compute_ksz2ksz2==1){
result_ttf = Fl_bl*Fl_bl*(cl_ksz_fft[i])*sqrt(2./(2.*_PI_)/(2.*_PI_)); // divide by 2pi^2 and multiply by 2 -- see formula
result_ttf = Fl_bl*Fl_bl2*(cl_ksz_fft[i])*sqrt(2./(2.*_PI_)/(2.*_PI_)); // divide by 2pi^2 and multiply by 2 -- see formula

}
else{
result_ttf = Fl_bl*Fl_bl*(cl_tt_lensed_fft[i]+cl_ksz_fft[i]+cl_noise_fft[i])*sqrt(2./(2.*_PI_)/(2.*_PI_)); // divide by 2pi^2 and multiply by 2 -- see formula
result_ttf = Fl_bl*Fl_bl2*(cl_tt_lensed_fft[i]+cl_ksz_fft[i]+cl_noise_fft[i])*sqrt(2./(2.*_PI_)/(2.*_PI_)); // divide by 2pi^2 and multiply by 2 -- see formula
}

if (l[i]<l_lensed[0]){
Expand Down Expand Up @@ -2512,9 +2513,16 @@ double lprime = exp(pclass_sz->ell_kSZ2_gal_multipole_grid[index_lprime]);
double theta = pclass_sz->theta_kSZ2_gal_theta_grid[index_theta];
double l = pclass_sz->ell[i];

double flprime = get_ksz_filter_at_l(lprime,pclass_sz);
double abs_l_plus_lprime = sqrt(l*l+lprime*lprime+2.*l*lprime*cos(theta));
double fl_plus_lprime = get_ksz_filter_at_l(abs_l_plus_lprime,pclass_sz);
/* Symmetrise over the two filters: average both filter assignments so that
swapping ksz_filter_file and ksz_filter_file2 gives identical results. */
// See Ferraro et al, Eqs.17-20: "while the second and third terms are the lowest
// order CMB lensing contribution and are equal in magnitude by symmetry."--> this assumes filters are the same
double f1_lprime = get_ksz_filter_at_l (lprime, pclass_sz);
double f2_l_plus_lprime = get_ksz_filter_at_l2(abs_l_plus_lprime,pclass_sz);
double f2_lprime = get_ksz_filter_at_l2(lprime, pclass_sz);
double f1_l_plus_lprime = get_ksz_filter_at_l (abs_l_plus_lprime,pclass_sz);
double filter_sym = 0.5*(f1_lprime*f2_l_plus_lprime + f2_lprime*f1_l_plus_lprime);

double clprime_tt_unlensed = pwl_value_1d(ppt->l_scalar_max-2,l_unlensed,cl_tt_unlensed,lprime);
if (lprime<l_unlensed[0])
Expand All @@ -2523,7 +2531,7 @@ if (lprime>l_unlensed[ppt->l_scalar_max-3])
clprime_tt_unlensed = 0.;


integrand_l_lprime_phi[index_lprime_theta] = lprime*lprime*flprime*clprime_tt_unlensed*fl_plus_lprime*cos(theta);
integrand_l_lprime_phi[index_lprime_theta] = lprime*lprime*filter_sym*clprime_tt_unlensed*cos(theta);
integrand_l_lprime_phi[index_lprime_theta] *= lprime; //for integration in log(ell)
// printf("l = %.4e lp = %.4e th = %.3e flprime = %.5e fl_plus_lprime = %.5e clprime = %.5e integrand_l_lprime_phi = %.5e\n",
// l,lprime,theta,flprime,fl_plus_lprime,clprime_tt_unlensed,
Expand Down Expand Up @@ -3284,6 +3292,8 @@ if(pclass_sz->has_kSZ_kSZ_gal_1h
// free(pclass_sz->theta_kSZ2_gal_theta_grid);
free(pclass_sz->l_unwise_filter);
free(pclass_sz->f_unwise_filter);
free(pclass_sz->l_unwise_filter2);
free(pclass_sz->f_unwise_filter2);
}
if(pclass_sz->has_kSZ_kSZ_gal_1h
|| pclass_sz->has_kSZ_kSZ_gal_2h
Expand Down Expand Up @@ -7959,6 +7969,7 @@ if ((int) pvectsz[pclass_sz->index_part_id_cov_hsv] == 6) {
double lnk[N],lnpk[N];
int ik;
double fl;
double fl2;
double taul;
double l;
double m = exp(logM);
Expand All @@ -7968,20 +7979,22 @@ if ((int) pvectsz[pclass_sz->index_part_id_cov_hsv] == 6) {
lnk[ik] = log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min));
l = k[ik];
fl = get_ksz_filter_at_l(l,pclass_sz);
fl2 = get_ksz_filter_at_l2(l,pclass_sz);

// pvectsz[pclass_sz->index_multipole_for_tau_profile] = l;
evaluate_tau_profile((l+0.5)/chi,pvecback,pvectsz,pba,pclass_sz);
taul = pvectsz[pclass_sz->index_tau_profile];
Pk[ik] = fl*taul;
Pko[ik] = fl*taul;
Pko[ik] = fl2*taul;
}

double r[N], xi[N];
double r[N], xi[N], xio[N];

xi2pk(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pko,r,xio,pclass_sz);
for (ik=0; ik<N; ik++){
// convolution:
xi[ik] = pow(xi[ik],2.);
// cross-convolution of filter1*tau and filter2*tau:
xi[ik] = xi[ik]*xio[ik];
}

pk2xi(N,r,xi,k,Pk,pclass_sz);
Expand Down Expand Up @@ -8024,6 +8037,7 @@ if ((int) pvectsz[pclass_sz->index_part_id_cov_hsv] == 6) {
double lnk[N],lnpk[N];
int ik;
double fl;
double fl2;
double taul;
double l;
double m = exp(logM);
Expand All @@ -8033,20 +8047,22 @@ if ((int) pvectsz[pclass_sz->index_part_id_cov_hsv] == 6) {
lnk[ik] = log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min));
l = k[ik];
fl = get_ksz_filter_at_l(l,pclass_sz);
fl2 = get_ksz_filter_at_l2(l,pclass_sz);

// pvectsz[pclass_sz->index_multipole_for_tau_profile] = l;
evaluate_tau_profile((l+0.5)/chi,pvecback,pvectsz,pba,pclass_sz);
taul = pvectsz[pclass_sz->index_tau_profile];
Pk[ik] = fl*taul;
Pko[ik] = fl*taul;
Pko[ik] = fl2*taul;
}

double r[N], xi[N];
double r[N], xi[N], xio[N];

xi2pk(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pko,r,xio,pclass_sz);
for (ik=0; ik<N; ik++){
// convolution:
xi[ik] = pow(xi[ik],2.);
// cross-convolution of filter1*tau and filter2*tau:
xi[ik] = xi[ik]*xio[ik];
}

pk2xi(N,r,xi,k,Pk,pclass_sz);
Expand Down Expand Up @@ -8088,6 +8104,7 @@ if ((int) pvectsz[pclass_sz->index_part_id_cov_hsv] == 6) {
double lnk[N],lnpk[N];
int ik;
double fl;
double fl2;
double taul;
double l;
double m = exp(logM);
Expand All @@ -8097,21 +8114,22 @@ if ((int) pvectsz[pclass_sz->index_part_id_cov_hsv] == 6) {
lnk[ik] = log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min));
l = k[ik];
fl = get_ksz_filter_at_l(l,pclass_sz);

fl2 = get_ksz_filter_at_l2(l,pclass_sz);
// pvectsz[pclass_sz->index_multipole_for_tau_profile] = l;
double kp = (l+0.5)/chi;
evaluate_galaxy_profile_2h(kp,m_delta_gal,r_delta_gal,c_delta_gal,pvecback,pvectsz,pba,pclass_sz);
taul = pvectsz[pclass_sz->index_galaxy_profile];
Pk[ik] = fl*taul;
Pko[ik] = fl*taul;
Pko[ik] = fl2*taul;
}

double r[N], xi[N];
double r[N], xi[N], xio[N];

xi2pk(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pko,r,xio,pclass_sz);
for (ik=0; ik<N; ik++){
// convolution:
xi[ik] = pow(xi[ik],2.);
// cross-convolution of filter1*tau and filter2*tau:
xi[ik] = xi[ik]*xio[ik];
}

pk2xi(N,r,xi,k,Pk,pclass_sz);
Expand Down Expand Up @@ -8236,6 +8254,7 @@ double k[N], Pk[N],Pko[N];
double lnk[N],lnpk[N];
int ik;
double fl;
double fl2;
double taul;
double l;
// double m = exp(logM);
Expand All @@ -8245,22 +8264,26 @@ k[ik] = exp(log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min)));
lnk[ik] = log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min));
l = k[ik];
fl = get_ksz_filter_at_l(l,pclass_sz);

fl2 = get_ksz_filter_at_l2(l,pclass_sz);
// pvectsz[pclass_sz->index_multipole_for_tau_profile] = l;//pclass_sz->ell_kSZ2_gal_multipole_grid[index_l_2];
evaluate_tau_profile((l+0.5)/chi,pvecback,pvectsz,pba,pclass_sz);
taul = pvectsz[pclass_sz->index_tau_profile];//get_tau_profile_at_z_m_l(z,m,l,pclass_sz,pba);
Pk[ik] = fl*taul;
Pko[ik] = fl*taul;
Pko[ik] = fl2*taul;
// if(l>3e3)
// printf("k = %.5e pk = %.5e\n",l,Pk[ik]);
}

double r[N], xi[N];
double r[N], xi[N], xio[N];

// printf("\n\n####################\n");
// printf("To position space:\n");
// printf("####################\n\n");
// pk2xi(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pko,r,xio,pclass_sz);
for (ik=0; ik<N; ik++){
// cross-convolution of filter1*tau and filter2*tau:
xi[ik] = xi[ik]*xio[ik];
}

pk2xi(N,r,xi,k,Pk,pclass_sz);
xi2pk(N,k,Pk,r,xi,pclass_sz);
for (ik=0; ik<N; ik++){
// printf("r = %.5e xi = %.5e\n",r[ik],xi[ik]);
Expand Down Expand Up @@ -8335,6 +8358,7 @@ double k[N], Pk[N],Pko[N];
double lnk[N],lnpk[N];
int ik;
double fl;
double fl2;
double taul;
double l;
// double m = exp(logM);
Expand All @@ -8344,7 +8368,7 @@ k[ik] = exp(log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min)));
lnk[ik] = log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min));
l = k[ik];
fl = get_ksz_filter_at_l(l,pclass_sz);

fl2 = get_ksz_filter_at_l2(l,pclass_sz);
// pvectsz[pclass_sz->index_multipole_for_tau_profile] = l;//pclass_sz->ell_kSZ2_gal_multipole_grid[index_l_2];
// evaluate_tau_profile((l+0.5)/chi,pvecback,pvectsz,pba,pclass_sz);

Expand All @@ -8354,17 +8378,21 @@ fl = get_ksz_filter_at_l(l,pclass_sz);

taul = pvectsz[pclass_sz->index_galaxy_profile];//get_tau_profile_at_z_m_l(z,m,l,pclass_sz,pba);
Pk[ik] = fl*taul;
Pko[ik] = fl*taul;
Pko[ik] = fl2*taul;
// if(l>3e3)
// printf("k = %.5e pk = %.5e\n",l,Pk[ik]);
}

double r[N], xi[N];
double r[N], xi[N], xio[N];

// printf("\n\n####################\n");
// printf("To position space:\n");
// printf("####################\n\n");
// pk2xi(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pko,r,xio,pclass_sz);
for (ik=0; ik<N; ik++){
// cross-convolution of filter1*tau and filter2*tau:
xi[ik] = xi[ik]*xio[ik];
}

pk2xi(N,r,xi,k,Pk,pclass_sz);
xi2pk(N,k,Pk,r,xi,pclass_sz);
for (ik=0; ik<N; ik++){
// printf("r = %.5e xi = %.5e\n",r[ik],xi[ik]);
Expand Down Expand Up @@ -8439,6 +8467,7 @@ double k[N], Pk[N],Pko[N];
double lnk[N],lnpk[N];
int ik;
double fl;
double fl2;
double taul;
double l;
// double m = exp(logM);
Expand All @@ -8448,22 +8477,26 @@ k[ik] = exp(log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min)));
lnk[ik] = log(l_min)+ik/(N-1.)*(log(l_max)-log(l_min));
l = k[ik];
fl = get_ksz_filter_at_l(l,pclass_sz);

fl2 = get_ksz_filter_at_l2(l,pclass_sz);
// pvectsz[pclass_sz->index_multipole_for_tau_profile] = l;//pclass_sz->ell_kSZ2_gal_multipole_grid[index_l_2];
evaluate_tau_profile((l+0.5)/chi,pvecback,pvectsz,pba,pclass_sz);
taul = pvectsz[pclass_sz->index_tau_profile];//get_tau_profile_at_z_m_l(z,m,l,pclass_sz,pba);
Pk[ik] = fl*taul;
Pko[ik] = fl*taul;
Pko[ik] = fl2*taul;
// if(l>3e3)
// printf("k = %.5e pk = %.5e\n",l,Pk[ik]);
}

double r[N], xi[N];
double r[N], xi[N], xio[N];

// printf("\n\n####################\n");
// printf("To position space:\n");
// printf("####################\n\n");
// pk2xi(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pk,r,xi,pclass_sz);
xi2pk(N,k,Pko,r,xio,pclass_sz);
for (ik=0; ik<N; ik++){
// cross-convolution of filter1*tau and filter2*tau:
xi[ik] = xi[ik]*xio[ik];
}

pk2xi(N,r,xi,k,Pk,pclass_sz);
xi2pk(N,k,Pk,r,xi,pclass_sz);
for (ik=0; ik<N; ik++){
// printf("r = %.5e xi = %.5e\n",r[ik],xi[ik]);
Expand Down Expand Up @@ -11075,6 +11108,18 @@ else
ell);
return fl;
}
double get_ksz_filter_at_l2(double ell,
struct class_sz_structure * pclass_sz){
double fl = 1.;
if (ell <= pclass_sz->l_unwise_filter2[0] || ell >= pclass_sz->l_unwise_filter2[pclass_sz->unwise_filter2_size-1])
fl = 0.;
else
fl = pwl_value_1d(pclass_sz->unwise_filter2_size,
pclass_sz->l_unwise_filter2,
pclass_sz->f_unwise_filter2,
ell);
return fl;
}

double get_M_min_of_z(double ell,
struct class_sz_structure * pclass_sz){
Expand Down
2 changes: 2 additions & 0 deletions class-sz/source/input.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6557,6 +6557,8 @@ class_read_int("no_tt_noise_in_kSZ2X_cov",pclass_sz->no_tt_noise_in_kSZ2X_cov);
// class_read_string("sBBN_file",ppr->sBBN_file);
class_read_string("ksz_filter_file",pclass_sz->ksz_filter_file);
class_read_string("projected_field_filter_file",pclass_sz->ksz_filter_file);
class_read_string("ksz_filter_file2",pclass_sz->ksz_filter_file2);
class_read_string("projected_field_filter_file2",pclass_sz->ksz_filter_file2);
class_read_string("full_path_to_dndz_gal",pclass_sz->full_path_to_dndz_gal);
class_read_string("full_path_and_prefix_to_dndz_ngal",pclass_sz->full_path_and_prefix_to_dndz_ngal);
class_read_string("full_path_to_redshift_dependent_M_min",pclass_sz->full_path_to_redshift_dependent_M_min);
Expand Down
Loading