Changeset 3330 in Sophya for trunk/Cosmo/SimLSS

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

Location:

trunk/Cosmo/SimLSS

Files:

• cmvconcherr.cc (modified) (1 diff)
• cmvobserv3d.cc (modified) (4 diffs)
• genefluct3d.cc (modified) (6 diffs)
• genefluct3d.h (modified) (3 diffs)
• geneutils.cc (modified) (3 diffs)
• geneutils.h (modified) (1 diff)

trunk/Cosmo/SimLSS/cmvconcherr.cc

@@ -32 +32 @@
 int main(int narg,char *arg[])
 {
- // "hpkgen" "hpkgenf" "hpkrec" "hpkgen2" "hpkgenf2" "hpkrec2"
+ // "hpkgen"  "hpkgenfb"  "hpkgenf"  "hpkrecb"  "hpkrec"
+ // "hpkgen2" "hpkgenfb2" "hpkgenf2" "hpkrecb2" "hpkrec2"
  string nameh = "";
  vector<string> ppfname;

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -419 +419 @@
 
    cout<<"\n--- Checking realization spectra after pixel shape convol."<<endl;
-   HistoErr hpkgenf(0.,knyqmax,nherr);
-   hpkgenf.ReCenterBin(); hpkgenf.Zero();
-   hpkgenf.Show();
-   fluct3d.ComputeSpectrum(hpkgenf);
+   HistoErr hpkgenfb(0.,knyqmax,nherr);
+   hpkgenfb.ReCenterBin(); hpkgenfb.Zero();
+   hpkgenfb.Show();
+   fluct3d.ComputeSpectrum(hpkgenfb);
+   {
+   tagobs = "hpkgenfb"; posobs.PutObject(hpkgenfb,tagobs);
+   }
+   PrtTim(">>>> End Checking realization spectra");
+
+   cout<<"\n--- Checking realization spectra after pixel shape convol. with pixel correc."<<endl;
+   HistoErr hpkgenf(hpkgenfb); hpkgenf.Zero();
+   fluct3d.ComputeSpectrum(hpkgenf,0.,filter_by_pixel);
    {
    tagobs = "hpkgenf"; posobs.PutObject(hpkgenf,tagobs);
    }
-   PrtTim(">>>> End Checking realization spectra");
+   PrtTim(">>>> End Checking realization spectra with pixel correc.");
 
    if(comp2dspec) {
      cout<<"\n--- Checking realization 2D spectra after pixel shape convol."<<endl;
-     Histo2DErr hpkgenf2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
-     hpkgenf2.ReCenterBin(); hpkgenf2.Zero();
-     hpkgenf2.Show();
-     fluct3d.ComputeSpectrum2D(hpkgenf2);
+     Histo2DErr hpkgenfb2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+     hpkgenfb2.ReCenterBin(); hpkgenfb2.Zero();
+     hpkgenfb2.Show();
+     fluct3d.ComputeSpectrum2D(hpkgenfb2);
+     {
+     tagobs = "hpkgenfb2"; posobs.PutObject(hpkgenfb2,tagobs);
+     }
+     PrtTim(">>>> End Checking realization 2D spectra");
+
+     cout<<"\n--- Checking realization 2D spectra after pixel shape convol. with pixel correc."<<endl;
+     Histo2DErr hpkgenf2(hpkgenfb2); hpkgenf2.Zero();
+     fluct3d.ComputeSpectrum2D(hpkgenf2,0.,filter_by_pixel);
      {
      tagobs = "hpkgenf2"; posobs.PutObject(hpkgenf2,tagobs);
      }
-     PrtTim(">>>> End Checking realization 2D spectra");
+     PrtTim(">>>> End Checking realization 2D spectra with pixel correc.");
    }
  }
@@ -589 +605 @@
 
  //-----------------------------------------------------------------
+ double snoisesave = 0.;
  if(snoise>0.) {
    cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<", evolution="<<isnoise_evol<<endl;
    fluct3d.AddNoise2Real(snoise,(isnoise_evol>0? true:false));
+   snoisesave = snoise;
      nm = fluct3d.MeanSigma2(rm,rs2);
    PrtTim(">>>> End Add noise");
@@ -603 +621 @@
    }
    cout<<"...scale="<<scalecube<<" offset="<<offsetcube<<endl;
-   if(scalecube>0.) fluct3d.ScaleOffset(scalecube,offsetcube);
+   if(scalecube>0.) {
+     fluct3d.ScaleOffset(scalecube,offsetcube);
+     snoisesave *= scalecube;
+   }
    PrtTim(">>>> End Scale cube");
  }
@@ -637 +658 @@
 
  cout<<endl<<"\n--- Computing final spectrum"<<endl;
- HistoErr hpkrec(0.,knyqmax,nherr);
- hpkrec.ReCenterBin();
- hpkrec.Show();
- fluct3d.ComputeSpectrum(hpkrec);
+ HistoErr hpkrecb(0.,knyqmax,nherr); hpkrecb.Zero();
+ hpkrecb.ReCenterBin();
+ hpkrecb.Show();
+ fluct3d.ComputeSpectrum(hpkrecb);
+ {
+ tagobs = "hpkrecb"; posobs.PutObject(hpkrecb,tagobs);
+ }
+ PrtTim(">>>> End Computing final spectrum");
+
+ cout<<endl<<"\n--- Computing final spectrum with pixel deconv."<<endl;
+ HistoErr hpkrec(hpkrecb); hpkrec.Zero();
+ fluct3d.ComputeSpectrum(hpkrec,snoisesave,filter_by_pixel);
+ {
  tagobs = "hpkrec"; posobs.PutObject(hpkrec,tagobs);
- PrtTim(">>>> End Computing final spectrum");
+ }
+ PrtTim(">>>> End Computing final spectrum with pixel deconv.");
 
  if(comp2dspec) {
    cout<<"\n--- Computing final 2D spectrum"<<endl;
-   Histo2DErr hpkrec2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
-   hpkrec2.ReCenterBin(); hpkrec2.Zero();
-   hpkrec2.Show();
-   fluct3d.ComputeSpectrum2D(hpkrec2);
+   Histo2DErr hpkrecb2(0.,ktnyqmax,nherrt,0.,kznyqmax,nherrz);
+   hpkrecb2.ReCenterBin(); hpkrecb2.Zero();
+   hpkrecb2.Show();
+   fluct3d.ComputeSpectrum2D(hpkrecb2);
+   {
+   tagobs = "hpkrecb2"; posobs.PutObject(hpkrecb2,tagobs);
+   }
+   PrtTim(">>>> End Computing final 2D spectrum");
+
+   cout<<"\n--- Computing final 2D spectrum with pixel deconv."<<endl;
+   Histo2DErr hpkrec2(hpkrecb2); hpkrec2.Zero();
+   fluct3d.ComputeSpectrum2D(hpkrec2,snoisesave,filter_by_pixel);
    {
    tagobs = "hpkrec2"; posobs.PutObject(hpkrec2,tagobs);
    }
-   PrtTim(">>>> End Computing final 2D spectrum");
+   PrtTim(">>>> End Computing final 2D spectrum with pixel deconv.");
+
  }
 

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);
      }
    }

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -51 +51 @@
   // Pour adressage fdata_[ip][0-1]
   inline int_8 IndexC(long i,long j,long k) {return (int_8)(k+NCz_*(j+Ny_*i));}
-  // - On peut aussi adresser:
-  // TArray< complex<r_8> >& pk = gf3d.GetComplexArray(); pk(i,j,k) = ...;
-  // TArray<r_8>& rgen = gf3d.GetRealArray(); rgen(i,j,k) = ...;
+  // On peut aussi adresser:
+  // TArray< complex<r_8> >& pk = gf3d.GetComplexArray();
+  //    pk(k,j,i) avec k=[0,NCz_[  j=[0,Ny_[   i=[0,Nx_[
+  //    pk[IndexC(i,j,k)]
+  // TArray<r_8>& rgen = gf3d.GetRealArray();
+  //    rgen(k,j,i) avec k=[0,NTz_[  j=[0,Ny_[   i=[0,Nx_[
+  //                mais seul k=[0,Nz_[ est utile
+  //    rgen[IndexR(i,j,k)]
+  // ATTENTION: TArray adresse en memoire a l'envers du C !
+  //            Tarray(n1,n2,n3) == Carray[n3][n2][n1]
 
   vector<long> GetNpix(void) {return N_;}
   int_8 NPix(void) {return NRtot_;}
+  long GetNx(void) {return Nx_;}
+  long GetNy(void) {return Ny_;}
+  long GetNz(void) {return Nz_;}
 
   // Return |K_i| module relative to pixel indices
@@ -83 +93 @@
 
   void ComputeReal(void);
+  void ApplyGrowthFactor(void);
 
   void ReComputeFourier(void);
@@ -88 +99 @@
   int  ComputeSpectrum(HistoErr& herr);
   int  ComputeSpectrum2D(Histo2DErr& herr);
-  void ApplyGrowthFactor(void);
+  int  ComputeSpectrum(HistoErr& herr,double sigma,bool pixcor);
+  int  ComputeSpectrum2D(Histo2DErr& herr,double sigma,bool pixcor);
 
   int_8 VarianceFrReal(double R,double& var);

trunk/Cosmo/SimLSS/geneutils.cc

@@ -55 +55 @@
 //-------------------------------------------------------------------
 // Classe d'inversion d'une fonction STRICTEMENT MONOTONE CROISSANTE
-// - Le vecteur y a "Nin" elements y_i tels que "y_i = f(x_i)"
-// - On a x(i) < x(i+1) et y(i) < y(i+1)
-// - La classe renvoie ymin=y(0) , ymax=y(Nin -1)
-//   et le vecteur x = f^-1(y) de "Nout" elements
-//   Les y_i sont regulierement espaces et ymin et ymax
-//   La re-interpolation inverse est faite par lineaire
+//
+// - On part de deux vecteurs x,y de "Nin" elements tels que "y_i = f[x_i]"
+//   ou la fonction "f" est strictement monotone croissante:
+//          x(i) < x(i+1)   et   y(i) < y(i+1)
+// - Le but de la classe est de remplir un vecteur X de "Nout" elements
+//   tels que:    X_j = f^-1[Y_j]    avec   j=[0,Nout[
+//   avec les Y_j regulierement espaces entre ymin=y(0) , ymax=y(Nin -1)
+//   cad:     X_j = f^-1[ ymin+j*(ymax-ymin)/(Nout-1) ]
+// - La construction du vecteur X est realisee
+//   par interpolation lineaire (ComputeLinear) ou parabolique (ComputeParab)
 InverseFunc::InverseFunc(vector<double>& x,vector<double>& y)
   : _ymin(0.) , _ymax(0.) , _x(x) , _y(y)
@@ -92 +96 @@
 }
 
-int InverseFunc::ComputeLinear(long n,vector<double>& xfcty)
+int InverseFunc::ComputeLinear(long nout,vector<double>& xfcty)
 // Compute table "xfcty" by linear interpolation of "x" versus "y"
-//   on "n" points from "ymin" to "ymax":
-// xfcty[i] = interpolation of function "x" for "ymin+i*(ymax-ymin)/(n-1.)"
-{
-  if(n<3) return -1;
-
-  xfcty.resize(n);
+//   on "nout" points from "ymin" to "ymax":
+// xfcty[i] = interpolation of function "x" for "ymin+i*(ymax-ymin)/(nout-1)"
+{
+  if(nout<3) return -1;
+
+  xfcty.resize(nout);
 
   long i1,i2;
   double x;
-  for(int_4 i=0;i<n;i++) {
-    double y = _ymin + i*(_ymax-_ymin)/(n-1.);
+  for(int_4 i=0;i<nout;i++) {
+    double y = _ymin + i*(_ymax-_ymin)/(nout-1.);
     find_in_y(y,i1,i2);
     double dy = _y[i2]-_y[i1];
@@ -118 +122 @@
 }
 
-int InverseFunc::ComputeParab(long n,vector<double>& xfcty)
-{
-  if(n<3) return -1;
-
-  xfcty.resize(n);
+int InverseFunc::ComputeParab(long nout,vector<double>& xfcty)
+{
+  if(nout<3) return -1;
+
+  xfcty.resize(nout);
 
   long i1,i2,i3;
   double x;
-  for(int_4 i=0;i<n;i++) {
-    double y = _ymin + i*(_ymax-_ymin)/(n-1.);
+  for(int_4 i=0;i<nout;i++) {
+    double y = _ymin + i*(_ymax-_ymin)/(nout-1.);
     find_in_y(y,i1,i2);
     // On cherche le 3ieme point selon la position de y / au 2 premiers

trunk/Cosmo/SimLSS/geneutils.h

@@ -71 +71 @@
     }
 
+  double _ymin,_ymax;
   vector<double>& _x;
   vector<double>& _y;
-  double _ymin,_ymax,_dy;
 };