Changeset 3358 in Sophya for trunk/Cosmo/SimLSS/genefluct3d.cc

Timestamp:

Oct 23, 2007, 12:19:27 PM (18 years ago)

Author:

cmv

Message:

suppression et rename de TurnMass2MeanNumber , cmv 23/10/2007

File:

• trunk/Cosmo/SimLSS/genefluct3d.cc (modified) (3 diffs)

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -1089 +1089 @@
  }
 
- if(lp_>0) cout<<"--- NumberOfBad "<<nbad<<" px out of ]"<<vmin<<","<<vmax<<"["<<endl;
+ if(lp_>0) cout<<"--- NumberOfBad "<<nbad<<" px out of ]"<<vmin<<","<<vmax
+               <<"[ i.e. frac="<<nbad/(double)NRtot_<<endl;
  return nbad;
 }
@@ -1206 +1207 @@
 }
 
-double GeneFluct3D::TurnMass2MeanNumber(double val_by_mpc3)
-{
- if(lp_>0) cout<<"--- TurnMass2MeanNumber : "<<val_by_mpc3<<" quantity (gal or mass)/Mpc^3"<<endl;
-
- double mall=0., mgood=0.;
- int_8  nall=0,  ngood=0;
+double GeneFluct3D::TurnFluct2MeanNumber(double val_by_mpc3)
+{
+ if(lp_>0) cout<<"--- TurnFluct2MeanNumber : "<<val_by_mpc3<<" quantity (gal or mass)/Mpc^3"<<endl;
+
+ // First convert dRho/Rho into 1+dRho/Rho
+ int_8  nball = 0; double sumall = 0., sumall2 = 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 (ou Msol) 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 n*Vol galaxies (ou Msol)
- //   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)
- // nb de gal (ou Msol) a mettre ds 1 px:
- double dn = val_by_mpc3*Vol_/ (mgood/mall) /(double)ngood;
- if(lp_>0) cout<<"...quantity density move from "
-               <<val_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
+   data_[ip] += 1.;
+   nball++; sumall += data_[ip]; sumall2 += data_[ip]*data_[ip];
+ }
+ if(nball>2) {
+   sumall /= (double)nball;
+   sumall2 = sumall2/(double)nball - sumall*sumall;
+   if(lp_>0) cout<<"1+dRho/Rho: mean="<<sumall<<" variance="<<sumall2
+                 <<" -> "<<sqrt(fabs(sumall2))<<endl;
+ }
+
+ // Find contribution for positive pixels
+ int_8  nbpos = 0; double sumpos = 0. , sumpos2 = 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);
+   double v = data_[ip];
+   if(data_[ip]>0.) {nbpos++; sumpos += v; sumpos2 += v*v;}
+ }
+ if(nbpos<1) {
+   cout<<"TurnFluct2MeanNumber_Error: nbpos<1"<<endl;
+   throw RangeCheckError("TurnFluct2MeanNumber_Error: nbpos<1");
+ }
+ sumpos2 = sumpos2/nball - sumpos*sumpos/(nball*nball);
+ if(lp_>0)
+   cout<<"1+dRho/Rho with v<0 set to zero: mean="<<sumpos/nball
+       <<" variance="<<sumpos2<<" -> "<<sqrt(fabs(sumpos2))<<endl;
+   cout<<"Sum of positive values: sumpos="<<sumpos
+       <<" (n(v>0) = "<<nbpos<<" frac(v>0)="<<nbpos/(double)NRtot_<<")"<<endl;
+
+ // - Mettre exactement val_by_mpc3*Vol galaxies (ou Msol) dans notre survey
+ // - Uniquement dans les pixels de masse >0.
+ // - Mise a zero des pixels <0
+ double dn = val_by_mpc3 * Vol_ / sumpos;
+ if(lp_>0) cout<<"...density move from "
+               <<val_by_mpc3*dVol_<<" to "<<dn<<" / pixel"<<endl;
 
  double sum = 0.;
@@ -1248 +1260 @@
  }
 
- if(lp_>0) cout<<sum<<"...quantity put into survey / "<<val_by_mpc3*Vol_<<endl;
+ if(lp_>0) cout<<"...quantity put into survey "<<sum<<" / "<<val_by_mpc3*Vol_<<endl;
 
  return sum;