Changeset 3320 in Sophya for trunk/Cosmo/SimLSS/schechter.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/schechter.cc (modified) (6 diffs)

trunk/Cosmo/SimLSS/schechter.cc

@@ -8 +8 @@
 
 #include "pexceptions.h"
+#include "histos.h"
+#include "perandom.h"
+#include "tvector.h"
+#include "fioarr.h"
 
 #include "constcosmo.h"
+#include "geneutils.h"
 #include "schechter.h"
 
@@ -27 +32 @@
 Schechter::Schechter(Schechter& f)
   : nstar_(f.nstar_) , mstar_(f.mstar_) , alpha_(f.alpha_) , outvalue_(f.outvalue_)
+{
+}
+
+Schechter::Schechter(void)
+  : nstar_(0.) , mstar_(0.) , alpha_(0.) , outvalue_(0)
 {
 }
@@ -45 +55 @@
 }
 
+unsigned short Schechter::GetOutValue(void)
+{
+  return outvalue_;
+}
+
+void Schechter::SetParam(double nstar,double mstar,double alpha)
+{
+  nstar_ = nstar; mstar_ = mstar; alpha_ = alpha;
+}
+
+void Schechter::GetParam(double& nstar,double& mstar,double& alpha)
+{
+  nstar = nstar_; mstar = mstar_; alpha = alpha_;
+}
+
 double Schechter::operator() (double m)
 // Return : "dn/dm = f(m)" or "m*dn/dm = f(m)"
@@ -52 +77 @@
   if(outvalue_) dndm *= m;  // on veut m*dn/dm
   return dndm;
+}
+
+
+double Schechter::Integrate(double massmin,double massmax,int npt)
+// Integrate from massmin to massmax with at least npt points log10 spaced
+{
+ if(massmin<=0. || massmax<=0. || massmax<=massmin) return 0.;
+ if(npt<1) npt = 100;
+ double lmassmin = log10(massmin), lmassmax = log10(massmax);
+ double perc=0.01, dlxinc=(lmassmax-lmassmin)/npt, dlxmax=10.*dlxinc; unsigned short glorder=4;
+ double sum = IntegrateFuncLog(*this,lmassmin,lmassmax,perc,dlxinc,dlxmax,glorder);
+ return sum;
 }
 
@@ -63 +100 @@
   cout<<endl;
 }
+
+///////////////////////////////////////////////////////////
+//******* Les facilites pour tirer sur Schechter ********//
+///////////////////////////////////////////////////////////
+
+SchechterMassDist::SchechterMassDist(Schechter sch,double massmin,double massmax,int nbinmass)
+// On veut une fonction de tirage aleatoire de la masse de N galaxies sur une Schechter
+// Pour optimiser on cree (eventuellement) des histos de tirage
+// ATTENTION: on donne les limites en masses pour les histos mais leurs abscisses
+//   sont en log10(masse): histo [log10(massmin),log10(massmax)] avec nbinmass bins
+// Si nbinmass<0 alors c'est le nombre de points par decade
+  : sch_(sch) , sch_outvalue_(sch.GetOutValue())
+  , massmin_(massmin) , massmax_(massmax) , nbinmass_(nbinmass)
+  , ngalmin_(0) , ngalmax_(0) , nvalngal_(0)
+  , ntrial_dir(0) , ntrial_tab(0)
+  , hmdndm_(NULL) , tirhmdndm_(NULL)
+{
+  if(massmin_>massmax_  || massmin_<=0. || massmax_<=0.|| nbinmass_==0) {
+    cout<<"SchechterMassDist::SchechterMassDist: error in input values"<<endl;
+    throw ParmError("SchechterMassDist::SchechterMassDist: error in input values");
+  }
+
+  // Creation du tirage aleatoire sur la Schechter
+  sch_outvalue_ = sch.GetOutValue();  // on sauve la configuration de la schechter
+  sch.SetOutValue(1);  // on veut m*dN/dm
+  double lnx1 = log10(massmin_), lnx2 = log10(massmax_); // Binning en log10 de la masse
+  if(nbinmass_<0) nbinmass_ = int((-nbinmass_)*(lnx2-lnx1+1.));
+  hmdndm_ = new Histo(lnx1,lnx2,nbinmass_); hmdndm_->ReCenterBin();
+  FuncToHisto(sch,*hmdndm_,true);  // true -> bin en log10(x)
+  tirhmdndm_ = new FunRan(*hmdndm_,true);  // true -> histo est pdf
+  sch.SetOutValue(sch_outvalue_);  // on remet a la valeur initiale
+}
+
+SchechterMassDist::SchechterMassDist(void)
+  : sch_outvalue_(0)
+  , massmin_(0.) , massmax_(0.) , nbinmass_(0)
+  , ngalmin_(0) , ngalmax_(0) , nvalngal_(0)
+  , ntrial_dir(0) , ntrial_tab(0)
+  , hmdndm_(NULL) , tirhmdndm_(NULL)
+{
+}
+
+SchechterMassDist::~SchechterMassDist(void)
+{
+ Delete();
+}
+
+void SchechterMassDist::Delete(void)
+{
+ if(hmdndm_) delete hmdndm_;
+ if(tirhmdndm_) delete tirhmdndm_;
+ hmass_.resize(0);
+ tmass_.resize(0);
+ ntrial_dir = ntrial_tab = 0;
+}
+
+int SchechterMassDist::GetMassLim(double& massmin,double& massmax)
+{
+ massmin = massmin_;
+ massmax = massmax_;
+ return nbinmass_;
+}
+
+int SchechterMassDist::SetNgalLim(int ngalmax,int ngalmin,unsigned long nalea)
+// Creation des Histos de tirage pour ngalmin a ngalmax galaxies
+{
+ int lp=2;
+ ngalmin_=ngalmax_=nvalngal_=0;
+ if(ngalmin<=0) ngalmin=1;
+ if(ngalmax<ngalmin || ngalmax==1) return 0;
+ ngalmin_ = ngalmin;
+ ngalmax_ = ngalmax;
+ nvalngal_ = ngalmax-ngalmin+1;
+  
+ if(nalea<1) nalea = 100000;
+
+ if(lp>0) cout<<"SchechterMassDist::SetNgalLim: ngal=["
+              <<ngalmin_<<","<<ngalmax_<<"] n="<<nvalngal_
+              <<" filling with "<<nalea<<" trials"<<endl;
+
+ //------- Construct histo
+ double lnx1 = log10(massmin_), lnx2 = log10(massmax_);
+ if(lp>0) cout<<"> Creating "<<nvalngal_<<" histos ["<<lnx1<<","<<lnx2<<"] n="<<nbinmass_<<endl;
+ for(int i=ngalmin_;i<=ngalmax_;i++) {
+   Histo h(*hmdndm_); h.Zero();
+   hmass_.push_back(h);
+ }
+ if(lp>1) cout<<"...number of histos is "<<hmass_.size()<<endl;
+
+ //------- Remplissage random
+ sch_.SetOutValue(1);  // on veut m*dN/dm
+ int lpmod = nalea/20; if(lpmod<=0) lpmod=1;
+ double s1=0., sc1=0.; unsigned long ns1=0;
+ double sax=0., scax=0.; unsigned long nsax=0;
+ for(unsigned long ia=0;ia<nalea;ia++) {
+   if(lp>1 && ia%lpmod==0) cout<<"...tirage "<<ia<<endl;
+   double sum = 0.;
+   for(int i=1;i<=ngalmax_;i++) {
+     //double l10m = tirhmdndm_->Random();
+     double l10m = tirhmdndm_->RandomInterp();
+     double m = pow(10.,l10m);
+     sum += m;
+     s1 += m; sc1 += m*m; ns1++;
+     int ipo = i-ngalmin_;
+     if(ipo<0) continue;
+     // ATTENTION: l'histo de tirage stoque le log10(moyenne=somme_masses/ngal).
+     //            Ca permet d'avoir le meme binning quelque soit ngal
+     double v = log10(sum/(double)i);
+     hmass_[ipo].Add(v);
+     if(i==ngalmax) {sax += sum/(double)i; scax += sum/(double)i*sum/(double)i; nsax++;}
+   }
+ }
+ sch_.SetOutValue(sch_outvalue_);  // on remet a la valeur initiale
+
+ if(ns1>1) {
+   s1 /= ns1; sc1 = sc1/ns1 - s1*s1;
+   cout<<"...Mean mass for ngal=1: "<<s1<<" ("<<log10(fabs(s1))<<")"
+       <<" s="<<sqrt(fabs(sc1))<<" (ntrials "<<ns1<<")"<<endl;
+ }
+ if(nsax>1) {
+   sax /= nsax; scax = scax/nsax - sax*sax;
+   cout<<"...Mean mass for ngal="<<ngalmax_<<": "<<sax<<" ("<<log10(fabs(sax))<<")"
+       <<" s="<<sqrt(fabs(scax))<<" (ntrials "<<nsax<<")"<<endl;
+ }
+
+ //------- Generation des classes de tirage aleatoire et des histos de verif
+ if(lp>0) cout<<"> Creating "<<nvalngal_<<" FunRan"<<endl;
+ for(unsigned int i=0;i<hmass_.size();i++) {
+   FunRan t(hmass_[i],true);
+   tmass_.push_back(t);
+ }
+ if(lp>1) cout<<"...number of funran is "<<tmass_.size()<<endl;
+
+ return nvalngal_;
+}
+
+int SchechterMassDist::GetNgalLim(int& ngalmax,int& ngalmin)
+{
+ ngalmax = ngalmax_;
+ ngalmin = ngalmin_;
+ return nvalngal_;
+}
+
+Histo SchechterMassDist::GetHmDnDm(void) const
+{
+  return *hmdndm_;
+}
+
+FunRan SchechterMassDist::GetTmDnDm(void) const
+{
+  return *tirhmdndm_;
+}
+
+Histo SchechterMassDist::GetHisto(int i) const
+{
+  if(i<0 || i>=nvalngal_) {
+    cout<<"SchechterMassDist::GetHisto: error in input values"<<endl;
+    throw ParmError("SchechterMassDist::GetHisto: error in input values");
+  }
+  return hmass_[i];
+}
+
+FunRan SchechterMassDist::GetFunRan(int i) const
+{
+  if(i<0 || i>=nvalngal_) {
+    cout<<"SchechterMassDist::GetFunRan: error in input values"<<endl;
+    throw ParmError("SchechterMassDist::GetFunRan: error in input values");
+  }
+  return tmass_[i];
+}
+
+double SchechterMassDist::TirMass(int ngal)
+{
+ if(ngal<1) return 0.;
+
+ int ipo = IndexFrNGal(ngal);
+ double masse_des_ngal = 0.;
+ if(ipo<0) {  // Pas d'histos de tirage pour ce nombre de galaxies
+   for(long i=0;i<ngal;i++) {  // On tire donc ngal fois sur la Schechter
+     // double lm = tirhmdndm_->Random();
+     double lm = tirhmdndm_->RandomInterp();
+     masse_des_ngal += pow(10.,lm);  // ATTENTION abscisse en log10(masse)
+   }
+   ntrial_dir++;
+ } else {
+   // ATTENTION l'histo de tirage stoque le log10(moyenne=somme_masses/ngal)
+   //double lmngal = tmass_[ipo].Random();
+   double lmngal = tmass_[ipo].RandomInterp();
+   masse_des_ngal = pow(10.,lmngal) * ngal;
+   ntrial_tab++;
+ }
+
+ return masse_des_ngal;
+}
+
+void SchechterMassDist::Print(void)
+{
+ cout<<"SchechterMassDist::Print: mass=["<<massmin_<<","<<massmax_<<"] n="<<nbinmass_<<endl
+     <<"                 ngal=["<<ngalmin_<<","<<ngalmax_<<"] n="<<nvalngal_<<endl;
+ sch_.Print();
+}
+
+void SchechterMassDist::PrintStatus(void)
+{
+ cout<<"SchechterMassDist::PrintStatus: number of trials: direct="<<ntrial_dir
+     <<" tabulated="<<ntrial_tab<<endl;
+}
+
+void SchechterMassDist::WritePPF(string ppfname)
+{
+ char str[64];
+ cout<<"SchechterMassDist::WritePPF into "<<ppfname<<endl;
+ POutPersist pos(ppfname.c_str());
+
+ double nstar,mstar,alpha;
+ sch_.GetParam(nstar,mstar,alpha);
+ TVector<r_8> tdum(20); tdum = 0.;
+ tdum(0) = nstar;
+ tdum(1) = mstar;
+ tdum(2) = alpha;
+ tdum(3) = sch_outvalue_;
+ tdum(4) = massmin_;
+ tdum(5) = massmax_;
+ tdum(6) = nbinmass_;
+ tdum(7) = ngalmin_;
+ tdum(8) = ngalmax_;
+ tdum(9) = nvalngal_;
+ tdum(10) = hmass_.size();
+ tdum(11) = tmass_.size();
+ pos << PPFNameTag("SMDparam") << tdum;
+
+ pos << PPFNameTag("SMDhmdndm") << *hmdndm_;
+ pos << PPFNameTag("SMDtirhmdndm") << *tirhmdndm_;
+
+ if(hmass_.size()>0) {
+   for(unsigned int i=0;i<hmass_.size();i++) {
+     sprintf(str,"SMDh%d",NGalFrIndex(i));
+     pos << PPFNameTag(str) << hmass_[i];
+   }
+ }
+
+ if(tmass_.size()>0) {
+   for(unsigned int i=0;i<tmass_.size();i++) {
+     sprintf(str,"SMDt%d",NGalFrIndex(i));
+     Histo hdum(tmass_[i]);
+     pos << PPFNameTag(str) << hdum;
+   }
+ }
+
+}
+
+void SchechterMassDist::ReadPPF(string ppfname)
+{
+ Delete(); // On des-alloue si deja remplit!
+
+ char str[64];
+ cout<<"SchechterMassDist::ReadPPF from "<<ppfname<<endl;
+ PInPersist pis(ppfname.c_str());
+
+ TVector<r_8> tdum;
+ pis >> PPFNameTag("SMDparam") >> tdum;
+ sch_.SetParam(tdum(0),tdum(1),tdum(2));
+ sch_.SetOutValue((unsigned short)(tdum(3)+0.1));
+ massmin_ = tdum(4);
+ massmax_ = tdum(5);
+ nbinmass_ = int(tdum(6)+0.1);
+ ngalmin_ = int(tdum(7)+0.1);
+ ngalmax_ = int(tdum(8)+0.1);
+ nvalngal_ = int(tdum(9)+0.1);
+ unsigned int nhmass = (unsigned int)(tdum(10)+0.1);
+ unsigned int ntmass = (unsigned int)(tdum(11)+0.1);
+
+ {
+ Histo hdum;
+ pis >> PPFNameTag("SMDhmdndm") >> hdum;
+ hmdndm_ = new Histo(hdum);
+ pis >> PPFNameTag("SMDtirhmdndm") >> hdum;
+ tirhmdndm_ = new FunRan(hdum,false);
+ }
+
+ if(nhmass>0) {
+   for(unsigned int i=0;i<nhmass;i++) {
+     sprintf(str,"SMDh%d",NGalFrIndex(i));
+     Histo hdum;
+     pis >> PPFNameTag(str) >> hdum;
+     hmass_.push_back(hdum);
+   }
+ }
+
+ if(ntmass>0) {
+   for(unsigned int i=0;i<ntmass;i++) {
+     sprintf(str,"SMDt%d",NGalFrIndex(i));
+     Histo hdum;
+     pis >> PPFNameTag(str) >> hdum;
+     FunRan fdum(hdum,false);
+     tmass_.push_back(fdum);
+   }
+ }
+
+}
+
+
+///////////////////////////////////////////////////////////
+//***************** Fonctions de Check ******************//
+///////////////////////////////////////////////////////////
+
+bool IsCompatible(Schechter& sch1,Schechter& sch2,double eps)
+// on compare les differences a eps pres
+{
+  if(eps<=0.) eps=1.e-4;
+  double nstar1,mstar1,alpha1;
+  sch1.GetParam(nstar1,mstar1,alpha1);
+  double nstar2,mstar2,alpha2;
+  sch2.GetParam(nstar2,mstar2,alpha2);
+
+  if(fabs(nstar1-nstar2)>fabs(nstar1+nstar2)/2.*eps) return false;
+  if(fabs(mstar1-mstar2)>fabs(mstar1+mstar2)/2.*eps) return false;
+  if(fabs(alpha1-alpha2)>fabs(alpha1+alpha2)/2.*eps) return false;
+  return true;
+}
+
 
 ///////////////////////////////////////////////////////////
@@ -98 +456 @@
   return f *d*d / 2.0e-28;
 }
+