Changeset 3262 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Jun 5, 2007, 7:52:59 PM (18 years ago)

Author:

cmv

Message:

affinage du calcul des meansigma cmv 5/6/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvobserv3d.cc (modified) (7 diffs)
• genefluct3d.cc (modified) (10 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -100 +100 @@
 
  // --- Decodage arguments
+ if(narg>0) {
+   for(int i=0;i<narg;i++) cout<<arg[i]<<" ";
+   cout<<endl;
+ }
 
  char c;
@@ -414 +418 @@
  if(1) {
    cout<<"\n--- Check mean and variance in real space"<<endl;
-   int_8 nlowone = fluct3d.NumberOfBad(-1.,1e+200);
+   fluct3d.NumberOfBad(-1.,1e+200);
    nm = fluct3d.MeanSigma2(rm,rs2);
-   cout<<"rgen:("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-       <<rs2<<" -> "<<sqrt(rs2)<<",  n(<-1)="<<nlowone<<endl;
    PrtTim(">>>> End Check mean and variance in real space");
  }
@@ -424 +426 @@
    cout<<"\n--- Check variance sigmaR in real space"<<endl;
    double varr;
-   int_8 nvarr = fluct3d.VarianceFrReal(R,varr);
-   cout<<"R="<<R<<" : sigmaR^2="<<varr<<" -> "<<sqrt(varr)<<",   n="<<nvarr<<endl;
+   fluct3d.VarianceFrReal(R,varr);
    PrtTim(">>>> End Check variance sigmaR in real space");
  }
@@ -432 +433 @@
  cout<<endl<<"\n--- Converting fluctuations into mass"<<endl;
  fluct3d.TurnFluct2Mass();
- nm = fluct3d.MeanSigma2(rm,rs2);
- cout<<"1+rgen: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
- cout<<"1+rgen with_neg_a_zero: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+   nm = fluct3d.MeanSigma2(rm,rs2);
  PrtTim(">>>> End Converting fluctuations into mass");
 
  cout<<"\n--- Converting mass into galaxy number"<<endl;
  rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
- cout<<rm<<" galaxies put into survey"<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2,0.);
- cout<<"galaxy mass: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
- cout<<"galaxy mass with_neg_a_zero: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
- PrtTim(">>>> End Converting mass into galaxy number");
+   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.);
- cout<<nm<<" negative in survey set to BAD"<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2,-998.);
- cout<<"galaxy mass: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
  PrtTim(">>>> End Set negative pixels to BAD etc...");
 
  cout<<"\n--- Apply poisson on galaxy number"<<endl;
- nm = fluct3d.ApplyPoisson();
- cout<<nm<<" galaxies into survey after poisson"<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2,-998.);
- cout<<"galaxy number: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
- cout<<"galaxy number with_neg_a_zero: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
+ fluct3d.ApplyPoisson();
+   nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
+   nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
  PrtTim(">>>> End Apply poisson on galaxy number");
 
@@ -483 +465 @@
  double mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
  cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2,-998.);
- cout<<"HI mass: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
- cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
- cout<<"HI mass with_neg_a_zero: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<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.);
- cout<<nm<<" BAD in survey set to zero"<<endl;
- nm = fluct3d.MeanSigma2(rm,rs2);
- cout<<"galaxy: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-     <<rs2<<" -> "<<sqrt(rs2)<<endl;
  PrtTim(">>>> End Set BAD pixels to Zero etc...");
 
@@ -513 +487 @@
        <<" , powlaw="<<powlaw_agn<<endl;
    fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
-   nm = fluct3d.MeanSigma2(rm,rs2);
-   cout<<"HI mass with AGN: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-       <<rs2<<" -> "<<sqrt(rs2)<<endl;
+     nm = fluct3d.MeanSigma2(rm,rs2);
    PrtTim(">>>> End Add AGN");
  }
@@ -522 +494 @@
    cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<endl;
    fluct3d.AddNoise2Real(snoise);
-   nm = fluct3d.MeanSigma2(rm,rs2);
-   cout<<"HI mass with noise: ("<<nm<<") Mean = "<<rm<<", Sigma^2 = "
-       <<rs2<<" -> "<<sqrt(rs2)<<endl;
+     nm = fluct3d.MeanSigma2(rm,rs2);
    PrtTim(">>>> End Add noise");
  }

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -381 +381 @@
    bool found_tag_k = pis.GotoNameTag("pkgen");
    if(found_tag_k) {
-     cout<<"           ...reading spectrun into TArray<complex <r_8> >"<<endl;
+     cout<<"           ...reading spectrum into TArray<complex <r_8> >"<<endl;
      pis >> PPFNameTag("pkgen")  >> T_;
      from_real = false;
@@ -792 +792 @@
 // Recompute MASS variance in spherical top-hat (rayon=R)
 {
- if(lp_>0) cout<<"--- VarianceFrReal ---"<<endl;
+  if(lp_>0) cout<<"--- VarianceFrReal R="<<R<<endl;
  check_array_alloc();
 
@@ -829 +829 @@
  sum /= nsum;
  sum2 = sum2/nsum - sum*sum;
- if(lp_>0) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
+ if(lp_>0) cout<<"...nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
 
  var = sum2/(sum*sum);  // <dM>^2/<M>^2
- if(lp_>0) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
+ if(lp_>0) cout<<"...sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
 
  return nsum;
@@ -851 +851 @@
  }
 
+ if(lp_>0) cout<<"--- NumberOfBad "<<nbad<<" px out of ]"<<vmin<<","<<vmax<<"["<<endl;
  return nbad;
 }
@@ -863 +864 @@
  bool tstval = (vmax>vmin)? true: false;
  if(lp_>0) {
-   cout<<"--- MeanSigma2: ";
-   if(tstval) cout<<"range=]"<<vmin<<","<<vmax<<"[";
+   cout<<"--- MeanSigma2";
+   if(tstval) cout<<"  range=]"<<vmin<<","<<vmax<<"[";
    if(useout) cout<<", useout="<<useout<<" vout="<<vout;
    cout<<endl;
@@ -898 +899 @@
 // cad set to "val0" if in [vmin,vmax] -> vmin and vmax are set to val0
 {
- if(lp_>0) cout<<"--- SetToVal set to="<<val0
-               <<" when in range=["<<vmin<<","<<vmax<<"]"<<endl;
  check_array_alloc();
 
@@ -909 +908 @@
  }
 
+ if(lp_>0) cout<<"--- SetToVal "<<nbad<<" px set to="<<val0
+               <<" because out of range=]"<<vmin<<","<<vmax<<"["<<endl;
  return nbad;
 }
@@ -931 +932 @@
  if(lp_>0) cout<<"--- TurnMass2MeanNumber ---"<<endl;
 
- double m,s2;
- int_8 ngood = MeanSigma2(m,s2,0.,1e+200);
- if(lp_>0) cout<<"...ngood="<<ngood
-               <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<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<<"TurnMass2MeanNumber_Error: ngood="<<ngood<<" <=0 || mall="<<mall<<" <=0"<<endl;
+   throw RangeCheckError("TurnMass2MeanNumber_Error: ngood<=0 || mall<=0");
+ }
 
  // On doit mettre n*Vol galaxies dans notre survey
@@ -943 +955 @@
  //   a la moyenne de <1+d_rho/rho> sur ces pixels
  //   (rappel sur tout les pixels <1+d_rho/rho>=1)
- double dn = n_by_mpc3*Vol_/m /(double)ngood;  // nb de gal a mettre ds 1 px
+ // nb de gal a mettre ds 1 px:
+ double dn = n_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
  if(lp_>0) cout<<"...galaxy density move from "
                <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / 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++) {
@@ -952 +966 @@
    if(data_[ip]>0.) sum += data_[ip];  // mais on ne les compte pas
  }
+
  if(lp_>0) cout<<sum<<"...galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;