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

Timestamp:

Jul 24, 2010, 6:15:07 PM (15 years ago)

Author:

cmv

Message:

suite de la mise au point pour lecture fichiers CAMB, cmv 24/07/2010

File:

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

trunk/Cosmo/SimLSS/cmvginit3d.cc

@@ -20 +20 @@
 #include "geneutils.h"
 #include "genefluct3d.h"
+
+string decodecambarg(string arg,double& h, double& z);
 
 void usage(void);
@@ -36 +38 @@
      <<" -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
      <<" -1 : compute 1D spectrum (default: no)"<<endl
      <<" -2 : compute 2D spectrum (default: no)"<<endl
@@ -69 +72 @@
  unsigned short flat = 0;
  double ob0 = 0.0444356;
- double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
+ double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73, w0=-1.;
  double zref = 0.5;
  double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
@@ -78 +81 @@
  double sigmaR = 1.;
 
- double kmin=1e-5,kmax=1000.;
- int npt = 10000;
- double lkmin=log10(kmin), lkmax=log10(kmax);
+ double kmin=1e-4,kmax=100.;
+ int npt = 1000;
  double eps=1.e-3;
+
+ // *** Spectrum read from CAMB file
+ string cambfread = "";
 
  // *** type de generation
@@ -115 +120 @@
  // --- Decodage des arguments
  char c;
- while((c = getopt(narg,arg,"haG:F:x:y:z:Z:128:v:n:CO:So:VT:")) != -1) {
+ while((c = getopt(narg,arg,"haG:F:x:y:z:Z:128:v:n:CO:So:VT:f:")) != -1) {
   int nth = 0;
   switch (c) {
@@ -175 +180 @@
     nthread = (nth<1)? 0: nth;
     break;
+  case 'f' :
+    cambfread = optarg;
+    break;
   case 'h' :
   default :
@@ -189 +197 @@
  cout<<"zref="<<zref<<endl;
  cout<<"nx="<<nx<<" dx="<<dx<<" ny="<<ny<<" dy="<<dy<<" nz="<<nz<<" dz="<<dz<<endl;
- cout<<"kmin="<<kmin<<" ("<<lkmin<<"), kmax="<<kmax<<" ("<<lkmax<<") Mpc^-1"<<", npt="<<npt<<endl;
+ cout<<"kmin="<<kmin<<" ("<<log10(kmin)<<"), kmax="<<kmax<<" ("<<log10(kmax)<<") Mpc^-1"<<", npt="<<npt<<endl;
  cout<<"Filter by pixel = "<<filter_by_pixel<<endl;
  if(comptransveloc) cout<<"Tansverse velocity generation requested"<<endl;
+ if(cambfread.size()>0) cout<<"use CAMB file: cambfread="<<cambfread<<endl;
  cout<<"R="<<R<<" Rg="<<Rg<<" Mpc, sigmaR="<<sigmaR<<endl;
  cout<<"Use_growth_factor = "<<use_growth_factor<<endl;
@@ -216 +225 @@
 
  //-----------------------------------------------------------------
- cout<<endl<<"\n--- Create Spectrum"<<endl;
-
- InitialPowerLaw pkini(ns,as);
-
- TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
- TransfertEisenstein tfnos(h100,om0-ob0,ob0,T_CMB_Par,false);
-   tfnos.SetNoOscEnv(2);
+ cout<<endl<<"\n--- Create GrowthFactor"<<endl;
 
  GrowthEisenstein growth(om0,ol0);
@@ -229 +232 @@
  cout<<"...Growth factor at z="<<zref<<" = "<<growth_at_z<<endl;
 
- PkSpecCalc pkz(pkini,tf,growth,zref);
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create Spectrum"<<endl;
+
+ InitialPowerLaw pkini(ns,as);
+
+ TransfertEisenstein tfnos(h100,om0-ob0,ob0,T_CMB_Par,false);
+   tfnos.SetNoOscEnv(2);
+ TransfertEisenstein tf(h100,om0-ob0,ob0,T_CMB_Par,false);
+
  PkSpecCalc pkznos(pkini,tfnos,growth,zref);
-
- //-----------------------------------------------------------------
- pkz.SetZ(0.);
+ PkSpectrum* pkz = NULL;
+ if(cambfread.size()>0) {  // read pk in CAMB file
+   double h100tab, ztab;
+   string fn = decodecambarg(cambfread,h100tab,ztab);
+   if(h100tab<=0.) h100tab = h100;
+   PkTabulate* pkzt = new PkTabulate;
+   pkzt->ReadCAMB(fn,h100tab,ztab);
+   pkzt->SetGrowthFactor(&growth);
+   pkzt->SetZ(zref);
+   pkz = pkzt;
+   // change k limits
+   kmin = pkzt->KMin();
+   kmax = pkzt->KMax();
+   cout<<"Change k limits to kmin="<<kmin<<" kmax="<<kmax<<endl;
+ } else {   // use Eisenstein pk
+   PkSpecCalc* pkze = new PkSpecCalc(pkini,tf,growth,zref);
+   pkz = pkze;
+ }
+
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Normalize Spectrum"<<endl;
+
+ pkz->SetZ(0.);
  pkznos.SetZ(0.);
- cout<<endl<<"\n--- Compute variance for top-hat R="<<R<<" at z="<<pkz.GetZ()<<endl;
- VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
- double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
- double pkmax_th = varpk_th(kfind_th);
- cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
- double k1=kmin, k2=kmax;
- int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
- cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
-     <<k2<<" ("<<log10(k2)<<")"<<endl;
-
- double ldlk = (log10(k2)-log10(k1))/npt;
- varpk_th.SetInteg(0.01,ldlk,-1.,4);
- double sr2 = varpk_th.Variance(k1,k2);
- cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;
-
- double normpkz = sigmaR*sigmaR/sr2;
- pkz.SetScale(normpkz);
- pkznos.SetScale(normpkz);
- cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
+
+ double normpkz[2] = {0.,0.};
+ for(int i=0;i<2;i++) {
+   PkSpectrum* pkvar = NULL; string nam = "";
+   if(i==0) pkvar = pkz; else {pkvar = &pkznos; nam = "NoOsc";}
+   cout<<endl<<"\n--- Compute variance for Pk "<<nam<<" with top-hat R="<<R<<" at z="<<pkvar->GetZ()<<endl;
+   VarianceSpectrum varpk_th(*pkvar,R,VarianceSpectrum::TOPHAT);
+   double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
+   double pkmax_th = varpk_th(kfind_th);
+   cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
+   double k1=kmin, k2=kmax;
+   int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
+   cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
+       <<k2<<" ("<<log10(k2)<<")"<<endl;
+   double ldlk = (log10(k2)-log10(k1))/npt;
+   varpk_th.SetInteg(0.01,ldlk,-1.,4);
+   double sr2 = varpk_th.Variance(k1,k2);
+   normpkz[i] = sigmaR*sigmaR/sr2;
+   cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<" normpkz="<<normpkz[i]<<endl;
+ }
+ cout<<endl;
+
+ if(cambfread.size()>0) {
+   pkz->SetScale(1.);
+   pkznos.SetScale(normpkz[1]);
+   cout<<"Warning: CAMB spectrum, no normalisation applied, scale="<<pkz->GetScale()<<endl;
+ } else {
+   pkz->SetScale(normpkz[0]);
+   pkznos.SetScale(normpkz[0]);
+   cout<<"Spectrum normalisation (osc+noosc), scale = "<<pkz->GetScale()<<endl;
+ }
+
  {
- Histo hpkz0(lkmin,lkmax,npt); hpkz0.ReCenterBin();
- FuncToHisto(pkz,hpkz0,true);
+ Histo hpkz0(log10(kmin),log10(kmax),npt); hpkz0.ReCenterBin();
+ FuncToHisto(*pkz,hpkz0,true);
  posobs.PutObject(hpkz0,"hpkz0");
  Histo hpkz0nos(hpkz0);
@@ -263 +307 @@
  }
 
- pkz.SetZ(zref);
+ pkz->SetZ(zref);
  pkznos.SetZ(zref);
+
  {
- Histo hpkz(lkmin,lkmax,npt); hpkz.ReCenterBin();
- FuncToHisto(pkz,hpkz,true);
+ Histo hpkz(log10(kmin),log10(kmax),npt); hpkz.ReCenterBin();
+ FuncToHisto(*pkz,hpkz,true);
  posobs.PutObject(hpkz,"hpkz");
  Histo hpkznos(hpkz);
@@ -275 +320 @@
 
  //-----------------------------------------------------------------
- cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
- VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
+ cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz->GetZ()<<endl;
+ VarianceSpectrum varpk_int(*pkz,R,VarianceSpectrum::NOFILTER);
  double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
  double pkmax_int = varpk_int(kfind_int);
@@ -325 +370 @@
 
  //-----------------------------------------------------------------
- cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
+ cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz->GetZ()<<endl;
  // En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
  // sphere: Vs = 4Pi/3 k^3 , cube inscrit (cote k*sqrt(2)): Vc = (k*sqrt(2))^3
@@ -339 +384 @@
  //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in Fourier space"<<endl;
- if(use_growth_factor>0) pkz.SetZ(0.); else pkz.SetZ(zref);
- cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
- fluct3d.ComputeFourier(pkz);
+ if(use_growth_factor>0) pkz->SetZ(0.); else pkz->SetZ(zref);
+ cout<<"Power spectrum set at redshift: "<<pkz->GetZ()<<endl;
+ fluct3d.ComputeFourier(*pkz);
  fluct3d.NTupleCheck(posobs,string("ntpkgen"),ntnent);
  PrtTim(">>>> End Computing a realization in Fourier space");
@@ -411 +456 @@
  //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in real space"<<endl;
+ double vdum1,vdum2;
  fluct3d.ComputeReal();
- double rmin,rmax; fluct3d.MinMax(rmin,rmax);
- cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+ fluct3d.MinMax(vdum1,vdum2);
+ fluct3d.MeanSigma2(vdum1,vdum2);
  fluct3d.NTupleCheck(posobs,string("ntreal"),ntnent);
  PrtTim(">>>> End Computing a realization in real space");
@@ -421 +467 @@
    cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
    fluct3d.ApplyGrowthFactor(use_growth_factor);
-   fluct3d.MinMax(rmin,rmax);
-   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   fluct3d.MinMax(vdum1,vdum2);
+   fluct3d.MeanSigma2(vdum1,vdum2);
    fluct3d.NTupleCheck(posobs,string("ntgrow"),ntnent);
    PrtTim(">>>> End Applying growth factor");
  }
-
- int_8 nm;
- double rmref,rs2ref;
- cout<<"\n--- Computing reference variance in real space"<<endl;
- nm = fluct3d.MeanSigma2(rmref,rs2ref);
- cout<<" rs2ref= "<<rs2ref<<" , rmref="<<rmref<<" ("<<nm<<")"<<endl;
- PrtTim(">>>> End Computing reference variance in real space");
 
  if(whattowrt[1]&1) {
@@ -561 +600 @@
    cout<<"\n--- Computing a realization in real space for Transverse velocity"<<endl;
    fluct3d.ComputeReal();
-   double rmin,rmax; fluct3d.MinMax(rmin,rmax);
-   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   fluct3d.MinMax(vdum1,vdum2);
+   fluct3d.MeanSigma2(vdum1,vdum2);
    fluct3d.NTupleCheck(posobs,string("nvreal"),ntnent);
    PrtTim(">>>> End Computing a realization in real space");
 
     if(use_growth_factor>0) {
-      cout<<"\n--- Apply Growth factor"<<endl;
+      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.MinMax(rmin,rmax);
-      cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+      fluct3d.MinMax(vdum1,vdum2);
+      fluct3d.MeanSigma2(vdum1,vdum2);
       fluct3d.NTupleCheck(posobs,string("nvgrow"),ntnent);
       PrtTim(">>>> End Applying growth factor");
     }
-
-    int_8 nmv;
-    double vmref,vs2ref;
-    cout<<"\n--- Computing reference variance in real space"<<endl;
-    nmv = fluct3d.MeanSigma2(vmref,vs2ref);
-    cout<<" vs2ref= "<<vs2ref<<" , vmref="<<vmref<<" ("<<nmv<<")"<<endl;
-    PrtTim(">>>> End Computing reference variance in real space");
 
     if(whattowrt[1]&1) {
@@ -610 +642 @@
  dvl("s8") = sigmaR;
  dvl("Dtype") = (int_4)use_growth_factor; dvl("Pxfilt") = (int_4)filter_by_pixel; 
+ if(cambfread.size()>0) dvl("fCAMB") = cambfread;
  posobs.PutObject(dvl,"siminfo");
  }
 
  delete RandGen;
+ if(pkz != NULL) delete pkz;
  PrtTim(">>>> End Of Job");
 
@@ -634 +668 @@
 }
 
+#include "strutilxx.h"
+string decodecambarg(string arg,double& h, double& z)
+// - decode argument for CAMB generation:
+// input: arg = "cambfilename,h100tab,ztab"
+// output: h = h100 use to generate CAMB file
+//         z = redshift use to generate CAMB file
+// return: CAMB file name
+{
+  h = -1.;
+  z = 0.;
+  vector<string> vs;
+  FillVStringFrString(arg,vs,',');
+  if(vs.size()>1) h = atof(vs[1].c_str());
+  if(vs.size()>2) z = atof(vs[2].c_str());
+  cout<<"decodecambarg: "<<arg<<" fn=\""<<vs[0]<<"\" h="<<h<<" z="<<z<<endl;
+  return vs[0];
+}