Changeset 3349 in Sophya for trunk/Cosmo/SimLSS/cmvobserv3d.cc


Ignore:
Timestamp:
Oct 11, 2007, 4:39:58 PM (18 years ago)
Author:
cmv
Message:
  • gros changements dans la structure de la classe GeneFluct3D (constructeur et logique d'aloocation memoire, init_fftw etc...)
  • suppression des valeurs de masse<0 mises a -999. directement mises a zero
  • suppression de TurnMass2HIMass qui fait maintenant la meme chose que TurnMass2MeanNumber
  • legere restructuration de cmvobserv3d.cc pour compat.

cmv 11/10/2007

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/cmvobserv3d.cc

    r3345 r3349  
    4545     <<" -R schmassdist.ppf : read mass distribution for trials from file"<<endl
    4646     <<"               instead of computing it (ONLY if \"-Q\" option is activated)"<<endl
    47      <<" -A <log10(S_agn in Jy at 1.4 GHz)>,sigma,powlaw :"<<endl
    48      <<"    AGN mean and sigma gaussian equiv. distrib. for solid angle of centeral pixel"<<endl
    49      <<"    powlaw: apply S_agn evolution as (Nu/1.4)^powlaw"<<endl
    5047     <<" -W : write cube in FITS format (complex cube is coded as real cube)"<<endl
    5148     <<" -P : write cube in PPF format"<<endl
    … …  
    5451     <<" -T nth : nombre de threads (si compil multi-thread, default: 0)"<<endl
    5552     <<endl;
     53     ////<<" -A <log10(S_agn in Jy at 1.4 GHz)>,sigma,powlaw :"<<endl
     54     ////<<"    AGN mean and sigma gaussian equiv. distrib. for solid angle of centeral pixel"<<endl
     55     ////<<"    powlaw: apply S_agn evolution as (Nu/1.4)^powlaw"<<endl
    5656}
    5757
    … …  
    101101 int noise_evol = 0;
    102102
    103  // *** AGN
    104  bool do_agn = false;
    105  double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
    106  double powlaw_agn = 0.;
     103 //// *** AGN
     104 ////bool do_agn = false;
     105 ////double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
     106 ////double powlaw_agn = 0.;
    107107
    108108 // *** type de generation
    … …  
    177177    schmassdistfile = optarg;
    178178    break;
    179   case 'A' :
    180     do_agn = true;
    181     sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
    182     break;
     179    ////  case 'A' :
     180    ////do_agn = true;
     181    ////sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
     182    ////break;
    183183  case 'V' :
    184184    compvarreal = true;
    … …  
    224224 cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<noise_evol<<endl;
    225225 cout<<"scalecube="<<scalecube<<", offsetcube="<<offsetcube<<endl;
    226  if(do_agn)
    227    cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
    228        <<" , powlaw="<<powlaw_agn<<endl;
     226 ////if(do_agn)
     227 ////  cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
     228 ////      <<" , powlaw="<<powlaw_agn<<endl;
    229229
    230230 string tagobs = "cmvobserv3d.ppf";
    … …  
    264264 cout<<endl<<"\n--- Compute variance for top-hat R="<<R
    265265     <<" at z="<<pkz.GetZ()<<endl;
    266  VarianceSpectrum varpk_th(pkz,0);
    267  double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
     266 VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
     267 double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
    268268 double pkmax_th = varpk_th(kfind_th);
    269269 cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
    270270 double k1=kmin, k2=kmax;
    271  int rc = varpk_th.FindLimits(R,pkmax_th/1.e4,k1,k2,eps);
     271 int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
    272272 cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
    273273     <<k2<<" ("<<log10(k2)<<")"<<endl;
    … …  
    275275 double ldlk = (log10(k2)-log10(k1))/npt;
    276276 varpk_th.SetInteg(0.01,ldlk,-1.,4);
    277  double sr2 = varpk_th.Variance(R,k1,k2);
     277 double sr2 = varpk_th.Variance(k1,k2);
    278278 cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;
    279279
    … …  
    292292 //-----------------------------------------------------------------
    293293 cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
    294  VarianceSpectrum varpk_int(pkz,2);
    295 
    296  double kfind_int = varpk_int.FindMaximum(R,kmin,kmax,eps);
     294 VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
     295
     296 double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
    297297 double pkmax_int = varpk_int(kfind_int);
    298298 cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
    299299 double k1int=kmin, k2int=kmax;
    300  int rcint = varpk_int.FindLimits(R,pkmax_int/1.e4,k1int,k2int,eps);
     300 int rcint = varpk_int.FindLimits(pkmax_int/1.e4,k1int,k2int,eps);
    301301 cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
    302302     <<k2int<<" ("<<log10(k2int)<<")"<<endl;
    … …  
    304304 double ldlkint = (log10(k2int)-log10(k1int))/npt;
    305305 varpk_int.SetInteg(0.01,ldlkint,-1.,4);
    306  double sr2int = varpk_int.Variance(R,k1int,k2int);
     306 double sr2int = varpk_int.Variance(k1int,k2int);
    307307 cout<<"varpk_int="<<sr2int<<"  ->  sigma="<<sqrt(sr2int)<<endl;
    308308 
    … …  
    355355 cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
    356356
    357  TArray< complex<r_8> > pkgen;
    358  GeneFluct3D fluct3d(pkgen);
    359  fluct3d.SetPrtLevel(2);
    360  fluct3d.SetNThread(nthread);
    361  fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
     357 GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
    362358 fluct3d.SetObservator(zref,-nz/2.);
    363359 fluct3d.SetCosmology(univ);
    364360 fluct3d.SetGrowthFactor(growth);
    365361 fluct3d.LosComRedshift(0.001,-1);
    366  //TArray<r_8>& rgen = fluct3d.GetRealArray();
     362 TArray< complex<r_8> >& pkgen = fluct3d.GetComplexArray();
     363 TArray<r_8>& rgen = fluct3d.GetRealArray();
    367364 cout<<endl; fluct3d.Print();
    368365 cout<<"\nMean number of galaxies per pixel = "<<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
    … …  
    391388 ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
    392389 varpk_int.SetInteg(0.01,ldlkint,-1.,4);
    393  double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax_mod);
     390 double sr2int_kmax = varpk_int.Variance(k1int,knyqmax_mod);
    394391 cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
    395392
    … …  
    538535
    539536 if(no_poisson) {
    540 
    541537   cout<<"\n--- Converting !!!DIRECTLY!!! mass into HI mass: mass per pixel ="
    542538       <<mass_by_mpc3*fluct3d.GetDVol()<<endl;
    543    rm = fluct3d.TurnMass2HIMass(mass_by_mpc3);
    544      nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
    545      nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
    546    PrtTim(">>>> End Converting !!!DIRECTLY!!! mass into HI mass");
    547 
     539   rm = fluct3d.TurnMass2MeanNumber(mass_by_mpc3);
    548540 } else {
    549 
    550541   cout<<"\n--- Converting mass into galaxy number: gal per pixel ="
    551542       <<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
    552543   rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
    553      nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
    554      nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
    555    PrtTim(">>>> End Converting mass into galaxy number");
    556 
    557    cout<<"\n--- Set negative and null pixels to BAD"<<endl;
    558    nm = fluct3d.SetToVal(0.,1e+200,-999.);
    559    PrtTim(">>>> End Set negative pixels to BAD etc...");
     544 }
     545 nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
     546 nm = fluct3d.MeanSigma2(rm,rs2);
     547 PrtTim(">>>> End Converting mass into galaxy number or mass");
     548
     549 if( !no_poisson ) {
    560550
    561551   cout<<"\n--- Apply poisson on galaxy number"<<endl;
    562552   fluct3d.ApplyPoisson();
    563      nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
    564      nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
     553     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
     554     nm = fluct3d.MeanSigma2(rm,rs2);
    565555   double xgalmin,xgalmax; fluct3d.MinMax(xgalmin,xgalmax,0.1,1.e50);
    566556   PrtTim(">>>> End Apply poisson on galaxy number");
    … …  
    584574   }
    585575   cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
    586      nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
     576     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
    587577     cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
    588      nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
     578     nm = fluct3d.MeanSigma2(rm,rs2);
    589579   PrtTim(">>>> End Convert Galaxy number to HI mass");
    590 
    591    cout<<"\n--- Set BAD pixels to Zero"<<endl;
    592    nm = fluct3d.SetToVal(-998.,1e+200,0.);
    593    PrtTim(">>>> End Set BAD pixels to Zero etc...");
    594580
    595581 }
    … …  
    609595
    610596 //-----------------------------------------------------------------
    611  if(do_agn) {
    612    cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
    613        <<" , powlaw="<<powlaw_agn<<endl;
    614    fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
    615      nm = fluct3d.MeanSigma2(rm,rs2);
    616    PrtTim(">>>> End Add AGN");
    617  }
     597 ////if(do_agn) {
     598 ////  cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
     599 ////      <<" , powlaw="<<powlaw_agn<<endl;
     600 ////  fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
     601 ////    nm = fluct3d.MeanSigma2(rm,rs2);
     602 ////   PrtTim(">>>> End Add AGN");
     603 ////}
    618604
    619605 //-----------------------------------------------------------------
    … …  
    627613 }
    628614
     615
     616 //-----------------------------------------------------------------
    629617 if(scalecube!=0.) {  // Si scalecube==0 pas de normalisation
    630618   cout<<"\n--- Scale cube rs2ref="<<rs2ref<<endl;
Note: See TracChangeset for help on using the changeset viewer.