CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
StackedDensityProfile.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2020 by Giorgio Lesci and Federico Marulli *
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 
35 #include "StackedDensityProfile.h"
36 #include "Data1D_extra.h"
37 #include "CCfits/CCfits"
38 #include<time.h>
39 
40 using namespace cbl;
41 using namespace measure::stackprofile;
42 
43 
44 // ============================================================================
45 
46 void cbl::measure::stackprofile::StackedDensityProfile::m_check_catalogue_variables ()
47 {
48  if (m_galData->isSetVar(cbl::catalogue::Var::_RA_) != true || m_galData->isSetVar(cbl::catalogue::Var::_Dec_) != true || m_galData->isSetVar(cbl::catalogue::Var::_Redshift_) != true || m_galData->isSetVar(cbl::catalogue::Var::_Shear1_) != true || m_galData->isSetVar(cbl::catalogue::Var::_Shear2_) != true || m_galData->isSetVar(cbl::catalogue::Var::_RedshiftMin_) != true || m_galData->isSetVar(cbl::catalogue::Var::_LensingWeight_) != true)
49  cbl::ErrorCBL("You must set the following galaxy catalogue variables: RA, Dec, redshift, redshift_min, shear1, shear2, lensing_weight","m_check_catalogue_variables","StackedDensityProfile.cpp");
50  if (m_cluData->isSetVar(cbl::catalogue::Var::_RA_) !=true || m_cluData->isSetVar(cbl::catalogue::Var::_Dec_) !=true || m_cluData->isSetVar(cbl::catalogue::Var::_Redshift_) !=true || m_cluData->isSetVar(cbl::catalogue::Var::_SN_) !=true || m_cluData->isSetVar(cbl::catalogue::Var::_MassProxy_) !=true)
51  cbl::ErrorCBL("You must set the following cluster catalogue variables: RA, Dec, redshift, S/N, Mass Proxy","m_check_catalogue_variables","StackedDensityProfile.cpp");
52 }
53 
54 
55 // ============================================================================
56 
57 void cbl::measure::stackprofile::StackedDensityProfile::m_resize (const double rad_min, const double rad_max, const int nRad, const bool log_rad)
58 {
59  if (m_z_binEdges.size()-1 != m_proxy_binEdges.size())
60  cbl::ErrorCBL("The number of proxy binnings must be equal to the number of redshift bins!","StackedDensityProfile","StackedDensityProfile.cpp");
61 
62  m_background_idx_z.resize(m_z_binEdges.size()-1, std::vector<int>());
63  m_background_idx_colour.resize(m_z_binEdges.size()-1, std::vector<int>());
64 
65  m_colourSel.resize(m_z_binEdges.size()-1);
66  m_photzSel.resize(m_z_binEdges.size()-1);
67  m_logicSel.resize(m_z_binEdges.size()-1);
68 
69  m_isSet_colourSel.resize(m_z_binEdges.size()-1, false);
70  m_isSet_photzSel.resize(m_z_binEdges.size()-1, false);
71  m_isSet_logicSel.resize(m_z_binEdges.size()-1, false);
72 
73  m_colour_sel_pars.resize(m_z_binEdges.size()-1);
74  m_zphot_sel_pars.resize(m_z_binEdges.size()-1);
75  m_logic_sel_par.resize(m_z_binEdges.size()-1);
76 
77  m_rad_arr.resize(nRad+1);
78  if (log_rad == false){
79  for (int i=0; i < nRad+1; i++)
80  m_rad_arr[i] = rad_min + (rad_max-rad_min)/double(nRad)*i;
81  }else{
82  for (int i = 0; i<nRad+1; i++)
83  m_rad_arr[i] = exp(log(rad_min) + (log(rad_max)-log(rad_min))/double(nRad)*i);
84  }
85 
86  m_ngal_arr.resize(m_rad_arr.size()-1,0);
87  m_deltasigma_t.resize(m_rad_arr.size()-1,0.);
88  m_deltasigma_x.resize(m_rad_arr.size()-1,0.);
89  m_deltasigma_err.resize(m_rad_arr.size()-1,0.);
90  m_wetasquareSum.resize(m_rad_arr.size()-1,0.);
91  m_rad_eff_arr.resize(m_rad_arr.size()-1,0.);
92  m_rad_sigma_arr.resize(m_rad_arr.size()-1,0.);
93  m_wSum.resize(m_rad_arr.size()-1,0.);
94  m_wetaSum.resize(m_rad_arr.size()-1,0.);
95  m_deltasigmaSum.resize(m_rad_arr.size()-1,0.);
96 
97  m_stacked_ngal_arr.resize(m_z_binEdges.size()-1);
98  m_stacked_deltasigma_t.resize(m_z_binEdges.size()-1);
99  m_stacked_deltasigma_x.resize(m_z_binEdges.size()-1);
100  m_stacked_deltasigma_t_err.resize(m_z_binEdges.size()-1);
101  m_stacked_deltasigma_x_err.resize(m_z_binEdges.size()-1);
102  m_stacked_deltasigma_err.resize(m_z_binEdges.size()-1);
103  m_stacked_wetasquareSum.resize(m_z_binEdges.size()-1);
104  m_stacked_wSum.resize(m_z_binEdges.size()-1);
105  m_stacked_wetaSum.resize(m_z_binEdges.size()-1);
106  m_stacked_deltasigmaSum.resize(m_z_binEdges.size()-1);
107  m_stacked_deltasigma_wei.resize(m_z_binEdges.size()-1);
108  m_proxy_eff.resize(m_z_binEdges.size()-1);
109  m_proxy_sigma.resize(m_z_binEdges.size()-1);
110  m_stacked_rad_eff_arr.resize(m_z_binEdges.size()-1);
111  m_stacked_rad_sigma_arr.resize(m_z_binEdges.size()-1);
112  m_single_deltasigma_t.resize(m_z_binEdges.size()-1);
113  m_single_deltasigma_x.resize(m_z_binEdges.size()-1);
114  m_single_deltasigma_err.resize(m_z_binEdges.size()-1);
115  m_single_ngal_arr.resize(m_z_binEdges.size()-1);
116  m_nClu_inBin.resize(m_z_binEdges.size()-1);
117  m_z_eff.resize(m_z_binEdges.size()-1);
118  m_z_sigma.resize(m_z_binEdges.size()-1);
119  for (size_t i=0; i<m_z_binEdges.size()-1; i++){
120  m_stacked_ngal_arr[i].resize(m_proxy_binEdges[i].size()-1);
121  m_stacked_deltasigma_t[i].resize(m_proxy_binEdges[i].size()-1);
122  m_stacked_deltasigma_x[i].resize(m_proxy_binEdges[i].size()-1);
123  m_stacked_deltasigma_t_err[i].resize(m_proxy_binEdges[i].size()-1);
124  m_stacked_deltasigma_x_err[i].resize(m_proxy_binEdges[i].size()-1);
125  m_stacked_deltasigma_err[i].resize(m_proxy_binEdges[i].size()-1);
126  m_stacked_wetasquareSum[i].resize(m_proxy_binEdges[i].size()-1);
127  m_stacked_wSum[i].resize(m_proxy_binEdges[i].size()-1);
128  m_stacked_wetaSum[i].resize(m_proxy_binEdges[i].size()-1);
129  m_stacked_deltasigmaSum[i].resize(m_proxy_binEdges[i].size()-1);
130  m_stacked_deltasigma_wei[i].resize(m_proxy_binEdges[i].size()-1);
131  m_proxy_eff[i].resize(m_proxy_binEdges[i].size()-1);
132  m_proxy_sigma[i].resize(m_proxy_binEdges[i].size()-1);
133  m_stacked_rad_eff_arr[i].resize(m_proxy_binEdges[i].size()-1);
134  m_stacked_rad_sigma_arr[i].resize(m_proxy_binEdges[i].size()-1);
135  m_single_deltasigma_t[i].resize(m_proxy_binEdges[i].size()-1);
136  m_single_deltasigma_x[i].resize(m_proxy_binEdges[i].size()-1);
137  m_single_deltasigma_err[i].resize(m_proxy_binEdges[i].size()-1);
138  m_single_ngal_arr[i].resize(m_proxy_binEdges[i].size()-1);
139  m_nClu_inBin[i].resize(m_proxy_binEdges[i].size()-1,0);
140  m_z_eff[i].resize(m_proxy_binEdges[i].size()-1);
141  m_z_sigma[i].resize(m_proxy_binEdges[i].size()-1);
142  for (size_t j=0; j<m_proxy_binEdges[i].size()-1; j++){
143  m_stacked_ngal_arr[i][j].resize(m_rad_arr.size()-1);
144  m_stacked_deltasigma_t[i][j].resize(m_rad_arr.size()-1);
145  m_stacked_deltasigma_x[i][j].resize(m_rad_arr.size()-1);
146  m_stacked_deltasigma_t_err[i][j].resize(m_rad_arr.size()-1);
147  m_stacked_deltasigma_x_err[i][j].resize(m_rad_arr.size()-1);
148  m_stacked_deltasigma_err[i][j].resize(m_rad_arr.size()-1);
149  m_stacked_wetasquareSum[i][j].resize(m_rad_arr.size()-1);
150  m_stacked_wSum[i][j].resize(m_rad_arr.size()-1);
151  m_stacked_wetaSum[i][j].resize(m_rad_arr.size()-1);
152  m_stacked_deltasigmaSum[i][j].resize(m_rad_arr.size()-1);
153  m_stacked_deltasigma_wei[i][j].resize(m_rad_arr.size()-1);
154  m_stacked_rad_eff_arr[i][j].resize(m_rad_arr.size()-1);
155  m_single_deltasigma_t[i][j].resize(m_rad_arr.size()-1);
156  m_single_deltasigma_x[i][j].resize(m_rad_arr.size()-1);
157  m_single_deltasigma_err[i][j].resize(m_rad_arr.size()-1);
158  m_single_ngal_arr[i][j].resize(m_rad_arr.size()-1);
159  m_stacked_rad_sigma_arr[i][j].resize(m_rad_arr.size()-1);
160  }}
161 
162  m_deltasigma_cov_matr.resize(m_z_binEdges.size()-1);
163  for (size_t i=0; i<m_z_binEdges.size()-1; i++){
164  m_deltasigma_cov_matr[i].resize(m_proxy_binEdges[i].size()-1);
165  for (size_t j=0; j<m_proxy_binEdges[i].size()-1; j++){
166  m_deltasigma_cov_matr[i][j].resize(m_rad_arr.size()-1);
167  for (size_t k=0; k<m_rad_arr.size()-1; k++){
168  m_deltasigma_cov_matr[i][j][k].resize(m_rad_arr.size()-1);
169  }}}
170 }
171 
172 // ============================================================================
173 
174 void cbl::measure::stackprofile::StackedDensityProfile::m_set_multiplicative_calib ()
175 {
176  if (m_m_calib.size() > 0) {
177 
178  if (m_m_calib.size() != m_m_calib_zEdges.size())
179  cbl::ErrorCBL("multiplicative_calibration_stats and multiplicative_calibration_zEdges must have the same size!","m_set_multiplicative_calib","StackedDensityProfile.cpp");
180 
181  for (size_t i=0; i<m_m_calib.size(); i++) {
182  if (m_m_calib[i].size() != 2)
183  cbl::ErrorCBL("All the vectors in multiplicative_calibration_stats must have size = 2","m_set_multiplicative_calib","StackedDensityProfile.cpp");
184  if (m_m_calib_zEdges[i].size() != 2)
185  cbl::ErrorCBL("All the vectors in multiplicative_calibration_zEdges must have size = 2","m_set_multiplicative_calib","StackedDensityProfile.cpp");
186  }
187 
188  coutCBL<<"I'm setting the m parameter for each galaxy by extracting values from Gaussian distributions..."<<std::endl;
189  std::vector<double> m_values((int)(m_galData->nObjects()));
190  for (size_t j=0; j<m_galData->nObjects(); j++) {
191  for (size_t k=0; k<m_m_calib_zEdges.size(); k++) {
192  if (m_galData->redshift(j) >= m_m_calib_zEdges[k][0] && m_galData->redshift(j) <= m_m_calib_zEdges[k][1]) {
193 
194  srand(time(0));
195  cbl::random::NormalRandomNumbers random(m_m_calib[k][0], m_m_calib[k][1], rand());
196  m_values[j] = random();
197 
198  }
199  }
200  }
201 
202  m_galData->set_var(cbl::catalogue::Var::_LensingCalib_, m_values);
203 
204  } else {
205  if (m_galData->isSetVar(cbl::catalogue::Var::_LensingCalib_) != true)
206  cbl::ErrorCBL("You must set the galaxy catalogue variable LensingCalib","m_set_multiplicative_calib","StackedDensityProfile.cpp");
207  }
208 }
209 
210 // ============================================================================
211 
212 bool cbl::measure::stackprofile::StackedDensityProfile::m_colourSelection_gri (const int i_gal, const int clu_zbin)
213 {
214  bool cond1 = m_galData->magnitudeG(i_gal)-m_galData->magnitudeR(i_gal) < m_colour_sel_pars[clu_zbin][0];
215  bool cond2 = m_galData->magnitudeR(i_gal)-m_galData->magnitudeI(i_gal) > m_colour_sel_pars[clu_zbin][1];
216  bool cond3 = m_galData->magnitudeR(i_gal)-m_galData->magnitudeI(i_gal) > m_colour_sel_pars[clu_zbin][2]*(m_galData->magnitudeG(i_gal)-m_galData->magnitudeR(i_gal))+m_colour_sel_pars[clu_zbin][3];
217 
218  bool is_good_colour = (cond1 || cond2 || cond3) && m_galData->redshift(i_gal)>=m_colour_sel_pars[clu_zbin][4];
219 
220  return is_good_colour;
221 }
222 
223 // ============================================================================
224 
225 bool cbl::measure::stackprofile::StackedDensityProfile::m_photzSelection (const int i_clu, const int i_gal, const int clu_zbin)
226 {
227  bool cond1 = m_galData->redshiftMin(i_gal) > m_cluData->redshift(i_clu)+m_zphot_sel_pars[clu_zbin][0];
228  bool cond2 = m_galData->redshift(i_gal)>=m_zphot_sel_pars[clu_zbin][1];
229  bool cond3 = m_galData->redshift(i_gal)<=m_zphot_sel_pars[clu_zbin][2];
230  bool cond4 = m_galData->odds(i_gal) > m_zphot_sel_pars[clu_zbin][3];
231 
232  bool is_good_photz = (cond1 && cond2 && cond3 && cond4);
233 
234  return is_good_photz;
235 }
236 
237 // ============================================================================
238 
239 void cbl::measure::stackprofile::StackedDensityProfile::m_linked_list ()
240 {
241  std::vector<double> ra_limits = {1000,-1000}, dec_limits = {1000,-1000};
242  for (size_t i=0; i<m_galData->nObjects(); i++){
243  if (m_galData->ra(i) < ra_limits[0]) ra_limits[0] = m_galData->ra(i);
244  if (m_galData->ra(i) > ra_limits[1]) ra_limits[1] = m_galData->ra(i);
245  if (m_galData->dec(i) < dec_limits[0]) dec_limits[0] = m_galData->dec(i);
246  if (m_galData->dec(i) > dec_limits[1]) dec_limits[1] = m_galData->dec(i);
247  }
248  m_ra_start = ra_limits[0]; m_dec_start = dec_limits[0];
249  m_nPix.resize(2);
250  m_nPix[0] = (ra_limits[1]-m_ra_start)/m_pix_size +1;
251  m_nPix[1] = (dec_limits[1]-m_dec_start)/m_pix_size +1;
252 
253  m_last.resize(m_nPix[0], std::vector<int>(m_nPix[1]));
254  m_first.resize(m_nPix[0], std::vector<int>(m_nPix[1],-1));
255  m_next.resize(m_galData->nObjects(),-1);
256 
257  int pix[2];
258  for (size_t i=0; i<m_next.size(); i++){
259  // Pixel position
260  pix[0] = std::floor((m_galData->ra(i)-m_ra_start)/m_pix_size);
261  pix[1] = std::floor((m_galData->dec(i)-m_dec_start)/m_pix_size);
262  if (m_first[pix[0]][pix[1]]==-1){ // If m_first==-1 then this is the first particle of the list
263  m_first[pix[0]][pix[1]] = i;
264  m_last[pix[0]][pix[1]] = i;
265  }
266  else{ // otherwise I search for the end of the list and put it there.
267  m_next[m_last[pix[0]][pix[1]]] = i;
268  m_last[pix[0]][pix[1]] = i;
269  }
270  }
271 }
272 
273 // ============================================================================
274 
275 void cbl::measure::stackprofile::StackedDensityProfile::m_add_galaxy (const int i_gal, const int clu_index, const double coscludec, const double clu_dist)
276 {
277  // Consider a spherical triangle with sides a, b, c. We want to find c, the angular distance between two sources.
278  // The law of haversines is: cos(c) = cos(a)cos(b) + sin(a)sin(b)cos(C), where C is the angle of the corner opposite c.
279  // Suppose also that a and b intersect in the north pole of a unit sphere, so that a = 90°-dec1, b = 90°-dec2, and C = ra1-ra2. Thus the formula becomes:
280  // cos(c) = sin(dec1)sin(dec2) + cos(dec1)cos(dec2)cos(ra1-ra2)
281 
282  std::vector<double> ang_dist(2);
283  ang_dist[0] = (m_cluData->ra(clu_index)-m_galData->ra(i_gal))*coscludec;
284  ang_dist[1] = m_galData->dec(i_gal)-m_cluData->dec(clu_index); // If ra1=ra2 -> cos(c) = sin(dec1)sin(dec2) + cos(dec1)cos(dec2) = cos(dec1-dec2)
285 
286  double alpha1 = m_cluData->ra(clu_index);
287  double delta1 = m_cluData->dec(clu_index);
288  double alpha2 = m_galData->ra(i_gal);
289  double delta2 = m_galData->dec(i_gal);
290 
291  double dist = acos(sin(delta1)*sin(delta2)+cos(delta1)*cos(delta2)*cos(alpha1-alpha2))*clu_dist/1.e+6;
292  double dist2 = dist*dist;
293 
294  if (dist2 < m_rad_arr[m_rad_arr.size()-1]*m_rad_arr[m_rad_arr.size()-1] && dist2 >= m_rad_arr[0]*m_rad_arr[0] && m_galData->redshift(i_gal) > m_cluData->redshift(clu_index)+m_delta_redshift){
295 
296  // Find the cluster redshift bin
297  int candidate_cluster_zbin = 0;
298  for (size_t z_clu_idx=0; z_clu_idx<m_z_binEdges.size()-1; z_clu_idx++)
299  if (m_cluData->redshift(clu_index) == m_z_binEdges[0])
300  candidate_cluster_zbin = z_clu_idx;
301  else if (m_z_binEdges[z_clu_idx] < m_cluData->redshift(clu_index) && m_cluData->redshift(clu_index) <= m_z_binEdges[z_clu_idx+1])
302  candidate_cluster_zbin = z_clu_idx;
303 
304  const int cluster_zbin = candidate_cluster_zbin;
305 
306  // Check if the galaxy is a good background candidate
307  std::vector<bool> single_selections = { (this->*m_colourSel[cluster_zbin])(i_gal, cluster_zbin), (this->*m_photzSel[cluster_zbin])(clu_index, i_gal, cluster_zbin) };
308  bool selection = m_logicSel[cluster_zbin](single_selections);
309 
310  if (selection){
311  // Fill the vector of indices of the selected lensed galaxies (used to write files)
312  if (single_selections[0])
313  m_background_idx_colour[cluster_zbin].emplace_back(i_gal);
314  if (single_selections[1])
315  m_background_idx_z[cluster_zbin].emplace_back(i_gal);
316 
317  // locate correct annulus
318  std::vector<double>::iterator it = upper_bound(m_rad_arr.begin(), m_rad_arr.end(), dist);
319  int p = it-m_rad_arr.begin()-1;
320  m_ngal_arr[p]++;
321 
322  double ds = 2.99792e+09*m_cosm->D_A(0.,m_galData->redshift(i_gal))/cbl::par::cc*100; // in pc/h
323  double dds = 2.99792e+09*m_cosm->D_A(m_cluData->redshift(clu_index),m_galData->redshift(i_gal))/cbl::par::cc*100; // in pc/h
324  double eta = clu_dist*dds/ds; // (sigma_crit / (c*c/4/pi/G))^-1
325  std::vector<double> ell = {m_galData->shear1(i_gal),m_galData->shear2(i_gal)};
326 
327  double phi = {atan2(ang_dist[1],ang_dist[0])}; // azimuthal angle between cluster and lensed source
328  double cos2phi = cos(2*phi); double sin2phi = sin(2*phi);
329  std::vector<double> new_ell = {-ell[0]*cos2phi-ell[1]*sin2phi, ell[0]*sin2phi-ell[1]*cos2phi}; // tang. shear = - (shear1*cos(2phi) + shear2*sin(2phi))
330 
331  m_wSum[p] += m_galData->lensingWeight(i_gal);
332  m_wetaSum[p] += m_galData->lensingWeight(i_gal)*eta;
333  m_wetasquareSum[p] += m_galData->lensingWeight(i_gal)*eta*eta;
334  m_deltasigmaSum[p] += m_galData->lensingWeight(i_gal)*eta*eta*m_galData->lensingCalib(i_gal); // uncertainty on m_wetasquareSum[p]
335  m_deltasigma_t[p] += new_ell[0]*m_galData->lensingWeight(i_gal)*eta;
336  m_deltasigma_x[p] += new_ell[1]*m_galData->lensingWeight(i_gal)*eta;
337  m_rad_eff_arr[p] += m_galData->lensingWeight(i_gal)*eta*eta*std::pow(dist,-m_rad_alpha);
338  m_rad_sigma_arr[p] += m_galData->lensingWeight(i_gal)*eta*eta*std::pow(dist,2.);
339  }
340  }
341 }
342 
343 
344 // ============================================================================
345 
346 void cbl::measure::stackprofile::StackedDensityProfile::m_profile (const int clu_index)
347 {
348  for (size_t i=0; i<m_rad_arr.size()-1; i++){
349  m_ngal_arr[i]=0; m_deltasigma_t[i]=0; m_deltasigma_x[i]=0; m_deltasigma_err[i]=0; m_wetasquareSum[i]=0; m_rad_eff_arr[i]=0; m_rad_sigma_arr[i]=0; m_wSum[i]=0; m_wetaSum[i]=0; m_deltasigmaSum[i]=0;
350  }
351 
352  const double clu_dist = m_cosm->D_A(0.,m_cluData->redshift(clu_index)) * 1.e6; // in pc/h
353  const double coscludec = cos(m_cluData->dec(clu_index));
354  std::vector<double> clu_pix = {std::floor((m_cluData->ra(clu_index)-m_ra_start)/m_pix_size), std::floor((m_cluData->dec(clu_index)-m_dec_start)/m_pix_size)};
355  double max_rad = m_rad_arr[m_rad_arr.size()-1]/clu_dist*1.e+6;
356 
357  // limits of relevant pixels
358  std::vector<int> pix_ra_lim(2); std::vector<int> pix_dec_lim(2);
359  pix_ra_lim[0] = std::min(std::max(0,int(clu_pix[0] - max_rad/coscludec/m_pix_size)),int(m_nPix[0]-1));
360  pix_ra_lim[1] = std::max(std::min(int(m_nPix[0]-1),int(clu_pix[0] + max_rad/coscludec/m_pix_size + 1)),0);
361  pix_dec_lim[0] = std::min(std::max(0,int(clu_pix[1] - max_rad/m_pix_size)),int(m_nPix[1]-1));
362  pix_dec_lim[1] = std::max(std::min(int(m_nPix[1]-1),int(clu_pix[1] + max_rad/m_pix_size + 1)),0);
363 
364  std::vector<int> pix(2);
365  for (pix[0] = pix_ra_lim[0]; pix[0] <= pix_ra_lim[1]; pix[0]++){
366  for (pix[1] = pix_dec_lim[0]; pix[1] <= pix_dec_lim[1]; pix[1]++){
367  int igal = m_first[pix[0]][pix[1]];
368  while ( igal!=-1 ){
369  m_add_galaxy(igal,clu_index,coscludec,clu_dist);
370  igal = m_next[igal];
371  }
372  }
373  }
374 
375  // Normalize
376  for (size_t i=0; i<m_rad_arr.size()-1; i++){
377  double deltasigmaMean = m_deltasigmaSum[i]/m_wetasquareSum[i];
378  m_deltasigma_t[i] /= (1.+deltasigmaMean)*(m_wetasquareSum[i]/m_sigma_fac); // final units: h*M_sun/(pc^2)
379  m_deltasigma_x[i] /= (1.+deltasigmaMean)*(m_wetasquareSum[i]/m_sigma_fac);
380  m_deltasigma_err[i] = 1./sqrt(m_wetasquareSum[i])*m_sigma_fac; // statistical uncertainty from weighted mean
381  m_rad_eff_arr[i] = std::pow(m_rad_eff_arr[i]/m_wetasquareSum[i],-1./m_rad_alpha);
382  m_rad_sigma_arr[i] = std::pow(m_rad_sigma_arr[i]/m_wetasquareSum[i] - std::pow(m_rad_eff_arr[i],2.), 0.5);
383  }
384 }
385 
386 // ============================================================================
387 
388 void cbl::measure::stackprofile::StackedDensityProfile::m_stacker ()
389 {
390 
391  m_set_multiplicative_calib();
392 
393  std::cout<<std::endl; coutCBL<<"Performing the stacking..."<<std::endl;
394 
395  // This is the inverse of Fac_eta_ToInv_Sigma_Cr in Mauro Sereno's code
396  m_sigma_fac = 1.6628e+12; // c^2/(4piG) in Msun/pc
397 
398  // linked_list pix_size could depend on cluster min z and rad_max
399  m_linked_list();
400 
401  // --------------------------------------------------------------------------------------------
402  // ---------------------------------- Perform the stacking ------------------------------------
403  // --------------------------------------------------------------------------------------------
404  std::vector<std::vector<double>> wetasquareSum_j(m_z_binEdges.size()-1,std::vector<double>(m_proxy_binEdges[0].size()-1)); // Array used for the stacked proxy and redshift
405  std::vector<std::vector<double>> wetasquareSum_tot(m_z_binEdges.size()-1,std::vector<double>(m_proxy_binEdges[0].size()-1)); // Array used for the stacked proxy and redshift
406 
407  double nClu = (double)(m_cluData->nObjects());
408  for (int i=0; i<nClu; i++){
409  if (m_cluData->redshift(i) >= m_z_binEdges[0] && m_cluData->redshift(i) < m_z_binEdges[m_z_binEdges.size()-1] && m_cluData->sn(i) > m_SN_min){
410  std::vector<double>::iterator it_p = std::upper_bound(m_z_binEdges.begin(), m_z_binEdges.end(), m_cluData->redshift(i)); // find redshift and mass proxy bin where to put cluster
411  int p = it_p-m_z_binEdges.begin()-1;
412 
413  if (m_cluData->mass_proxy(i) >= m_proxy_binEdges[p][0] && m_cluData->mass_proxy(i) < m_proxy_binEdges[p][m_proxy_binEdges[p].size()-1]){
414  std::vector<double>::iterator it_q = std::upper_bound(m_proxy_binEdges[p].begin(), m_proxy_binEdges[p].end(), m_cluData->mass_proxy(i));
415  int q = it_q-m_proxy_binEdges[p].begin()-1;
416 
417  if (p>=0 && p < (int)(m_z_binEdges.size()) && q>=0 && q < (int)(m_proxy_binEdges[p].size())) {
418 
419  bool cluster_has_background = false;
420  m_profile(i);
421 
422  for (size_t j=0; j<m_rad_arr.size()-1; j++){
423  if (m_ngal_arr[j] > 0){
424 
425  cluster_has_background = true;
426 
427  double err = m_deltasigma_err[j]; double wei = 1./err/err;
428  m_stacked_deltasigma_t[p][q][j] += m_deltasigma_t[j]*wei; m_stacked_deltasigma_x[p][q][j] += m_deltasigma_x[j]*wei;
429  m_stacked_deltasigma_wei[p][q][j] += wei; m_stacked_ngal_arr[p][q][j] += m_ngal_arr[j];
430  m_single_deltasigma_t[p][q][j].emplace_back(m_deltasigma_t[j]); m_single_deltasigma_x[p][q][j].emplace_back(m_deltasigma_x[j]); // For bootstrap/resampling in general
431  m_single_deltasigma_err[p][q][j].emplace_back(m_deltasigma_err[j]); m_single_ngal_arr[p][q][j].emplace_back(m_ngal_arr[j]); // For bootstrap/resampling in general
432  wetasquareSum_j[p][q] += m_wetasquareSum[j];
433  m_stacked_wetasquareSum[p][q][j] += m_wetasquareSum[j];
434  m_stacked_rad_eff_arr[p][q][j] += pow(m_rad_eff_arr[j],m_obs_gamma)*m_wetasquareSum[j];
435  if (m_ngal_arr[j] > 1) m_stacked_rad_sigma_arr[p][q][j] += pow(m_rad_sigma_arr[j],m_obs_gamma*2.)*m_wetasquareSum[j];
436  }
437  }
438 
439  if (cluster_has_background) {
440  m_nClu_inBin[p][q] ++;
441 
442  m_proxy_eff[p][q] += pow(m_cluData->mass_proxy(i),m_obs_gamma)*wetasquareSum_j[p][q]; // wetasquaresum = \Sum (weight * (sigma_crit / (c*c/4/pi/G))^-2)
443  m_proxy_sigma[p][q] += pow(m_cluData->mass_proxy(i),2*m_obs_gamma)*wetasquareSum_j[p][q];
444  m_z_eff[p][q] += pow(m_cluData->redshift(i),m_obs_gamma)*wetasquareSum_j[p][q];
445  m_z_sigma[p][q] += pow(m_cluData->redshift(i),2*m_obs_gamma)*wetasquareSum_j[p][q];
446  wetasquareSum_tot[p][q] += wetasquareSum_j[p][q];
447  wetasquareSum_j[p][q] = 0;
448  }
449 
450  }
451  }
452  }
453  const int progress = (int)((i+1)/nClu*100);
454  coutCBL << progress << "% \r"; std::cout.flush();
455  }
456  for (size_t p=0; p<m_z_binEdges.size()-1; p++){
457  for (size_t q=0; q<m_proxy_binEdges[p].size()-1; q++){
458  m_proxy_eff[p][q] = pow(m_proxy_eff[p][q]/wetasquareSum_tot[p][q],1./m_obs_gamma);
459  m_proxy_sigma[p][q] = pow((m_proxy_sigma[p][q]/wetasquareSum_tot[p][q] - pow(m_proxy_eff[p][q],2*m_obs_gamma))/m_nClu_inBin[p][q],1./2./m_obs_gamma);
460  m_z_eff[p][q] = pow(m_z_eff[p][q]/wetasquareSum_tot[p][q],1./m_obs_gamma);
461  m_z_sigma[p][q] = pow((m_z_sigma[p][q]/wetasquareSum_tot[p][q] - pow(m_z_eff[p][q],2*m_obs_gamma))/m_nClu_inBin[p][q],1./2./m_obs_gamma);
462  for (size_t j=0; j<m_rad_arr.size()-1; j++){
463  m_stacked_deltasigma_t[p][q][j] /= m_stacked_deltasigma_wei[p][q][j];
464  m_stacked_deltasigma_x[p][q][j] /= m_stacked_deltasigma_wei[p][q][j];
465  m_stacked_deltasigma_t_err[p][q][j] = sqrt(1./m_stacked_deltasigma_wei[p][q][j]);
466  m_stacked_deltasigma_x_err[p][q][j] = sqrt(1./m_stacked_deltasigma_wei[p][q][j]);
467  // Make m_rad_arr
468  // for radius we consider total scatter of galaxy data in bin (from weighted average of scatters), and not sample-average error
469  m_stacked_rad_eff_arr[p][q][j] = pow(m_stacked_rad_eff_arr[p][q][j]/m_stacked_wetasquareSum[p][q][j],1./m_obs_gamma);;
470  m_stacked_rad_sigma_arr[p][q][j] = pow(m_stacked_rad_sigma_arr[p][q][j]/m_stacked_wetasquareSum[p][q][j],1./2./m_obs_gamma);
471  }}}
472 }
473 
474 
475 // ============================================================================
476 
477 std::shared_ptr<data::Data> cbl::measure::stackprofile::StackedDensityProfile::m_make_bootstrap (const std::vector<int> z_proxy_bin)
478 {
479  m_stacker(); // Perform the stacking
480 
481  std::cout<<std::endl; coutCBL<<"Computing the bootstrap covariance matrix..."<<std::endl<<std::endl;
482  std::default_random_engine generator;
483  for (size_t p=0; p<m_z_binEdges.size()-1; p++){
484  for (size_t q=0; q<m_proxy_binEdges[p].size()-1; q++){
485  std::uniform_int_distribution<int> distribution(0,m_single_ngal_arr[p][q][0].size()-1);
486  std::vector< std::vector<double> > deltasigma_t_boot(m_n_resampling,std::vector<double>(m_rad_arr.size()-1,0.));
487  std::vector< std::vector<double> > deltasigma_x_boot(m_n_resampling,std::vector<double>(m_rad_arr.size()-1,0.));
488  std::vector< std::vector<double> > deltasigma_wei_boot(m_n_resampling,std::vector<double>(m_rad_arr.size()-1,0.));
489  std::vector<double> deltasigma_t_boot_mean(m_rad_arr.size()-1,0.);
490 
491  for (int i=0; i<m_n_resampling; i++){
492  for (size_t j=0; j<m_single_ngal_arr[p][q][0].size(); j++){
493  int clu_num = distribution(generator); // random extraction of a cluster
494  for (size_t k=0; k<m_rad_arr.size()-1; k++){
495  if (m_single_ngal_arr[p][q][k][clu_num] > 0){
496  deltasigma_t_boot[i][k] += m_single_deltasigma_t[p][q][k][clu_num]/m_single_deltasigma_err[p][q][k][clu_num]/m_single_deltasigma_err[p][q][k][clu_num];
497  deltasigma_x_boot[i][k] += m_single_deltasigma_x[p][q][k][clu_num]/m_single_deltasigma_err[p][q][k][clu_num]/m_single_deltasigma_err[p][q][k][clu_num];
498  deltasigma_wei_boot[i][k] += 1./m_single_deltasigma_err[p][q][k][clu_num]/m_single_deltasigma_err[p][q][k][clu_num];
499  }}}
500  for (size_t k=0; k<m_rad_arr.size()-1; k++){
501  deltasigma_t_boot[i][k] /= deltasigma_wei_boot[i][k];
502  deltasigma_x_boot[i][k] /= deltasigma_wei_boot[i][k];
503  }
504  }
505  for (size_t k=0; k<m_rad_arr.size()-1; k++){
506  for (int i=0; i<m_n_resampling; i++){
507  deltasigma_t_boot_mean[k] += deltasigma_t_boot[i][k];
508  }
509  deltasigma_t_boot_mean[k] /= double(m_n_resampling);
510  }
511  for (size_t k=0; k<m_rad_arr.size()-1; k++){
512  for (size_t l=0; l<m_rad_arr.size()-1; l++){
513  for (int i=0; i<m_n_resampling; i++){
514  m_deltasigma_cov_matr[p][q][k][l] += (deltasigma_t_boot[i][k]-deltasigma_t_boot_mean[k])*(deltasigma_t_boot[i][l]-deltasigma_t_boot_mean[l]);
515  }
516  m_deltasigma_cov_matr[p][q][k][l] /= double(m_n_resampling);
517  }}
518  }}
519 
520  std::vector<double> bins = m_stacked_rad_eff_arr[z_proxy_bin[0]][z_proxy_bin[1]];
521  std::vector<double> profile = m_stacked_deltasigma_t[z_proxy_bin[0]][z_proxy_bin[1]];
522  std::vector<double> error; // Set the diagonal of the covariance matrix as the error
523  std::vector<std::vector<double>> extra_info (4);
524  for (size_t i=0; i<m_rad_arr.size()-1; i++){
525  error.emplace_back(sqrt(m_deltasigma_cov_matr[z_proxy_bin[0]][z_proxy_bin[1]][i][i]));
526  extra_info[0].emplace_back(m_z_eff[z_proxy_bin[0]][z_proxy_bin[1]]);
527  extra_info[1].emplace_back(m_z_sigma[z_proxy_bin[0]][z_proxy_bin[1]]);
528  extra_info[2].emplace_back(m_proxy_eff[z_proxy_bin[0]][z_proxy_bin[1]]);
529  extra_info[3].emplace_back(m_proxy_sigma[z_proxy_bin[0]][z_proxy_bin[1]]);
530  }
531  auto dataset = std::make_shared<data::Data1D_extra> (data::Data1D_extra(bins, profile, error, extra_info));
532 
533  return dataset;
534 }
535 
536 // ============================================================================
537 
538 bool cbl::measure::stackprofile::StackedDensityProfile::m_check_file (const std::string checked_file, const std::vector<int> z_proxy_bin)
539 {
540  bool out = false;
541 
542  m_inputs_to_str.clear();
543 
544  m_inputs_to_str.emplace_back(cbl::conv(m_cosm->OmegaM(),cbl::par::fDP5));
545  m_inputs_to_str.emplace_back(cbl::conv(m_cosm->Omega_baryon(),cbl::par::fDP5));
546  m_inputs_to_str.emplace_back(cbl::conv(m_cosm->Omega_k(),cbl::par::fDP5));
547  m_inputs_to_str.emplace_back(cbl::conv(m_cosm->Omega_neutrinos(),cbl::par::fDP5));
548  m_inputs_to_str.emplace_back(cbl::conv(pow(10,9)*m_cosm->scalar_amp(),cbl::par::fDP3));
549  m_inputs_to_str.emplace_back(cbl::conv(m_cosm->tau(),cbl::par::fDP5));
550  m_inputs_to_str.emplace_back(cbl::conv(m_cosm->n_spec(),cbl::par::fDP5));
551  m_inputs_to_str.emplace_back(cbl::conv(m_cosm->hh(),cbl::par::fDP5));
552 
553  m_inputs_to_str.emplace_back(cbl::conv(m_galData->nObjects(),cbl::par::fINT));
554  m_inputs_to_str.emplace_back(cbl::conv(m_cluData->nObjects(),cbl::par::fINT));
555 
556  m_inputs_to_str.emplace_back(cbl::conv(m_delta_redshift,cbl::par::fDP5));
557 
558  for (size_t i=0; i<m_m_calib.size(); i++) {
559  m_inputs_to_str.emplace_back(cbl::conv(m_m_calib[i][0],cbl::par::fDP5));
560  m_inputs_to_str.emplace_back(cbl::conv(m_m_calib[i][1],cbl::par::fDP5));
561  m_inputs_to_str.emplace_back(cbl::conv(m_m_calib_zEdges[i][0],cbl::par::fDP5));
562  m_inputs_to_str.emplace_back(cbl::conv(m_m_calib_zEdges[i][1],cbl::par::fDP5));
563  }
564 
565  for (size_t j=0; j<m_zphot_sel_pars.size(); j++)
566  for (size_t jj=0; jj<m_zphot_sel_pars[j].size(); jj++)
567  m_inputs_to_str.emplace_back(cbl::conv(m_zphot_sel_pars[j][jj],cbl::par::fDP5));
568 
569  for (size_t j=0; j<m_colour_sel_pars.size(); j++)
570  for (size_t jj=0; jj<m_colour_sel_pars[j].size(); jj++)
571  m_inputs_to_str.emplace_back(cbl::conv(m_colour_sel_pars[j][jj],cbl::par::fDP5));
572 
573  for (size_t j=0; j<m_logic_sel_par.size(); j++)
574  for (size_t jj=0; jj<m_logic_sel_par[j].size(); jj++)
575  m_inputs_to_str.emplace_back(m_logic_sel_par[j][jj]);
576 
577  for (size_t j=0; j<m_z_binEdges.size(); j++){
578  m_inputs_to_str.emplace_back(cbl::conv(m_z_binEdges[j],cbl::par::fDP5));
579  }
580  for (size_t i=0; i<m_proxy_binEdges.size(); i++){
581  for (size_t j=0; j<m_proxy_binEdges[i].size(); j++){
582  m_inputs_to_str.emplace_back(cbl::conv(m_proxy_binEdges[i][j],cbl::par::fDP5));
583  }
584  }
585  for (size_t j=0; j<m_rad_arr.size(); j++){
586  m_inputs_to_str.emplace_back(cbl::conv(m_rad_arr[j],cbl::par::fDP5));
587  }
588  m_inputs_to_str.emplace_back(cbl::conv(m_SN_min,cbl::par::fDP5));
589  m_inputs_to_str.emplace_back(cbl::conv(m_pix_size,cbl::par::fDP5));
590  m_inputs_to_str.emplace_back(cbl::conv(m_n_resampling,cbl::par::fINT));
591  m_inputs_to_str.emplace_back(cbl::conv(m_rad_alpha,cbl::par::fDP5));
592  m_inputs_to_str.emplace_back(cbl::conv(m_obs_gamma,cbl::par::fDP5));
593 
594  std::string check_string = "# ";
595  for (size_t i=0; i<m_inputs_to_str.size(); i++){
596  check_string += m_inputs_to_str[i];
597  check_string += " ";
598  }
599  std::ifstream file_input(checked_file.c_str(),std::ios::in);
600  std::vector<double> rad, DSigma, error;
601  std::vector<std::vector<double>> extra_info (4);
602  if (file_input.std::ifstream::is_open()){
603  std::string line;
604  std::getline(file_input,line);
605  if (line == check_string){
606  out = true;
607  while (std::getline(file_input,line)){
608  if (line.at(0) != '#'){
609  std::stringstream ss(line);
610  std::vector<double> num; double NUM;
611  while (ss>>NUM) num.emplace_back(NUM);
612  if ((int)(num[0])==z_proxy_bin[0] && (int)(num[1])==z_proxy_bin[1]){
613  rad.emplace_back(num[6]);
614  DSigma.emplace_back(num[8]);
615  error.emplace_back(num[14]);
616  extra_info[0].emplace_back(num[2]);
617  extra_info[1].emplace_back(num[3]);
618  extra_info[2].emplace_back(num[4]);
619  extra_info[3].emplace_back(num[5]);
620  }
621  }
622  }
623  if (rad.size()==0)
624  cbl::ErrorCBL("The chosen redshift and proxy bins are not valid!","StackedDensityProfile","StackedDensityProfile.cpp");
625  }
626  }
627  file_input.clear(); file_input.close();
628 
629  m_dataset = std::make_shared<data::Data1D_extra> (data::Data1D_extra(rad, DSigma, error, extra_info));
630 
631  return out;
632 }
633 
634 // ============================================================================
635 
636 void cbl::measure::stackprofile::StackedDensityProfile::m_write(const std::string output_dir, const std::string output_file)
637 {
638  const std::string mkdir = "mkdir -p "+output_dir; if (system(mkdir.c_str())) {}
639 
640  std::ofstream out_stack (output_dir+"/"+output_file+".dat",std::ios::out);
641 
642  // Write the inputs in the header
643  out_stack << "# ";
644  for (size_t i=0; i<m_inputs_to_str.size(); i++)
645  out_stack << m_inputs_to_str[i] << " ";
646  out_stack << std::endl;
647 
648  // Write the rest of the file
649  out_stack << "# Redshift array:" << std::endl << "# ";
650  for (size_t i=0; i<m_z_binEdges.size(); i++)
651  out_stack << m_z_binEdges[i] << " ";
652  out_stack << std::endl;
653  out_stack << "# Proxy arrays:";
654  for (size_t i=0; i<m_proxy_binEdges.size(); i++){
655  out_stack << std::endl << "# ";
656  for (size_t j=0; j<m_proxy_binEdges[i].size(); j++){
657  out_stack << m_proxy_binEdges[i][j] << " ";
658  if (i==m_proxy_binEdges.size()-1 && j==m_proxy_binEdges[i].size()-1)
659  out_stack << std::endl;
660  }
661  }
662  out_stack << "# 1) redshift bin" << std::endl;
663  out_stack << "# 2) proxy bin" << std::endl;
664  out_stack << "# 3) mean redshift" << std::endl;
665  out_stack << "# 4) redshift sigma" << std::endl;
666  out_stack << "# 5) mean proxy" << std::endl;
667  out_stack << "# 6) proxy sigma" << std::endl;
668  out_stack << "# 7) mean radius" << std::endl;
669  out_stack << "# 8) radius sigma" << std::endl;
670  out_stack << "# 9) deltasigma_+ [h*MSun/pc^2]" << std::endl;
671  out_stack << "# 10) deltasigma_+ error [h*MSun/pc^2]" << std::endl;
672  out_stack << "# 11) deltasigma_x [h*MSun/pc^2]" << std::endl;
673  out_stack << "# 12) deltasigma_x error [h*MSun/pc^2]" << std::endl;
674  out_stack << "# 13) number of clusters in the bins of z and proxy" << std::endl;
675  out_stack << "# 14) total weight" << std::endl;
676  out_stack << "# 15) bootstrap error on deltasigma_+" << std::endl;
677  out_stack << "# 16) rad min" << std::endl;
678  out_stack << "# 17) rad max" << std::endl;
679  out_stack << "# 18) number of background galaxies in the bins of z, proxy and radius" << std::endl;
680 
681  for (size_t i=0; i<m_z_binEdges.size()-1; i++){
682  for (size_t j=0; j<m_proxy_binEdges[i].size()-1; j++){
683  for (size_t k=0; k<m_rad_arr.size()-1; k++){
684  out_stack << i << " " << j << " " << m_z_eff[i][j] << " " << m_z_sigma[i][j] << " " << m_proxy_eff[i][j] << " " << m_proxy_sigma[i][j] << " " << m_stacked_rad_eff_arr[i][j][k] << " " << m_stacked_rad_sigma_arr[i][j][k] << " " << m_stacked_deltasigma_t[i][j][k] << " " << m_stacked_deltasigma_t_err[i][j][k] << " " << m_stacked_deltasigma_x[i][j][k] << " " << m_stacked_deltasigma_x_err[i][j][k] << " " << m_nClu_inBin[i][j] << " " << m_stacked_wetasquareSum[i][j][k] << " " << sqrt(m_deltasigma_cov_matr[i][j][k][k]) << " " << m_rad_arr[k] << " " << m_rad_arr[k+1] << " " << m_stacked_ngal_arr[i][j][k] << std::endl;
685  }
686  }
687  }
688  out_stack.clear(); out_stack.close(); coutCBL << "I wrote the file: " << output_dir+"/"+output_file+".dat" << std::endl << std::endl;
689 
690 
691  // Write the covariance matrix files
692 
693  const std::string cov_dir = output_dir+"/covariance/";
694  const std::string mkdir_cov = "mkdir -p "+cov_dir; if (system(mkdir_cov.c_str())) {}
695 
696  for (size_t i=0; i<m_z_binEdges.size()-1; i++) {
697  for (size_t j=0; j<m_proxy_binEdges[i].size()-1; j++) {
698 
699  const std::string file_path = cov_dir+"covariance_"+cbl::conv(i,cbl::par::fINT)+cbl::conv(j,cbl::par::fINT)+".dat";
700  std::ofstream cov_file (file_path,std::ios::out);
701 
702  for (size_t r=0; r<m_rad_arr.size()-1; r++)
703  for (size_t r2=0; r2<m_rad_arr.size()-1; r2++)
704  if (r2 < m_rad_arr.size()-2)
705  cov_file << m_deltasigma_cov_matr[i][j][r][r2] << " ";
706  else
707  cov_file << m_deltasigma_cov_matr[i][j][r][r2] << std::endl;
708 
709  cov_file.clear(); cov_file.close(); coutCBL << "I wrote the file: " << file_path << std::endl << std::endl;
710  }
711  }
712 
713 
714  // Write the files containing the indices of the background galaxies
715 
716  if (m_write_background) {
717 
718  for (size_t i=0; i<m_z_binEdges.size()-1; i++) {
719 
720  const std::string path = output_dir+"/background_galaxies/zbin_"+cbl::conv(i,cbl::par::fINT)+"/";
721  const std::string mkdir_bkg = "mkdir -p "+path; if (system(mkdir_bkg.c_str())) {}
722 
723  // Since multiple clusters may share the same sources, build vectors of unique elements
724  const std::vector<int> idx_z = cbl::different_elements(m_background_idx_z[i]);
725  const std::vector<int> idx_c = cbl::different_elements(m_background_idx_colour[i]);
726 
727  const std::string z_selection_file_name = path+"redshift_selection.dat";
728  std::ofstream z_selection_file (z_selection_file_name,std::ios::out);
729  for (size_t j=0; j<idx_z.size(); j++) {
730  z_selection_file << idx_z[j] << std::endl;
731  }
732  z_selection_file.clear(); z_selection_file.close(); coutCBL << "I wrote the file: " << z_selection_file_name << std::endl << std::endl;
733 
734  const std::string colour_selection_file_name = path+"colour_selection.dat";
735  std::ofstream colour_selection_file (colour_selection_file_name,std::ios::out);
736  for (size_t j=0; j<idx_c.size(); j++) {
737  colour_selection_file << idx_c[j] << std::endl;
738  }
739  colour_selection_file.clear(); colour_selection_file.close(); coutCBL << "I wrote the file: " << colour_selection_file_name << std::endl << std::endl;
740 
741  }
742 
743  }
744 
745 }
746 
747 // ============================================================================
748 
749 void cbl::measure::stackprofile::StackedDensityProfile::set_colour_selection (const std::string colour_space, const int z_bin, const double C1, const double C2, const double C3, const double C4, const double z_min)
750 {
751  if (z_bin > (int)(m_z_binEdges.size()-1))
752  cbl::ErrorCBL("The chosen cluster z bin index is outside the z bin index range!","set_colour_selection_gri","StackedDensityProfile.cpp");
753 
754  if (colour_space == "ri_gr")
755  m_colourSel[z_bin] = &cbl::measure::stackprofile::StackedDensityProfile::m_colourSelection_gri;
756  else
757  cbl::ErrorCBL("The selected colour selection is not available!","set_colour_selection","StackedDensityProfile.cpp");
758 
759  m_colour_sel_pars[z_bin] = {C1, C2, C3, C4, z_min};
760  m_isSet_colourSel[z_bin] = true;
761 
762  coutCBL << "I've set the colour selection in the redshift bin " << z_bin << std::endl;
763 }
764 
765 // ============================================================================
766 
767 void cbl::measure::stackprofile::StackedDensityProfile::set_colour_selection (const std::string colour_space, const double C1, const double C2, const double C3, const double C4, const double z_min)
768 {
769  for (size_t i=0; i<m_z_binEdges.size()-1; i++)
770  this->set_colour_selection(colour_space, i, C1, C2, C3, C4, z_min);
771 }
772 
773 // ============================================================================
774 
775 void cbl::measure::stackprofile::StackedDensityProfile::set_zphot_selection (const int z_bin, const double deltaz, const double zgal_min, const double zgal_max, const double ODDS_min)
776 {
777  if (z_bin > (int)(m_z_binEdges.size()-1))
778  cbl::ErrorCBL("The chosen cluster z bin index is outside the z bin index range!","set_zphot_selection","StackedDensityProfile.cpp");
779 
780  m_zphot_sel_pars[z_bin] = {deltaz, zgal_min, zgal_max, ODDS_min};
781  m_photzSel[z_bin] = &cbl::measure::stackprofile::StackedDensityProfile::m_photzSelection;
782  m_isSet_photzSel[z_bin] = true;
783 
784  coutCBL << "I've set the redshift selection in the redshift bin " << z_bin << std::endl;
785 }
786 
787 // ============================================================================
788 
789 void cbl::measure::stackprofile::StackedDensityProfile::set_zphot_selection (const double deltaz, const double zgal_min, const double zgal_max, const double ODDS_min)
790 {
791  for (size_t i=0; i<m_z_binEdges.size()-1; i++)
792  this->set_zphot_selection(i, deltaz, zgal_min, zgal_max, ODDS_min);
793 }
794 
795 // ============================================================================
796 
797 void cbl::measure::stackprofile::StackedDensityProfile::set_logic_selection (const int z_bin, const std::string sel)
798 {
799  if (z_bin > (int)(m_z_binEdges.size()-1))
800  cbl::ErrorCBL("The chosen cluster z bin index is outside the z bin index range!","set_logic_selection","StackedDensityProfile.cpp");
801 
802  std::function<bool(const std::vector<bool>)> logicSelection_and = [] (const std::vector<bool> x){
803  // x[0]->colour selection, x[1]->photz selection
804  bool thebool = (x[0] && x[1]);
805  return thebool;
806  };
807 
808  std::function<bool(const std::vector<bool>)> logicSelection_or = [] (const std::vector<bool> x){
809  // x[0]->colour selection, x[1]->photz selection
810  bool thebool = (x[0] || x[1]);
811  return thebool;
812  };
813  if (sel=="and")
814  m_logicSel[z_bin] = logicSelection_and;
815  else if (sel=="or")
816  m_logicSel[z_bin] = logicSelection_or;
817  else
818  cbl::ErrorCBL("The chosen logic operator between the selection criteria is not valid! Choose \'and\' or \'or\'.","set_logic_selection","StackedDensityProfile.cpp");
819 
820  m_isSet_logicSel[z_bin] = true;
821  m_logic_sel_par[z_bin] = {sel};
822 
823  coutCBL << "I've set the logic selection in the redshift bin " << z_bin << std::endl;
824 }
825 
826 // ============================================================================
827 
828 void cbl::measure::stackprofile::StackedDensityProfile::set_logic_selection (const std::string sel)
829 {
830  for (size_t i=0; i<m_z_binEdges.size()-1; i++)
831  this->set_logic_selection(i, sel);
832 }
833 
834 // ============================================================================
835 
836 void cbl::measure::stackprofile::StackedDensityProfile::measure (const std::vector<int> z_proxy_bin, const std::string output_dir, const std::string output_file_root, const ErrorType errorType, const int n_resampling, const bool write_background)
837 {
838  m_write_background = write_background;
839 
840  for (size_t i=0; i<m_z_binEdges.size()-1; i++) {
841  if (m_isSet_colourSel[i] == false)
842  cbl::ErrorCBL("The colour selection for the redshift bin "+cbl::conv(i,cbl::par::fINT)+" is not set!","measure","StackedDensityProfile.cpp");
843  if (m_isSet_photzSel[i] == false)
844  cbl::ErrorCBL("The redshift selection for the redshift bin "+cbl::conv(i,cbl::par::fINT)+" is not set!","measure","StackedDensityProfile.cpp");
845  if (m_isSet_logicSel[i] == false)
846  cbl::ErrorCBL("The logic selection for the redshift bin "+cbl::conv(i,cbl::par::fINT)+" is not set!","measure","StackedDensityProfile.cpp");
847  }
848 
849  m_n_resampling = n_resampling;
850 
851  bool file_exists = m_check_file(output_dir+"/"+output_file_root+".dat", z_proxy_bin);
852 
853  if (file_exists){
854  m_measure_is_read = true;
855  coutCBL << "I've just read the already existing stacking file!" << std::endl;
856  }else{
857  switch (errorType){
858  case (ErrorType::_Bootstrap_):
859  m_dataset = m_make_bootstrap(z_proxy_bin);
860  break;
861  default:
862  ErrorCBL("The input ErrorType is not allowed!", "measure", "StackedDensityProfile.cpp");
863  }
864  m_write(output_dir, output_file_root);
865  }
866 }
867 
868 // ============================================================================
869 
870 void cbl::measure::stackprofile::StackedDensityProfile::write (const std::string dir, const std::string file)
871 {
872  std::string mkdir = "mkdir -p "+dir;
873  if (system(mkdir.c_str())) {}
874 
875  std::string header = "bin_center profile error mean_redshift mean_redshift_error mean_mass_proxy mean_mass_proxy_error";
876  m_dataset->write(dir, file, header, 4, 8);
877 }
878 
879 
880 // ============================================================================
881 
882 cbl::measure::stackprofile::StackedDensityProfile::StackedDensityProfile (cosmology::Cosmology cosm, const catalogue::Catalogue gal_cat, const catalogue::Catalogue clu_cat, const double delta_redshift, std::vector<double> z_binEdges, std::vector<std::vector<double>> proxy_binEdges, const double rad_min, const double rad_max, const int nRad, const bool log_rad, const double SN_min, const double pix_size, const std::vector<std::vector<double>> multiplicative_calibration_stats, const std::vector<std::vector<double>> multiplicative_calibration_zEdges, const double rad_alpha, const double obs_gamma)
883 {
884  // WARNING: if you change the arguments of this constructor, you must edit m_check_file !!!
885 
886  m_cosm = std::make_shared<cosmology::Cosmology>(cosmology::Cosmology(std::move(cosm)));
887  m_cosm->set_unit(true); // Force cosmological units
888 
889  m_galData = std::make_shared<catalogue::Catalogue>(catalogue::Catalogue(std::move(gal_cat)));
890  m_cluData = std::make_shared<catalogue::Catalogue>(catalogue::Catalogue(std::move(clu_cat)));
891 
892  m_delta_redshift = delta_redshift;
893  m_z_binEdges = z_binEdges;
894  m_proxy_binEdges = proxy_binEdges;
895  m_SN_min = SN_min;
896  m_pix_size = pix_size * (cbl::par::pi/180.);
897  m_m_calib = multiplicative_calibration_stats;
898  m_m_calib_zEdges = multiplicative_calibration_zEdges;
899  m_rad_alpha = rad_alpha;
900  m_obs_gamma = obs_gamma;
901 
902  m_check_catalogue_variables();
903  m_resize(rad_min,rad_max,nRad,log_rad);
904 
905  m_measure_is_read = false;
906 }
Data1D_extra.h
The class Data1D_extra.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
StackedDensityProfile.h
The class StackedDensityProfile.
cbl::catalogue::Catalogue
The class Catalogue.
Definition: Catalogue.h:654
cbl::cosmology::Cosmology
The class Cosmology.
Definition: Cosmology.h:277
cbl::data::Data1D_extra
The class Data1D_extra.
Definition: Data1D_extra.h:52
cbl::measure::stackprofile::StackedDensityProfile::m_colourSelection_gri
bool m_colourSelection_gri(const int i_gal, const int clu_zbin)
colour selection function in the colour-colour space ( ) vs ( )
Definition: StackedDensityProfile.cpp:212
cbl::measure::stackprofile::StackedDensityProfile::set_logic_selection
void set_logic_selection(const int z_bin, const std::string sel)
Set the logic operator linking redshift and colour selections in the i-th cluster redshift bin.
Definition: StackedDensityProfile.cpp:797
cbl::measure::stackprofile::StackedDensityProfile::StackedDensityProfile
StackedDensityProfile()=default
default constructor
cbl::measure::stackprofile::StackedDensityProfile::m_add_galaxy
void m_add_galaxy(const int i_gal, const int clu_index, const double coscludec, const double clu_dist)
add galaxies to profiles
Definition: StackedDensityProfile.cpp:275
cbl::measure::stackprofile::StackedDensityProfile::m_profile
void m_profile(const int clu_index)
compute a single cluster profile
Definition: StackedDensityProfile.cpp:346
cbl::measure::stackprofile::StackedDensityProfile::m_resize
void m_resize(const double rad_min, const double rad_max, const int nRad, const bool log_rad)
resize the private member arrays
Definition: StackedDensityProfile.cpp:57
cbl::measure::stackprofile::StackedDensityProfile::set_colour_selection
void set_colour_selection(const std::string colour_space, const int z_bin, const double C1, const double C2, const double C3, const double C4, const double z_min)
Set the colour selection in the i-th cluster redshift bin.
Definition: StackedDensityProfile.cpp:749
cbl::measure::stackprofile::StackedDensityProfile::m_make_bootstrap
std::shared_ptr< data::Data > m_make_bootstrap(const std::vector< int > z_proxy_bin)
bootstrap resampling to retrieve the covariance matrix of the stacked signal, providing the stacked p...
Definition: StackedDensityProfile.cpp:477
cbl::measure::stackprofile::StackedDensityProfile::m_photzSelection
bool m_photzSelection(const int i_clu, const int i_gal, const int clu_zbin)
redshift selection function
Definition: StackedDensityProfile.cpp:225
cbl::measure::stackprofile::StackedDensityProfile::m_check_file
bool m_check_file(const std::string checked_file, const std::vector< int > z_proxy_bin)
check if the results obtained from the stacking have already been written to a file
Definition: StackedDensityProfile.cpp:538
cbl::measure::stackprofile::StackedDensityProfile::measure
void measure(const std::vector< int > z_proxy_bin, const std::string output_dir, const std::string output_file_root, const ErrorType errorType=ErrorType::_Bootstrap_, const int n_resampling=10000, const bool write_background=false)
measure the stacked profiles in all the bins of redshift and mass proxy, providing in output the stac...
Definition: StackedDensityProfile.cpp:836
cbl::measure::stackprofile::StackedDensityProfile::m_write
void m_write(const std::string output_dir, const std::string output_file)
write the results obtained from the stacking on file
Definition: StackedDensityProfile.cpp:636
cbl::measure::stackprofile::StackedDensityProfile::write
void write(const std::string dir, const std::string file)
write on file the measure in a given bin of redshift and mass proxy
Definition: StackedDensityProfile.cpp:870
cbl::measure::stackprofile::StackedDensityProfile::m_check_catalogue_variables
void m_check_catalogue_variables()
check if all the necessary variables in the galaxy and cluster catalogues are set
Definition: StackedDensityProfile.cpp:46
cbl::measure::stackprofile::StackedDensityProfile::set_zphot_selection
void set_zphot_selection(const int z_bin, const double deltaz, const double zgal_min, const double zgal_max, const double ODDS_min)
Set the redshift selection in the i-th cluster redshift bin.
Definition: StackedDensityProfile.cpp:775
cbl::measure::stackprofile::StackedDensityProfile::m_set_multiplicative_calib
void m_set_multiplicative_calib()
set the multiplicative shear calibration, m
Definition: StackedDensityProfile.cpp:174
cbl::random::NormalRandomNumbers
The class NormalRandomNumbers.
Definition: RandomNumbers.h:393
cbl::par::fDP3
static const char fDP3[]
conversion symbol for: double -> std::string
Definition: Constants.h:136
cbl::par::fDP5
static const char fDP5[]
conversion symbol for: double -> std::string
Definition: Constants.h:142
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::cc
static const double cc
: the speed of light in vacuum (the value is exact) [km/sec]
Definition: Constants.h:221
cbl::catalogue::Var::_LensingCalib_
@ _LensingCalib_
lensing calibration factor
cbl::catalogue::Var::_MassProxy_
@ _MassProxy_
mass proxy
cbl::catalogue::Var::_Shear2_
@ _Shear2_
second component of the shear signal
cbl::catalogue::Var::_SN_
@ _SN_
Signal-to-noise ratio.
cbl::catalogue::Var::_RA_
@ _RA_
Right Ascension.
cbl::catalogue::Var::_RedshiftMin_
@ _RedshiftMin_
minimum redshift
cbl::catalogue::Var::_Shear1_
@ _Shear1_
first component of the shear signal
cbl::catalogue::Var::_LensingWeight_
@ _LensingWeight_
lensing weight
cbl::catalogue::Var::_Dec_
@ _Dec_
Declination.
cbl::measure::ErrorType
ErrorType
the two-point correlation function error type
Definition: Measure.h:59
cbl::measure::ErrorType::_Bootstrap_
@ _Bootstrap_
Bootstrap resampling.
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
cbl::distribution
void distribution(std::vector< double > &xx, std::vector< double > &fx, std::vector< double > &err, const std::vector< double > FF, const std::vector< double > WW, const int nbin, const bool linear=true, const std::string file_out=par::defaultString, const double fact=1., const double V1=par::defaultDouble, const double V2=par::defaultDouble, const std::string bin_type="Linear", const bool conv=false, const double sigma=0.)
derive and store the number distribution of a given std::vector
Definition: Func.cpp:1654
cbl::conv
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
cbl::different_elements
std::vector< T > different_elements(const std::vector< T > vect_input)
get the unique elements of a std::vector
Definition: Kernel.h:1349
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