46 using namespace catalogue;
47 using namespace measure;
48 using namespace triplets;
49 using namespace threept;
61 ThreePointCorrelation_comoving_connected::measure(dir_output_triplets, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
62 m_scale = ThreePointCorrelation_comoving_connected::m_scale;
69 double rMin = m_ddd->r12();
70 double rMax = m_ddd->r12()+m_ddd->r13()+m_ddd->r13_binSize();
71 double binSize = 0.05;
74 TwoP.measure(ErrorType::_Poisson_, dir_output_triplets, {},
par::defaultString, 0, 1, 1, 1, tcount,
cbl::measure::twopt::Estimator::_LandySzalay_, fact);
76 vector<double> log_r(TwoP.dd()->nbins()), log_xi(TwoP.dd()->nbins());
77 for (
int i=0; i<TwoP.dd()->nbins(); i++) {
78 log_r[i] = log10(TwoP.xx()[i]);
79 log_xi[i] = log10(1.+TwoP.xi1D()[i]);
82 vector<double> values_interp(m_ddd->nbins()+2, 0.), xi_real_lin(m_ddd->nbins()+2, 0.);
84 values_interp[0] = log10(m_ddd->r12());
85 values_interp[1] = log10(m_ddd->r13());
87 vector<double> theta(m_ddd->nbins(),0.);
89 for (
int i=0; i<m_ddd->nbins(); i++) {
90 theta[i]=(i+0.5)*m_ddd->binSize();
91 double tmp_value = sqrt(m_ddd->r12()*m_ddd->r12()+m_ddd->r13()*m_ddd->r13()-2.*m_ddd->r12()*m_ddd->r13()*cos(theta[i]));
93 values_interp[i+2] = log10(tmp_value);
96 string file_2pt = dir_output_2pt+
"2ptCorrelation_3pt.dat";
97 ofstream fout(file_2pt.c_str());
checkIO(fout, file_2pt);
99 for (
size_t i=0; i<values_interp.size(); i++) {
100 xi_real_lin[i] = pow(10.,
interpolated(values_interp[i], log_r, log_xi,
"Linear"))-1.;
102 fout << pow(10.,values_interp[i]) <<
" " << xi_real_lin[i] << endl;
105 fout.clear(); fout.close();
coutCBL <<
"I wrote the file "<<file_2pt<<endl;
110 m_QQ.resize(m_ddd->nbins()); m_error.resize(m_ddd->nbins());
112 for (
int i=0; i<m_ddd->nbins(); i++) {
113 m_QQ[i] = m_zeta[i]/((xi_real_lin[0]*xi_real_lin[1])+(xi_real_lin[0]*xi_real_lin[i+2])+(xi_real_lin[1]*xi_real_lin[i+2]));
117 m_dataset = move(unique_ptr<data::Data1D>(
new data::Data1D(theta, m_QQ, m_error)));
124 void cbl::measure::threept::ThreePointCorrelation_comoving_reduced::measure (
const std::vector<std::vector<double>> weight,
const bool doJK,
const std::string dir_output_triplets,
const std::string dir_output_2pt,
const std::vector<std::string> dir_input_triplets,
const bool count_ddd,
const bool count_rrr,
const bool count_ddr,
const bool count_drr,
const bool tcount,
const double fact,
const int seed)
130 ThreePointCorrelation_comoving_connected::measure(weight, doJK, dir_output_triplets, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
131 m_scale = ThreePointCorrelation_comoving_connected::m_scale;
138 double rMin = m_ddd->r12();
139 double rMax = m_ddd->r12()+m_ddd->r13()+m_ddd->r13_binSize();
140 double binSize = 0.05;
143 TwoP.measure(ErrorType::_Poisson_,dir_output_triplets, {},
par::defaultString, 0, 1, 1, 1, tcount,
cbl::measure::twopt::Estimator::_LandySzalay_, fact);
145 vector<double> log_r(TwoP.dd()->nbins()), log_xi(TwoP.dd()->nbins());
146 for (
int i=0; i<TwoP.dd()->nbins(); i++) {
147 log_r[i] = log10(TwoP.xx()[i]);
148 log_xi[i] = log10(1.+TwoP.xi1D()[i]);
151 vector<double> values_interp(m_ddd->nbins()+2, 0.), xi_real_lin(m_ddd->nbins()+2, 0.);
153 values_interp[0] = log10(m_ddd->r12());
154 values_interp[1] = log10(m_ddd->r13());
156 vector<double> theta(m_ddd->nbins(),0.);
158 for (
int i=0; i<m_ddd->nbins(); i++) {
159 theta[i]=(i+0.5)*m_ddd->binSize();
160 double tmp_value = sqrt(m_ddd->r12()*m_ddd->r12()+m_ddd->r13()*m_ddd->r13()-2.*m_ddd->r12()*m_ddd->r13()*cos(theta[i]));
162 values_interp[i+2] = log10(tmp_value);
165 string file_2pt = dir_output_2pt+
"2ptCorrelation_3pt.dat";
166 ofstream fout(file_2pt.c_str());
checkIO(fout, file_2pt);
168 for (
size_t i=0; i<values_interp.size(); i++) {
169 xi_real_lin[i] = pow(10.,
interpolated(values_interp[i], log_r, log_xi,
"Linear"))-1.;
171 fout << pow(10.,values_interp[i]) <<
" " << xi_real_lin[i] << endl;
174 fout.clear(); fout.close();
coutCBL <<
"I wrote the file "<<file_2pt<<endl;
179 m_QQ.resize(m_ddd->nbins()); m_error.resize(m_ddd->nbins());
181 for (
int i=0; i<m_ddd->nbins(); i++) {
182 m_QQ[i] = m_zeta[i]/((xi_real_lin[0]*xi_real_lin[1])+(xi_real_lin[0]*xi_real_lin[i+2])+(xi_real_lin[1]*xi_real_lin[i+2]));
188 vector<vector<double>> resampling_threept(weight.size(), vector<double>(m_ddd->nbins(), 0));
189 vector<long> region_list = m_data->region_list();
191 vector<double> nData_reg_weighted, nRandom_reg_weighted;
193 const int nRegions = m_data->nRegions();
195 for (
int i=0; i<nRegions; i++) {
196 nData_reg_weighted.push_back(m_data->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
197 nRandom_reg_weighted.push_back(m_random->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
200 for (
size_t i=0; i<weight.size(); i++) {
205 for (
size_t j=0; j<weight[i].size(); j++) {
206 nData += weight[i][j]*nData_reg_weighted[j];
207 nRan += weight[i][j]*nRandom_reg_weighted[j];
210 double norm1 = (double(nData)*double(nData-1)*double(nData-2))/6.;
211 double norm2 = (double(nData)*double(nData-1)*double(nRan))*0.5;
212 double norm3 = (double(nData)*double(nRan)*double(nRan-1))*0.5;
213 double norm4 = (double(nRan)*double(nRan-1)*double(nRan-2))/6.;
215 for (
int j=0; j<m_ddd->nbins(); j++)
216 if (m_ddd_regions[i]->TT1D(j)>0 && m_rrr_regions[i]->TT1D(j)>0) {
217 resampling_threept[i][j] = ((m_ddd_regions[i]->TT1D(j)/norm1)/(m_rrr_regions[i]->TT1D(j)/norm4))-3.*((m_ddr_regions[i]->TT1D(j)/norm2)/(m_rrr_regions[i]->TT1D(j)/norm4))+3.*(((m_drr_regions[i]->TT1D(j)/norm3)/(m_rrr_regions[i]->TT1D(j)/norm4)))-1.;
218 resampling_threept[i][j] /= ((xi_real_lin[0]*xi_real_lin[1])+(xi_real_lin[0]*xi_real_lin[j+2])+(xi_real_lin[1]*xi_real_lin[j+2]));
222 vector<vector<double>> cov_mat;
224 m_dataset = move(unique_ptr<data::Data1D>(
new data::Data1D(m_scale, m_zeta, cov_mat)));
226 m_error = m_dataset->error();
234 void cbl::measure::threept::ThreePointCorrelation_comoving_reduced::measure (
const cbl::measure::ErrorType errorType,
const std::string dir_output_triplets,
const std::string dir_output_2pt,
const std::vector<std::string> dir_input_triplets,
const int nResamplings,
const bool count_ddd,
const bool count_rrr,
const bool count_ddr,
const bool count_drr,
const bool tcount,
const double fact,
const int seed)
241 measure(dir_output_triplets, dir_output_2pt, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
247 const int nRegions = m_data->nRegions();
249 vector<vector<double>> weight(nRegions, vector<double>(nRegions, 1));
250 for (
int i=0; i<nRegions; i++)
253 measure(weight,
true, dir_output_triplets, dir_output_2pt, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
259 const int nRegions = m_data->nRegions();
265 vector<vector<double>> weight(nResamplings, vector<double>(nRegions, 0));
266 for (
int i=0; i<nResamplings; i++)
267 for (
int j=0; j<val*nRegions; j++)
270 measure(weight,
false, dir_output_triplets, dir_output_2pt, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
275 ErrorCBL(
"no such kind of error type!",
"measure",
"hreePointCorrelation_comoving_reduced.cpp");
286 checkDim(m_scale, m_ddd->nbins(),
"scale");
288 string file_out = dir+file;
289 ofstream fout(file_out.c_str());
checkIO(fout, file_out);
292 fout <<
"# scale Q error(work in progress)" << endl;
293 for (
size_t i=0; i<m_scale.size(); i++)
294 fout << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_scale[i]
295 <<
" " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_QQ[i]
296 <<
" " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_error[i] << endl;
300 fout <<
"# scale z error(work in progrss) Q error(work in progress)" << endl;
301 for (
size_t i=0; i<m_scale.size(); i++)
302 fout << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_scale[i]
303 <<
" " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_zeta[i]
304 <<
" " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << ThreePointCorrelation_comoving_connected::m_error[i]
305 <<
" " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_QQ[i]
306 <<
" " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_error[i] << endl;
309 fout.close();
coutCBL << endl <<
"I wrote the file: " << file_out << endl << endl;
318 m_dataset->write_covariance(dir, file);
#define coutCBL
CBL print message.
The class ThreePointCorrelation_comoving_reduced.
The class TwoPointCorrelation1D_monopole.
void measure(const std::string dir_output_triplets, const std::string dir_output_2pt, const std::vector< std::string > dir_input_triplets={}, const bool count_ddd=true, const bool count_rrr=true, const bool count_ddr=true, const bool count_drr=true, const bool tcount=false, const double fact=0.1, const int seed=3213) override
method to measure the three-point correlation function
void write_covariance(const std::string dir, const std::string file) const override
write the measured three-point correlation covariance
void write(const std::string dir, const std::string file, const bool connected) const override
write the monopole of the two-point correlation function
The class TwoPointCorrelation1D_monopole.
static const std::string defaultString
default std::string value
@ _LandySzalay_
Landy&Szalay estimator.
ErrorType
the two-point correlation function error type
@ _Bootstrap_
Bootstrap resampling.
@ _None_
No error computed.
@ _Jackknife_
Jackknife resampling.
The global namespace of the CosmoBolognaLib
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
void checkIO(const std::ifstream &fin, const std::string file="NULL")
check if an input file can be opened
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
double interpolated(const double _xx, const std::vector< double > xx, const std::vector< double > yy, const std::string type)
1D interpolation
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