43 using namespace catalogue;
44 using namespace chainmesh;
46 using namespace pairs;
47 using namespace measure;
48 using namespace twopt;
54 void cbl::measure::twopt::TwoPointCorrelation2D_cartesian::set_parameters (
const BinType binType_rp,
const double rpMin,
const double rpMax,
const int nbins_rp,
const double shift_rp,
const BinType binType_pi,
const double piMin,
const double piMax,
const int nbins_pi,
const double shift_pi,
const CoordinateUnits angularUnits, std::function<
double(
double)> angularWeight,
const bool compute_extra_info)
56 if (!compute_extra_info) {
57 if (binType_rp==BinType::_logarithmic_) {
58 if (binType_pi==BinType::_logarithmic_) {
59 m_dd = move(Pair::Create(PairType::_comovingCartesian_loglog_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits, angularWeight));
62 m_dd = move(Pair::Create(PairType::_comovingCartesian_loglin_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits, angularWeight));
66 if (binType_pi==BinType::_logarithmic_) {
67 m_dd = move(Pair::Create(PairType::_comovingCartesian_linlog_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits, angularWeight));
70 m_dd = move(Pair::Create(PairType::_comovingCartesian_linlin_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits, angularWeight));
75 if (binType_rp==BinType::_logarithmic_) {
76 if (binType_pi==BinType::_logarithmic_) {
77 m_dd = move(Pair::Create(PairType::_comovingCartesian_loglog_, PairInfo::_extra_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits, angularWeight));
80 m_dd = move(Pair::Create(PairType::_comovingCartesian_loglin_, PairInfo::_extra_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits, angularWeight));
84 if (binType_pi==BinType::_logarithmic_) {
85 m_dd = move(Pair::Create(PairType::_comovingCartesian_linlog_, PairInfo::_extra_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits, angularWeight));
88 m_dd = move(Pair::Create(PairType::_comovingCartesian_linlin_, PairInfo::_extra_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits, angularWeight));
93 if (binType_rp==BinType::_logarithmic_) {
94 if (binType_pi==BinType::_logarithmic_) {
95 m_rr = move(Pair::Create(PairType::_comovingCartesian_loglog_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits));
96 m_dr = move(Pair::Create(PairType::_comovingCartesian_loglog_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits));
99 m_rr = move(Pair::Create(PairType::_comovingCartesian_loglin_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits));
100 m_dr = move(Pair::Create(PairType::_comovingCartesian_loglin_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits));
104 if (binType_pi==BinType::_logarithmic_) {
105 m_rr = move(Pair::Create(PairType::_comovingCartesian_linlog_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits));
106 m_dr = move(Pair::Create(PairType::_comovingCartesian_linlog_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits));
109 m_rr = move(Pair::Create(PairType::_comovingCartesian_linlin_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits));
110 m_dr = move(Pair::Create(PairType::_comovingCartesian_linlin_, PairInfo::_standard_, rpMin, rpMax, nbins_rp, shift_rp, piMin, piMax, nbins_pi, shift_pi, angularUnits));
119 void cbl::measure::twopt::TwoPointCorrelation2D_cartesian::set_parameters (
const BinType binType_rp,
const double rpMin,
const double rpMax,
const double binSize_rp,
const double shift_rp,
const BinType binType_pi,
const double piMin,
const double piMax,
const double binSize_pi,
const double shift_pi,
const CoordinateUnits angularUnits, std::function<
double(
double)> angularWeight,
const bool compute_extra_info)
121 if (!compute_extra_info) {
122 if (binType_rp==BinType::_logarithmic_) {
123 if (binType_pi==BinType::_logarithmic_) {
124 m_dd = move(Pair::Create(PairType::_comovingCartesian_loglog_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits, angularWeight));
127 m_dd = move(Pair::Create(PairType::_comovingCartesian_loglin_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits, angularWeight));
131 if (binType_pi==BinType::_logarithmic_) {
132 m_dd = move(Pair::Create(PairType::_comovingCartesian_linlog_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits, angularWeight));
135 m_dd = move(Pair::Create(PairType::_comovingCartesian_linlin_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits, angularWeight));
140 if (binType_rp==BinType::_logarithmic_) {
141 if (binType_pi==BinType::_logarithmic_) {
142 m_dd = move(Pair::Create(PairType::_comovingCartesian_loglog_, PairInfo::_extra_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits, angularWeight));
145 m_dd = move(Pair::Create(PairType::_comovingCartesian_loglin_, PairInfo::_extra_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits, angularWeight));
149 if (binType_pi==BinType::_logarithmic_) {
150 m_dd = move(Pair::Create(PairType::_comovingCartesian_linlog_, PairInfo::_extra_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits, angularWeight));
153 m_dd = move(Pair::Create(PairType::_comovingCartesian_linlin_, PairInfo::_extra_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits, angularWeight));
159 if (binType_rp==BinType::_logarithmic_) {
160 if (binType_pi==BinType::_logarithmic_) {
161 m_rr = move(Pair::Create(PairType::_comovingCartesian_loglog_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits));
162 m_dr = move(Pair::Create(PairType::_comovingCartesian_loglog_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits));
165 m_rr = move(Pair::Create(PairType::_comovingCartesian_loglin_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits));
166 m_dr = move(Pair::Create(PairType::_comovingCartesian_loglin_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits));
170 if (binType_pi==BinType::_logarithmic_) {
171 m_rr = move(Pair::Create(PairType::_comovingCartesian_linlog_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits));
172 m_dr = move(Pair::Create(PairType::_comovingCartesian_linlog_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits));
175 m_rr = move(Pair::Create(PairType::_comovingCartesian_linlin_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits));
176 m_dr = move(Pair::Create(PairType::_comovingCartesian_linlin_, PairInfo::_standard_, rpMin, rpMax, binSize_rp, shift_rp, piMin, piMax, binSize_pi, shift_pi, angularUnits));
187 m_dataset->read(dir+file);
196 vector<double> xx = m_dataset->xx(), yy = m_dataset->yy();
200 string header =
"[1] perpendicular separation at the bin centre # [2] parallel separation at the bin centre # [3] 2D two-point correlation function # [4] error";
201 if (m_compute_extra_info) header +=
" # [5] mean perpendicular separation # [6] standard deviation of the distribution of perpendicular separations # [7] mean parallel separation # [8] standard deviation of the distribution of parallel separations # [9] mean redshift # [10] standard deviation of the redshift distribution";
203 m_dataset->write(dir, file, header, full, 5, rank);
210 void cbl::measure::twopt::TwoPointCorrelation2D_cartesian::measure (
const ErrorType errorType,
const std::string dir_output_pairs,
const std::vector<std::string> dir_input_pairs,
const std::string dir_output_resample,
const int nMocks,
const bool count_dd,
const bool count_rr,
const bool count_dr,
const bool tcount,
const Estimator estimator,
const double fact,
const int seed)
213 case (ErrorType::_Poisson_) :
214 measurePoisson(dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
216 case (ErrorType::_Jackknife_) :
217 measureJackknife(dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact);
219 case (ErrorType::_Bootstrap_) :
220 measureBootstrap(nMocks, dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact, seed);
223 ErrorCBL(
"unknown type of error!",
"measure",
"TwoPointCorrelation2D_cartesian.cpp");
235 count_allPairs(m_twoPType, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
240 if (estimator==Estimator::_natural_)
241 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
242 else if (estimator==Estimator::_LandySzalay_) {
243 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
246 ErrorCBL(
"the chosen estimator is not implemented!",
"measurePoisson",
"TwoPointCorrelation2D_cartesian.cpp");
257 string mkdir =
"mkdir -p "+dir_output_resample;
258 if (system(mkdir.c_str())) {}
261 vector<long> region_list = m_data->region_list();
262 size_t nRegions = region_list.size();
264 vector< vector<double > > xi_SubSamples, covariance;
266 vector<shared_ptr<Pair> > dd_regions, rr_regions, dr_regions;
267 count_allPairs_region(dd_regions, rr_regions, dr_regions, m_twoPType, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
269 vector<shared_ptr<Data> > data_SS = (estimator==Estimator::_natural_) ? XiJackknife(dd_regions, rr_regions) : XiJackknife(dd_regions, rr_regions, dr_regions);
271 for (
size_t i=0; i<nRegions; i++) {
275 string header =
"[1] perpendicular separation at the bin centre # [2] parallel separation at the bin centre # [3] 2D two-point correlation function # [4] error";
276 if (m_compute_extra_info) header +=
" # [5] mean perpendicular separation # [6] standard deviation of the distribution of perpendicular separations # [7] mean parallel separation # [8] standard deviation of the distribution of parallel separations # [9] mean redshift # [10] standard deviation of the redshift distribution";
277 data_SS[i]->write(dir_output_resample, file, header,
false, 10, 0);
280 xi_SubSamples.push_back(data_SS[i]->data());
285 if (estimator==Estimator::_natural_)
286 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
287 else if (estimator==Estimator::_LandySzalay_)
288 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
290 ErrorCBL(
"the chosen estimator is not implemented!",
"measureJackknife",
"TwoPointCorrelation2D_cartesian.cpp");
292 m_dataset->set_covariance(covariance);
298 void cbl::measure::twopt::TwoPointCorrelation2D_cartesian::measureBootstrap (
const int nMocks,
const std::string dir_output_pairs,
const std::vector<std::string> dir_input_pairs,
const std::string dir_output_resample,
const bool count_dd,
const bool count_rr,
const bool count_dr,
const bool tcount,
const Estimator estimator,
const double fact,
const int seed)
301 ErrorCBL(
"the number of mocks must be >0",
"measureBootstrap",
"TwoPointCorrelation2D_cartesian.cpp");
304 string mkdir =
"mkdir -p "+dir_output_resample;
305 if (system(mkdir.c_str())) {}
308 vector< vector<double > > xi_SubSamples, covariance;
310 vector<shared_ptr<Pair> > dd_regions, rr_regions, dr_regions;
311 count_allPairs_region(dd_regions, rr_regions, dr_regions, m_twoPType, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
313 vector<shared_ptr<Data> > data_SS = (estimator==Estimator::_natural_) ? XiBootstrap(nMocks, dd_regions, rr_regions, seed) : XiBootstrap(nMocks, dd_regions, rr_regions, dr_regions, seed);
315 for (
int i=0; i<nMocks; i++) {
319 string header =
"[1] perpendicular separation at the bin centre # [2] parallel separation at the bin centre # [3] 2D two-point correlation function # [4] error";
320 if (m_compute_extra_info) header +=
" # [5] mean perpendicular separation # [6] standard deviation of the distribution of perpendicular separations # [7] mean parallel separation # [8] standard deviation of the distribution of parallel separations # [9] mean redshift # [10] standard deviation of the redshift distribution";
321 data_SS[i]->write(dir_output_resample, file, header,
false, 10, 0);
324 xi_SubSamples.push_back(data_SS[i]->data());
329 if (estimator==Estimator::_natural_)
330 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
331 else if (estimator==Estimator::_LandySzalay_)
332 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
334 ErrorCBL(
"the chosen estimator is not implemented!",
"measureBootstrap",
"TwoPointCorrelation2D_cartesian.cpp");
336 m_dataset->set_covariance(covariance);
The class TwoPointCorrelation2D_cartesian.
void measureJackknife(const std::string dir_output_pairs=par::defaultString, const std::vector< std::string > dir_input_pairs={}, const std::string dir_output_resample=par::defaultString, const bool count_dd=true, const bool count_rr=true, const bool count_dr=true, const bool tcount=true, const Estimator estimator=Estimator::_LandySzalay_, const double fact=0.1) override
measure the 2D two-point correlation function in Cartesian coordinates, estimating the covariance wit...
void measurePoisson(const std::string dir_output_pairs=par::defaultString, const std::vector< std::string > dir_input_pairs={}, const bool count_dd=true, const bool count_rr=true, const bool count_dr=true, const bool tcount=true, const Estimator estimator=Estimator::_LandySzalay_, const double fact=0.1) override
measure the 2D two-point correlation function in Cartesian coordinates, with Poisson errors
void measure(const ErrorType errorType=ErrorType::_Poisson_, const std::string dir_output_pairs=par::defaultString, const std::vector< std::string > dir_input_pairs={}, const std::string dir_output_resample=par::defaultString, const int nMocks=0, const bool count_dd=true, const bool count_rr=true, const bool count_dr=true, const bool tcount=true, const Estimator estimator=Estimator::_LandySzalay_, const double fact=0.1, const int seed=3213) override
measure the 2D two-point correlation function in Cartesian coordinates
void set_parameters(const BinType binType_rp, const double rpMin, const double rpMax, const int nbins_rp, const double shift_rp, const BinType binType_pi, const double piMin, const double piMax, const int nbins_pi, const double shift_pi, const CoordinateUnits angularUnits=CoordinateUnits::_radians_, std::function< double(double)> angularWeight=nullptr, const bool compute_extra_info=false)
set the binning parameters
void write(const std::string dir, const std::string file, const bool full, const int rank=0) const override
write the 2D two-point correlation function
void read(const std::string dir, const std::string file) override
read the 2D two-point correlation function
void measureBootstrap(const int nMocks, const std::string dir_output_pairs=par::defaultString, const std::vector< std::string > dir_input_pairs={}, const std::string dir_output_resample=par::defaultString, const bool count_dd=true, const bool count_rr=true, const bool count_dr=true, const bool tcount=true, const Estimator estimator=Estimator::_LandySzalay_, const double fact=0.1, const int seed=3213) override
measure the 2D two-point correlation function in Cartesian coordinates, estimating the covariance wit...
static const char fINT[]
conversion symbol for: int -> std::string
static const std::string defaultString
default std::string value
Estimator
the two-point correlation estimator
ErrorType
the two-point correlation function error type
The global namespace of the CosmoBolognaLib
std::string conv(const T val, const char *fact)
convert a number to a std::string
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
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
void covariance_matrix(const std::vector< std::vector< double >> mat, std::vector< std::vector< double >> &cov, const bool JK=false)
compute the covariance matrix from an input dataset
CoordinateUnits
the coordinate units