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

/* Distributions de masse pour n=1,2,3,... tirages */
#include "sopnamsp.h"
#include "machdefs.h"
// Pour faire les histogrammes des distributions de masses
//  pour ntirages:
// 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+13 -n -100 -r 200 -N 100000
// compare with:
// > cmvschdist -a -v -m 1e+7,1e+13 -n -100 -r 200 -N 100000 -0
// Fabriquer un fichier pour la prod
// > cmvschdist -a -m 1e+7,1e+13 -n -100 -r 2000 -N 5000000 -W schdist_2000.ppf
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "sophyainit.h"
#include "timing.h"
#include "histos.h"
#include "ntuple.h"
//#include "randr48.h"
#include "randfmt.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 -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
     <<"  -0 : only use ngal=1 histogram"<<endl
     <<"  -a : init auto de l\'aleatoire"<<endl
     <<"  -R readfile.ppf : read SchechterMassDist from ppf file"<<endl;
}

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

 double xmin=1.e7, xmax=1.e13;
 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 = "";

 // --- Choix du generateur aleatoire (a faire ICI imperativement avant AutoInitRand)
 FMTRandGen *RandGen = new FMTRandGen;
 RandGen->SelectGaussianAlgo(C_Gaussian_RandLibSNorm);
 RandGen->SelectPoissonAlgo(C_Poisson_Ahrens);
 RandomGeneratorInterface::SetGlobalRandGenP(RandGen);

 // --- decodage des arguments
 char c;
 while((c = getopt(narg,arg,"h0avm:N:n:r:R:W:")) != -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 'v' :
     verif = true;
     break;
   case '0' :
     useonly1 = true;
     break;
   case 'R' :
     readfrppf = true;
     namefrppf = optarg;
     break;
   case 'W' :
     nametoppf = optarg;
     break;
   case 'a' :
     AutoInitRand(5);
     break;
   case 'h' :
   default :
     usage(); return -1;
   }
 }

 double lnx1=log10(xmin), lnx2=log10(xmax);
 cout<<"xmin="<<xmin<<" ("<<lnx1<<") xmax="<<xmax<<" ("<<lnx2<<")"<<endl;
 cout<<"npt="<<npt<<",  nalea="<<nalea<<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;


 //-------m*dn/dm
 cout<<"> Schechter m*dn/dm nstar="<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
 Schechter sch(nstar,mstar,alpha);

 //------- 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();
 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 "<<ngmax-ngmin+1<<" Histos for back-verif"<<endl;
   for(int i=ngmin;i<=ngmax;i++) {
     Histo h(hmdndm); h.Zero();
     vthisto.push_back(h);
     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++) {
     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");
 }

 //------- 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 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]);
 }

 delete RandGen;

 PrtTim("End of Job");

 return 0;
}

/*
delobjs *
openppf cmvschdist.ppf

set n 70
echo ${h70.vmax}

set h h$n
set t t$n
set th th$n

disp tirhmdndm
disp $t "same red"

disp $h
disp $th "same red"

n/plot $h.val%x ! ! "connectpoints"
n/plot $th.val%x ! ! "nsta connectpoints same red"

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 \
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; \
}

# Compare evolution
disp h2000               
disp h1000 "same red"        
disp h500 "same blue"        
disp h200 "same green"     
disp h100 "same red"
disp h50 "same blue"
disp h10 "same green"
disp h1 "same yellow"
*/