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

/* 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+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>
#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 -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[])
{
 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;
 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,"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' :
     Auto_Ini_Ranf(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]);
 }

 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 h200                
disp h100 "same red"
disp h50 "same blue"
disp h10 "same green"
disp h1 "same yellow"
*/