Changeset 3924 in Sophya for trunk/Cosmo

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

Location:

trunk/Cosmo/SimLSS

Files:

• cmvginit3d.cc (modified) (10 diffs)
• genefluct3d.cc (modified) (6 diffs)
• genefluct3d.h (modified) (1 diff)

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
 
  //-----------------------------------------------------------------

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -904 +904 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ToRedshiftSpace(void)
+void GeneFluct3D::ToRedshiftSpace(int type_beta,double fixed_beta)
 // Apply redshift distortion corrections
-// ---- ATTENTION: Le spectre Pk doit etre (dRho/Rho)(k)
+// P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
+// type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
+//            >0 compute factor with fix beta=Omega^fixed_beta
+//            <0 compute factor with fix beta=fixed_beta
+// ---- ATTENTION:
+// 1- Le spectre Pk doit etre (dRho/Rho)(k)
+// 2- pas d'evolution en redhsuft possible,
+//    le spectre doit etre calcule au redshift final de la boite
 {
  double zpk = compute_pk_redsh_ref_;
- double beta = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
- if(lp_>0) cout<<"--- ToRedshiftSpace --- at z="<<zpk<<", Beta="<<beta
-               <<" (comp. for dz="<<good_dzinc_<<")"<<endl;
+ if(lp_>0) cout<<"--- ToRedshiftSpace --- at z="<<zpk<<", type_beta="<<type_beta<<endl;
+
+ double beta = 0.;
+ if(type_beta>0) {
+   double omegam = cosmo_->Omatter(zpk); beta = pow(omegam,fixed_beta);
+   if(lp_>0) cout<<"...beta = OmegaM^expo = "<<omegam<<"^"<<fixed_beta<<" = "<<beta<<endl;
+ } else if(type_beta<0) {
+   beta = fixed_beta;
+   if(lp_>0) cout<<"...beta = "<<beta<<endl;
+ } else {
+   beta = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+   if(lp_>0) cout<<"...beta = -(1+z)*D'/D = "<<beta<<endl;
+ }
+
  check_array_alloc();
 
@@ -929 +947 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ToVelLoS(void)
+void GeneFluct3D::ToVelLoS(int type_beta,double fixed_beta)
 /*
+// P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
+// type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
+//            >0 compute factor with fix beta=Omega^fixed_beta
+//            <0 compute factor with fix beta=fixed_beta
 ---- Vitesse particuliere
 Un atome au redshift "z" emet a la longueur d'onde "Le" si il est au repos
@@ -948 +970 @@
 {
  double zpk = compute_pk_redsh_ref_;
- double dpsd = -cosmo_->H(zpk) * (1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
- if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", -H*(1+z)*D'/D="<<dpsd
-               <<" (km/s)/Mpc (comp. for dz="<<good_dzinc_<<")"<<endl;
+ if(lp_>0) cout<<"--- ToVelLoS --- at z="<<zpk<<", type_beta="<<type_beta<<endl;
+
+ double beta = 0;
+ if(type_beta>0) {
+   double omegam = cosmo_->Omatter(zpk);
+   beta = pow(omegam,fixed_beta);
+   if(lp_>0) cout<<"...beta = OmegaM^expo = "<<omegam<<"^"<<fixed_beta<<"="<<beta<<endl;
+ } else if(type_beta<0) {
+   beta = fixed_beta;
+   if(lp_>0) cout<<"...beta="<<beta<<endl;
+ } else {
+   beta = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+   if(lp_>0) cout<<"...beta = -(1+z)*D'/D = "<<beta<<endl;
+ }
+ double dpsd = cosmo_->H(zpk) * beta;
+ if(lp_>0) cout<<"......H*beta = "<<dpsd<<" (km/s)/Mpc"<<endl;
+
  check_array_alloc();
 
@@ -1019 +1055 @@
 
 //-------------------------------------------------------------------
-void GeneFluct3D::ApplyVelLosGrowthFactor(int type_evol)
+void GeneFluct3D::ApplyVelLosGrowthFactor(int type_evol,int type_beta,double fixed_beta)
 // Apply Growth(z) to transverse velocity real space cube
-{
- if(lp_>0) cout<<"--- ApplyVelLosGrowthFactor: evol="<<type_evol<<endl;
+// P(k)_z = P(k) * ( 1 + beta * mu^2 )^2
+// type_beta: =0 compute factor automatically with beta=-(1+z)/D*dD/z
+//            >0 compute factor with fix beta=Omega^fixed_beta
+//            <0 compute factor with fix beta=fixed_beta
+{
+ if(lp_>0) cout<<"--- ApplyVelLosGrowthFactor: evol="<<type_evol<<", type_beta="<<type_beta<<endl;
+
  check_array_alloc();
 
@@ -1036 +1077 @@
  double zpk = compute_pk_redsh_ref_;
  double grw_orig = (*growth_)(zpk);
- double dpsd_orig = - cosmo_->H(zpk) * (1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / grw_orig;
- if(lp_>0) cout<<"    original growth="<<grw_orig<<" dpsd="<<dpsd_orig<<" computed at z="<<zpk<<endl;
+ double beta_orig = 0.;
+ if(type_beta>0) {double omegam = cosmo_->Omatter(zpk); beta_orig = pow(omegam,fixed_beta);}
+ else if(type_beta<0) beta_orig = fixed_beta;
+ else beta_orig = -(1.+zpk) * growth_->DsDz(zpk,good_dzinc_) / (*growth_)(zpk);
+ double dpsd_orig = cosmo_->H(zpk) * beta_orig;
+ if(lp_>0) cout<<"...original at z="<<zpk<<": growth="<<grw_orig<<" beta="<<beta_orig
+               <<" H*beta="<<dpsd_orig<<" (km/s)/Mpc"<<endl;
+
+ if(type_beta<0) cout<<"***WARNING: type_beta<0 , no evolution used for beta !!!"<<endl;
 
  InterpFunc interpinv(loscom2zred_min_,loscom2zred_max_,loscom2zred_);
@@ -1052 +1100 @@
        double z = interpinv(dz);   // interpolation par morceau
        double grw = (*growth_)(z);
-       double dpsd = interdpd(z);
+       double dpsd = 0;
+       if(type_beta>0) {
+         dpsd = cosmo_->H(z) * pow(cosmo_->Omatter(z),fixed_beta);
+       } else if(type_beta<0) {
+         dpsd = cosmo_->H(z) * fixed_beta;
+       } else {
+         dpsd = interdpd(z); // interpole -H(z)*(1+z)*D'(z)/D(z)
+       }
        int_8 ip = IndexR(i,j,l);
        // on remet le beta au bon z

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -118 +118 @@
   void ComputeFourier(PkSpectrum& pk_at_z);
   void FilterByPixel(void);
-  void ToRedshiftSpace(void);
-  void ToVelLoS(void);
+  void ToRedshiftSpace(int type_beta=0,double fixed_beta=0.);
+  void ToVelLoS(int type_beta=0,double fixed_beta=0.);
 
   void ComputeReal(void);
   void ApplyGrowthFactor(int type_evol=1);
   void ApplyRedshiftSpaceGrowthFactor(void);
-  void ApplyVelLosGrowthFactor(int type_evol=1);
+  void ApplyVelLosGrowthFactor(int type_evol=1,int type_beta=0,double fixed_beta=0.);
 
   void ReComputeFourier(void);