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

Timestamp:

Sep 5, 2007, 6:16:45 PM (18 years ago)

Author:

cmv

Message:

Intro SchechterMassDist + lecture sur fichier cmv 05/09/2007

File:

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

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -34 +34 @@
      <<"                  if evol>0 noise evolved with distance (def no)"<<endl
      <<" -2 : compute also 2D spectrum (default: no)"<<endl
+     <<" -N scalecube,offsetcube: normalize cube before doing final spectrum (default: automatic)"<<endl
      <<" -M schmin,schmax,nsch : min,max mass and nb points for schechter HI"<<endl
+     <<"                         If nsch<0 alors no,bre de points par decade"<<endl
+     <<" -Q naleagal : use quick method for turning ngal to mass"<<endl
+     <<" -R schmassdist.ppf : read mass distribution for trials from file"<<endl
+     <<"               instead of computing it (ONLY if \"-Q\" option is activated)"<<endl
      <<" -A <log10(S_agn in Jy at 1.4 GHz)>,sigma,powlaw :"<<endl
      <<"    AGN mean and sigma gaussian equiv. distrib. for solid angle of centeral pixel"<<endl
@@ -79 +84 @@
  double alpha = -1.37;
 
- double schmin=1e8, schmax=1e12;
- int schnpt = 1000;
+ double schmin=1e8, schmax=1e13;
+ int schnpt = -100;
+ bool use_schmassdist = false;
+ long naleagal = 100000;
+ bool recompute_schmassdist = true;
+ string schmassdistfile = "";
 
  // *** Niveau de bruit
@@ -103 +112 @@
  bool wslice = false;
  bool compvarreal = false;
+ double scalecube = -1., offsetcube = 0.;
 
  // --- Decodage arguments
  if(narg>0) {
+   cout<<"\n--- Arguments: "<<endl;
    for(int i=0;i<narg;i++) cout<<arg[i]<<" ";
    cout<<endl;
@@ -111 +122 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWSV2GF:x:y:z:s:Z:M:A:T:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWSV2GF:x:y:z:s:Z:M:A:T:N:Q:R:")) != -1) {
   int nth = 0;
   switch (c) {
@@ -144 +155 @@
     comp2dspec = true;
     break;
+  case 'N' :
+    sscanf(optarg,"%lf,%lf",&scalecube,&offsetcube);
+    break;
   case 'M' :
     sscanf(optarg,"%lf,%lf,%d",&schmin,&schmax,&schnpt);
+    break;
+  case 'Q' :
+    use_schmassdist = true;
+    sscanf(optarg,"%ld",&naleagal);
+    break;
+  case 'R' :
+    schmassdistfile = optarg;
     break;
   case 'A' :
@@ -173 +194 @@
  }
 
- double lschmin=log10(schmin), lschmax=log10(schmax), dlsch=(lschmax-lschmin)/schnpt;
-
- string tagobs = "cmvobserv3d.ppf";
- POutPersist posobs(tagobs);
+ double lschmin=log10(schmin), lschmax=log10(schmax);
+ if(schnpt<=0) {  // alors c'est un nombre de points par decade
+   schnpt = long( (-schnpt)*(lschmax-lschmin+1.) + 0.5 );
+   if(schnpt<=2) schnpt = 1000;
+ }
+ if(naleagal<=2) naleagal = 100000;
 
  cout<<"zref="<<zref<<endl;
@@ -189 +212 @@
      <<", schnpt="<<schnpt<<endl;
  cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<isnoise_evol<<endl;
+ cout<<"scalecube="<<scalecube<<", offsetcube="<<offsetcube<<endl;
  if(do_agn)
    cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
        <<" , powlaw="<<powlaw_agn<<endl;
+
+ string tagobs = "cmvobserv3d.ppf";
+ POutPersist posobs(tagobs);
 
  //-----------------------------------------------------------------
@@ -221 +248 @@
 
  PkSpectrumZ pkz(pk0,growth,zref);
- 
- //-----------------------------------------------------------------
- cout<<endl<<"\n--- Create mass function"<<endl;
-
- Schechter sch(nstar,mstar,alpha);
- sch.Print();
 
  //-----------------------------------------------------------------
@@ -274 +295 @@
  double sr2int = varpk_int.Variance(R,k1int,k2int);
  cout<<"varpk_int="<<sr2int<<"  ->  sigma="<<sqrt(sr2int)<<endl;
-
- //-----------------------------------------------------------------
- cout<<endl<<"\n--- Compute galaxy density number"<<endl;
+ 
+ //-----------------------------------------------------------------
+ cout<<endl<<"\n--- Create mass function, compute number/mass density, init mass trials"<<endl;
+
+ Schechter sch(nstar,mstar,alpha);
+ sch.Print();
 
  sch.SetOutValue(0);
  cout<<"sch(mstar) = "<<sch(mstar)<<" /Mpc^3/Msol"<<endl;
- double ngal_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
+ double ngal_by_mpc3 = sch.Integrate(schmin,schmax,schnpt);
  cout<<"Galaxy density number = "<<ngal_by_mpc3<<" /Mpc^3 between limits"<<endl;
 
  sch.SetOutValue(1);
  cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
- double mass_by_mpc3 = IntegrateFuncLog(sch,lschmin,lschmax,0.01,dlsch,10.*dlsch,4);
+ double mass_by_mpc3 = sch.Integrate(schmin,schmax,schnpt);
  cout<<"Galaxy mass density= "<<mass_by_mpc3<<" Msol/Mpc^3 between limits"<<endl;
  cout<<"Omega_HI at z=0 is "<<mass_by_mpc3/(univ.Rhoc(0.)*GCm3toMsolMpc3_Cst)<<endl
      <<"         at z="<<zref<<" is "<<mass_by_mpc3/(univ.Rhoc(zref)*GCm3toMsolMpc3_Cst)<<endl;
+
+ SchechterMassDist schmdist(sch,schmin,schmax,schnpt);
+ if(use_schmassdist && schmassdistfile.size()>0) {
+   cout<<"\nWARNING: SchechterMassDist read from "<<schmassdistfile<<endl
+       <<"      PLEASE CHECK CONSISTENCY WITH REQUESTED PARAMETERS"<<endl;
+   schmdist.ReadPPF(schmassdistfile);
+ }
+ schmdist.Print();
+ Histo hmdndm = schmdist.GetHmDnDm();
+ FunRan tirhmdndm = schmdist.GetTmDnDm();
+ {
+ tagobs = "hmdndm"; posobs.PutObject(hmdndm,tagobs);
+ Histo hdum1(tirhmdndm);
+ tagobs = "tirhmdndm"; posobs.PutObject(hdum1,tagobs);
+ }
 
  PrtTim(">>>> End of definition");
@@ -305 +344 @@
  fluct3d.SetGrowthFactor(growth);
  fluct3d.LosComRedshift(0.001,-1);
- TArray<r_8>& rgen = fluct3d.GetRealArray();
+ //TArray<r_8>& rgen = fluct3d.GetRealArray();
  cout<<endl; fluct3d.Print();
+ cout<<"\nMean number of galaxies per pixel = "<<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
+ cout<<"Mean mass per pixel = "<<mass_by_mpc3*fluct3d.GetDVol()<<endl;
 
  double dkmin = fluct3d.GetKincMin();
  double knyqmax = fluct3d.GetKmax();
  long nherr = long(knyqmax/dkmin+0.5);
- cout<<"For HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
+ cout<<"\nFor HistoErr: d="<<dkmin<<" max="<<knyqmax<<" n="<<nherr<<endl;
 
  double dktmin = fluct3d.GetKTincMin();
@@ -421 +462 @@
  cout<<"\n--- Computing a realization in real space"<<endl;
  fluct3d.ComputeReal();
- double rmin,rmax; rgen.MinMax(rmin,rmax);
+ double rmin,rmax; fluct3d.MinMax(rmin,rmax);
  cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
  PrtTim(">>>> End Computing a realization in real space");
@@ -429 +470 @@
    cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
    fluct3d.ApplyGrowthFactor();
-   rgen.MinMax(rmin,rmax);
+   fluct3d.MinMax(rmin,rmax);
    cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
    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(wfits) {
@@ -448 +495 @@
  }
 
- int_8 nm;
- double rm,rs2;
  cout<<"\n--- Check mean and variance in real space"<<endl;
  fluct3d.NumberOfBad(-1.,1e+200);
+ double rm,rs2;
  nm = fluct3d.MeanSigma2(rm,rs2);
  PrtTim(">>>> End Check mean and variance in real space");
@@ -465 +511 @@
  }
 
- //return -41; // CMVTEST
+ //STOP_HERE_FOR QUICK_DEBUG
+ // return -41;
 
  //-----------------------------------------------------------------
@@ -473 +520 @@
  PrtTim(">>>> End Converting fluctuations into mass");
 
- cout<<"\n--- Converting mass into galaxy number"<<endl;
+ cout<<"\n--- Converting mass into galaxy number: gal per pixel ="
+     <<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
  rm = fluct3d.TurnMass2MeanNumber(ngal_by_mpc3);
    nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
@@ -487 +535 @@
    nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
    nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+ double xgalmin,xgalmax; fluct3d.MinMax(xgalmin,xgalmax,0.1,1.e50);
  PrtTim(">>>> End Apply poisson on galaxy number");
  if(wslice) {
@@ -494 +543 @@
 
  cout<<"\n--- Convert Galaxy number to HI mass"<<endl;
- long nhmdndm = long( (lschmax-lschmin+1.)*100. + 0.5);
- Histo hmdndm(lschmin,lschmax,nhmdndm);
- sch.SetOutValue(1);
- FuncToHisto(sch,hmdndm,true);
- FunRan tirhmdndm(hmdndm,true);
- {
- tagobs = "hmdndm"; posobs.PutObject(hmdndm,tagobs);
- Histo hdum1(tirhmdndm);
- tagobs = "tirhmdndm"; posobs.PutObject(hdum1,tagobs);
- }
- double mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
+ double mhi = 0.;
+ if(use_schmassdist) {
+   if(recompute_schmassdist) {
+     int ngalmax = int(xgalmax+0.5);
+     schmdist.SetNgalLim(ngalmax,1,naleagal);
+     PrtTim(">>>> End creating tabulated histograms for trials");
+   }
+   mhi = fluct3d.TurnNGal2MassQuick(schmdist);
+   schmdist.PrintStatus();
+ } else {
+   mhi = fluct3d.TurnNGal2Mass(tirhmdndm,true);
+ }
  cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
    nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
@@ -536 +586 @@
  }
 
- double scalecube = -1., offsetcube=0.;
  if(snoise>0.) {
    cout<<"\n--- Add noise to HI Flux snoise="<<snoise<<", evolution="<<isnoise_evol<<endl;
    fluct3d.AddNoise2Real(snoise,(isnoise_evol>0? true:false));
      nm = fluct3d.MeanSigma2(rm,rs2);
-     if(rs2>0. && sr2int>0.) scalecube = sqrt(sr2int)/sqrt(rs2);
-     offsetcube = -rm;
    PrtTim(">>>> End Add noise");
  }
 
- if(scalecube>0.) {
-   cout<<"\n--- Scale and offset scale="<<scalecube
-       <<" off="<<offsetcube<<endl;
-   fluct3d.ScaleOffset(scalecube,offsetcube);
+ if(scalecube!=0.) {  // Si scalecube==0 pas de normalisation
+   cout<<"\n--- Scale cube rs2ref="<<rs2ref<<endl;
+   if(scalecube<0. && rs2ref>0.) {  // si negatif on scale automatiquement
+     nm = fluct3d.MeanSigma2(rm,rs2);
+     if(rs2>0.) {scalecube=sqrt(rs2ref)/sqrt(rs2); offsetcube=-rm;}
+   }
+   cout<<"...scale="<<scalecube<<" offset="<<offsetcube<<endl;
+   if(scalecube>0.) fluct3d.ScaleOffset(scalecube,offsetcube);
+   PrtTim(">>>> End Scale cube");
  }