44 using namespace catalogue;
45 using namespace chainmesh;
46 using namespace pairs;
47 using namespace measure;
48 using namespace twopt;
58 vector<double> weightTOT(dd2D->nbins_D1(), 0.), scale_mean(dd2D->nbins_D1(), 0.), scale_sigma(dd2D->nbins_D1(), 0.), z_mean(dd2D->nbins_D1(), 0.), z_sigma(dd2D->nbins_D1(), 0.);
60 double fact_err, fact_scale, fact_z;
62 for (
int i=0; i<dd2D->nbins_D1(); ++i) {
64 for (
int j=0; j<dd2D->nbins_D2(); ++j)
65 weightTOT[i] += dd2D->PP2D_weighted(i, j);
67 for (
int j=0; j<dd2D->nbins_D2(); ++j) {
68 scale_mean[i] += dd2D->scale_D1_mean(i, j)*dd2D->PP2D_weighted(i, j)/weightTOT[i];
69 z_mean[i] += dd2D->z_mean(i, j)*dd2D->PP2D_weighted(i, j)/weightTOT[i];
72 scale_sigma[i] = pow(dd2D->scale_D1_sigma(i, 0), 2)*dd2D->PP2D_weighted(i, 0);
73 z_sigma[i] = pow(dd2D->z_sigma(i, 0), 2)*dd2D->PP2D_weighted(i, 0);
75 for (
int j=1; j<dd2D->nbins_D2(); ++j) {
76 if (dd2D->PP2D_weighted(i, j)>0) {
77 fact_err = dd2D->PP2D_weighted(i, j)*dd2D->PP2D_weighted(i, j-1)/(dd2D->PP2D_weighted(i, j)+dd2D->PP2D_weighted(i, j-1));
78 fact_scale = pow(dd2D->scale_D1_mean(i, j)-dd2D->scale_D1_mean(i, j-1), 2)*fact_err;
79 fact_z = pow(dd2D->z_mean(i, j)-dd2D->z_mean(i, j-1), 2)*fact_err;
80 scale_sigma[i] += pow(dd2D->scale_D1_sigma(i, j), 2)*dd2D->PP2D_weighted(i, j)+fact_scale;
81 z_sigma[i] += pow(dd2D->z_sigma(i, j),2)*weightTOT[i]+fact_z;
86 vector<vector<double>> extra(4);
88 for (
int i=0; i<dd2D->nbins_D1(); ++i) {
89 extra[0].push_back(scale_mean[i]);
90 extra[1].push_back(sqrt(scale_sigma[i]/weightTOT[i]));
91 extra[2].push_back(z_mean[i]);
92 extra[3].push_back(sqrt(z_sigma[i]/weightTOT[i]));
95 return move(unique_ptr<data::Data1D_extra>(
new data::Data1D_extra(rp, ww, error, extra)));
104 vector<double> ww, error;
106 ww.resize(0); ww.resize(rp.size(), 0.);
107 error.resize(0); error.resize(rp.size(), 0.);
109 const double binSize = 1./m_dd->binSize_inv_D2();
111 const int pim =
nint((m_piMax_integral-
Min(
pi))/binSize);
113 for (
size_t i=0; i<rp.size(); i++) {
117 for (
size_t j=0; j<pim*1.000001; j++) {
118 if (j<xi[i].size() and j<error_xi[i].size()) {
119 ww[i] = ww[i]+2.*binSize*xi[i][j];
120 if (ww[i]>-1.) error[i] += pow(2.*binSize*error_xi[i][j], 2);
125 for_each( error.begin(), error.end(), [] (
double &vv) { vv = sqrt(vv);} );
127 return (!m_compute_extra_info) ? move(unique_ptr<data::Data1D>(
new data::Data1D(rp, ww, error))) : data_with_extra_info(rp, ww, error);
137 m_dataset->read(dir+file);
146 vector<double> xx = m_dataset->xx();
148 checkDim(xx, m_dd->nbins_D1(),
"rp");
150 string header =
"[1] perpendicular separation at the bin centre # [2] projected two-point correlation function # [3] error";
151 if (m_compute_extra_info) header +=
" # [4] mean perpendicular separation # [5] standard deviation of the distribution of perpendicular separations # [6] mean redshift # [7] standard deviation of the redshift distribution";
153 m_dataset->write(dir, file, header, 5, rank);
160 void cbl::measure::twopt::TwoPointCorrelation_projected::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)
163 case (ErrorType::_Poisson_) :
164 measurePoisson(dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
166 case (ErrorType::_Jackknife_) :
167 measureJackknife(dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact);
169 case (ErrorType::_Bootstrap_) :
170 measureBootstrap(nMocks, dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact, seed);
173 ErrorCBL(
"unknown type of error!",
"measure",
"TwoPointCorrelation_projected.cpp");
185 TwoPointCorrelation2D_cartesian::measurePoisson(dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
190 m_dataset = Projected(TwoPointCorrelation2D_cartesian::xx(), TwoPointCorrelation2D_cartesian::yy(), TwoPointCorrelation2D_cartesian::xi2D(), TwoPointCorrelation2D_cartesian::error2D());
201 string mkdir =
"mkdir -p "+dir_output_resample;
202 if (system(mkdir.c_str())) {}
205 vector<shared_ptr<data::Data>> data;
206 vector<shared_ptr<pairs::Pair>> dd_regions, rr_regions, dr_regions;
207 count_allPairs_region(dd_regions, rr_regions, dr_regions, TwoPType::_2D_Cartesian_, dir_output_pairs,dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
209 auto data_cart = (estimator==Estimator::_natural_) ? correlation_NaturalEstimator(m_dd, m_rr) : correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
211 if (estimator==Estimator::_natural_)
212 data = XiJackknife(dd_regions, rr_regions);
213 else if (estimator==Estimator::_LandySzalay_)
214 data = XiJackknife(dd_regions, rr_regions, dr_regions);
216 ErrorCBL(
"the chosen estimator is not implemented!",
"measureJackknife",
"TwoPointCorrelation_projected.cpp");
218 vector<vector<double>> ww, covariance;
219 for (
size_t i=0; i<data.size(); i++) {
220 ww.push_back(data[i]->data());
222 string file =
"xi_projected_Jackkknife_"+
conv(i,
par::fINT)+
".dat";
223 string header =
"[1] perpendicular separation at the bin centre # [2] projected two-point correlation function # [3] error";
224 if (m_compute_extra_info) header +=
" # [4] mean perpendicular separation # [5] standard deviation of the distribution of perpendicular separations # [6] mean redshift # [7] standard deviation of the redshift distribution";
225 data[i]->write(dir_output_resample, file, header, 10, 0);
231 vector<double> xx_cart = data_cart->xx(), yy_cart = data_cart->yy();
232 vector<vector<double> > dd_cart, error_cart;
233 data_cart->get_data(dd_cart); data_cart->get_error(error_cart);
235 m_dataset = Projected(xx_cart, yy_cart, dd_cart, error_cart);
236 m_dataset->set_covariance(covariance);
244 void cbl::measure::twopt::TwoPointCorrelation_projected::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)
247 ErrorCBL(
"the number of mocks must be >0!",
"measureBootstrap",
"TwoPointCorrelation1D_monopole.cpp");
250 string mkdir =
"mkdir -p "+dir_output_resample;
251 if (system(mkdir.c_str())) {}
254 vector<shared_ptr<data::Data>> data;
255 vector<shared_ptr<pairs::Pair>> dd_regions, rr_regions, dr_regions;
256 count_allPairs_region(dd_regions, rr_regions, dr_regions, TwoPType::_2D_Cartesian_, dir_output_pairs,dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
258 auto data_cart = (estimator==Estimator::_natural_) ? correlation_NaturalEstimator(m_dd, m_rr) : correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
260 if (estimator==Estimator::_natural_)
261 data = XiBootstrap(nMocks, dd_regions, rr_regions, seed);
262 else if (estimator==Estimator::_LandySzalay_)
263 data = XiBootstrap(nMocks, dd_regions, rr_regions, dr_regions, seed);
265 ErrorCBL(
"the chosen estimator is not implemented!",
"measureBootstrap",
"TwoPointCorrelation_projected.cpp");
267 vector<vector<double>> ww, covariance;
268 for (
size_t i=0; i<data.size(); i++) {
269 ww.push_back(data[i]->data());
271 string file =
"xi_projected_Bootstrap_"+
conv(i,
par::fINT)+
".dat";
272 string header =
"[1] perpendicular separation at the bin centre # [2] projected two-point correlation function # [3] error";
273 if (m_compute_extra_info) header +=
" # [4] mean perpendicular separation # [5] standard deviation of the distribution of perpendicular separations # [6] mean redshift # [7] standard deviation of the redshift distribution";
274 data[i]->write(dir_output_resample, file, header, 10, 0);
280 vector<double> xx_cart = data_cart->xx(), yy_cart = data_cart->yy();
281 vector<vector<double> > dd_cart, error_cart;
282 data_cart->get_data(dd_cart); data_cart->get_error(error_cart);
284 m_dataset = Projected(xx_cart, yy_cart, dd_cart, error_cart);
285 m_dataset->set_covariance(covariance);
294 vector<shared_ptr<data::Data>> data;
296 auto data2d = TwoPointCorrelation2D_cartesian::XiJackknife(dd, rr);
298 for (
size_t i=0; i<data2d.size(); i++) {
299 vector<double> xx_cart = data2d[i]->xx(), yy_cart = data2d[i]->yy();
300 vector<vector<double>> dd_cart, error_cart;
301 data2d[i]->get_data(dd_cart); data2d[i]->get_error(error_cart);
302 data.push_back(move(Projected(xx_cart, yy_cart, dd_cart, error_cart)));
314 vector<shared_ptr<data::Data>> data;
316 auto data2d = TwoPointCorrelation2D_cartesian::XiJackknife(dd, rr, dr);
318 for (
size_t i=0; i<data2d.size(); i++) {
319 vector<double> xx_cart = data2d[i]->xx(), yy_cart = data2d[i]->yy();
320 vector<vector<double>> dd_cart, error_cart;
321 data2d[i]->get_data(dd_cart); data2d[i]->get_error(error_cart);
322 data.push_back(move(Projected(xx_cart, yy_cart, dd_cart, error_cart)));
334 vector<shared_ptr<data::Data>> data;
336 auto data2d = TwoPointCorrelation2D_cartesian::XiBootstrap(nMocks, dd, rr, seed);
338 for (
size_t i=0; i<data2d.size(); i++) {
339 vector<double> xx_cart = data2d[i]->xx(), yy_cart = data2d[i]->yy();
340 vector<vector<double>> dd_cart, error_cart;
341 data2d[i]->get_data(dd_cart); data2d[i]->get_error(error_cart);
342 data.push_back(move(Projected(xx_cart, yy_cart, dd_cart, error_cart)));
354 vector<shared_ptr<data::Data>> data;
356 auto data2d = TwoPointCorrelation2D_cartesian::XiBootstrap(nMocks, dd, rr, dr, seed);
358 for (
size_t i=0; i<data2d.size(); i++) {
359 vector<double> xx_cart = data2d[i]->xx(), yy_cart = data2d[i]->yy();
360 vector<vector<double>> dd_cart, error_cart;
361 data2d[i]->get_data(dd_cart); data2d[i]->get_error(error_cart);
362 data.push_back(move(Projected(xx_cart, yy_cart, dd_cart, error_cart)));
374 m_dataset->set_covariance(dir+file);
383 m_dataset->write_covariance(dir, file);
392 vector<vector<double>> Xi;
394 for (
size_t i=0; i<xi.size(); i++)
395 Xi.push_back(xi[i]->data());
397 vector<vector<double>> cov_mat;
400 m_dataset->set_covariance(cov_mat);
409 vector<double> rad, mean;
410 vector<vector<double>> cov_mat;
413 m_dataset->set_covariance(cov_mat);
The class TwoPointCorrelation_projected.
void write_covariance(const std::string dir, const std::string file) const override
write the measured two-point correlation
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 projected two-point correlation function with Poisson errors
std::shared_ptr< data::Data > Projected(const std::vector< double > rp, const std::vector< double > pi, const std::vector< std::vector< double >> xi, const std::vector< std::vector< double >> error_xi) override
measure projected correlation function
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 projected two-point correlation function
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 projected two-point correlation function estimating the covariance with Jackknife resampl...
std::shared_ptr< data::Data > data_with_extra_info(const std::vector< double > rp, const std::vector< double > ww, const std::vector< double > error) const
return a data object with extra info
std::vector< std::shared_ptr< data::Data > > XiBootstrap(const int nMocks, const std::vector< std::shared_ptr< pairs::Pair >> dd, const std::vector< std::shared_ptr< pairs::Pair >> rr, const int seed=3213) override
measure the bootstrap resampling of the two-point correlation function, ξ(r)
void read(const std::string dir, const std::string file) override
read the projected two-point correlation function
void write(const std::string dir=par::defaultString, const std::string file=par::defaultString, const int rank=0) const override
write the projected two-point correlation function
void compute_covariance(const std::vector< std::shared_ptr< data::Data >> xi, const bool JK) override
compute the covariance matrix
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=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 projected two-point correlation function estimating the covariance with Bootstrap resampl...
void read_covariance(const std::string dir, const std::string file) override
read the measured covariance matrix
std::vector< std::shared_ptr< data::Data > > XiJackknife(const std::vector< std::shared_ptr< pairs::Pair >> dd, const std::vector< std::shared_ptr< pairs::Pair >> rr) override
measure the jackknife resampling of the two-point correlation function, ξ(r)
std::shared_ptr< pairs::Pair > m_dd
number of data-data pairs
static const char fINT[]
conversion symbol for: int -> std::string
static const std::string defaultString
default std::string value
static const double pi
: the ratio of a circle's circumference to its diameter
Estimator
the two-point correlation estimator
ErrorType
the two-point correlation function error type
The global namespace of the CosmoBolognaLib
T Min(const std::vector< T > vect)
minimum element of a std::vector
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
int nint(const T val)
the nearest integer