Changeset 3320 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Sep 5, 2007, 9:58:05 AM (18 years ago)

Author:

cmv

Message:

intro SchechterMassDist pour accelerer la simulations cmv 05/09/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvschdist.cc (modified) (9 diffs)
• cmvtstsch.cc (modified) (2 diffs)
• genefluct3d.cc (modified) (2 diffs)
• genefluct3d.h (modified) (2 diffs)
• geneutils.cc (modified) (1 diff)
• schechter.cc (modified) (6 diffs)
• schechter.h (modified) (3 diffs)

trunk/Cosmo/SimLSS/cmvschdist.cc

@@ -6 +6 @@
 // ex: on a un pixel qui a N galaxies, quelle est la distribution
 //     statistique de masse dans un tel pixel ?
-// > cmvschdist -a -v -m 1e+7,1e+14 -n 140 -r 200 -N 100000
+// > cmvschdist -a -v -m 1e+8,1e+13 -n -100 -r 200 -N 100000
+// compare with:
+// > cmvschdist -a -v -m 1e+8,1e+13 -n -100 -r 200 -N 100000 -0
+// Fabriquer un fichier
+// > cmvschdist -a -m 1e+8,1e+13 -n -100 -r 2000 -N 5000000 -W schdist_2000.ppf
 #include <iostream>
 #include <stdlib.h>
@@ -28 +32 @@
 void usage(void)
 {
- cout<<"cmvschdist -m xmin,xmax -a -n nbin -r ngmax,ngmin -N nalea"<<endl
-     <<"  xmin,xmax : limites en masse"<<endl
-     <<"  nbin : nombre de bins"<<endl
-     <<"  ngmax,ngmin : distribution pour tirage de ngmin a ngmax galaxies"<<endl
-     <<"  nalea : nombre de tirages aleatoires"<<endl
+ cout<<"cmvschdist -m xmin,xmax -a -0 -n nbin -r ngmax,ngmin -N nalea "
+     <<"-R readfile.ppf -W writefile.ppf"<<endl
+     <<"  -m xmin,xmax : limites en masse"<<endl
+     <<"  -n nbin : nombre de bins (si <0 alors nb par decade)"<<endl
+     <<"  -r ngmax,ngmin : distribution pour tirage de ngmin a ngmax galaxies"<<endl
+     <<"  -N nalea : nombre de tirages aleatoires"<<endl
      <<"  -v : back verification of random trials"<<endl
-     <<"  -a : init auto de l\'aleatoire"<<endl;
+     <<"  -0 : only use ngal=1 histogram"<<endl
+     <<"  -a : init auto de l\'aleatoire"<<endl
+     <<"  -R readfile.ppf : read SchechterMassDist from ppf file"<<endl;
 }
 
@@ -46 +53 @@
 
  double xmin=1e8, xmax=1e13;
- int npt = 100;
+ int npt = -100;
  int ngmax = 200, ngmin = 1;
+ bool useonly1 = false;
  unsigned long long nalea = 10000;
  bool verif = false;
+ bool readfrppf = false;
+ string namefrppf = "";
+ string nametoppf = "";
 
  char c;
- while((c = getopt(narg,arg,"havm:N:n:r:")) != -1) {
+ while((c = getopt(narg,arg,"h0avm:N:n:r:R:W:")) != -1) {
    switch (c) {
    case 'm' :
@@ -68 +79 @@
    case 'n' :
      sscanf(optarg,"%d",&npt);
-     if(npt<=0) npt = 100;
+     if(npt<=0) npt = -100;
      break;
    case 'N' :
@@ -77 +88 @@
      verif = true;
      break;
+   case '0' :
+     useonly1 = true;
+     break;
+   case 'R' :
+     readfrppf = true;
+     namefrppf = optarg;
+     break;
+   case 'W' :
+     nametoppf = optarg;
+     break;
    case 'a' :
      Auto_Ini_Ranf(5);
@@ -86 +107 @@
  }
 
- double lnx1=log10(xmin), lnx2=log10(xmax), dlnx = (lnx2-lnx1)/npt;
- cout<<"xmin="<<xmin<<" ("<<lnx1<<") xmax="<<xmax<<" ("<<lnx2<<"), dlnx="<<dlnx<<endl;
+ double lnx1=log10(xmin), lnx2=log10(xmax);
+ cout<<"xmin="<<xmin<<" ("<<lnx1<<") xmax="<<xmax<<" ("<<lnx2<<")"<<endl;
  cout<<"npt="<<npt<<",  nalea="<<nalea<<endl;
- cout<<"ngmin="<<ngmin<<" ngmax="<<ngmax<<endl;
+ cout<<"useonly1="<<useonly1<<" ngmin="<<ngmin<<" ngmax="<<ngmax<<endl;
  cout<<"verif="<<verif<<endl;
+ if(readfrppf) cout<<"SchechterMassDist will be read from file "<<namefrppf<<endl;
+ if(nametoppf.size()>0) cout<<"SchechterMassDist will be written to file "<<nametoppf<<endl;
 
 
@@ -96 +119 @@
  cout<<"> Schechter m*dn/dm nstar="<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
  Schechter sch(nstar,mstar,alpha);
- sch.SetOutValue(1);  // on veut m*dN/dm
- Histo hmdndm(lnx1,lnx2,npt); hmdndm.ReCenterBin();
- FuncToHisto(sch,hmdndm,true);
- FunRan tirhmdndm = FunRan(hmdndm,true);
-
- //------- Construct histo
- int nbhist = ngmax-ngmin+1;
- cout<<"> Creating "<<nbhist<<" histos"<<endl;
- vector<Histo> vhisto;  vector<string> vhname;
- for(int i=ngmin;i<=ngmax;i++) {
-   Histo h(hmdndm); h.Zero();
-   vhisto.push_back(h);
-   char str[32]; sprintf(str,"h%d",i);
-   vhname.push_back(string(str));
- }
+
+ //------- Construct Mass Distribution
+ cout<<"> Creating  Mass Distribution"<<endl;
+ SchechterMassDist schdmass(sch,xmin,xmax,npt);
+ if(readfrppf) {
+   cout<<"> Mass Distribution read from "<<namefrppf<<endl;
+   schdmass.ReadPPF(namefrppf);
+ }
+ Histo hmdndm = schdmass.GetHmDnDm();
+ FunRan tirhmdndm = schdmass.GetTmDnDm();
 
  //------- Random
  InitTim();
- cout<<"> Random: "<<nalea<<" trials"<<endl;
- int lpmod = nalea/25; if(lpmod<=0) lpmod=1;
- double bigsum = 0.; long nbigsum=0;
- for(unsigned long long ia=0;ia<nalea;ia++) {
-   if(ia%lpmod==0) cout<<"... "<<ia<<endl;
-   double sum = 0.;
-   for(int i=1;i<=ngmax;i++) {
-     //double l10m = tirhmdndm.Random();
-     double l10m = tirhmdndm.RandomInterp();
-     double m = pow(10.,l10m);
-     sum += m;
-     bigsum += m; nbigsum++;
-     int ipo = i-ngmin;
-     if(ipo<0) continue;
-     double v = log10(sum/(double)i);
-     vhisto[ipo].Add(v);
-   }
-   if(ia%lpmod==0) PrtTim(" ");
- }
- PrtTim("End Random loop");
- if(nbigsum>0) {
-   bigsum /= (double)nbigsum;
-   cout<<"Mean mass : "<<bigsum<<" ("<<log10(fabs(bigsum))<<")"<<endl;
- }
-
- //------- Generation des classes de tirage aleatoire et des histos de verif
- vector<FunRan> vtir;  vector<string> vtname;
- vector<Histo> vthisto;  vector<string> vthname;
+ if(!useonly1 && !readfrppf) {
+   cout<<"> Creating "<<ngmax-ngmin+1<<" histos and filling with  "<<nalea<<" trials"<<endl;
+   schdmass.SetNgalLim(ngmax,ngmin,nalea);
+   PrtTim("End of filling histos for trials");
+ }
+ schdmass.Print();
+
+ //------- Generation des histos de verif
+ vector<Histo> vthisto; vector<string> vthname;
  if(verif) {
-   cout<<"> Creating "<<nbhist<<" FunRan and Histos for back-verif"<<endl;
+   cout<<"> Creating "<<ngmax-ngmin+1<<" Histos for back-verif"<<endl;
    for(int i=ngmin;i<=ngmax;i++) {
-     FunRan t(vhisto[i-ngmin],true);
-     vtir.push_back(t);
-     char str[32]; sprintf(str,"t%d",i);
-     vtname.push_back(string(str));
-     Histo h(vhisto[i-ngmin]); h.Zero();
+     Histo h(hmdndm); h.Zero();
      vthisto.push_back(h);
-     sprintf(str,"th%d",i);
+     char str[32]; sprintf(str,"th%d",i);
      vthname.push_back(string(str));
    }
+   cout<<"   "<<vthisto.size()<<" Histos created"<<endl;
+   vthisto[20].Show();
 
    cout<<"> Checking tirage"<<endl;
+   int lpmod = nalea/20; if(lpmod<=0) lpmod=1;
    for(unsigned long long ia=0;ia<nalea;ia++) {
-     for(unsigned int i=0;i<vtir.size();i++) {
-       //double v = vtir[i].Random();
-       double v = vtir[i].RandomInterp();
-       vthisto[i].Add(v);
+     if(ia%lpmod==0) cout<<"...filling tirage "<<ia<<endl;
+     for(unsigned int i=0;i<vthisto.size();i++) {
+       int ng = ngmin+i;
+       double v = schdmass.TirMass(ng)/(double)ng;
+       vthisto[i].Add(log10(v));
      }
    }
+   schdmass.PrintStatus();
    PrtTim("End Random loop back-tirage");
  }
 
- //------- Ecriture ppf
- cout<<"Ecriture"<<endl;
+ //------- Ecritures ppf
+ if(nametoppf.size()>0) {
+   cout<<"Ecriture de l\'objet SchechterMassDist dans "<<nametoppf<<endl;
+   schdmass.WritePPF(nametoppf);
+ }
+
+ cout<<"Ecriture pour verification"<<endl;
  string tag = "cmvschdist.ppf";
  POutPersist pos(tag);
+ {
+ cout<<"  writing hmdndm tirhmdndm"<<endl;
  tag = "hmdndm"; pos.PutObject(hmdndm,tag);
- Histo hdum2(tirhmdndm);
- tag = "tirhmdndm"; pos.PutObject(hdum2,tag);
- if(vhisto.size()>0)
-   for(unsigned int i=0;i<vhisto.size();i++)
-     if(vhisto[i].NEntries()>0) pos.PutObject(vhisto[i],vhname[i]);
- if(vthisto.size()>0)
-   for(unsigned int i=0;i<vthisto.size();i++)
+ Histo hdum(tirhmdndm);
+ tag = "tirhmdndm"; pos.PutObject(hdum,tag);
+ }
+ if(schdmass.GetNgalLim()>0) {
+   cout<<"  writing h t"<<endl;
+   for(int i=0;i<schdmass.GetNgalLim();i++) {
+     int ng = schdmass.NGalFrIndex(i);
+     if(ng<=0) continue;
+     char str[32];
+     sprintf(str,"h%d",ng);
+     Histo hdum = schdmass.GetHisto(i);
+     tag = str; pos.PutObject(hdum,tag);
+     sprintf(str,"t%d",ng);
+     Histo hdum2(schdmass.GetFunRan(i));
+     tag = str; pos.PutObject(hdum2,tag);
+   }
+ }
+ if(vthisto.size()>0) {
+  cout<<"  writing th"<<endl;
+  for(unsigned int i=0;i<vthisto.size();i++)
      if(vthisto[i].NEntries()>0) pos.PutObject(vthisto[i],vthname[i]);
+ }
+
+ PrtTim("End of Job");
 
  return 0;
@@ -191 +214 @@
 
 set h h$n
+set t t$n
 set th th$n
+
+disp tirhmdndm
+disp $t "same red"
 
 disp $h
@@ -201 +228 @@
 n/plot $h.log10(val)%x val>0. ! "connectpoints"
 n/plot $th.log10(val)%x val>0. ! "nsta connectpoints same red"
+
+n/plot $h.val%pow(10.,x) ! ! "connectpoints"
+n/plot $th.val%pow(10.,x) ! ! "nsta connectpoints same red"
 
 c++exec \

trunk/Cosmo/SimLSS/cmvtstsch.cc

@@ -93 +93 @@
  cout<<"......"<<endl;
  for(double lnx=lnx1-3.;lnx<=lnx2-1.;lnx+=1.) {
-   double num = IntegrateFuncLog(sch,lnx,lnx2,perc,dlxinc,dlxmax,glorder);
-   cout<<"["<<lnx<<","<<lnx2<<"] integrated number = "<<num<<" Mpc^-3"<<endl;
+   double num = sch.Integrate(pow(10.,lnx),pow(10.,lnx2),npt);
+   double numc = IntegrateFuncLog(sch,lnx,lnx2,perc,dlxinc,dlxmax,glorder);
+   cout<<"["<<lnx<<","<<lnx2<<"] integrated number = "<<num<<" Mpc^-3 , check: "<<numc<<endl;
  }
  cout<<"......"<<endl;
  for(double lnx=lnx1+1.;lnx<=lnx2+3.;lnx+=1.) {
-   double num = IntegrateFuncLog(sch,lnx1,lnx,perc,dlxinc,dlxmax,glorder);
-   cout<<"["<<lnx1<<","<<lnx<<"] integrated number = "<<num<<" Mpc^-3"<<endl;
+   double num = sch.Integrate(pow(10.,lnx1),pow(10.,lnx),npt);
+   double numc = IntegrateFuncLog(sch,lnx1,lnx,perc,dlxinc,dlxmax,glorder);
+   cout<<"["<<lnx1<<","<<lnx<<"] integrated number = "<<num<<" Mpc^-3 , check: "<<numc<<endl;
  }
  cout<<"......"<<endl;
  for(double lnx=lnx1-3.;lnx<=lnx2+3.;lnx+=1.) {
-   double num = IntegrateFuncLog(sch,lnx,lnx+1.,perc,dlxinc,dlxmax,glorder);
-   cout<<"["<<lnx<<","<<lnx+1.<<"] integrated number = "<<num<<" Mpc^-3"<<endl;
+   double num = sch.Integrate(pow(10.,lnx),pow(10.,lnx+1.),npt);
+   double numc = IntegrateFuncLog(sch,lnx,lnx+1.,perc,dlxinc,dlxmax,glorder);
+   cout<<"["<<lnx<<","<<lnx+1.<<"] integrated number = "<<num<<" Mpc^-3 , check: "<<numc<<endl;
  }
  Histo hdndm(lnx1,lnx2,npt); hdndm.ReCenterBin();
@@ -115 +118 @@
  cout<<"......"<<endl;
  for(double lnx=lnx1-3.;lnx<=lnx2-1.;lnx+=1.) {
-   double sum = IntegrateFuncLog(sch,lnx,lnx2,perc,dlxinc,dlxmax,glorder);
-   cout<<"["<<lnx<<","<<lnx2<<"] integrated mass = "<<sum<<" Msol.Mpc^-3"<<endl;
+   double sum = sch.Integrate(pow(10.,lnx),pow(10.,lnx2),npt);
+   double sumc = IntegrateFuncLog(sch,lnx,lnx2,perc,dlxinc,dlxmax,glorder);
+   cout<<"["<<lnx<<","<<lnx2<<"] integrated mass = "<<sum<<" Msol.Mpc^-3 , check: "<<sumc<<endl;
  }
  cout<<"......"<<endl;
  for(double lnx=lnx1+1.;lnx<=lnx2+3.;lnx+=1.) {
-   double sum = IntegrateFuncLog(sch,lnx1,lnx,perc,dlxinc,dlxmax,glorder);
-   cout<<"["<<lnx1<<","<<lnx<<"] integrated mass = "<<sum<<" Msol.Mpc^-3"<<endl;
+   double sum = sch.Integrate(pow(10.,lnx1),pow(10.,lnx),npt);
+   double sumc = IntegrateFuncLog(sch,lnx1,lnx,perc,dlxinc,dlxmax,glorder);
+   cout<<"["<<lnx1<<","<<lnx<<"] integrated mass = "<<sum<<" Msol.Mpc^-3 , check: "<<sumc<<endl;
  }
  cout<<"......"<<endl;
  for(double lnx=lnx1-3.;lnx<=lnx2+3.;lnx+=1.) {
-   double sum = IntegrateFuncLog(sch,lnx,lnx+1.,perc,dlxinc,dlxmax,glorder);
-   cout<<"["<<lnx<<","<<lnx+1.<<"] integrated mass = "<<sum<<" Msol.Mpc^-3"<<endl;
+   double sum = sch.Integrate(pow(10.,lnx),pow(10.,lnx+1.),npt);
+   double sumc = IntegrateFuncLog(sch,lnx,lnx+1.,perc,dlxinc,dlxmax,glorder);
+   cout<<"["<<lnx<<","<<lnx+1.<<"] integrated mass = "<<sum<<" Msol.Mpc^-3 , check: "<<sumc<<endl;
  }
  Histo hmdndm(lnx1,lnx2,npt); hmdndm.ReCenterBin();

trunk/Cosmo/SimLSS/genefluct3d.cc

@@ -940 +940 @@
 }
 
+int_8 GeneFluct3D::MinMax(double& xmin,double& xmax,double vmin,double vmax)
+// Calcul des valeurs xmin et xmax dans le cube reel avec valeurs ]vmin,vmax[ extremites exclues
+{
+ bool tstval = (vmax>vmin)? true: false;
+ if(lp_>0) {
+   cout<<"--- MinMax";
+   if(tstval) cout<<"  range=]"<<vmin<<","<<vmax<<"[";
+   cout<<endl;
+ }
+ check_array_alloc();
+
+ int_8 n = 0;
+ xmin = xmax = data_[0];
+
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
+   int_8 ip = IndexR(i,j,l);
+   double x = data_[ip];
+   if(tstval && (x<=vmin || x>=vmax)) continue;
+   if(x<xmin) xmin = x;
+   if(x>xmax) xmax = x;
+   n++;
+ }
+
+ if(lp_>0) cout<<"  n="<<n<<" min="<<xmin<<" max="<<xmax<<endl;
+
+ return n;
+}
+
 int_8 GeneFluct3D::MeanSigma2(double& rm,double& rs2,double vmin,double vmax
                              ,bool useout,double vout)
@@ -1114 +1142 @@
        data_[ip] += m;
      }
+     sum += data_[ip];
+   }
+ }
+ if(lp_>0) cout<<sum<<" MSol HI mass put into survey"<<endl;
+
+ return sum;
+}
+
+double GeneFluct3D::TurnNGal2MassQuick(SchechterMassDist& schmdist)
+// idem TurnNGal2Mass mais beaucoup plus rapide
+{
+ if(lp_>0) cout<<"--- TurnNGal2MassQuick ---"<<endl;
+ check_array_alloc();
+
+ double sum = 0.;
+ for(long i=0;i<Nx_;i++) for(long j=0;j<Ny_;j++) for(long l=0;l<Nz_;l++) {
+   int_8 ip = IndexR(i,j,l);
+   double v = data_[ip];
+   if(v>0.) {
+     long ngal = long(v+0.1);
+     data_[ip] = schmdist.TirMass(ngal);
      sum += data_[ip];
    }

trunk/Cosmo/SimLSS/genefluct3d.h

@@ -93 +93 @@
   int_8 MeanSigma2(double& rm,double& rs2,double vmin=1.,double vmax=-1.
                   ,bool useout=false,double vout=0.);
+  int_8 MinMax(double& xmin,double& xmax,double vmin=1.,double vmax=-1.);
   int_8 NumberOfBad(double vmin=-1.e+150,double vmax=1.e+150);
   int_8 SetToVal(double vmin, double vmax,double val0=0.);
@@ -101 +102 @@
   double ApplyPoisson(void);
   double TurnNGal2Mass(FunRan& massdist,bool axeslog=false);
+  double TurnNGal2MassQuick(SchechterMassDist& schmdist);
   double TurnMass2Flux(void);
   void AddAGN(double lfjy,double lsigma,double powlaw=0.);

trunk/Cosmo/SimLSS/geneutils.cc

@@ -366 +366 @@
 ////////////////////////////////////////////////////////////////////////////////////
 double IntegrateFuncLog(GenericFunc& func,double lxmin,double lxmax
-         ,double perc,double dlxinc,double dlxmax,unsigned short  glorder)
+         ,double perc,double dlxinc,double dlxmax,unsigned short glorder)
 // --- Integration adaptative ---
 // Idem IntegrateFunc but integrate on logarithmic (base 10) intervals:

trunk/Cosmo/SimLSS/schechter.cc

@@ -8 +8 @@
 
 #include "pexceptions.h"
+#include "histos.h"
+#include "perandom.h"
+#include "tvector.h"
+#include "fioarr.h"
 
 #include "constcosmo.h"
+#include "geneutils.h"
 #include "schechter.h"
 
@@ -27 +32 @@
 Schechter::Schechter(Schechter& f)
   : nstar_(f.nstar_) , mstar_(f.mstar_) , alpha_(f.alpha_) , outvalue_(f.outvalue_)
+{
+}
+
+Schechter::Schechter(void)
+  : nstar_(0.) , mstar_(0.) , alpha_(0.) , outvalue_(0)
 {
 }
@@ -45 +55 @@
 }
 
+unsigned short Schechter::GetOutValue(void)
+{
+  return outvalue_;
+}
+
+void Schechter::SetParam(double nstar,double mstar,double alpha)
+{
+  nstar_ = nstar; mstar_ = mstar; alpha_ = alpha;
+}
+
+void Schechter::GetParam(double& nstar,double& mstar,double& alpha)
+{
+  nstar = nstar_; mstar = mstar_; alpha = alpha_;
+}
+
 double Schechter::operator() (double m)
 // Return : "dn/dm = f(m)" or "m*dn/dm = f(m)"
@@ -52 +77 @@
   if(outvalue_) dndm *= m;  // on veut m*dn/dm
   return dndm;
+}
+
+
+double Schechter::Integrate(double massmin,double massmax,int npt)
+// Integrate from massmin to massmax with at least npt points log10 spaced
+{
+ if(massmin<=0. || massmax<=0. || massmax<=massmin) return 0.;
+ if(npt<1) npt = 100;
+ double lmassmin = log10(massmin), lmassmax = log10(massmax);
+ double perc=0.01, dlxinc=(lmassmax-lmassmin)/npt, dlxmax=10.*dlxinc; unsigned short glorder=4;
+ double sum = IntegrateFuncLog(*this,lmassmin,lmassmax,perc,dlxinc,dlxmax,glorder);
+ return sum;
 }
 
@@ -63 +100 @@
   cout<<endl;
 }
+
+///////////////////////////////////////////////////////////
+//******* Les facilites pour tirer sur Schechter ********//
+///////////////////////////////////////////////////////////
+
+SchechterMassDist::SchechterMassDist(Schechter sch,double massmin,double massmax,int nbinmass)
+// On veut une fonction de tirage aleatoire de la masse de N galaxies sur une Schechter
+// Pour optimiser on cree (eventuellement) des histos de tirage
+// ATTENTION: on donne les limites en masses pour les histos mais leurs abscisses
+//   sont en log10(masse): histo [log10(massmin),log10(massmax)] avec nbinmass bins
+// Si nbinmass<0 alors c'est le nombre de points par decade
+  : sch_(sch) , sch_outvalue_(sch.GetOutValue())
+  , massmin_(massmin) , massmax_(massmax) , nbinmass_(nbinmass)
+  , ngalmin_(0) , ngalmax_(0) , nvalngal_(0)
+  , ntrial_dir(0) , ntrial_tab(0)
+  , hmdndm_(NULL) , tirhmdndm_(NULL)
+{
+  if(massmin_>massmax_  || massmin_<=0. || massmax_<=0.|| nbinmass_==0) {
+    cout<<"SchechterMassDist::SchechterMassDist: error in input values"<<endl;
+    throw ParmError("SchechterMassDist::SchechterMassDist: error in input values");
+  }
+
+  // Creation du tirage aleatoire sur la Schechter
+  sch_outvalue_ = sch.GetOutValue();  // on sauve la configuration de la schechter
+  sch.SetOutValue(1);  // on veut m*dN/dm
+  double lnx1 = log10(massmin_), lnx2 = log10(massmax_); // Binning en log10 de la masse
+  if(nbinmass_<0) nbinmass_ = int((-nbinmass_)*(lnx2-lnx1+1.));
+  hmdndm_ = new Histo(lnx1,lnx2,nbinmass_); hmdndm_->ReCenterBin();
+  FuncToHisto(sch,*hmdndm_,true);  // true -> bin en log10(x)
+  tirhmdndm_ = new FunRan(*hmdndm_,true);  // true -> histo est pdf
+  sch.SetOutValue(sch_outvalue_);  // on remet a la valeur initiale
+}
+
+SchechterMassDist::SchechterMassDist(void)
+  : sch_outvalue_(0)
+  , massmin_(0.) , massmax_(0.) , nbinmass_(0)
+  , ngalmin_(0) , ngalmax_(0) , nvalngal_(0)
+  , ntrial_dir(0) , ntrial_tab(0)
+  , hmdndm_(NULL) , tirhmdndm_(NULL)
+{
+}
+
+SchechterMassDist::~SchechterMassDist(void)
+{
+ Delete();
+}
+
+void SchechterMassDist::Delete(void)
+{
+ if(hmdndm_) delete hmdndm_;
+ if(tirhmdndm_) delete tirhmdndm_;
+ hmass_.resize(0);
+ tmass_.resize(0);
+ ntrial_dir = ntrial_tab = 0;
+}
+
+int SchechterMassDist::GetMassLim(double& massmin,double& massmax)
+{
+ massmin = massmin_;
+ massmax = massmax_;
+ return nbinmass_;
+}
+
+int SchechterMassDist::SetNgalLim(int ngalmax,int ngalmin,unsigned long nalea)
+// Creation des Histos de tirage pour ngalmin a ngalmax galaxies
+{
+ int lp=2;
+ ngalmin_=ngalmax_=nvalngal_=0;
+ if(ngalmin<=0) ngalmin=1;
+ if(ngalmax<ngalmin || ngalmax==1) return 0;
+ ngalmin_ = ngalmin;
+ ngalmax_ = ngalmax;
+ nvalngal_ = ngalmax-ngalmin+1;
+  
+ if(nalea<1) nalea = 100000;
+
+ if(lp>0) cout<<"SchechterMassDist::SetNgalLim: ngal=["
+              <<ngalmin_<<","<<ngalmax_<<"] n="<<nvalngal_
+              <<" filling with "<<nalea<<" trials"<<endl;
+
+ //------- Construct histo
+ double lnx1 = log10(massmin_), lnx2 = log10(massmax_);
+ if(lp>0) cout<<"> Creating "<<nvalngal_<<" histos ["<<lnx1<<","<<lnx2<<"] n="<<nbinmass_<<endl;
+ for(int i=ngalmin_;i<=ngalmax_;i++) {
+   Histo h(*hmdndm_); h.Zero();
+   hmass_.push_back(h);
+ }
+ if(lp>1) cout<<"...number of histos is "<<hmass_.size()<<endl;
+
+ //------- Remplissage random
+ sch_.SetOutValue(1);  // on veut m*dN/dm
+ int lpmod = nalea/20; if(lpmod<=0) lpmod=1;
+ double s1=0., sc1=0.; unsigned long ns1=0;
+ double sax=0., scax=0.; unsigned long nsax=0;
+ for(unsigned long ia=0;ia<nalea;ia++) {
+   if(lp>1 && ia%lpmod==0) cout<<"...tirage "<<ia<<endl;
+   double sum = 0.;
+   for(int i=1;i<=ngalmax_;i++) {
+     //double l10m = tirhmdndm_->Random();
+     double l10m = tirhmdndm_->RandomInterp();
+     double m = pow(10.,l10m);
+     sum += m;
+     s1 += m; sc1 += m*m; ns1++;
+     int ipo = i-ngalmin_;
+     if(ipo<0) continue;
+     // ATTENTION: l'histo de tirage stoque le log10(moyenne=somme_masses/ngal).
+     //            Ca permet d'avoir le meme binning quelque soit ngal
+     double v = log10(sum/(double)i);
+     hmass_[ipo].Add(v);
+     if(i==ngalmax) {sax += sum/(double)i; scax += sum/(double)i*sum/(double)i; nsax++;}
+   }
+ }
+ sch_.SetOutValue(sch_outvalue_);  // on remet a la valeur initiale
+
+ if(ns1>1) {
+   s1 /= ns1; sc1 = sc1/ns1 - s1*s1;
+   cout<<"...Mean mass for ngal=1: "<<s1<<" ("<<log10(fabs(s1))<<")"
+       <<" s="<<sqrt(fabs(sc1))<<" (ntrials "<<ns1<<")"<<endl;
+ }
+ if(nsax>1) {
+   sax /= nsax; scax = scax/nsax - sax*sax;
+   cout<<"...Mean mass for ngal="<<ngalmax_<<": "<<sax<<" ("<<log10(fabs(sax))<<")"
+       <<" s="<<sqrt(fabs(scax))<<" (ntrials "<<nsax<<")"<<endl;
+ }
+
+ //------- Generation des classes de tirage aleatoire et des histos de verif
+ if(lp>0) cout<<"> Creating "<<nvalngal_<<" FunRan"<<endl;
+ for(unsigned int i=0;i<hmass_.size();i++) {
+   FunRan t(hmass_[i],true);
+   tmass_.push_back(t);
+ }
+ if(lp>1) cout<<"...number of funran is "<<tmass_.size()<<endl;
+
+ return nvalngal_;
+}
+
+int SchechterMassDist::GetNgalLim(int& ngalmax,int& ngalmin)
+{
+ ngalmax = ngalmax_;
+ ngalmin = ngalmin_;
+ return nvalngal_;
+}
+
+Histo SchechterMassDist::GetHmDnDm(void) const
+{
+  return *hmdndm_;
+}
+
+FunRan SchechterMassDist::GetTmDnDm(void) const
+{
+  return *tirhmdndm_;
+}
+
+Histo SchechterMassDist::GetHisto(int i) const
+{
+  if(i<0 || i>=nvalngal_) {
+    cout<<"SchechterMassDist::GetHisto: error in input values"<<endl;
+    throw ParmError("SchechterMassDist::GetHisto: error in input values");
+  }
+  return hmass_[i];
+}
+
+FunRan SchechterMassDist::GetFunRan(int i) const
+{
+  if(i<0 || i>=nvalngal_) {
+    cout<<"SchechterMassDist::GetFunRan: error in input values"<<endl;
+    throw ParmError("SchechterMassDist::GetFunRan: error in input values");
+  }
+  return tmass_[i];
+}
+
+double SchechterMassDist::TirMass(int ngal)
+{
+ if(ngal<1) return 0.;
+
+ int ipo = IndexFrNGal(ngal);
+ double masse_des_ngal = 0.;
+ if(ipo<0) {  // Pas d'histos de tirage pour ce nombre de galaxies
+   for(long i=0;i<ngal;i++) {  // On tire donc ngal fois sur la Schechter
+     // double lm = tirhmdndm_->Random();
+     double lm = tirhmdndm_->RandomInterp();
+     masse_des_ngal += pow(10.,lm);  // ATTENTION abscisse en log10(masse)
+   }
+   ntrial_dir++;
+ } else {
+   // ATTENTION l'histo de tirage stoque le log10(moyenne=somme_masses/ngal)
+   //double lmngal = tmass_[ipo].Random();
+   double lmngal = tmass_[ipo].RandomInterp();
+   masse_des_ngal = pow(10.,lmngal) * ngal;
+   ntrial_tab++;
+ }
+
+ return masse_des_ngal;
+}
+
+void SchechterMassDist::Print(void)
+{
+ cout<<"SchechterMassDist::Print: mass=["<<massmin_<<","<<massmax_<<"] n="<<nbinmass_<<endl
+     <<"                 ngal=["<<ngalmin_<<","<<ngalmax_<<"] n="<<nvalngal_<<endl;
+ sch_.Print();
+}
+
+void SchechterMassDist::PrintStatus(void)
+{
+ cout<<"SchechterMassDist::PrintStatus: number of trials: direct="<<ntrial_dir
+     <<" tabulated="<<ntrial_tab<<endl;
+}
+
+void SchechterMassDist::WritePPF(string ppfname)
+{
+ char str[64];
+ cout<<"SchechterMassDist::WritePPF into "<<ppfname<<endl;
+ POutPersist pos(ppfname.c_str());
+
+ double nstar,mstar,alpha;
+ sch_.GetParam(nstar,mstar,alpha);
+ TVector<r_8> tdum(20); tdum = 0.;
+ tdum(0) = nstar;
+ tdum(1) = mstar;
+ tdum(2) = alpha;
+ tdum(3) = sch_outvalue_;
+ tdum(4) = massmin_;
+ tdum(5) = massmax_;
+ tdum(6) = nbinmass_;
+ tdum(7) = ngalmin_;
+ tdum(8) = ngalmax_;
+ tdum(9) = nvalngal_;
+ tdum(10) = hmass_.size();
+ tdum(11) = tmass_.size();
+ pos << PPFNameTag("SMDparam") << tdum;
+
+ pos << PPFNameTag("SMDhmdndm") << *hmdndm_;
+ pos << PPFNameTag("SMDtirhmdndm") << *tirhmdndm_;
+
+ if(hmass_.size()>0) {
+   for(unsigned int i=0;i<hmass_.size();i++) {
+     sprintf(str,"SMDh%d",NGalFrIndex(i));
+     pos << PPFNameTag(str) << hmass_[i];
+   }
+ }
+
+ if(tmass_.size()>0) {
+   for(unsigned int i=0;i<tmass_.size();i++) {
+     sprintf(str,"SMDt%d",NGalFrIndex(i));
+     Histo hdum(tmass_[i]);
+     pos << PPFNameTag(str) << hdum;
+   }
+ }
+
+}
+
+void SchechterMassDist::ReadPPF(string ppfname)
+{
+ Delete(); // On des-alloue si deja remplit!
+
+ char str[64];
+ cout<<"SchechterMassDist::ReadPPF from "<<ppfname<<endl;
+ PInPersist pis(ppfname.c_str());
+
+ TVector<r_8> tdum;
+ pis >> PPFNameTag("SMDparam") >> tdum;
+ sch_.SetParam(tdum(0),tdum(1),tdum(2));
+ sch_.SetOutValue((unsigned short)(tdum(3)+0.1));
+ massmin_ = tdum(4);
+ massmax_ = tdum(5);
+ nbinmass_ = int(tdum(6)+0.1);
+ ngalmin_ = int(tdum(7)+0.1);
+ ngalmax_ = int(tdum(8)+0.1);
+ nvalngal_ = int(tdum(9)+0.1);
+ unsigned int nhmass = (unsigned int)(tdum(10)+0.1);
+ unsigned int ntmass = (unsigned int)(tdum(11)+0.1);
+
+ {
+ Histo hdum;
+ pis >> PPFNameTag("SMDhmdndm") >> hdum;
+ hmdndm_ = new Histo(hdum);
+ pis >> PPFNameTag("SMDtirhmdndm") >> hdum;
+ tirhmdndm_ = new FunRan(hdum,false);
+ }
+
+ if(nhmass>0) {
+   for(unsigned int i=0;i<nhmass;i++) {
+     sprintf(str,"SMDh%d",NGalFrIndex(i));
+     Histo hdum;
+     pis >> PPFNameTag(str) >> hdum;
+     hmass_.push_back(hdum);
+   }
+ }
+
+ if(ntmass>0) {
+   for(unsigned int i=0;i<ntmass;i++) {
+     sprintf(str,"SMDt%d",NGalFrIndex(i));
+     Histo hdum;
+     pis >> PPFNameTag(str) >> hdum;
+     FunRan fdum(hdum,false);
+     tmass_.push_back(fdum);
+   }
+ }
+
+}
+
+
+///////////////////////////////////////////////////////////
+//***************** Fonctions de Check ******************//
+///////////////////////////////////////////////////////////
+
+bool IsCompatible(Schechter& sch1,Schechter& sch2,double eps)
+// on compare les differences a eps pres
+{
+  if(eps<=0.) eps=1.e-4;
+  double nstar1,mstar1,alpha1;
+  sch1.GetParam(nstar1,mstar1,alpha1);
+  double nstar2,mstar2,alpha2;
+  sch2.GetParam(nstar2,mstar2,alpha2);
+
+  if(fabs(nstar1-nstar2)>fabs(nstar1+nstar2)/2.*eps) return false;
+  if(fabs(mstar1-mstar2)>fabs(mstar1+mstar2)/2.*eps) return false;
+  if(fabs(alpha1-alpha2)>fabs(alpha1+alpha2)/2.*eps) return false;
+  return true;
+}
+
 
 ///////////////////////////////////////////////////////////
@@ -98 +456 @@
   return f *d*d / 2.0e-28;
 }
+

trunk/Cosmo/SimLSS/schechter.h

@@ -5 +5 @@
 #include "genericfunc.h"
 
+namespace SOPHYA {
 
-namespace SOPHYA {
+class Histo;
+class FunRan;
 
 //-----------------------------------------------------------------------------------
@@ -13 +15 @@
   Schechter(double nstar,double mstar,double alpha);
   Schechter(Schechter& f);
+  Schechter(void);
   virtual ~Schechter(void);
+
   void SetOutValue(unsigned short outvalue=0);
+  unsigned short GetOutValue(void);
+  void SetParam(double nstar,double mstar,double alpha);
+  void GetParam(double& nstar,double& mstar,double& alpha);
+
   virtual double operator() (double m);
+
+  double Integrate(double massmin,double massmax,int npt=100);
+
   virtual void Print(void);
+
 protected:
   double nstar_,mstar_,alpha_;
@@ -22 +34 @@
 };
 
+//-----------------------------------------------------------------------------------
+class SchechterMassDist : public AnyDataObj {
+  friend class ObjFileIO<SchechterMassDist>;
+public:
+  SchechterMassDist(Schechter sch,double massmin,double massmax,int nbinmass);
+  SchechterMassDist(void);
+  virtual ~SchechterMassDist(void);
+
+  int GetMassLim(double& massmin,double& massmax);
+  int SetNgalLim(int ngalmax,int ngalmin=1,unsigned long nalea=10000);
+  int GetNgalLim(int& ngalmax,int& ngalmin);
+  int GetNgalLim(void) {return nvalngal_;}
+  Schechter GetSchechter(void) {return sch_;}
+
+  inline int IndexFrNGal(int ngal) {
+    int i = ngal-ngalmin_;
+    if(nvalngal_<1 || i<0) return -1;
+    if(i>=nvalngal_) return -2; else return i;
+  }
+  inline int NGalFrIndex(int i) {
+    if(nvalngal_<1 || i<0 || i>=nvalngal_) return -1;
+    return ngalmin_+i;
+  }
+
+  Histo GetHmDnDm(void) const;
+  FunRan GetTmDnDm(void) const;
+
+  Histo GetHisto(int i) const;
+  FunRan GetFunRan(int i) const;
+
+  double TirMass(int ngal);
+
+  void Print(void);
+  void PrintStatus(void);
+
+  void WritePPF(string ppfname);
+  void ReadPPF(string ppfname);
+
+protected:
+ void Delete(void);
+
+ Schechter sch_;
+ unsigned short sch_outvalue_;
+
+ double massmin_,massmax_;   int nbinmass_;
+ int ngalmin_,ngalmax_,nvalngal_;
+ unsigned long ntrial_dir, ntrial_tab;
+
+ Histo* hmdndm_;
+ FunRan* tirhmdndm_;
+
+ vector<Histo> hmass_;
+ vector<FunRan> tmass_;
+};
+
 } // Fin du namespace
+
+//-----------------------------------------------------------------------------------
+bool IsCompatible(Schechter& sch1,Schechter& sch2,double eps=1.e-4);
 
 //-----------------------------------------------------------------------------------