forked from EiffL/Angpow4CCL
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathangpow.cc
308 lines (246 loc) · 8.72 KB
/
angpow.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
#include <iostream>
#include "Angpow/walltimer.h" //profiling
#include "Angpow/angpow_exceptions.h" //exceptions
#include "Angpow/angpow_numbers.h" //r_8... def
#include "Angpow/angpow_parameters.h" //control parameters
#include "Angpow/angpow_powspec.h" //power spectrum IMPLEMENTATION
#include "Angpow/angpow_integrand.h" //f_ell(k,z) integrand functions
#include "Angpow/angpow_radial.h" //radial window
#include "Angpow/angpow_utils.h" //utility functions
#include "Angpow/angpow_pk2cl.h" //utility class that produce the Cl (auto/cross corr)
#include "Angpow/angpow_clbase.h" //class to manage Clbase structure (ell, val)
#include "Angpow/angpow_ctheta.h" //class function to manage CTheta values
#include "Angpow/angpow_tools.h"
#ifdef _OPENMP
#include<omp.h>
#endif
//------------------------------------------------------
namespace Angpow {
//------------------------------
// Exemple of processing from P(k) to Cl
//------------------------------
void process() {
//Get the pointer to the job processing user parameters
Parameters para = Param::Instance().GetParam();
int Lmax = para.Lmax; //ell in [0, Lmax-1]
tstack_push("Processing....");
//Radial (redshift) selection windows
RadSelectBase* Z1win = 0;
switch(para.wtype1) {
case Parameters::GaussGal:
{
r_8 zmean = para.mean1;
r_8 zsigma = para.width1;
r_8 zmin = zmean - para.n_sigma_cut * zsigma;
r_8 zmax = zmean + para.n_sigma_cut * zsigma;
std::cout << "Galaxy density + Gaussian window Z1: mean/sigma/min/max: "
<< zmean << "/"
<< zsigma << "/"
<< zmin << "/"
<< zmax
<< std::endl;
Z1win = new RadModifiedGaussSelect(zmean, zsigma, zmin, zmax);
break;
}
case Parameters::Gauss:
{
r_8 zmean = para.mean1;
r_8 zsigma = para.width1;
r_8 zmin = zmean - para.n_sigma_cut * zsigma;
r_8 zmax = zmean + para.n_sigma_cut * zsigma;
std::cout << "Gaussian window Z1: mean/sigma/min/max: "
<< zmean << "/"
<< zsigma << "/"
<< zmin << "/"
<< zmax
<< std::endl;
Z1win = new RadGaussSelect(zmean, zsigma, zmin, zmax);
break;
}
case Parameters::Dirac:
{
r_8 zmin = 0; //not used
r_8 zmax = 0; //not used
r_8 zmean = para.mean1;
std::cout << "Dirac window Z1: mean: " << zmean << std::endl;
Z1win = new DiracSelect(zmean, zmin, zmax);
break;
}
case Parameters::TopHat:
{
r_8 zmean = para.mean1;
r_8 zwidth = para.width1; //note 1/2 full width
r_8 smooth_edges = para.smooth_edges;
r_8 nSigmaCut = para.n_sigma_cut;
r_8 zmin = zmean - (1. + nSigmaCut*smooth_edges) * zwidth;
r_8 zmax = zmean + (1. + nSigmaCut*smooth_edges) * zwidth;
std::cout << "TopHas window Z1: mean, width, min, max: " <<
zmean << " " << zwidth << " " << zmin << " " << zmax << std::endl;
Z1win = new RadTopHatSmoothSelect(zmin,zmax,zmean,zwidth,smooth_edges);
break;
}
default:
throw AngpowError("Unknown Select 1 type");
break;
}//sw
if(Z1win == 0) throw AngpowError("process: FATAL ERROR: unknown Z1 selection function");
RadSelectBase* Z2win = 0;
switch(para.wtype2) {
case Parameters::GaussGal:
{
r_8 zmean = para.mean2;
r_8 zsigma = para.width2;
r_8 zmin = zmean - para.n_sigma_cut * zsigma;
r_8 zmax = zmean + para.n_sigma_cut * zsigma;
std::cout << "Galaxy density + Gaussian window Z2: mean/sigma/min/max: "
<< zmean << "/"
<< zsigma << "/"
<< zmin << "/"
<< zmax
<< std::endl;
Z2win = new RadModifiedGaussSelect(zmean, zsigma, zmin, zmax);
break;
}
case Parameters::Gauss:
{
r_8 zmean = para.mean2;
r_8 zsigma = para.width2;
r_8 zmin = zmean - para.n_sigma_cut * zsigma;
r_8 zmax = zmean + para.n_sigma_cut * zsigma;
std::cout << "Gaussian window Z2: mean/sigma/min/max: "
<< zmean << "/"
<< zsigma << "/"
<< zmin << "/"
<< zmax
<< std::endl;
Z2win = new RadGaussSelect(zmean, zsigma, zmin, zmax);
break;
}
case Parameters::Dirac:
{
r_8 zmin = 0; //not used
r_8 zmax = 0; //not used
r_8 zmean = para.mean2;
std::cout << "Dirac window Z2: mean: " << zmean << std::endl;
Z2win = new DiracSelect(zmean, zmin, zmax);
break;
}
case Parameters::TopHat:
{
r_8 zmean = para.mean2;
r_8 zwidth = para.width2; //note 1/2 full width
r_8 smooth_edges = para.smooth_edges;
r_8 nSigmaCut = para.n_sigma_cut;
r_8 zmin = zmean - (1. + nSigmaCut*smooth_edges) * zwidth;
r_8 zmax = zmean + (1. + nSigmaCut*smooth_edges) * zwidth;
std::cout << "TopHas window Z2: mean, width, min, max: " <<
zmean << " " << zwidth << " " << zmin << " " << zmax << std::endl;
Z2win = new RadTopHatSmoothSelect(zmin,zmax,zmean,zwidth,smooth_edges);
}
break;
default:
throw AngpowError("Unknown Select 2 type");
break;
}//sw
if(Z2win == 0) throw AngpowError("process: FATAL ERROR: unknown Z2 selection function");
//The cosmological distance tool
CosmoCoord cosmo(para.cosmo_zmin, para.cosmo_zmax, para.cosmo_npts, para.cosmo_precision);
//Usage of an external file of P(k) and use of an internal Growth function
//If k unit is h/Mpc and P(k) unit is [Mpc/h]^3 then set use_h_rescaling = true
//If the P(k) is defined at zref=0 then no growth rescaling is needed.
//
std::string pwName = para.power_spectrum_input_dir + para.power_spectrum_input_file;
r_8 pw_kmin = para.pw_kmin;
r_8 pw_kmax = para.pw_kmax;
r_8 zref = 0.;
bool use_h_rescaling = true;
bool use_growth_rescaling= false;
PowerSpecFile pws(cosmo,pwName,zref,pw_kmin,pw_kmax,use_h_rescaling,use_growth_rescaling);
// JN 20/04/2017
r_8 bias = para.bias;
bool include_rsd = para.include_rsd;
Integrand int1(&pws, &cosmo, bias, include_rsd);
Integrand int2(&pws, &cosmo, bias, include_rsd);
//Initialize the Cl with parameters to select the ell set which is interpolated after the processing
Clbase clout(Lmax,para.linearStep, para.logStep);
//Main class
Pk2Cl pk2cl; //Default: the user parameters are used in the Constructor
pk2cl.PrintParam();
pk2cl.Compute(int1, int2, cosmo, Z1win, Z2win, Lmax, clout);
//pk2cl.Compute(pws, cosmo, Z1win, Z2win, Lmax, clout);
tstack_pop("Processing....");
tstack_report("Processing....");
{//save the Cls
std::fstream ofs;
std::string outName = para.output_dir + para.common_file_tag + "cl.txt";
ofs.open(outName, std::fstream::out);
for(int index_l=0;index_l<clout.Size();index_l++){
ofs << setprecision(20) << clout[index_l].first << " " << clout[index_l].second << endl;
}
ofs.close();
}
{//save ctheta
CTheta ct(clout,para.apod);
std::fstream ofs;
std::string outName = para.output_dir + para.common_file_tag + "ctheta.txt";
//define theta values
const int Npts=100;
const double theta_max=para.theta_max*M_PI/180;
double step=theta_max/(Npts-1);
ofs.open(outName, std::fstream::out);
for (size_t i=0;i<Npts;i++){
double t=i*step;
ofs << setprecision(20) << t << " " << ct(t) << endl;
}
ofs.close();
outName = para.output_dir + para.common_file_tag + "apod_cl.txt";
ct.WriteApodCls(outName);
}
//clear
pws.ExplicitDestroy();
if(Z1win) delete Z1win; Z1win = 0;
if(Z2win) delete Z2win; Z2win = 0;
std::cout << "End process......" << std::endl;
}//process
}//namespace
//------------------------------------------------------
int main(int narg, char* argv[]) {
using namespace Angpow;
//Default return code
int rc = 0;
try {
//Default Initialization
if(narg<2) {
throw AngpowError("usage: angpow <init-file>");
}
std::string fileIni = argv[1];
#ifdef _OPENMP
cout << "using OpenMP with max threads=" << omp_get_max_threads() << endl;
cout << "dynamic adjustment is: "<< omp_get_dynamic()<<endl;
#endif
Param::Instance().ReadParam(fileIni);
{
Parameters para = Param::Instance().GetParam();
std::fstream ofs;
std::string outName = para.output_dir + para.common_file_tag + "used-param.txt";
ofs.open(outName, std::fstream::out);
Param::Instance().WriteParam(ofs);
}
//Go to processing
process();
}
catch (std::exception& sex) {
cerr << "\n angpow.cc std::exception :" << (string)typeid(sex).name()
<< "\n msg= " << sex.what() << endl;
rc = 78;
}
catch ( string str ) {
cerr << "angpow.cc Exception raised: " << str << endl;
}
catch (...) {
cerr << " angpow.cc catched unknown (...) exception " << endl;
rc = 79;
}
cout << ">>>> angpow.cc ------- END ----------- RC=" << rc << endl;
return rc;
}//main