Context Navigation

source: Sophya/trunk/Cosmo/SimLSS/cmvobserv3d.cc@ 3743

Visit:

Last change on this file since 3743 was 3615, checked in by cmv, 16 years ago
Modifs relatives a l'introduction de RandomGeneratorInterface + delete de srandgen.c, cmv 01/05/2009
File size: 27.3 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	#include "sophyainit.h"
10	#include "timing.h"
11	#include "ntuple.h"
12	#include "matharr.h"
13	#include "randfmt.h"
14	//#include "randr48.h"
15	#include "srandgen.h"
16
17	#include "constcosmo.h"
18	#include "cosmocalc.h"
19	#include "schechter.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	<<" -0 : use ComputeFourier0 method (defaut: no, use normal way)"<<endl
29	<<" -G typevol: compute Pk(z=0) and apply growth factor in real space"<<endl
30	<<" typevol=1 evolved with distance / observateur (def)"<<endl
31	<<" typevol=2 evolved with distance to middle of Z planes"<<endl
32	<<" else : no evol, spectrum Pk(z=z_median) for all cube (def)"<<endl
33	<<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl
34	<<" -U typ =0 compute <NGal>, poisson fluctuate then convert to HI mass (def)"<<endl
35	<<" >0 compute directly <HI mass>, do NOT poisson fluctuate with <Ngal>"<<endl
36	<<" <0 just multiply dRho/Rho by mass_by_pixel (possible negative pixel values)"<<endl
37	<<" -x nx,dx : size along x axis (npix,Mpc)"<<endl
38	<<" -y ny,dy : size along y axis (npix,Mpc)"<<endl
39	<<" if ny or dy <=0 take same value as for x"<<endl
40	<<" -z nz,dz : size along z axis (redshift axis, npix,Mpc)"<<endl
41	<<" -Z zref : redshift for the center of the simulation cube"<<endl
42	<<" -s snoise,typevol : gaussian noise sigma in equivalent Msol"<<endl
43	<<" typevol=0 no evolution (def)"<<endl
44	<<" typevol=1 noise evolved with distance / observateur"<<endl
45	<<" typevol=2 noise evolved with distance to middle of Z planes"<<endl
46	<<" -2 : compute also 2D spectrum (default: no)"<<endl
47	<<" -N scalemass: facteur d\'unite pour la masse (default: -1)"<<endl
48	<<" ex: si on veut unites 10^8 Msol -> scalemass=1.e-8"<<endl
49	<<" si <0 alors facteur=-scalemass*Mpix"<<endl
50	<<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
51	<<" If nsch<0 alors no,bre de points par decade"<<endl
52	<<" -Q naleagal : use quick method for turning ngal to mass"<<endl
53	<<" -R schmassdist.ppf : read mass distribution for trials from file"<<endl
54	<<" instead of computing it (ONLY if \"-Q\" option is activated)"<<endl
55	<<" -8 sigmaR,R : normalisation du spectre de puissance, R en Mpc"<<endl
56	<<" (default sigmaR=1, R=8/h100 Mpc)"<<endl
57	<<" -W : write cube in FITS format (complex cube is coded as real cube)"<<endl
58	<<" -P : write cube in PPF format"<<endl
59	<<" -O a,b,c,d,e : tell what you want to write (with the -W and -P options)"<<endl
60	<<" a=1 : write generated fourier cube (_k0)"<<endl
61	<<" b=1 : write real space cube dRho/Rho at z (_r0)"<<endl
62	<<" c=1 : write real space cube conv to gal,poiss,mass etc... (_r)"<<endl
63	<<" d=1 : write previous real space cube + noise (_rf)"<<endl
64	<<" e=1 : write fourier transform of previous cube (_k)"<<endl
65	<<" -S : write cube slices in PPF format"<<endl
66	<<" -o root_name_out : root string for output file name (def: cmvobserv3d)"<<endl
67	<<" -V : compute variance from real space (for check, default: no)"<<endl
68	<<" -T nth : nombre de threads (si compil multi-thread, default: 0)"<<endl
69	<<endl;
70	////<<" -A <log10(S_agn in Jy at 1.4 GHz)>,sigma,powlaw :"<<endl
71	////<<" AGN mean and sigma gaussian equiv. distrib. for solid angle of centeral pixel"<<endl
72	////<<" powlaw: apply S_agn evolution as (Nu/1.4)^powlaw"<<endl
73	}
74
75	int main(int narg,char *arg[])
76	{
77	SophyaInit();
78	InitTim();
79
80	//-----------------------------------------------------------------
81	// *** Survey definition
82	long nx=360, ny=-1, nz=64; double dx=1., dy=-1., dz=-1.;
83	//long nx=1000, ny=-1, nz=128; double dx=3., dy=-1., dz=6.;
84	//long nx=1200, ny=-1, nz=128; double dx=1., dy=-1., dz=3;
85
86	// *** Cosmography definition (WMAP)
87	unsigned short flat = 0;
88	double ob0 = 0.0444356;
89	double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
90	double zref = 0.5;
91	double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
92
93	// *** Spectrum and variance definition
94	double ns = 1., as = 1.;
95	double R=8./h100, Rg=R/sqrt(5.);
96	double sigmaR = 1.;
97
98	double kmin=1e-5,kmax=1000.;
99	int npt = 10000;
100	double lkmin=log10(kmin), lkmax=log10(kmax);
101	double eps=1.e-3;
102
103	// *** Schechter mass function definition
104	double h75 = h100 / 0.75;
105	double nstar = 0.006*pow(h75,3.);
106	double mstar = pow(10.,9.8); // MSol
107	double alpha = -1.37;
108
109	double schmin=1.e7, schmax=1.e13;
110	int schnpt = -100;
111	bool use_schmassdist = false;
112	long naleagal = 100000;
113	bool recompute_schmassdist = true;
114	string schmassdistfile = "";
115	int no_poisson_type = 0;
116
117	double scalemass = -1.;
118
119	// *** Niveau de bruit
120	double snoise= 0.; // en equivalent MSol
121	int noise_evol = 0;
122
123	//// *** AGN
124	////bool do_agn = false;
125	////double lfjy_agn=-99., lsigma_agn=0.; // en Jy
126	////double powlaw_agn = 0.;
127
128	// *** type de generation
129	bool computefourier0=false;
130	int use_growth_factor = 0;
131	unsigned short nthread=0;
132	int filter_by_pixel = 1;
133
134	// *** What to do
135	bool comp2dspec = false;
136	bool wfits = false;
137	bool wppf = false;
138	bool wslice = false;
139	bool compvarreal = false;
140	unsigned short whattowrt[5] = {1,1,1,1,1};
141	string rootnameout = "cmvobserv3d";
142
143	unsigned long ntnent = 10000; // 0 = do not fill NTuple
144
145	// --- Decodage arguments
146	if(narg>0) {
147	cout<<"\n--- Arguments: "<<endl;
148	for(int i=0;i<narg;i++) cout<<arg[i]<<" ";
149	cout<<endl;
150	}
151	system("date -u");
152
153	// --- Choix du generateur aleatoire (a faire ICI imperativement avant AutoInitRand)
154	FMTRandGen *RandGen = new FMTRandGen;
155	RandGen->SelectGaussianAlgo(C_Gaussian_RandLibSNorm);
156	RandGen->SelectPoissonAlgo(C_Poisson_Ahrens);
157	RandomGeneratorInterface::SetGlobalRandGenP(RandGen);
158
159	// --- Decodage des arguments
160	char c;
161	while((c = getopt(narg,arg,"ha0PWSV2U:G:F:x:y:z:s:Z:M:A:T:N:Q:R:8:O:o:")) != -1) {
162	int nth = 0;
163	switch (c) {
164	case 'a' :
165	AutoInitRand(5);
166	break;
167	case '0' :
168	computefourier0 = true;
169	break;
170	case 'G' :
171	sscanf(optarg,"%d",&use_growth_factor);
172	break;
173	case 'U' :
174	sscanf(optarg,"%d",&no_poisson_type);
175	break;
176	case 'F' :
177	sscanf(optarg,"%d",&filter_by_pixel);
178	break;
179	case 'x' :
180	sscanf(optarg,"%ld,%lf",&nx,&dx);
181	break;
182	case 'y' :
183	sscanf(optarg,"%ld,%lf",&ny,&dy);
184	break;
185	case 'z' :
186	sscanf(optarg,"%ld,%lf",&nz,&dz);
187	break;
188	case 's' :
189	sscanf(optarg,"%lf,%d",&snoise,&noise_evol);
190	break;
191	case 'Z' :
192	sscanf(optarg,"%lf",&zref);
193	break;
194	case '2' :
195	comp2dspec = true;
196	break;
197	case 'N' :
198	sscanf(optarg,"%lf",&scalemass);
199	if(scalemass==0.) scalemass = 1.;
200	break;
201	case 'M' :
202	sscanf(optarg,"%lf,%lf,%d",&schmin,&schmax,&schnpt);
203	break;
204	case 'Q' :
205	use_schmassdist = true;
206	sscanf(optarg,"%ld",&naleagal);
207	break;
208	case 'R' :
209	schmassdistfile = optarg;
210	break;
211	//// case 'A' :
212	////do_agn = true;
213	////sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
214	////break;
215	case '8' :
216	sscanf(optarg,"%lf,%lf",&sigmaR,&R);
217	break;
218	case 'V' :
219	compvarreal = true;
220	break;
221	case 'W' :
222	wfits = true;
223	break;
224	case 'P' :
225	wppf = true;
226	break;
227	case 'O' :
228	sscanf(optarg,"%hu,%hu,%hu,%hu,%hu"
229	,&whattowrt[0],&whattowrt[1],&whattowrt[2],&whattowrt[3],&whattowrt[4]);
230	break;
231	case 'S' :
232	wslice = true;
233	break;
234	case 'o' :
235	rootnameout = optarg;
236	break;
237	case 'T' :
238	sscanf(optarg,"%d",&nth);
239	nthread = (nth<1)? 0: nth;
240	break;
241	case 'h' :
242	default :
243	usage(); return -1;
244	}
245	}
246
247	double lschmin=log10(schmin), lschmax=log10(schmax);
248	if(schnpt<=0) { // alors c'est un nombre de points par decade
249	schnpt = long( (-schnpt)*(lschmax-lschmin+1.) + 0.5 );
250	if(schnpt<=2) schnpt = 1000;
251	}
252	if(naleagal<=2) naleagal = 100000;
253
254	cout<<"zref="<<zref<<endl;
255	cout<<"nx="<<nx<<" dx="<<dx<<" ny="<<ny<<" dy="<<dy<<" nz="<<nz<<" dz="<<dz<<endl;
256	cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<") Mpc^-1"
257	<<", npt="<<npt<<endl;
258	cout<<"Filter by pixel = "<<filter_by_pixel<<endl;
259	cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
260	cout<<"Use_growth_factor = "<<use_growth_factor<<endl;
261	cout<<"nstar= "<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
262	cout<<"schmin="<<schmin<<" ("<<lschmin
263	<<"), schmax="<<schmax<<" ("<<lschmax<<") Msol"
264	<<", schnpt="<<schnpt<<endl;
265	if(no_poisson_type!=0) cout<<"No poisson fluctuation, direct conversion to HI mass, typ="
266	<<no_poisson_type<<endl;
267	cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<noise_evol<<endl;
268	cout<<"scalemass="<<scalemass<<endl;
269	////if(do_agn)
270	//// cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
271	//// <<" , powlaw="<<powlaw_agn<<endl;
272	cout<<"wfits="<<wfits<<" wppf="<<wppf<<" wslice="<<wslice<<" what?="
273	<<whattowrt[0]<<","<<whattowrt[1]<<","<<whattowrt[2]<<","<<whattowrt[3]<<","
274	<<whattowrt[4]<<endl;
275	cout<<"rootnameout="<<rootnameout<<endl;
276	ShowRandom();
277	cout<<" First random is: "<<drand01()<<endl;
278
279	string tagobs = rootnameout + ".ppf";
280	POutPersist posobs(tagobs);
281
282	//-----------------------------------------------------------------
283	cout<<endl<<"\n--- Create Cosmology"<<endl;
284
285	CosmoCalc univ(flat,true,zref+1.);
286	univ.SetInteg(perc,dzinc,dzmax,glorder);
287	univ.SetDynParam(h100,om0,or0,ol0,w0);
288	univ.PrtInteg();
289	univ.Print();
290	double loscomref = univ.Dloscom(zref);
291	cout<<"\nzref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
292	univ.Print(zref);
293
294	//-----------------------------------------------------------------
295	cout<<endl<<"\n--- Create Spectrum"<<endl;
296
297	InitialSpectrum pkini(ns,as);
298
299	TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
300	//tf.SetNoOscEnv(2);
301
302	GrowthFactor growth(om0,ol0);
303	// GrowthFactor growth(1.,0.); // D(z) = 1/(1+z)
304	double growth_at_z = growth(zref);
305	cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
306
307	PkSpectrum0 pk0(pkini,tf);
308
309	PkSpectrumZ pkz(pk0,growth,zref);
310
311	//-----------------------------------------------------------------
312	pkz.SetZ(0.);
313	cout<<endl<<"\n--- Compute variance for top-hat R="<<R
314	<<" at z="<<pkz.GetZ()<<endl;
315	VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
316	double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
317	double pkmax_th = varpk_th(kfind_th);
318	cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
319	double k1=kmin, k2=kmax;
320	int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
321	cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
322	<<k2<<" ("<<log10(k2)<<")"<<endl;
323
324	double ldlk = (log10(k2)-log10(k1))/npt;
325	varpk_th.SetInteg(0.01,ldlk,-1.,4);
326	double sr2 = varpk_th.Variance(k1,k2);
327	cout<<"varpk_th="<<sr2<<" -> sigma="<<sqrt(sr2)<<endl;
328
329	double normpkz = sigmaR*sigmaR/sr2;
330	pkz.SetScale(normpkz);
331	cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
332
333	{
334	Histo hpkz0(lkmin,lkmax,npt); hpkz0.ReCenterBin();
335	FuncToHisto(pkz,hpkz0,true);
336	tagobs = "hpkz0"; posobs.PutObject(hpkz0,tagobs);
337	}
338
339	pkz.SetZ(zref);
340
341	{
342	Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
343	FuncToHisto(pkz,hpkz,true);
344	tagobs = "hpkz"; posobs.PutObject(hpkz,tagobs);
345	}
346
347	//-----------------------------------------------------------------
348	cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
349	VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
350
351	double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
352	double pkmax_int = varpk_int(kfind_int);
353	cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
354	double k1int=kmin, k2int=kmax;
355	int rcint = varpk_int.FindLimits(pkmax_int/1.e4,k1int,k2int,eps);
356	cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
357	<<k2int<<" ("<<log10(k2int)<<")"<<endl;
358
359	double ldlkint = (log10(k2int)-log10(k1int))/npt;
360	varpk_int.SetInteg(0.01,ldlkint,-1.,4);
361	double sr2int = varpk_int.Variance(k1int,k2int);
362	cout<<"varpk_int="<<sr2int<<" -> sigma="<<sqrt(sr2int)<<endl;
363
364	//-----------------------------------------------------------------
365	cout<<endl<<"\n--- Create mass function, compute number/mass density, init mass trials"<<endl;
366
367	Schechter sch(nstar,mstar,alpha);
368	sch.Print();
369
370	sch.SetOutValue(0);
371	cout<<"sch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
372	double ngal_by_mpc3 = sch.Integrate(schmin,schmax,schnpt);
373	cout<<"Galaxy density number = "<<ngal_by_mpc3<<" /Mpc^3 between limits"<<endl;
374
375	sch.SetOutValue(1);
376	cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
377	double mass_by_mpc3 = sch.Integrate(schmin,schmax,schnpt);
378	cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl;
379	cout<<"Omega_HI at z=0 is "<<mass_by_mpc3/(univ.Rhoc(0.)*GCm3toMsolMpc3_Cst)<<endl
380	<<" at z="<<zref<<" is "<<mass_by_mpc3/(univ.Rhoc(zref)*GCm3toMsolMpc3_Cst)<<endl;
381
382	SchechterMassDist schmdist(sch,schmin,schmax,schnpt);
383	if(use_schmassdist && schmassdistfile.size()>0) {
384	cout<<"\nSchechterMassDist read from "<<schmassdistfile<<endl;
385	schmdist.ReadPPF(schmassdistfile);
386	bool compsch = IsCompatible(sch,schmdist.GetSchechter());
387	double m1,m2; schmdist.GetMassLim(m1,m2);
388	if( !compsch \|\| fabs(m1-schmin)>1e-4schmin \|\| fabs(m2-schmax)>1e-4schmax) {
389	cout<<"FATAL_ERROR: INCONSISTENT SchechterMassDist file / Schechter or limits"<<endl;
390	cout<<"Schechter:"<<endl; sch.Print();
391	cout<<"schmin="<<schmin<<" schmax="<<schmax<<endl;
392	cout<<"SchechterMassDist:"<<endl; schmdist.Print();
393	return -10;
394	}
395	recompute_schmassdist = false;
396	}
397	schmdist.Print();
398	Histo hmdndm = schmdist.GetHmDnDm();
399	FunRan tirhmdndm = schmdist.GetTmDnDm();
400	{
401	tagobs = "hmdndm"; posobs.PutObject(hmdndm,tagobs);
402	Histo hdum1(tirhmdndm);
403	tagobs = "tirhmdndm"; posobs.PutObject(hdum1,tagobs);
404	}
405
406	PrtTim(">>>> End of definition");
407
408	//-----------------------------------------------------------------
409	// FFTW3 (p26): faster if sizes 2^a 3^b 5^c 7^d 11^e 13^f with e+f=0 ou 1
410	cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
411
412	GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
413	fluct3d.SetObservator(zref,-nz/2.);
414	fluct3d.SetCosmology(univ);
415	fluct3d.SetGrowthFactor(growth);
416	fluct3d.LosComRedshift(0.001,-1);
417	//TArray< complex<GEN3D_TYPE> >& pkgen = fluct3d.GetComplexArray();
418	//TArray<GEN3D_TYPE>& rgen = fluct3d.GetRealArray();
419	cout<<endl; fluct3d.Print();
420	cout<<"\nMean number of galaxies per pixel = "<<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
421	double mass_by_pixel = mass_by_mpc3 * fluct3d.GetDVol();
422	cout<<"Mean mass per pixel = "<<mass_by_pixel<<endl;
423
424	double dkmin = fluct3d.GetKincMin();
425	double knyqmax = fluct3d.GetKmax();
426	long nherr = long(knyqmax/dkmin+0.5);
427	cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
428
429	double dktmin = fluct3d.GetKTincMin();
430	double ktnyqmax = fluct3d.GetKTmax();
431	long nherrt = long(ktnyqmax/dktmin+0.5);
432	double dkzmin = fluct3d.GetKinc()[2];
433	double kznyqmax = fluct3d.GetKnyq()[2];
434	long nherrz = long(kznyqmax/dkzmin+0.5);
435	cout<<"For Histo2DErr: d="<<dktmin<<","<<dkzmin
436	<<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
437
438	//-----------------------------------------------------------------
439	cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
440	// En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
441	// sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote ksqrt(2)): Vc = (ksqrt(2))^3
442	// Vc/Vs = 0.675 -> keff = kmax * (0.675)^(1/3) = kmax * 0.877
443	double knyqmax_mod = 0.877*knyqmax;
444	ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
445	varpk_int.SetInteg(0.01,ldlkint,-1.,4);
446	double sr2int_kmax = varpk_int.Variance(k1int,knyqmax_mod);
447	cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<" -> sigma="<<sqrt(sr2int_kmax)<<endl;
448
449	PrtTim(">>>> End Initialisation de GeneFluct3D");
450
451	//-----------------------------------------------------------------
452	cout<<"\n--- Computing a realization in Fourier space"<<endl;
453	if(use_growth_factor>0) pkz.SetZ(0.); else pkz.SetZ(zref);
454	cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
455	if(computefourier0) fluct3d.ComputeFourier0(pkz);
456	else fluct3d.ComputeFourier(pkz);
457	fluct3d.NTupleCheck(posobs,string("ntpkgen"),ntnent);
458	PrtTim(">>>> End Computing a realization in Fourier space");
459
460	cout<<"\n--- Checking realization spectra"<<endl;
461	HistoErr hpkgen(0.,knyqmax,nherr);
462	hpkgen.ReCenterBin(); hpkgen.Zero();
463	hpkgen.Show();
464	fluct3d.ComputeSpectrum(hpkgen);
465	{
466	tagobs = "hpkgen"; posobs.PutObject(hpkgen,tagobs);
467	}
468	PrtTim(">>>> End Checking realization spectra");
469
470	if(comp2dspec) {
471	cout<<"\n--- Checking realization 2D spectra"<<endl;
472	Histo2DErr hpkgen2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
473	hpkgen2.ReCenterBin(); hpkgen2.Zero();
474	hpkgen2.Show();
475	fluct3d.ComputeSpectrum2D(hpkgen2);
476	{
477	tagobs = "hpkgen2"; posobs.PutObject(hpkgen2,tagobs);
478	}
479	PrtTim(">>>> End Checking realization 2D spectra");
480	}
481
482	if(filter_by_pixel!=0) {
483	cout<<"\n--- Computing convolution by pixel shape"<<endl;
484	fluct3d.FilterByPixel();
485	fluct3d.NTupleCheck(posobs,string("ntpkgenf"),ntnent);
486	PrtTim(">>>> End Computing convolution by pixel shape");
487
488	cout<<"\n--- Checking realization spectra after pixel shape convol."<<endl;
489	HistoErr hpkgenfb(0.,knyqmax,nherr);
490	hpkgenfb.ReCenterBin(); hpkgenfb.Zero();
491	hpkgenfb.Show();
492	fluct3d.ComputeSpectrum(hpkgenfb);
493	{
494	tagobs = "hpkgenfb"; posobs.PutObject(hpkgenfb,tagobs);
495	}
496	PrtTim(">>>> End Checking realization spectra");
497
498	cout<<"\n--- Checking realization spectra after pixel shape convol. with pixel correc."<<endl;
499	HistoErr hpkgenf(hpkgenfb); hpkgenf.Zero();
500	fluct3d.ComputeSpectrum(hpkgenf,0.,filter_by_pixel);
501	{
502	tagobs = "hpkgenf"; posobs.PutObject(hpkgenf,tagobs);
503	}
504	PrtTim(">>>> End Checking realization spectra with pixel correc.");
505
506	if(comp2dspec) {
507	cout<<"\n--- Checking realization 2D spectra after pixel shape convol."<<endl;
508	Histo2DErr hpkgenfb2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
509	hpkgenfb2.ReCenterBin(); hpkgenfb2.Zero();
510	hpkgenfb2.Show();
511	fluct3d.ComputeSpectrum2D(hpkgenfb2);
512	{
513	tagobs = "hpkgenfb2"; posobs.PutObject(hpkgenfb2,tagobs);
514	}
515	PrtTim(">>>> End Checking realization 2D spectra");
516
517	cout<<"\n--- Checking realization 2D spectra after pixel shape convol. with pixel correc."<<endl;
518	Histo2DErr hpkgenf2(hpkgenfb2); hpkgenf2.Zero();
519	fluct3d.ComputeSpectrum2D(hpkgenf2,0.,filter_by_pixel);
520	{
521	tagobs = "hpkgenf2"; posobs.PutObject(hpkgenf2,tagobs);
522	}
523	PrtTim(">>>> End Checking realization 2D spectra with pixel correc.");
524	}
525	}
526
527	if(whattowrt[0]==1) {
528	if(wfits) {
529	tagobs = "!" + rootnameout + "_k0.fits";
530	fluct3d.WriteFits(tagobs);
531	PrtTim(">>>> End WriteFits");
532	}
533	if(wppf) {
534	tagobs = rootnameout + "_k0.ppf";
535	fluct3d.WritePPF(tagobs,false);
536	PrtTim(">>>> End WritePPF");
537	}
538	}
539
540	//-----------------------------------------------------------------
541	cout<<"\n--- Computing a realization in real space"<<endl;
542	fluct3d.ComputeReal();
543	double rmin,rmax; fluct3d.MinMax(rmin,rmax);
544	cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
545	fluct3d.NTupleCheck(posobs,string("ntreal"),ntnent);
546	PrtTim(">>>> End Computing a realization in real space");
547
548	if(use_growth_factor>0) {
549	cout<<"\n--- Apply Growth factor"<<endl;
550	cout<<"...D(z=0)="<<growth(0.)<<" D(z="<<zref<<")="<<growth(zref)<<endl;
551	fluct3d.ApplyGrowthFactor(use_growth_factor);
552	fluct3d.MinMax(rmin,rmax);
553	cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
554	fluct3d.NTupleCheck(posobs,string("ntgrow"),ntnent);
555	PrtTim(">>>> End Applying growth factor");
556	}
557
558	int_8 nm;
559	double rmref,rs2ref;
560	cout<<"\n--- Computing reference variance in real space"<<endl;
561	nm = fluct3d.MeanSigma2(rmref,rs2ref);
562	cout<<" rs2ref= "<<rs2ref<<" , rmref="<<rmref<<" ("<<nm<<")"<<endl;
563	PrtTim(">>>> End Computing reference variance in real space");
564
565	if(whattowrt[1]==1) {
566	if(wfits) {
567	tagobs = "!" + rootnameout + "_r0.fits";
568	fluct3d.WriteFits(tagobs);
569	PrtTim(">>>> End WriteFits");
570	}
571	if(wppf) {
572	tagobs = rootnameout + "_r0.ppf";
573	fluct3d.WritePPF(tagobs,true);
574	PrtTim(">>>> End WritePPF");
575	}
576	if(wslice) {
577	tagobs = rootnameout + "_s_r0.ppf";
578	fluct3d.WriteSlicePPF(tagobs);
579	PrtTim(">>>> End WriteSlicePPF");
580	}
581	}
582
583	cout<<"\n--- Check mean and variance in real space"<<endl;
584	fluct3d.NumberOfBad(-1.,1e+200);
585	double rm,rs2;
586	nm = fluct3d.MeanSigma2(rm,rs2);
587	PrtTim(">>>> End Check mean and variance in real space");
588
589	if(compvarreal) {
590	cout<<"\n--- Check variance sigmaR in real space"<<endl;
591	double varr;
592	fluct3d.VarianceFrReal(R,varr);
593	cout<<"...Computed variance = "<<varr
594	<<" , Theorical variance at (z=0) = "<<pow(sigmaR,2.)
595	<<" , at (z="<<zref<<") = "<<pow(sigmaR*growth_at_z,2.)<<endl;
596	PrtTim(">>>> End Check variance sigmaR in real space");
597	}
598
599	//-----------------------------------------------------------------
600	if(no_poisson_type!=0) {
601	cout<<"\n--- Converting !!!DIRECTLY!!! mass into HI mass: mass per pixel ="
602	<<mass_by_pixel<<endl;
603	if(no_poisson_type<0) {
604	cout<<"!!!WARNING!!! dRho/Rho is just multiplied by "<<mass_by_pixel<<endl;
605	fluct3d.ScaleOffset(mass_by_pixel);
606	} else {
607	rm = fluct3d.TurnFluct2MeanNumber(mass_by_mpc3); // ici on doit donner Msol/Mpc^3
608	}
609	fluct3d.NTupleCheck(posobs,string("ntmhi"),ntnent);
610	} else {
611	cout<<"\n--- Converting mass into galaxy number: gal per pixel ="
612	<<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
613	rm = fluct3d.TurnFluct2MeanNumber(ngal_by_mpc3); // ici on doit donner Ngal/Mpc^3
614	fluct3d.NTupleCheck(posobs,string("ntmeang"),ntnent);
615	}
616	nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
617	nm = fluct3d.MeanSigma2(rm,rs2);
618	if(rm>0.) cout<<"normalised sigma(dM/M) is "<<sqrt(rs2)/rm<<endl;
619	PrtTim(">>>> End Converting mass into galaxy number or mass");
620
621	if( no_poisson_type==0 ) {
622
623	cout<<"\n--- Apply poisson on galaxy number"<<endl;
624	fluct3d.ApplyPoisson();
625	nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
626	nm = fluct3d.MeanSigma2(rm,rs2);
627	double xgalmin,xgalmax; fluct3d.MinMax(xgalmin,xgalmax,0.1,1.e50);
628	fluct3d.NTupleCheck(posobs,string("ntpois"),ntnent);
629	PrtTim(">>>> End Apply poisson on galaxy number");
630	if(wslice) {
631	tagobs = rootnameout + "_s_rn.ppf";
632	fluct3d.WriteSlicePPF(tagobs);
633	PrtTim(">>>> End WriteSlicePPF");
634	}
635
636	cout<<"\n--- Convert Galaxy number to HI mass"<<endl;
637	double mhi = 0.;
638	if(use_schmassdist) {
639	if(recompute_schmassdist) {
640	int ngalmax = int(xgalmax+0.5);
641	schmdist.SetNgalLim(ngalmax,1,naleagal);
642	PrtTim(">>>> End creating tabulated histograms for trials");
643	}
644	mhi = fluct3d.TurnNGal2MassQuick(schmdist);
645	schmdist.PrintStatus();
646	} else {
647	mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
648	}
649	cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
650	nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
651	cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
652	nm = fluct3d.MeanSigma2(rm,rs2);
653	fluct3d.NTupleCheck(posobs,string("ntmhi"),ntnent);
654	PrtTim(">>>> End Convert Galaxy number to HI mass");
655
656	}
657
658	//-----------------------------------------------------------------
659	if(whattowrt[2]==1) {
660	if(wfits) {
661	tagobs = "!" + rootnameout + "_r.fits";
662	fluct3d.WriteFits(tagobs);
663	PrtTim(">>>> End WriteFits");
664	}
665	if(wppf) {
666	tagobs = rootnameout + "_r.ppf";
667	fluct3d.WritePPF(tagobs,true);
668	PrtTim(">>>> End WritePPF");
669	}
670	if(wslice) {
671	tagobs = rootnameout + "_s_r.ppf";
672	fluct3d.WriteSlicePPF(tagobs);
673	PrtTim(">>>> End WriteSlicePPF");
674	}
675	}
676
677	//-----------------------------------------------------------------
678	////if(do_agn) {
679	//// cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
680	//// <<" , powlaw="<<powlaw_agn<<endl;
681	//// fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
682	//// nm = fluct3d.MeanSigma2(rm,rs2);
683	//// fluct3d.NTupleCheck(posobs,string("ntagn"),ntnent);
684	//// PrtTim(">>>> End Add AGN");
685	////}
686
687	//-----------------------------------------------------------------
688	double snoisesave = 0.;
689	if(snoise>0.) {
690	cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<", evolution="<<noise_evol<<endl;
691	fluct3d.AddNoise2Real(snoise,noise_evol);
692	snoisesave = snoise;
693	nm = fluct3d.MeanSigma2(rm,rs2);
694	fluct3d.NTupleCheck(posobs,string("ntnois"),ntnent);
695	PrtTim(">>>> End Add noise");
696	}
697
698	//-----------------------------------------------------------------
699	if(scalemass!=0. && scalemass!=1.) { // Si scalemass==0 pas de normalisation
700	if(scalemass<0.) scalemass = 1. / (-scalemass * mass_by_pixel);
701	cout<<"\n--- Scale cube scale="<<scalemass<<endl;
702	fluct3d.ScaleOffset(scalemass);
703	snoisesave *= scalemass;
704	nm = fluct3d.MeanSigma2(rm,rs2);
705	PrtTim(">>>> End Scale cube");
706	}
707	fluct3d.NTupleCheck(posobs,string("ntfin"),ntnent);
708
709	//-----------------------------------------------------------------
710	if(whattowrt[3]==1) {
711	if(wfits) {
712	tagobs = "!" + rootnameout + "_rf.fits";
713	fluct3d.WriteFits(tagobs);
714	PrtTim(">>>> End WriteFits");
715	}
716	if(wppf) {
717	tagobs = rootnameout + "_rf.ppf";
718	fluct3d.WritePPF(tagobs,true);
719	PrtTim(">>>> End WritePPF");
720	}
721	if(wslice) {
722	tagobs = rootnameout + "_s_rf.ppf";
723	fluct3d.WriteSlicePPF(tagobs);
724	PrtTim(">>>> End WriteSlicePPF");
725	}
726	}
727
728	//-----------------------------------------------------------------
729	// -- NE PAS FAIRE CA SI ON VEUT CONTINUER LA SIMULATION -> d_rho/rho ecrase
730
731	cout<<endl<<"\n--- ReComputing spectrum from real space"<<endl;
732	fluct3d.ReComputeFourier();
733	fluct3d.NTupleCheck(posobs,string("ntpkrec"),ntnent);
734	PrtTim(">>>> End ReComputing spectrum");
735
736	if(whattowrt[4]==1) {
737	if(wfits) {
738	tagobs = "!" + rootnameout + "_k.fits";
739	fluct3d.WriteFits(tagobs);
740	PrtTim(">>>> End WriteFits");
741	}
742	if(wppf) {
743	tagobs = rootnameout + "_k.ppf";
744	fluct3d.WritePPF(tagobs,false);
745	PrtTim(">>>> End WritePPF");
746	}
747	}
748
749	cout<<endl<<"\n--- Computing final spectrum"<<endl;
750	HistoErr hpkrecb(0.,knyqmax,nherr); hpkrecb.Zero();
751	hpkrecb.ReCenterBin();
752	hpkrecb.Show();
753	fluct3d.ComputeSpectrum(hpkrecb);
754	{
755	tagobs = "hpkrecb"; posobs.PutObject(hpkrecb,tagobs);
756	}
757	PrtTim(">>>> End Computing final spectrum");
758
759	cout<<endl<<"\n--- Computing final spectrum with pixel deconv."<<endl;
760	HistoErr hpkrec(hpkrecb); hpkrec.Zero();
761	fluct3d.ComputeSpectrum(hpkrec,snoisesave,filter_by_pixel);
762	{
763	tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
764	}
765	PrtTim(">>>> End Computing final spectrum with pixel deconv.");
766
767	if(comp2dspec) {
768	cout<<"\n--- Computing final 2D spectrum"<<endl;
769	Histo2DErr hpkrecb2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
770	hpkrecb2.ReCenterBin(); hpkrecb2.Zero();
771	hpkrecb2.Show();
772	fluct3d.ComputeSpectrum2D(hpkrecb2);
773	{
774	tagobs = "hpkrecb2"; posobs.PutObject(hpkrecb2,tagobs);
775	}
776	PrtTim(">>>> End Computing final 2D spectrum");
777
778	cout<<"\n--- Computing final 2D spectrum with pixel deconv."<<endl;
779	Histo2DErr hpkrec2(hpkrecb2); hpkrec2.Zero();
780	fluct3d.ComputeSpectrum2D(hpkrec2,snoisesave,filter_by_pixel);
781	{
782	tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs);
783	}
784	PrtTim(">>>> End Computing final 2D spectrum with pixel deconv.");
785
786	}
787
788	//-----------------------------------------------------------------
789	delete RandGen;
790	PrtTim(">>>> End Of Job");
791
792	return 0;
793	}
794

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

Download in other formats: