source: Sophya/trunk/Cosmo/SimLSS/cmvginit3d.cc@ 3805

Last change on this file since 3805 was 3805, checked in by cmv, 15 years ago

Refonte totale de la structure des classes InitialSpectrum, TransfertFunction, GrowthFactor, PkSpectrum et classes tabulate pour lecture des fichiers CAMB, cmv 24/07/2010
File size: 22.0 KB
Line 
1#include "sopnamsp.h"
2#include "machdefs.h"
3#include <iostream>
4#include <stdlib.h>
5#include <stdio.h>
6#include <string.h>
7#include <math.h>
8#include <unistd.h>
9
10#include "sophyainit.h"
11#include "timing.h"
12#include "dvlist.h"
13#include "ntuple.h"
14#include "matharr.h"
15#include "randfmt.h"
16//#include "randr48.h"
17#include "srandgen.h"
18
19#include "constcosmo.h"
20#include "geneutils.h"
21#include "genefluct3d.h"
22
23void usage(void);
24void usage(void)
25{
26 cout<<"cmvginit3d [...options...]"<<endl
27 <<" -a : auto init random seed (needed for multiple simul)"<<endl
28 <<" -G typevol: compute Pk(z=0) and apply growth factor in real space"<<endl
29 <<" typevol=1 evolved with distance / observateur (def)"<<endl
30 <<" typevol=2 evolved with distance to middle of Z planes"<<endl
31 <<" else : no evol, spectrum Pk(z=z_median) for all cube (def)"<<endl
32 <<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl
33 <<" -x nx,dx : size along x axis (npix,Mpc)"<<endl
34 <<" -y ny,dy : size along y axis (npix,Mpc)"<<endl
35 <<" if ny or dy <=0 take same value as for x"<<endl
36 <<" -z nz,dz : size along z axis (redshift axis, npix,Mpc)"<<endl
37 <<" -Z zref : redshift for the center of the simulation cube"<<endl
38 <<" -1 : compute 1D spectrum (default: no)"<<endl
39 <<" -2 : compute 2D spectrum (default: no)"<<endl
40 <<" -C : go back to cube in FT space and compute P(k) for check (def: do nothing)"<<endl
41 <<" -8 sigmaR,R : normalisation du spectre de puissance, R en Mpc"<<endl
42 <<" (default sigmaR=1, R=8/h100 Mpc)"<<endl
43 <<" -v temp_file.ppf: generate los velocity cube"<<endl
44 <<" temporary save cube in temp_file.ppf"<<endl
45 <<" -n nent : fill control ntuple with nent entries (def: do not fill)"<<endl
46 <<" -O a,b : tell what you want to write"<<endl
47 <<" a : write generated fourier cube (_k0)"<<endl
48 <<" b : write real space cube dRho/Rho at z (_r0)"<<endl
49 <<" a,b= 0 no write, 1 ppf write, 2 fits write, 3 ppf+fits write"<<endl
50 <<" -S : write cube slices in PPF format"<<endl
51 <<" -o out_base_name : base string for output file name (def: cmvginit3d)"<<endl
52 <<" -V : compute variance from real space (for check, default: no)"<<endl
53 <<" -T nth : nombre de threads (si compil multi-thread, default: 0)"<<endl
54 <<endl;
55}
56
57int main(int narg,char *arg[])
58{
59 SophyaInit();
60 InitTim();
61
62 //-----------------------------------------------------------------
63 // *** Survey definition
64 long nx=360, ny=-1, nz=64; double dx=1., dy=-1., dz=-1.;
65 //long nx=1000, ny=-1, nz=128; double dx=3., dy=-1., dz=6.;
66 //long nx=1200, ny=-1, nz=128; double dx=1., dy=-1., dz=3;
67
68 // *** Cosmography definition (WMAP)
69 unsigned short flat = 0;
70 double ob0 = 0.0444356;
71 double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
72 double zref = 0.5;
73 double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
74
75 // *** Spectrum and variance definition
76 double ns = 1., as = 1.;
77 double R=8./h100, Rg=R/sqrt(5.);
78 double sigmaR = 1.;
79
80 double kmin=1e-5,kmax=1000.;
81 int npt = 10000;
82 double lkmin=log10(kmin), lkmax=log10(kmax);
83 double eps=1.e-3;
84
85 // *** type de generation
86 int use_growth_factor = 0;
87 unsigned short nthread=0;
88 int filter_by_pixel = 1;
89
90 // *** What to do
91 bool comptransveloc = false;
92 string temporay_file = "cmvginit3d_tmp.ppf";
93 // compute: 1=1D spectra, 2=2D spectra,
94 // 4=recompute spectra after real space generation
95 unsigned short compdspec = 0;
96 long ntnent = 0; // fill control NTuple
97 bool wslice = false;
98 bool compvarreal = false;
99 unsigned short whattowrt[2] = {1,1};
100 string rootnameout = "cmvobserv3d";
101
102 // --- Decodage arguments
103 if(narg>0) {
104 cout<<"\n--- Arguments: "<<endl;
105 for(int i=0;i<narg;i++) cout<<arg[i]<<" ";
106 cout<<endl;
107 }
108 system("date -u");
109
110 // --- Choix du generateur aleatoire (a faire ICI imperativement avant AutoInitRand)
111 FMTRandGen *RandGen = new FMTRandGen;
112 RandGen->SelectGaussianAlgo(C_Gaussian_RandLibSNorm);
113 RandomGeneratorInterface::SetGlobalRandGenP(RandGen);
114
115 // --- Decodage des arguments
116 char c;
117 while((c = getopt(narg,arg,"haG:F:x:y:z:Z:128:v:n:CO:So:VT:")) != -1) {
118 int nth = 0;
119 switch (c) {
120 case 'a' :
121 AutoInitRand(5);
122 break;
123 case 'G' :
124 use_growth_factor = atoi(optarg);
125 break;
126 case 'F' :
127 filter_by_pixel = atoi(optarg);
128 break;
129 case 'x' :
130 sscanf(optarg,"%ld,%lf",&nx,&dx);
131 break;
132 case 'y' :
133 sscanf(optarg,"%ld,%lf",&ny,&dy);
134 break;
135 case 'z' :
136 sscanf(optarg,"%ld,%lf",&nz,&dz);
137 break;
138 case 'Z' :
139 zref = atof(optarg);
140 break;
141 case 'v' :
142 comptransveloc = true;
143 temporay_file = optarg;
144 break;
145 case '1' :
146 compdspec |= 1;
147 break;
148 case '2' :
149 compdspec |= 2;
150 break;
151 case 'C' :
152 compdspec |= 4;
153 break;
154 case '8' :
155 sscanf(optarg,"%lf,%lf",&sigmaR,&R);
156 break;
157 case 'V' :
158 compvarreal = true;
159 break;
160 case 'n':
161 ntnent = atol(optarg);
162 if(ntnent<0) ntnent = 0;
163 break;
164 case 'O' :
165 sscanf(optarg,"%hu,%hu",&whattowrt[0],&whattowrt[1]);
166 break;
167 case 'S' :
168 wslice = true;
169 break;
170 case 'o' :
171 rootnameout = optarg;
172 break;
173 case 'T' :
174 nth = atoi(optarg);
175 nthread = (nth<1)? 0: nth;
176 break;
177 case 'h' :
178 default :
179 usage(); return -1;
180 }
181 }
182
183 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
184 try {
185 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
186
187 //-----------------------------------------------------------------
188 cout<<endl<<"\n--- Initialisation and parameters"<<endl;
189 cout<<"zref="<<zref<<endl;
190 cout<<"nx="<<nx<<" dx="<<dx<<" ny="<<ny<<" dy="<<dy<<" nz="<<nz<<" dz="<<dz<<endl;
191 cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<") Mpc^-1"<<", npt="<<npt<<endl;
192 cout<<"Filter by pixel = "<<filter_by_pixel<<endl;
193 if(comptransveloc) cout<<"Tansverse velocity generation requested"<<endl;
194 cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
195 cout<<"Use_growth_factor = "<<use_growth_factor<<endl;
196 cout<<"wslice="<<wslice<<" what?="<<whattowrt[0]<<","<<whattowrt[1]<<endl;
197 cout<<"rootnameout="<<rootnameout<<endl;
198 ShowRandom();
199 cout<<" First random is: "<<drand01()<<endl;
200
201 string tagobs = rootnameout + ".ppf";
202 cout<<"Writing result in "<<tagobs<<endl;
203 POutPersist posobs(tagobs);
204
205 //-----------------------------------------------------------------
206 cout<<endl<<"\n--- Create Cosmology"<<endl;
207
208 CosmoCalc univ(flat,true,zref+1.);
209 univ.SetInteg(perc,dzinc,dzmax,glorder);
210 univ.SetDynParam(h100,om0,or0,ol0,w0);
211 univ.PrtInteg();
212 univ.Print();
213 double loscomref = univ.Dloscom(zref);
214 cout<<"\nzref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
215 univ.Print(zref);
216
217 //-----------------------------------------------------------------
218 cout<<endl<<"\n--- Create Spectrum"<<endl;
219
220 InitialPowerLaw pkini(ns,as);
221
222 TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
223 TransfertEisenstein tfnos(h100,om0-ob0,ob0,T_CMB_Par,false);
224 tfnos.SetNoOscEnv(2);
225
226 GrowthEisenstein growth(om0,ol0);
227 // GrowthEisenstein growth(1.,0.); // D(z) = 1/(1+z)
228 double growth_at_z = growth(zref);
229 cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
230
231 PkSpecCalc pkz(pkini,tf,growth,zref);
232 PkSpecCalc pkznos(pkini,tfnos,growth,zref);
233
234 //-----------------------------------------------------------------
235 pkz.SetZ(0.);
236 pkznos.SetZ(0.);
237 cout<<endl<<"\n--- Compute variance for top-hat R="<<R<<" at z="<<pkz.GetZ()<<endl;
238 VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
239 double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
240 double pkmax_th = varpk_th(kfind_th);
241 cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
242 double k1=kmin, k2=kmax;
243 int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
244 cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
245 <<k2<<" ("<<log10(k2)<<")"<<endl;
246
247 double ldlk = (log10(k2)-log10(k1))/npt;
248 varpk_th.SetInteg(0.01,ldlk,-1.,4);
249 double sr2 = varpk_th.Variance(k1,k2);
250 cout<<"varpk_th="<<sr2<<" -> sigma="<<sqrt(sr2)<<endl;
251
252 double normpkz = sigmaR*sigmaR/sr2;
253 pkz.SetScale(normpkz);
254 pkznos.SetScale(normpkz);
255 cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
256 {
257 Histo hpkz0(lkmin,lkmax,npt); hpkz0.ReCenterBin();
258 FuncToHisto(pkz,hpkz0,true);
259 posobs.PutObject(hpkz0,"hpkz0");
260 Histo hpkz0nos(hpkz0);
261 FuncToHisto(pkznos,hpkz0nos,true);
262 posobs.PutObject(hpkz0nos,"hpkz0nos");
263 }
264
265 pkz.SetZ(zref);
266 pkznos.SetZ(zref);
267 {
268 Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
269 FuncToHisto(pkz,hpkz,true);
270 posobs.PutObject(hpkz,"hpkz");
271 Histo hpkznos(hpkz);
272 FuncToHisto(pkznos,hpkznos,true);
273 posobs.PutObject(hpkznos,"hpkznos");
274 }
275
276 //-----------------------------------------------------------------
277 cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
278 VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
279 double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
280 double pkmax_int = varpk_int(kfind_int);
281 cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
282 double k1int=kmin, k2int=kmax;
283 int rcint = varpk_int.FindLimits(pkmax_int/1.e4,k1int,k2int,eps);
284 cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
285 <<k2int<<" ("<<log10(k2int)<<")"<<endl;
286
287 double ldlkint = (log10(k2int)-log10(k1int))/npt;
288 varpk_int.SetInteg(0.01,ldlkint,-1.,4);
289 double sr2int = varpk_int.Variance(k1int,k2int);
290 cout<<"varpk_int="<<sr2int<<" -> sigma="<<sqrt(sr2int)<<endl;
291
292 //-----------------------------------------------------------------
293
294 PrtTim(">>>> End of definition");
295
296 //-----------------------------------------------------------------
297 // Le cube et sa cosmlogie
298 // FFTW3 (p26): faster if sizes 2^a 3^b 5^c 7^d 11^e 13^f with e+f=0 ou 1
299 cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
300
301 GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
302 fluct3d.SetCosmology(univ);
303 fluct3d.SetGrowthFactor(growth);
304 fluct3d.SetObservator(zref,-nz/2.);
305 fluct3d.LosComRedshift(0.001,-1);
306 //TArray< complex<GEN3D_TYPE> >& pkgen = fluct3d.GetComplexArray();
307 //TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
308 cout<<endl; fluct3d.Print();
309
310 //-----------------------------------------------------------------
311 // Les bins et ranges en k pour les histos 1D et 2D
312 double dkmin = fluct3d.GetKincMin();
313 double knyqmax = fluct3d.GetKmax();
314 long nherr = long(knyqmax/dkmin+0.5);
315 cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
316
317 double dktmin = fluct3d.GetKTincMin();
318 double ktnyqmax = fluct3d.GetKTmax();
319 long nherrt = long(ktnyqmax/dktmin+0.5);
320 double dkzmin = fluct3d.GetKinc()[2];
321 double kznyqmax = fluct3d.GetKnyq()[2];
322 long nherrz = long(kznyqmax/dkzmin+0.5);
323 cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
324 <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
325
326 //-----------------------------------------------------------------
327 cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
328 // En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
329 // sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote k*sqrt(2)): Vc = (k*sqrt(2))^3
330 // Vc/Vs = 0.675 -> keff = kmax * (0.675)^(1/3) = kmax * 0.877
331 double knyqmax_mod = 0.877*knyqmax;
332 ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
333 varpk_int.SetInteg(0.01,ldlkint,-1.,4);
334 double sr2int_kmax = varpk_int.Variance(k1int,knyqmax_mod);
335 cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<" -> sigma="<<sqrt(sr2int_kmax)<<endl;
336
337 PrtTim(">>>> End Initialisation of GeneFluct3D");
338
339 //-----------------------------------------------------------------
340 cout<<"\n--- Computing a realization in Fourier space"<<endl;
341 if(use_growth_factor>0) pkz.SetZ(0.); else pkz.SetZ(zref);
342 cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
343 fluct3d.ComputeFourier(pkz);
344 fluct3d.NTupleCheck(posobs,string("ntpkgen"),ntnent);
345 PrtTim(">>>> End Computing a realization in Fourier space");
346
347 HistoErr hpkgen(0.,knyqmax,nherr);
348 if(compdspec&1) {
349 cout<<"\n--- Checking realization spectra"<<endl;
350 hpkgen.ReCenterBin(); hpkgen.Zero(); hpkgen.Show();
351 fluct3d.ComputeSpectrum(hpkgen);
352 posobs.PutObject(hpkgen,"hpkgen");
353 PrtTim(">>>> End Checking realization 1D spectra");
354 }
355 Histo2DErr hpkgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
356 if(compdspec&2) {
357 cout<<"\n--- Checking realization 2D spectra"<<endl;
358 hpkgen2.ReCenterBin(); hpkgen2.Zero(); hpkgen2.Show();
359 fluct3d.ComputeSpectrum2D(hpkgen2);
360 posobs.PutObject(hpkgen2,"hpkgen2");
361 PrtTim(">>>> End Checking realization 2D spectra");
362 }
363
364 //-----------------------------------------------------------------
365 if(filter_by_pixel!=0) {
366 cout<<"\n--- Computing convolution by pixel shape"<<endl;
367 fluct3d.FilterByPixel();
368 fluct3d.NTupleCheck(posobs,string("ntpkgenf"),ntnent);
369 PrtTim(">>>> End Computing convolution by pixel shape");
370 }
371
372 if(compdspec&1) {
373 HistoErr hpkgenf(hpkgen);
374 if(filter_by_pixel!=0) {
375 cout<<"\n--- Checking realization spectra"<<endl;
376 hpkgenf.Zero(); hpkgenf.Show();
377 fluct3d.ComputeSpectrum(hpkgenf);
378 PrtTim(">>>> End Checking realization 1D spectra");
379 }
380 posobs.PutObject(hpkgenf,"hpkgenf");
381 }
382 if(compdspec&2) {
383 Histo2DErr hpkgenf2(hpkgen2);
384 if(filter_by_pixel!=0) {
385 cout<<"\n--- Checking realization 2D spectra"<<endl;
386 hpkgenf2.Zero(); hpkgenf2.Show();
387 fluct3d.ComputeSpectrum2D(hpkgenf2);
388 PrtTim(">>>> End Checking realization 2D spectra");
389 }
390 posobs.PutObject(hpkgenf2,"hpkgenf2");
391 }
392
393 //-----------------------------------------------------------------
394 if(whattowrt[0]&1) {
395 tagobs = rootnameout + "_k.ppf";
396 fluct3d.WritePPF(tagobs,false);
397 PrtTim(">>>> End WritePPF");
398 temporay_file = tagobs;
399 }
400 if(whattowrt[0]&2) {
401 tagobs = "!" + rootnameout + "_k.fits";
402 fluct3d.WriteFits(tagobs);
403 PrtTim(">>>> End WriteFits");
404 }
405 if(comptransveloc && !(whattowrt[0]&1)) {
406 cout<<">>> Writing FT cube (for transv veloc.) in "<<temporay_file<<endl;
407 fluct3d.WritePPF(temporay_file,false);
408 PrtTim(">>>> End Writing FT cube");
409 }
410
411 //-----------------------------------------------------------------
412 cout<<"\n--- Computing a realization in real space"<<endl;
413 fluct3d.ComputeReal();
414 double rmin,rmax; fluct3d.MinMax(rmin,rmax);
415 cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
416 fluct3d.NTupleCheck(posobs,string("ntreal"),ntnent);
417 PrtTim(">>>> End Computing a realization in real space");
418
419 if(use_growth_factor>0) {
420 cout<<"\n--- Apply Growth factor"<<endl;
421 cout<<"...D(z=0)="<<growth(0.)<<" D(z="<<zref<<")="<<growth(zref)<<endl;
422 fluct3d.ApplyGrowthFactor(use_growth_factor);
423 fluct3d.MinMax(rmin,rmax);
424 cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
425 fluct3d.NTupleCheck(posobs,string("ntgrow"),ntnent);
426 PrtTim(">>>> End Applying growth factor");
427 }
428
429 int_8 nm;
430 double rmref,rs2ref;
431 cout<<"\n--- Computing reference variance in real space"<<endl;
432 nm = fluct3d.MeanSigma2(rmref,rs2ref);
433 cout<<" rs2ref= "<<rs2ref<<" , rmref="<<rmref<<" ("<<nm<<")"<<endl;
434 PrtTim(">>>> End Computing reference variance in real space");
435
436 if(whattowrt[1]&1) {
437 tagobs = rootnameout + "_r.ppf";
438 fluct3d.WritePPF(tagobs,true);
439 PrtTim(">>>> End WritePPF");
440 }
441 if(whattowrt[1]&2) {
442 tagobs = "!" + rootnameout + "_r.fits";
443 fluct3d.WriteFits(tagobs);
444 PrtTim(">>>> End WriteFits");
445 }
446 if(wslice) {
447 tagobs = rootnameout + "_s_r.ppf";
448 fluct3d.WriteSlicePPF(tagobs);
449 PrtTim(">>>> End WriteSlicePPF");
450 }
451
452 //-----------------------------------------------------------------
453 if(compvarreal) {
454 cout<<"\n--- Check variance sigmaR in real space"<<endl;
455 double varr;
456 fluct3d.VarianceFrReal(R,varr);
457 cout<<"...Computed variance = "<<varr
458 <<" , Theorical variance at (z=0) = "<<pow(sigmaR,2.)
459 <<" , at (z="<<zref<<") = "<<pow(sigmaR*growth_at_z,2.)<<endl;
460 PrtTim(">>>> End Check variance sigmaR in real space");
461 }
462
463 //-----------------------------------------------------------------
464 if(compdspec&4) { // ATTENTION: d_rho/rho ecrase
465 cout<<endl<<"\n--- ReComputing spectrum from real space"
466 <<"\n Warning: REAL SPACE CUBE OVERWRITTEN"<<endl;
467 fluct3d.ReComputeFourier();
468 fluct3d.NTupleCheck(posobs,string("ntpkrec"),ntnent);
469 PrtTim(">>>> End ReComputing spectrum");
470
471 cout<<endl<<"\n--- Computing final spectrum"<<endl;
472 HistoErr hpkrecb(0.,knyqmax,nherr); hpkrecb.Zero();
473 hpkrecb.ReCenterBin(); hpkrecb.Show();
474 fluct3d.ComputeSpectrum(hpkrecb);
475 posobs.PutObject(hpkrecb,"hpkrecb");
476 PrtTim(">>>> End Computing final spectrum");
477
478 HistoErr hpkrec(hpkrecb);
479 if(filter_by_pixel!=0) {
480 cout<<endl<<"\n--- Computing final spectrum with pixel deconv."<<endl;
481 hpkrec.Zero();
482 fluct3d.ComputeSpectrum(hpkrec,0.,filter_by_pixel);
483 PrtTim(">>>> End Computing final spectrum with pixel deconv.");
484 }
485 posobs.PutObject(hpkrec,"hpkrec");
486
487 cout<<"\n--- Computing final 2D spectrum"<<endl;
488 Histo2DErr hpkrecb2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
489 hpkrecb2.ReCenterBin(); hpkrecb2.Zero(); hpkrecb2.Show();
490 fluct3d.ComputeSpectrum2D(hpkrecb2);
491 posobs.PutObject(hpkrecb2,"hpkrecb2");
492 PrtTim(">>>> End Computing final 2D spectrum");
493
494 Histo2DErr hpkrec2(hpkrecb2);
495 if(filter_by_pixel!=0) {
496 cout<<"\n--- Computing final 2D spectrum with pixel deconv."<<endl;
497 hpkrec2.Zero();
498 fluct3d.ComputeSpectrum2D(hpkrec2,0.,filter_by_pixel);
499 PrtTim(">>>> End Computing final 2D spectrum with pixel deconv.");
500 }
501 posobs.PutObject(hpkrec2,"hpkrec2");
502 }
503
504 //-----------------------------------------------------------------
505 //-----------------------------------------------------------------
506 //-----------------------------------------------------------------
507
508 if(comptransveloc) {
509 cout<<"\n\n\n"<<"---------------------------------------------\n"
510 <<"--- Transverse velocity space computation ---\n"
511 <<"---------------------------------------------\n"<<endl;
512
513 //-----------------------------------------------------------------
514 cout<<"\n--- Reading back the Pk(vec(k)) cube from "<<temporay_file<<endl;
515 fluct3d.ReadPPF(temporay_file);
516 PrtTim(">>>> End Reading the Pk(vec(k)) cube");
517
518 if(compdspec&1) {
519 cout<<"\n--- Check realization spectra that has been re-read"<<endl;
520 HistoErr hpkgenf_rr(hpkgen); hpkgenf_rr.Zero();
521 fluct3d.ComputeSpectrum(hpkgenf_rr);
522 PrtTim(">>>> End Checking re-read realization 1D spectra");
523 posobs.PutObject(hpkgenf_rr,"hpkgenf_rr");
524 }
525
526 //-----------------------------------------------------------------
527 cout<<"\n--- Modifying cube for Transverse velocity"<<endl;
528 fluct3d.ToVelLoS();
529 fluct3d.NTupleCheck(posobs,string("ntpvgenf"),ntnent);
530 PrtTim(">>>> End Modifying cube for Transverse velocity");
531
532 if(compdspec&1) {
533 cout<<"\n--- Checking realization spectra"<<endl;
534 HistoErr hpvgen(0.,knyqmax,nherr);
535 hpvgen.ReCenterBin(); hpvgen.Zero(); hpvgen.Show();
536 fluct3d.ComputeSpectrum(hpvgen);
537 posobs.PutObject(hpvgen,"hpvgen");
538 PrtTim(">>>> End Checking realization 1D spectra");
539 }
540 if(compdspec&2) {
541 cout<<"\n--- Checking realization 2D spectra"<<endl;
542 Histo2DErr hpvgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
543 hpvgen2.ReCenterBin(); hpvgen2.Zero(); hpvgen2.Show();
544 fluct3d.ComputeSpectrum2D(hpvgen2);
545 posobs.PutObject(hpvgen2,"hpvgen2");
546 PrtTim(">>>> End Checking realization 2D spectra");
547 }
548
549 if(whattowrt[0]&1) {
550 tagobs = rootnameout + "_kv.ppf";
551 fluct3d.WritePPF(tagobs,false);
552 PrtTim(">>>> End WritePPF");
553 }
554 if(whattowrt[0]&2) {
555 tagobs = "!" + rootnameout + "_kv.fits";
556 fluct3d.WriteFits(tagobs);
557 PrtTim(">>>> End WriteFits");
558 }
559
560 //-----------------------------------------------------------------
561 cout<<"\n--- Computing a realization in real space for Transverse velocity"<<endl;
562 fluct3d.ComputeReal();
563 double rmin,rmax; fluct3d.MinMax(rmin,rmax);
564 cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
565 fluct3d.NTupleCheck(posobs,string("nvreal"),ntnent);
566 PrtTim(">>>> End Computing a realization in real space");
567
568 if(use_growth_factor>0) {
569 cout<<"\n--- Apply Growth factor"<<endl;
570 cout<<"...D(z=0)="<<growth(0.)<<" D(z="<<zref<<")="<<growth(zref)<<endl;
571 fluct3d.ApplyVelLosGrowthFactor(use_growth_factor);
572 fluct3d.MinMax(rmin,rmax);
573 cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
574 fluct3d.NTupleCheck(posobs,string("nvgrow"),ntnent);
575 PrtTim(">>>> End Applying growth factor");
576 }
577
578 int_8 nmv;
579 double vmref,vs2ref;
580 cout<<"\n--- Computing reference variance in real space"<<endl;
581 nmv = fluct3d.MeanSigma2(vmref,vs2ref);
582 cout<<" vs2ref= "<<vs2ref<<" , vmref="<<vmref<<" ("<<nmv<<")"<<endl;
583 PrtTim(">>>> End Computing reference variance in real space");
584
585 if(whattowrt[1]&1) {
586 tagobs = rootnameout + "_rv.ppf";
587 fluct3d.WritePPF(tagobs,true);
588 PrtTim(">>>> End WritePPF");
589 }
590 if(whattowrt[1]&2) {
591 tagobs = "!" + rootnameout + "_rv.fits";
592 fluct3d.WriteFits(tagobs);
593 PrtTim(">>>> End WriteFits");
594 }
595 if(wslice) {
596 tagobs = rootnameout + "_s_rv.ppf";
597 fluct3d.WriteSlicePPF(tagobs);
598 PrtTim(">>>> End WriteSlicePPF");
599 }
600
601 }
602
603 //-----------------------------------------------------------------
604 {
605 DVList dvl;
606 dvl("Nx") = (int_4)fluct3d.Nx(); dvl("Ny") = (int_4)fluct3d.Ny(); dvl("Nz") = (int_4)fluct3d.Nz();
607 dvl("Dx") = fluct3d.Dx(); dvl("Dy") = fluct3d.Dy(); dvl("Dz") = fluct3d.Dz();
608 dvl("Z") = fluct3d.Zref(); dvl("Zpk") = fluct3d.ZrefPk();
609 dvl("D") = fluct3d.Dref(); dvl("Dpk") = fluct3d.DrefPk();
610 dvl("s8") = sigmaR;
611 dvl("Dtype") = (int_4)use_growth_factor; dvl("Pxfilt") = (int_4)filter_by_pixel;
612 posobs.PutObject(dvl,"siminfo");
613 }
614
615 delete RandGen;
616 PrtTim(">>>> End Of Job");
617
618 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
619 } catch (PException& exc) {
620 cerr<<"cmvginit3d.cc catched PException"<<exc.Msg()<<endl;
621 return 77;
622 } catch (std::exception& sex) {
623 cerr << "cmvginit3d.cc std::exception :"
624 << (string)typeid(sex).name() << "\n msg= "
625 << sex.what() << endl;
626 return 78;
627 } catch (...) {
628 cerr << "cmvginit3d.cc catched unknown (...) exception " << endl;
629 return 79;
630 }
631 //----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
632
633 return 0;
634}
635
Note: See TracBrowser for help on using the repository browser.