Changeset 3150 in Sophya

Timestamp:

Jan 18, 2007, 7:23:56 PM (19 years ago)

Author:

cmv

Message:

suite modifs structure code + bug HistoErr/2D dans ComputeSpectrum cmv 18/01/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvconcherr.cc (modified) (3 diffs)
• cmvobserv3d.cc (modified) (4 diffs)
• genefluct3d.cc (modified) (34 diffs)
• genefluct3d.h (modified) (3 diffs)

trunk/Cosmo/SimLSS/cmvconcherr.cc

@@ -9 +9 @@
 #include "timing.h"
 #include "histerr.h"
-
+#include "hist2err.h"
+#include "fitshisterr.h"
+
+int ObjectType(string nameh,string nameppf);
+int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp);
+int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp);
 void usage(void);
 void usage(void)
 {
- cout<<"cmvconchprof -n name_histoerr file1.ppf file2.ppf ..."<<endl;
-}
-
+cout
+<<"cmvconcherr -w wtyp -n name_histoerr file1.ppf file2.ppf ..."<<endl
+<<"  name_histoerr: object name"<<endl
+<<"  wtyp: bit0 write ASCII, bit1 write PPF, bit2 write FITS (all=7)"<<endl;
+}
+
+
+//---------------------------------------------------------------
 int main(int narg,char *arg[])
 {
- string nameh = ""; // "hpkgen" "hpkgenf" "hpkrec"
+ // "hpkgen" "hpkgenf" "hpkrec" "hpkgen2" "hpkgenf2" "hpkrec2"
+ string nameh = "";
  vector<string> ppfname;
+ int Write_Type = 3;
 
  // --- Decodage arguments
  char c;
- while((c = getopt(narg,arg,"hn:")) != -1) {
+ while((c = getopt(narg,arg,"hn:w:")) != -1) {
   switch (c) {
   case 'n' :
     nameh = optarg;
+    break;
+  case 'w' :
+    sscanf(optarg,"%d",&Write_Type);
     break;
   case 'h' :
@@ -35 +50 @@
  if(nameh.size()<=0 || optind>=narg) {usage(); return -1;}
  for (int i=optind;i<narg;i++) ppfname.push_back(arg[i]);
- cout<<"nameh = "<<nameh<<" , number of file to be read = "<<ppfname.size()<<endl;
-
- // --- lecture et moyenne des HistoErr
+ cout<<"nameh = "<<nameh
+     <<" , number of file to be read = "<<ppfname.size()
+     <<" , write_type="<<Write_Type<<endl;
+
+ // --- Quel type d'objet ?
+ int obtyp = ObjectType(nameh,ppfname[0]);
+
+ // --- lecture et moyenne des HistoErr ou Histo2DErr
+ int nread=0;
+ if(obtyp==1) nread = ConcatHistoErr(nameh,ppfname,Write_Type);
+ else if(obtyp==2) nread = ConcatHisto2DErr(nameh,ppfname,Write_Type);
+ else {
+   cout<<"Type of object "<<nameh<<" not decoded"<<endl;
+   return -1;
+ }
+ cout<<"Number of file read: "<<nread<<endl;
+ return 0;
+}
+
+//---------------------------------------------------------------
+int ObjectType(string nameh,string nameppf)
+// retourne: 1 si HistoErr
+//           2 si Histo2DErr,
+//           0 si autre objet
+//          -1 si il n'y a pas d'objet de nom "nameh"
+{
+ PInPersist pis(nameppf.c_str());
+
+ bool found_tag = pis.GotoNameTag(nameh.c_str());
+ if(!found_tag) return -1;
+
+ PPersist *pps = pis.ReadObject();
+ AnyDataObj *obj = pps->DataObj();
+ cout<<"Object type read from input PPF file : "<<typeid(*obj).name()<< endl;
+
+ HistoErr *herr = dynamic_cast<HistoErr *>(obj);
+ if(herr) return 1;
+
+ Histo2DErr *herr2 = dynamic_cast<Histo2DErr *>(obj);
+ if(herr2) return 2;
+
+ return 0;
+}
+
+//---------------------------------------------------------------
+int ConcatHistoErr(string nameh,vector<string> ppfname,int wrtyp)
+{
+ if(ppfname.size()<=0) return 0;
+
  HistoErr *herrconc = NULL;
- int nherr = 0;
- int itest = 0;
+ int nherr=0, nread=0, itest=0;
  double sum=0., sum2=0., nsum=0;
  for (int ifile=0;ifile<ppfname.size();ifile++) {
@@ -46 +106 @@
    HistoErr herr;
    pis.GetObject(herr,nameh);
-   int n = herr.NBins();
-   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<n<<endl;
+   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="<<herr.NBins()<<endl;
    if(ifile==0) { 
      herrconc = new HistoErr(herr.XMin(),herr.XMax(),herr.NBins());
-     nherr = n;
+     nherr = herr.NBins();
      itest = int(herrconc->NBins()/5.+0.5);
+   } else if(nherr!=herr.NBins()) {
+     cout<<"BAD NUMBER OF BINS"<<endl;
+     continue;
    }
-   if(n!=nherr) {cout<<"BAD NUMBER OF BINS"<<endl; continue;}
+   if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
+                             <<" nmean="<<herr.NMean()<<endl;
+   nread++;
    for(int i=0;i<nherr;i++) herrconc->AddBin(i,herr(i));
    sum+=herr(itest); sum2+=herr(itest)*herr(itest), nsum++;
  }
- herrconc->ToCorrel();
+ herrconc->ToVariance();
  if(nsum>0) {sum/=nsum; sum2/=nsum; sum2-=sum*sum;}
  cout<<"Test bin "<<itest<<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
- cout<<"         "<<" mean="<<(*herrconc)(itest)<<" sigma^2="<<herrconc->Error2(itest)<<endl;
-
- // --- ecriture ascii
- if(1) {
-   FILE * fdata = fopen("cmvconchprof.data","w");
-   fprintf(fdata,"# k(Mpc^-1) <pkz_reconstruit> variance_pkz\n");
+ cout<<"         "<<" mean="<<(*herrconc)(itest)
+                  <<" sigma^2="<<herrconc->Error2(itest)<<endl;
+
+ // --- ecriture
+ if(wrtyp&1) {   // ecriture ASCII
+   FILE * fdata = fopen("cmvconcherr.data","w");
+   fprintf(fdata,"# x <mean> variance\n");
+   fprintf(fdata,"%.17e %.17e %d %d\n",herrconc->XMin(),herrconc->XMax()
+          ,herrconc->NBins(),herrconc->NMean());
    for(int i=0;i<herrconc->NBins();i++) {
-     double k = herrconc->BinCenter(i);
-     fprintf(fdata,"%e %e %e\n",k,(*herrconc)(i),herrconc->Error2(i));
+     double x = herrconc->BinCenter(i);
+     fprintf(fdata,"%e %e %e\n",x,(*herrconc)(i),herrconc->Error2(i));
    }
    fclose(fdata);
  }
-
- // --- ecriture ppf
- string tagobs = "cmvconchprof.ppf";
- POutPersist posobs(tagobs);
- tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
+ if(wrtyp&2) {   // ecriture PPF
+   string tagobs = "cmvconcherr.ppf";
+   POutPersist posobs(tagobs);
+   tagobs = "herrconc"; posobs.PutObject(*herrconc,tagobs);
+ }
+ if(wrtyp&4) {   // ecriture FITS
+   FitsInOutFile fio("!cmvconcherr.fits",FitsInOutFile::Fits_Create);
+   fio << *herrconc;
+ }
 
  delete herrconc;
 
- return 0;
+ return nread;
+}
+
+//---------------------------------------------------------------
+int ConcatHisto2DErr(string nameh,vector<string> ppfname,int wrtyp)
+{
+ if(ppfname.size()<=0) return 0;
+
+ Histo2DErr *herrconc = NULL;
+ int nherrx=0, nherry=0, nread=0, itestx=0, itesty=0;
+ double sum=0., sum2=0., nsum=0;
+ for (int ifile=0;ifile<ppfname.size();ifile++) {
+   PInPersist pis(ppfname[ifile].c_str());
+   Histo2DErr herr;
+   pis.GetObject(herr,nameh);
+   cout<<ifile<<" : "<<ppfname[ifile]<<"  nbin="
+       <<herr.NBinX()<<","<<herr.NBinY()<<endl;
+   if(ifile==0) { 
+     herrconc = new Histo2DErr(herr.XMin(),herr.XMax(),herr.NBinX()
+                              ,herr.YMin(),herr.YMax(),herr.NBinY());
+     nherrx = herr.NBinX(); nherry = herr.NBinY();
+     itestx = int(herrconc->NBinX()/5.+0.5);
+     itesty = int(herrconc->NBinY()/5.+0.5);
+   } else if(nherrx!=herr.NBinX() || nherry!=herr.NBinY()) {
+     cout<<"BAD NUMBER OF BINS"<<endl;
+     continue;
+   }
+   if(herr.NMean()!=1) cout<<"ATTENTION: file="<<ppfname[ifile]
+                             <<" nmean="<<herr.NMean()<<endl;
+   nread++;
+   for(int i=0;i<nherrx;i++)
+     for(int j=0;j<nherrx;j++) herrconc->AddBin(i,j,herr(i,j));
+   sum+=herr(itestx,itesty); sum2+=herr(itestx,itesty)*herr(itestx,itesty), nsum++;
+ }
+ herrconc->ToVariance();
+ cout<<"Test bin "<<itestx<<","<<itesty
+     <<" mean="<<sum<<" sigma^2="<<sum2<<"  nsum="<<nsum<<endl;
+ cout<<"         "<<" mean="<<(*herrconc)(itestx,itesty)
+                  <<" sigma^2="<<herrconc->Error2(itestx,itesty)<<endl;
+
+ // --- ecriture
+ if(wrtyp&1) {   // ecriture ASCII
+   cout<<"ERROR: No ASCII writing implemented"<<endl;
+ }
+ if(wrtyp&2) {   // ecriture PPF
+   string tagobs = "cmvconcherr2.ppf";
+   POutPersist posobs(tagobs);
+   tagobs = "herrconc2"; posobs.PutObject(*herrconc,tagobs);
+ }
+ if(wrtyp&4) {   // ecriture FITS
+   FitsInOutFile fio("!cmvconcherr2.fits",FitsInOutFile::Fits_Create);
+   fio << *herrconc;
+ }
+
+ delete herrconc;
+
+ return nread;
 }
 
 /*
-openppf cmvconchprof.ppf
-
-disp herrconc
+#### HistoErr 1D
+openppf cmvconcherr.ppf
+
+disp herrconc "hbincont err"
+disp herrconc "hbinerr"
+disp herrconc "hbinent"
 
 n/plot herrconc.val%log10(x) x>0 ! "connectpoints"
-n/plot herrconc.err%log10(x) x>0 ! "connectpoints same red"
-
-n/plot herrconc.err/val%log10(x) x>0&&val>0 ! "connectpoints"
+n/plot herrconc.sqrt(err2)%log10(x) x>0&&err2>0. ! "connectpoints same red"
+
+n/plot herrconc.sqrt(err2)/val%log10(x) x>0&&err2>0.&&val>0 ! "connectpoints"
+
+#### Histo2DErr 2D
+openppf cmvconcherr2.ppf
+
+disp herrconc2 "hbincont"
+disp herrconc2 "hbinerr"
+disp herrconc2 "hbinent"
 */

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -204 +204 @@
  TArray< complex<r_8> > pkgen;
  GeneFluct3D fluct3d(pkgen);
+ fluct3d.SetPrtLevel(3);
  fluct3d.SetNThread(nthread);
  fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
@@ -431 +432 @@
 
 zone
+set k pow(10.,x)
+n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
+
+zone
 n/plot hpkz.val%x ! ! "nsta connectpoints"
 n/plot hpkgen.val%log10(x) x>0 ! "nsta same red connectpoints"
@@ -436 +441 @@
 n/plot hpkrec.val%log10(x) x>0 ! "nsta same blue connectpoints"
 
-set k pow(10.,x)
-n/plot hpkz.val*$k*$k/(2*M_PI*M_PI)%x ! "connectpoints"
+disp hpkgen "hbincont err"
+disp hpkgenf "hbincont err"
+disp hpkrec "hbincont err"
 
 zone 2 2
@@ -444 +450 @@
 imag hpkrec2
 
-zone
-disp hmdndm
-disp tirhmdndm
+zone 2 1
+disp hmdndm "nsta"
+disp tirhmdndm "nsta"
 addline 0 1 20 1 "red"
 

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -32 +32 @@
 //-------------------------------------------------------
 GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
-  : T_(T) , array_allocated_(false)
+  : T_(T) , array_allocated_(false) , lp_(0)
 {
  SetNThread();
@@ -257 +257 @@
 //    (attention on fait une fft pour realiser le spectre)
 {
- int lp=2;
-
  // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
- if(lp>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: before fftw_plan ---");
 #ifdef FFTW_THREAD
  if(nthread_>0) {
@@ -270 +268 @@
  pf_ = fftw_plan_dft_r2c_3d(Nx_,Ny_,Nz_,data_,fdata_,FFTW_ESTIMATE);
  pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
- if(lp>1) PrtTim("--- ComputeFourier0: after fftw_plan ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after fftw_plan ---");
 
  // --- realisation d'un tableau de tirage gaussiens
- if(lp>1) PrtTim("--- ComputeFourier0: before gaussian filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: before gaussian filling ---");
  // On tient compte du pb de normalisation de FFTW3
  double sntot = sqrt((double)NRtot_);
@@ -280 +278 @@
    data_[ip] = NorRand()/sntot;
  }
- if(lp>1) PrtTim("--- ComputeFourier0: after gaussian filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after gaussian filling ---");
 
  // --- realisation d'un tableau de tirage gaussiens
- if(lp>1) PrtTim("--- ComputeFourier0: before fft real ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: before fft real ---");
  fftw_execute(pf_);
- if(lp>1) PrtTim("--- ComputeFourier0: after fft real ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after fft real ---");
 
  // --- On remplit avec une realisation
- if(lp>1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: before Fourier realization filling ---");
  T_(0,0,0) = complex<r_8>(0.);  // on coupe le continue et on l'initialise
  long lmod = Nx_/10; if(lmod<1) lmod=1;
@@ -294 +292 @@
    long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_;  kx *= kx;
-   if(lp>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   if(lp_>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
    for(long j=0;j<Ny_;j++) {
      long jj = (j>Ny_/2) ? Ny_-j : j;
@@ -309 +307 @@
    }
  }
- if(lp>1) PrtTim("--- ComputeFourier0: after Fourier realization filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after Fourier realization filling ---");
 
  double p = compute_power_carte();
  cout<<"Puissance dans la realisation: "<<p<<endl;
- if(lp>1) PrtTim("--- ComputeFourier0: after Computing power ---");
+ if(lp_>1) PrtTim("--- ComputeFourier0: after Computing power ---");
 
 }
@@ -328 +326 @@
 //                                          sum(fb(x_i)^2) = S*N^2
 {
- int lp=2;
-
  // --- Initialisation de fftw3 (attention data est sur-ecrit a l'init)
- if(lp>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: before fftw_plan ---");
 #ifdef FFTW_THREAD
  if(nthread_>0) {
@@ -342 +338 @@
  pb_ = fftw_plan_dft_c2r_3d(Nx_,Ny_,Nz_,fdata_,data_,FFTW_ESTIMATE);
  //fftw_print_plan(pb_); cout<<endl;
- if(lp>1) PrtTim("--- ComputeFourier: after fftw_plan ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: after fftw_plan ---");
 
  // --- RaZ du tableau
@@ -348 +344 @@
 
  // --- On remplit avec une realisation
- if(lp>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: before Fourier realization filling ---");
  long lmod = Nx_/10; if(lmod<1) lmod=1;
  for(long i=0;i<Nx_;i++) {
    long ii = (i>Nx_/2) ? Nx_-i : i;
    double kx = ii*Dkx_;  kx *= kx;
-   if(lp>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
+   if(lp_>1 && i%lmod==0) cout<<"i="<<i<<" ii="<<ii<<endl;
    for(long j=0;j<Ny_;j++) {
      long jj = (j>Ny_/2) ? Ny_-j : j;
@@ -369 +365 @@
    }
  }
- if(lp>1) PrtTim("--- ComputeFourier: after Fourier realization filling ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: after Fourier realization filling ---");
 
  manage_coefficients();   // gros effet pour les spectres que l'on utilise !
- if(lp>1) PrtTim("--- ComputeFourier: after managing coefficients ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: after managing coefficients ---");
 
  double p = compute_power_carte();
  cout<<"Puissance dans la realisation: "<<p<<endl;
- if(lp>1) PrtTim("--- ComputeFourier: after before Computing power ---");
+ if(lp_>1) PrtTim("--- ComputeFourier: after before Computing power ---");
 
 }
@@ -490 +486 @@
 //                          avec y = k_x*dx/2
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- FilterByPixel: before ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- FilterByPixel: before ---");
 
  for(long i=0;i<Nx_;i++) {
@@ -511 +506 @@
  }
 
- if(lp>1) PrtTim("--- FilterByPixel: after ---");
+ if(lp_>1) PrtTim("--- FilterByPixel: after ---");
 }
 
@@ -518 +513 @@
 // Calcule une realisation dans l'espace reel
 {
- int lp=2;
  check_array_alloc();
 
  // On fait la FFT
- if(lp>1) PrtTim("--- ComputeReal: before fftw backward ---");
+ if(lp_>1) PrtTim("--- ComputeReal: before fftw backward ---");
  fftw_execute(pb_);
- if(lp>1) PrtTim("--- ComputeReal: after fftw backward ---");
+ if(lp_>1) PrtTim("--- ComputeReal: after fftw backward ---");
 }
 
@@ -530 +524 @@
 void GeneFluct3D::ReComputeFourier(void)
 {
- int lp=2;
  check_array_alloc();
 
  // On fait la FFT
- if(lp>1) PrtTim("--- ComputeReal: before fftw forward ---");
+ if(lp_>1) PrtTim("--- ComputeReal: before fftw forward ---");
  fftw_execute(pf_);
  // On corrige du pb de la normalisation de FFTW3
@@ -542 +535 @@
      for(long l=0;l<NCz_;l++) T_(l,j,i) /= complex<r_8>(v);
 
- if(lp>1) PrtTim("--- ComputeReal: after fftw forward ---");
+ if(lp_>1) PrtTim("--- ComputeReal: after fftw forward ---");
 }
 
@@ -552 +545 @@
 {
  check_array_alloc();
+ if(lp_>1) PrtTim("--- ComputeSpectrum: before ---");
 
  if(herr.NBins()<0) return -1;
@@ -571 +565 @@
    }
  }
- herr.ToCorrel();
+ herr.ToVariance();
 
  // renormalize to directly compare to original spectrum
@@ -577 +571 @@
  herr *= norm;
 
+ if(lp_>1) PrtTim("--- ComputeSpectrum: after ---");
+
  return 0;
 }
@@ -583 +579 @@
 {
  check_array_alloc();
+ if(lp_>1) PrtTim("--- ComputeSpectrum2D: before ---");
 
  if(herr.NBinX()<0 || herr.NBinY()<0) return -1;
@@ -602 +599 @@
    }
  }
- herr.ToCorrel();
+ herr.ToVariance();
 
  // renormalize to directly compare to original spectrum
@@ -608 +605 @@
  herr *= norm;
 
+ if(lp_>1) PrtTim("--- ComputeSpectrum2D: after ---");
+
  return 0;
 }
@@ -615 +614 @@
 // Recompute MASS variance in spherical top-hat (rayon=R)
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- VarianceFrReal: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- VarianceFrReal: before computation ---");
 
  long dnx = long(R/Dx_+0.5); if(dnx<=0) dnx = 1;
@@ -654 +652 @@
  sum /= nsum;
  sum2 = sum2/nsum - sum*sum;
- if(lp) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
+ if(lp_) cout<<"VarianceFrReal: nsum="<<nsum<<" <M>="<<sum<<" <(M-<M>)^2>="<<sum2<<endl;
 
  var = sum2/(sum*sum);  // <dM>^2/<M>^2
- if(lp) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
-
- if(lp>1) PrtTim("--- VarianceFrReal: after computation ---");
+ if(lp_) cout<<"sigmaR^2="<<var<<" -> "<<sqrt(var)<<endl;
+
+ if(lp_>1) PrtTim("--- VarianceFrReal: after computation ---");
  return nsum;
 }
@@ -726 +724 @@
 // d_rho/rho -> Mass  (add one!)
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- TurnFluct2Mass: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- TurnFluct2Mass: before computation ---");
 
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
@@ -736 +733 @@
  }
 
- if(lp>1) PrtTim("--- TurnFluct2Mass: after computation ---");
+ if(lp_>1) PrtTim("--- TurnFluct2Mass: after computation ---");
 }
 
@@ -742 +739 @@
 // do NOT treate negative or nul values
 {
- int lp=2;
- if(lp>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");
+ if(lp_>1) PrtTim("--- TurnMass2MeanNumber: before computation ---");
 
  double m,s2;
  int_8 ngood = MeanSigma2(m,s2,0.,1e+200);
- if(lp) cout<<"TurnMass2MeanNumber: ngood="<<ngood
+ if(lp_) cout<<"TurnMass2MeanNumber: ngood="<<ngood
             <<" m="<<m<<" s2="<<s2<<" -> "<<sqrt(s2)<<endl;
 
@@ -758 +754 @@
  //   (rappel sur tout les pixels <1+d_rho/rho>=1)
  double dn = n_by_mpc3*Vol_/m /(double)ngood;  // nb de gal a mettre ds 1 px
- if(lp) cout<<"galaxy density move from "
+ if(lp_) cout<<"galaxy density move from "
             <<n_by_mpc3*Vol_/double(NRtot_)<<" to "<<dn<<" / pixel"<<endl;
  double sum = 0.;
@@ -766 +762 @@
    if(data_[ip]>0.) sum += data_[ip];  // mais on ne les compte pas
  }
- if(lp) cout<<sum<<" galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
-
- if(lp>1) PrtTim("--- TurnMass2MeanNumber: after computation ---");
+ if(lp_) cout<<sum<<" galaxies put into survey / "<<n_by_mpc3*Vol_<<endl;
+
+ if(lp_>1) PrtTim("--- TurnMass2MeanNumber: after computation ---");
  return sum;
 }
@@ -775 +771 @@
 // do NOT treate negative or nul mass  -> let it as it is
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");
 
  double sum = 0.;
@@ -790 +785 @@
    }
  }
- if(lp) cout<<sum<<" galaxies put into survey"<<endl;
-
- if(lp>1) PrtTim("--- ApplyPoisson: before computation ---");
+ if(lp_) cout<<sum<<" galaxies put into survey"<<endl;
+
+ if(lp_>1) PrtTim("--- ApplyPoisson: before computation ---");
  return sum;
 }
@@ -804 +799 @@
 // RETURN la masse totale
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- TurnNGal2Mass: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- TurnNGal2Mass: before computation ---");
 
  double sum = 0.;
@@ -824 +818 @@
    }
  }
- if(lp) cout<<sum<<" MSol HI mass put into survey"<<endl;
-
- if(lp>1) PrtTim("--- TurnNGal2Mass: after computation ---");
+ if(lp_) cout<<sum<<" MSol HI mass put into survey"<<endl;
+
+ if(lp_>1) PrtTim("--- TurnNGal2Mass: after computation ---");
  return sum;
 }
@@ -833 +827 @@
 // add noise to every pixels (meme les <=0 !)
 {
- int lp=2;
- check_array_alloc();
-
- if(lp>1) PrtTim("--- AddNoise2Real: before computation ---");
+ check_array_alloc();
+
+ if(lp_>1) PrtTim("--- AddNoise2Real: before computation ---");
 
  for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
@@ -843 +836 @@
  }
 
- if(lp>1) PrtTim("--- AddNoise2Real: after computation ---");
-}
+ if(lp_>1) PrtTim("--- AddNoise2Real: after computation ---");
+}

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -85 +85 @@
   void ReadPPF(string cfname);
 
+  void SetPrtLevel(int lp=0) {lp_ = lp;}
   void Print(void);
 
@@ -96 +97 @@
   inline double pixelfilter(double x)
     {return (x<0.025) ? 1.-x*x/6.*(1.-x*x/20.): sin(x)/x;}
-
 
   long Nx_,Ny_,Nz_;  vector<long> N_;
@@ -111 +111 @@
   fftw_plan pf_,pb_;
   unsigned short nthread_;
+  int lp_;
 
   bool array_allocated_;  // true if array has been allocated