Changeset 3329 in Sophya

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

Location:

trunk/Cosmo/SimLSS

Files:

• Makefile (modified) (1 diff)
• cmvobserv3d.cc (modified) (10 diffs)
• cmvtvarspec.cc (modified) (1 diff)
• genefluct3d.cc (modified) (1 diff)
• genefluct3d.h (modified) (1 diff)
• geneutils.cc (modified) (4 diffs)
• pkspectrum.cc (modified) (1 diff)

trunk/Cosmo/SimLSS/Makefile

@@ -14 +14 @@
 MYLIB = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvsimbao -lfftw3_threads -lfftw3 -lm
 
-#---- Les programmes utilitaires de calcul de cartes
 PROGS = \
        $(EXE)cmvtuniv $(EXE)cmvtransf $(EXE)cmvtgrowth $(EXE)cmvtstpk \

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;

trunk/Cosmo/SimLSS/cmvtvarspec.cc

@@ -24 +24 @@
  double ob0 = 0.0444356;
  // -- WMAP
- double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ double h100=0.71, om0=0.267804, ol0=0.73;
  // -- ouvert matter only
- //double h100=0.71, om0=0.3, or0=0., ol0=0.,w0=-1.;
+ //double h100=0.71, om0=0.3, ol0=0.;
  // -- plat matter only
- //double h100=0.71, om0=1., or0=0., ol0=0.,w0=-1.;
+ //double h100=0.71, om0=1., ol0=0.;
 
  double ns = 1., as = 1.;

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -1054 +1054 @@
 }
 
+double  GeneFluct3D::TurnMass2HIMass(double m_by_mpc3)
+{
+ if(lp_>0) cout<<"--- TurnMass2HIMass : "<<m_by_mpc3<<" Msol/Mpc^3"<<endl;
+
+ double mall=0., mgood=0.;
+ int_8  nall=0,  ngood=0;
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
+   int_8 ip = IndexR(i,j,l);
+   mall += data_[ip]; nall++;
+   if(data_[ip]>0.) {mgood += data_[ip]; ngood++;}
+ }
+ if(ngood>0) mgood /= (double)ngood;
+ if(nall>0)  mall  /= (double)nall;
+ if(lp_>0) cout<<"...ngood="<<ngood<<" mgood="<<mgood
+               <<",  nall="<<nall<<" mall="<<mall<<endl;
+ if(ngood<=0 || mall<=0.) {
+   cout<<"TurnMass2HIMass_Error: ngood="<<ngood<<" <=0 || mall="<<mall<<" <=0"<<endl;
+   throw RangeCheckError("TurnMass2HIMass_Error: ngood<=0 || mall<=0");
+ }
+
+ // On doit mettre m*Vol masse de HI dans notre survey
+ // On en met uniquement dans les pixels de masse >0.
+ // On MET a zero les pixels <0
+ // On renormalise sur les pixels>0 pour qu'on ait m*Vol masse de HI
+ //   comme on ne prend que les pixels >0, on doit normaliser
+ //   a la moyenne de <1+d_rho/rho> sur ces pixels
+ //   (rappel sur tout les pixels <1+d_rho/rho>=1)
+ // masse de HI a mettre ds 1 px:
+ double dm = m_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
+ if(lp_>0) cout<<"...HI mass density move from "
+               <<m_by_mpc3*Vol_/double(NRtot_)<<" to "<<dm<<" / pixel"<<endl;
+
+ double sum = 0.;
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
+   int_8 ip = IndexR(i,j,l);
+   if(data_[ip]<=0.) data_[ip] = 0.;
+   else {
+     data_[ip] *= dm; // par coherence on multiplie aussi les <=0
+     sum += data_[ip];  // mais on ne les compte pas
+   }
+ }
+
+ if(lp_>0) cout<<sum<<"...mass HI put into survey / "<<m_by_mpc3*Vol_<<endl;
+
+ return sum;
+}
+
 double GeneFluct3D::TurnMass2MeanNumber(double n_by_mpc3)
 // do NOT treate negative or nul values
 {
- if(lp_>0) cout<<"--- TurnMass2MeanNumber ---"<<endl;
+ if(lp_>0) cout<<"--- TurnMass2MeanNumber : "<<n_by_mpc3<<" gal/Mpc^3"<<endl;
 
  double mall=0., mgood=0.;

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -99 +99 @@
 
   void TurnFluct2Mass(void);
+  double TurnMass2HIMass(double m_by_mpc3);
   double TurnMass2MeanNumber(double n_by_mpc3);
   double ApplyPoisson(void);

trunk/Cosmo/SimLSS/geneutils.cc

@@ -64 +64 @@
   : _ymin(0.) , _ymax(0.) , _x(x) , _y(y)
 {
-  int_4 ns = _x.size();
+  uint_4 ns = _x.size();
   if(ns<3 || _y.size()<=0 || ns!=_y.size())
     throw ParmError("InverseFunc::InverseFunc_Error: bad array size");
 
   // Check "x" strictement monotone croissant
-  for(int_4 i=0;i<ns-1;i++)
+  for(uint_4 i=0;i<ns-1;i++)
     if(_x[i+1]<=_x[i]) {
       cout<<"InverseFunc::InverseFunc_Error: _x array not stricly growing"<<endl;
@@ -76 +76 @@
 
   // Check "y" monotone croissant
-  for(int_4 i=0;i<ns-1;i++)
+  for(uint_4 i=0;i<ns-1;i++)
     if(_y[i+1]<_y[i]) {
       cout<<"InverseFunc::InverseFunc_Error: _y array not growing"<<endl;
@@ -134 +134 @@
     // Protection
     if(i1<0) {i1++; i2++; i3++;}
-    if(i3==_y.size()) {i1--; i2--; i3--;}
+    if(i3==(long)_y.size()) {i1--; i2--; i3--;}
     // Interpolation parabolique
     double dy = _y[i3]-_y[i1];
@@ -163 +163 @@
 {
  long n = X.size();
- if(Y.size()<n || n<2)
+ if(n>(long)Y.size() || n<2)
    throw ParmError("InterpTab_Error :  incompatible size between X and Y tables!");
 

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -78 +78 @@
   O0_ = Oc_ + Ob_;
   if(nobaryon_) {O0_ = Oc_; Ob_ = 0.;}
-  if(lp_) cout<<"h100="<<h_<<" Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
-
   double H0 = 100. * h_, h2 = h_*h_;
+  if(lp_) cout<<"h100="<<h_<<" H0="<<H0<<") Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
+
 
   if(tcmb_<0.) tcmb_ = T_CMB_Par;