CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Catalogue.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2015 by Federico Marulli and Alfonso Veropalumbo *
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 
34 #include "Field3D.h"
35 #include "Catalogue.h"
36 #include "Object.h"
37 
38 using namespace std;
39 
40 using namespace cbl;
41 using namespace catalogue;
42 using namespace chainmesh;
43 
44 
45 // ============================================================================
46 
48 
49 template cbl::catalogue::Catalogue::Catalogue (vector<cbl::catalogue::RandomObject>);
50 template cbl::catalogue::Catalogue::Catalogue (vector<cbl::catalogue::Mock>);
51 template cbl::catalogue::Catalogue::Catalogue (vector<cbl::catalogue::Halo>);
52 template cbl::catalogue::Catalogue::Catalogue (vector<cbl::catalogue::Galaxy>);
53 template cbl::catalogue::Catalogue::Catalogue (vector<cbl::catalogue::Cluster>);
54 template cbl::catalogue::Catalogue::Catalogue (vector<cbl::catalogue::Void>);
55 template cbl::catalogue::Catalogue::Catalogue (vector<cbl::catalogue::HostHalo>);
56 template cbl::catalogue::Catalogue::Catalogue (vector<cbl::catalogue::ChainMeshCell>);
57 
58 template void cbl::catalogue::Catalogue::add_object (cbl::catalogue::RandomObject);
59 template void cbl::catalogue::Catalogue::add_object (cbl::catalogue::Mock);
60 template void cbl::catalogue::Catalogue::add_object (cbl::catalogue::Halo);
61 template void cbl::catalogue::Catalogue::add_object (cbl::catalogue::Galaxy);
62 template void cbl::catalogue::Catalogue::add_object (cbl::catalogue::Cluster);
63 template void cbl::catalogue::Catalogue::add_object (cbl::catalogue::Void);
64 template void cbl::catalogue::Catalogue::add_object (cbl::catalogue::HostHalo);
65 template void cbl::catalogue::Catalogue::add_object (cbl::catalogue::ChainMeshCell);
66 
67 template void cbl::catalogue::Catalogue::add_objects (vector<cbl::catalogue::RandomObject>);
68 template void cbl::catalogue::Catalogue::add_objects (vector<cbl::catalogue::Mock>);
69 template void cbl::catalogue::Catalogue::add_objects (vector<cbl::catalogue::Halo>);
70 template void cbl::catalogue::Catalogue::add_objects (vector<cbl::catalogue::Galaxy>);
71 template void cbl::catalogue::Catalogue::add_objects (vector<cbl::catalogue::Cluster>);
72 template void cbl::catalogue::Catalogue::add_objects (vector<cbl::catalogue::Void>);
73 template void cbl::catalogue::Catalogue::add_objects (vector<cbl::catalogue::HostHalo>);
74 template void cbl::catalogue::Catalogue::add_objects (vector<cbl::catalogue::ChainMeshCell>);
75 
76 template void cbl::catalogue::Catalogue::replace_objects (vector<cbl::catalogue::RandomObject>);
77 template void cbl::catalogue::Catalogue::replace_objects (vector<cbl::catalogue::Mock>);
78 template void cbl::catalogue::Catalogue::replace_objects (vector<cbl::catalogue::Halo>);
79 template void cbl::catalogue::Catalogue::replace_objects (vector<cbl::catalogue::Galaxy>);
80 template void cbl::catalogue::Catalogue::replace_objects (vector<cbl::catalogue::Cluster>);
81 template void cbl::catalogue::Catalogue::replace_objects (vector<cbl::catalogue::Void>);
82 template void cbl::catalogue::Catalogue::replace_objects (vector<cbl::catalogue::HostHalo>);
83 template void cbl::catalogue::Catalogue::replace_objects (vector<cbl::catalogue::ChainMeshCell>);
84 
86 
87 
88 // ============================================================================
89 
90 
91 cbl::catalogue::Catalogue::Catalogue (const ObjectType objectType, const CoordinateType coordinateType, const std::vector<double> coord1, const std::vector<double> coord2, const std::vector<double> coord3, const std::vector<double> weight, const cosmology::Cosmology &cosm, const CoordinateUnits inputUnits)
92 {
93  // check the vector dimensions
94  if (!(coord1.size()==coord2.size() && coord2.size()==coord3.size()))
95  ErrorCBL("coordinates with different dimensions!", "Catalogue", "Catalogue.cpp");
96 
97  // weight[]=1, if are not used
98  vector<double> _weight = weight;
99  if (_weight.size()==0) _weight.resize(coord1.size(), 1);
100 
101 
102  // include the objects in the catalogue
103 
104  for (size_t i=0; i<coord1.size(); ++i) {
105 
106  if (coordinateType==cbl::CoordinateType::_comoving_) { // comoving coordinates (x, y, z)
107  comovingCoordinates coord = {coord1[i], coord2[i], coord3[i]};
108  m_object.push_back(move(Object::Create(objectType, coord, _weight[i])));
109  }
110  else if (coordinateType==cbl::CoordinateType::_observed_) { // observed coordinates (R.A., Dec, redshift)
111  observedCoordinates coord = {coord1[i], coord2[i], coord3[i]};
112  m_object.push_back(move(Object::Create(objectType, coord, inputUnits, cosm, _weight[i])));
113  }
114  else ErrorCBL("CoordinateType is not valid!", "Catalogue", "Catalogue.cpp");
115 
116  }
117 
118 }
119 
120 
121 // ============================================================================
122 
123 
124 cbl::catalogue::Catalogue::Catalogue (const ObjectType objectType, const CoordinateType coordinateType, const std::vector<std::string> file, const int col1, const int col2, const int col3, const int colWeight, const int colRegion, const double nSub, const double fact, const cosmology::Cosmology &cosm, const CoordinateUnits inputUnits, const CharEncode charEncode, const std::string comment, const int seed)
125 {
126  // parameters for random numbers used in case nSub!=1
127  random::UniformRandomNumbers ran(0., 1., seed);
128 
129  // read the input catalogue files
130 
131  for (size_t dd=0; dd<file.size(); ++dd) {
132 
133  string line, file_in = file[dd];
134 
135  if (charEncode==CharEncode::_ascii_) {
136 
137  coutCBL << "Reading the catalogue: " << file_in << endl;
138  ifstream finr(file_in.c_str()); checkIO(finr, file_in);
139 
140  double Weight, Value; long Region;
141 
142  long maxRegion = -100000000; // check!
143 
144  while (getline(finr, line)) { // read the lines
145 
146  if (line.find(comment)==string::npos) { // skip a comment
147 
148  if (ran()<nSub) { // extract a subsample
149 
150  stringstream ss(line);
151  vector<double> value; while (ss>>Value) value.emplace_back(Value);
152  checkDim(value, col1, "value", false); checkDim(value, col2, "value", false);
153 
154  if (coordinateType==cbl::CoordinateType::_comoving_) { // comoving coordinates (x, y, z)
155  checkDim(value, col3, "value", false);
156  comovingCoordinates coord;
157  coord.xx = value[col1-1]*fact;
158  coord.yy = value[col2-1]*fact;
159  coord.zz = value[col3-1]*fact;
160  Weight = (colWeight!=-1 && colWeight-1<(int)value.size()) ? value[colWeight-1] : 1.;
161  Region = (colRegion!=-1 && colRegion-1<(int)value.size()) ? (long)value[colRegion-1] : 0;
162  m_object.push_back(move(Object::Create(objectType, coord, Weight, Region)));
163  }
164 
165  else if (coordinateType==cbl::CoordinateType::_observed_) { // observed coordinates (R.A., Dec (redshift))
166  observedCoordinates coord;
167  coord.ra = value[col1-1]*fact;
168  coord.dec = value[col2-1]*fact;
169  coord.redshift = ((int)value.size()>=col3) ? value[col3-1] : 1.;
170  Weight = (colWeight!=-1 && colWeight-1<(int)value.size()) ? value[colWeight-1] : 1.;
171  Region = (colRegion!=-1 && colRegion-1<(int)value.size()) ? (long)value[colRegion-1] : 0;
172  m_object.push_back(move(Object::Create(objectType, coord, inputUnits, cosm, Weight, Region)));
173  }
174 
175  else ErrorCBL("CoordinateType is not valid!", "Catalogue", "Catalogue.cpp");
176 
177  maxRegion = max(maxRegion, Region);
178  }
179  }
180  }
181 
182  if (colRegion>0) {
183  m_nRegions = maxRegion+1;
184  WarningMsgCBL("The total number of region is "+conv(m_nRegions, par::fINT)+", deducted from the maximum value of the input regions; if needed, this number can be changed with the function set_region_number()", "Catalogue", "Catalogue.cpp");
185  }
186 
187  finr.clear(); finr.close();
188  }
189 
190  else if (charEncode==CharEncode::_binary_) {
191 
192  // read the input catalogue files
193 
194  coutCBL << "Reading the catalogue: " << file_in << endl;
195 
196  short num_bin;
197  float val;
198 
199  ifstream finr(file_in.c_str(), ios::in|ios::binary|ios::ate); checkIO(finr, file_in);
200 
201  comovingCoordinates coord;
202 
203  if (finr.is_open())
204  finr.seekg(0, ios::beg);
205 
206  finr.read((char*)(&num_bin), 2);
207 
208  int n_blocks = num_bin;
209 
210  for (int i=1; i<=n_blocks; ++i) {
211  finr.read((char*)(&num_bin), 4);
212  int n_objs = num_bin;
213 
214  for (int j=1; j<=n_objs; ++j) {
215  finr.read((char*)(&val), 4);
216  coord.xx = (val)*fact;
217  finr.read((char*)(&val), 4);
218  coord.yy = (val)*fact;
219  finr.read((char*)(&val), 4);
220  coord.zz = (val)*fact;
221  // Weight = (colWeight!=-1 && colWeight-1<value.size()) ? value[colWeight-1] : 1.;
222  // Region = (colRegion!=-1 && colRegion-1<value.size()) ? (long)value[colRegion-1] : 0;
223  if (ran()<nSub) m_object.push_back(move(Object::Create(objectType, coord)));
224  // if (ran()<nSub) m_object.push_back(move(Object::Create(objectType, coord, Weight, Region)));
225  }
226 
227  finr.read((char*)(&num_bin), 4);
228 
229  int n_objs2 = num_bin;
230 
231  if (n_objs2!=n_objs) ErrorCBL("wrong reading of input binary file", "Catalogue", "Catalogue.cpp");
232  }
233 
234  finr.clear(); finr.close();
235  }
236 
237  else ErrorCBL("charEncode is not valid!", "Catalogue", "Catalogue.cpp");
238  }
239 }
240 
241 
242 // ============================================================================
243 
244 
245 cbl::catalogue::Catalogue::Catalogue (const ObjectType objectType, const CoordinateType coordinateType, const std::vector<Var> attribute, const std::vector<int> column, const std::vector<std::string> file, const int comments, const double nSub, const double fact, const cosmology::Cosmology &cosm, const CoordinateUnits inputUnits, const char delimiter, const int seed)
246 {
247 
248  // preliminary check on vector sizes
249  size_t nvar;
250  if (attribute.size()==column.size()) nvar = attribute.size();
251  else ErrorCBL("Column vector and attribute vector must have equal size!", "Catalogue", "Catalogue.cpp");
252 
253  if (!(std::is_sorted(column.begin(), column.end()))) ErrorCBL("Column vector must be sort in ascending order!", "Catalogue", "Catalogue.cpp");
254 
255  // std::unordered_map to map columns and attributes
256  unordered_map<int, Var> varMap;
257  for (size_t ii=0; ii<nvar; ii++)
258  varMap.insert({column[ii], attribute[ii]});
259 
260  // parameters for random numbers used in case nSub!=1
261  random::UniformRandomNumbers ran(0., 1., seed);
262 
263  // read the input catalogue files
264 
265  for (size_t dd=0; dd<file.size(); ++dd) {
266 
267  string line, file_in = file[dd];
268 
269  coutCBL << "Reading the catalogue: " << file_in << endl;
270  ifstream finr(file_in.c_str()); checkIO(finr, file_in);
271 
272  // prepare default coordinates
273  comovingCoordinates defaultComovingCoord = {par::defaultDouble, par::defaultDouble, par::defaultDouble};
274  observedCoordinates defaultObservedCoord = {par::defaultDouble, -1., 0.1};
275 
276  for (int cc=0; cc<comments; cc++) getline(finr, line); // ignore commented lines at the beginning of file
277 
278  while (getline(finr, line)) { // read the lines
279 
280  if (ran()<nSub) { // extract a subsample
281 
282  if (coordinateType==cbl::CoordinateType::_comoving_)
283  m_object.push_back(move(Object::Create(objectType, defaultComovingCoord, 1.)));
284 
285  else if (coordinateType==cbl::CoordinateType::_observed_)
286  m_object.push_back(move(Object::Create(objectType, defaultObservedCoord, inputUnits, cosm, 1.)));
287 
288  else ErrorCBL("CoordinateType is not valid!", "Catalogue", "Catalogue.cpp");
289 
290  stringstream ss(line);
291 
292  if (delimiter!='\t') {
293  string temp_string;
294  stringstream temp(line);
295  ss.clear();
296  while (getline(temp, temp_string, delimiter))
297  ss << temp_string+" ";
298  }
299 
300  double Value_d;
301  int Value_i;
302  //std::string Value_s;
303  size_t ii = nObjects()-1;
304  int index = 0;
305 
306  for (int jj=1; jj<=cbl::Max(column); jj++) {
307  if (varMap[column[index]]==Var::_ID_ || varMap[column[index]]==Var::_IDHOST_) {
308 
309  ss>>Value_i;
310  if (std::find(column.begin(), column.end(), jj)!=column.end()) {
311  set_var(ii, varMap[column[index]], Value_i);
312  index++;
313 
314  }
315 
316  }
317  else if (varMap[column[index]]==Var::_GalaxyTag_) {
318  ss>>Value_d;
319  if (std::find(column.begin(), column.end(), jj)!=column.end()) {
320  set_var(ii, varMap[column[index]], Value_d);
321  index++;
322 
323  }
324  }
325  else {
326  ss>>Value_d;
327  if (std::find(column.begin(), column.end(), jj)!=column.end()) {
328  if ((varMap[column[index]]==Var::_RA_) || (varMap[column[index]]==Var::_Dec_)) Value_d = radians(Value_d, inputUnits);
329  set_var(ii, varMap[column[index]],
330  ((varMap[column[index]]==Var::_X_) || (varMap[column[index]]==Var::_Y_) || (varMap[column[index]]==Var::_Z_)) ?
331  Value_d*fact : Value_d,
332  cosm);
333  index++;
334  }
335  }
336 
337  }
338  }
339  }
340  finr.clear(); finr.close();
341  }
342 }
343 
344 // ============================================================================
345 
346 
347 cbl::catalogue::Catalogue::Catalogue (const ObjectType objectType, const std::vector<Var> attribute, const std::vector<int> column, const std::vector<std::string> file, const int comments, const double nSub, const int seed)
348 {
349  // parameters for random numbers used in case nSub!=1
350  random::UniformRandomNumbers ran(0., 1., seed);
351 
352  // read the input catalogue files
353 
354  vector<vector<double>> vars;
355 
356  for (size_t dd=0; dd<file.size(); ++dd) {
357 
358  string line, file_in = file[dd];
359 
360  coutCBL << "Reading the catalogue: " << file_in << endl;
361  ifstream finr(file_in.c_str()); checkIO(finr, file_in);
362 
363  // start reading catalogue
364  for (int cc=0; cc<comments; cc++) getline(finr, line); // ignore commented lines at the beginning of file
365  while (getline(finr, line)) { // read the lines
366 
367  if (ran()<nSub) { // extract a subsample
368 
369  stringstream ss(line);
370  vector<double> num; double NN = -1.e30;
371  while (ss>>NN) num.push_back(NN);
372 
373  m_object.push_back(move(Object::Create(objectType)));
374 
375  vector<double> vv;
376  for (size_t i=0; i<column.size(); i++)
377  vv.push_back(num[column[i]-1]);
378 
379  vars.push_back(vv);
380  }
381  }
382 
383  finr.clear(); finr.close();
384  }
385 
386  vars = transpose(vars);
387 
388  for (size_t i=0; i<vars.size(); i++)
389  set_var(attribute[i], vars[i]);
390 }
391 
392 
393 // ============================================================================
394 
395 
396 size_t cbl::catalogue::Catalogue::nRegions ()
397 {
398  if (m_nRegions==0) {
399  m_nRegions = Max(Var::_Region_)+1;
400 
401  WarningMsgCBL("The total number of region is "+conv(m_nRegions, par::fINT)+", deducted from the maximum value of the input regions; if needed, this number can be changed with the function set_region_number()", "Catalogue", "Catalogue.cpp");
402  }
403 
404  return m_nRegions;
405 }
406 
407 
408 // ============================================================================
409 
410 
411 std::vector<long> cbl::catalogue::Catalogue::region_list () const
412 {
413  vector<long> vv(m_nRegions);
414 
415  for (size_t i=0; i<m_nRegions; ++i) vv[i] = static_cast<long>(i);
416 
417  return vv;
418 }
419 
420 // ============================================================================
421 
422 
423 std::vector<long> cbl::catalogue::Catalogue::region () const
424 {
425  vector<long> vv(m_object.size());
426 
427  for (size_t i=0; i<nObjects(); ++i) vv[i] = m_object[i]->region();
428 
429  return vv;
430 }
431 
432 
433 // ============================================================================
434 
435 
436 std::vector<std::string> cbl::catalogue::Catalogue::field () const
437 {
438  vector<string> vv(m_object.size());
439 
440  for (size_t i=0; i<nObjects(); ++i) vv[i] = m_object[i]->field();
441 
442  return vv;
443 }
444 
445 
446 // ============================================================================
447 
448 
449 double cbl::catalogue::Catalogue::var (int index, Var var_name) const
450 {
451  double vv = 0;
452 
453  switch (var_name) {
454 
455  case Var::_X_:
456  vv = m_object[index]->xx();
457  break;
458 
459  case Var::_Y_:
460  vv = m_object[index]->yy();
461  break;
462 
463  case Var::_Z_:
464  vv = m_object[index]->zz();
465  break;
466 
467  case Var::_RA_:
468  vv = m_object[index]->ra();
469  break;
470 
471  case Var::_Dec_:
472  vv = m_object[index]->dec();
473  break;
474 
475  case Var::_TileRA_:
476  vv = m_object[index]->ra_tile();
477  break;
478 
479  case Var::_TileDec_:
480  vv = m_object[index]->dec_tile();
481  break;
482 
483  case Var::_SN_:
484  vv = m_object[index]->sn();
485  break;
486 
487  case Var::_Redshift_:
488  vv = m_object[index]->redshift();
489  break;
490 
491  case Var::_RedshiftMin_:
492  vv = m_object[index]->redshiftMin();
493  break;
494 
495  case Var::_RedshiftMax_:
496  vv = m_object[index]->redshiftMax();
497  break;
498 
499  case Var::_Shear1_:
500  vv = m_object[index]->shear1();
501  break;
502 
503  case Var::_Shear2_:
504  vv = m_object[index]->shear2();
505  break;
506 
507  case Var::_ODDS_:
508  vv = m_object[index]->odds();
509  break;
510 
511  case Var::_LensingWeight_:
512  vv = m_object[index]->lensingWeight();
513  break;
514 
515  case Var::_LensingCalib_:
516  vv = m_object[index]->lensingCalib();
517  break;
518 
519  case Var::_Dc_:
520  vv = m_object[index]->dc();
521  break;
522 
523  case Var::_Weight_:
524  vv = m_object[index]->weight();
525  break;
526 
527  case Var::_Mass_:
528  vv = m_object[index]->mass();
529  break;
530 
531  case Var::_Magnitude_:
532  vv = m_object[index]->magnitude();
533  break;
534 
535  case Var::_MagnitudeU_:
536  vv = m_object[index]->magnitudeU();
537  break;
538 
539  case Var::_MagnitudeG_:
540  vv = m_object[index]->magnitudeG();
541  break;
542 
543  case Var::_MagnitudeR_:
544  vv = m_object[index]->magnitudeR();
545  break;
546 
547  case Var::_MagnitudeI_:
548  vv = m_object[index]->magnitudeI();
549  break;
550 
551  case Var::_SFR_:
552  vv = m_object[index]->SFR();
553  break;
554 
555  case Var::_sSFR_:
556  vv = m_object[index]->sSFR();
557  break;
558 
559  case Var::_MassProxy_:
560  vv = m_object[index]->mass_proxy();
561  break;
562 
563  case Var::_MassProxyError_:
564  vv = m_object[index]->mass_proxy_error();
565  break;
566 
567  case Var::_Mstar_:
568  vv = m_object[index]->mstar();
569  break;
570 
571  case Var::_MassInfall_:
572  vv = m_object[index]->massinfall();
573  break;
574 
575  case Var::_Vx_:
576  vv = m_object[index]->vx();
577  break;
578 
579  case Var::_Vy_:
580  vv = m_object[index]->vy();
581  break;
582 
583  case Var::_Vz_:
584  vv = m_object[index]->vz();
585  break;
586 
587  case Var::_Region_:
588  vv = m_object[index]->region();
589  break;
590 
591  case Var::_Generic_:
592  vv = m_object[index]->generic();
593  break;
594 
595  case Var::_Radius_:
596  vv = m_object[index]->radius();
597  break;
598 
599  case Var::_DensityContrast_:
600  vv = m_object[index]->densityContrast();
601  break;
602 
603  case Var::_CentralDensity_:
604  vv = m_object[index]->centralDensity();
605  break;
606 
607  case Var::_X_displacement_:
608  vv = m_object[index]->x_displacement();
609  break;
610 
611  case Var::_Y_displacement_:
612  vv = m_object[index]->y_displacement();
613  break;
614 
615  case Var::_Z_displacement_:
616  vv = m_object[index]->z_displacement();
617  break;
618 
619  case Var::_MassEstimate_:
620  vv = m_object[index]->mass_estimate();
621  break;
622 
623  case Var::_RadiusEstimate_:
624  vv = m_object[index]->radius_estimate();
625  break;
626 
627  case Var::_VeldispEstimate_:
628  vv = m_object[index]->veldisp_estimate();
629  break;
630 
631  case Var::_XCM_:
632  vv = m_object[index]->xcm();
633  break;
634 
635  case Var::_YCM_:
636  vv = m_object[index]->ycm();
637  break;
638 
639  case Var::_ZCM_:
640  vv = m_object[index]->zcm();
641  break;
642 
643  case Var::_XSpin_:
644  vv = m_object[index]->spin_x();
645  break;
646 
647  case Var::_YSpin_:
648  vv = m_object[index]->spin_y();
649  break;
650 
651  case Var::_ZSpin_:
652  vv = m_object[index]->spin_z();
653  break;
654 
655  case Var::_VelDisp_:
656  vv = m_object[index]->veldisp();
657  break;
658 
659  case Var::_Vmax_:
660  vv = m_object[index]->vmax();
661  break;
662 
663  case Var::_VmaxRad_:
664  vv = m_object[index]->vmax_rad();
665  break;
666 
667  case Var::_TotMass_:
668  vv = m_object[index]->tot_mass();
669  break;
670 
671  case Var::_ID_:
672  vv = m_object[index]->ID();
673  break;
674 
675  case Var::_Nsub_:
676  vv = m_object[index]->nsub();
677  break;
678 
679  case Var::_Parent_:
680  vv = m_object[index]->parent();
681  break;
682 
683  case Var::_GalaxyTag_:
684  vv = m_object[index]->galaxyTag();
685  break;
686 
687  default:
688  ErrorCBL("no such a variable in the list!", "var", "Catalogue.cpp");
689  }
690 
691  return vv;
692 }
693 
694 
695 // ============================================================================
696 
697 int cbl::catalogue::Catalogue::var_int (int index, Var var_name) const
698 {
699  int vv;
700 
701  switch (var_name) {
702 
703  case Var::_ID_:
704  vv = m_object[index]->ID();
705  break;
706 
707  default:
708  ErrorCBL("no such a variable in the list!", "var", "Catalogue.cpp");
709  }
710 
711  return vv;
712 }
713 
714 // ============================================================================
715 
716 std::vector<int> cbl::catalogue::Catalogue::var_int (Var var_name) const
717 {
718  vector<int> vv(m_object.size(), 0.);
719 
720  for (size_t i=0; i<nObjects(); ++i) vv[i] = var(i, var_name);
721 
722  return vv;
723 }
724 
725 // ============================================================================
726 
727 std::vector<double> cbl::catalogue::Catalogue::var (Var var_name) const
728 {
729  vector<double> vv(m_object.size(), 0.);
730 
731  for (size_t i=0; i<nObjects(); ++i) vv[i] = var(i, var_name);
732 
733  return vv;
734 }
735 
736 
737 // ============================================================================
738 
739 
740 bool cbl::catalogue::Catalogue::isSetVar (int index, Var var_name) const
741 {
742  if (var_name==Var::_X_)
743  return m_object[index]->isSet_xx();
744 
745  else if (var_name==Var::_Y_)
746  return m_object[index]->isSet_yy();
747 
748  else if (var_name==Var::_Z_)
749  return m_object[index]->isSet_zz();
750 
751  else if (var_name==Var::_RA_)
752  return m_object[index]->isSet_ra();
753 
754  else if (var_name==Var::_Dec_)
755  return m_object[index]->isSet_dec();
756 
757  else if (var_name==Var::_TileRA_)
758  return m_object[index]->isSet_ra_tile();
759 
760  else if (var_name==Var::_TileDec_)
761  return m_object[index]->isSet_dec_tile();
762 
763  else if (var_name==Var::_SN_)
764  return m_object[index]->isSet_sn();
765 
766  else if (var_name==Var::_Redshift_)
767  return m_object[index]->isSet_redshift();
768 
769  else if (var_name==Var::_RedshiftMin_)
770  return m_object[index]->isSet_redshiftMin();
771 
772  else if (var_name==Var::_RedshiftMax_)
773  return m_object[index]->isSet_redshiftMax();
774 
775  else if (var_name==Var::_Shear1_)
776  return m_object[index]->isSet_shear1();
777 
778  else if (var_name==Var::_Shear2_)
779  return m_object[index]->isSet_shear2();
780 
781  else if (var_name==Var::_ODDS_)
782  return m_object[index]->isSet_odds();
783 
784  else if (var_name==Var::_LensingWeight_)
785  return m_object[index]->isSet_lensingWeight();
786 
787  else if (var_name==Var::_LensingCalib_)
788  return m_object[index]->isSet_lensingCalib();
789 
790  else if (var_name==Var::_Dc_)
791  return m_object[index]->isSet_dc();
792 
793  else if (var_name==Var::_Weight_)
794  return m_object[index]->isSet_weight();
795 
796  else if (var_name==Var::_Mass_)
797  return m_object[index]->isSet_mass();
798 
799  else if (var_name==Var::_Magnitude_)
800  return m_object[index]->isSet_magnitude();
801 
802  else if (var_name==Var::_MagnitudeU_)
803  return m_object[index]->isSet_magnitudeU();
804 
805  else if (var_name==Var::_MagnitudeG_)
806  return m_object[index]->isSet_magnitudeG();
807 
808  else if (var_name==Var::_MagnitudeR_)
809  return m_object[index]->isSet_magnitudeR();
810 
811  else if (var_name==Var::_MagnitudeI_)
812  return m_object[index]->isSet_magnitudeI();
813 
814  else if (var_name==Var::_SFR_)
815  return m_object[index]->isSet_SFR();
816 
817  else if (var_name==Var::_sSFR_)
818  return m_object[index]->isSet_sSFR();
819 
820  else if (var_name==Var::_MassProxy_)
821  return m_object[index]->isSet_mass_proxy();
822 
823  else if (var_name==Var::_MassProxyError_)
824  return m_object[index]->isSet_mass_proxy_error();
825 
826  else if (var_name==Var::_Mstar_)
827  return m_object[index]->isSet_mstar();
828 
829  else if (var_name==Var::_MassInfall_)
830  return m_object[index]->isSet_massinfall();
831 
832  else if (var_name==Var::_Vx_)
833  return m_object[index]->isSet_vx();
834 
835  else if (var_name==Var::_Vy_)
836  return m_object[index]->isSet_vy();
837 
838  else if (var_name==Var::_Vz_)
839  return m_object[index]->isSet_vz();
840 
841  else if (var_name==Var::_Region_)
842  return m_object[index]->isSet_region();
843 
844  else if (var_name==Var::_Generic_)
845  return m_object[index]->isSet_generic();
846 
847  else if (var_name==Var::_Radius_)
848  return m_object[index]->isSet_radius();
849 
850  else if (var_name==Var::_DensityContrast_)
851  return m_object[index]->isSet_densityContrast();
852 
853  else if (var_name==Var::_CentralDensity_)
854  return m_object[index]->isSet_centralDensity();
855 
856  else if (var_name==Var::_X_displacement_)
857  return m_object[index]->isSet_x_displacement();
858 
859  else if (var_name==Var::_Y_displacement_)
860  return m_object[index]->isSet_y_displacement();
861 
862  else if (var_name==Var::_Z_displacement_)
863  return m_object[index]->isSet_z_displacement();
864 
865  else if (var_name==Var::_MassEstimate_)
866  return m_object[index]->isSet_mass_estimate();
867 
868  else if (var_name==Var::_RadiusEstimate_)
869  return m_object[index]->isSet_radius_estimate();
870 
871  else if (var_name==Var::_VeldispEstimate_)
872  return m_object[index]->isSet_veldisp_estimate();
873 
874  else if (var_name==Var::_XCM_)
875  return m_object[index]->isSet_xcm();
876 
877  else if (var_name==Var::_YCM_)
878  return m_object[index]->isSet_ycm();
879 
880  else if (var_name==Var::_ZCM_)
881  return m_object[index]->isSet_zcm();
882 
883  else if (var_name==Var::_XSpin_)
884  return m_object[index]->isSet_spin_x();
885 
886  else if (var_name==Var::_YSpin_)
887  return m_object[index]->isSet_spin_y();
888 
889  else if (var_name==Var::_ZSpin_)
890  return m_object[index]->isSet_spin_z();
891 
892  else if (var_name==Var::_VelDisp_)
893  return m_object[index]->isSet_veldisp();
894 
895  else if (var_name==Var::_Vmax_)
896  return m_object[index]->isSet_vmax();
897 
898  else if (var_name==Var::_VmaxRad_)
899  return m_object[index]->isSet_vmax_rad();
900 
901  else if (var_name==Var::_TotMass_)
902  return m_object[index]->isSet_tot_mass();
903 
904  else if (var_name==Var::_ID_)
905  return m_object[index]->isSet_ID();
906 
907  if (var_name==Var::_GalaxyTag_)
908  return m_object[index]->isSet_galaxyTag();
909 
910  else
911  return ErrorCBL("no such a variable in the list!", "isSetVar", "Catalogue.cpp");
912 }
913 
914 
915 // ============================================================================
916 
917 
918 bool cbl::catalogue::Catalogue::isSetVar (Var var_name) const
919 {
920  bool ret = true;
921 
922  size_t i = 0;
923  while (ret==true && i<nObjects())
924  ret = isSetVar(i++, var_name);
925 
926  return ret;
927 }
928 
929 
930 // ============================================================================
931 
932 
933 void cbl::catalogue::Catalogue::set_region (const std::vector<long> region, const int nRegions)
934 {
935  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_region(region[i]);
936 
937  set_region_number( static_cast<size_t>( (nRegions<0) ? cbl::Max<long>(region)+1 : nRegions) );
938 }
939 
940 
941 // ============================================================================
942 
943 
944 void cbl::catalogue::Catalogue::set_region_number (const size_t nRegions)
945 {
946  for (size_t i=0; i<m_object.size(); i++)
947  if (m_object[i]->region()>=static_cast<int>(nRegions))
948  ErrorCBL("region index for object "+conv(i, par::fINT)+" is larger than input number of regions! "+conv(m_object[i]->region(), par::fINT)+" >= "+conv(nRegions, par::fINT), "set_region_number", "Catalogue.cpp");
949 
950  m_nRegions = nRegions;
951 }
952 
953 
954 // ============================================================================
955 
956 
957 void cbl::catalogue::Catalogue::set_ra_dec_tile (const std::vector<double> RA_tile, const std::vector<double> Dec_tile, const CoordinateUnits inputUnits)
958 {
959  for (size_t i=0; i<m_object.size(); i++) {
960  m_object[i]->set_ra_tile(RA_tile[i], inputUnits);
961  m_object[i]->set_dec_tile(Dec_tile[i], inputUnits);
962  }
963 }
964 
965 
966 // ============================================================================
967 
968 
969 void cbl::catalogue::Catalogue::set_field (const std::vector<std::string> field)
970 {
971  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_field(field[i]);
972 }
973 
974 
975 // ============================================================================
976 
977 
978 void cbl::catalogue::Catalogue::set_var (const int index, const Var var_name, const double value, const cosmology::Cosmology cosmology, const bool update_coordinates)
979 {
980  switch (var_name) {
981 
982  case Var::_X_:
983  m_object[index]->set_xx(value);
984  break;
985 
986  case Var::_Y_:
987  m_object[index]->set_yy(value);
988  break;
989 
990  case Var::_Z_:
991  m_object[index]->set_zz(value);
992  break;
993 
994  case Var::_RA_:
995  m_object[index]->set_ra(value);
996  break;
997 
998  case Var::_Dec_:
999  m_object[index]->set_dec(value);
1000  break;
1001 
1002  case Var::_SN_:
1003  m_object[index]->set_sn(value);
1004  break;
1005 
1006  case Var::_Redshift_:
1007  m_object[index]->set_redshift(value, cosmology, update_coordinates);
1008  break;
1009 
1010  case Var::_RedshiftMin_:
1011  m_object[index]->set_redshiftMin(value);
1012  break;
1013 
1014  case Var::_RedshiftMax_:
1015  m_object[index]->set_redshiftMax(value);
1016  break;
1017 
1018  case Var::_Shear1_:
1019  m_object[index]->set_shear1(value);
1020  break;
1021 
1022  case Var::_Shear2_:
1023  m_object[index]->set_shear2(value);
1024  break;
1025 
1026  case Var::_ODDS_:
1027  m_object[index]->set_odds(value);
1028  break;
1029 
1030  case Var::_LensingWeight_:
1031  m_object[index]->set_lensingWeight(value);
1032  break;
1033 
1034  case Var::_LensingCalib_:
1035  m_object[index]->set_lensingCalib(value);
1036  break;
1037 
1038  case Var::_Dc_:
1039  m_object[index]->set_dc(value);
1040  break;
1041 
1042  case Var::_Weight_:
1043  m_object[index]->set_weight(value);
1044  break;
1045 
1046  case Var::_Mass_:
1047  m_object[index]->set_mass(value);
1048  break;
1049 
1050  case Var::_Magnitude_:
1051  m_object[index]->set_magnitude(value);
1052  break;
1053 
1054  case Var::_MagnitudeU_:
1055  m_object[index]->set_magnitudeU(value);
1056  break;
1057 
1058  case Var::_MagnitudeG_:
1059  m_object[index]->set_magnitudeG(value);
1060  break;
1061 
1062  case Var::_MagnitudeR_:
1063  m_object[index]->set_magnitudeR(value);
1064  break;
1065 
1066  case Var::_MagnitudeI_:
1067  m_object[index]->set_magnitudeI(value);
1068  break;
1069 
1070  case Var::_SFR_:
1071  m_object[index]->set_SFR(value);
1072  break;
1073 
1074  case Var::_sSFR_:
1075  m_object[index]->set_sSFR(value);
1076  break;
1077 
1078  case Var::_MassProxy_:
1079  m_object[index]->set_mass_proxy(value);
1080  break;
1081 
1082  case Var::_MassProxyError_:
1083  m_object[index]->set_mass_proxy_error(value);
1084  break;
1085 
1086  case Var::_Mstar_:
1087  m_object[index]->set_mstar(value);
1088  break;
1089 
1090  case Var::_MassInfall_:
1091  m_object[index]->set_massinfall(value);
1092  break;
1093 
1094  case Var::_Vx_:
1095  m_object[index]->set_vx(value);
1096  break;
1097 
1098  case Var::_Vy_:
1099  m_object[index]->set_vy(value);
1100  break;
1101 
1102  case Var::_Vz_:
1103  m_object[index]->set_vz(value);
1104  break;
1105 
1106  case Var::_Region_:
1107  m_object[index]->set_region(value);
1108  break;
1109 
1110  case Var::_Generic_:
1111  m_object[index]->set_generic(value);
1112  break;
1113 
1114  case Var::_Radius_:
1115  m_object[index]->set_radius(value);
1116  break;
1117 
1118  case Var::_CentralDensity_:
1119  m_object[index]->set_centralDensity(value);
1120  break;
1121 
1122  case Var::_DensityContrast_:
1123  m_object[index]->set_densityContrast(value);
1124  break;
1125 
1126  case Var::_X_displacement_:
1127  m_object[index]->set_x_displacement(value);
1128  break;
1129 
1130  case Var::_Y_displacement_:
1131  m_object[index]->set_y_displacement(value);
1132  break;
1133 
1134  case Var::_Z_displacement_:
1135  m_object[index]->set_z_displacement(value);
1136  break;
1137 
1138  case Var::_MassEstimate_:
1139  m_object[index]->set_mass_estimate(value);
1140  break;
1141 
1142  case Var::_RadiusEstimate_:
1143  m_object[index]->set_radius_estimate(value);
1144  break;
1145 
1146  case Var::_VeldispEstimate_:
1147  m_object[index]->set_veldisp_estimate(value);
1148  break;
1149 
1150  case Var::_XCM_:
1151  m_object[index]->set_xcm(value);
1152  break;
1153 
1154  case Var::_YCM_:
1155  m_object[index]->set_ycm(value);
1156  break;
1157 
1158  case Var::_ZCM_:
1159  m_object[index]->set_zcm(value);
1160  break;
1161 
1162  case Var::_XSpin_:
1163  m_object[index]->set_spin_x(value);
1164  break;
1165 
1166  case Var::_YSpin_:
1167  m_object[index]->set_spin_y(value);
1168  break;
1169 
1170  case Var::_ZSpin_:
1171  m_object[index]->set_spin_z(value);
1172  break;
1173 
1174  case Var::_VelDisp_:
1175  m_object[index]->set_veldisp(value);
1176  break;
1177 
1178  case Var::_Vmax_:
1179  m_object[index]->set_vmax(value);
1180  break;
1181 
1182  case Var::_VmaxRad_:
1183  m_object[index]->set_vmax_rad(value);
1184  break;
1185 
1186  case Var::_TotMass_:
1187  m_object[index]->set_tot_mass(value);
1188  break;
1189 
1190  case Var::_GalaxyTag_:
1191  m_object[index]->set_galaxyTag(value);
1192  break;
1193 
1194  default:
1195  ErrorCBL("no such a variable in the list!", "set_var", "Catalogue.cpp");
1196  }
1197 
1198 }
1199 
1200 
1201 // ============================================================================
1202 
1203 
1204 void cbl::catalogue::Catalogue::set_var (const int index, const Var var_name, const int value, const cosmology::Cosmology cosmology)
1205 {
1206  switch (var_name) {
1207 
1208  case Var::_ID_:
1209  m_object[index]->set_ID(value);
1210  (void)cosmology;
1211  break;
1212 
1213  case Var::_IDHOST_:
1214  m_object[index]->set_IDHost(value);
1215  (void)cosmology;
1216  break;
1217 
1218  default:
1219  ErrorCBL("no such a variable in the list!", "set_var", "Catalogue.cpp");
1220  }
1221 
1222 }
1223 
1224 // ============================================================================
1225 
1226 
1227 // void cbl::catalogue::Catalogue::set_var (const int index, const Var var_name, const double value, const cosmology::Cosmology cosmology)
1228 // {
1229 // switch (var_name) {
1230 
1231 // case Var::_GalaxyTag_:
1232 // m_object[index]->set_galaxyTag(value);
1233 // (void)cosmology;
1234 // break;
1235 
1236 
1237 // default:
1238 // ErrorCBL("no such a variable in the list!", "set_var", "Catalogue.cpp");
1239 // }
1240 
1241 // }
1242 
1243 // ============================================================================
1244 
1245 
1246 void cbl::catalogue::Catalogue::set_var (const Var var_name, const std::vector<double> var, const cosmology::Cosmology cosmology, const bool update_coordinates)
1247 {
1248  if (m_object.size()!=var.size()) ErrorCBL("m_object.size()!=var.size()!", "set_var", "Catalogue.cpp");
1249 
1250  switch (var_name) {
1251 
1252  case Var::_X_:
1253  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_xx(var[i]);
1254  break;
1255 
1256  case Var::_Y_:
1257  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_yy(var[i]);
1258  break;
1259 
1260  case Var::_Z_:
1261  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_zz(var[i]);
1262  break;
1263 
1264  case Var::_RA_:
1265  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_ra(var[i]);
1266  break;
1267 
1268  case Var::_Dec_:
1269  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_dec(var[i]);
1270  break;
1271 
1272  case Var::_SN_:
1273  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_sn(var[i]);
1274  break;
1275 
1276  case Var::_Redshift_:
1277  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_redshift(var[i], cosmology, update_coordinates);
1278  break;
1279 
1280  case Var::_RedshiftMin_:
1281  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_redshiftMin(var[i]);
1282  break;
1283 
1284  case Var::_RedshiftMax_:
1285  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_redshiftMax(var[i]);
1286  break;
1287 
1288  case Var::_Shear1_:
1289  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_shear1(var[i]);
1290  break;
1291 
1292  case Var::_Shear2_:
1293  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_shear2(var[i]);
1294  break;
1295 
1296  case Var::_ODDS_:
1297  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_odds(var[i]);
1298  break;
1299 
1300  case Var::_LensingWeight_:
1301  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_lensingWeight(var[i]);
1302  break;
1303 
1304  case Var::_LensingCalib_:
1305  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_lensingCalib(var[i]);
1306  break;
1307 
1308  case Var::_Dc_:
1309  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_dc(var[i]);
1310  break;
1311 
1312  case Var::_Weight_:
1313  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_weight(var[i]);
1314  break;
1315 
1316  case Var::_Mass_:
1317  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_mass(var[i]);
1318  break;
1319 
1320  case Var::_Magnitude_:
1321  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_magnitude(var[i]);
1322  break;
1323 
1324  case Var::_MagnitudeU_:
1325  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_magnitudeU(var[i]);
1326  break;
1327 
1328  case Var::_MagnitudeG_:
1329  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_magnitudeG(var[i]);
1330  break;
1331 
1332  case Var::_MagnitudeR_:
1333  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_magnitudeR(var[i]);
1334  break;
1335 
1336  case Var::_MagnitudeI_:
1337  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_magnitudeI(var[i]);
1338  break;
1339 
1340  case Var::_SFR_:
1341  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_SFR(var[i]);
1342  break;
1343 
1344  case Var::_sSFR_:
1345  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_sSFR(var[i]);
1346  break;
1347 
1348  case Var::_MassProxy_:
1349  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_mass_proxy(var[i]);
1350  break;
1351 
1352  case Var::_MassProxyError_:
1353  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_mass_proxy_error(var[i]);
1354  break;
1355 
1356  case Var::_Mstar_:
1357  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_mstar(var[i]);
1358  break;
1359 
1360  case Var::_MassInfall_:
1361  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_massinfall(var[i]);
1362  break;
1363 
1364  case Var::_Vx_:
1365  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_vx(var[i]);
1366  break;
1367 
1368  case Var::_Vy_:
1369  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_vy(var[i]);
1370  break;
1371 
1372  case Var::_Vz_:
1373  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_vz(var[i]);
1374  break;
1375 
1376  case Var::_Region_:
1377  {
1378  vector<long> regList(nObjects());
1379  for (size_t i=0; i<nObjects(); ++i) regList[i] = static_cast<long>(var[i]);
1380  set_region(regList);
1381  break;
1382  }
1383 
1384  case Var::_Generic_:
1385  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_generic(var[i]);
1386  break;
1387 
1388  case Var::_Radius_:
1389  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_radius(var[i]);
1390  break;
1391 
1392  case Var::_CentralDensity_:
1393  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_centralDensity(var[i]);
1394  break;
1395 
1396  case Var::_DensityContrast_:
1397  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_densityContrast(var[i]);
1398  break;
1399 
1400  case Var::_X_displacement_:
1401  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_x_displacement(var[i]);
1402  break;
1403 
1404  case Var::_Y_displacement_:
1405  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_y_displacement(var[i]);
1406  break;
1407 
1408  case Var::_Z_displacement_:
1409  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_z_displacement(var[i]);
1410  break;
1411 
1412  case Var::_MassEstimate_:
1413  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_mass_estimate(var[i]);
1414  break;
1415 
1416  case Var::_RadiusEstimate_:
1417  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_radius_estimate(var[i]);
1418  break;
1419 
1420  case Var::_VeldispEstimate_:
1421  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_veldisp_estimate(var[i]);
1422  break;
1423 
1424  case Var::_XCM_:
1425  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_xcm(var[i]);
1426  break;
1427 
1428  case Var::_YCM_:
1429  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_ycm(var[i]);
1430  break;
1431 
1432  case Var::_ZCM_:
1433  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_zcm(var[i]);
1434  break;
1435 
1436  case Var::_XSpin_:
1437  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_spin_x(var[i]);
1438  break;
1439 
1440  case Var::_YSpin_:
1441  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_spin_y(var[i]);
1442  break;
1443 
1444  case Var::_ZSpin_:
1445  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_spin_z(var[i]);
1446  break;
1447 
1448  case Var::_VelDisp_:
1449  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_veldisp(var[i]);
1450  break;
1451 
1452  case Var::_Vmax_:
1453  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_vmax(var[i]);
1454  break;
1455 
1456  case Var::_VmaxRad_:
1457  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_vmax_rad(var[i]);
1458  break;
1459 
1460  case Var::_TotMass_:
1461  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_tot_mass(var[i]);
1462  break;
1463 
1464  case Var::_GalaxyTag_:
1465  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_galaxyTag(var[i]);
1466  break;
1467 
1468  default:
1469  ErrorCBL("no such a variable in the list!", "set_var", "Catalogue.cpp");
1470  }
1471 
1472 }
1473 
1474 // ============================================================================
1475 
1476 
1477 void cbl::catalogue::Catalogue::set_var (const Var var_name, const std::vector<int> var, const cosmology::Cosmology cosmology)
1478 {
1479  if (m_object.size()!=var.size()) ErrorCBL("m_object.size()!=var.size()!", "set_var", "Catalogue.cpp");
1480 
1481  switch (var_name) {
1482 
1483  case Var::_ID_:
1484  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_ID(var[i]);
1485  (void)cosmology;
1486  break;
1487 
1488  case Var::_IDHOST_:
1489  for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_IDHost(var[i]);
1490  (void)cosmology;
1491 
1492  break;
1493  default:
1494  ErrorCBL("no such a variable in the list!", "set_var", "Catalogue.cpp");
1495  }
1496 }
1497 
1498 //==============================================================================
1499 
1500 
1501 // void cbl::catalogue::Catalogue::set_var (const Var var_name, const std::vector<std::string> var, const cosmology::Cosmology cosmology)
1502 // {
1503 // if (m_object.size()!=var.size()) ErrorCBL("m_object.size()!=var.size()!", "set_var", "Catalogue.cpp");
1504 
1505 // switch (var_name) {
1506 
1507 // case Var::_GalaxyTag_:
1508 // for (size_t i=0; i<nObjects(); ++i) m_object[i]->set_galaxyTag(var[i]);
1509 // (void)cosmology;
1510 // break;
1511 
1512 // default:
1513 // ErrorCBL("no such a variable in the list!", "set_var", "Catalogue.cpp");
1514 // }
1515 
1516 // }
1517 
1518 // ============================================================================
1519 
1520 
1521 std::vector<double> cbl::catalogue::Catalogue::stats_var (const Var var_name) const
1522 {
1523  return { Average(var(var_name)), Sigma(var(var_name)), Quartile(var(var_name))[1], Quartile(var(var_name))[2]-Quartile(var(var_name))[0] };
1524 }
1525 
1526 
1527 // ============================================================================
1528 
1529 
1530 std::vector<std::vector<double>> cbl::catalogue::Catalogue::stats_var (const std::vector<Var> var_name) const
1531 {
1532  vector<vector<double>> stats;
1533 
1534  for (auto &&vv : var_name)
1535  stats.emplace_back(stats_var(vv));
1536 
1537  return stats;
1538 }
1539 
1540 
1541 // ============================================================================
1542 
1543 
1544 void cbl::catalogue::Catalogue::var_distr (const Var var_name, std::vector<double> &_var, std::vector<double> &dist, std::vector<double> &err, const int nbin, const bool linear, const std::string file_out, const double Volume, const bool norm, const double V1, const double V2, const std::string bin_type, const bool convolution, const double sigma) const
1545 {
1546  distribution(_var, dist, err, var(var_name), var(Var::_Weight_), nbin, linear, file_out, (norm) ? Volume*weightedN() : Volume, V1, V2, bin_type, convolution, sigma);
1547 }
1548 
1549 
1550 // ============================================================================
1551 
1552 
1553 void cbl::catalogue::Catalogue::computeComovingCoordinates (const cosmology::Cosmology &cosm, const CoordinateUnits inputUnits)
1554 {
1555  double red, xx, yy, zz;
1556 
1557 
1558  // ----- unit conversion -----
1559 
1560  vector<double> RA(nObjects()), DEC(nObjects());
1561 
1562  for (size_t i=0; i<nObjects(); ++i) {
1563  RA[i] = (inputUnits==CoordinateUnits::_radians_) ? ra(i) : radians(ra(i), inputUnits);
1564  DEC[i] = (inputUnits==CoordinateUnits::_radians_) ? dec(i) : radians(dec(i), inputUnits);
1565  }
1566 
1567 
1568  // ----- compute comoving coordinates -----
1569 
1570  for (size_t i=0; i<nObjects(); ++i) {
1571 
1572  red = redshift(i);
1573  m_object[i]->set_dc(cosm.D_C(red));
1574 
1575  cartesian_coord(RA[i], DEC[i], dc(i), xx, yy, zz);
1576 
1577  m_object[i]->set_xx(xx);
1578  m_object[i]->set_yy(yy);
1579  m_object[i]->set_zz(zz);
1580  }
1581 }
1582 
1583 
1584 // ============================================================================
1585 
1586 
1587 void cbl::catalogue::Catalogue::computePolarCoordinates (const CoordinateUnits outputUnits)
1588 {
1589  double ra, dec, dc;
1590 
1591 
1592  // ----- compute polar coordinates in radians -----
1593 
1594  for (size_t i=0; i<nObjects(); ++i) {
1595  polar_coord(xx(i), yy(i), zz(i), ra, dec, dc);
1596  m_object[i]->set_ra(ra);
1597  m_object[i]->set_dec(dec);
1598  m_object[i]->set_dc(dc);
1599  }
1600 
1601 
1602  // ----- unit conversion -----
1603 
1604  if (outputUnits!=CoordinateUnits::_radians_) {
1605 
1606  if (outputUnits==CoordinateUnits::_degrees_)
1607  for (size_t i=0; i<nObjects(); ++i) {
1608  ra = m_object[i]->ra();
1609  dec = m_object[i]->dec();
1610  m_object[i]->set_ra(degrees(ra));
1611  m_object[i]->set_dec(degrees(dec));
1612  }
1613 
1614  else if (outputUnits==CoordinateUnits::_arcseconds_)
1615  for (size_t i=0; i<nObjects(); ++i) {
1616  ra = m_object[i]->ra();
1617  dec = m_object[i]->dec();
1618  m_object[i]->set_ra(arcseconds(ra));
1619  m_object[i]->set_dec(arcseconds(dec));
1620  }
1621 
1622  else if (outputUnits==CoordinateUnits::_arcminutes_)
1623  for (size_t i=0; i<nObjects(); ++i) {
1624  ra = m_object[i]->ra();
1625  dec = m_object[i]->dec();
1626  m_object[i]->set_ra(arcminutes(ra));
1627  m_object[i]->set_dec(arcminutes(dec));
1628  }
1629 
1630  else ErrorCBL("outputUnits type not allowed!", "computePolarCoordinates", "Catalogue.cpp");
1631  }
1632 
1633 }
1634 
1635 // ============================================================================
1636 
1637 
1638 void cbl::catalogue::Catalogue::computePolarCoordinates (const cosmology::Cosmology &cosmology, const double z1, const double z2, const CoordinateUnits outputUnits)
1639 {
1640  double ra, dec, dc;
1641 
1642 
1643  // ----- compute polar coordinates in radians -----
1644 
1645  for (size_t i=0; i<nObjects(); ++i) {
1646  polar_coord(xx(i), yy(i), zz(i), ra, dec, dc);
1647  m_object[i]->set_ra(ra);
1648  m_object[i]->set_dec(dec);
1649  m_object[i]->set_dc(dc);
1650  m_object[i]->set_redshift(cosmology.Redshift(dc, z1, z2), cosmology);
1651  }
1652 
1653 
1654  // ----- unit conversion -----
1655 
1656  if (outputUnits!=CoordinateUnits::_radians_) {
1657 
1658  if (outputUnits==CoordinateUnits::_degrees_)
1659  for (size_t i=0; i<nObjects(); ++i) {
1660  m_object[i]->set_ra(degrees(ra));
1661  m_object[i]->set_dec(degrees(dec));
1662  }
1663 
1664  else if (outputUnits==CoordinateUnits::_arcseconds_)
1665  for (size_t i=0; i<nObjects(); ++i) {
1666  m_object[i]->set_ra(arcseconds(ra));
1667  m_object[i]->set_dec(arcseconds(dec));
1668  }
1669 
1670  else if (outputUnits==CoordinateUnits::_arcminutes_)
1671  for (size_t i=0; i<nObjects(); ++i) {
1672  m_object[i]->set_ra(arcminutes(ra));
1673  m_object[i]->set_dec(arcminutes(dec));
1674  }
1675 
1676  else ErrorCBL("outputUnits type not allowed!", "computePolarCoordinates", "Catalogue.cpp");
1677  }
1678 
1679 }
1680 
1681 // ============================================================================
1682 
1683 
1684 void cbl::catalogue::Catalogue::normalizeComovingCoordinates ()
1685 {
1686  for (size_t i=0; i<nObjects(); ++i) {
1687  m_object[i]->set_xx(xx(i)/dc(i));
1688  m_object[i]->set_yy(yy(i)/dc(i));
1689  m_object[i]->set_zz(zz(i)/dc(i));
1690  }
1691 }
1692 
1693 
1694 // ============================================================================
1695 
1696 
1697 void cbl::catalogue::Catalogue::restoreComovingCoordinates ()
1698 {
1699  for (size_t i=0; i<nObjects(); ++i) {
1700  m_object[i]->set_xx(xx(i)*dc(i));
1701  m_object[i]->set_yy(yy(i)*dc(i));
1702  m_object[i]->set_zz(zz(i)*dc(i));
1703  }
1704 }
1705 
1706 
1707 // ============================================================================
1708 
1709 
1710 void cbl::catalogue::Catalogue::Order (const std::vector<int> vv)
1711 {
1712  const int nObj = m_object.size();
1713 
1714  if (int(vv.size())!=nObj)
1715  ErrorCBL("vv.size()="+conv(vv.size(), par::fINT)+" and nObj="+conv(nObj, par::fINT)+" must be equal!", "Order", "Catalogue.cpp");
1716 
1717  vector<shared_ptr<Object>> obj(nObj);
1718 
1719  m_index.resize(nObj);
1720 
1721  for (size_t i=0; i<vv.size(); i++) {
1722  m_index[i] = vv[i];
1723  obj[i] = m_object[vv[i]];
1724  }
1725 
1726  m_object = obj;
1727 }
1728 
1729 
1730 // ============================================================================
1731 
1732 
1733 void cbl::catalogue::Catalogue::Order ()
1734 {
1735  if (m_index.size()>0) {
1736 
1737  const size_t nObj = m_object.size();
1738 
1739  if (m_index.size()!=nObj)
1740  ErrorCBL("m_index.size()="+conv(m_index.size(), par::fINT)+" and nObj="+conv(nObj, par::fINT)+" must be equal!", "Order", "Catalogue.cpp");
1741 
1742  vector<shared_ptr<Object>> obj(nObj);
1743 
1744  obj = m_object;
1745 
1746  for (size_t i=0; i<nObj; i++)
1747  m_object[i] = obj[m_index[i]];
1748 
1749  }
1750 }
1751 
1752 
1753 double cbl::catalogue::Catalogue::weightedN () const
1754 {
1755  double nn = 0.;
1756  for (size_t i=0; i<m_object.size(); i++)
1757  nn += m_object[i]->weight();
1758  return nn;
1759 }
1760 
1761 
1762 // ============================================================================
1763 
1764 
1765 void cbl::catalogue::Catalogue::write_comoving_coordinates (const std::string outputFile) const
1766 {
1767  if (m_object.size()==0) ErrorCBL("m_object.size()=0!", "write_comoving_coordinates", "Catalogue.cpp");
1768 
1769  coutCBL << "I'm writing the file: " << outputFile << "..." << endl;
1770 
1771  ofstream fout(outputFile.c_str()); checkIO(fout, outputFile);
1772 
1773  for (size_t i=0; i<nObjects(); ++i)
1774  fout << xx(i) << " " << yy(i) << " " << zz(i) << endl;
1775 
1776  fout.clear(); fout.close();
1777 }
1778 
1779 
1780 // ============================================================================
1781 
1782 
1783 void cbl::catalogue::Catalogue::write_obs_coordinates (const std::string outputFile) const
1784 {
1785  if (m_object.size()==0) ErrorCBL("m_object.size()=0!", "write_obs_coordinates", "Catalogue.cpp");
1786 
1787  coutCBL << "I'm writing the file: " << outputFile << "..." << endl;
1788 
1789  ofstream fout(outputFile.c_str()); checkIO(fout, outputFile);
1790 
1791  if (!isSet(ra(0)) || !isSet(dec(0)) || !isSet(redshift(0)))
1792  ErrorCBL("polar coordinates are not set!", "write_obs_coordinates", "Catalogue.cpp");
1793 
1794  for (size_t i=0; i<nObjects(); ++i)
1795  fout << ra(i) << " " << dec(i) << " " << redshift(i) << endl;
1796 
1797  fout.clear(); fout.close();
1798 }
1799 
1800 
1801 // ============================================================================
1802 
1803 
1804 void cbl::catalogue::Catalogue::write_data (const std::string outputFile, const std::vector<Var> var_name, const std::string sep, const std::string header) const
1805 {
1806  coutCBL << "Writing the file: " << outputFile << "..." << endl;
1807 
1808  ofstream fout(outputFile.c_str()); checkIO(fout, outputFile);
1809 
1810  if (header != "") {
1811  fout << header << endl;
1812  }
1813 
1814  if (var_name.size()==0) {
1815 
1816  if (!isSet(ra(0)) || !isSet(dec(0)) || !isSet(redshift(0)) || !isSet(dc(0)))
1817  ErrorCBL("Polar coordinates are not set!", "write_data", "Catalogue.cpp");
1818 
1819  if (!isSet(region(0)))
1820  for (size_t i=0; i<nObjects(); ++i)
1821  fout << xx(i) << sep << yy(i) << sep << zz(i) << sep << ra(i) << sep << dec(i) << sep << redshift(i) << sep << dc(i) << endl;
1822  else
1823  for (size_t i=0; i<nObjects(); ++i)
1824  fout << xx(i) << sep << yy(i) << sep << zz(i) << sep << ra(i) << sep << dec(i) << sep << redshift(i) << sep << dc(i) << sep << region(i) << endl;
1825 
1826  }
1827 
1828  else {
1829  vector<vector<double>> data;
1830  for (size_t j=0; j<var_name.size(); j++)
1831  data.push_back(var(var_name[j]));
1832 
1833  for (size_t i=0; i<nObjects(); ++i) {
1834  for (size_t j=0; j<data.size(); j++)
1835  fout << data[j][i] << sep;
1836  fout << endl;
1837  }
1838 
1839  }
1840 
1841  fout.clear(); fout.close();
1842 }
1843 
1844 
1845 // ============================================================================
1846 
1847 
1848 Catalogue cbl::catalogue::Catalogue::sub_catalogue (const mask_function mask, const bool excl) const
1849 {
1850  vector<shared_ptr<Object>> objects;
1851  vector<bool> w(m_object.size());
1852 
1853  bool fact = (excl) ? false : true;
1854 
1855  for (size_t i=0; i<m_object.size(); i++)
1856  w[i] = mask(m_object[i])&&fact;
1857 
1858  for (size_t i=0; i<m_object.size(); i++)
1859  if (w[i])
1860  objects.push_back(m_object[i]);
1861 
1862  return Catalogue{objects};
1863 }
1864 
1865 
1866 // ============================================================================
1867 
1868 
1869 Catalogue cbl::catalogue::Catalogue::sub_catalogue (const cbl::catalogue::MaskObject &mask, const bool excl) const
1870 {
1871  vector<shared_ptr<Object>> objects;
1872  vector<bool> w(m_object.size());
1873 
1874  bool fact = (excl) ? false : true;
1875 
1876  for (size_t i=0; i<m_object.size(); i++)
1877  w[i] = mask(m_object[i])&&fact;
1878 
1879  for (size_t i=0; i<m_object.size(); i++)
1880  if (w[i])
1881  objects.push_back(m_object[i]);
1882 
1883  return Catalogue{objects};
1884 }
1885 
1886 
1887 // ============================================================================
1888 
1889 
1890 Catalogue cbl::catalogue::Catalogue::sub_catalogue (const Var var_name, const double down, const double up, const bool excl) const
1891 {
1892  vector<shared_ptr<Object>> objects;
1893  vector<double> vvar = var(var_name);
1894  vector<int> w(vvar.size());
1895 
1896  for (size_t i=0; i<m_object.size(); i++) {
1897  w[i] = (excl) ? 1 : 0;
1898  if (vvar[i] >= down && vvar[i] < up)
1899  w[i] = (excl) ? 0 : 1;
1900  }
1901 
1902  for (size_t i=0; i<m_object.size(); i++)
1903  if (w[i]==1)
1904  objects.push_back(m_object[i]);
1905 
1906  return Catalogue{objects};
1907 }
1908 
1909 
1910 // ============================================================================
1911 
1912 
1913 Catalogue cbl::catalogue::Catalogue::sub_catalogue (catalogue::Catalogue &data, const double tile_width_RA, const double tile_width_Dec, const bool write_tiles, const std::string dir_tiles, const std::string file_tiles)
1914 {
1915  const double RA_hw = 0.5 * tile_width_RA * (par::pi/180.); // half tile width in radians along R.A.
1916  const double Dec_hw = 0.5 * tile_width_Dec * (par::pi/180.); // along Dec
1917 
1918 
1919  // Set the vector of tile numbers
1920  std::vector<long int> dummy_tiles = data.region();
1921  std::sort(dummy_tiles.begin(), dummy_tiles.end());
1922  std::vector<long int> unique_tile_numbers = cbl::different_elements (dummy_tiles);
1923  const int n_tiles = (int)(unique_tile_numbers.size());
1924 
1925 
1926  // Define the vectors of min/max R.A./Dec,
1927  // and check if the internal variables are properly set
1928  std::vector<double> RA_min (n_tiles);
1929  std::vector<double> RA_max (n_tiles);
1930  std::vector<double> Dec_min (n_tiles);
1931  std::vector<double> Dec_max (n_tiles);
1932 
1933  std::vector<bool> isSet_region (n_tiles);
1934 
1935  for (size_t i=0; i<data.nObjects(); i++) {
1936 
1937  if (data.isSetVar(i, catalogue::Var::_TileRA_) == false)
1938  ErrorCBL("The tile central R.A. for the object "+cbl::conv(i,cbl::par::fINT)+" in the data catalogue is not set.", "sub_catalogue", "Catalogue.cpp");
1939  if (data.isSetVar(i, catalogue::Var::_TileDec_) == false)
1940  ErrorCBL("The tile central Dec for the object "+cbl::conv(i,cbl::par::fINT)+" in the data catalogue is not set.", "sub_catalogue", "Catalogue.cpp");
1941  if (data.isSetVar(i, catalogue::Var::_Region_) == false)
1942  ErrorCBL("The tile number for the object "+cbl::conv(i,cbl::par::fINT)+" in the data catalogue is not set.", "sub_catalogue", "Catalogue.cpp");
1943  if (data.var(i, catalogue::Var::_Region_) < 0)
1944  ErrorCBL("The tile number for the object "+cbl::conv(i,cbl::par::fINT)+" in the data catalogue is <0. The tile numbers must be all the integers between 0 and N, where N is the highest tile number.", "sub_catalogue", "Catalogue.cpp");
1945 
1946  if (data.region(i) < n_tiles)
1947  isSet_region[data.region(i)] = true;
1948  else
1949  ErrorCBL("The tile number "+cbl::conv(i,cbl::par::fINT)+" in the data catalogue is greater than the number of tiles. The tile numbers must be all the integers between 0 and N, where N is the highest tile number.", "sub_catalogue", "Catalogue.cpp");
1950 
1951  // Set min/max R.A./Dec for each tile
1952  double candidate_RA_min = data.ra_tile(i) - RA_hw/cos(std::abs(data.dec_tile(i)));
1953  if (candidate_RA_min < 0)
1954  candidate_RA_min = 2*cbl::par::pi + candidate_RA_min;
1955 
1956  double candidate_RA_max = data.ra_tile(i) + RA_hw/cos(std::abs(data.dec_tile(i)));
1957  if (candidate_RA_max > 2*cbl::par::pi)
1958  candidate_RA_max = candidate_RA_max - 2*cbl::par::pi;
1959 
1960  double candidate_Dec_min = data.dec_tile(i) - Dec_hw;
1961  double candidate_Dec_max = data.dec_tile(i) + Dec_hw;
1962 
1963  RA_min[data.region(i)] = candidate_RA_min;
1964  RA_max[data.region(i)] = candidate_RA_max;
1965  Dec_min[data.region(i)] = candidate_Dec_min;
1966  Dec_max[data.region(i)] = candidate_Dec_max;
1967 
1968  }
1969 
1970  for (int i=0; i<n_tiles; i++)
1971  if (isSet_region[i] == false)
1972  ErrorCBL("The tile number "+cbl::conv(i,cbl::par::fINT)+" in the original catalogue (data) is missing. The tile numbers must be all the integers between 0 and N, where N is the highest tile number.", "sub_catalogue", "Catalogue.cpp");
1973 
1974  std::vector<bool>().swap(isSet_region);
1975 
1976 
1977  // Write the tiles file
1978  if (write_tiles) {
1979  std::string mkdir = "mkdir -p "+dir_tiles; if (system(mkdir.c_str())) {}
1980  std::ofstream myfile; myfile.open(dir_tiles+file_tiles);
1981  myfile << "# RA_min RA_max Dec_min Dec_max [all in degrees]" <<std::endl;
1982  for (int i=0; i<n_tiles; i++)
1983  myfile << RA_min[i]*(180./cbl::par::pi) << std::setw(20) << RA_max[i]*(180./cbl::par::pi) << std::setw(20) << Dec_min[i]*(180./cbl::par::pi) << std::setw(20) << Dec_max[i]*(180./cbl::par::pi) << std::endl;
1984  myfile.close();
1985  coutCBL<<"I wrote the file "+dir_tiles+file_tiles<<std::endl;
1986  }
1987 
1988 
1989  // Collect the objects falling in the tiles
1990  std::vector<shared_ptr<Object>> objects;
1991 
1992  for (size_t i=0; i<this->nObjects(); i++) {
1993 
1994  for (int j=0; j<n_tiles; j++) {
1995 
1996  if (Dec_min[j] < this->dec(i) && Dec_max[j] > this->dec(i)) {
1997 
1998  if (RA_max[j] > RA_min[j]) { // To manage the cases in which a tile lies at the zero of R.A.
1999 
2000  if (RA_min[j] < this->ra(i) && RA_max[j] > this->ra(i)) {
2001 
2002  objects.push_back(m_object[i]);
2003  break;
2004  }
2005 
2006  } else {
2007 
2008  if (RA_min[j] < this->ra(i) || RA_max[j] > this->ra(i)) { // To manage the cases in which a tile lies at the zero of R.A.
2009 
2010  objects.push_back(m_object[i]);
2011  break;
2012  }
2013 
2014  }
2015  }
2016  }
2017  }
2018 
2019  return Catalogue{objects};
2020 
2021 }
2022 
2023 
2024 // ============================================================================
2025 
2026 
2027 Catalogue cbl::catalogue::Catalogue::mangle_cut (const std::string mangle_mask, const bool excl) const
2028 {
2029 
2030  vector<shared_ptr<Object>> objects;
2031 
2032  cbl::Path path;
2033  string mangle_dir = path.DirCosmo()+"/External/mangle/";
2034 
2035  string mangle_working_dir = mangle_dir+"output/";
2036  string mkdir = "mkdir -p "+mangle_working_dir;
2037  if (system(mkdir.c_str())) {}
2038 
2039  string input = mangle_working_dir+"temporary.dat";
2040  string output = mangle_working_dir+"temporary_output.dat";
2041 
2042  write_obs_coordinates(input);
2043 
2044  string polyid = mangle_dir+"/bin/polyid -ur "+mangle_mask+" "+input+" "+output;
2045  if (system(polyid.c_str())) {}
2046 
2047  ifstream fin(output);
2048  string line;
2049 
2050  int nn=0;
2051  getline(fin, line);
2052 
2053  while(getline(fin, line))
2054  {
2055  stringstream ss(line);
2056  vector<double> num; double NN = par::defaultDouble;
2057  while (ss>>NN) num.push_back(NN);
2058 
2059  if( !excl && num.size()>2)
2060  objects.push_back(m_object[nn]);
2061  if( excl && num.size()<3)
2062  objects.push_back(m_object[nn]);
2063  nn++;
2064  }
2065  fin.clear(); fin.close();
2066 
2067  string rm = "rm "+input+" "+output;
2068  if (system(rm.c_str())) {}
2069 
2070  return Catalogue{objects};
2071 }
2072 
2073 // ============================================================================
2074 
2075 
2076 Catalogue cbl::catalogue::Catalogue::diluted_catalogue (const double nSub, const int seed) const
2077 {
2078  if (nSub<=0 || nSub>1 || !isfinite(nSub)) ErrorCBL("nSub must be in the range (0,1] !", "diluted_catalogue", "Catalogue.cpp");
2079 
2080  // copy the catalogue into a new one
2081  auto diluted_catalogue = *this;
2082 
2083  // set the index vector that will be used to remove the objects
2084  vector<bool> index(nObjects(), false);
2085 
2086  // nObjects()*(1-nSub) will be removed
2087  for (size_t i=0; i<nObjects()*(1-nSub); ++i) index[i] = true;
2088 
2089  // shuffle the indexes of the objects that will be removed
2090  default_random_engine engine(seed);
2091  std::shuffle(begin(index), end(index), engine);
2092 
2093  // dilute the new catalogue
2094  diluted_catalogue.remove_objects(index);
2095 
2096  return diluted_catalogue;
2097 }
2098 
2099 
2100 // ============================================================================
2101 
2102 
2103 double cbl::catalogue::Catalogue::distance (const int i, shared_ptr<Object> obj) const
2104 {
2105  return sqrt((m_object[i]->xx()-obj->xx())*(m_object[i]->xx()-obj->xx())+
2106  (m_object[i]->yy()-obj->yy())*(m_object[i]->yy()-obj->yy())+
2107  (m_object[i]->zz()-obj->zz())*(m_object[i]->zz()-obj->zz()));
2108 }
2109 
2110 
2111 // ============================================================================
2112 
2113 
2114 double cbl::catalogue::Catalogue::angsep_xyz (const int i, shared_ptr<Object> obj) const
2115 {
2116  return 2.*asin(0.5*sqrt((m_object[i]->xx()-obj->xx())*(m_object[i]->xx()-obj->xx())+
2117  (m_object[i]->yy()-obj->yy())*(m_object[i]->yy()-obj->yy())+
2118  (m_object[i]->zz()-obj->zz())*(m_object[i]->zz()-obj->zz())));
2119 }
2120 
2121 
2122 // ============================================================================
2123 
2124 
2125 shared_ptr<Catalogue> cbl::catalogue::Catalogue::smooth (const double gridsize, const cosmology::Cosmology cosmology, const std::vector<Var> vars, const int SUB)
2126 {
2127  (void)vars;
2128 
2129  shared_ptr<Catalogue> cat {new Catalogue(*this)};
2130 
2131  if (gridsize<1.e-30) return cat;
2132 
2133  double rMAX = 0.;
2134 
2135  vector<shared_ptr<Object>> sample;
2136 
2137 
2138  // ----------------------------------------------------------------------------
2139  // ----- subdivide the catalogue in SUB^3 sub-catalogues, to avoid ------------
2140  // ----- memory problems possibly caused by too large chain-mesh vectors ------
2141  // ----------------------------------------------------------------------------
2142 
2143  coutCBL <<"Please wait, I'm subdividing the catalogue in "<<pow(SUB, 3)<<" sub-catalogues..."<<endl;
2144 
2145  int nx = SUB, ny = SUB, nz = SUB;
2146 
2147  double Cell_X = (Max(Var::_X_)-Min(Var::_X_))/nx;
2148  double Cell_Y = (Max(Var::_Y_)-Min(Var::_Y_))/ny;
2149  double Cell_Z = (Max(Var::_Z_)-Min(Var::_Z_))/nz;
2150  double MinX = Min(Var::_X_);
2151  double MinY = Min(Var::_Y_);
2152  double MinZ = Min(Var::_Z_);
2153 
2154  vector<long> regions(cat->nObjects());
2155 
2156  for (size_t i=0; i<cat->nObjects(); i++) {
2157  int i1 = min(int((cat->xx(i)-MinX)/Cell_X), nx-1);
2158  int j1 = min(int((cat->yy(i)-MinY)/Cell_Y), ny-1);
2159  int z1 = min(int((cat->zz(i)-MinZ)/Cell_Z), nz-1);
2160  long index = z1+nz*(j1+ny*i1);
2161  regions[i] = index;
2162  }
2163 
2164  cat->set_region(regions);
2165 
2166  size_t nRegions = cat->nRegions();
2167 
2168  vector<Catalogue> subSamples(nRegions);
2169 
2170  for (size_t i=0; i<nRegions; ++i) {
2171  double start = (double)cat->region_list()[i];
2172  double stop = start+1;
2173  subSamples[i] = sub_catalogue(Var::_Region_, start, stop);
2174  }
2175 
2176 
2177  // -----------------------------------------
2178  // ----- smooth all the sub-catalogues -----
2179  // -----------------------------------------
2180 
2181  for (size_t rr=0; rr<nRegions; ++rr) {
2182 
2183  vector<double> _xx = subSamples[rr].var(Var::_X_), _yy = subSamples[rr].var(Var::_Y_), _zz = subSamples[rr].var(Var::_Z_);
2184 
2185  ChainMesh3D ll(gridsize, _xx, _yy, _zz, rMAX, (long)-1.e5, (long)1.e5);
2186 
2187 
2188  for (int i=0; i<ll.nCell(); i++) {
2189  vector<long> list = ll.get_list(i);
2190 
2191  int nObj = list.size();
2192  if (nObj>0) {
2193 
2194  double XX = 0., YY = 0., ZZ = 0., RA = 0., DEC = 0., REDSHIFT = 0., WEIGHT = 0.;
2195 
2196  for (size_t j=0; j<list.size(); j++) {
2197  XX += subSamples[rr].xx(list[j]);
2198  YY += subSamples[rr].yy(list[j]);
2199  ZZ += subSamples[rr].zz(list[j]);
2200  RA += subSamples[rr].ra(list[j]);
2201  DEC += subSamples[rr].dec(list[j]);
2202  REDSHIFT += subSamples[rr].redshift(list[j]);
2203  WEIGHT += subSamples[rr].weight(list[j]);
2204  }
2205 
2206  shared_ptr<Object> obj{new Object()};
2207  XX /= nObj; obj->set_xx(XX);
2208  YY /= nObj; obj->set_yy(YY);
2209  ZZ /= nObj; obj->set_zz(ZZ);
2210  RA /= nObj; obj->set_ra(RA);
2211  DEC /= nObj; obj->set_dec(DEC);
2212  REDSHIFT /= nObj; obj->set_redshift(REDSHIFT, cosmology);
2213  obj->set_weight(WEIGHT);
2214  sample.push_back(obj);
2215 
2216  }
2217  }
2218 
2219  }
2220 
2221  shared_ptr<Catalogue> cat_new(new Catalogue{sample});
2222  return cat_new;
2223 }
2224 
2225 
2226 // ============================================================================
2227 
2228 
2229 int cbl::catalogue::Catalogue::nObjects_condition (const Var var_name, const double down, const double up, const bool excl)
2230 {
2231  int nObjw = 0;
2232  vector<double> vvar = var(var_name);
2233 
2234  for (size_t i=0; i<m_object.size(); i++)
2235 
2236  if (vvar[i] >= down && vvar[i] < up)
2237  nObjw ++;
2238 
2239  nObjw = (excl) ? weightedN()-nObjw : nObjw;
2240 
2241  return nObjw;
2242 }
2243 
2244 
2245 // ============================================================================
2246 
2247 
2248 double cbl::catalogue::Catalogue::weightedN_condition (const Var var_name, const double down, const double up, const bool excl)
2249 {
2250  double nObjw = 0;
2251  vector<double> vvar = var(var_name);
2252 
2253  for (size_t i=0; i<m_object.size(); i++)
2254  if (vvar[i] >= down && vvar[i] < up)
2255  nObjw += weight(i);
2256 
2257  nObjw = (excl) ? weightedN()-nObjw : nObjw;
2258 
2259  return nObjw;
2260 }
2261 
2262 
2263 // ============================================================================
2264 
2265 
2266 data::ScalarField3D cbl::catalogue::Catalogue::counts_in_cell (const double cell_size, const int interpolation_type, const bool useMass, const double minX, const double maxX, const double minY, const double maxY, const double minZ, const double maxZ) const
2267 {
2268  double _minX = (minX>par::defaultDouble) ? minX : Min(Var::_X_)-3*cell_size;
2269  double _maxX = (maxX>par::defaultDouble) ? maxX : Max(Var::_X_)+3*cell_size;
2270  double _minY = (minY>par::defaultDouble) ? minY : Min(Var::_Y_)-3*cell_size;
2271  double _maxY = (maxY>par::defaultDouble) ? maxY : Max(Var::_Y_)+3*cell_size;
2272  double _minZ = (minZ>par::defaultDouble) ? minZ : Min(Var::_Z_)-3*cell_size;
2273  double _maxZ = (maxZ>par::defaultDouble) ? maxZ : Max(Var::_Z_)+3*cell_size;
2274 
2275  data::ScalarField3D density(cell_size, _minX, _maxX, _minY, _maxY, _minZ, _maxZ);
2276 
2277  double deltaX = density.deltaX();
2278  double deltaY = density.deltaY();
2279  double deltaZ = density.deltaZ();
2280  int nx = density.nx();
2281  int ny = density.ny();
2282  int nz = density.nz();
2283 
2284  for (size_t i=0; i<nObjects(); ++i) {
2285  int i1 = min(int((xx(i)-density.MinX())/deltaX),nx-1);
2286  int j1 = min(int((yy(i)-density.MinY())/deltaY),ny-1);
2287  int k1 = min(int((zz(i)-density.MinZ())/deltaZ),nz-1);
2288 
2289  double w = (useMass) ? mass(i)*weight(i) : weight(i);
2290 
2291  if (interpolation_type==0) {
2292  density.set_ScalarField(w, i1, j1, k1, 1);
2293  }
2294  else if (interpolation_type==1) {
2295 
2296  double dx = (xx(i)-density.XX(i1))/deltaX;
2297  double dy = (yy(i)-density.YY(j1))/deltaY;
2298  double dz = (zz(i)-density.ZZ(k1))/deltaZ;
2299 
2300  int i2 = (dx<0) ? i1-1 : i1+1;
2301  int j2 = (dy<0) ? j1-1 : j1+1;
2302  int k2 = (dz<0) ? k1-1 : k1+1;
2303 
2304  dx = fabs(dx);
2305  dy = fabs(dy);
2306  dz = fabs(dz);
2307 
2308  double tx = 1-dx;
2309  double ty = 1-dy;
2310  double tz = 1-dz;
2311 
2312  double ww = 0.;
2313  ww += w*tx*ty*tz;
2314  ww += w*dx*ty*tz;
2315  ww += w*tx*dy*tz;
2316  ww += w*tx*ty*dz;
2317  ww += w*dx*dy*tz;
2318  ww += w*dx*ty*dz;
2319  ww += w*tx*dy*dz;
2320  ww += w*dx*dy*dz;
2321 
2322  density.set_ScalarField(w*tx*ty*tz, i1, j1, k1, 1);
2323  density.set_ScalarField(w*dx*ty*tz, i2, j1, k1, 1);
2324  density.set_ScalarField(w*tx*dy*tz, i1, j2, k1, 1);
2325  density.set_ScalarField(w*tx*ty*dz, i1, j1, k2, 1);
2326  density.set_ScalarField(w*dx*dy*tz, i2, j2, k1, 1);
2327  density.set_ScalarField(w*dx*ty*dz, i2, j1, k2, 1);
2328  density.set_ScalarField(w*tx*dy*dz, i1, j2, k2, 1);
2329  density.set_ScalarField(w*dx*dy*dz, i2, j2, k2, 1);
2330 
2331  }
2332  }
2333 
2334  return density;
2335 }
2336 
2337 
2338 // ============================================================================
2339 
2340 
2341 data::ScalarField3D cbl::catalogue::Catalogue::density_field (const double cell_size, const Catalogue mask_catalogue, const int interpolation_type, const double kernel_radius, const bool useMass) const
2342 {
2343 
2344  data::ScalarField3D data_cic = counts_in_cell(cell_size, interpolation_type, useMass);
2345 
2346  data::ScalarField3D mask_cic = mask_catalogue.counts_in_cell(cell_size, 0 /*interpolation_type*/, useMass, data_cic.MinX(), data_cic.MaxX(), data_cic.MinY(), data_cic.MaxY(), data_cic.MinZ(), data_cic.MaxZ());
2347 
2348 
2349  data::ScalarField3D density(cell_size, data_cic.MinX(), data_cic.MaxX(), data_cic.MinY(), data_cic.MaxY(), data_cic.MinZ(), data_cic.MaxZ());
2350 
2351  double data_tot=0, random_tot=0;
2352  int nrandom = 0;
2353 
2354  for (int i=0; i<density.nx(); i++)
2355  for (int j=0; j<density.ny(); j++)
2356  for (int k=0; k<density.nz(); k++) {
2357  if (mask_cic.ScalarField(i, j, k)>0) {
2358  random_tot += mask_cic.ScalarField(i, j, k);
2359  nrandom ++;
2360  }
2361  }
2362 
2363  double mean_random = random_tot/nrandom;
2364  coutCBL << "Mean random objects " << mean_random << " in " << nrandom << " cells " << endl;
2365 
2366  random_tot=0;
2367  int masked_cells=0;
2368 
2369  for (int i=0; i<density.nx(); i++)
2370  for (int j=0; j<density.ny(); j++)
2371  for (int k=0; k<density.nz(); k++) {
2372  if (mask_cic.ScalarField(i, j, k)>0.) {
2373  data_tot += data_cic.ScalarField(i, j, k);
2374  random_tot += mask_cic.ScalarField(i, j, k);
2375  }
2376  else if (mask_cic.ScalarField(i, j, k)>0 && mask_cic.ScalarField(i, j, k)<0.1*mean_random)
2377  {
2378  masked_cells ++;
2379  data_cic.set_ScalarField(0, i, j, k, 0);
2380  mask_cic.set_ScalarField(0, i, j, k, 0);
2381  }
2382  }
2383 
2384  coutCBL << "Masked " << masked_cells << "/" << nrandom << " for bad random coverage " << endl;
2385 
2386  double norm = int(random_tot)/data_tot;
2387  for (int i=0; i<density.nx(); i++)
2388  for (int j=0; j<density.ny(); j++)
2389  for (int k=0; k<density.nz(); k++) {
2390  double val=0;
2391  val = (mask_cic.ScalarField(i, j, k)>0) ? data_cic.ScalarField(i,j,k)/mask_cic.ScalarField(i,j,k)*norm-1 : 0;
2392  density.set_ScalarField(val, i, j, k);
2393  }
2394 
2395  if (kernel_radius>0)
2396  density.GaussianConvolutionField(kernel_radius);
2397 
2398  return density;
2399 }
2400 
2401 
2402 // ============================================================================
2403 
2404 
2405 cbl::catalogue::Catalogue::Catalogue (const Catalogue input_catalogue, const Catalogue target_catalogue, const Var var_name, const int nbin, const int seed)
2406 {
2407  vector<double> fvar_input(nbin, 0);
2408  vector<double> fvar_target(nbin, 0);
2409 
2410  vector<double> input_var = input_catalogue.var(var_name);
2411  vector<double> target_var = target_catalogue.var(var_name);
2412 
2413  double Vmin = target_catalogue.Min(var_name);
2414  double Vmax = target_catalogue.Max(var_name);
2415 
2416  double binSize_inv = pow((Vmax-Vmin)/nbin, -1);
2417 
2418  for (size_t i=0; i<input_var.size(); i++)
2419  if (input_var[i] < Vmax && Vmin < input_var[i]) {
2420  int occ = max(0, min(int((input_var[i]-Vmin)*binSize_inv), nbin));
2421  fvar_input[occ] ++;
2422  }
2423 
2424  for (size_t i=0; i<target_var.size(); i++)
2425  if (target_var[i] < Vmax && Vmin < target_var[i]) {
2426  int occ = max(0, min(int((target_var[i]-Vmin)*binSize_inv), nbin));
2427  fvar_target[occ] ++;
2428  }
2429 
2430  random::UniformRandomNumbers ran(0., 1., seed);
2431 
2432  for (size_t i=0; i<input_var.size(); i++)
2433  if (input_var[i] < Vmax && Vmin < input_var[i]) {
2434  int occ = max(0, min(int((input_var[i]-Vmin)*binSize_inv), nbin));
2435  if (ran() < fvar_target[occ]/fvar_input[occ])
2436  m_object.push_back(shared_ptr<Object>(input_catalogue.catalogue_object(i)));
2437  }
2438 }
2439 
2440 
2441 // ============================================================================
2442 
2443 
2444 cbl::catalogue::Catalogue::Catalogue (const Catalogue input_catalogue, const Catalogue target_catalogue, const Var var_name1, const int nbin1, const Var var_name2, const int nbin2, const int seed)
2445 {
2446  vector< vector<double> > fvars_input(nbin1, vector<double>(nbin2, 0));
2447  vector< vector<double> > fvars_target(nbin1, vector<double>(nbin2, 0));
2448 
2449  vector<double> input_var1 = input_catalogue.var(var_name1);
2450  vector<double> target_var1 = target_catalogue.var(var_name1);
2451  vector<double> input_var2 = input_catalogue.var(var_name2);
2452  vector<double> target_var2 = target_catalogue.var(var_name2);
2453 
2454  double V1min = target_catalogue.Min(var_name1);
2455  double V1max = target_catalogue.Max(var_name1);
2456 
2457  double V2min = target_catalogue.Min(var_name2);
2458  double V2max = target_catalogue.Max(var_name2);
2459 
2460  double binSize1_inv = pow((V1max-V1min)/nbin1,-1);
2461  double binSize2_inv = pow((V2max-V2min)/nbin2,-1);
2462 
2463  for (size_t i=0; i<input_var1.size(); i++)
2464  if ( (input_var1[i] < V1max && V1min < input_var1[i]) && (input_var2[i] < V2max && V2min < input_var2[i])) {
2465  int occ1 = max(0, min(int((input_var1[i]-V1min)*binSize1_inv), nbin1));
2466  int occ2 = max(0, min(int((input_var2[i]-V2min)*binSize2_inv), nbin2));
2467  fvars_input[occ1][occ2] ++;
2468  }
2469 
2470  for (size_t i=0; i<target_var1.size(); i++)
2471  if ( (target_var1[i] < V1max && V1min < target_var2[i]) && (target_var2[i] < V2max && V2min < target_var2[i])) {
2472  int occ1 = max(0, min(int((target_var1[i]-V1min)*binSize1_inv), nbin1));
2473  int occ2 = max(0, min(int((target_var2[i]-V2min)*binSize2_inv), nbin2));
2474  fvars_target[occ1][occ2] ++;
2475  }
2476 
2477 
2478  random::UniformRandomNumbers ran(0., 1., seed);
2479 
2480  for (size_t i=0; i<input_var1.size(); i++)
2481  if ((input_var1[i] < V1max && V1min < input_var1[i]) && (input_var2[i] < V2max && V2min < input_var2[i])) {
2482  const int occ1 = max(0, min(int((input_var1[i]-V1min)*binSize1_inv), nbin1));
2483  const int occ2 = max(0, min(int((input_var2[i]-V2min)*binSize2_inv), nbin2));
2484  if (ran() < fvars_target[occ1][occ2]/fvars_input[occ1][occ2])
2485  m_object.push_back(shared_ptr<Object>(input_catalogue.catalogue_object(i)));
2486  }
2487 
2488 }
2489 
2490 
2491 // ============================================================================
2492 
2493 
2494 void cbl::catalogue::Catalogue::remove_objects (const std::vector<bool> index)
2495 {
2496  if (index.size() != m_object.size()) ErrorCBL ("argument size not valid!", "remove_objects", "Catalogue.cpp");
2497 
2498  decltype(m_object) object_temp;
2499 
2500  for (size_t ii = 0; ii<index.size(); ii++)
2501  if (!index[ii]) object_temp.emplace_back(m_object[ii]);
2502 
2503  m_object.swap(object_temp);
2504 }
2505 
2506 
2507 // ============================================================================
2508 
2509 
2510 void cbl::catalogue::Catalogue::swap_objects (const int ind1, const int ind2)
2511 {
2512  shared_ptr<Object> temp = m_object[ind1];
2513  m_object[ind1] = m_object[ind2];
2514  m_object[ind2] = temp;
2515 }
2516 
2517 
2518 // ============================================================================
2519 
2520 
2521 void cbl::catalogue::Catalogue::sort (const Var var_name, const bool increasing)
2522 {
2523  coutCBL << "I'm sorting the catalogue..." << endl;
2524 
2525  vector<double> variable = var(var_name);
2526  bool swap = true;
2527  while (swap) {
2528  swap = false;
2529  for (size_t i = 0; i<nObjects()-1; ++i) {
2530  if (increasing) {
2531  if (variable[i] > variable[i+1]) {
2532  double temp = variable[i];
2533  variable[i] = variable[i+1];
2534  variable[i+1] = temp;
2535  swap_objects(i, i+1);
2536  swap = true;
2537  }
2538  }
2539  else {
2540  if (variable[i] < variable[i+1]) {
2541  double temp = variable[i];
2542  variable[i] = variable[i+1];
2543  variable[i+1] = temp;
2544  swap_objects(i, i+1);
2545  swap = true;
2546  }
2547  }
2548  }
2549  }
2550 }
2551 
2552 
2553 // ============================================================================
2554 
2555 
2556 void cbl::catalogue::Catalogue::shuffle (const int seed)
2557 {
2558  coutCBL << "I'm shuffling objects in the catalogue..." << endl;
2559 
2561  std::mt19937_64 generator;
2562  generator.seed(seed);
2563 
2564  std::shuffle(m_object.begin(), m_object.end(), generator);
2565 }
2566 
2567 
2568 
2569 // ============================================================================
2570 
2571 std::vector<double> cbl::catalogue::Catalogue::compute_catalogueProperties_box (const double boxside)
2572 {
2573  std::vector<double> prop(5);
2574  double vol = volume(boxside);
2575  double ndensity = m_object.size()/vol;
2576  double mps = pow(ndensity, -1./3.);
2577  double numdensity_error = pow(m_object.size(), 1./2)/vol;
2578  double mps_error = mps*(1./3)*(numdensity_error/ndensity);
2579 
2580  coutCBL << "Sample volume = " << vol << " (Mpc/h)^3" << endl;
2581  coutCBL << "Sample density = " << ndensity << " \u00b1 " << numdensity_error << " (Mpc/h)^-3" << endl;
2582  coutCBL << "Sample mps = " << mps << " \u00b1 " << mps_error << " Mpc/h" << endl;
2583 
2584  prop[0]=vol;
2585  prop[1]=ndensity;
2586  prop[2]=mps;
2587  prop[3]=numdensity_error;
2588  prop[4]=mps_error;
2589 
2590  return prop;
2591 }
2592 
2593 // ============================================================================
2594 
2595 std::vector<std::vector<double>> cbl::catalogue::Catalogue::compute_catalogueProperties_lightCone (cbl::cosmology::Cosmology cosmology, const std::vector<double> RA_range, const std::vector<double> DEC_range, const unsigned int nbin)
2596 {
2597  std::vector<vector<double>> prop(7, vector<double>(nbin));
2598  prop[0] = z_bins(nbin);
2599  std::vector<double> bin_limits = linear_bin_vector(nbin+1, cbl::Min(var(Var::_Redshift_)), cbl::Max(var(Var::_Redshift_)));
2600 
2601  double THETA_min = -DEC_range[1]+par::pi/2, THETA_max = -DEC_range[0]+par::pi/2;
2602  double delta_PHI = RA_range[1]-RA_range[0];
2603  double vol=0.;
2604 
2605  for (int i=0; i<(int)nbin; i++) {
2606  auto temp_cat = sub_catalogue(Var::_Redshift_, bin_limits[i], bin_limits[i+1]);
2607  prop[1][i]=temp_cat.nObjects(); //number of objects
2608  prop[2][i]=(-(std::cos(THETA_max)-std::cos(THETA_min))*delta_PHI)/(4*cbl::par::pi)*
2609  (volume_sphere(cosmology.D_C(bin_limits[i+1]))-volume_sphere(cosmology.D_C(bin_limits[i]))); //volume
2610  prop[3][i]=temp_cat.nObjects()/prop[2][i]; //numdensity
2611  prop[4][i]=pow(prop[3][i], -1./3.); //mps
2612  vol = vol + prop[2][i];
2613  prop[5][i]=sqrt(prop[1][0])/prop[2][i]; //numdensity error
2614  prop[6][i]=prop[4][i]*(1./3)*(prop[5][i]/prop[3][i]); //mps error
2615  }
2616 
2617  double numdensity = m_object.size()/vol;
2618  double mps = pow(numdensity, -1./3.);
2619  double numdensity_error = pow(m_object.size(), 1./2)/vol;
2620  double mps_error = mps*(1./3)*(numdensity_error/numdensity);
2621  cout << endl;
2622  coutCBL << "Sample volume = " << vol << " (Mpc/h)^3" << endl;
2623  coutCBL << "Sample density = " << numdensity << " \u00b1 " << numdensity_error << " (Mpc/h)^-3" << endl;
2624  coutCBL << "Sample mps = " << mps << " \u00b1 " << mps_error << " Mpc/h" << endl;
2625  coutCBL << "In addiction, volume, density and mps (with errors) have been calculated for " << nbin << " slices in redshift, from z = " << bin_limits[0] << " to z = " << bin_limits[nbin] << endl;
2626  cout << endl;
2627 
2628  return prop;
2629 }
2630 
2631 
Catalogue.h
The class Catalogue
Field3D.h
The class field3D.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
Object.h
The class Object.
cbl::Path
The class Path.
Definition: Path.h:145
cbl::Path::DirCosmo
std::string DirCosmo()
get the directory where the CosmoBolognaLbi are stored
Definition: Path.h:178
cbl::catalogue::Catalogue
The class Catalogue.
Definition: Catalogue.h:654
cbl::catalogue::Catalogue::Catalogue
Catalogue()=default
default constructor
cbl::catalogue::Catalogue::nRegions
size_t nRegions()
get the private member m_nRegions
Definition: Catalogue.cpp:396
cbl::catalogue::Catalogue::remove_objects
void remove_objects()
remove all objects remove all objects of the catalogue
Definition: Catalogue.h:2693
cbl::catalogue::Catalogue::add_objects
void add_objects(std::vector< std::shared_ptr< Object > > sample)
add some objects to the catalogue
Definition: Catalogue.h:2654
cbl::catalogue::Catalogue::restoreComovingCoordinates
void restoreComovingCoordinates()
restore comoving coordinates
Definition: Catalogue.cpp:1697
cbl::catalogue::Catalogue::Min
double Min(const Var var_name) const
get the minimum value of a variable of the catalogue objects
Definition: Catalogue.h:2287
cbl::catalogue::Catalogue::normalizeComovingCoordinates
void normalizeComovingCoordinates()
normalize comoving coordinates
Definition: Catalogue.cpp:1684
cbl::catalogue::Catalogue::sub_catalogue
Catalogue sub_catalogue(const Var var_name, const double down, const double up, const bool excl=false) const
create a sub-catalogue
Definition: Catalogue.cpp:1890
cbl::catalogue::Catalogue::weightedN_condition
double weightedN_condition(const Var var_name, const double down, const double up, const bool excl=false)
return the weighted number of objectes following a condition on the variable VAR
Definition: Catalogue.cpp:2248
cbl::catalogue::Catalogue::mangle_cut
Catalogue mangle_cut(const std::string mangle_mask, const bool excl=false) const
create a sub-catalogue
Definition: Catalogue.cpp:2027
cbl::catalogue::Catalogue::nObjects
size_t nObjects() const
get the number of objects of the catalogue
Definition: Catalogue.h:2264
cbl::catalogue::Catalogue::set_region_number
void set_region_number(const size_t nRegions)
set the private variable m_nRegion
Definition: Catalogue.cpp:944
cbl::catalogue::Catalogue::set_ra_dec_tile
void set_ra_dec_tile(const std::vector< double > RA_tile, const std::vector< double > Dec_tile, const CoordinateUnits inputUnits=CoordinateUnits::_degrees_)
set the central R.A. and Dec of the tiles
Definition: Catalogue.cpp:957
cbl::catalogue::Catalogue::write_comoving_coordinates
void write_comoving_coordinates(const std::string outputFile) const
write the comoving coordinates of the catalogue to an output file
Definition: Catalogue.cpp:1765
cbl::catalogue::Catalogue::isSetVar
bool isSetVar(const int index, const Var var_name) const
check if the given variable of the i-th object is set
Definition: Catalogue.cpp:740
cbl::catalogue::Catalogue::ra_tile
double ra_tile(const int i) const
get the private member Catalogue::m_object[i]->m_ra_tile
Definition: Catalogue.h:1870
cbl::catalogue::Catalogue::add_object
void add_object(std::shared_ptr< Object > object)
add one single object to the catalogue
Definition: Catalogue.h:2641
cbl::catalogue::Catalogue::catalogue_object
std::shared_ptr< Object > catalogue_object(const int i) const
get the i-th object of the catalogue
Definition: Catalogue.h:2236
cbl::catalogue::Catalogue::var_distr
void var_distr(const Var var_name, std::vector< double > &_var, std::vector< double > &dist, std::vector< double > &err, const int nbin, const bool linear=true, const std::string file_out=par::defaultString, const double Volume=1., const bool norm=false, const double V1=par::defaultDouble, const double V2=par::defaultDouble, const std::string bin_type="Linear", const bool convolution=false, const double sigma=0.) const
get the distribution of a variable
Definition: Catalogue.cpp:1544
cbl::catalogue::Catalogue::set_field
void set_field(const std::vector< std::string > field)
set a private variable
Definition: Catalogue.cpp:969
cbl::catalogue::Catalogue::write_obs_coordinates
void write_obs_coordinates(const std::string outputFile) const
write the polar coordinates of the catalogue to an output file
Definition: Catalogue.cpp:1783
cbl::catalogue::Catalogue::var
double var(const int index, const Var var_name) const
get the value of the i-th object variable
Definition: Catalogue.cpp:449
cbl::catalogue::Catalogue::compute_catalogueProperties_lightCone
std::vector< std::vector< double > > compute_catalogueProperties_lightCone(cbl::cosmology::Cosmology cosmology, const std::vector< double > RA_range, const std::vector< double > DEC_range, const unsigned int nbin)
compute catalogue volume, number density and mean particle separation in light cones
Definition: Catalogue.cpp:2595
cbl::catalogue::Catalogue::sort
void sort(const Var var_name, const bool increasing=false)
bubble sort of a catalogue wrt a variable
Definition: Catalogue.cpp:2521
cbl::catalogue::Catalogue::density_field
data::ScalarField3D density_field(const double cell_size, const Catalogue mask_catalogue, const int interpolation_type=0, const double kernel_radius=0., const bool useMass=false) const
return the density field from object position
Definition: Catalogue.cpp:2341
cbl::catalogue::Catalogue::angsep_xyz
double angsep_xyz(const int i, std::shared_ptr< Object > obj) const
get the angular distrance between the i-th object of the catalogue and another object
Definition: Catalogue.cpp:2114
cbl::catalogue::Catalogue::computeComovingCoordinates
void computeComovingCoordinates(const cosmology::Cosmology &cosm, const CoordinateUnits inputUnits=CoordinateUnits::_radians_)
compute the comoving coordinates (x, y, z) from the observed coordinates (R.A., Dec,...
Definition: Catalogue.cpp:1553
cbl::catalogue::Catalogue::nObjects_condition
int nObjects_condition(const Var var_name, const double down, const double up, const bool excl=false)
return the number of objectes following a condition on the variable VAR
Definition: Catalogue.cpp:2229
cbl::catalogue::Catalogue::Max
double Max(const Var var_name) const
get the maximum value of a variable of the catalogue objects
Definition: Catalogue.h:2295
cbl::catalogue::Catalogue::weightedN
double weightedN() const
get the total weight of the objects of the catalogue
Definition: Catalogue.cpp:1753
cbl::catalogue::Catalogue::replace_objects
void replace_objects(std::vector< T > sample)
replace existing objects with new ones
Definition: Catalogue.h:2674
cbl::catalogue::Catalogue::region
long region(const int i) const
get the private member Catalogue::m_object[i]->m_region
Definition: Catalogue.h:1954
cbl::catalogue::Catalogue::computePolarCoordinates
void computePolarCoordinates(const CoordinateUnits outputUnits=CoordinateUnits::_radians_)
compute the polar coordinates (R.A., Dec, dc) from the comoving coordinates (x, y,...
Definition: Catalogue.cpp:1587
cbl::catalogue::Catalogue::distance
double distance(const int i, std::shared_ptr< Object > obj) const
get the distrance between the i-th object of the catalogue and another object
Definition: Catalogue.cpp:2103
cbl::catalogue::Catalogue::set_region
void set_region(const std::vector< long > region, const int nRegions=-1)
set a private variable
Definition: Catalogue.cpp:933
cbl::catalogue::Catalogue::set_var
void set_var(const int index, const Var var_name, const double value, const cosmology::Cosmology cosmology={}, const bool update_coordinates=true)
set a private variable
Definition: Catalogue.cpp:978
cbl::catalogue::Catalogue::compute_catalogueProperties_box
std::vector< double > compute_catalogueProperties_box(const double boxside)
compute catalogue volume, number density and mean particle separation
Definition: Catalogue.cpp:2571
cbl::catalogue::Catalogue::field
std::vector< std::string > field() const
get the values of the object fields
Definition: Catalogue.cpp:436
cbl::catalogue::Catalogue::stats_var
std::vector< double > stats_var(const Var var_name) const
get the mean, the median, the standard deviation, and the difference between the third and first quar...
Definition: Catalogue.cpp:1521
cbl::catalogue::Catalogue::write_data
void write_data(const std::string outputFile, const std::vector< Var > var_name={}, const std::string sep=" ", const std::string header="") const
write both the comoving and polar coordinates, and the regions (if present) of the catalogue to an ou...
Definition: Catalogue.cpp:1804
cbl::catalogue::Catalogue::smooth
std::shared_ptr< Catalogue > smooth(const double gridsize, const cosmology::Cosmology cosmology, const std::vector< Var > vars={}, const int SUB=1)
create a smoothed version of the catalogue averaging quantities on a X, Y, Z grid
Definition: Catalogue.cpp:2125
cbl::catalogue::Catalogue::swap_objects
void swap_objects(const int ind1, const int ind2)
swap two existing objects
Definition: Catalogue.cpp:2510
cbl::catalogue::Catalogue::dec_tile
double dec_tile(const int i) const
get the private member Catalogue::m_object[i]->m_dec_tile
Definition: Catalogue.h:1877
cbl::catalogue::Catalogue::var_int
int var_int(const int index, const Var var_name) const
get the value of the i-th object integer variable
Definition: Catalogue.cpp:697
cbl::catalogue::Catalogue::region_list
std::vector< long > region_list() const
get the list of regions in which the catalogue is divided
Definition: Catalogue.cpp:411
cbl::catalogue::Catalogue::Order
void Order()
restore the original vector (i.e. the opposite of Order(std::vector<int>))
Definition: Catalogue.cpp:1733
cbl::catalogue::Catalogue::shuffle
void shuffle(const int seed)
shuffle objects in the catalogue
Definition: Catalogue.cpp:2556
cbl::catalogue::Catalogue::region
std::vector< long > region() const
get the values of the object regions
Definition: Catalogue.cpp:423
cbl::catalogue::Catalogue::counts_in_cell
data::ScalarField3D counts_in_cell(const double cell_size, const int interpolation_type=0, const bool useMass=false, const double minX=par::defaultDouble, const double maxX=par::defaultDouble, const double minY=par::defaultDouble, const double maxY=par::defaultDouble, const double minZ=par::defaultDouble, const double maxZ=par::defaultDouble) const
return the density field from object positions
Definition: Catalogue.cpp:2266
cbl::catalogue::Catalogue::diluted_catalogue
Catalogue diluted_catalogue(const double nSub, const int seed=3213) const
create a diluted catalogue
Definition: Catalogue.cpp:2076
cbl::catalogue::ChainMeshCell
The class ChainMeshCell.
Definition: ChainMeshCell.h:53
cbl::catalogue::Cluster
The class Cluster.
Definition: Cluster.h:55
cbl::catalogue::Galaxy
The class Galaxy.
Definition: Galaxy.h:54
cbl::catalogue::Halo
The class CatalogueChainMesh.
Definition: Halo.h:53
cbl::catalogue::HostHalo
The class HostHalo.
Definition: HostHalo.h:53
cbl::catalogue::Mock
The class Mock.
Definition: Mock.h:53
cbl::catalogue::Object
The class Object.
Definition: Object.h:132
cbl::catalogue::Object::set_xx
void set_xx(const double xx)
set the member m_xx
Definition: Object.h:1261
cbl::catalogue::RandomObject
The class RandomObject.
Definition: RandomObject.h:55
cbl::catalogue::Void
The class Void.
Definition: Void.h:52
cbl::chainmesh::ChainMesh3D
The class ChainMesh3D.
Definition: ChainMesh.h:392
cbl::chainmesh::ChainMesh::nCell
long nCell() const
get the private member ChainMesh::m_nCell_tot
Definition: ChainMesh.h:156
cbl::chainmesh::ChainMesh::get_list
std::vector< long > get_list(const long cell_index) const
get the index of the object inside a cell
Definition: ChainMesh.cpp:526
cbl::cosmology::Cosmology
The class Cosmology.
Definition: Cosmology.h:277
cbl::cosmology::Cosmology::Redshift
double Redshift(const double d_c=1., const double z1_guess=0., const double z2_guess=10., const double prec=0.0001) const
redshift at a given comoving distance
Definition: Cosmology.cpp:1045
cbl::cosmology::Cosmology::D_C
double D_C(const double redshift) const
the comoving line-of-sight distance at a given redshift
Definition: Cosmology.cpp:741
cbl::data::Field3D::deltaX
double deltaX() const
get the private member m_deltaX
Definition: Field3D.h:319
cbl::data::Field3D::MaxY
double MaxY() const
get the private member m_MaxY
Definition: Field3D.h:307
cbl::data::Field3D::MinY
double MinY() const
get the private member m_MinY
Definition: Field3D.h:289
cbl::data::Field3D::XX
double XX(const int i) const
get the value of the X coordinates at the i-th cell
Definition: Field3D.h:346
cbl::data::Field3D::deltaZ
double deltaZ() const
get the private member m_deltaZ
Definition: Field3D.h:331
cbl::data::Field3D::ny
int ny() const
get the private member m_nY
Definition: Field3D.h:253
cbl::data::Field3D::MinX
double MinX() const
get the private member m_MinX
Definition: Field3D.h:283
cbl::data::Field3D::nx
int nx() const
get the private member m_nX
Definition: Field3D.h:247
cbl::data::Field3D::MaxZ
double MaxZ() const
get the private member m_MaxZ
Definition: Field3D.h:313
cbl::data::Field3D::ZZ
double ZZ(const int i) const
get the value of the Z coordinates at the i-th cell
Definition: Field3D.h:364
cbl::data::Field3D::deltaY
double deltaY() const
get the private member m_deltaY
Definition: Field3D.h:325
cbl::data::Field3D::nz
int nz() const
get the private member m_nZ
Definition: Field3D.h:259
cbl::data::Field3D::MinZ
double MinZ() const
get the private member m_MinZ
Definition: Field3D.h:295
cbl::data::Field3D::YY
double YY(const int i) const
get the value of the Y coordinates at the i-th cell
Definition: Field3D.h:355
cbl::data::Field3D::MaxX
double MaxX() const
get the private member m_MaxX
Definition: Field3D.h:301
cbl::data::ScalarField3D
The class ScalarField3D.
Definition: Field3D.h:648
cbl::data::ScalarField3D::ScalarField
double ScalarField(const int i, const int j, const int k) const
get the value of the scalar field
Definition: Field3D.cpp:316
cbl::data::ScalarField3D::GaussianConvolutionField
void GaussianConvolutionField(const double kernel_size)
perform a smoothing of the field with a gaussian kernel
Definition: Field3D.cpp:257
cbl::data::ScalarField3D::set_ScalarField
void set_ScalarField(const double value, const int i, const int j, const int k, const bool add=0)
set the value of the scalar field
Definition: Field3D.cpp:289
cbl::random::UniformRandomNumbers
The class UniformRandomNumbers.
Definition: RandomNumbers.h:254
cbl::par::fINT
static const char fINT[]
conversion symbol for: int -> std::string
Definition: Constants.h:121
cbl::par::defaultDouble
static const double defaultDouble
default double value
Definition: Constants.h:348
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::mask_function
std::function< bool(const std::shared_ptr< Object > obj)> mask_function
Definition of a new type to manage mask function.
Definition: Catalogue.h:266
cbl::catalogue::Var
Var
the catalogue variables
Definition: Catalogue.h:70
cbl::catalogue::ObjectType
ObjectType
the object types
Definition: Object.h:51
cbl::catalogue::CharEncode
CharEncode
character encoding of input file
Definition: Catalogue.h:476
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::Min
T Min(const std::vector< T > vect)
minimum element of a std::vector
Definition: Kernel.h:1324
cbl::conv
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
cbl::Average
double Average(const std::vector< double > vect)
the average of a std::vector
Definition: Func.cpp:870
cbl::CoordinateType
CoordinateType
the coordinate type
Definition: Kernel.h:624
cbl::CoordinateType::_observed_
@ _observed_
observed coordinates (R.A., Dec, redshift)
cbl::CoordinateType::_comoving_
@ _comoving_
comoving coordinates (x, y, z)
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::volume_sphere
T volume_sphere(const T RR)
the volume of a sphere of a given radius
Definition: Func.h:1231
cbl::j2
double j2(const double xx)
the l=2 spherical Bessel function
Definition: Func.cpp:2048
cbl::isSet
bool isSet(const std::string var)
check if the value of a [string] variable has already been set
Definition: Kernel.h:803
cbl::linear_bin_vector
std::vector< T > linear_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with linearly spaced values
Definition: Kernel.h:1604
cbl::checkDim
void checkDim(const std::vector< T > vect, const int val, const std::string vector, bool equal=true)
check if the dimension of a std::vector is equal/lower than an input value
Definition: Kernel.h:1532
cbl::checkIO
void checkIO(const std::ifstream &fin, const std::string file="NULL")
check if an input file can be opened
Definition: Kernel.cpp:160
cbl::cartesian_coord
void cartesian_coord(const double ra, const double dec, const double dd, double &XX, double &YY, double &ZZ)
conversion from polar coordinates to Cartesian coordinates
Definition: Func.cpp:221
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::Max
T Max(const std::vector< T > vect)
maximum element of a std::vector
Definition: Kernel.h:1336
cbl::transpose
std::vector< std::vector< T > > transpose(std::vector< std::vector< T >> matrix)
transpose a matrix
Definition: Kernel.h:1729
cbl::radians
double radians(const double angle, const CoordinateUnits inputUnits=CoordinateUnits::_degrees_)
conversion to radians
Definition: Func.cpp:126
cbl::Sigma
double Sigma(const std::vector< double > vect)
the standard deviation of a std::vector
Definition: Func.cpp:925
cbl::Quartile
std::vector< double > Quartile(const std::vector< double > Vect)
the first, second and third quartiles of a std::vector
Definition: Func.cpp:1002
cbl::volume
double volume(const double boxSize, const int frac, const double Bord, const double mean_redshift, cosmology::Cosmology &real_cosm)
get the volume of a simulation box
Definition: Func.cpp:78
cbl::degrees
double degrees(const double angle, const CoordinateUnits inputUnits=CoordinateUnits::_radians_)
conversion to degrees
Definition: Func.cpp:104
cbl::CoordinateUnits
CoordinateUnits
the coordinate units
Definition: Kernel.h:562
cbl::WarningMsgCBL
void WarningMsgCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL)
internal CBL warning message
Definition: Kernel.h:747
cbl::convolution
void convolution(const std::vector< double > f1, const std::vector< double > f2, std::vector< double > &res, const double deltaX)
convolution of two functions
Definition: Func.cpp:1580
cbl::arcseconds
double arcseconds(const double angle, const CoordinateUnits inputUnits=CoordinateUnits::_radians_)
conversion to arcseconds
Definition: Func.cpp:148
cbl::polar_coord
void polar_coord(const double XX, const double YY, const double ZZ, double &ra, double &dec, double &dd)
conversion from Cartesian coordinates to polar coordinates
Definition: Func.cpp:214
cbl::arcminutes
double arcminutes(const double angle, const CoordinateUnits inputUnits=CoordinateUnits::_radians_)
conversion to arcminutes
Definition: Func.cpp:170
cbl::catalogue::MaskObject
object that encapsulate the mask function.
Definition: Catalogue.h:277