Changeset 3157 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Jan 25, 2007, 6:04:48 PM (19 years ago)

Author:

cmv

Message:

intro du facteur de croissance dans la simul cmv 25/01/2007

Location:

trunk/Cosmo/SimLSS

Files:

• Makefile (modified) (3 diffs)
• cmvobserv3d.cc (modified) (13 diffs)
• cmvtinterp.cc (added)
• genefluct3d.cc (modified) (5 diffs)
• genefluct3d.h (modified) (4 diffs)
• geneutils.cc (modified) (4 diffs)
• geneutils.h (modified) (5 diffs)

trunk/Cosmo/SimLSS/Makefile

@@ -19 +19 @@
        $(EXE)cmvtstsch $(EXE)cmvtstblack $(EXE)cmvtvarspec \
        $(EXE)cmvdefsurv $(EXE)cmvobserv3d $(EXE)cmvtintfun \
-       $(EXE)cmvtpoisson $(EXE)cmvconcherr
+       $(EXE)cmvtpoisson $(EXE)cmvconcherr $(EXE)cmvtinterp
 #$(EXE)cmvtluc
  
@@ -25 +25 @@
           $(OBJ)cmvtuniv.o $(OBJ)cmvtransf.o $(OBJ)cmvtgrowth.o $(OBJ)cmvtstpk.o \
           $(OBJ)cmvtstsch.o $(OBJ)cmvtstblack.o $(OBJ)cmvtvarspec.o $(OBJ)cmvdefsurv.o \
-          $(OBJ)cmvobserv3d.o $(OBJ)cmvtintfun.o \
+          $(OBJ)cmvobserv3d.o $(OBJ)cmvtintfun.o $(OBJ)cmvtinterp.o \
           $(OBJ)cmvtpoisson.o $(OBJ)cmvconcherr.o $(OBJ)cmvtluc.o
  
@@ -165 +165 @@
 
 ##############################################################################
+cmvtinterp: $(EXE)cmvtinterp
+        echo $@ " done"
+$(EXE)cmvtinterp: $(OBJ)cmvtinterp.o $(LIB)libcmvsimbao.a
+        $(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvtinterp.o $(MYLIB) 
+$(OBJ)cmvtinterp.o: cmvtinterp.cc
+        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtinterp.cc
+
+##############################################################################
 cmvtluc: $(EXE)cmvtluc
         echo $@ " done"

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -12 +12 @@
 
 #include "constcosmo.h"
+#include "cosmocalc.h"
 #include "schechter.h"
 #include "geneutils.h"
@@ -23 +24 @@
      <<" -a : init auto de l aleatoire"<<endl
      <<" -0 : use methode ComputeFourier0"<<endl
+     <<" -G : compute Pk(z=0) and apply growth factor in real space"<<endl
      <<" -x nx,dx : taille en x (npix,Mpc)"<<endl
      <<" -y ny,dy : taille en y (npix,Mpc)"<<endl
@@ -46 +48 @@
 
  // *** Cosmography definition   (WMAP)
+ unsigned short flat = 0;
  double ob0 = 0.0444356;
  double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
  double zref = 0.5;
+ double perc=0.01,dzinc=-1.,dzmax=5.; unsigned short glorder=4;
 
  // *** Spectrum and variance definition
@@ -75 +79 @@
  // *** type de generation
  bool computefourier0=false;
+ bool use_growth_factor = false;
  unsigned short nthread=4;
 
@@ -86 +91 @@
 
  char c;
- while((c = getopt(narg,arg,"ha0PWV2x:y:z:s:Z:")) != -1) {
+ while((c = getopt(narg,arg,"ha0PWV2Gx:y:z:s:Z:")) != -1) {
   switch (c) {
   case 'a' :
@@ -93 +98 @@
   case '0' :
     computefourier0 = true;
+    break;
+  case 'G' :
+    use_growth_factor = true;
     break;
   case 'x' :
@@ -142 +150 @@
 
  //-----------------------------------------------------------------
- cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
+ cout<<endl<<"\n--- Create Cosmology"<<endl;
+
+ CosmoCalc univ(flat,true,zref+1.);
+ univ.SetInteg(perc,dzinc,dzmax,glorder);
+ univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.Print();
+ double loscomref = univ.Dloscom(zref);
+ cout<<"zref = "<<zref<<" -> dloscom = "<<loscomref<<" Mpc"<<endl;
+
+ //-----------------------------------------------------------------
+  cout<<endl<<"\n--- Create Spectrum and mass function"<<endl;
 
  InitialSpectrum pkini(ns,as);
@@ -149 +167 @@
  //tf.SetNoOscEnv(2);
 
- GrowthFactor d1(om0,ol0);
+ GrowthFactor growth(om0,ol0);
 
  PkSpectrum0 pk0(pkini,tf);
 
- PkSpectrumZ pkz(pk0,d1,zref);
+ PkSpectrumZ pkz(pk0,growth,zref);
  
  Schechter sch(nstar,mstar,alpha);
@@ -223 +241 @@
  cout<<endl<<"\n--- Initialisation de GeneFluct3D"<<endl;
 
- pkz.SetZ(zref);
  TArray< complex<r_8> > pkgen;
  GeneFluct3D fluct3d(pkgen);
@@ -229 +246 @@
  fluct3d.SetNThread(nthread);
  fluct3d.SetSize(nx,ny,nz,dx,dy,dz);
+ fluct3d.SetObservator(zref,nz/2.);
+ fluct3d.SetCosmology(univ);
+ fluct3d.SetGrowthFactor(growth);
+ fluct3d.LosComRedshift(0.001);
  TArray<r_8>& rgen = fluct3d.GetRealArray();
- fluct3d.Print();
+ cout<<endl; fluct3d.Print();
 
  double dkmin = fluct3d.GetKincMin();
@@ -246 +267 @@
      <<" max="<<ktnyqmax<<","<<kznyqmax<<" n="<<nherrt<<","<<nherrz<<endl;
 
+ //-----------------------------------------------------------------
  cout<<"\n--- Computing spectra variance up to Kmax at z="<<pkz.GetZ()<<endl;
  // En fait on travaille sur un cube inscrit dans la sphere de rayon kmax:
@@ -258 +280 @@
  PrtTim(">>>> End Initialisation de GeneFluct3D");
 
+ //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in Fourier space"<<endl;
+ if(use_growth_factor) pkz.SetZ(0.); else pkz.SetZ(zref);
+ cout<<"Power spectrum set at redshift: "<<pkz.GetZ()<<endl;
  if(computefourier0) fluct3d.ComputeFourier0(pkz);
    else              fluct3d.ComputeFourier(pkz);
@@ -326 +351 @@
  }
 
+ //-----------------------------------------------------------------
  cout<<"\n--- Computing a realization in real space"<<endl;
  fluct3d.ComputeReal();
  double rmin,rmax; rgen.MinMax(rmin,rmax);
  cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
-   PrtTim(">>>> End Computing a realization in real space");
+ PrtTim(">>>> End Computing a realization in real space");
+
+ if(use_growth_factor) {
+   cout<<"\n--- Apply Growth factor"<<endl;
+   cout<<"...D(z=0)="<<growth(0.)<<"  D(z="<<zref<<")="<<growth(zref)<<endl;
+   fluct3d.ApplyGrowthFactor(-1);
+   rmin,rmax; rgen.MinMax(rmin,rmax);
+   cout<<"rgen.Min = "<<rmin<<" , Max="<<rmax<<endl;
+   PrtTim(">>>> End Applying growth factor");
+ }
 
  if(wfits) {

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -32 +32 @@
 GeneFluct3D::GeneFluct3D(TArray< complex<r_8 > >& T)
   : T_(T) , Nx_(0) , Ny_(0) , Nz_(0) , array_allocated_(false) , lp_(0)
-  , redshref_(-999.) , kredshref_(-999.)
-{
+  , redshref_(-999.) , kredshref_(0.) , cosmo_(NULL) , growth_(NULL)
+  , loscom_ref_(-999.), loscom_min_(-999.), loscom_max_(-999.)
+
+
+{
+ xobs_[0] = xobs_[1] = xobs_[2] = 0.;
+ zred_.resize(0);
+ loscom_.resize(0);
  SetNThread();
 }
@@ -64 +70 @@
 //     kredshref = indice (en double) correspondant a ce redshift
 //                 Si kredshref<0 alors kredshref=0
+// Exemple: redshref=1.5 kredshref=250.75
+//    -> Le pixel i=nx/2 j=ny/2 k=250.75 est au redshift 1.5
 {
  if(redshref<0.) redshref = 0.;
  if(kredshref<0.) kredshref = 0.;
- redshref_ = redshref;
+ redshref_  = redshref;
  kredshref_ = kredshref;
+}
+
+void GeneFluct3D::SetCosmology(CosmoCalc& cosmo)
+{
+ cosmo_ = &cosmo;
+ if(lp_>1) cosmo_->Print();
+}
+
+void GeneFluct3D::SetGrowthFactor(GrowthFactor& growth)
+{
+ growth_ = &growth;
 }
 
@@ -164 +183 @@
 }
 
+//-------------------------------------------------------
+long GeneFluct3D::LosComRedshift(double zinc)
+// Given a position of the cube relative to the observer
+// and a cosmology
+// (SetObservator() and SetCosmology() should have been called !)
+// This routine filled:
+//   the vector "zred_" of scanned redshift (by zinc increments)
+//   the vector "loscom_" of corresponding los comoving distance
+// 
+{
+ if(zinc<=0.) zinc = 0.01;
+ if(lp_>0) cout<<"--- LosComRedshift: zinc="<<zinc<<endl;
+
+ if(cosmo_ == NULL || redshref_<0.) {
+   cout<<"GeneFluct3D::LosComRedshift_Error: set Observator and Cosmology first"<<endl;
+   throw ParmError("");
+ }
+
+ // On calcule les coordonnees de l'observateur
+ // Il est sur un axe centre sur le milieu de la face Oxy
+ double loscom_ref_ = cosmo_->Dloscom(redshref_);
+ xobs_[0] = Nx_/2.*Dx_;
+ xobs_[1] = Ny_/2.*Dy_;
+ xobs_[2] = kredshref_*Dz_ - loscom_ref_;
+
+ // L'observateur est-il dans le cube?
+ bool obs_in_cube = false;
+ if(xobs_[2]>=0. && xobs_[2]<=Nz_*Dz_) obs_in_cube = true;
+
+ // Find MINIMUM los com distance to the observer:
+ // c'est le centre de la face a k=0
+ // (ou zero si l'observateur est dans le cube)
+ loscom_min_ = 0.;
+ if(!obs_in_cube) loscom_min_ = -xobs_[2];
+
+ // Find MAXIMUM los com distance to the observer:
+ // ou que soit positionne l'observateur, la distance
+ // maximal est sur un des coins du cube
+ loscom_max_ = 0.;
+ for(long i=0;i<=1;i++) {
+   double dx2 = xobs_[0] - i*Nx_*Dx_; dx2 *= dx2;
+   for(long j=0;j<=1;j++) {
+     double dy2 = xobs_[1] - j*Ny_*Dy_; dy2 *= dy2;
+     for(long k=0;k<=1;k++) {
+       double dz2 = xobs_[2] - k*Nz_*Dz_; dz2 *= dz2;
+       dz2 = sqrt(dx2+dy2+dz2);
+       if(dz2>loscom_max_) loscom_max_ = dz2;
+     }
+   }
+ }
+ if(lp_>0) {
+   cout<<"...zref="<<redshref_<<" kzref="<<kredshref_<<" losref="<<loscom_ref_<<" Mpc\n"
+       <<"   xobs="<<xobs_[0]<<" , "<<xobs_[1]<<" , "<<xobs_[2]<<" Mpc "
+       <<" in_cube="<<obs_in_cube
+       <<" loscom_min="<<loscom_min_<<" loscom_max="<<loscom_max_<<" Mpc "<<endl;
+ }
+
+ // Fill the corresponding vectors
+ for(double z=0.; ; z+=zinc) {
+   double dlc = cosmo_->Dloscom(z);
+   if(dlc<loscom_min_) {zred_.resize(0); loscom_.resize(0);}
+   zred_.push_back(z);
+   loscom_.push_back(dlc);
+   z += zinc;
+   if(dlc>loscom_max_) break; // on break apres avoir stoque un dlc>dlcmax
+ }
+
+ long n = zred_.size();
+ if(lp_>0) {
+   cout<<"...n="<<n;
+   if(n>0) cout<<" z="<<zred_[0]<<" -> d="<<loscom_[0];
+   if(n>1) cout<<" , z="<<zred_[n-1]<<" -> d="<<loscom_[n-1];
+   cout<<endl;
+ }
+
+ return n;
+}
 
 //-------------------------------------------------------
@@ -323 +419 @@
 
  // --- On remplit avec une realisation
- if(lp_>0) cout<<"...before Fourier realization filling ---"<<endl;
+ if(lp_>0) cout<<"...before Fourier realization filling"<<endl;
  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;
@@ -529 +625 @@
 
 //-------------------------------------------------------------------
+void GeneFluct3D::ApplyGrowthFactor(long npoints)
+// Apply Growth to real space
+// Using the correspondance between redshift and los comoving distance
+// describe in vector "zred_" "loscom_"
+{
+ if(npoints<3) npoints = zred_.size();
+ if(lp_>0) cout<<"--- ApplyGrowthFactor ---  npoints="<<npoints<<endl;
+ check_array_alloc();
+
+ if(growth_ == NULL) {
+   cout<<"GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first"<<endl;
+   throw ParmError("GeneFluct3D::ApplyGrowthFactor_Error: set GrowthFactor first");
+ }
+
+ long n = zred_.size();
+ InverseFunc invfun(zred_,loscom_);
+ vector<double> invlos;
+ invfun.ComputeParab(npoints,invlos);
+
+ InterpFunc interpinv(invfun.YMin(),invfun.YMax(),invlos);
+ unsigned short ok;
+
+ //CHECK: Histo hgr(0.9*zred_[0],1.1*zred_[n-1],1000);
+ for(long i=0;i<Nx_;i++) {
+   double dx2 = xobs_[0] - i*Dx_; dx2 *= dx2;
+   for(long j=0;j<Ny_;j++) {
+     double dy2 = xobs_[1] - j*Dy_; dy2 *= dy2;
+     for(long l=0;l<Nz_;l++) {
+       double dz2 = xobs_[2] - l*Dz_; dz2 *= dz2;
+       dz2 = sqrt(dx2+dy2+dz2);
+       double z = interpinv(dz2);
+       //CHECK: hgr.Add(z);
+       double dzgr = (*growth_)(z);   // interpolation par morceau
+       //double dzgr = growth_->Linear(z,ok);  // interpolation lineaire
+       //double dzgr = growth_->Parab(z,ok);  // interpolation parabolique
+       int_8 ip = IndexR(i,j,l);
+       data_[ip] *= dzgr;
+     }
+   }
+ }
+
+ //CHECK: {POutPersist pos("applygrowth.ppf"); string tag="hgr"; pos.PutObject(hgr,tag);}
+
+}
+
+//-------------------------------------------------------------------
 void GeneFluct3D::ComputeReal(void)
 // Calcule une realisation dans l'espace reel

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -14 +14 @@
 #include <vector>
 
+#include "cosmocalc.h"
 #include "pkspectrum.h"
 
@@ -27 +28 @@
   void SetNThread(unsigned short nthread=0) {nthread_ = nthread;}
   void SetSize(long nx,long ny,long nz,double dx,double dy,double dz);  // Mpc
-  void SetObservator(double redshref=0.,double kredshref=-1.);
+  void SetObservator(double redshref=0.,double kredshref=0.);
+  void SetCosmology(CosmoCalc& cosmo);
+  void SetGrowthFactor(GrowthFactor& growth);
+  long LosComRedshift(double zinc=0.001);
 
   TArray< complex<r_8> >& GetComplexArray(void) {return T_;}
@@ -67 +71 @@
   int  ComputeSpectrum(HistoErr& herr);
   int  ComputeSpectrum2D(Histo2DErr& herr);
+  void ApplyGrowthFactor(long npoints=-1);
 
   int_8 VarianceFrReal(double R,double& var);
@@ -127 +132 @@
   // l'observateur
   double redshref_,kredshref_;
+  CosmoCalc *cosmo_;
+  GrowthFactor *growth_;
+  double xobs_[3];
+  double loscom_ref_, loscom_min_, loscom_max_;
+  vector<double> zred_, loscom_;
 };
 

trunk/Cosmo/SimLSS/geneutils.cc

@@ -19 +19 @@
 // Classe d'interpolation lineaire:
 // Le vecteur y a n elements y_i tels que y_i = f(x_i)
-// Les x_i sont regulierement espaces
-//   et x_0=xmin et x_{n-1}=xmax
+//   ou les x_i sont regulierement espaces
+//   et x_0=xmin et x_{n-1}=xmax   (xmax inclus!)
 InterpFunc::InterpFunc(double xmin,double xmax,vector<double>& y)
   : _xmin(xmin), _xmax(xmax), _y(y)
@@ -29 +29 @@
   }
   _nm1   = _y.size()-1;
-  _xlarg = _xmax-_xmin;
-  _dx    = _xlarg/(double)_nm1;
+  _dx    = (_xmax-_xmin)/(double)_nm1;
 }
 
@@ -37 +36 @@
   ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2;
   x -= _xmin;
-  long i = long(x/_xlarg*_nm1);  // On prend le "i" juste en dessous
+  long i = long(x/_dx);  // On prend le "i" juste en dessous
   if(i<0) i=0; else if(i>=_nm1) i=_nm1-1;
   return _y[i] + (_y[i+1]-_y[i])/_dx*(x-i*_dx);
@@ -46 +45 @@
   ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2;
   x -= _xmin;
-  long i = long(x/_xlarg*_nm1+0.5);  // On prend le "i" le + proche
+  long i = long(x/_dx+0.5);  // On prend le "i" le + proche
   if(i<1) i=1; else if(i>=_nm1-1) i=_nm1-2;
   double a = (_y[i+1]-2.*_y[i]+_y[i-1])/(2.*_dx*_dx);
   double b = (_y[i+1]-_y[i-1])/(2.*_dx);
   return _y[i] + (x-i*_dx)*(a*(x-i*_dx)+b);
+}
+
+//-------------------------------------------------------------------
+// Classe d'inversion d'une fonction STRICTEMENT MONOTONE CROISSANTE
+// - Le vecteur y a "Nin" elements y_i tels que "y_i = f(x_i)"
+// - On a x(i) < x(i+1) et y(i) < y(i+1)
+// - La classe renvoie ymin=y(0) , ymax=y(Nin -1)
+//   et le vecteur x = f^-1(y) de "Nout" elements
+//   Les y_i sont regulierement espaces et ymin et ymax
+//   La re-interpolation inverse est faite par lineaire
+InverseFunc::InverseFunc(vector<double>& x,vector<double>& y)
+  : _ymin(0.) , _ymax(0.) , _x(x) , _y(y)
+{
+  int_4 ns = _x.size();
+  if(ns<3 || _y.size()<=0 || ns!=_y.size())
+    throw ParmError("InverseFunc::InverseFunc_Error: bad array size");
+
+  // Check "x" strictement monotone croissant
+  for(int_4 i=0;i<ns-1;i++)
+    if(_x[i+1]<=_x[i]) {
+      cout<<"InverseFunc::InverseFunc_Error: _x array not stricly growing"<<endl;
+      throw ParmError("InverseFunc::InverseFunc_Error: _x array not stricly growing");
+    }
+
+  // Check "y" monotone croissant
+  for(int_4 i=0;i<ns-1;i++)
+    if(_y[i+1]<_y[i]) {
+      cout<<"InverseFunc::InverseFunc_Error: _y array not growing"<<endl;
+      throw ParmError("InverseFunc::InverseFunc_Error: _y array not growing");
+    }
+
+  // define limits
+  _ymin = _y[0];
+  _ymax = _y[ns-1];
+
+}
+
+InverseFunc::~InverseFunc(void)
+{
+}
+
+int InverseFunc::ComputeLinear(long n,vector<double>& xfcty)
+{
+  if(n<3) return -1;
+
+  xfcty.resize(n);
+
+  long i1,i2;
+  double x;
+  for(int_4 i=0;i<n;i++) {
+    double y = _ymin + i*(_ymax-_ymin)/(n-1.);
+    find_in_y(y,i1,i2);
+    double dy = _y[i2]-_y[i1];
+    if(dy==0.) {
+      x = (_x[i2]+_x[i1])/2.; // la fct a inverser est plate!
+    } else {
+      x = _x[i1] + (_x[i2]-_x[i1])/dy * (y-_y[i1]);
+    }
+    xfcty[i] = x;
+  }
+
+  return 0;
+}
+
+int InverseFunc::ComputeParab(long n,vector<double>& xfcty)
+{
+  if(n<3) return -1;
+
+  xfcty.resize(n);
+
+  long i1,i2,i3;
+  double x;
+  for(int_4 i=0;i<n;i++) {
+    double y = _ymin + i*(_ymax-_ymin)/(n-1.);
+    find_in_y(y,i1,i2);
+    // On cherche le 3ieme point selon la position de y / au 2 premiers
+    double my = (_y[i1]+_y[i2])/2.;
+    if(y<my) {i3=i2; i2=i1; i1--;} else {i3=i2+1;}
+    // Protection
+    if(i1<0) {i1++; i2++; i3++;}
+    if(i3==_y.size()) {i1--; i2--; i3--;}
+    // Interpolation parabolique
+    double dy = _y[i3]-_y[i1];
+    if(dy==0.) {
+      x = (_x[i3]+_x[i1])/2.; // la fct a inverser est plate!
+    } else {
+      double X1=_x[i1]-_x[i2], X3=_x[i3]-_x[i2];
+      double Y1=_y[i1]-_y[i2], Y3=_y[i3]-_y[i2];
+      double den = Y1*Y3*dy;
+      double a = (X3*Y1-X1*Y3)/den;
+      double b = (X1*Y3*Y3-X3*Y1*Y1)/den;
+      y -= _y[i2];
+      x = (a*y+b)*y + _x[i2];
+    }
+    xfcty[i] = x;
+  }
+
+  return 0;
 }
 

trunk/Cosmo/SimLSS/geneutils.h

@@ -7 +7 @@
 #include "histos.h"
 #include "tvector.h"
+#include "cspline.h"
 
 #include <vector>
@@ -12 +13 @@
 namespace SOPHYA {
 
+//----------------------------------------------------
 class InterpFunc {
 public:
@@ -19 +21 @@
   double XMin(void) {return _xmin;}
   double XMax(void) {return _xmax;}
+  inline double X(long i) {return _xmin + i*_dx;}
 
   //! Retourne l'element le plus proche de f donnant y=f(x)
@@ -24 +27 @@
          {
          x -= _xmin;
-         long i = long(x/_xlarg*_nm1+0.5);  // On prend le "i" le plus proche
+         long i = long(x/_dx+0.5);  // On prend le "i" le plus proche
          if(i<0) i=0; else if(i>=_nm1) i=_nm1-1;
          return _y[i];
@@ -43 +46 @@
 
 protected:
-  double _xmin,_xmax,_xlarg,_dx;
+  double _xmin,_xmax,_dx;
+  long _nm1;  // n-1
   vector<double>& _y;
-  long _nm1;  // n-1
+};
+
+//----------------------------------------------------
+class InverseFunc {
+public:
+  InverseFunc(vector<double>& x,vector<double>& y);
+  virtual ~InverseFunc(void);
+  int ComputeLinear(long n,vector<double>& xfcty);
+  int ComputeParab(long n,vector<double>& xfcty);
+  double YMin(void) {return _ymin;}
+  double YMax(void) {return _ymax;}
+protected:
+  inline void find_in_y(double x,long& klo,long& khi)
+    {
+      long k;
+      klo=0, khi=_y.size()-1;
+      while (khi-klo > 1) {
+        k = (khi+klo) >> 1;
+        if (_y[k] > x) khi=k; else klo=k;
+      }
+    }
+
+  vector<double>& _x;
+  vector<double>& _y;
+  double _ymin,_ymax,_dy;
 };
 
 }
 
+//----------------------------------------------------
 int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex=false);
 int FuncToVec(GenericFunc& func,TVector<r_8>& h,double xmin,double xmax,bool logaxex=false);
 
+//----------------------------------------------------
 double AngSol(double dtheta,double dphi,double theta0=M_PI/2.);
 double AngSol(double dtheta);
 
+//----------------------------------------------------
 unsigned long PoissRandLimit(double mu,double mumax=10.);