Context Navigation

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

Visit:

Last change on this file since 3771 was 3770, checked in by cmv, 15 years ago
relecture des redshift-distorsion + calcul des spectres, cmv 07/05/2010
File size: 21.4 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
23	void usage(void);
24	void usage(void)
25	{
26	cout<<"cmvobserv3d [...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 (with the -W and -P options)"<<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 root_name_out : root string for output file name (def: cmvobserv3d)"<<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
57	int 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	InitialSpectrum pkini(ns,as);
221
222	TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
223	//tf.SetNoOscEnv(2);
224
225	GrowthFactor growth(om0,ol0);
226	// GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
227	double growth_at_z = growth(zref);
228	cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
229
230	PkSpectrum0 pk0(pkini,tf);
231
232	PkSpectrumZ pkz(pk0,growth,zref);
233
234	//-----------------------------------------------------------------
235	pkz.SetZ(0.);
236	cout<<endl<<"\n--- Compute variance for top-hat R="<<R<<" at z="<<pkz.GetZ()<<endl;
237	VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
238	double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
239	double pkmax_th = varpk_th(kfind_th);
240	cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
241	double k1=kmin, k2=kmax;
242	int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
243	cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
244	<<k2<<" ("<<log10(k2)<<")"<<endl;
245
246	double ldlk = (log10(k2)-log10(k1))/npt;
247	varpk_th.SetInteg(0.01,ldlk,-1.,4);
248	double sr2 = varpk_th.Variance(k1,k2);
249	cout<<"varpk_th="<<sr2<<" -> sigma="<<sqrt(sr2)<<endl;
250
251	double normpkz = sigmaR*sigmaR/sr2;
252	pkz.SetScale(normpkz);
253	cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
254	{
255	Histo hpkz0(lkmin,lkmax,npt); hpkz0.ReCenterBin();
256	FuncToHisto(pkz,hpkz0,true);
257	posobs.PutObject(hpkz0,"hpkz0");
258	}
259
260	pkz.SetZ(zref);
261	{
262	Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
263	FuncToHisto(pkz,hpkz,true);
264	posobs.PutObject(hpkz,"hpkz");
265	}
266
267	//-----------------------------------------------------------------
268	cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
269	VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
270	double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
271	double pkmax_int = varpk_int(kfind_int);
272	cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
273	double k1int=kmin, k2int=kmax;
274	int rcint = varpk_int.FindLimits(pkmax_int/1.e4,k1int,k2int,eps);
275	cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
276	<<k2int<<" ("<<log10(k2int)<<")"<<endl;
277
278	double ldlkint = (log10(k2int)-log10(k1int))/npt;
279	varpk_int.SetInteg(0.01,ldlkint,-1.,4);
280	double sr2int = varpk_int.Variance(k1int,k2int);
281	cout<<"varpk_int="<<sr2int<<" -> sigma="<<sqrt(sr2int)<<endl;
282
283	//-----------------------------------------------------------------
284
285	PrtTim(">>>> End of definition");
286
287	//-----------------------------------------------------------------
288	// Le cube et sa cosmlogie
289	// FFTW3 (p26): faster if sizes 2^a 3^b 5^c 7^d 11^e 13^f with e+f=0 ou 1
290	cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
291
292	GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
293	fluct3d.SetCosmology(univ);
294	fluct3d.SetGrowthFactor(growth);
295	fluct3d.SetObservator(zref,-nz/2.);
296	fluct3d.LosComRedshift(0.001,-1);
297	//TArray< complex<GEN3D_TYPE> >& pkgen = fluct3d.GetComplexArray();
298	//TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
299	cout<<endl; fluct3d.Print();
300
301	//-----------------------------------------------------------------
302	// Les bins et ranges en k pour les histos 1D et 2D
303	double dkmin = fluct3d.GetKincMin();
304	double knyqmax = fluct3d.GetKmax();
305	long nherr = long(knyqmax/dkmin+0.5);
306	cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
307
308	double dktmin = fluct3d.GetKTincMin();
309	double ktnyqmax = fluct3d.GetKTmax();
310	long nherrt = long(ktnyqmax/dktmin+0.5);
311	double dkzmin = fluct3d.GetKinc()[2];
312	double kznyqmax = fluct3d.GetKnyq()[2];
313	long nherrz = long(kznyqmax/dkzmin+0.5);
314	cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
315	<<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
316
317	//-----------------------------------------------------------------
318	cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
319	// En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
320	// sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote ksqrt(2)): Vc = (ksqrt(2))^3
321	// Vc/Vs = 0.675 -> keff = kmax * (0.675)^(1/3) = kmax * 0.877
322	double knyqmax_mod = 0.877*knyqmax;
323	ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
324	varpk_int.SetInteg(0.01,ldlkint,-1.,4);
325	double sr2int_kmax = varpk_int.Variance(k1int,knyqmax_mod);
326	cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<" -> sigma="<<sqrt(sr2int_kmax)<<endl;
327
328	PrtTim(">>>> End Initialisation of GeneFluct3D");
329
330	//-----------------------------------------------------------------
331	cout<<"\n--- Computing a realization in Fourier space"<<endl;
332	if(use_growth_factor>0) pkz.SetZ(0.); else pkz.SetZ(zref);
333	cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
334	fluct3d.ComputeFourier(pkz);
335	fluct3d.NTupleCheck(posobs,string("ntpkgen"),ntnent);
336	PrtTim(">>>> End Computing a realization in Fourier space");
337
338	HistoErr hpkgen(0.,knyqmax,nherr);
339	if(compdspec&1) {
340	cout<<"\n--- Checking realization spectra"<<endl;
341	hpkgen.ReCenterBin(); hpkgen.Zero(); hpkgen.Show();
342	fluct3d.ComputeSpectrum(hpkgen);
343	posobs.PutObject(hpkgen,"hpkgen");
344	PrtTim(">>>> End Checking realization 1D spectra");
345	}
346	Histo2DErr hpkgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
347	if(compdspec&2) {
348	cout<<"\n--- Checking realization 2D spectra"<<endl;
349	hpkgen2.ReCenterBin(); hpkgen2.Zero(); hpkgen2.Show();
350	fluct3d.ComputeSpectrum2D(hpkgen2);
351	posobs.PutObject(hpkgen2,"hpkgen2");
352	PrtTim(">>>> End Checking realization 2D spectra");
353	}
354
355	//-----------------------------------------------------------------
356	if(filter_by_pixel!=0) {
357	cout<<"\n--- Computing convolution by pixel shape"<<endl;
358	fluct3d.FilterByPixel();
359	fluct3d.NTupleCheck(posobs,string("ntpkgenf"),ntnent);
360	PrtTim(">>>> End Computing convolution by pixel shape");
361	}
362
363	if(compdspec&1) {
364	HistoErr hpkgenf(hpkgen);
365	if(filter_by_pixel!=0) {
366	cout<<"\n--- Checking realization spectra"<<endl;
367	hpkgenf.Zero(); hpkgenf.Show();
368	fluct3d.ComputeSpectrum(hpkgenf);
369	PrtTim(">>>> End Checking realization 1D spectra");
370	}
371	posobs.PutObject(hpkgenf,"hpkgenf");
372	}
373	if(compdspec&2) {
374	Histo2DErr hpkgenf2(hpkgen2);
375	if(filter_by_pixel!=0) {
376	cout<<"\n--- Checking realization 2D spectra"<<endl;
377	hpkgenf2.Zero(); hpkgenf2.Show();
378	fluct3d.ComputeSpectrum2D(hpkgenf2);
379	PrtTim(">>>> End Checking realization 2D spectra");
380	}
381	posobs.PutObject(hpkgenf2,"hpkgenf2");
382	}
383
384	//-----------------------------------------------------------------
385	if(whattowrt[0]&1) {
386	tagobs = rootnameout + "_k.ppf";
387	fluct3d.WritePPF(tagobs,false);
388	PrtTim(">>>> End WritePPF");
389	temporay_file = tagobs;
390	}
391	if(whattowrt[0]&2) {
392	tagobs = "!" + rootnameout + "_k.fits";
393	fluct3d.WriteFits(tagobs);
394	PrtTim(">>>> End WriteFits");
395	}
396	if(comptransveloc && !(whattowrt[0]&1)) {
397	cout<<">>> Writing FT cube (for transv veloc.) in "<<temporay_file<<endl;
398	fluct3d.WritePPF(temporay_file,false);
399	PrtTim(">>>> End Writing FT cube");
400	}
401
402	//-----------------------------------------------------------------
403	cout<<"\n--- Computing a realization in real space"<<endl;
404	fluct3d.ComputeReal();
405	double rmin,rmax; fluct3d.MinMax(rmin,rmax);
406	cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
407	fluct3d.NTupleCheck(posobs,string("ntreal"),ntnent);
408	PrtTim(">>>> End Computing a realization in real space");
409
410	if(use_growth_factor>0) {
411	cout<<"\n--- Apply Growth factor"<<endl;
412	cout<<"...D(z=0)="<<growth(0.)<<" D(z="<<zref<<")="<<growth(zref)<<endl;
413	fluct3d.ApplyGrowthFactor(use_growth_factor);
414	fluct3d.MinMax(rmin,rmax);
415	cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
416	fluct3d.NTupleCheck(posobs,string("ntgrow"),ntnent);
417	PrtTim(">>>> End Applying growth factor");
418	}
419
420	int_8 nm;
421	double rmref,rs2ref;
422	cout<<"\n--- Computing reference variance in real space"<<endl;
423	nm = fluct3d.MeanSigma2(rmref,rs2ref);
424	cout<<" rs2ref= "<<rs2ref<<" , rmref="<<rmref<<" ("<<nm<<")"<<endl;
425	PrtTim(">>>> End Computing reference variance in real space");
426
427	if(whattowrt[1]&1) {
428	tagobs = rootnameout + "_r.ppf";
429	fluct3d.WritePPF(tagobs,true);
430	PrtTim(">>>> End WritePPF");
431	}
432	if(whattowrt[1]&2) {
433	tagobs = "!" + rootnameout + "_r.fits";
434	fluct3d.WriteFits(tagobs);
435	PrtTim(">>>> End WriteFits");
436	}
437	if(wslice) {
438	tagobs = rootnameout + "_s_r.ppf";
439	fluct3d.WriteSlicePPF(tagobs);
440	PrtTim(">>>> End WriteSlicePPF");
441	}
442
443	//-----------------------------------------------------------------
444	if(compvarreal) {
445	cout<<"\n--- Check variance sigmaR in real space"<<endl;
446	double varr;
447	fluct3d.VarianceFrReal(R,varr);
448	cout<<"...Computed variance = "<<varr
449	<<" , Theorical variance at (z=0) = "<<pow(sigmaR,2.)
450	<<" , at (z="<<zref<<") = "<<pow(sigmaR*growth_at_z,2.)<<endl;
451	PrtTim(">>>> End Check variance sigmaR in real space");
452	}
453
454	//-----------------------------------------------------------------
455	if(compdspec&4) { // ATTENTION: d_rho/rho ecrase
456	cout<<endl<<"\n--- ReComputing spectrum from real space"
457	<<"\n Warning: REAL SPACE CUBE OVERWRITTEN"<<endl;
458	fluct3d.ReComputeFourier();
459	fluct3d.NTupleCheck(posobs,string("ntpkrec"),ntnent);
460	PrtTim(">>>> End ReComputing spectrum");
461
462	cout<<endl<<"\n--- Computing final spectrum"<<endl;
463	HistoErr hpkrecb(0.,knyqmax,nherr); hpkrecb.Zero();
464	hpkrecb.ReCenterBin(); hpkrecb.Show();
465	fluct3d.ComputeSpectrum(hpkrecb);
466	posobs.PutObject(hpkrecb,"hpkrecb");
467	PrtTim(">>>> End Computing final spectrum");
468
469	HistoErr hpkrec(hpkrecb);
470	if(filter_by_pixel!=0) {
471	cout<<endl<<"\n--- Computing final spectrum with pixel deconv."<<endl;
472	hpkrec.Zero();
473	fluct3d.ComputeSpectrum(hpkrec,0.,filter_by_pixel);
474	PrtTim(">>>> End Computing final spectrum with pixel deconv.");
475	}
476	posobs.PutObject(hpkrec,"hpkrec");
477
478	cout<<"\n--- Computing final 2D spectrum"<<endl;
479	Histo2DErr hpkrecb2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
480	hpkrecb2.ReCenterBin(); hpkrecb2.Zero(); hpkrecb2.Show();
481	fluct3d.ComputeSpectrum2D(hpkrecb2);
482	posobs.PutObject(hpkrecb2,"hpkrecb2");
483	PrtTim(">>>> End Computing final 2D spectrum");
484
485	Histo2DErr hpkrec2(hpkrecb2);
486	if(filter_by_pixel!=0) {
487	cout<<"\n--- Computing final 2D spectrum with pixel deconv."<<endl;
488	hpkrec2.Zero();
489	fluct3d.ComputeSpectrum2D(hpkrec2,0.,filter_by_pixel);
490	PrtTim(">>>> End Computing final 2D spectrum with pixel deconv.");
491	}
492	posobs.PutObject(hpkrec2,"hpkrec2");
493	}
494
495	//-----------------------------------------------------------------
496	//-----------------------------------------------------------------
497	//-----------------------------------------------------------------
498
499	if(comptransveloc) {
500	cout<<"\n\n\n"<<"---------------------------------------------\n"
501	<<"--- Transverse velocity space computation ---\n"
502	<<"---------------------------------------------\n"<<endl;
503
504	//-----------------------------------------------------------------
505	cout<<"\n--- Reading back the Pk(vec(k)) cube from "<<temporay_file<<endl;
506	fluct3d.ReadPPF(temporay_file);
507	PrtTim(">>>> End Reading the Pk(vec(k)) cube");
508
509	//-----------------------------------------------------------------
510	cout<<"\n--- Modifying cube for Transverse velocity"<<endl;
511	fluct3d.ToVelLoS();
512	fluct3d.NTupleCheck(posobs,string("ntpvgenf"),ntnent);
513	PrtTim(">>>> End Modifying cube for Transverse velocity");
514
515	if(compdspec&1) {
516	cout<<"\n--- Checking realization spectra"<<endl;
517	HistoErr hpvgen(0.,knyqmax,nherr);
518	hpvgen.ReCenterBin(); hpvgen.Zero(); hpvgen.Show();
519	fluct3d.ComputeSpectrum(hpvgen);
520	posobs.PutObject(hpvgen,"hpvgen");
521	PrtTim(">>>> End Checking realization 1D spectra");
522	}
523	if(compdspec&2) {
524	cout<<"\n--- Checking realization 2D spectra"<<endl;
525	Histo2DErr hpvgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
526	hpvgen2.ReCenterBin(); hpvgen2.Zero(); hpvgen2.Show();
527	fluct3d.ComputeSpectrum2D(hpvgen2);
528	posobs.PutObject(hpvgen2,"hpvgen2");
529	PrtTim(">>>> End Checking realization 2D spectra");
530	}
531
532	if(whattowrt[0]&1) {
533	tagobs = rootnameout + "_kv.ppf";
534	fluct3d.WritePPF(tagobs,false);
535	PrtTim(">>>> End WritePPF");
536	}
537	if(whattowrt[0]&2) {
538	tagobs = "!" + rootnameout + "_kv.fits";
539	fluct3d.WriteFits(tagobs);
540	PrtTim(">>>> End WriteFits");
541	}
542
543	//-----------------------------------------------------------------
544	cout<<"\n--- Computing a realization in real space for Transverse velocity"<<endl;
545	fluct3d.ComputeReal();
546	double rmin,rmax; fluct3d.MinMax(rmin,rmax);
547	cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
548	fluct3d.NTupleCheck(posobs,string("nvreal"),ntnent);
549	PrtTim(">>>> End Computing a realization in real space");
550
551	if(use_growth_factor>0) {
552	cout<<"\n--- Apply Growth factor"<<endl;
553	cout<<"...D(z=0)="<<growth(0.)<<" D(z="<<zref<<")="<<growth(zref)<<endl;
554	fluct3d.ApplyDerGrowthFactor(use_growth_factor);
555	fluct3d.MinMax(rmin,rmax);
556	cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
557	fluct3d.NTupleCheck(posobs,string("nvgrow"),ntnent);
558	PrtTim(">>>> End Applying growth factor");
559	}
560
561	int_8 nmv;
562	double vmref,vs2ref;
563	cout<<"\n--- Computing reference variance in real space"<<endl;
564	nmv = fluct3d.MeanSigma2(vmref,vs2ref);
565	cout<<" vs2ref= "<<vs2ref<<" , vmref="<<vmref<<" ("<<nmv<<")"<<endl;
566	PrtTim(">>>> End Computing reference variance in real space");
567
568	if(whattowrt[1]&1) {
569	tagobs = rootnameout + "_rv.ppf";
570	fluct3d.WritePPF(tagobs,true);
571	PrtTim(">>>> End WritePPF");
572	}
573	if(whattowrt[1]&2) {
574	tagobs = "!" + rootnameout + "_rv.fits";
575	fluct3d.WriteFits(tagobs);
576	PrtTim(">>>> End WriteFits");
577	}
578	if(wslice) {
579	tagobs = rootnameout + "_s_rv.ppf";
580	fluct3d.WriteSlicePPF(tagobs);
581	PrtTim(">>>> End WriteSlicePPF");
582	}
583
584	}
585
586	//-----------------------------------------------------------------
587	{
588	DVList dvl;
589	dvl("Nx") = (int_4)fluct3d.Nx(); dvl("Ny") = (int_4)fluct3d.Ny(); dvl("Nz") = (int_4)fluct3d.Nz();
590	dvl("Dx") = fluct3d.Dx(); dvl("Dy") = fluct3d.Dy(); dvl("Dz") = fluct3d.Dz();
591	dvl("Z") = fluct3d.Zref(); dvl("Zpk") = fluct3d.ZrefPk();
592	dvl("D") = fluct3d.Dref(); dvl("Dpk") = fluct3d.DrefPk();
593	dvl("s8") = sigmaR;
594	dvl("Dtype") = (int_4)use_growth_factor; dvl("Pxfilt") = (int_4)filter_by_pixel;
595	posobs.PutObject(dvl,"siminfo");
596	}
597
598	delete RandGen;
599	PrtTim(">>>> End Of Job");
600
601	//----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
602	} catch (PException& exc) {
603	cerr<<"cmvginit3d.cc catched PException"<<exc.Msg()<<endl;
604	return 77;
605	} catch (std::exception& sex) {
606	cerr << "cmvginit3d.cc std::exception :"
607	<< (string)typeid(sex).name() << "\n msg= "
608	<< sex.what() << endl;
609	return 78;
610	} catch (...) {
611	cerr << "cmvginit3d.cc catched unknown (...) exception " << endl;
612	return 79;
613	}
614	//----TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH-TRY-CATCH
615
616	return 0;
617	}
618

Note: See TracBrowser for help on using the repository browser.

Download in other formats: