Changeset 3320 in Sophya for trunk/Cosmo/SimLSS/cmvschdist.cc

Timestamp:

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

Author:

cmv

Message:

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

File:

• trunk/Cosmo/SimLSS/cmvschdist.cc (modified) (9 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 \