Index: trunk/Cosmo/SimLSS/cmvobserv3d.cc
===================================================================
--- trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3345)
+++ trunk/Cosmo/SimLSS/cmvobserv3d.cc	(revision 3349)
@@ -45,7 +45,4 @@
      <<" -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
-     <<"    powlaw: apply S_agn evolution as (Nu/1.4)^powlaw"<<endl
      <<" -W : write cube in FITS format (complex cube is coded as real cube)"<<endl
      <<" -P : write cube in PPF format"<<endl
@@ -54,4 +51,7 @@
      <<" -T nth : nombre de threads (si compil multi-thread, default: 0)"<<endl
      <<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
+     ////<<"    powlaw: apply S_agn evolution as (Nu/1.4)^powlaw"<<endl
 }
 
@@ -101,8 +101,8 @@
  int noise_evol = 0;
 
- // *** AGN
- bool do_agn = false;
- double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
- double powlaw_agn = 0.;
+ //// *** AGN
+ ////bool do_agn = false;
+ ////double lfjy_agn=-99., lsigma_agn=0.;   // en Jy
+ ////double powlaw_agn = 0.;
 
  // *** type de generation
@@ -177,8 +177,8 @@
     schmassdistfile = optarg;
     break;
-  case 'A' :
-    do_agn = true;
-    sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
-    break;
+    ////  case 'A' :
+    ////do_agn = true;
+    ////sscanf(optarg,"%lf,%lf,%lf",&lfjy_agn,&lsigma_agn,&powlaw_agn);
+    ////break;
   case 'V' :
     compvarreal = true;
@@ -224,7 +224,7 @@
  cout<<"snoise="<<snoise<<" equivalent Msol, evolution="<<noise_evol<<endl;
  cout<<"scalecube="<<scalecube<<", offsetcube="<<offsetcube<<endl;
- if(do_agn)
-   cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
-       <<" , powlaw="<<powlaw_agn<<endl;
+ ////if(do_agn)
+ ////  cout<<"AGN: <log10(Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
+ ////      <<" , powlaw="<<powlaw_agn<<endl;
 
  string tagobs = "cmvobserv3d.ppf";
@@ -264,10 +264,10 @@
  cout<<endl<<"\n--- Compute variance for top-hat R="<<R
      <<" at z="<<pkz.GetZ()<<endl;
- VarianceSpectrum varpk_th(pkz,0);
- double kfind_th = varpk_th.FindMaximum(R,kmin,kmax,eps);
+ 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(R,pkmax_th/1.e4,k1,k2,eps);
+ int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
  cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
      <<k2<<" ("<<log10(k2)<<")"<<endl;
@@ -275,5 +275,5 @@
  double ldlk = (log10(k2)-log10(k1))/npt;
  varpk_th.SetInteg(0.01,ldlk,-1.,4);
- double sr2 = varpk_th.Variance(R,k1,k2);
+ double sr2 = varpk_th.Variance(k1,k2);
  cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;
 
@@ -292,11 +292,11 @@
  //-----------------------------------------------------------------
  cout<<endl<<"\n--- Compute variance for Pk at z="<<pkz.GetZ()<<endl;
- VarianceSpectrum varpk_int(pkz,2);
-
- double kfind_int = varpk_int.FindMaximum(R,kmin,kmax,eps);
+ VarianceSpectrum varpk_int(pkz,R,VarianceSpectrum::NOFILTER);
+
+ double kfind_int = varpk_int.FindMaximum(kmin,kmax,eps);
  double pkmax_int = varpk_int(kfind_int);
  cout<<"kfind_int = "<<kfind_int<<" ("<<log10(kfind_int)<<"), integrand="<<pkmax_int<<endl;
  double k1int=kmin, k2int=kmax;
- int rcint = varpk_int.FindLimits(R,pkmax_int/1.e4,k1int,k2int,eps);
+ int rcint = varpk_int.FindLimits(pkmax_int/1.e4,k1int,k2int,eps);
  cout<<"limit_int: rc="<<rcint<<" : "<<k1int<<" ("<<log10(k1int)<<") , "
      <<k2int<<" ("<<log10(k2int)<<")"<<endl;
@@ -304,5 +304,5 @@
  double ldlkint = (log10(k2int)-log10(k1int))/npt;
  varpk_int.SetInteg(0.01,ldlkint,-1.,4);
- double sr2int = varpk_int.Variance(R,k1int,k2int);
+ double sr2int = varpk_int.Variance(k1int,k2int);
  cout<<"varpk_int="<<sr2int<<"  ->  sigma="<<sqrt(sr2int)<<endl;
  
@@ -355,14 +355,11 @@
  cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
 
- TArray< complex<r_8> > pkgen;
- GeneFluct3D fluct3d(pkgen);
- fluct3d.SetPrtLevel(2);
- fluct3d.SetNThread(nthread);
- fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
+ GeneFluct3D fluct3d(nx,ny,nz,dx,dy,dz,nthread,2);
  fluct3d.SetObservator(zref,-nz/2.);
  fluct3d.SetCosmology(univ);
  fluct3d.SetGrowthFactor(growth);
  fluct3d.LosComRedshift(0.001,-1);
- //TArray<r_8>& rgen = fluct3d.GetRealArray();
+ TArray< complex<r_8> >& pkgen = fluct3d.GetComplexArray();
+ TArray<r_8>& rgen = fluct3d.GetRealArray();
  cout<<endl; fluct3d.Print();
  cout<<"\nMean number of galaxies per pixel = "<<ngal_by_mpc3*fluct3d.GetDVol()<<endl;
@@ -391,5 +388,5 @@
  ldlkint = (log10(knyqmax_mod)-log10(k1int))/npt;
  varpk_int.SetInteg(0.01,ldlkint,-1.,4);
- double sr2int_kmax = varpk_int.Variance(R,k1int,knyqmax_mod);
+ double sr2int_kmax = varpk_int.Variance(k1int,knyqmax_mod);
  cout<<"varpk_int(<"<<knyqmax_mod<<")="<<sr2int_kmax<<"  ->  sigma="<<sqrt(sr2int_kmax)<<endl;
 
@@ -538,29 +535,22 @@
 
  if(no_poisson) {
-
    cout<<"\n--- Converting !!!DIRECTLY!!! mass into HI mass: mass per pixel ="
        <<mass_by_mpc3*fluct3d.GetDVol()<<endl;
-   rm = fluct3d.TurnMass2HIMass(mass_by_mpc3);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
-   PrtTim(">>>> End Converting !!!DIRECTLY!!! mass into HI mass");
-
+   rm = fluct3d.TurnMass2MeanNumber(mass_by_mpc3);
  } else {
-
    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);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
-   PrtTim(">>>> End Converting mass into galaxy number");
-
-   cout<<"\n--- Set negative and null pixels to BAD"<<endl;
-   nm = fluct3d.SetToVal(0.,1e+200,-999.);
-   PrtTim(">>>> End Set negative pixels to BAD etc...");
+ }
+ nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
+ nm = fluct3d.MeanSigma2(rm,rs2);
+ PrtTim(">>>> End Converting mass into galaxy number or mass");
+
+ if( !no_poisson ) {
 
    cout<<"\n--- Apply poisson on galaxy number"<<endl;
    fluct3d.ApplyPoisson();
-     nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
+     nm = fluct3d.MeanSigma2(rm,rs2);
    double xgalmin,xgalmax; fluct3d.MinMax(xgalmin,xgalmax,0.1,1.e50);
    PrtTim(">>>> End Apply poisson on galaxy number");
@@ -584,12 +574,8 @@
    }
    cout<<mhi<<" MSol in survey / "<<mass_by_mpc3*fluct3d.GetVol()<<endl;
-     nm = fluct3d.MeanSigma2(rm,rs2,-998.,1e200);
+     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200);
      cout<<"Equivalent: "<<rm*nm/fluct3d.NPix()<<" Msol / pixels"<<endl;
-     nm = fluct3d.MeanSigma2(rm,rs2,0.,1e200,true,0.);
+     nm = fluct3d.MeanSigma2(rm,rs2);
    PrtTim(">>>> End Convert Galaxy number to HI mass");
-
-   cout<<"\n--- Set BAD pixels to Zero"<<endl;
-   nm = fluct3d.SetToVal(-998.,1e+200,0.);
-   PrtTim(">>>> End Set BAD pixels to Zero etc...");
 
  }
@@ -609,11 +595,11 @@
 
  //-----------------------------------------------------------------
- if(do_agn) {
-   cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
-       <<" , powlaw="<<powlaw_agn<<endl;
-   fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
-     nm = fluct3d.MeanSigma2(rm,rs2);
-   PrtTim(">>>> End Add AGN");
- }
+ ////if(do_agn) {
+ ////  cout<<"\n--- Add AGN: <log10(S Jy)>="<<lfjy_agn<<" , sigma="<<lsigma_agn
+ ////      <<" , powlaw="<<powlaw_agn<<endl;
+ ////  fluct3d.AddAGN(lfjy_agn,lsigma_agn,powlaw_agn);
+ ////    nm = fluct3d.MeanSigma2(rm,rs2);
+ ////   PrtTim(">>>> End Add AGN");
+ ////}
 
  //-----------------------------------------------------------------
@@ -627,4 +613,6 @@
  }
 
+
+ //-----------------------------------------------------------------
  if(scalecube!=0.) {  // Si scalecube==0 pas de normalisation
    cout<<"\n--- Scale cube rs2ref="<<rs2ref<<endl;
