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

Timestamp:

Oct 1, 2007, 7:10:50 PM (18 years ago)

Author:

cmv

Message:

intro ComputeSpectrum avec soustraction de bruit et deconvolution par la fct de transfert du pixel cmv 01/10/2007

File:

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

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -321 +321 @@
    if(lp_>1 && n>0)
      for(int i=0;i<n;i++)
-       if(i==0 || abs(i-n/2)<2 || i==n-1)
-         cout<<"    "<<i<<"  "<<loscom2zred_[i]<<endl;
+       if(i<2 || abs(i-n/2)<2 || i>=n-2)
+         cout<<"    i="<<i
+             <<"  d="<<loscom2zred_min_+i*(loscom2zred_max_-loscom2zred_min_)/(n-1.)
+             <<" Mpc   z="<<loscom2zred_[i]<<endl;
  }
 
@@ -374 +376 @@
    init_fftw();
    SetObservator(zref,kzref);
+   array_allocated_ = true;
  } catch (PThrowable & exc) {
    cout<<"Exception : "<<(string)typeid(exc).name()
@@ -438 +441 @@
    init_fftw();
    SetObservator(zref,kzref);
+   array_allocated_ = true;
  } catch (PThrowable & exc) {
    cout<<"Exception : "<<(string)typeid(exc).name()
@@ -724 +728 @@
    long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_ *Dx_/2;
-   double pk_x = pixelfilter(kx); 
+   double pk_x = pixelfilter(kx);
    for(long j=0;j<Ny_;j++) {
      long jj = (j>Ny_/2) ? Ny_-j : j;
      double ky = jj*Dky_ *Dy_/2;
-     double pk_y =  pixelfilter(ky); 
+     double pk_y = pixelfilter(ky);
      for(long l=0;l<NCz_;l++) {
        double kz = l*Dkz_ *Dz_/2;
@@ -826 +830 @@
      double ky = jj*Dky_;  ky *= ky;
      for(long l=0;l<NCz_;l++) {
-       double kz = l*Dkz_;  kz *= kz;
-       double k = sqrt(kx+ky+kz);
+       double kz = l*Dkz_;
+       double k = sqrt(kx+ky+kz*kz);
        double pk = MODULE2(T_(l,j,i));
        herr.Add(k,pk);
@@ -862 +866 @@
        double pk = MODULE2(T_(l,j,i));
        herr.Add(kt,kz,pk);
+     }
+   }
+ }
+ herr.ToVariance();
+
+ // renormalize to directly compare to original spectrum
+ double norm = Vol_;
+ herr *= norm;
+
+ return 0;
+}
+
+//-------------------------------------------------------------------
+int GeneFluct3D::ComputeSpectrum(HistoErr& herr,double sigma,bool pixcor)
+// Compute spectrum from "T" and fill HistoErr "herr"
+// AVEC la soustraction du niveau de bruit et la correction par filterpixel()
+// Si on ne fait pas ca, alors on obtient un spectre non-isotrope!
+//
+// T : dans le format standard de GeneFuct3D: T(nz,ny,nx)
+// cad T(kz,ky,kx) avec  0<kz<kz_nyq  -ky_nyq<ky<ky_nyq  -kx_nyq<kx<kx_nyq
+{
+  if(lp_>0) cout<<"--- ComputeSpectrum: sigma="<<sigma<<endl;
+ check_array_alloc();
+
+ if(sigma<=0.) sigma = 0.;
+ double sigma2 = sigma*sigma / (double)NRtot_;
+
+ if(herr.NBins()<0) return -1;
+ herr.Zero();
+
+ TVector<r_8> vfz(NCz_);
+ if(pixcor)  // kz = l*Dkz_
+   for(long l=0;l<NCz_;l++) {vfz(l)=pixelfilter(l*Dkz_ *Dz_/2); vfz(l)*=vfz(l);}
+
+ // Attention a l'ordre
+ for(long i=0;i<Nx_;i++) {
+   long ii = (i>Nx_/2) ? Nx_-i : i;
+   double kx = ii*Dkx_;
+   double fx = (pixcor) ? pixelfilter(kx*Dx_/2): 1.;
+   kx *= kx; fx *= fx;
+   for(long j=0;j<Ny_;j++) {
+     long jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_;
+     double fy = (pixcor) ? pixelfilter(ky*Dy_/2): 1.;
+     ky *= ky; fy *= fy;
+     for(long l=0;l<NCz_;l++) {
+       double kz = l*Dkz_;
+       double k = sqrt(kx+ky+kz*kz);
+       double pk = MODULE2(T_(l,j,i)) - sigma2;
+       double fz = (pixcor) ? vfz(l): 1.;
+       double f = fx*fy*fz;
+       if(f>0.) herr.Add(k,pk/f);
+     }
+   }
+ }
+ herr.ToVariance();
+
+ // renormalize to directly compare to original spectrum
+ double norm = Vol_;
+ herr *= norm;
+
+ return 0;
+}
+
+int GeneFluct3D::ComputeSpectrum2D(Histo2DErr& herr,double sigma,bool pixcor)
+// AVEC la soustraction du niveau de bruit et la correction par filterpixel()
+{
+ if(lp_>0) cout<<"--- ComputeSpectrum2D: sigma="<<sigma<<endl;
+ check_array_alloc();
+
+ if(sigma<=0.) sigma = 0.;
+ double sigma2 = sigma*sigma / (double)NRtot_;
+
+ if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
+ herr.Zero();
+
+ TVector<r_8> vfz(NCz_);
+ if(pixcor)  // kz = l*Dkz_
+   for(long l=0;l<NCz_;l++) {vfz(l)=pixelfilter(l*Dkz_ *Dz_/2); vfz(l)*=vfz(l);}
+
+ // Attention a l'ordre
+ for(long i=0;i<Nx_;i++) {
+   long ii = (i>Nx_/2) ? Nx_-i : i;
+   double kx = ii*Dkx_;
+   double fx = (pixcor) ? pixelfilter(kx*Dx_/2) : 1.;
+   kx *= kx; fx *= fx;
+   for(long j=0;j<Ny_;j++) {
+     long jj = (j>Ny_/2) ? Ny_-j : j;
+     double ky = jj*Dky_;
+     double fy = (pixcor) ? pixelfilter(ky*Dy_/2) : 1.;
+     ky *= ky; fy *= fy;
+     double kt = sqrt(kx+ky);
+     for(long l=0;l<NCz_;l++) {
+       double kz = l*Dkz_;
+       double pk = MODULE2(T_(l,j,i)) - sigma2;
+       double fz = (pixcor) ? vfz(l): 1.;
+       double f = fx*fy*fz;
+       if(f>0.) herr.Add(kt,kz,pk/f);
      }
    }