Changeset 3381 in Sophya for trunk

Timestamp:

Nov 14, 2007, 7:25:34 PM (18 years ago)

Author:

cmv

Message:

1-/ pkspectrum: GrowthFactor + PkSpectrum on enleve les memorisations

des calculs precedents, tout est recalcule

2-/ cmvfitpk: mise de Ob+Ol+Oc+h100+ns dans le fit

cmv 14/11/2007

Location:

trunk/Cosmo/SimLSS

Files:

• cmvfitpk.cc (modified) (13 diffs)
• pkspectrum.cc (modified) (9 diffs)
• pkspectrum.h (modified) (5 diffs)

trunk/Cosmo/SimLSS/cmvfitpk.cc

@@ -23 +23 @@
 using namespace ROOT::Minuit2;
 
-#include "Minuit2/FCNBase.h"
 namespace ROOT { namespace Minuit2 {
 class MyFCN : public FCNBase {
 public:
-  MyFCN(HistoErr& herr,GenericFunc& pk);
-  double operator()(const std::vector<double>& par) const;
+  MyFCN(HistoErr& herr,PkSpectrumZ& pk);
+  void Init(const vector<double>& par) const;
+  double Func(double x,const vector<double>& par) const;
+  double operator()(const vector<double>& par) const;
   double Up(void) const {return 1.;}
 private:
   HistoErr& herr_;
-  GenericFunc& pk_;
+  PkSpectrumZ& pk_;
 };
 } }  // namespace ROOT + Minuit2
@@ -59 +60 @@
  double klim[2]={0.,-1.};
  string fconc = "";
+
+ // --- les options de print/plot
+ int nprint = 3, nplot = -1;  // -1 = all
 
  // --- Reading arguments
@@ -84 +88 @@
  double ocdm0 = om0-ob0;
  TransfertEisenstein tf(h100,ocdm0,ob0,T_CMB_Par,nobaryon);
- GrowthFactor growth(om0,ol0);
+ GrowthFactor growth(om0,ol0); cout<<"...Growth="<<growth(zref)<<endl;
  PkSpectrum0 pk0(pkini,tf);
  PkSpectrumZ pkz(pk0,growth,zref);
@@ -125 +129 @@
 
  // --- NTuple et ppf
+ const int npar = 7;  // Number of parameters in the fit
+ TMatrix<r_8> Cov(npar,npar); Cov = 0.;
  POutPersist pos("cmvfitpk.ppf");
- const int nvar = 6;
- char *vname[nvar] = {"xi2","a","b","ea","eb","eab"};
+ const int nvar = 2*npar+1;
+ char *vname[nvar] = {
+       "xi2"
+       ,"b","a","oc","ob","ol","h","n"
+       ,"eb","ea","eoc","eob","eol","eh","en"
+       };
  NTuple nt(nvar,vname);
  double xnt[nvar];
@@ -134 +144 @@
  // --- Boucle sur les fichiers a fitter
  // ***
- int nfile = 0;
+ int nfile=0, ncov=0;
  for(int ifile=optind;ifile<narg;ifile++) {
 
@@ -158 +168 @@
    // --- Initialisation de minuit
     MnUserParameters upar;
-    const int npar = 2;
+    upar.Add("B",b_ini,b_ini/100.);
     upar.Add("A",1.,0.01);
-    upar.Add("B",b_ini,b_ini/100.);
+    upar.Add("Ocdm",ocdm0,0.001,ocdm0/2.,2.*ocdm0);
+    //upar.Fix("Ocdm");
+    upar.Add("Ob",ob0,0.001,ob0/10.,10.*ob0);
+    //upar.Fix("Ob");
+    upar.Add("Ol",ol0,0.001,0.1,2.);
+    upar.Fix("Ol");
+    upar.Add("h100",h100,0.01,10.,150.);
+    upar.Fix("h100");
+    upar.Add("ns",ns,0.001,0.7,1.5);
+    upar.Fix("ns");
     MyFCN fcn(herr,pkz);
     if(nbinok<upar.VariableParameters()) {
@@ -178 +197 @@
 
    MnUserParameterState ustate = min.UserState();
+   double xi2r = ustate.Fval()/(nbinok-ustate.VariableParameters());
+   vector<double> parfit  = ustate.Parameters().Params();
+   vector<double> eparfit = ustate.Parameters().Errors();
+   if(ustate.HasCovariance()) {
+     MnUserCovariance ucov = ustate.Covariance();
+     for(unsigned int i=0;i<upar.VariableParameters();i++)
+       for(unsigned int j=0;j<upar.VariableParameters();j++)
+         Cov(ustate.ExtOfInt(i),ustate.ExtOfInt(j)) += ucov(i,j);
+     ncov++;
+   }
 
    // --- 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] = 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);
+   for(int i=0;i<npar;i++) xnt[1+i] = parfit[i];
+   for(int i=0;i<npar;i++) xnt[npar+1+i] = eparfit[i];
    nt.Fill(xnt);
 
-   // --- stoquage des histos
-   if(nfile<10) {
+   // --- stoquage des histos de plots
+   if(nfile<nplot || nplot<0) {
      HistoErr herrf(herr); herrf.Zero();
      HistoErr herrd(herr); herrd.Zero();
+     fcn.Init(parfit);
      for(int i=1;i<herr.NBins();i++) {
-       double f = A*pkz(herr.BinCenter(i))+B;
+       double f = fcn.Func(herr.BinCenter(i),parfit);
        herrf.SetBin(i,f,herr.Error(i),1.);
        herrd.SetBin(i,herr(i)-f,herr.Error(i),1.);
@@ -208 +232 @@
 
    // --- Print de debug
-   if(nfile<3) {
+   if(nfile<nprint) {
      cout<<"minimum: "<<min<<endl;
      for(int i=0;i<nvar;i++) cout<<"["<<i<<"] = "<<xnt[i]<<endl;
@@ -225 +249 @@
 
  // --- Fin du job
- if(nfile>1) pos.PutObject(nt,"nt");
+ if(nfile>0) pos.PutObject(nt,"nt");
+ if(ncov>0) {Cov /= (double)ncov; pos.PutObject(Cov,"cov");}
 
  return 0;
@@ -231 +256 @@
 
 //--------------------------------------------------
-MyFCN::MyFCN(HistoErr& herr,GenericFunc& pk)
+MyFCN::MyFCN(HistoErr& herr,PkSpectrumZ& pk)
   : herr_(herr) , pk_(pk)
 {
 }
 
-double MyFCN::operator()(const std::vector<double>& par) const
-{
+void MyFCN::Init(const vector<double>& par) const
+{
+ double A=par[1], Ocdm0=par[2], Ob0=par[3], Ol0=par[4], h100=par[5], ns=par[6];
+ pk_.GetPk0().GetPkIni().SetSlope(ns);
+ pk_.GetPk0().GetPkIni().SetNorm(A);
+ pk_.GetPk0().GetTransfert().SetParTo(h100,Ocdm0,Ob0);
+ pk_.GetGrowthFactor().SetParTo(Ocdm0+Ob0,Ol0);
+}
+
+double MyFCN::Func(double x,const vector<double>& par) const
+// par: [0]=B,  [1]=A,  [2]=Ocdm0,  [3]=Ob0,  [4]=Ol0, [5]=h100, [6]=npow
+{
+ Init(par);
+ return pk_(x) + par[0];
+}
+
+double MyFCN::operator()(const vector<double>& par) const
+{
+ Init(par);
  double xi2 = 0.;
  for(int i=0;i<herr_.NBins();i++) {
@@ -244 +286 @@
    double x = herr_.BinCenter(i);
    double v = herr_(i);
-   double f = par[0]*pk_(x) + par[1];
+   double f = Func(x,par);
    xi2 += (v-f)*(v-f)/e2;
  }
@@ -252 +294 @@
 /*
 openppf cmvfitpk.ppf
+
+print nt
 
 set i 0
@@ -270 +314 @@
 n/plot nt.b%a ! ! "crossmarker5"
 
+n/plot nt.oc
+n/plot nt.ob
+n/plot nt.oc+ob
+n/plot nt.ob%oc ! ! "crossmarker5"
+
 n/plot nt.ea
 n/plot nt.eb
 n/plot nt.ea%eb ! ! "crossmarker5"
 
+disp cov "zoomx30"
+del cor
+c++exec \
+TMatrix<r_8> cor(cov,false); cor = 0.; \
+for(int i=0;i<cor.NRows();i++) { \
+  for(int j=0;j<cor.NCols();j++) { \
+    double v = cov(i,i)*cov(j,j); \
+    if(v<=0.) continue; \
+    cor(i,j) = cov(i,j)/sqrt(v); \
+} } \
+KeepObj(cor); \
+cout<<"Matrice cor cree "<<endl;
+
+disp cor "zoomx30"
+
  */

trunk/Cosmo/SimLSS/pkspectrum.cc

@@ -31 +31 @@
 InitialSpectrum::~InitialSpectrum(void)
 {
-}
-
-void InitialSpectrum::SetNorm(double a)
-{
-  A_ = a;
-}
-
-void InitialSpectrum::SetSlope(double n)
-{
-  n_ = n;
 }
 
@@ -177 +167 @@
 }
 
-bool TransfertEisenstein::SetParTo(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon)
+bool TransfertEisenstein::SetParTo(double h100,double OmegaCDM0,double OmegaBaryon0)
 // Changement des valeurs des parametres (suivi de re-init eventuel)
+// Si h100,Omega...<=0. alors pas de changement, on garde l'ancienne valeur
 {
  bool haschanged = false;
@@ -185 +176 @@
  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
@@ -371 +360 @@
 // Pour avoir D(z) = 1/(1+z) faire: OmegaMatter0=1 OmegaLambda0=0
 GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
-  : O0_(OmegaMatter0) , Ol_(OmegaLambda0) , Ok_(1.-OmegaMatter0-OmegaLambda0)
+  : O0_(OmegaMatter0) , Ol_(OmegaLambda0)
 {
   if(OmegaMatter0==0.) {
@@ -377 +366 @@
     throw ParmError("GrowthFactor::GrowthFactor: Error badOmegaMatter0  value");
   }
-  norm_ = 1.; // puisque (*this)(0.) a besoin de norm_
-  norm_ = (*this)(0.);
-  cout<<"GrowthFactor::GrowthFactor : norm="<<norm_<<endl;
 }
 
 GrowthFactor::GrowthFactor(GrowthFactor& d1)
-  : O0_(d1.O0_) , Ol_(d1.Ol_) , Ok_(d1.Ok_) , norm_(d1.norm_)
+  : O0_(d1.O0_) , Ol_(d1.Ol_)
 {
 }
@@ -396 +382 @@
  z += 1.;
  double z2 = z*z, z3 = z2*z;
- double den = Ol_ + Ok_*z2 + O0_*z3;
+
+ // Calcul de la normalisation (pour z=0 -> growth=1.)
+ double D1z0 = pow(O0_,4./7.) - Ol_ + (1.+O0_/2.)*(1.+Ol_/70.);
+ D1z0 = 2.5*O0_ / D1z0;
+
+ // Calcul du growthfactor pour z
+ double Ok = 1. - O0_ - Ol_;
+ double den = Ol_ + Ok*z2 + O0_*z3;
  double o0z = O0_ *z3 / den;
  double olz = Ol_ / den;
 
- // 4./7. = 0.571429
- double D1z = pow(o0z,0.571429) - olz + (1.+o0z/2.)*(1.+olz/70.);
+ double D1z = pow(o0z,4./7.) - olz + (1.+o0z/2.)*(1.+olz/70.);
  D1z = 2.5*o0z / z / D1z;
 
- 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;
+ return D1z / D1z0;
+}
+
+void GrowthFactor::SetParTo(double OmegaMatter0,double OmegaLambda0)
+{
+ if(OmegaMatter0>0.) O0_ = OmegaMatter0;
+ Ol_ = OmegaLambda0;
+}
+
+bool GrowthFactor::SetParTo(double OmegaMatter0)
+// idem precedent sans changer OmegaLambda0
+{
+ if(OmegaMatter0<=0.) return false;
+ O0_ = OmegaMatter0;
+ return true;
 }
 
@@ -453 +442 @@
 PkSpectrumZ::PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref)
   : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(PK)
-  , zold_(-1.) , d1old_(1.)
 {
 }
@@ -459 +447 @@
 PkSpectrumZ::PkSpectrumZ(PkSpectrumZ& pkz)
   : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(PK)
-  , zold_(pkz.zold_) , d1old_(pkz.d1old_)
 {
 }
@@ -481 +468 @@
 double PkSpectrumZ::operator() (double k,double z)
 {
-  double d1;
-  if(z == zold_) d1 = d1old_;
-    else {d1 = d1old_ = d1_(z); zold_ = z;}
+  double d1 = d1_(z);
 
   double v = pk0_(k) * d1*d1;

trunk/Cosmo/SimLSS/pkspectrum.h

@@ -14 +14 @@
   virtual ~InitialSpectrum(void);
   virtual double operator() (double k) {return A_ * pow(k,n_);}
-  void SetNorm(double a);
-  void SetSlope(double n);
+  void SetNorm(double a) {A_ = a;}
+  void SetSlope(double n) {n_ = n;}
 protected:
   double n_, A_;
@@ -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);
+  bool SetParTo(double h100,double OmegaCDM0,double OmegaBaryon0);
   virtual double operator() (double k);
   double KPeak(void);
@@ -77 +77 @@
   virtual ~GrowthFactor(void);
   virtual double operator() (double z);
-  bool SetParTo(double OmegaMatter0=-1.,double OmegaLambda0=-12345.);
+  void SetParTo(double OmegaMatter0,double OmegaLambda0);
+  bool SetParTo(double OmegaMatter0);
 protected:
-  double O0_,Ol_,Ok_;
-  double norm_;
+  double O0_,Ol_;
 };
 
@@ -107 +107 @@
   virtual double operator() (double k);
   virtual double operator() (double k,double z);
-  inline void   SetZ(double z) {zref_ = z;}
-  inline double GetZ(void) {return zref_;}
+  void   SetZ(double z) {zref_ = z;}
+  double GetZ(void) {return zref_;}
   void SetTypSpec(ReturnSpectrum typspec=PK);
-  inline void SetScale(double scale=1.) {scale_=scale; zold_=-1.;}
-  inline double GetScale(void) {return scale_;}
+  void SetScale(double scale=1.) {scale_=scale;}
+  double GetScale(void) {return scale_;}
   PkSpectrum0& GetPk0(void) {return pk0_;}
   GrowthFactor& GetGrowthFactor(void) {return d1_;}
@@ -119 +119 @@
   double zref_, scale_;
   ReturnSpectrum typspec_;
-  mutable double zold_, d1old_;
 };