CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Pair2D_extra.cpp
Go to the documentation of this file.
1 /*******************************************************************
2  * Copyright (C) 2015 by 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 "Pair2D_extra.h"
36 
37 using namespace std;
38 
39 using namespace cbl;
40 using namespace catalogue;
41 using namespace pairs;
42 
43 
44 // ============================================================================================
45 
46 
47 void cbl::pairs::Pair2D_comovingCartesian_linlin_extra::put (const std::shared_ptr<Object> obj1, const std::shared_ptr<Object> obj2)
48 {
49  const double rp = perpendicular_distance(obj1->ra(), obj2->ra(), obj1->dec(), obj2->dec(), obj1->dc(), obj2->dc());
50  const double pi = fabs(obj1->dc()-obj2->dc());
51 
52  if (m_rpMin<rp && rp<m_rpMax && m_piMin<pi && pi<m_piMax) {
53 
54  const int ir = max(0, min(int((rp-m_rpMin)*m_binSize_inv_D1), m_nbins_D1));
55  const int jr = max(0, min(int((pi-m_piMin)*m_binSize_inv_D2), m_nbins_D2));
56 
57  const double angWeight = (m_angularWeight==nullptr) ? 1.
58  : max(0., m_angularWeight(converted_angle(angular_distance(obj1->xx()/obj1->dc(), obj2->xx()/obj2->dc(), obj1->yy()/obj1->dc(), obj2->yy()/obj2->dc(), obj1->zz()/obj1->dc(), obj2->zz()/obj2->dc()), CoordinateUnits::_radians_, m_angularUnits)));
59 
60  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
61 
62  m_PP2D[ir][jr] ++;
63  m_PP2D_weighted[ir][jr] += WeightTOT;
64 
65  if (m_PP2D_weighted[ir][jr]>0) {
66 
67  const double scale_D1_mean_p = m_scale_D1_mean[ir][jr];
68  const double scale_D2_mean_p = m_scale_D2_mean[ir][jr];
69  m_scale_D1_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(rp-scale_D1_mean_p);
70  m_scale_D2_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pi-scale_D2_mean_p);
71  m_scale_D1_S[ir][jr] += WeightTOT*(rp-scale_D1_mean_p)*(rp-m_scale_D1_mean[ir][jr]);
72  m_scale_D2_S[ir][jr] += WeightTOT*(pi-scale_D2_mean_p)*(pi-m_scale_D2_mean[ir][jr]);
73 
74  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
75  const double z_mean_p = m_z_mean[ir][jr];
76  m_z_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pair_redshift-z_mean_p);
77  m_z_S[ir][jr] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[ir][jr]);
78 
79  }
80 
81  }
82 }
83 
84 // ============================================================================================
85 
86 
87 void cbl::pairs::Pair2D_comovingCartesian_linlog_extra::put (const std::shared_ptr<Object> obj1, const std::shared_ptr<Object> obj2)
88 {
89  const double rp = perpendicular_distance(obj1->ra(), obj2->ra(), obj1->dec(), obj2->dec(), obj1->dc(), obj2->dc());
90  const double pi = fabs(obj1->dc()-obj2->dc());
91 
92  if (m_rpMin<rp && rp<m_rpMax && m_piMin<pi && pi<m_piMax) {
93 
94  const int ir = max(0, min(int((rp-m_rpMin)*m_binSize_inv_D1), m_nbins_D1));
95  const int jr = max(0, min(int((log10(pi)-log10(m_piMin))*m_binSize_inv_D2), m_nbins_D2));
96 
97  const double angWeight = (m_angularWeight==nullptr) ? 1.
98  : max(0., m_angularWeight(converted_angle(angular_distance(obj1->xx()/obj1->dc(), obj2->xx()/obj2->dc(), obj1->yy()/obj1->dc(), obj2->yy()/obj2->dc(), obj1->zz()/obj1->dc(), obj2->zz()/obj2->dc()), CoordinateUnits::_radians_, m_angularUnits)));
99 
100  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
101 
102  m_PP2D[ir][jr] ++;
103  m_PP2D_weighted[ir][jr] += WeightTOT;
104 
105  if (m_PP2D_weighted[ir][jr]>0) {
106 
107  const double scale_D1_mean_p = m_scale_D1_mean[ir][jr];
108  const double scale_D2_mean_p = m_scale_D2_mean[ir][jr];
109  m_scale_D1_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(rp-scale_D1_mean_p);
110  m_scale_D2_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pi-scale_D2_mean_p);
111  m_scale_D1_S[ir][jr] += WeightTOT*(rp-scale_D1_mean_p)*(rp-m_scale_D1_mean[ir][jr]);
112  m_scale_D2_S[ir][jr] += WeightTOT*(pi-scale_D2_mean_p)*(pi-m_scale_D2_mean[ir][jr]);
113 
114  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
115  const double z_mean_p = m_z_mean[ir][jr];
116  m_z_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pair_redshift-z_mean_p);
117  m_z_S[ir][jr] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[ir][jr]);
118 
119  }
120 
121  }
122 }
123 
124 
125 // ============================================================================================
126 
127 
128 void cbl::pairs::Pair2D_comovingCartesian_loglin_extra::put (const std::shared_ptr<Object> obj1, const std::shared_ptr<Object> obj2)
129 {
130  const double rp = perpendicular_distance(obj1->ra(), obj2->ra(), obj1->dec(), obj2->dec(), obj1->dc(), obj2->dc());
131  const double pi = fabs(obj1->dc()-obj2->dc());
132 
133  if (m_rpMin<rp && rp<m_rpMax && m_piMin<pi && pi<m_piMax) {
134 
135  const int ir = max(0, min(int((log10(rp)-log10(m_rpMin))*m_binSize_inv_D1), m_nbins_D1));
136  const int jr = max(0, min(int((pi-m_piMin)*m_binSize_inv_D2), m_nbins_D2));
137 
138  const double angWeight = (m_angularWeight==nullptr) ? 1.
139  : max(0., m_angularWeight(converted_angle(angular_distance(obj1->xx()/obj1->dc(), obj2->xx()/obj2->dc(), obj1->yy()/obj1->dc(), obj2->yy()/obj2->dc(), obj1->zz()/obj1->dc(), obj2->zz()/obj2->dc()), CoordinateUnits::_radians_, m_angularUnits)));
140 
141  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
142 
143  m_PP2D[ir][jr] ++;
144  m_PP2D_weighted[ir][jr] += WeightTOT;
145 
146  if (m_PP2D_weighted[ir][jr]>0) {
147 
148  const double scale_D1_mean_p = m_scale_D1_mean[ir][jr];
149  const double scale_D2_mean_p = m_scale_D2_mean[ir][jr];
150  m_scale_D1_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(rp-scale_D1_mean_p);
151  m_scale_D2_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pi-scale_D2_mean_p);
152  m_scale_D1_S[ir][jr] += WeightTOT*(rp-scale_D1_mean_p)*(rp-m_scale_D1_mean[ir][jr]);
153  m_scale_D2_S[ir][jr] += WeightTOT*(pi-scale_D2_mean_p)*(pi-m_scale_D2_mean[ir][jr]);
154 
155  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
156  const double z_mean_p = m_z_mean[ir][jr];
157  m_z_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pair_redshift-z_mean_p);
158  m_z_S[ir][jr] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[ir][jr]);
159 
160  }
161 
162  }
163 }
164 
165 
166 // ============================================================================================
167 
168 
169 void cbl::pairs::Pair2D_comovingCartesian_loglog_extra::put (const std::shared_ptr<Object> obj1, const std::shared_ptr<Object> obj2)
170 {
171  const double rp = perpendicular_distance(obj1->ra(), obj2->ra(), obj1->dec(), obj2->dec(), obj1->dc(), obj2->dc());
172  const double pi = fabs(obj1->dc()-obj2->dc());
173 
174  if (m_rpMin<rp && rp<m_rpMax && m_piMin<pi && pi<m_piMax) {
175 
176  const int ir = max(0, min(int((log10(rp)-log10(m_rpMin))*m_binSize_inv_D1), m_nbins_D1));
177  const int jr = max(0, min(int((log10(pi)-log10(m_piMin))*m_binSize_inv_D2), m_nbins_D2));
178 
179  const double angWeight = (m_angularWeight==nullptr) ? 1.
180  : max(0., m_angularWeight(converted_angle(angular_distance(obj1->xx()/obj1->dc(), obj2->xx()/obj2->dc(), obj1->yy()/obj1->dc(), obj2->yy()/obj2->dc(), obj1->zz()/obj1->dc(), obj2->zz()/obj2->dc()), CoordinateUnits::_radians_, m_angularUnits)));
181 
182  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
183 
184  m_PP2D[ir][jr] ++;
185  m_PP2D_weighted[ir][jr] += WeightTOT;
186 
187  if (m_PP2D_weighted[ir][jr]>0) {
188 
189  const double scale_D1_mean_p = m_scale_D1_mean[ir][jr];
190  const double scale_D2_mean_p = m_scale_D2_mean[ir][jr];
191  m_scale_D1_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(rp-scale_D1_mean_p);
192  m_scale_D2_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pi-scale_D2_mean_p);
193  m_scale_D1_S[ir][jr] += WeightTOT*(rp-scale_D1_mean_p)*(rp-m_scale_D1_mean[ir][jr]);
194  m_scale_D2_S[ir][jr] += WeightTOT*(pi-scale_D2_mean_p)*(pi-m_scale_D2_mean[ir][jr]);
195 
196  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
197  const double z_mean_p = m_z_mean[ir][jr];
198  m_z_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pair_redshift-z_mean_p);
199  m_z_S[ir][jr] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[ir][jr]);
200 
201  }
202 
203  }
204 }
205 
206 
207 // ============================================================================================
208 
209 
210 void cbl::pairs::Pair2D_comovingPolar_linlin_extra::put (const std::shared_ptr<Object> obj1, const std::shared_ptr<Object> obj2)
211 {
212  double rr = Euclidean_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz());
213  double mu = fabs(obj1->dc()-obj2->dc())/rr;
214 
215  if (m_rMin<rr && rr<m_rMax && m_muMin<mu && mu<m_muMax) {
216 
217  const int ir = max(0, min(int((rr-m_rMin)*m_binSize_inv_D1), m_nbins_D1));
218  const int jr = max(0, min(int((mu-m_muMin)*m_binSize_inv_D2), m_nbins_D2));
219 
220  const double angWeight = (m_angularWeight==nullptr) ? 1.
221  : max(0., m_angularWeight(converted_angle(angular_distance(obj1->xx()/obj1->dc(), obj2->xx()/obj2->dc(), obj1->yy()/obj1->dc(), obj2->yy()/obj2->dc(), obj1->zz()/obj1->dc(), obj2->zz()/obj2->dc()), CoordinateUnits::_radians_, m_angularUnits)));
222 
223  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
224 
225  m_PP2D[ir][jr] ++;
226  m_PP2D_weighted[ir][jr] += WeightTOT;
227 
228  if (m_PP2D_weighted[ir][jr]>0) {
229 
230  const double scale_D1_mean_p = m_scale_D1_mean[ir][jr];
231  const double scale_D2_mean_p = m_scale_D2_mean[ir][jr];
232  m_scale_D1_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(rr-scale_D1_mean_p);
233  m_scale_D2_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(mu-scale_D2_mean_p);
234  m_scale_D1_S[ir][jr] += WeightTOT*(rr-scale_D1_mean_p)*(rr-m_scale_D1_mean[ir][jr]);
235  m_scale_D2_S[ir][jr] += WeightTOT*(mu-scale_D2_mean_p)*(mu-m_scale_D2_mean[ir][jr]);
236 
237  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
238  const double z_mean_p = m_z_mean[ir][jr];
239  m_z_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pair_redshift-z_mean_p);
240  m_z_S[ir][jr] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[ir][jr]);
241 
242  }
243 
244  }
245 }
246 
247 
248 // ============================================================================================
249 
250 
251 void cbl::pairs::Pair2D_comovingPolar_linlog_extra::put (const std::shared_ptr<Object> obj1, const std::shared_ptr<Object> obj2)
252 {
253  double rr = Euclidean_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz());
254  double mu = fabs(obj1->dc()-obj2->dc())/rr;
255 
256  if (m_rMin<rr && rr<m_rMax && m_muMin<mu && mu<m_muMax) {
257 
258  const int ir = max(0, min(int((rr-m_rMin)*m_binSize_inv_D1), m_nbins_D1));
259  const int jr = max(0, min(int((log10(mu)-log10(m_muMin))*m_binSize_inv_D2), m_nbins_D2));
260 
261  const double angWeight = (m_angularWeight==nullptr) ? 1.
262  : max(0., m_angularWeight(converted_angle(angular_distance(obj1->xx()/obj1->dc(), obj2->xx()/obj2->dc(), obj1->yy()/obj1->dc(), obj2->yy()/obj2->dc(), obj1->zz()/obj1->dc(), obj2->zz()/obj2->dc()), CoordinateUnits::_radians_, m_angularUnits)));
263 
264  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
265 
266  m_PP2D[ir][jr] ++;
267  m_PP2D_weighted[ir][jr] += WeightTOT;
268 
269  if (m_PP2D_weighted[ir][jr]>0) {
270 
271  const double scale_D1_mean_p = m_scale_D1_mean[ir][jr];
272  const double scale_D2_mean_p = m_scale_D2_mean[ir][jr];
273  m_scale_D1_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(rr-scale_D1_mean_p);
274  m_scale_D2_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(mu-scale_D2_mean_p);
275  m_scale_D1_S[ir][jr] += WeightTOT*(rr-scale_D1_mean_p)*(rr-m_scale_D1_mean[ir][jr]);
276  m_scale_D2_S[ir][jr] += WeightTOT*(mu-scale_D2_mean_p)*(mu-m_scale_D2_mean[ir][jr]);
277 
278  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
279  const double z_mean_p = m_z_mean[ir][jr];
280  m_z_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pair_redshift-z_mean_p);
281  m_z_S[ir][jr] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[ir][jr]);
282 
283  }
284 
285  }
286 }
287 
288 
289 // ============================================================================================
290 
291 
292 void cbl::pairs::Pair2D_comovingPolar_loglin_extra::put (const std::shared_ptr<Object> obj1, const std::shared_ptr<Object> obj2)
293 {
294  double rr = Euclidean_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz());
295  double mu = fabs(obj1->dc()-obj2->dc())/rr;
296 
297  if (m_rMin<rr && rr<m_rMax && m_muMin<mu && mu<m_muMax) {
298 
299  const int ir = max(0, min(int((log10(rr)-log10(m_rMin))*m_binSize_inv_D1), m_nbins_D1));
300  const int jr = max(0, min(int((mu-m_muMin)*m_binSize_inv_D2), m_nbins_D2));
301 
302  const double angWeight = (m_angularWeight==nullptr) ? 1.
303  : max(0., m_angularWeight(converted_angle(angular_distance(obj1->xx()/obj1->dc(), obj2->xx()/obj2->dc(), obj1->yy()/obj1->dc(), obj2->yy()/obj2->dc(), obj1->zz()/obj1->dc(), obj2->zz()/obj2->dc()), CoordinateUnits::_radians_, m_angularUnits)));
304 
305  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
306 
307  m_PP2D[ir][jr] ++;
308  m_PP2D_weighted[ir][jr] += WeightTOT;
309 
310  if (m_PP2D_weighted[ir][jr]>0) {
311 
312  const double scale_D1_mean_p = m_scale_D1_mean[ir][jr];
313  const double scale_D2_mean_p = m_scale_D2_mean[ir][jr];
314  m_scale_D1_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(rr-scale_D1_mean_p);
315  m_scale_D2_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(mu-scale_D2_mean_p);
316  m_scale_D1_S[ir][jr] += WeightTOT*(rr-scale_D1_mean_p)*(rr-m_scale_D1_mean[ir][jr]);
317  m_scale_D2_S[ir][jr] += WeightTOT*(mu-scale_D2_mean_p)*(mu-m_scale_D2_mean[ir][jr]);
318 
319  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
320  const double z_mean_p = m_z_mean[ir][jr];
321  m_z_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pair_redshift-z_mean_p);
322  m_z_S[ir][jr] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[ir][jr]);
323 
324  }
325 
326  }
327 }
328 
329 
330 // ============================================================================================
331 
332 
333 void cbl::pairs::Pair2D_comovingPolar_loglog_extra::put (const std::shared_ptr<Object> obj1, const std::shared_ptr<Object> obj2)
334 {
335  double rr = Euclidean_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz());
336  double mu = fabs(obj1->dc()-obj2->dc())/rr;
337 
338  if (m_rMin<rr && rr<m_rMax && m_muMin<mu && mu<m_muMax) {
339 
340  const int ir = max(0, min(int((log10(rr)-log10(m_rMin))*m_binSize_inv_D1), m_nbins_D1));
341  const int jr = max(0, min(int((log10(mu)-log10(m_muMin))*m_binSize_inv_D2), m_nbins_D2));
342 
343  const double angWeight = (m_angularWeight==nullptr) ? 1.
344  : max(0., m_angularWeight(converted_angle(angular_distance(obj1->xx()/obj1->dc(), obj2->xx()/obj2->dc(), obj1->yy()/obj1->dc(), obj2->yy()/obj2->dc(), obj1->zz()/obj1->dc(), obj2->zz()/obj2->dc()), CoordinateUnits::_radians_, m_angularUnits)));
345 
346  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
347 
348  m_PP2D[ir][jr] ++;
349  m_PP2D_weighted[ir][jr] += WeightTOT;
350 
351  if (m_PP2D_weighted[ir][jr]>0) {
352 
353  const double scale_D1_mean_p = m_scale_D1_mean[ir][jr];
354  const double scale_D2_mean_p = m_scale_D2_mean[ir][jr];
355  m_scale_D1_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(rr-scale_D1_mean_p);
356  m_scale_D2_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(mu-scale_D2_mean_p);
357  m_scale_D1_S[ir][jr] += WeightTOT*(rr-scale_D1_mean_p)*(rr-m_scale_D1_mean[ir][jr]);
358  m_scale_D2_S[ir][jr] += WeightTOT*(mu-scale_D2_mean_p)*(mu-m_scale_D2_mean[ir][jr]);
359 
360  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
361  const double z_mean_p = m_z_mean[ir][jr];
362  m_z_mean[ir][jr] += WeightTOT/m_PP2D_weighted[ir][jr]*(pair_redshift-z_mean_p);
363  m_z_S[ir][jr] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[ir][jr]);
364 
365  }
366 
367  }
368 }
369 
370 
371 // ============================================================================================
372 
373 
374 void cbl::pairs::Pair2D_extra::add_data2D (const int i, const int j, const std::vector<double> data)
375 {
376  /*
377  checkDim(m_PP2D, i, j, "m_PP2D", false);
378  checkDim(m_PP2D_weighted, i, j, "m_PP2D_weighted", false);
379  checkDim(m_scale_D1_mean, i, j, "m_scale_D1_mean", false);
380  checkDim(m_scale_D2_mean, i, j, "m_scale_D2_mean", false);
381  checkDim(m_scale_D1_sigma, i, j, "m_scale_D1_sigma", false);
382  checkDim(m_scale_D2_sigma, i, j, "m_scale_D2_sigma", false);
383  checkDim(m_z_mean, i, j, "m_z_mean", false);
384  checkDim(m_z_sigma, i, j, "m_z_sigma", false);
385  checkDim(data, 8, "data");
386  */
387 
388  // mean scale in the first dimension up to this computation step
389  const double scale_D1_mean_temp_p = m_scale_D1_mean[i][j];
390 
391  // mean scale in the second dimention up to this computation step
392  const double scale_D2_mean_temp_p = m_scale_D2_mean[i][j];
393 
394  // mean redshift up to this computation step
395  const double z_mean_temp_p = m_z_mean[i][j];
396 
397  // number of pairs
398  m_PP2D[i][j] += data[0];
399 
400  // number of weighted pairs
401  m_PP2D_weighted[i][j] += data[1];
402 
403 
404  if (m_PP2D_weighted[i][j]>0) {
405 
406  // mean separation in the first dimension
407  m_scale_D1_mean[i][j] += data[1]/m_PP2D_weighted[i][j]*(data[2]-scale_D1_mean_temp_p);
408 
409  // mean separation in the first dimension
410  m_scale_D2_mean[i][j] += data[1]/m_PP2D_weighted[i][j]*(data[4]-scale_D2_mean_temp_p);
411 
412  // mean redshift
413  m_z_mean[i][j] += data[1]/m_PP2D_weighted[i][j]*(data[6]-z_mean_temp_p);
414 
415  // compute the weighted standard deviation of the scale distribution in the first dimension
416  m_scale_D1_S[i][j] += data[3]+pow(data[2]-scale_D1_mean_temp_p, 2)*data[1]*(m_PP2D_weighted[i][j]-data[1])/m_PP2D_weighted[i][j];
417  m_scale_D1_sigma[i][j] = sqrt(m_scale_D1_S[i][j]/m_PP2D_weighted[i][j]);
418 
419  // compute the weighted standard deviation of the scale distribution in the second dimension
420  m_scale_D2_S[i][j] += data[5]+pow(data[4]-scale_D2_mean_temp_p, 2)*data[1]*(m_PP2D_weighted[i][j]-data[1])/m_PP2D_weighted[i][j];
421  m_scale_D2_sigma[i][j] = sqrt(m_scale_D2_S[i][j]/m_PP2D_weighted[i][j]);
422 
423  // compute the weighted standard deviation of the redshift distribution
424  m_z_S[i][j] += data[7]+pow(data[6]-z_mean_temp_p, 2)*data[1]*(m_PP2D_weighted[i][j]-data[1])/m_PP2D_weighted[i][j];
425  m_z_sigma[i][j] = sqrt(m_z_S[i][j]/m_PP2D_weighted[i][j]);
426 
427  }
428 
429 }
430 
431 
432 // ============================================================================================
433 
434 
435 void cbl::pairs::Pair2D_extra::add_data2D (const int i, const int j, const std::shared_ptr<pairs::Pair> pair, const double ww)
436 {
437  add_data2D(i, j, {ww*pair->PP2D(i, j), ww*pair->PP2D_weighted(i, j), pair->scale_D1_mean(i, j), pair->scale_D1_S(i, j), pair->scale_D2_mean(i, j), pair->scale_D2_S(i, j), pair->z_mean(i, j), pair->z_S(i, j)});
438 }
439 
440 
441 // ============================================================================================
442 
443 
444 void cbl::pairs::Pair2D_extra::Sum (const std::shared_ptr<Pair> pair, const double ww)
445 {
446  if (m_nbins_D1 != pair->nbins_D1() || m_nbins_D2 != pair->nbins_D2())
447  ErrorCBL("dimension problems!", "Sum", "Pair2D_extra.cpp");
448 
449  for (int i=0; i<m_nbins_D1; ++i)
450  for (int j=0; j<m_nbins_D2; ++j)
451  add_data2D(i, j, pair, ww);
452 }
453 
Pair2D_extra.h
The classes Pair2D_extra*.
cbl::pairs::Pair2D_comovingCartesian_linlin_extra::put
void put(const std::shared_ptr< catalogue::Object > obj1, const std::shared_ptr< catalogue::Object > obj2) override
estimate the distance between two objects and update the pair vector accordingly
Definition: Pair2D_extra.cpp:47
cbl::pairs::Pair2D_comovingCartesian_linlog_extra::put
void put(const std::shared_ptr< catalogue::Object > obj1, const std::shared_ptr< catalogue::Object > obj2) override
estimate the distance between two objects and update the pair vector accordingly
Definition: Pair2D_extra.cpp:87
cbl::pairs::Pair2D_comovingCartesian_loglin_extra::put
void put(const std::shared_ptr< catalogue::Object > obj1, const std::shared_ptr< catalogue::Object > obj2) override
estimate the distance between two objects and update the pair vector accordingly
Definition: Pair2D_extra.cpp:128
cbl::pairs::Pair2D_comovingCartesian_loglog_extra::put
void put(const std::shared_ptr< catalogue::Object > obj1, const std::shared_ptr< catalogue::Object > obj2) override
estimate the distance between two objects and update the pair vector accordingly
Definition: Pair2D_extra.cpp:169
cbl::pairs::Pair2D_comovingPolar_linlin_extra::put
void put(const std::shared_ptr< catalogue::Object > obj1, const std::shared_ptr< catalogue::Object > obj2) override
estimate the distance between two objects and update the pair vector accordingly
Definition: Pair2D_extra.cpp:210
cbl::pairs::Pair2D_comovingPolar_linlog_extra::put
void put(const std::shared_ptr< catalogue::Object > obj1, const std::shared_ptr< catalogue::Object > obj2) override
estimate the distance between two objects and update the pair vector accordingly
Definition: Pair2D_extra.cpp:251
cbl::pairs::Pair2D_comovingPolar_loglin_extra::put
void put(const std::shared_ptr< catalogue::Object > obj1, const std::shared_ptr< catalogue::Object > obj2) override
estimate the distance between two objects and update the pair vector accordingly
Definition: Pair2D_extra.cpp:292
cbl::pairs::Pair2D_comovingPolar_loglog_extra::put
void put(const std::shared_ptr< catalogue::Object > obj1, const std::shared_ptr< catalogue::Object > obj2) override
estimate the distance between two objects and update the pair vector accordingly
Definition: Pair2D_extra.cpp:333
cbl::pairs::Pair2D_extra::add_data2D
void add_data2D(const int i, const int j, const std::vector< double > data) override
set the protected members by adding new data
Definition: Pair2D_extra.cpp:374
cbl::pairs::Pair2D_extra::Sum
void Sum(const std::shared_ptr< Pair > pair, const double ww=1) override
sum the number of binned pairs
Definition: Pair2D_extra.cpp:444
cbl::par::pi
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
cbl::Euclidean_distance
double Euclidean_distance(const double x1, const double x2, const double y1, const double y2, const double z1, const double z2)
the Euclidean distance in 3D relative to the origin (0,0,0), i.e. the Euclidean norm
Definition: Func.cpp:250
cbl::converted_angle
double converted_angle(const double angle, const CoordinateUnits inputUnits=CoordinateUnits::_radians_, const CoordinateUnits outputUnits=CoordinateUnits::_degrees_)
conversion to angle units
Definition: Func.cpp:192
cbl::angular_distance
double angular_distance(const double x1, const double x2, const double y1, const double y2, const double z1, const double z2)
the angular separation in 3D
Definition: Func.cpp:277
cbl::ErrorCBL
int ErrorCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_)
throw an exception: it is used for handling exceptions inside the CosmoBolognaLib
Definition: Kernel.h:780
cbl::perpendicular_distance
double perpendicular_distance(const double ra1, const double ra2, const double dec1, const double dec2, const double d1, const double d2)
the perpendicular separation, rp
Definition: Func.cpp:259