Changeset 3924 in Sophya for trunk/Cosmo/SimLSS/cmvginit3d.cc

Timestamp:

Dec 10, 2010, 4:38:23 PM (15 years ago)

Author:

cmv

Message:

intro choix dabs calcul de beta pour les redshift-distorsions, cmv 10/12/2010

File:

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

trunk/Cosmo/SimLSS/cmvginit3d.cc

@@ -27 +27 @@
 {
  cout<<"cmvginit3d [...options...]"<<endl
-     <<" -a : auto init random seed (needed for multiple simul)"<<endl
+     <<"-- redshift"<<endl
+     <<" -Z zref : redshift for the center of the simulation cube"<<endl
+     <<"-- geometry"<<endl
+     <<" -x nx,dx : size along x axis (npix,Mpc)"<<endl
+     <<" -y ny,dy : size along y axis (npix,Mpc)"<<endl
+     <<"            if ny or dy <=0 take same value as for x"<<endl
+     <<" -z nz,dz : size along z axis (redshift axis, npix,Mpc)"<<endl
+     <<"-- input spectra"<<endl
+     <<" -8 sigmaR,R : normalisation du spectre de puissance, R en Mpc"<<endl
+     <<"               (default sigmaR=1, R=8/h100 Mpc)"<<endl
+     <<" -f cambspec.dat,h100tab,ztab : use CAMB file (def: Eisenstein spectrum)"<<endl
+     <<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl
+     <<"-- redshift evolution"<<endl
      <<" -G typevol: compute Pk(z=0) and apply growth factor in real space"<<endl
      <<"       typevol=1 evolved with distance / observateur (def)"<<endl
      <<"       typevol=2 evolved with distance to middle of Z planes"<<endl
      <<"       else : no evol, spectrum Pk(z=z_median) for all cube (def)"<<endl
-     <<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl
+     <<"-- redshift distorsions"<<endl
+     <<" -B type,fixed_beta : redshift distorted P(k)_z = P(k)*(1+beta*mu^2)^2"<<endl
+     <<"      type=0 : compute beta = -(1+z)*D'/D (default)"<<endl
+     <<"      type>0 : compute beta = OmegaM(z)^fixed_beta"<<endl
+     <<"      type<0 : use beta = fixed_beta"<<endl
      <<" -K : modify spectrum with Kaiser redshift distortion (default=no)"<<endl
-     <<" -x nx,dx : size along x axis (npix,Mpc)"<<endl
-     <<" -y ny,dy : size along y axis (npix,Mpc)"<<endl
-     <<"            if ny or dy <=0 take same value as for x"<<endl
-     <<" -z nz,dz : size along z axis (redshift axis, npix,Mpc)"<<endl
-     <<" -Z zref : redshift for the center of the simulation cube"<<endl
-     <<" -f cambspec.dat,h100tab,ztab : use CAMB file (def: Eisenstein spectrum)"<<endl
+     <<" -v temp_file.ppf: generate los velocity cube"<<endl
+     <<"                   temporary save cube in temp_file.ppf"<<endl
+     <<"-- program init and efficiency"<<endl
+     <<" -a : auto init random seed (needed for multiple simul)"<<endl
+     <<" -T nth : nombre de threads (si compil multi-thread, default: 0)"<<endl
+     <<"-- check by re-computing Pk spectra, ntuple"<<endl
      <<" -1 : compute 1D spectrum (default: no)"<<endl
      <<" -2 : compute 2D spectrum (default: no)"<<endl
      <<" -C : go back to cube in FT space and compute P(k) for check (def: do nothing)"<<endl
-     <<" -8 sigmaR,R : normalisation du spectre de puissance, R en Mpc"<<endl
-     <<"               (default sigmaR=1, R=8/h100 Mpc)"<<endl
-     <<" -v temp_file.ppf: generate los velocity cube"<<endl
-     <<"                   temporary save cube in temp_file.ppf"<<endl
      <<" -n nent : fill control ntuple with nent entries (def: do not fill)"<<endl
+     <<" -V : compute variance from real space (for check, default: no)"<<endl
+     <<"-- output type selection"<<endl
+     <<" -o out_base_name : base string for output file name (def: cmvginit3d)"<<endl
      <<" -O a,b : tell what you want to write"<<endl
      <<"      a : write generated fourier cube (_k0)"<<endl
@@ -53 +68 @@
      <<"      a,b= 0 no write, 1 ppf write, 2 fits write, 3 ppf+fits write"<<endl
      <<" -S : write cube slices in PPF format"<<endl
-     <<" -o out_base_name : base string for output file name (def: cmvginit3d)"<<endl
-     <<" -V : compute variance from real space (for check, default: no)"<<endl
-     <<" -T nth : nombre de threads (si compil multi-thread, default: 0)"<<endl
      <<endl;
 }
@@ -93 +105 @@
  unsigned short nthread=0;
  int filter_by_pixel = 1;
- bool kaiser_modify = false;
-
- // *** What to do
+
+ // *** traitement des distorsions de redshift
  bool comptransveloc = false;
  string temporay_file = "cmvginit3d_tmp.ppf";
+ bool kaiser_modify = false;
+ int type_beta = 0;
+ double fixed_beta = 0.6;
+
+ // *** What to do
  // compute: 1=1D spectra,  2=2D spectra,
  //          4=recompute spectra after real space generation
@@ -122 +138 @@
  // --- Decodage des arguments
  char c;
- while((c = getopt(narg,arg,"haG:F:Kx:y:z:Z:128:v:n:CO:So:VT:f:")) != -1) {
+ while((c = getopt(narg,arg,"haG:F:KB:x:y:z:Z:128:v:n:CO:So:VT:f:")) != -1) {
   int nth = 0;
   switch (c) {
@@ -135 +151 @@
     break;
   case 'K' :
-    kaiser_modify =  true;
+    kaiser_modify = true;
+    break;
+  case 'B' :
+    sscanf(optarg,"%d,%lf",&type_beta,&fixed_beta);
     break;
   case 'x' :
@@ -206 +225 @@
  if(kaiser_modify) cout<<"Modify spectrum with Kaiser redshift distorted formula"<<endl;
  if(comptransveloc) cout<<"Tansverse velocity generation requested"<<endl;
+ if(type_beta>0) cout<<"Beta: Pk factor is (1 + OmegaM^"<<fixed_beta<<" * mu^2)^2"<<endl;
+   else if(type_beta<0) cout<<"Beta: Pk factor is (1 + "<<fixed_beta<<" * mu^2)^2"<<endl;
+     else cout<<"Beta: automatically computed with dD/dEta"<<endl;
  if(cambfread.size()>0) cout<<"use CAMB file: cambfread="<<cambfread<<endl;
  cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
@@ -414 +436 @@
  if(kaiser_modify) {
    cout<<"\n--- Modify spectrum coeff with Kaiser redshift distorted formula"<<endl;
-   fluct3d.ToRedshiftSpace();
+   fluct3d.ToRedshiftSpace(type_beta,fixed_beta);
    spec_is_modified = true;
    PrtTim(">>>> End Computing Modify spectrum coeff with Kaiser");
@@ -571 +593 @@
    //-----------------------------------------------------------------
    cout<<"\n--- Modifying cube for Transverse velocity"<<endl;
-   fluct3d.ToVelLoS();
+   fluct3d.ToVelLoS(type_beta,fixed_beta);
    fluct3d.NTupleCheck(posobs,string("ntpvgenf"),ntnent);
    PrtTim(">>>> End Modifying cube for Transverse velocity");
@@ -614 +636 @@
       cout<<"\n--- Apply Growth factor for transverse velocity"<<endl;
       cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
-      fluct3d.ApplyVelLosGrowthFactor(use_growth_factor);
+      fluct3d.ApplyVelLosGrowthFactor(use_growth_factor,type_beta,fixed_beta);
       fluct3d.MinMax(vdum1,vdum2);
       fluct3d.MeanSigma2(vdum1,vdum2);
@@ -637 +659 @@
     }
 
- }
+ }  //// Fin de Transverse velocity space computation
 
  //-----------------------------------------------------------------