Context Navigation

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

Visit:

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

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

Download in other formats: