43 using namespace catalogue;
44 using namespace chainmesh;
46 using namespace pairs;
47 using namespace measure;
48 using namespace twopt;
56 if (!compute_extra_info)
57 m_dd = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_comoving_log_, PairInfo::_standard_, rMin, rMax, nbins, shift, angularUnits, angularWeight))
58 : move(Pair::Create(PairType::_comoving_lin_, PairInfo::_standard_, rMin, rMax, nbins, shift, angularUnits, angularWeight));
60 m_dd = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_comoving_log_, PairInfo::_extra_, rMin, rMax, nbins, shift, angularUnits, angularWeight))
61 : move(Pair::Create(PairType::_comoving_lin_, PairInfo::_extra_, rMin, rMax, nbins, shift, angularUnits, angularWeight));
63 m_rr = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_comoving_log_, PairInfo::_standard_, rMin, rMax, nbins, shift, angularUnits))
64 : move(Pair::Create(PairType::_comoving_lin_, PairInfo::_standard_, rMin, rMax, nbins, shift, angularUnits));
66 m_dr = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_comoving_log_, PairInfo::_standard_, rMin, rMax, nbins, shift, angularUnits))
67 : move(Pair::Create(PairType::_comoving_lin_, PairInfo::_standard_, rMin, rMax, nbins, shift, angularUnits));
76 if (!compute_extra_info)
77 m_dd = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_comoving_log_, PairInfo::_standard_, rMin, rMax, binSize, shift, angularUnits, angularWeight))
78 : move(Pair::Create(PairType::_comoving_lin_, PairInfo::_standard_, rMin, rMax, binSize, shift, angularUnits, angularWeight));
80 m_dd = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_comoving_log_, PairInfo::_extra_, rMin, rMax, binSize, shift, angularUnits, angularWeight))
81 : move(Pair::Create(PairType::_comoving_lin_, PairInfo::_extra_, rMin, rMax, binSize, shift, angularUnits, angularWeight));
83 m_rr = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_comoving_log_, PairInfo::_standard_, rMin, rMax, binSize, shift, angularUnits))
84 : move(Pair::Create(PairType::_comoving_lin_, PairInfo::_standard_, rMin, rMax, binSize, shift, angularUnits));
86 m_dr = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_comoving_log_, PairInfo::_standard_, rMin, rMax, binSize, shift, angularUnits))
87 : move(Pair::Create(PairType::_comoving_lin_, PairInfo::_standard_, rMin, rMax, binSize, shift, angularUnits));
96 m_dataset->read(dir+file);
105 vector<double> xx = m_dataset->xx();
109 string header =
"[1] separation at the bin centre # [2] spherically averagerded two-point correlation function # [3] error";
110 if (m_compute_extra_info) header +=
" # [4] mean separation # [5] standard deviation of the separation distribution # [6] mean redshift # [7] standard deviation of the redshift distribution";
112 m_dataset->write(dir, file, header, 5, rank);
119 void cbl::measure::twopt::TwoPointCorrelation1D_monopole::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)
122 case (ErrorType::_Poisson_) :
123 measurePoisson(dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
125 case (ErrorType::_Jackknife_) :
126 measureJackknife(dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact);
128 case (ErrorType::_Bootstrap_) :
129 measureBootstrap(nMocks, dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact, seed);
131 case (ErrorType::_JackknifeTest_) :
132 measureJackknifeTest(dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact);
136 ErrorCBL(
"unknown type of error",
"measure",
"TwoPointCorrelation1D_monopole.cpp");
148 count_allPairs(m_twoPType, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
153 if (estimator==Estimator::_natural_)
154 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
155 else if (estimator==Estimator::_LandySzalay_)
156 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
158 ErrorCBL(
"the chosen estimator is not implemented!",
"measurePoisson",
"TwoPointCorrelation1D_monopole.cpp");
169 string mkdir =
"mkdir -p "+dir_output_resample;
170 if (system(mkdir.c_str())) {}
173 vector<long> region_list = m_data->region_list();
174 size_t nRegions = region_list.size();
176 vector<vector<double> > xi_SubSamples, covariance;
178 vector<shared_ptr<Pair> > dd_regions, rr_regions, dr_regions;
180 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);
182 vector<shared_ptr<Data> > data_SS = (estimator==Estimator::_natural_) ? XiJackknife(dd_regions, rr_regions) : XiJackknife(dd_regions, rr_regions, dr_regions);
184 for (
size_t i=0; i<nRegions; i++) {
188 string header =
"[1] separation at the bin centre # [2] spherically averagerded two-point correlation function # [3] error";
189 if (m_compute_extra_info) header +=
" # [4] mean separation # [5] standard deviation of the separation distribution # [6] mean redshift # [7] standard deviation of the redshift distribution";
190 data_SS[i]->write(dir_output_resample, file, header, 10, 0);
193 vector<double> dd; data_SS[i]->get_data(dd);
195 xi_SubSamples.push_back(dd);
200 if (estimator==Estimator::_natural_)
201 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
202 else if (estimator==Estimator::_LandySzalay_)
203 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
205 ErrorCBL(
"the chosen estimator is not implemented!",
"measureJackknife",
"TwoPointCorrelation1D_monopole.cpp");
207 m_dataset->set_covariance(covariance);
218 string mkdir =
"mkdir -p "+dir_output_resample;
219 if (system(mkdir.c_str())) {}
222 vector<long> region_list = m_data->region_list();
223 size_t nRegions = region_list.size();
224 vector<double> weights(nRegions, 1.);
226 vector<vector<double> > xi_SubSamples, covariance;
228 count_allPairs_region_test(m_twoPType, weights, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
230 vector<shared_ptr<Data> > data_SS = (estimator==Estimator::_natural_) ? XiJackknifeTest(m_dd_res, m_rr_res) : XiJackknifeTest(m_dd_res, m_rr_res, m_dr_res);
232 for (
size_t i=0; i<nRegions; i++) {
236 string header =
"[1] separation at the bin centre # [2] spherically averagerded two-point correlation function # [3] error";
237 if (m_compute_extra_info) header +=
" # [4] mean separation # [5] standard deviation of the separation distribution # [6] mean redshift # [7] standard deviation of the redshift distribution";
238 data_SS[i]->write(dir_output_resample, file, header, 10, 0);
241 vector<double> dd; data_SS[i]->get_data(dd);
243 xi_SubSamples.push_back(dd);
248 if (estimator==Estimator::_natural_)
249 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
250 else if (estimator==Estimator::_LandySzalay_)
251 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
253 ErrorCBL(
"the chosen estimator is not implemented!",
"measureJackknifeTest",
"TwoPointCorrelation1D_monopole.cpp");
255 m_dataset->set_covariance(covariance);
262 void cbl::measure::twopt::TwoPointCorrelation1D_monopole::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)
265 ErrorCBL(
"number of mocks must be >0",
"measureBootstrap",
"TwoPointCorrelation1D_monopole.cpp");
268 string mkdir =
"mkdir -p "+dir_output_resample;
269 if (system(mkdir.c_str())) {}
272 vector<vector<double> > xi_SubSamples, covariance;
274 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_) ? XiBootstrap(nMocks, dd_regions, rr_regions, seed) : XiBootstrap(nMocks, dd_regions, rr_regions, dr_regions, seed);
280 for (
int i=0; i<nMocks; i++) {
284 string header =
"[1] separation at the bin centre # [2] spherically averagerded two-point correlation function # [3] error";
285 if (m_compute_extra_info) header +=
" # [4] mean separation # [5] standard deviation of the separation distribution # [6] mean redshift # [7] standard deviation of the redshift distribution";
286 data_SS[i]->write(dir_output_resample, file, header, 10, 0);
289 vector<double> dd; data_SS[i]->get_data(dd);
291 xi_SubSamples.push_back(dd);
296 if (estimator==Estimator::_natural_)
297 m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
298 else if (estimator==Estimator::_LandySzalay_)
299 m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
301 ErrorCBL(
"the chosen estimator is not implemented!",
"measureBootstrap",
"TwoPointCorrelation1D_monopole.cpp");
303 m_dataset->set_covariance(covariance);
The class TwoPointCorrelation1D_monopole.
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 monopole of the two-point correlation function with Poisson errors
void measureJackknifeTest(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 monopole of the two-point correlation function estimating the covariance with Jackknife r...
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 monopole of the two-point correlation function estimating the covariance with Jackknife r...
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 monopole of the two-point correlation function (with the direct estimator)
void write(const std::string dir=par::defaultString, const std::string file=par::defaultString, const int rank=0) const override
write the monopole of the two-point correlation function
void set_parameters(const BinType binType, const double rMin, const double rMax, const int nbins, const double shift, const CoordinateUnits angularUnits=CoordinateUnits::_radians_, std::function< double(double)> angularWeight=nullptr, const bool compute_extra_info=false)
set the binning parameters
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 monopole of the two-point correlation function estimating the covariance with Bootstrap r...
void read(const std::string dir, const std::string file) override
read the monopole of the two-point correlation
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