Changeset 3381 in Sophya for trunk/Cosmo/SimLSS/cmvfitpk.cc

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

File:

• trunk/Cosmo/SimLSS/cmvfitpk.cc (modified) (13 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"
+
  */