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

Timestamp:

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

Author:

cmv

Message:

affinage du calcul des meansigma cmv 5/6/2007

File:

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

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;