#include "sopnamsp.h" #include "machdefs.h" #include #include #include #include #include #include #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; namespace ROOT { namespace Minuit2 { class MyFCN : public FCNBase { public: MyFCN(HistoErr& herr,PkEisenstein& pk); void Init(const vector& par) const; double Func(double x,const vector& par) const; double operator()(const vector& par) const; double Up(void) const {return 1.;} private: HistoErr& herr_; PkEisenstein& 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 ..."<=narg) return -2; // --- Create spectrum cout< sigma="<> PPFNameTag("herrconc") >> herrconc; } cout<>> Reading["<> PPFNameTag("hpkrec") >> herr; if(herr.NBins()!=herrconc.NBins()) { cout<<"... bad number of bins: "<klim[1]) {herr.SetErr2(i,-1.); continue;} herr.SetErr2(i,herrconc.Error2(i)); nbinok++; } // --- Initialisation de minuit MnUserParameters upar; upar.Add("B",b_ini,b_ini/100.); upar.Add("A",1.,0.01); 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 parfit = ustate.Parameters().Params(); vector eparfit = ustate.Parameters().Errors(); if(ustate.HasCovariance()) { MnUserCovariance ucov = ustate.Covariance(); for(unsigned int i=0;i0) pos.PutObject(nt,"nt"); if(ncov>0) {Cov /= (double)ncov; pos.PutObject(Cov,"cov");} return 0; } //-------------------------------------------------- MyFCN::MyFCN(HistoErr& herr,PkEisenstein& pk) : herr_(herr) , pk_(pk) { } void MyFCN::Init(const vector& par) const { double A=par[1], Ocdm0=par[2], Ob0=par[3], Ol0=par[4], h100=par[5], ns=par[6]; pk_.GetPkIni().SetSlope(ns); pk_.GetPkIni().SetNorm(A); pk_.GetTransfert().SetParTo(h100,Ocdm0,Ob0); pk_.GetGrowthFactor().SetParTo(Ocdm0+Ob0,Ol0); } double MyFCN::Func(double x,const vector& 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& par) const { Init(par); double xi2 = 0.; for(int i=0;i cor(cov,false); cor = 0.; \ for(int i=0;i