Changeset 3329 in Sophya for trunk/Cosmo/SimLSS/cmvobserv3d.cc

Timestamp:

Sep 27, 2007, 6:21:55 PM (18 years ago)

Author:

cmv

Message:

possibilite de ne pas faire poisson sur Ngal cmv 27/09/2007

File:

• trunk/Cosmo/SimLSS/cmvobserv3d.cc (modified) (10 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -26 +26 @@
      <<"      (default: no, spectrum Pk(z=z_median) for all cube)"<<endl
      <<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl
+     <<" -U : do not poisson fluctuate with Ngal, convert directly to HI mass"<<endl
      <<" -x nx,dx : size along x axis (npix,Mpc)"<<endl
      <<" -y ny,dy : size along y axis (npix,Mpc)"<<endl
@@ -90 +91 @@
  bool recompute_schmassdist = true;
  string schmassdistfile = "";
+ bool no_poisson = false;
 
  // *** Niveau de bruit
@@ -122 +124 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWSV2GF:x:y:z:s:Z:M:A:T:N:Q:R:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWSV2GUF:x:y:z:s:Z:M:A:T:N:Q:R:")) != -1) {
   int nth = 0;
   switch (c) {
@@ -133 +135 @@
   case 'G' :
     use_growth_factor = true;
+    break;
+  case 'U' :
+    no_poisson = true;
     break;
   case 'F' :
@@ -211 +216 @@
      <<"), schmax="<<schmax<<" ("<<lschmax<<") Msol"
      <<", schnpt="<<schnpt<<endl;
+ if(no_poisson) cout<<"No poisson fluctuation, direct conversion to HI mass"<<endl;
  cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<isnoise_evol<<endl;
  cout<<"scalecube="<<scalecube<<", offsetcube="<<offsetcube<<endl;
@@ -385 +391 @@
  PrtTim(">>>> End Computing a realization in Fourier space");
 
- /*
- // CMVTEST
- for(int l=0;l<pkgen.SizeX();l++)
- for(int j=0;j<pkgen.SizeY();j++)
-   for(int i=0;i<pkgen.SizeZ();i++) {
-     double k2 = fluct3d.Kx(i)*fluct3d.Kx(i)
-                +fluct3d.Ky(j)*fluct3d.Ky(j)
-                +fluct3d.Kz(l)*fluct3d.Kz(l);
-     double pk = (k2>0.)? 1./k2: 0.;
-     //pkgen(l,j,i) = ComplexGaussRan(sqrt(pk/2.));
-     //pkgen(l,j,i) = sqrt(pk/2.);
-     pkgen(l,j,i) = 1.;
-   }
-   // CMVTEST
-   */
-
  cout<<"\n--- Checking realization spectra"<<endl;
  HistoErr hpkgen(0.,knyqmax,nherr);
@@ -512 +502 @@
  }
 
- //STOP_HERE_FOR QUICK_DEBUG
- // return -41;
-
  //-----------------------------------------------------------------
  cout<<endl<<"\n--- Converting fluctuations into mass"<<endl;
@@ -521 +508 @@
  PrtTim(">>>> End Converting fluctuations into mass");
 
- cout<<"\n--- Converting mass into galaxy number: gal per pixel ="
-     <<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
- rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
-   nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
-   nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
-  PrtTim(">>>> End Converting mass into galaxy number");
-
- cout<<"\n--- Set negative and null pixels to BAD"<<endl;
- nm = fluct3d.SetToVal(0.,1e+200,-999.);
- PrtTim(">>>> End Set negative pixels to BAD etc...");
-
- cout<<"\n--- Apply poisson on galaxy number"<<endl;
- fluct3d.ApplyPoisson();
-   nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
-   nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
- double xgalmin,xgalmax; fluct3d.MinMax(xgalmin,xgalmax,0.1,1.e50);
- PrtTim(">>>> End Apply poisson on galaxy number");
- if(wslice) {
-   fluct3d.WriteSlicePPF("cmvobserv3d_s_rn.ppf");
-   PrtTim(">>>> End WriteSlicePPF");
- }
-
- cout<<"\n--- Convert Galaxy number to HI mass"<<endl;
- double mhi = 0.;
- if(use_schmassdist) {
-   if(recompute_schmassdist) {
-     int ngalmax = int(xgalmax+0.5);
-     schmdist.SetNgalLim(ngalmax,1,naleagal);
-     PrtTim(">>>> End creating tabulated histograms for trials");
+ if(no_poisson) {
+
+   cout<<"\n--- Converting !!!DIRECTLY!!! mass into HI mass: mass per pixel ="
+       <<mass_by_mpc3*fluct3d.GetDVol()<<endl;
+   rm = fluct3d.TurnMass2HIMass(mass_by_mpc3);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+   PrtTim(">>>> End Converting !!!DIRECTLY!!! mass into HI mass");
+
+ } else {
+
+   cout<<"\n--- Converting mass into galaxy number: gal per pixel ="
+       <<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
+   rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+   PrtTim(">>>> End Converting mass into galaxy number");
+
+   cout<<"\n--- Set negative and null pixels to BAD"<<endl;
+   nm = fluct3d.SetToVal(0.,1e+200,-999.);
+   PrtTim(">>>> End Set negative pixels to BAD etc...");
+
+   cout<<"\n--- Apply poisson on galaxy number"<<endl;
+   fluct3d.ApplyPoisson();
+     nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+   double xgalmin,xgalmax; fluct3d.MinMax(xgalmin,xgalmax,0.1,1.e50);
+   PrtTim(">>>> End Apply poisson on galaxy number");
+   if(wslice) {
+     fluct3d.WriteSlicePPF("cmvobserv3d_s_rn.ppf");
+     PrtTim(">>>> End WriteSlicePPF");
    }
-   mhi = fluct3d.TurnNGal2MassQuick(schmdist);
-   schmdist.PrintStatus();
- } else {
-   mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
- }
- cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
-   nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
-   cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
-   nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
- PrtTim(">>>> End Convert Galaxy number to HI mass");
-
- cout<<"\n--- Set BAD pixels to Zero"<<endl;
- nm = fluct3d.SetToVal(-998.,1e+200,0.);
- PrtTim(">>>> End Set BAD pixels to Zero etc...");
+
+   cout<<"\n--- Convert Galaxy number to HI mass"<<endl;
+   double mhi = 0.;
+   if(use_schmassdist) {
+     if(recompute_schmassdist) {
+       int ngalmax = int(xgalmax+0.5);
+       schmdist.SetNgalLim(ngalmax,1,naleagal);
+       PrtTim(">>>> End creating tabulated histograms for trials");
+     }
+     mhi = fluct3d.TurnNGal2MassQuick(schmdist);
+     schmdist.PrintStatus();
+   } else {
+     mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
+   }
+   cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
+     nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
+     cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+   PrtTim(">>>> End Convert Galaxy number to HI mass");
+
+   cout<<"\n--- Set BAD pixels to Zero"<<endl;
+   nm = fluct3d.SetToVal(-998.,1e+200,0.);
+   PrtTim(">>>> End Set BAD pixels to Zero etc...");
+
+ }
 
  if(wfits) {
@@ -579 +579 @@
  }
 
+ //-----------------------------------------------------------------
  if(do_agn) {
    cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
@@ -587 +588 @@
  }
 
+ //-----------------------------------------------------------------
  if(snoise>0.) {
    cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<", evolution="<<isnoise_evol<<endl;