Changeset 3971 in Sophya for trunk/Cosmo
- Timestamp:
- Apr 1, 2011, 5:47:26 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Cosmo/SimLSS/cmvginit3d.cc
r3924 r3971 37 37 <<" -8 sigmaR,R : normalisation du spectre de puissance, R en Mpc"<<endl 38 38 <<" (default sigmaR=1, R=8/h100 Mpc)"<<endl 39 <<" -L : use no-oscillation (no-BAO) spectrum"<<endl 39 40 <<" -f cambspec.dat,h100tab,ztab : use CAMB file (def: Eisenstein spectrum)"<<endl 40 41 <<" -F : filter spectrum by pixel shape (0=no 1=yes(default)"<<endl … … 91 92 // *** Spectrum and variance definition 92 93 double ns = 1., as = 1.; 94 bool use_noosc = false; 93 95 double R=8./h100, Rg=R/sqrt(5.); 94 96 double sigmaR = 1.; … … 138 140 // --- Decodage des arguments 139 141 char c; 140 while((c = getopt(narg,arg,"haG:F: KB:x:y:z:Z:128:v:n:CO:So:VT:f:")) != -1) {142 while((c = getopt(narg,arg,"haG:F:LKB:x:y:z:Z:128:v:n:CO:So:VT:f:")) != -1) { 141 143 int nth = 0; 142 144 switch (c) { … … 149 151 case 'F' : 150 152 filter_by_pixel = atoi(optarg); 153 break; 154 case 'L' : 155 use_noosc = true; 151 156 break; 152 157 case 'K' : … … 223 228 cout<<"kmin="<<kmin<<" ("<<log10(kmin)<<"), kmax="<<kmax<<" ("<<log10(kmax)<<") Mpc^-1"<<", npt="<<npt<<endl; 224 229 cout<<"Filter by pixel = "<<filter_by_pixel<<endl; 230 if(use_noosc) cout<<"Use spectrum without BAO"<<endl; 225 231 if(kaiser_modify) cout<<"Modify spectrum with Kaiser redshift distorted formula"<<endl; 226 232 if(comptransveloc) cout<<"Tansverse velocity generation requested"<<endl; … … 270 276 271 277 PkSpecCalc pkznos(pkini,tfnos,growth,zref); 278 PkSpecCalc pkzwos(pkini,tf,growth,zref); 272 279 PkSpectrum* pkz = NULL; 273 280 if(cambfread.size()>0) { // read pk in CAMB file 281 cout<<"...using CAMB spectrum"<<endl; 274 282 double h100tab, ztab; 275 283 string fn = decodecambarg(cambfread,h100tab,ztab); … … 285 293 kmax = pkzt->KMax(); 286 294 cout<<"Change k limits to kmin="<<kmin<<" kmax="<<kmax<<endl; 287 } else { // use Eisenstein pk 295 } else if(use_noosc) { // use Eisenstein pk without BAO 296 cout<<"...using spectrum WITHOUT BAO"<<endl; 297 PkSpecCalc* pkze = new PkSpecCalc(pkini,tfnos,growth,zref); 298 pkz = pkze; 299 } else { // use Eisenstein pk with BAO 300 cout<<"...using spectrum with BAO"<<endl; 288 301 PkSpecCalc* pkze = new PkSpecCalc(pkini,tf,growth,zref); 289 302 pkz = pkze; … … 294 307 295 308 pkz->SetZ(0.); 309 pkzwos.SetZ(0.); 296 310 pkznos.SetZ(0.); 297 311 298 312 double normpkz[2] = {0.,0.}; 299 313 for(int i=0;i<2;i++) { 300 PkSpectrum* pkvar = NULL; string nam = " ";301 if(i==0) pkvar = pkz; else {pkvar = &pkznos; nam = "NoOsc";}314 PkSpectrum* pkvar = NULL; string nam = "With-BAO"; 315 if(i==0) pkvar = &pkzwos; else {pkvar = &pkznos; nam = "No-BAO";} 302 316 cout<<endl<<"\n--- Compute variance for Pk "<<nam<<" with top-hat R="<<R<<" at z="<<pkvar->GetZ()<<endl; 303 317 VarianceSpectrum varpk_th(*pkvar,R,VarianceSpectrum::TOPHAT); … … 317 331 cout<<endl; 318 332 333 pkzwos.SetScale(normpkz[0]); 334 pkznos.SetScale(normpkz[1]); 319 335 if(cambfread.size()>0) { 320 336 pkz->SetScale(1.); 321 pkznos.SetScale(normpkz[1]);322 337 cout<<"Warning: CAMB spectrum, no normalisation applied, scale="<<pkz->GetScale()<<endl; 323 338 } else { 324 pkz->SetScale(normpkz[0]); 325 pkznos.SetScale(normpkz[0]); 326 cout<<"Spectrum normalisation (osc+noosc), scale = "<<pkz->GetScale()<<endl; 339 pkz->SetScale( ((use_noosc) ? normpkz[1]: normpkz[0]) ); 340 cout<<"Spectrum normalisation with scale = "<<pkz->GetScale()<<endl; 327 341 } 328 342 … … 334 348 FuncToHisto(pkznos,hpkz0nos,true); 335 349 posobs.PutObject(hpkz0nos,"hpkz0nos"); 350 Histo hpkz0wos(hpkz0); 351 FuncToHisto(pkzwos,hpkz0wos,true); 352 posobs.PutObject(hpkz0wos,"hpkz0wos"); 336 353 } 337 354 338 355 pkz->SetZ(zref); 339 356 pkznos.SetZ(zref); 357 pkzwos.SetZ(zref); 340 358 341 359 { … … 346 364 FuncToHisto(pkznos,hpkznos,true); 347 365 posobs.PutObject(hpkznos,"hpkznos"); 366 Histo hpkzwos(hpkz); 367 FuncToHisto(pkzwos,hpkzwos,true); 368 posobs.PutObject(hpkzwos,"hpkzwos"); 348 369 } 349 370
Note:
See TracChangeset
for help on using the changeset viewer.