source: Sophya/trunk/Cosmo/SimLSS/cmvschdist.cc@ 3315

/* Distributions de masse pour n=1,2,3,... tirages */
#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "timing.h"
#include "histos.h"
#include "ntuple.h"
#include "srandgen.h"
#include "perandom.h"

#include "geneutils.h"
#include "cosmocalc.h"
#include "schechter.h"

#include <vector>

void usage(void);
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
     <<"  -a : init auto de l\'aleatoire"<<endl;
}

int main(int narg,char *arg[])
{
 double h75 = 0.71 / 0.75;
 double nstar = 0.006*pow(h75,3.);  //  
 double mstar = pow(10.,9.8/(h75*h75));  // MSol
 double alpha = -1.37;
 cout<<"h75= "<<h75<<" nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;

 double xmin=1e8, xmax=1e13;
 int npt = 100;
 int ngmax = 200, ngmin = 1;
 unsigned long long nalea = 10000;

 char c;
 while((c = getopt(narg,arg,"ham:N:n:r:")) != -1) {
   switch (c) {
   case 'm' :
     sscanf(optarg,"%lf,%lf",&xmin,&xmax);
     if(xmin<=0.) xmin=1e6;
     if(xmax<=0.) xmax=xmin;
     if(xmin>xmax) xmax=10.*xmin;
     break;
   case 'r' :
     sscanf(optarg,"%d,%d",&ngmax,&ngmin);
     if(ngmin<=0) ngmin=1;
     if(ngmax<=0) ngmax=ngmin;
     if(ngmin>ngmax) ngmin=ngmax;
     break;
   case 'n' :
     sscanf(optarg,"%d",&npt);
     if(npt<=0) npt = 100;
     break;
   case 'N' :
     sscanf(optarg,"%llu",&nalea);
     if(nalea<=0) nalea=10000;
     break;
   case 'a' :
     Auto_Ini_Ranf(5);
     break;
   case 'h' :
   default :
     usage(); return -1;
   }
 }

 double lnx1=log10(xmin), lnx2=log10(xmax), dlnx = (lnx2-lnx1)/npt;
 cout<<"xmin="<<xmin<<" ("<<lnx1<<") xmax="<<xmax<<" ("<<lnx2<<"), dlnx="<<dlnx<<endl;
 cout<<"npt="<<npt<<",  nalea="<<nalea<<endl;
 cout<<"ngmin="<<ngmin<<" ngmax="<<ngmax<<endl;


 //-------m*dn/dm
 cout<<"> Schechter m*dn/dm nstar="<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
 Schechter sch(nstar,mstar,alpha);
 sch.SetOutValue(1);
 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 = new Histo(hmdndm); h->Zero();
   vhisto.push_back(h);
   char str[32]; sprintf(str,"h%d",i);
   vhname.push_back(string(str));
 }

 //------- Random
 InitTim();
 cout<<"> Random: "<<nalea<<" trials"<<endl;
 int lpmod = nalea/25; if(lpmod<=0) lpmod=1;
 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++) {
     //sum += tirhmdndm.Random();
     sum += tirhmdndm.RandomInterp();
     int ipo = i-ngmin;
     if(ipo<0) continue;
     vhisto[ipo]->Add(sum/(double)i);
   }
   if(ia%lpmod==0) PrtTim(" ");
 }
 PrtTim("End Random loop");

 //------- Ecriture ppf
 cout<<"Ecriture"<<endl;
 string tag = "cmvschdist.ppf";
 POutPersist pos(tag);
 tag = "hmdndm"; pos.PutObject(hmdndm,tag);
 Histo hdum2(tirhmdndm);
 tag = "tirhmdndm"; pos.PutObject(hdum2,tag);
 for(unsigned int i=0;i<vhisto.size();i++)
   if(vhisto[i]->NEntries()>0) pos.PutObject(*vhisto[i],vhname[i]);

 //------- des-allocation
 for(unsigned int i=0;i<vhisto.size();i++) delete vhisto[i];

 return 0;
}

/*
delobjs *
openppf cmvschdist.ppf

set h h100
disp $h
n/plot $h.val%x ! ! "connectpoints"
n/plot $h.log10(val)%x val>0. ! "connectpoints"

echo ${h100.vmax}

c++exec \
for(int i=0;i<$h.NBins();i++) \
  if($h(i)>0.) { \
  cout<<i<<" m="<<$h.BinCenter(i)<<" ("<<pow(10.,$h.BinCenter(i))<<") h="<<$h(i)<<endl; \
  break; \
} \
for(int i=$h.NBins()-1;i>=0;i--) \
  if($h(i)>0.) { \
  cout<<i<<" m="<<$h.BinCenter(i)<<" ("<<pow(10.,$h.BinCenter(i))<<") h="<<$h(i)<<endl; \
  break; \
}

*/