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_polar::set_parameters (
const BinType binType_rad,
const double rMin,
const double rMax,
const int nbins_rad,
const double shift_rad,
const BinType binType_mu,
const double muMin,
const double muMax,
const int nbins_mu,
const double shift_mu,
const CoordinateUnits angularUnits, std::function<
double(
double)> angularWeight,
const bool compute_extra_info)
56 if (muMin<0.)
ErrorCBL(
"mMun must be >0 !",
"set_parameters",
"TwoPointCorrelation2D_polar.cpp");
57 if (muMin>1.)
ErrorCBL(
"mMun must be <1 !",
"set_parameters",
"TwoPointCorrelation2D_polar.cpp");
58 if (rMin<0.)
ErrorCBL(
"rMun must be >0 !",
"set_parameters",
"TwoPointCorrelation2D_polar.cpp");
60 if (!compute_extra_info) {
61 if (binType_rad==BinType::_logarithmic_) {
62 if (binType_mu==BinType::_logarithmic_) {
63 m_dd = move(Pair::Create(PairType::_comovingPolar_loglog_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits, angularWeight));
66 m_dd = move(Pair::Create(PairType::_comovingPolar_loglin_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits, angularWeight));
70 if (binType_mu==BinType::_logarithmic_) {
71 m_dd = move(Pair::Create(PairType::_comovingPolar_linlog_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits, angularWeight));
74 m_dd = move(Pair::Create(PairType::_comovingPolar_linlin_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits, angularWeight));
79 if (binType_rad==BinType::_logarithmic_) {
80 if (binType_mu==BinType::_logarithmic_) {
81 m_dd = move(Pair::Create(PairType::_comovingPolar_loglog_,PairInfo::_extra_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits, angularWeight));
84 m_dd = move(Pair::Create(PairType::_comovingPolar_loglin_, PairInfo::_extra_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits, angularWeight));
88 if (binType_mu==BinType::_logarithmic_) {
89 m_dd = move(Pair::Create(PairType::_comovingPolar_linlog_, PairInfo::_extra_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits, angularWeight));
92 m_dd = move(Pair::Create(PairType::_comovingPolar_linlin_, PairInfo::_extra_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits, angularWeight));
98 if (binType_rad==BinType::_logarithmic_) {
99 if (binType_mu==BinType::_logarithmic_) {
100 m_rr = move(Pair::Create(PairType::_comovingPolar_loglog_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits));
101 m_dr = move(Pair::Create(PairType::_comovingPolar_loglog_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits));
104 m_rr = move(Pair::Create(PairType::_comovingPolar_loglin_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits));
105 m_dr = move(Pair::Create(PairType::_comovingPolar_loglin_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits));
109 if (binType_mu==BinType::_logarithmic_) {
110 m_rr = move(Pair::Create(PairType::_comovingPolar_linlog_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits));
111 m_dr = move(Pair::Create(PairType::_comovingPolar_linlog_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits));
114 m_rr = move(Pair::Create(PairType::_comovingPolar_linlin_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits));
115 m_dr = move(Pair::Create(PairType::_comovingPolar_linlin_, PairInfo::_standard_, rMin, rMax, nbins_rad, shift_rad, muMin, muMax, nbins_mu, shift_mu, angularUnits));
125 void cbl::measure::twopt::TwoPointCorrelation2D_polar::set_parameters (
const BinType binType_rad,
const double rMin,
const double rMax,
const double binSize_rad,
const double shift_rad,
const BinType binType_mu,
const double muMin,
const double muMax,
const double binSize_mu,
const double shift_mu,
const CoordinateUnits angularUnits, std::function<
double(
double)> angularWeight,
const bool compute_extra_info)
127 if (muMin<0.)
ErrorCBL(
"mMun must be >0 !",
"set_parameters",
"TwoPointCorrelation2D_polar.cpp");
128 if (muMin>1.)
ErrorCBL(
"mMun must be <1 !",
"set_parameters",
"TwoPointCorrelation2D_polar.cpp");
129 if (rMin<0.)
ErrorCBL(
"rMun must be >0 !",
"set_parameters",
"TwoPointCorrelation2D_polar.cpp");
131 if (!compute_extra_info) {
132 if (binType_rad==BinType::_logarithmic_) {
133 if (binType_mu==BinType::_logarithmic_) {
134 m_dd = move(Pair::Create(PairType::_comovingPolar_loglog_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits, angularWeight));
137 m_dd = move(Pair::Create(PairType::_comovingPolar_loglin_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits, angularWeight));
141 if (binType_mu==BinType::_logarithmic_) {
142 m_dd = move(Pair::Create(PairType::_comovingPolar_linlog_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits, angularWeight));
145 m_dd = move(Pair::Create(PairType::_comovingPolar_linlin_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits, angularWeight));
150 if (binType_rad==BinType::_logarithmic_) {
151 if (binType_mu==BinType::_logarithmic_) {
152 m_dd = move(Pair::Create(PairType::_comovingPolar_loglog_,PairInfo::_extra_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits, angularWeight));
155 m_dd = move(Pair::Create(PairType::_comovingPolar_loglin_, PairInfo::_extra_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits, angularWeight));
159 if (binType_mu==BinType::_logarithmic_) {
160 m_dd = move(Pair::Create(PairType::_comovingPolar_linlog_, PairInfo::_extra_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits, angularWeight));
163 m_dd = move(Pair::Create(PairType::_comovingPolar_linlin_, PairInfo::_extra_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits, angularWeight));
169 if (binType_rad==BinType::_logarithmic_) {
170 if (binType_mu==BinType::_logarithmic_) {
171 m_rr = move(Pair::Create(PairType::_comovingPolar_loglog_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits));
172 m_dr = move(Pair::Create(PairType::_comovingPolar_loglog_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits));
175 m_rr = move(Pair::Create(PairType::_comovingPolar_loglin_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits));
176 m_dr = move(Pair::Create(PairType::_comovingPolar_loglin_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits));
180 if (binType_mu==BinType::_logarithmic_) {
181 m_rr = move(Pair::Create(PairType::_comovingPolar_linlog_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits));
182 m_dr = move(Pair::Create(PairType::_comovingPolar_linlog_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits));
185 m_rr = move(Pair::Create(PairType::_comovingPolar_linlin_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits));
186 m_dr = move(Pair::Create(PairType::_comovingPolar_linlin_, PairInfo::_standard_, rMin, rMax, binSize_rad, shift_rad, muMin, muMax, binSize_mu, shift_mu, angularUnits));
197 m_dataset->read(dir+file);
206 vector<double> xx = m_dataset->xx(), yy = m_dataset->yy();
210 string header =
"[1] separation at the bin centre # [2] angular separation at the bin centre # [3] 2D two-point correlation function # [4] error";
211 if (m_compute_extra_info) header +=
" # [5] mean separation # [6] standard deviation of the separation distribution # [7] mean angular separation # [8] standard deviation of the angular separation distribution # [9] mean redshift # [10] standard deviation of the redshift distribution";
213 m_dataset->write(dir, file, header, full, 5, rank);
220 void cbl::measure::twopt::TwoPointCorrelation2D_polar::measure (
const ErrorType errorType,
const std::string dir_output_pairs,
const std::vector<std::string> dir_input_pairs,
const std::string dir_output_resample,
int nMocks,
bool count_dd,
const bool count_rr,
const bool count_dr,
const bool tcount,
const Estimator estimator,
const double fact,
const int seed)
223 case (ErrorType::_Poisson_) :
224 measurePoisson(dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
226 case (ErrorType::_Jackknife_) :
227 measureJackknife(dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact);
229 case (ErrorType::_Bootstrap_) :
230 measureBootstrap(nMocks, dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact, seed);
233 ErrorCBL(
"unknown type of error!",
"measure",
"TwoPointCorrelation2D_polar.cpp");
245 count_allPairs(m_twoPType, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
250 if (estimator==Estimator::_natural_)
251 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
252 else if (estimator==Estimator::_LandySzalay_)
253 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
255 ErrorCBL(
"the chosen estimator is not implemented!",
"measurePoisson",
"TwoPointCorrelation2D_polar.cpp");
266 string mkdir =
"mkdir -p "+dir_output_resample;
267 if (system(mkdir.c_str())) {}
270 vector<long> region_list = m_data->region_list();
271 size_t nRegions = region_list.size();
273 vector< vector<double > > xi_SubSamples, covariance;
275 vector<shared_ptr<Pair> > dd_regions, rr_regions, dr_regions;
276 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);
278 vector<shared_ptr<Data> > data_SS = (estimator==Estimator::_natural_) ? XiJackknife(dd_regions, rr_regions) : XiJackknife(dd_regions, rr_regions, dr_regions);
280 for (
size_t i=0; i<nRegions; i++) {
284 string header =
"[1] separation at the bin centre # [2] angular separation at the bin centre # [3] 2D two-point correlation function # [4] error";
285 if (m_compute_extra_info) header +=
" # [5] mean separation # [6] standard deviation of the separation distribution # [7] mean angular separation # [8] standard deviation of the angular separation distribution # [9] mean redshift # [10] standard deviation of the redshift distribution";
286 data_SS[i]->write(dir_output_resample, file, header,
false, 10, 0);
289 xi_SubSamples.push_back(data_SS[i]->data());
294 if (estimator==Estimator::_natural_)
295 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
296 else if (estimator==Estimator::_LandySzalay_)
297 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
299 ErrorCBL(
"the chosen estimator is not implemented!",
"measureJackknife",
"TwoPointCorrelation2D_polar.cpp");
301 m_dataset->set_covariance(covariance);
308 void cbl::measure::twopt::TwoPointCorrelation2D_polar::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)
311 ErrorCBL(
"the number of mocks must be >0!",
"measureBootstrap",
"TwoPointCorrelation2D_polar.cpp");
314 string mkdir =
"mkdir -p "+dir_output_resample;
315 if (system(mkdir.c_str())) {}
318 vector< vector<double > > xi_SubSamples, covariance;
320 vector<shared_ptr<Pair> > dd_regions, rr_regions, dr_regions;
321 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);
323 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);
325 for (
int i=0; i<nMocks; i++) {
329 string header =
"[1] separation at the bin centre # [2] angular separation at the bin centre # [3] 2D two-point correlation function # [4] error";
330 if (m_compute_extra_info) header +=
" # [5] mean separation # [6] standard deviation of the separation distribution # [7] mean angular separation # [8] standard deviation of the angular separation distribution # [9] mean redshift # [10] standard deviation of the redshift distribution";
331 data_SS[i]->write(dir_output_resample, file, header,
false, 10, 0);
334 xi_SubSamples.push_back(data_SS[i]->data());
340 if (estimator==Estimator::_natural_)
341 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
342 else if (estimator==Estimator::_LandySzalay_)
343 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
345 ErrorCBL(
"the chosen estimator is not implemented!",
"measureBootstrap",
"TwoPointCorrelation2D_polar.cpp");
347 m_dataset->set_covariance(covariance);
The class TwoPointCorrelation2D_polar.
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 polar coordinates
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 polar coordinates, with Poisson errors
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 measureJackknife(const std::string dir_output_pairs, const std::vector< std::string > dir_input_pairs={}, const std::string dir_output_resample="NULL", 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 polar coordinates, estimating the covariance with Ja...
void measureBootstrap(const int nMocks, const std::string dir_output_pairs, const std::vector< std::string > dir_input_pairs={}, const std::string dir_output_resample="NULL", 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 polar coordinates, estimating the covariance with Bo...
void set_parameters(const BinType binType_rad, const double rMin, const double rMax, const int nbins_rad, const double shift_rad, const BinType binType_mu, const double muMin, const double muMax, const int nbins_mu, const double shift_mu, const CoordinateUnits angularUnits=CoordinateUnits::_radians_, std::function< double(double)> angularWeight=nullptr, const bool compute_extra_info=false)
set the binning parameters
void read(const std::string dir, const std::string file) override
read the 2D two-point correlation function
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