CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Pair1D_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 "Pair1D_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::Pair1D_angular_lin_extra::put (const shared_ptr<Object> obj1, const shared_ptr<Object> obj2)
48 {
49  const double dist = (m_angularUnits==CoordinateUnits::_radians_) ? angular_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz())
50  : converted_angle(angular_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz()), CoordinateUnits::_radians_, m_angularUnits);
51 
52  if (m_thetaMin < dist && dist < m_thetaMax) {
53 
54  const int kk = max(0, min(int((dist-m_thetaMin)*m_binSize_inv), m_nbins));
55 
56  const double WeightTOT = obj1->weight()*obj2->weight();
57 
58  m_PP1D[kk] ++;
59  m_PP1D_weighted[kk] += WeightTOT;
60 
61  if (m_PP1D_weighted[kk]>0) {
62 
63  const double scale_mean_p = m_scale_mean[kk];
64  m_scale_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
65  m_scale_S[kk] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
66 
67  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
68  const double z_mean_p = m_z_mean[kk];
69  m_z_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(pair_redshift-z_mean_p);
70  m_z_S[kk] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[kk]);
71 
72  }
73 
74  }
75 }
76 
77 
78 // ============================================================================================
79 
80 
81 void cbl::pairs::Pair1D_angular_log_extra::put (const shared_ptr<Object> obj1, const shared_ptr<Object> obj2)
82 {
83  const double dist = (m_angularUnits==CoordinateUnits::_radians_) ? angular_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz())
84  : converted_angle(angular_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz()), CoordinateUnits::_radians_, m_angularUnits);
85 
86  if (m_thetaMin < dist && dist < m_thetaMax) {
87 
88  const int kk = max(0, min(int((log10(dist)-log10(m_thetaMin))*m_binSize_inv), m_nbins));
89 
90  const double WeightTOT = obj1->weight()*obj2->weight();
91 
92  m_PP1D[kk] ++;
93  m_PP1D_weighted[kk] += WeightTOT;
94 
95  if (m_PP1D_weighted[kk]>0) {
96 
97  const double scale_mean_p = m_scale_mean[kk];
98  m_scale_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
99  m_scale_S[kk] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
100 
101  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
102  const double z_mean_p = m_z_mean[kk];
103  m_z_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(pair_redshift-z_mean_p);
104  m_z_S[kk] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[kk]);
105 
106  }
107 
108  }
109 }
110 
111 
112 // ============================================================================================
113 
114 
115 void cbl::pairs::Pair1D_comoving_lin_extra::put (const shared_ptr<Object> obj1, const shared_ptr<Object> obj2)
116 {
117  const double dist = Euclidean_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz());
118 
119  if (m_rMin < dist && dist < m_rMax) {
120 
121  const int kk = max(0, min(int((dist-m_rMin)*m_binSize_inv), m_nbins));
122 
123  const double angWeight = (m_angularWeight==nullptr) ? 1.
124  : 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)));
125 
126  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
127 
128  m_PP1D[kk] ++;
129  m_PP1D_weighted[kk] += WeightTOT;
130 
131  if (m_PP1D_weighted[kk]>0) {
132 
133  const double scale_mean_p = m_scale_mean[kk];
134  m_scale_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
135  m_scale_S[kk] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
136 
137  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
138  const double z_mean_p = m_z_mean[kk];
139  m_z_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(pair_redshift-z_mean_p);
140  m_z_S[kk] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[kk]);
141 
142  }
143 
144  }
145 }
146 
147 
148 // ============================================================================================
149 
150 
151 void cbl::pairs::Pair1D_comoving_log_extra::put (const shared_ptr<Object> obj1, const shared_ptr<Object> obj2)
152 {
153  const double dist = Euclidean_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz());
154 
155  if (m_rMin < dist && dist < m_rMax) {
156 
157  const int kk = max(0, min(int((log10(dist)-log10(m_rMin))*m_binSize_inv), m_nbins));
158 
159  const double angWeight = (m_angularWeight==nullptr) ? 1.
160  : 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)));
161 
162  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
163 
164  m_PP1D[kk] ++;
165  m_PP1D_weighted[kk] += WeightTOT;
166 
167  if (m_PP1D_weighted[kk]>0) {
168 
169  const double scale_mean_p = m_scale_mean[kk];
170  m_scale_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
171  m_scale_S[kk] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
172 
173  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
174  const double z_mean_p = m_z_mean[kk];
175  m_z_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(pair_redshift-z_mean_p);
176  m_z_S[kk] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[kk]);
177 
178  }
179 
180  }
181 }
182 
183 
184 // ============================================================================================
185 
186 
187 void cbl::pairs::Pair1D_comoving_multipoles_lin_extra::put (const shared_ptr<Object> obj1, const shared_ptr<Object> obj2)
188 {
189  const double dist = Euclidean_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz());
190 
191  if (m_rMin < dist && dist < m_rMax) {
192 
193  const int kk = max(0, min(int((dist-m_rMin)*m_binSize_inv), m_nbins));
194 
195  const double angWeight = (m_angularWeight==nullptr) ? 1.
196  : 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)));
197 
198  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
199  const double cosmu2 = pow((obj2->dc()-obj1->dc())/dist, 2);
200 
201  m_PP1D[kk] ++;
202  m_PP1D_weighted[kk] += WeightTOT;
203 
204  double leg_pol_2 = 0.5*(3.*cosmu2-1.);
205  m_PP1D[(m_nbins+1)+kk] += 5.*leg_pol_2 ;
206  m_PP1D_weighted[(m_nbins+1)+kk] += 5.*WeightTOT*leg_pol_2;
207 
208  double leg_pol_4 = 0.125*(35.*cosmu2*cosmu2-30.*cosmu2+3.);
209  m_PP1D[2.*(m_nbins+1)+kk] += 9.*leg_pol_4;
210  m_PP1D_weighted[2.*(m_nbins+1)+kk] += 9.*WeightTOT*leg_pol_4;
211 
212  if (m_PP1D_weighted[kk]>0) {
213 
214  const double scale_mean_p = m_scale_mean[kk];
215  m_scale_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
216  m_scale_S[kk] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
217 
218  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
219  const double z_mean_p = m_z_mean[kk];
220  m_z_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(pair_redshift-z_mean_p);
221  m_z_S[kk] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[kk]);
222 
223  m_scale_mean[kk+(m_nbins+1)] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
224  m_scale_S[kk+(m_nbins+1)] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
225 
226  m_scale_mean[kk+2*(m_nbins+1)] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
227  m_scale_S[kk+2*(m_nbins+1)] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
228  }
229  }
230 }
231 
232 
233 // ============================================================================================
234 
235 
236 void cbl::pairs::Pair1D_comoving_multipoles_log_extra::put (const shared_ptr<Object> obj1, const shared_ptr<Object> obj2)
237 {
238  const double dist = Euclidean_distance(obj1->xx(), obj2->xx(), obj1->yy(), obj2->yy(), obj1->zz(), obj2->zz());
239 
240  if (m_rMin < dist && dist < m_rMax) {
241 
242  const int kk = max(0, min(int((log10(dist)-log10(m_rMin))*m_binSize_inv), m_nbins));
243 
244  const double angWeight = (m_angularWeight==nullptr) ? 1.
245  : 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)));
246 
247  const double WeightTOT = obj1->weight()*obj2->weight()*angWeight;
248  const double cosmu2 = pow((obj2->dc()-obj1->dc())/dist, 2);
249 
250  m_PP1D[kk] ++;
251  m_PP1D_weighted[kk] += WeightTOT;
252 
253  double leg_pol_2 = 0.5*(3.*cosmu2-1.);
254  m_PP1D[(m_nbins+1)+kk] += 5.*leg_pol_2 ;
255  m_PP1D_weighted[(m_nbins+1)+kk] += 5.*WeightTOT*leg_pol_2;
256 
257  double leg_pol_4 = 0.125*(35.*cosmu2*cosmu2-30.*cosmu2+3.);
258  m_PP1D[2.*(m_nbins+1)+kk] += 9.*leg_pol_4;
259  m_PP1D_weighted[2.*(m_nbins+1)+kk] += 9.*WeightTOT*leg_pol_4;
260 
261  if (m_PP1D_weighted[kk]>0) {
262 
263  const double scale_mean_p = m_scale_mean[kk];
264  m_scale_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
265  m_scale_S[kk] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
266 
267  const double pair_redshift = (obj1->redshift()>0 && obj2->redshift()>0) ? (obj1->redshift()+obj2->redshift())*0.5 : -1.;
268  const double z_mean_p = m_z_mean[kk];
269  m_z_mean[kk] += WeightTOT/m_PP1D_weighted[kk]*(pair_redshift-z_mean_p);
270  m_z_S[kk] += WeightTOT*(pair_redshift-z_mean_p)*(pair_redshift-m_z_mean[kk]);
271 
272  m_scale_mean[kk+(m_nbins+1)] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
273  m_scale_S[kk+(m_nbins+1)] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
274 
275  m_scale_mean[kk+2*(m_nbins+1)] += WeightTOT/m_PP1D_weighted[kk]*(dist-scale_mean_p);
276  m_scale_S[kk+2*(m_nbins+1)] += WeightTOT*(dist-scale_mean_p)*(dist-m_scale_mean[kk]);
277  }
278 
279  }
280 }
281 
282 
283 // ============================================================================================
284 
285 
286 void cbl::pairs::Pair1D_extra::add_data1D (const int i, const std::vector<double> data)
287 {
288  /*
289  checkDim(m_PP1D, i, "m_PP1D", false);
290  checkDim(m_PP1D_weighted, i, "m_PP1D_weighted", false);
291  checkDim(m_scale_mean, i, "m_scale_mean", false);
292  checkDim(m_scale_sigma, i, "m_scale_sigma", false);
293  checkDim(m_z_mean, i, "m_z_mean", false);
294  checkDim(m_z_sigma, i, "m_z_sigma", false);
295  checkDim(data, 6, "data");
296  */
297 
298  // mean scale up to this computation step
299  const double scale_mean_temp_p = m_scale_mean[i];
300 
301  // mean redshift up to this computation step
302  const double z_mean_temp_p = m_z_mean[i];
303 
304  // number of pairs
305  m_PP1D[i] += data[0];
306 
307  // number of weighted pairs
308  m_PP1D_weighted[i] += data[1];
309 
310 
311  if (m_PP1D_weighted[i]>0) {
312 
313  // mean separation
314  m_scale_mean[i] += data[1]/m_PP1D_weighted[i]*(data[2]-scale_mean_temp_p);
315 
316  // mean redshift
317  m_z_mean[i] += data[1]/m_PP1D_weighted[i]*(data[4]-z_mean_temp_p);
318 
319  // compute the weighted standard deviation of the scale distribution
320  m_scale_S[i] += data[3]+pow(data[2]-scale_mean_temp_p, 2)*data[1]*(m_PP1D_weighted[i]-data[1])/m_PP1D_weighted[i];
321  m_scale_sigma[i] = sqrt(m_scale_S[i]/m_PP1D_weighted[i]);
322 
323  // compute the weighted standard deviation of the redshift distribution
324  m_z_S[i] += data[5]+pow(data[4]-z_mean_temp_p, 2)*data[1]*(m_PP1D_weighted[i]-data[1])/m_PP1D_weighted[i];
325  m_z_sigma[i] = sqrt(m_z_S[i]/m_PP1D_weighted[i]);
326 
327  }
328 
329 }
330 
331 
332 // ============================================================================================
333 
334 
335 void cbl::pairs::Pair1D_extra::add_data1D (const int i, const std::shared_ptr<pairs::Pair> pair, const double ww)
336 {
337  add_data1D(i, {ww*pair->PP1D(i), ww*pair->PP1D_weighted(i), pair->scale_mean(i), pair->scale_S(i), pair->z_mean(i), pair->z_S(i)});
338 }
339 
340 
341 // ============================================================================================
342 
343 
344 void cbl::pairs::Pair1D_extra::Sum (const std::shared_ptr<Pair> pair, const double ww)
345 {
346  if (m_nbins != pair->nbins())
347  ErrorCBL("dimension problems!", "Sum", "Pair1D_extra.cpp");
348 
349  for (int i=0; i<m_nbins; ++i)
350  add_data1D(i, pair, ww);
351 }
352 
353 
354 // ============================================================================================
355 
356 
357 void cbl::pairs::Pair1D_comoving_multipoles_lin_extra::Sum (const std::shared_ptr<Pair> pair, const double ww)
358 {
359  if (m_nbins != pair->nbins())
360  ErrorCBL("dimension problems!", "Sum", "Pair1D_extra.cpp");
361 
362  for (int l=0; l<3; ++l)
363  for (int i=0; i<m_nbins; ++i)
364  add_data1D(l*(m_nbins+1)+i, pair, ww);
365 }
366 
367 
368 // ============================================================================================
369 
370 
371 void cbl::pairs::Pair1D_comoving_multipoles_log_extra::Sum (const std::shared_ptr<Pair> pair, const double ww)
372 {
373  if (m_nbins != pair->nbins())
374  ErrorCBL("dimension problems!", "Sum", "Pair1D_extra.cpp");
375 
376  for (int l=0; l<3; ++l)
377  for (int i=0; i<m_nbins; ++i)
378  add_data1D(l*(m_nbins+1)+i, pair, ww);
379 }
Pair1D_extra.h
The classes Pair1D_extra*.
cbl::pairs::Pair1D_angular_lin_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: Pair1D_extra.cpp:47
cbl::pairs::Pair1D_angular_log_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: Pair1D_extra.cpp:81
cbl::pairs::Pair1D_comoving_lin_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: Pair1D_extra.cpp:115
cbl::pairs::Pair1D_comoving_log_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: Pair1D_extra.cpp:151
cbl::pairs::Pair1D_comoving_multipoles_lin_extra::Sum
virtual void Sum(const std::shared_ptr< Pair > pair, const double ww=1) override
sum the number of binned pairs
Definition: Pair1D_extra.cpp:357
cbl::pairs::Pair1D_comoving_multipoles_lin_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: Pair1D_extra.cpp:187
cbl::pairs::Pair1D_comoving_multipoles_log_extra::Sum
virtual void Sum(const std::shared_ptr< Pair > pair, const double ww=1) override
sum the number of binned pairs
Definition: Pair1D_extra.cpp:371
cbl::pairs::Pair1D_comoving_multipoles_log_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: Pair1D_extra.cpp:236
cbl::pairs::Pair1D_extra::Sum
virtual void Sum(const std::shared_ptr< Pair > pair, const double ww=1) override
sum the number of binned pairs
Definition: Pair1D_extra.cpp:344
cbl::pairs::Pair1D_extra::add_data1D
void add_data1D(const int i, const std::vector< double > data) override
set the protected members by adding new data
Definition: Pair1D_extra.cpp:286
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