source: Sophya/trunk/Cosmo/SimLSS/cmvfitpk.cc@ 3380

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

#include "ntuple.h"
#include "histerr.h"

#include "constcosmo.h"
#include "pkspectrum.h"

#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MinosError.h"
#include "Minuit2/FCNBase.h"
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;
  double Up(void) const {return 1.;}
private:
  HistoErr& herr_;
  GenericFunc& pk_;
};
} }  // namespace ROOT + Minuit2

void usage(void);
void usage(void)
{
cout<<"usage: cmvfitpk [-k kmin,kmax] -e conchpkrec.ppf observ3d_1.ppf observ3d_2.ppf ..."<<endl;
}

int main(int narg,char *arg[])
{
 // --- Cosmography definition (WMAP)
 double zref = 0.175;
 double ob0 = 0.0444356;
 double h100=0.71, om0=0.267804, ol0=0.73;

 // --- Spectrum and variance definition
 double ns = 1., as = 1.;
 double R=8./h100, sigmaR = 1.;
 double kmin=1e-5,kmax=1000.;
 int npt = 10000;
 double eps=1.e-3;

 // --- les fichiers
 double klim[2]={0.,-1.};
 string fconc = "";

 // --- Reading arguments
 char c;
 while((c = getopt(narg,arg,"hk:e:")) != -1) {
  switch (c) {
  case 'k' :
    sscanf(optarg,"%lf,%lf",&klim[0],&klim[1]);
    break;
  case 'e' :
    fconc = optarg;
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }

 if(optind>=narg) return -2;

 // --- Create spectrum
 cout<<endl<<"\n--- Create Spectrum"<<endl;
 InitialSpectrum pkini(ns,as);
 bool nobaryon = false;
 double ocdm0 = om0-ob0;
 TransfertEisenstein tf(h100,ocdm0,ob0,T_CMB_Par,nobaryon);
 GrowthFactor growth(om0,ol0);
 PkSpectrum0 pk0(pkini,tf);
 PkSpectrumZ pkz(pk0,growth,zref);

 // --- Compute variance and normalize spectrum
 pkz.SetZ(0.);
 cout<<endl<<"\n--- Compute variance for top-hat R="<<R
     <<" at z="<<pkz.GetZ()<<endl;
 VarianceSpectrum varpk_th(pkz,R,VarianceSpectrum::TOPHAT);
 double kfind_th = varpk_th.FindMaximum(kmin,kmax,eps);
 double pkmax_th = varpk_th(kfind_th);
 cout<<"kfind_th = "<<kfind_th<<" ("<<log10(kfind_th)<<"), integrand="<<pkmax_th<<endl;
 double k1=kmin, k2=kmax;
 int rc = varpk_th.FindLimits(pkmax_th/1.e4,k1,k2,eps);
 cout<<"limit_th: rc="<<rc<<" : "<<k1<<" ("<<log10(k1)<<") , "
     <<k2<<" ("<<log10(k2)<<")"<<endl;

 double ldlk = (log10(k2)-log10(k1))/npt;
 varpk_th.SetInteg(0.01,ldlk,-1.,4);
 double sr2 = varpk_th.Variance(k1,k2);
 cout<<"varpk_th="<<sr2<<"  ->  sigma="<<sqrt(sr2)<<endl;

 double normpkz = sigmaR*sigmaR/sr2;
 pkz.SetScale(normpkz);
 cout<<"Spectrum normalisation = "<<pkz.GetScale()<<endl;
 pkz.SetZ(zref);

 // --- Reading reference mean spectrum
 cout<<endl<<"\n--- Reading spectrum form file "<<fconc<<endl;
 HistoErr herrconc;
 {
 PInPersist pis(fconc.c_str());
 pis >> PPFNameTag("herrconc") >> herrconc;
 }
 cout<<herrconc;
 if(herrconc.NBins()<1) return -3;
 double b_ini = herrconc(2*herrconc.NBins()/3);
 cout<<"b_ini = "<<b_ini<<endl;
 if(klim[1]<=klim[0]) klim[1] = herrconc.XMax();

 // --- NTuple et ppf
 POutPersist pos("cmvfitpk.ppf");
 const int nvar = 6;
 char *vname[nvar] = {"xi2","a","b","ea","eb","eab"};
 NTuple nt(nvar,vname);
 double xnt[nvar];

 // ***
 // --- Boucle sur les fichiers a fitter
 // ***
 int nfile = 0;
 for(int ifile=optind;ifile<narg;ifile++) {

   // --- lecture du fichier
   cout<<">>> Reading["<<nfile<<","<<ifile<<"]: "<<arg[ifile]<<endl;
   HistoErr herr;
   PInPersist pis(arg[ifile]);
   pis >> PPFNameTag("hpkrec") >> herr;
   if(herr.NBins()!=herrconc.NBins()) {
     cout<<"... bad number of bins: "<<herr.NBins()<<endl;
     continue;
   }

   // --- Mise en forme
   unsigned int nbinok = 0;
   for(int i=0;i<herr.NBins();i++) {
     double k = herr.BinCenter(i);
     if(i==0 || k<klim[0] || k>klim[1]) {herr.SetErr2(i,-1.); continue;}
     herr.SetErr2(i,herrconc.Error2(i));
     nbinok++;
   }

   // --- Initialisation de minuit
    MnUserParameters upar;
    const int npar = 2;
    upar.Add("A",1.,0.01);
    upar.Add("B",b_ini,b_ini/100.);
    MyFCN fcn(herr,pkz);
    if(nbinok<upar.VariableParameters()) {
      cout<<"... nombre de degres de liberte <0: nbinok="<<nbinok<<endl;
      continue;
    }

   // --- Fit minuit
   unsigned int nfcall = 1000; double tolerance = 0.05;
   MnMigrad migrad(fcn,upar,MnStrategy(1));
   FunctionMinimum min = migrad(nfcall,tolerance);
   if(!min.IsValid()) {
     cout<<"!!!!!!!!! Echec Minuit"<<endl;
     //cout<<"minimum: "<<min<<endl;
     continue;
   }

   MnUserParameterState ustate = min.UserState();

   // --- 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);
   nt.Fill(xnt);

   // --- stoquage des histos
   if(nfile<10) {
     HistoErr herrf(herr); herrf.Zero();
     HistoErr herrd(herr); herrd.Zero();
     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.);
     }
     char str[64];
     sprintf(str,"h%d",nfile); pos.PutObject(herr,str);
     sprintf(str,"hf%d",nfile); pos.PutObject(herrf,str);
     sprintf(str,"hd%d",nfile); pos.PutObject(herrd,str);
   }

   // --- Print de debug
   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;
     }
   }

   nfile++;
 }

 // --- Fin du job
 if(nfile>1) pos.PutObject(nt,"nt");

 return 0;
}

//--------------------------------------------------
MyFCN::MyFCN(HistoErr& herr,GenericFunc& pk)
  : herr_(herr) , pk_(pk)
{
}

double MyFCN::operator()(const std::vector<double>& par) const
{
 double xi2 = 0.;
 for(int i=0;i<herr_.NBins();i++) {
   double e2 = herr_.Error2(i);
   if(e2<=0.) continue;
   double x = herr_.BinCenter(i);
   double v = herr_(i);
   double f = par[0]*pk_(x) + par[1];
   xi2 += (v-f)*(v-f)/e2;
 }
 return xi2;
}

/*
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"

 */