Changeset 3524 in Sophya for trunk/Cosmo

Timestamp:

Sep 22, 2008, 4:07:02 PM (17 years ago)

Author:

cmv

Message:

intro du remplissage du NTuple de debug cmv 22/09/2008

Location:

trunk/Cosmo/SimLSS

Files:

• cmvobserv3d.cc (modified) (18 diffs)
• genefluct3d.cc (modified) (8 diffs)
• genefluct3d.h (modified) (2 diffs)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -28 +28 @@
      <<"      (default: no, spectrum Pk(z=z_median) for all cube)"<<endl
      <<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl
-     <<" -U typ : do not poisson fluctuate with Ngal, convert directly to HI mass"<<endl
-     <<"          (typ<0 just multiply dRho/Rho by mass_by_pixel)"<<endl
+     <<" -U typ =0 compute <NGal>, poisson fluctuate then convert to HI mass (def)"<<endl
+     <<"        >0 compute directly <HI mass>, do NOT poisson fluctuate with <Ngal>"<<endl
+     <<"        <0 just multiply dRho/Rho by mass_by_pixel (possible negative pixel values)"<<endl
      <<" -x nx,dx : size along x axis (npix,Mpc)"<<endl
      <<" -y ny,dy : size along y axis (npix,Mpc)"<<endl
@@ -100 +101 @@
  bool recompute_schmassdist = true;
  string schmassdistfile = "";
- bool no_poisson = false;
- int no_poisson_type = 1;
+ int no_poisson_type = 0;
 
  double scalemass = -1.;
@@ -126 +126 @@
  bool wslice = false;
  bool compvarreal = false;
+ unsigned long ntnent = 10000;  // 0 = do not fill NTuple
 
  // --- Decodage arguments
@@ -148 +149 @@
     break;
   case 'U' :
-    no_poisson = true;
     sscanf(optarg,"%d",&no_poisson_type);
     break;
@@ -233 +233 @@
      <<"), schmax="<<schmax<<" ("<<lschmax<<") Msol"
      <<", schnpt="<<schnpt<<endl;
- if(no_poisson) cout<<"No poisson fluctuation, direct conversion to HI mass, typ="
-                    <<no_poisson_type<<endl;
+ if(no_poisson_type!=0) cout<<"No poisson fluctuation, direct conversion to HI mass, typ="
+                            <<no_poisson_type<<endl;
  cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<noise_evol<<endl;
  cout<<"scalemass="<<scalemass<<endl;
@@ -419 +419 @@
  if(computefourier0) fluct3d.ComputeFourier0(pkz);
    else              fluct3d.ComputeFourier(pkz);
+ fluct3d.NTupleCheck(posobs,string("ntpkgen"),ntnent);
  PrtTim(">>>> End Computing a realization in Fourier space");
 
@@ -446 +447 @@
    cout<<"\n--- Computing convolution by pixel shape"<<endl;
    fluct3d.FilterByPixel();
+   fluct3d.NTupleCheck(posobs,string("ntpkgenf"),ntnent);
    PrtTim(">>>> End Computing convolution by pixel shape");
 
@@ -501 +503 @@
  double rmin,rmax; fluct3d.MinMax(rmin,rmax);
  cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+ fluct3d.NTupleCheck(posobs,string("ntreal"),ntnent);
  PrtTim(">>>> End Computing a realization in real space");
 
@@ -509 +512 @@
    fluct3d.MinMax(rmin,rmax);
    cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   fluct3d.NTupleCheck(posobs,string("ntgrow"),ntnent);
    PrtTim(">>>> End Applying growth factor");
  }
@@ -549 +553 @@
 
  //-----------------------------------------------------------------
- if(no_poisson) {
+ if(no_poisson_type!=0) {
    cout<<"\n--- Converting !!!DIRECTLY!!! mass into HI mass: mass per pixel ="
        <<mass_by_pixel<<endl;
@@ -558 +562 @@
      rm = fluct3d.TurnFluct2MeanNumber(mass_by_mpc3); // ici on doit donner Msol/Mpc^3
    }
+   fluct3d.NTupleCheck(posobs,string("ntmhi"),ntnent);
  } else {
    cout<<"\n--- Converting mass into galaxy number: gal per pixel ="
        <<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
    rm = fluct3d.TurnFluct2MeanNumber(ngal_by_mpc3); // ici on doit donner Ngal/Mpc^3
+   fluct3d.NTupleCheck(posobs,string("ntmeang"),ntnent);
  }
  nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
@@ -568 +574 @@
  PrtTim(">>>> End Converting mass into galaxy number or mass");
 
- if( !no_poisson ) {
+ if( no_poisson_type==0 ) {
 
    cout<<"\n--- Apply poisson on galaxy number"<<endl;
@@ -575 +581 @@
      nm = fluct3d.MeanSigma2(rm,rs2);
    double xgalmin,xgalmax; fluct3d.MinMax(xgalmin,xgalmax,0.1,1.e50);
+   fluct3d.NTupleCheck(posobs,string("ntpois"),ntnent);
    PrtTim(">>>> End Apply poisson on galaxy number");
    if(wslice) {
@@ -598 +605 @@
      cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
      nm = fluct3d.MeanSigma2(rm,rs2);
+   fluct3d.NTupleCheck(posobs,string("ntmhi"),ntnent);
    PrtTim(">>>> End Convert Galaxy number to HI mass");
 
@@ -622 +630 @@
  ////  fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
  ////    nm = fluct3d.MeanSigma2(rm,rs2);
- ////   PrtTim(">>>> End Add AGN");
+ ////  fluct3d.NTupleCheck(posobs,string("ntagn"),ntnent);
+ ////  PrtTim(">>>> End Add AGN");
  ////}
 
@@ -632 +641 @@
    snoisesave = snoise;
      nm = fluct3d.MeanSigma2(rm,rs2);
+   fluct3d.NTupleCheck(posobs,string("ntnois"),ntnent);
    PrtTim(">>>> End Add noise");
  }
@@ -644 +654 @@
    PrtTim(">>>> End Scale cube");
  }
+ fluct3d.NTupleCheck(posobs,string("ntfin"),ntnent);
 
  //----------------------------------------------------------------- 
@@ -664 +675 @@
  cout<<endl<<"\n--- ReComputing spectrum from real space"<<endl;
  fluct3d.ReComputeFourier();
+ fluct3d.NTupleCheck(posobs,string("ntpkrec"),ntnent);
  PrtTim(">>>> End ReComputing spectrum");
 

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -15 +15 @@
 #include "fabtwriter.h"
 #include "fioarr.h"
+#include "ntuple.h"
 
 #include "arrctcast.h"
@@ -83 +84 @@
  nthread_ = 0;
  lp_ = 0;
- array_allocated_ = false;
+ array_allocated_ = false; array_type = 0;
  cosmo_ = NULL;
  growth_ = NULL;
@@ -153 +154 @@
  try {
    T_.ReSize(3,SzK_);
-   array_allocated_ = true;
+   array_allocated_ = true; array_type=0;
    if(lp_>1) cout<<"  allocating: "<<T_.Size()*sizeof(complex<GEN3D_TYPE>)/1.e6<<" Mo"<<endl;
  } catch (...) {
@@ -552 +553 @@
 
 //-------------------------------------------------------
+void GeneFluct3D::NTupleCheck(POutPersist &pos,string ntname,unsigned long nent)
+// Remplit le NTuple "ntname" avec "nent" valeurs du cube (reel ou complex) et l'ecrit dans "pos"
+{
+  if(ntname.size()<=0 || nent==0) return;
+  int nvar = 0;
+  if(array_type==1) nvar = 3;
+  else if(array_type==2) nvar = 4;
+  else return;
+  char *vname[4] = {"t","z","re","im"};
+  float xnt[4];
+  NTuple nt(nvar,vname);
+
+  if(array_type==1) {
+    unsigned long nmod = Nx_*Ny_*Nz_/nent; if(nmod==0) nmod=1;
+    unsigned long n=0;
+    for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
+      if(n==nmod) {
+        int_8 ip = IndexR(i,j,l);
+        xnt[0]=sqrt(i*i+j*j); xnt[1]=l; xnt[2]=data_[ip];
+        nt.Fill(xnt);
+        n=0;
+      }
+      n++;
+    }
+  } else {
+    unsigned long nmod = Nx_*Ny_*NCz_/nent; if(nmod==0) nmod=1;
+    unsigned long n=0;
+    for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<NCz_;l++) {
+      if(n==nmod) {
+        xnt[0]=sqrt(i*i+j*j); xnt[1]=l; xnt[2]=T_(l,j,i).real(); xnt[3]=T_(l,j,i).imag();
+        nt.Fill(xnt);
+        n=0;
+      }
+      n++;
+    }
+  }
+
+  pos.PutObject(nt,ntname);
+}
+
+//-------------------------------------------------------
 void GeneFluct3D::Print(void)
 {
@@ -611 +653 @@
  }
 
+ array_type = 2;
+
  if(lp_>0) cout<<"...computing power"<<endl;
  double p = compute_power_carte();
@@ -654 +698 @@
  }
 
+ array_type = 2;
  manage_coefficients();   // gros effet pour les spectres que l'on utilise !
 
@@ -851 +896 @@
  // On fait la FFT
  GEN3D_FFTW_EXECUTE(pb_);
+ array_type = 1;
 }
 
@@ -861 +907 @@
  // On fait la FFT
  GEN3D_FFTW_EXECUTE(pf_);
+ array_type = 2;
+
  // On corrige du pb de la normalisation de FFTW3
  double v = (double)NRtot_;

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -137 +137 @@
   void ReadPPF(string cfname);
   void WriteSlicePPF(string cfname);
+  void NTupleCheck(POutPersist &pos,string ntname,unsigned long nent);
 
   void SetPrtLevel(int lp=0) {lp_ = lp;}
@@ -177 +178 @@
   // le stockage du Cube de donnees et les pointeurs
   bool array_allocated_;  // true if array has been allocated
+  unsigned short array_type; // 0=empty, 1=real, 2=complex
   TArray< complex<GEN3D_TYPE> > T_;
   GEN3D_FFTW_COMPLEX *fdata_;