Skip to content

GitLab

Commit 477d8c4f authored by Reza ANSARI's avatar Reza ANSARI
Browse files

Suite implementation fit avec Sophya::GeneraFit() ou Fit avec Minuit, presque... 

Suite implementation fit avec Sophya::GeneraFit() ou Fit avec Minuit, presque fait completement , makefile MAJ, Reza 26/02/2019
parent c1c138a4
Hide whitespace changes
Inline Side-by-side
Showing with 290 additions and 149 deletions
gacfit.h
View file @ 477d8c4f
......@@ -11,6 +11,8 @@
#include "acbeam.h"
#include "tkfit_minuit.h"
static int debug_gacfit=0;
//-----------------------------------------------------------------
......@@ -29,16 +31,18 @@ public:
v_lambdas.resize(v_freqs.size());
for(size_t k=0; k<v_freqs.size(); k++)
v_lambdas[k] = clight/(v_freqs[k]*1.e6);
ndatapoints_=0;
for(size_t j=0; j<v_time_data.size(); j++) ndatapoints_+=v_time_data[j].size();
}
virtual TVector<double> getTimeVec(size_t j)
virtual TVector<double> getTimeVec(size_t j) const
{
vector<double> & tvec = v_time_data[j];
TVector<double> tmv(tvec);
return tmv;
}
virtual int getDataSignal(size_t j, vector<double> & sig)
virtual int getDataSignal(size_t j, vector<double> & sig) const
{
//DEL sig.resize(v_time_data[j].size());
vector<double> & vdata = vv_data[j][nac_];
......@@ -47,7 +51,7 @@ public:
return (int)sig.size();
}
virtual TVector<double> getDataSignal(size_t j)
virtual TVector<double> getDataSignal(size_t j) const
{
vector<double> sig;
getDataSignal(j, sig);
......@@ -57,7 +61,8 @@ public:
}
virtual int getExpectedSignal(size_t j, vector<double> & sig, double D_dish, double theta_b, double phi_b, double A, double B)
virtual int getExpectedSignal(size_t j, vector<double> & sig, double D_dish, double theta_b, double phi_b,
double A, double B) const
{
//DBG if (debug_gacfit > 0) cout<<"*DBG*A*getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<endl;
sig.resize(v_time_data[j].size());
......@@ -84,7 +89,7 @@ public:
return (int)sig.size();
}
virtual int getExpectedSignal(size_t j, vector<double> & sig, const double* parm)
virtual int getExpectedSignal(size_t j, vector<double> & sig, const double* parm) const
{
/*
param[0] = D (Dish diameter)
......@@ -101,7 +106,7 @@ public:
return getExpectedSignal(j, sig, D_dish, theta_b, phi_b, A, B);
}
virtual TVector<double> getExpectedSignal(size_t j, double D_dish, double theta_b, double phi_b, double A, double B)
virtual TVector<double> getExpectedSignal(size_t j, double D_dish, double theta_b, double phi_b, double A, double B) const
{
vector<double> sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
......@@ -112,12 +117,12 @@ public:
return rvs;
}
virtual double getXi2(const double* parm, int & npts)
virtual double getXi2(const double* parm, int & npts) const
{
double xi2=0.;
npts=0;
for(size_t j=0; j<v_time_data.size(); j++) {
if (v_noAC[j]) continue;
// if (v_noAC[j]) continue; On ne peut pas faire ca , il faut modifier sinon la liste des parametres
vector<double> sig;
npts += getExpectedSignal(j, sig, parm);
vector<double> & vdata = vv_data[j][nac_];
......@@ -140,6 +145,7 @@ public:
vector< SLinInterp1D > & v_interp_sazim;
size_t nac_;
bool fggauss_beam_;
size_t ndatapoints_;
};
//-----------------------------------------------------------------
......@@ -149,14 +155,18 @@ public:
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_, vector<bool> & v_noAC_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t nac=0, bool fggauss_beam=true)
: GeneralXi2(3+2*v_time_data_.size()) ,
: GeneralXi2(3+2*v_time_data_.size()) ,
MyACSignal(v_time_data_, vv_data_, vv_err_, v_freqs_, v_noAC_,
v_interp_elev_, v_interp_sazim_, nac, fggauss_beam)
{
}
size_t getNbParams() const { return NPar(); }
size_t getNbDataPoints() const { return ndatapoints_; }
virtual double Value(GeneralFitData& data, double* parm, int& ndataused)
/*
nparam = 3+2*v_time_data_.size()
param[0] = D (Dish diameter)
parm[1] = delta_elev = theta
parm[2] = azim = phi
......@@ -180,6 +190,7 @@ public:
v_interp_elev_, v_interp_sazim_, nac, fggauss_beam)
{
}
virtual double Value(double const parm[])
{
int ndataused;
......@@ -188,16 +199,18 @@ public:
};
#ifdef TKF_AVEC_MINUIT
#ifndef TKF_AVEC_MINUIT
typedef MyACGenXi2 TkF_ACXi2 ;
#else
// Ajustement avec Minuit
class MyAC_Xi2_Minuit : public ROOT::Minuit2::FCNBase , public MyACSignal {
class MyAC_Xi2_Minuit : public Trk_FCNBase , public MyACSignal {
public:
MyAC_Xi2_Minuit(vector< vector<double> > & v_time_data_, vector< vector< vector<double> > > & vv_data_,
vector< vector< vector<double> > > & vv_err_, vector<double> & v_freqs_, vector<bool> & v_noAC_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
size_t nac=0, bool fggauss_beam=true)
: MyACSignal(v_time_data_, vv_data_, vv_err_, v_freqs_, v_noAC_,
v_interp_elev_, v_interp_sazim_, nac, fggauss_beam) , errDef_(1)
v_interp_elev_, v_interp_sazim_, nac, fggauss_beam) , errDef_(1.)
{
nbparam_=3+2*v_time_data_.size();
myparam_ = new double[nbparam_];
......@@ -207,6 +220,8 @@ public:
{
delete[] myparam_;
}
virtual size_t getNbParams() const { return nbparam_; }
virtual size_t getNbDataPoints() const { return ndatapoints_; }
/*
param[0] = D (Dish diameter)
parm[1] = delta_elev = theta
......@@ -229,6 +244,8 @@ public:
double errDef_;
};
typedef MyAC_Xi2_Minuit TkF_ACXi2 ;
#endif /* TKF_AVEC_MINUIT */
#endif
gcxfit.h
View file @ 477d8c4f
......@@ -11,6 +11,8 @@
#include "cxbeam.h"
#include "tkfit_minuit.h"
class MyCxSignal {
public:
MyCxSignal(vector< vector<double> > & v_time_data_, vector< vector< vector< complex<double> > > > & vv_cxdata_,
......@@ -22,16 +24,19 @@ public:
ncx_(ncx), beam0_(beam), phlinfac_(1./250.)
{
//DBG cout << " MyCxSignal() CxCorr=ncx "<<ncx_<< endl;
ndatapoints_=0;
for(size_t j=0; j<v_time_data.size(); j++) ndatapoints_+=v_time_data[j].size();
ndatapoints_*=2;
}
virtual TVector<double> getTimeVec(size_t j)
virtual TVector<double> getTimeVec(size_t j) const
{
vector<double> & tvec = v_time_data[j];
TVector<double> tmv(tvec);
return tmv;
}
virtual int getDataSignal(size_t j, vector< complex<double> > & sig)
virtual int getDataSignal(size_t j, vector< complex<double> > & sig) const
{
//DEL sig.resize(v_time_data[j].size());
vector< complex<double> > & vcxdata = vv_cxdata[j][ncx_];
......@@ -40,7 +45,7 @@ public:
return (int)sig.size();
}
virtual TVector< complex<double> > getDataSignal(size_t j)
virtual TVector< complex<double> > getDataSignal(size_t j) const
{
vector< complex<double> >sig;
getDataSignal(j, sig);
......@@ -51,13 +56,13 @@ public:
//---- This method should be identical to the one in My6CxSignalsB (gcxfitbaseline.h)
// Model for linearly varying phase with frequency : Phase = Phi0 + a_phi * (freq-1250)/250
inline double getPhase4Freq(double phi0, double a_phi, double freq)
inline double getPhase4Freq(double phi0, double a_phi, double freq) const
{
return (phi0+a_phi*(freq-1250.)*phlinfac_);
}
virtual int getExpectedSignal(size_t j, vector< complex<double> > & sig, double phase,
double A, complex<double> B=complex<double>(0.,0.))
double A, complex<double> B=complex<double>(0.,0.)) const
{
//DBG if (debug_gacfit > 0) cout<<"*DBG*A*getExpectedSignal() j="<<j<<" size="<<v_time_data[j].size()<<endl;
sig.resize(v_time_data[j].size());
......@@ -83,7 +88,7 @@ public:
return (int)sig.size();
}
virtual int getExpectedSignal(size_t j, vector< complex<double> > & sig, const double* parm)
virtual int getExpectedSignal(size_t j, vector< complex<double> > & sig, const double* parm) const
{
double phase = getPhase4Freq(parm[0], parm[1], v_freqs[j]);
double A = parm[2+3*j];
......@@ -92,7 +97,7 @@ public:
return getExpectedSignal(j, sig, phase, A);
}
virtual TVector< complex<double> > getExpectedSignal(size_t j, double phase, double A, complex<double> B=complex<double>(0.,0.))
virtual TVector< complex<double> > getExpectedSignal(size_t j, double phase, double A, complex<double> B=complex<double>(0.,0.)) const
{
vector< complex<double> > sig;
//DBG cout << "*DBG*getExpectedSignal() : calling getExpectedSignal(j="<<j<<" ,sig,...)"<<endl;
......@@ -103,7 +108,7 @@ public:
return rvs;
}
virtual double getXi2(const double* parm, int & npts)
virtual double getXi2(const double* parm, int & npts) const
{
double xi2=0.;
npts=0;
......@@ -132,6 +137,7 @@ public:
CxBeam beam0_;
size_t ncx_;
double phlinfac_; // 1./FreqBandWidth = 1./250.
size_t ndatapoints_;
};
//-----------------------------------------------------------------
......@@ -146,6 +152,9 @@ public:
{
}
size_t getNbParams() const { return NPar(); }
size_t getNbDataPoints() const { return ndatapoints_; }
virtual double Value(GeneralFitData& data, double* parm, int& ndataused)
/*
phase = phi_0 + a_phi*(freq-1250.)/250. freq in MHz
......@@ -160,5 +169,46 @@ public:
}
};
#ifndef TKF_AVEC_MINUIT
typedef MyCxGenXi2 TkF_CxXi2 ;
#else
// Ajustement avec Minuit
class MyCx_Xi2_Minuit : public Trk_FCNBase , public MyCxSignal {
public:
MyCx_Xi2_Minuit(vector< vector<double> > & v_time_data_, vector< vector< vector< complex<double> > > > & vv_cxdata_,
vector< vector< vector<double> > > & vv_cxerr_, vector<double> & v_freqs_,
vector< SLinInterp1D > & v_interp_elev_, vector< SLinInterp1D > & v_interp_sazim_,
CxBeam & beam, size_t ncx=0)
: MyCxSignal(v_time_data_, vv_cxdata_, vv_cxerr_, v_freqs_, v_interp_elev_, v_interp_sazim_, beam, ncx), errDef_(1.)
{
nbparam_= 2+3*v_time_data_.size();
myparam_ = new double[nbparam_];
}
virtual ~MyCx_Xi2_Minuit()
{
delete[] myparam_;
}
virtual size_t getNbParams() const { return nbparam_; }
virtual size_t getNbDataPoints() const { return ndatapoints_; }
virtual double operator() (const std::vector< double > & p) const
{
if (p.size() != nbparam_) throw ParmError("MyCx_Xi2_Minuit::operator()(p) p.size() != nbparam_");
for(size_t i=0; i<nbparam_; i++) myparam_[i]=p[i];
int ndataused=0;
return getXi2(myparam_, ndataused);
}
virtual r_8 Up() const { return errDef_;}
size_t nbparam_;
double * myparam_;
double errDef_;
};
typedef MyCx_Xi2_Minuit TkF_CxXi2 ;
#endif /* TKF_AVEC_MINUIT */
#endif
makefile
View file @ 477d8c4f
......@@ -5,21 +5,30 @@ include $(SOPHYABASE)/include/sophyamake.inc
OBJ = ./Objs/
EXE = ./Objs/
#######################
### Pour les programme qui ont besoin de Minuit (mn_trkacxfit) - Definir la variable MINUITDIR
MINUITLIB = -L${MINUITDIR}/lib/
MINUITINC = -I${MINUITDIR}/include/
# List of include files of this package, and .o files to handle dependencies
MYINCLISTHERE = p4autils.h visip4reader.h visp4winreader.h p4gnugain.h p4gvcor.h p4phcor.h p4freqselmgr.h
MYINCLISTHEREFIT = acbeam.h gacfit.h cxbeam.h gcxfit.h gcxfitbaseline.h trkfit.h
MYOLISTHERE = $(OBJ)/p4autils.o $(OBJ)/visip4reader.o $(OBJ)/visp4winreader.o $(OBJ)/p4gnugain.o $(OBJ)/p4gvcor.o $(OBJ)/p4phcor.o $(OBJ)/p4freqselmgr.o
MYOLISTHEREFIT = $(OBJ)/trkfit.o
MYOLISTHEREFIT = $(OBJ)/trkfit.o
MYOLISTHEREFITMN = $(OBJ)/mn_trkfit.o
# Define our target list
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm tfm2dt trkacxfit
all : rdvisip4 visi2ntac visi2dtacx visi2tmfreq tstutls p4conv2fits msvis2dt visiavg filt_blind ckconvphases vis2ra p4vdblist rdthermrfi visamm tfm2dt trkacxfit
clean :
rm -f $(EXE)/rdvisip4 $(EXE)/visi2ntac $(EXE)/visi2dtacx $(EXE)/visi2tmfreq $(EXE)/p4conv2fits $(EXE)/msvis2dt $(EXE)/visiavg \
$(EXE)/filt_blind $(EXE)/ckconvphases $(EXE)/vis2ra $(EXE)/p4vdblist $(EXE)/rdthermrfi $(EXE)/visamm $(EXE)/tfm2dt $(EXE)/trkacxfit
rm -f $(OBJ)/rdvisip4.o $(OBJ)/visi2ntac.o $(OBJ)/visi2dtacx.o $(OBJ)/visi2tmfreq.o $(OBJ)/p4conv2fits.o $(OBJ)/msvis2dt.o $(OBJ)/visiavg.o \
$(OBJ)/filt_blind.o $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o $(OBJ)/p4vdblist.o $(OBJ)/rdthermrfi.o $(OBJ)/visamm.o $(OBJ)/tfm2dt.o $(OBJ)/trkacxfit.o
$(OBJ)/filt_blind.o $(OBJ)/ckconvphases.o $(OBJ)/vis2ra.o $(OBJ)/p4vdblist.o $(OBJ)/rdthermrfi.o $(OBJ)/visamm.o $(OBJ)/tfm2dt.o $(OBJ)/trkacxfit.o $(OBJ)/mn_trkacxfit.o
rm -f $(MYOLISTHERE)
rm -f $(MYOLISTHEREFIT)
rm -f $(MYOLISTHEREFITMN)
depend :
mkdir Objs/
......@@ -98,6 +107,20 @@ $(OBJ)/trkacxfit.o : trkacxfit.cc $(MYINCLISTHERE) $(MYINCLISTHEREFIT)
$(OBJ)/trkfit.o : trkfit.cc $(MYINCLISTHERE) $(MYINCLISTHEREFIT)
$(CXXCOMPILE) -o $(OBJ)/trkfit.o trkfit.cc
## Le meme programme de fit utilisant Minuit
mn_trkacxfit : $(EXE)/mn_trkacxfit
echo '---mn_trkacxfit made'
$(EXE)/mn_trkacxfit : $(OBJ)/mn_trkacxfit.o $(MYOLISTHERE) $(MYOLISTHEREFITMN)
$(CXXLINK) -o $(EXE)/mn_trkacxfit $(OBJ)/mn_trkacxfit.o $(MYOLISTHERE) $(MYOLISTHEREFITMN) $(SOPHYAEXTSLBLIST) ${MINUITLIB} -lMinuit2
$(OBJ)/mn_trkacxfit.o : trkacxfit.cc $(MYINCLISTHERE) $(MYINCLISTHEREFIT)
$(CXXCOMPILE) -DTKF_AVEC_MINUIT ${MINUITINC} -o $(OBJ)/mn_trkacxfit.o trkacxfit.cc
$(OBJ)/mn_trkfit.o : trkfit.cc $(MYINCLISTHERE) $(MYINCLISTHEREFIT)
$(CXXCOMPILE) -DTKF_AVEC_MINUIT ${MINUITINC} -o $(OBJ)/mn_trkfit.o trkfit.cc
## programme pour extraire un DataTable a partir des TFM (fichier PPF), input programme trkacfit
tfm2dt : $(EXE)/tfm2dt
echo '---tfm2dt made'
......
tkfit_minuit.h
View file @ 477d8c4f
......@@ -4,7 +4,19 @@
#ifndef TKFIT_MINUIT_H_SEEN
#define TKFIT_MINUIT_H_SEEN
/Minuit
#include "generalfit.h"
#ifndef TKF_AVEC_MINUIT
class TkF_Fitter : public GeneralFit {
public:
TkF_Fitter(GeneralXi2 & gxi2) : GeneralFit(&gxi2)
{
}
int doFit() { return Fit(); }
};
#else
#include "Math/MinimizerOptions.h"
#include "Minuit2/FCNBase.h"
......@@ -21,138 +33,161 @@
using namespace ROOT::Minuit2;
void DumpMinuitParam(const MnUserParameters& upar) {
int pr = cout.precision();
#define PRECISION 13
#define WIDTH 12
for(std::vector<MinuitParameter>::const_iterator ipar = upar.Parameters().begin();
ipar != upar.Parameters().end(); ipar++) {
cout << std::setw(4) << (*ipar).Number() << std::setw(5) << "||";
cout << std::setw(10) << (*ipar).Name() << std::setw(3) << "||";
if((*ipar).IsConst()) {
cout << " const ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() << " ||" << std::endl;
} else if((*ipar).IsFixed()) {
cout << " fixed ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() << " ||" << std::endl;
} else if((*ipar).HasLimits()) {
if((*ipar).Error() > 0.) {
cout << " limited ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value();
cout << " || " << std::setw(12) << (*ipar).Error() << " ||";
cout << " [" << ((*ipar).HasLowerLimit()?Num2String((*ipar).LowerLimit()):" - ")
<< ", " << ( (*ipar).HasUpperLimit()?Num2String((*ipar).UpperLimit()):" - ") << "]" << std::endl;
} else
cout << " free ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value()
<< " || " << std::setw(12) << "no" << std::endl;
} else {
if((*ipar).Error() > 0.)
cout << " free ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value()
<< " || " << std::setw(12) << (*ipar).Error() << std::endl;
else
cout << " free ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value()
<< " || " << std::setw(12) << "no" << std::endl;
class Trk_FCNBase : public ROOT::Minuit2::FCNBase {
public :
virtual size_t getNbParams() const = 0;
virtual size_t getNbDataPoints() const = 0;
};
class TkF_Fitter {
public:
TkF_Fitter(Trk_FCNBase & gxi2) : gxi2_(gxi2), minxi2_p_(NULL)
{
v_pname.resize(gxi2_.getNbParams());
char buff[32];
for(size_t i=0; i<gxi2_.getNbParams(); i++) {
sprintf(buff,"P%d",(int)i);
v_pname[i]=buff;
upar_.Add(v_pname[i],0.);
}
}//loop on parameter
cout << std::endl;
cout.precision(pr);
#undef PRECISION
#undef WIDTH
max_step_=1000;
tolerance_=0.05;
bestxi2red_=-1.;
rcfit_=-9999;
nstepdone_=0;
}
void SetData(GeneralFitData * gd)
{
return;
}
void SetParam(size_t num, string const & pname, double value, double err, double min=1., double max=-1.)
{
v_pname[num]=pname;
return SetParam(num, value, err, min, max);
}
}
4) Classe de minimisation qui herrite de FCNBase. Par exemple
class AutoCorrChi2 : public FCNBase {
public:
AutoCorrChi2(GeneralFunction* func, GeneralFitData* data): func_(func), errDef_(1.) {
npts_ = data->NData();
for(size_t i=0;i<npts_;i++){
x_.push_back(data->X1(i));
y_.push_back(data->Val(i));
erry_.push_back(data->EVal(i));
bool flag = data->IsValid(i) == 1 ? true: false;
valid_.push_back(flag);
}
void SetParam(size_t num, double value, double err, double min=1., double max=-1.)
{
upar_.SetValue((unsigned int) num, value);
upar_.SetError((unsigned int) num, err);
if (max>min) upar_.SetLimits((unsigned int) num, min, max);
return;
}
//Dtor
virtual ~AutoCorrChi2() {}
//Specialization of the baseclass
virtual r_8 operator() (const vector< r_8 > &par) const {
void SetFix(size_t num, double value)
{
upar_.SetValue((unsigned int) num, value);
upar_.Fix((unsigned int)num);
return;
}
r_8 chi2 =0.;
void SetMaxStep(int maxstep) { max_step_=maxstep; }
for (size_t i=0;i<npts_;i++){
if(valid_[i]) {
static string Num2String(double v)
{
char buff[48];
sprintf(buff,"%lg",v);
return buff;
}
r_8 xp[1]; xp[0] = x_[i];
void PrintFit()
{
int pr = cout.precision();
int PRECISION=13;
int WIDTH=12;
r_8 val = (func_->Value(xp, par.data()) - y_[i])/erry_[i];
chi2 += val*val;
for(std::vector<MinuitParameter>::const_iterator ipar = upar_.Parameters().begin();
ipar != upar_.Parameters().end(); ipar++) {
cout << std::setw(4) << (*ipar).Number() << std::setw(5) << "||";
cout << std::setw(10) << (*ipar).Name() << std::setw(3) << "||";
if((*ipar).IsConst()) {
cout << " const ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() << " ||" << std::endl;
} else if((*ipar).IsFixed()) {
cout << " fixed ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() << " ||" << std::endl;
} else if((*ipar).HasLimits()) {
if((*ipar).Error() > 0.) {
cout << " limited ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value();
cout << " || " << std::setw(12) << (*ipar).Error() << " ||";
cout << " [" << ((*ipar).HasLowerLimit()?Num2String((*ipar).LowerLimit()):" - ")
<< ", " << ( (*ipar).HasUpperLimit()?Num2String((*ipar).UpperLimit()):" - ") << "]" << std::endl;
} else
cout << " free ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value()
<< " || " << std::setw(12) << "no" << std::endl;
} else {
if((*ipar).Error() > 0.)
cout << " free ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value()
<< " || " << std::setw(12) << (*ipar).Error() << std::endl;
else
cout << " free ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value()
<< " || " << std::setw(12) << "no" << std::endl;
}
}
return chi2;
}//loop on parameter
cout << std::endl;
cout.precision(pr);
return;
}
virtual r_8 Up() const { return errDef_;}
private:
size_t npts_; //nbre of points in the fit
GeneralFunction* func_; //function to fit
vector<r_8> x_; //position values
vector<r_8> y_; //measurements
vector<r_8> erry_; //error on measurements
vector<bool> valid_; // if the measurement is valid or not
r_8 errDef_; // 1: 1-sigma Chi2, n*n for n-sigma errors
}; //AutoCorrChi2
5) Ensuite il faut instancier la fonction qui sera donner a Minuit
AutoCorrChi2 fitf(&funcAuto,&dataFit);
funcAuto est la GeneralFunction au sens de Sophya
dataFit sont les donnes au sens de GeneralfitData de Sophya
6) La declartion des parametres
MnUserParameters upar;
upar.Add(nom, initial, err, min, max)
exemple upar.Add("Deff",4.5,0.1,3.,6.);
si on veut fixer un parametre > upar.Fix(nom)
au contraire upar.Release(nom) pour le reliberer
7) Migrad instantiation et minimisation
MnMigrad migrad(fitf, upar);
unsigned int nparFree = migrad.VariableParameters();
cout << "Nbre of free parameters = " << nparFree << endl;
unsigned int maxfcn = 200+100*nparFree + 100*nparFree*nparFree;
FunctionMinimum min = migrad(maxfcn,0.001);
cout << "Minimizer..." << endl;
cout << min << endl;
int doFit()
{
MnMigrad migrad(gxi2_, upar_);
// unsigned int nparFree = migrad.VariableParameters();
// cout << "Nbre of free parameters = " << nparFree << endl;
minxi2_p_ = new FunctionMinimum(migrad(max_step_,tolerance_));
nstepdone_ = rcfit_= minxi2_p_->NFcn();
bestxi2red_=minxi2_p_->Fval()/(double)gxi2_.getNbDataPoints();
// cout << "Minimizer..." << endl;
// cout << min << endl;
if(!minxi2_p_->IsValid()) {
rcfit_ = -rcfit_;
}
return rcfit_;
}
8) recuperation des valeurs fittees
void PrintFitErr(int rcfit)
{
cout << "TkF_Fitter::PrintFitErr() rcfit="<< rcfit << " +/- Nb of iterations , negative -> error "<<endl;
return;
}
Test si le fit a converge
if(!min.IsValid()) {
rc = -1; //TODO be more explicit
cout << "Minuit FAIL: Freq " << freqCentre << endl;
}
ensuite
min.UserParameters().Value("Deff");
double GetParm(unsigned int num) const
{
if ( minxi2_p_)
return minxi2_p_->UserParameters().Value(num);
else return 0.;
}