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

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