Changeset 3378 in Sophya for trunk/Cosmo

Timestamp:

Nov 9, 2007, 7:04:12 PM (18 years ago)

Author:

cmv

Message:

modif sur remplissage NTuple/Minos, possibilite de changer les parametres O0,Om,Ol etc... , cmv 09/11/2007

Location:

trunk/Cosmo/SimLSS

Files:

• Makefile (modified) (6 diffs)
• cmvfitpk.cc (modified) (6 diffs)
• pkspectrum.cc (modified) (12 diffs)
• pkspectrum.h (modified) (3 diffs)

trunk/Cosmo/SimLSS/Makefile

@@ -14 +14 @@
 MYLIB = $(SOPHYAEXTSLBLIST) -L$(LIB) -lcmvsimbao -lfftw3_threads -lfftw3 -lm
 
+#--------------------------------------------------------------------------
+# ---- Les programmes de simulation
 PROGS = \
        $(EXE)cmvtuniv $(EXE)cmvtransf $(EXE)cmvtgrowth $(EXE)cmvtstpk \
@@ -34 +36 @@
 LIBR = $(LIB)libcmvsimbao.a
 
+#--------------------------------------------------------------------------
+# ---- Les programmes de test
 PROGSTEST = \
         $(EXE)cmvtluc $(EXE)cmvchkwhu $(EXE)hu_sigma8
@@ -41 +45 @@
         $(OBJ)cmvtluc.o $(OBJ)cmvchkwhu.o $(OBJ)hu_sigma8.o
 
+#--------------------------------------------------------------------------
+#---- Les programmes utilisant des librairies non standard SOPHYA
+SPROGS = \
+       $(EXE)cmvfitpk
+
+SPROGSOBJ = \
+       $(OBJ)cmvfitpk.o
+
 #----
 all: lib prog
@@ -49 +61 @@
 
 progtest: $(PROGSTEST)
+
+sprog: $(SPROGS)
+
+allprog: all progtest sprog
 
 clean:
@@ -56 +72 @@
         rm -rf $(OBJ)/CmvBAO_cxxrep/
         rm -f $(PROGSTEST) $(PROGSTESTOBJ)
+        rm -f $(SPROGS) $(SPROGSOBJ)
 
 ##############################################################################
@@ -215 +232 @@
 $(OBJ)cmvtluc.o: cmvtluc.cc
         $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvtluc.cc
+
+##############################################################################
+cmvfitpk: $(EXE)cmvfitpk
+        echo $@ " done"
+$(EXE)cmvfitpk: $(OBJ)cmvfitpk.o $(LIBR)
+        $(CXXLINK) $(CXXREP) -o $@ $(OBJ)cmvfitpk.o $(MYLIB) -lMinuit2Base
+$(OBJ)cmvfitpk.o: cmvfitpk.cc
+        $(CXXCOMPILE) $(CXXREP) -I$(MYEXTINC) -o $@ cmvfitpk.cc

trunk/Cosmo/SimLSS/cmvfitpk.cc

@@ -124 +124 @@
  // --- NTuple et ppf
  POutPersist pos("cmvfitpk.ppf");
- const int nvar = 10;
- char *vname[nvar] = {"xi2","a","b","ea","eb","eab","ea1","ea2","eb1","eb2"};
+ const int nvar = 6;
+ char *vname[nvar] = {"xi2","a","b","ea","eb","eab"};
  NTuple nt(nvar,vname);
  double xnt[nvar];
@@ -156 +156 @@
    // --- Initialisation de minuit
     MnUserParameters upar;
+    const int npar = 2;
     upar.Add("A",1.,0.01);
     upar.Add("B",b_ini,b_ini/100.);
@@ -176 +177 @@
    MnUserParameterState ustate = min.UserState();
 
-   MnMinos Minos(fcn,min);
-   pair<double,double> eminos[2];
-   for(int ip=0;ip<2;ip++) {
-     eminos[ip].first = eminos[ip].second = 0.;
-     if(!ustate.Parameter(ip).IsFixed()) eminos[ip] = Minos(ip);
-   }
-
    // --- Recuperation des informations et remplissage NTuple
    for(int i=0;i<nvar;i++) xnt[i] = 0.;
    double xi2r = ustate.Fval()/(nbinok-upar.VariableParameters());
+   double A = ustate.Value((unsigned int)0);
+   double B = ustate.Value((unsigned int)1);
    xnt[0] = xi2r;
-   xnt[1] = ustate.Value((unsigned int)0);
-   xnt[2] = ustate.Value((unsigned int)1);
+   xnt[1] = A;
+   xnt[2] = B;
    xnt[3] = ustate.Error((unsigned int)0);
    xnt[4] = ustate.Error((unsigned int)1);
    if(ustate.HasCovariance()) xnt[5] = ustate.Covariance()(0,1);
-   xnt[6] = eminos[0].first;
-   xnt[7] = eminos[0].second;
-   xnt[8] = eminos[1].first;
-   xnt[9] = eminos[1].second;
    nt.Fill(xnt);
 
@@ -202 +194 @@
      HistoErr herrf(herr); herrf.Zero();
      HistoErr herrd(herr); herrd.Zero();
-     for(int i=0;i<herr.NBins();i++) {
-       double f = pkz(herr.BinCenter(i));
+     for(int i=1;i<herr.NBins();i++) {
+       double f = A*pkz(herr.BinCenter(i))+B;
        herrf.SetBin(i,f,herr.Error(i),1.);
        herrd.SetBin(i,herr(i)-f,herr.Error(i),1.);
@@ -214 +206 @@
 
    // --- Print de debug
-   if(nfile==0) {
+   if(nfile<3) {
      cout<<"minimum: "<<min<<endl;
      for(int i=0;i<nvar;i++) cout<<"["<<i<<"] = "<<xnt[i]<<endl;
+
+     MnMinos Minos(fcn,min);
+     pair<double,double> eminos[npar];
+     for(int ip=0;ip<npar;ip++) {
+       eminos[ip].first = eminos[ip].second = 0.;
+       if(!ustate.Parameter(ip).IsFixed()) eminos[ip] = Minos(ip);
+       cout<<"Minos: "<<ip<<" err=["<<eminos[ip].first<<","<<eminos[ip].second<<"]"<<endl;
+     }
    }
 
@@ -248 +248 @@
 }
 
+/*
+openppf cmvfitpk.ppf
+
+set i 0
+zone
+disp h$i "nsta hbincont err"
+disp hf$i "nsta same red"
+
+zone
+disp hd$i "nsta hbincont err red"
+disp hd$i "nsta same"
+
+zone
+n/plot nt.xi2
+
+zone
+n/plot nt.a
+n/plot nt.b
+n/plot nt.b%a ! ! "crossmarker5"
+
+n/plot nt.ea
+n/plot nt.eb
+n/plot nt.ea%eb ! ! "crossmarker5"
+
+ */

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -51 +51 @@
 TransfertEisenstein::TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon,int lp)
   : lp_(lp)
-  , Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , tcmb_(tcmb)
+  , Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h100_(h100) , tcmb_(tcmb)
   , nobaryon_(nobaryon) , nooscenv_(0), retpart_(ALL)
 {
@@ -60 +60 @@
 TransfertEisenstein::TransfertEisenstein(TransfertEisenstein& tf)
   : lp_(tf.lp_)
-  ,Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , tcmb_(tf.tcmb_)
+  ,Oc_(tf.Oc_) , Ob_(tf.Ob_) , h100_(tf.h100_) , tcmb_(tf.tcmb_)
   , nobaryon_(tf.nobaryon_) , nooscenv_(tf.nooscenv_), retpart_(tf.retpart_)
 {
@@ -67 +67 @@
 }
 
+TransfertEisenstein::~TransfertEisenstein(void)
+{
+}
+
 void TransfertEisenstein::zero_(void)
 {
@@ -78 +82 @@
   O0_ = Oc_ + Ob_;
   if(nobaryon_) {O0_ = Oc_; Ob_ = 0.;}
-  double H0 = 100. * h_, h2 = h_*h_;
-  if(lp_) cout<<"h100="<<h_<<" H0="<<H0<<") Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
+  double H0 = 100. * h100_, h2 = h100_*h100_;
+  if(lp_) cout<<"h100="<<h100_<<" H0="<<H0<<") Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;
 
 
@@ -173 +177 @@
 }
 
-TransfertEisenstein::~TransfertEisenstein(void)
-{
+bool TransfertEisenstein::SetParTo(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon)
+// Changement des valeurs des parametres (suivi de re-init eventuel)
+{
+ bool haschanged = false;
+
+ if(h100>0.) {h100_ = h100; haschanged = true;}
+ if(OmegaCDM0>0.) {Oc_ = OmegaCDM0; haschanged = true;}
+ if(OmegaBaryon0>0.) {Ob_ = OmegaBaryon0; haschanged = true;}
+ if(tcmb>0.) {tcmb_ = tcmb; haschanged = true;}
+ if(nobaryon!=nobaryon_) {nobaryon_ = nobaryon; haschanged = true;}
+
+ // et on recalcule les initialisations
+ if(haschanged) Init_();
+
+ return haschanged;
 }
 
@@ -202 +219 @@
 }
 
+void TransfertEisenstein::SetPrintLevel(int lp)
+{
+ lp_ = lp;
+}
+
+
 double TransfertEisenstein::T0tild(double k,double alphac,double betac)
 {
   // Formule 10 p 608
-  //double q = k*th2p7_*th2p7_/(O0_*h_*h_);
+  //double q = k*th2p7_*th2p7_/(O0_*h100_*h100_);
   double q = k/(13.41*keq_);
   // Formule 20 p 610
@@ -220 +243 @@
   //  OU Pour function lissee sans oscillation baryon
   if(nobaryon_  || nooscenv_ == 2) {
-    double gamma = O0_*h_;
+    double gamma = O0_*h100_;
     // Calcul de Gamma_eff, formule 30 page 612 (pour fct lissee)
     if( nobaryon_==false && nooscenv_ == 2 )
-      gamma = O0_*h_*(alphag_ + (1.-alphag_)/(1.+pow(0.43*k*sfit_,4.))); // Gamma_eff
+      gamma = O0_*h100_*(alphag_ + (1.-alphag_)/(1.+pow(0.43*k*sfit_,4.))); // Gamma_eff
     // Formule 28 page 612 : qui est est equivalent a:
-    //    q = k / h_ * th2p7_*th2p7_ / gamma;
+    //    q = k / h100_ * th2p7_*th2p7_ / gamma;
     // qui est est equivalent a:
     //    q = k / (13.41 * keq)                   pour Ob=0
@@ -231 +254 @@
     // Les resultats sont legerement differents a cause des valeurs approx.
     // des constantes numeriques: on prend comme W.Hu (tf_fit.c)
-    //double q = k / h_ * th2p7_*th2p7_ / gamma;  // Mpc^-1
-    double q = k/(13.41*keq_) * (O0_*h_/gamma);  // Mpc^-1
+    //double q = k / h100_ * th2p7_*th2p7_ / gamma;  // Mpc^-1
+    double q = k/(13.41*keq_) * (O0_*h100_/gamma);  // Mpc^-1
     // Formules 29 page 612
     double l0 = log(2.*M_E + 1.8*q);
@@ -284 +307 @@
 
 TransfertTabulate::TransfertTabulate(double h100,double OmegaCDM0,double OmegaBaryon0)
-: Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , kmin_(1.) , kmax_(-1.)
+: Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h100_(h100) , kmin_(1.) , kmax_(-1.)
 , interptyp_(0)
 {
@@ -290 +313 @@
 
 TransfertTabulate::TransfertTabulate(TransfertTabulate& tf)
-: Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , kmin_(tf.kmin_) , kmax_(tf.kmax_)
+: Oc_(tf.Oc_) , Ob_(tf.Ob_) , h100_(tf.h100_) , kmin_(tf.kmin_) , kmax_(tf.kmax_)
 , interptyp_(tf.interptyp_) , k_(tf.k_) , tf_(tf.tf_)
 {
@@ -324 +347 @@
    double k,tc,tb,tf;
    sscanf(line,"%lf %lf %lf",&k,&tc,&tb);
-   k *= h_;     // convert h^-1 Mpc  ->  Mpc
+   k *= h100_;     // convert h^-1 Mpc  ->  Mpc
    tf = (Oc_*tc+Ob_*tb)/(Oc_+Ob_);
    if(tf>tmax) tmax = tf;
@@ -382 +405 @@
 
  return D1z / norm_;
+}
+
+bool GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
+{
+ bool haschanged = false;
+
+ if(OmegaMatter0>0.) {O0_ = OmegaMatter0; haschanged = true;}
+ if(fabs(OmegaLambda0+12345.)>1e-6) {Ol_ = OmegaLambda0; haschanged = true;}
+
+ // et on recalcule les initialisations
+ if(haschanged) {
+   Ok_ = 1. - O0_ - Ol_;
+   norm_ = 1.; // puisque (*this)(0.) a besoin de norm_
+   norm_ = (*this)(0.);
+ }
+
+ return haschanged;
 }
 

trunk/Cosmo/SimLSS/pkspectrum.h

@@ -29 +29 @@
   TransfertEisenstein(TransfertEisenstein& tf);
   virtual ~TransfertEisenstein(void);
+  bool SetParTo(double h100=-1.,double OmegaCDM0=-1.,double OmegaBaryon0=-1.,double tcmb=-1.,bool nobaryon=false);
   virtual double operator() (double k);
   double KPeak(void);
   void SetNoOscEnv(unsigned short nooscenv=0);
   void SetReturnPart(ReturnPart retpart=ALL);
+  void SetPrintLevel(int lp=0);
 protected:
   int lp_;
-  double O0_,Oc_,Ob_,h_,tcmb_;
+  double O0_,Oc_,Ob_,h100_,tcmb_;
   double th2p7_;
   double zeq_,keq_,zd_,Req_,Rd_,s_,ksilk_,alphac_,betac_,bnode_,alphab_,betab_;
@@ -61 +63 @@
   int ReadCMBFast(string filename);
 protected:
-  double Oc_,Ob_,h_;
+  double Oc_,Ob_,h100_;
   double kmin_,kmax_;
   int interptyp_;
@@ -75 +77 @@
   virtual ~GrowthFactor(void);
   virtual double operator() (double z);
+  bool SetParTo(double OmegaMatter0=-1.,double OmegaLambda0=-12345.);
 protected:
   double O0_,Ol_,Ok_;