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

Timestamp:

Sep 3, 2007, 12:49:21 PM (18 years ago)

Author:

cmv

Message:

commit pour memoire avant intro class SchechterMassDist cmv 03/09/2007

File:

• trunk/Cosmo/SimLSS/cmvschdist.cc (modified) (11 diffs)

trunk/Cosmo/SimLSS/cmvschdist.cc

@@ -2 +2 @@
 #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+14 -n 140 -r 200 -N 100000
 #include <iostream>
 #include <stdlib.h>
@@ -28 +33 @@
      <<"  ngmax,ngmin : distribution pour tirage de ngmin a ngmax galaxies"<<endl
      <<"  nalea : nombre de tirages aleatoires"<<endl
+     <<"  -v : back verification of random trials"<<endl
      <<"  -a : init auto de l\'aleatoire"<<endl;
 }
@@ -43 +49 @@
  int ngmax = 200, ngmin = 1;
  unsigned long long nalea = 10000;
+ bool verif = false;
 
  char c;
- while((c = getopt(narg,arg,"ham:N:n:r:")) != -1) {
+ while((c = getopt(narg,arg,"havm:N:n:r:")) != -1) {
    switch (c) {
    case 'm' :
@@ -67 +74 @@
      if(nalea<=0) nalea=10000;
      break;
+   case 'v' :
+     verif = true;
+     break;
    case 'a' :
      Auto_Ini_Ranf(5);
@@ -80 +90 @@
  cout<<"npt="<<npt<<",  nalea="<<nalea<<endl;
  cout<<"ngmin="<<ngmin<<" ngmax="<<ngmax<<endl;
+ cout<<"verif="<<verif<<endl;
 
 
@@ -85 +96 @@
  cout<<"> Schechter m*dn/dm nstar="<<nstar<<" mstar="<<mstar<<" alpha="<<alpha<<endl;
  Schechter sch(nstar,mstar,alpha);
- sch.SetOutValue(1);
+ sch.SetOutValue(1);  // on veut m*dN/dm
  Histo hmdndm(lnx1,lnx2,npt); hmdndm.ReCenterBin();
  FuncToHisto(sch,hmdndm,true);
@@ -93 +104 @@
  int nbhist = ngmax-ngmin+1;
  cout<<"> Creating "<<nbhist<<" histos"<<endl;
- vector<Histo *> vhisto;
- vector<string> vhname;
+ vector<Histo> vhisto;  vector<string> vhname;
  for(int i=ngmin;i<=ngmax;i++) {
-   Histo* h = new Histo(hmdndm); h->Zero();
+   Histo h(hmdndm); h.Zero();
    vhisto.push_back(h);
    char str[32]; sprintf(str,"h%d",i);
@@ -106 +116 @@
  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++) {
-     //sum += tirhmdndm.Random();
-     sum += tirhmdndm.RandomInterp();
+     //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;
-     vhisto[ipo]->Add(sum/(double)i);
+     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(verif) {
+   cout<<"> Creating "<<nbhist<<" FunRan and 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();
+     vthisto.push_back(h);
+     sprintf(str,"th%d",i);
+     vthname.push_back(string(str));
+   }
+
+   cout<<"> Checking tirage"<<endl;
+   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);
+     }
+   }
+   PrtTim("End Random loop back-tirage");
+ }
 
  //------- Ecriture ppf
@@ -127 +173 @@
  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];
+ 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++)
+     if(vthisto[i].NEntries()>0) pos.PutObject(vthisto[i],vthname[i]);
 
  return 0;
@@ -140 +187 @@
 openppf cmvschdist.ppf
 
-set h h100
+set n 70
+echo ${h70.vmax}
+
+set h h$n
+set th th$n
+
 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"
-
-echo ${h100.vmax}
+n/plot $th.log10(val)%x val>0. ! "nsta connectpoints same red"
 
 c++exec \
@@ -159 +214 @@
 }
 
+# Compare evolution
+disp h200                
+disp h100 "same red"
+disp h50 "same blue"
+disp h10 "same green"
+disp h1 "same yellow"
 */