source: Sophya/trunk/SophyaLib/NTools/objfitter.cc@ 1516

#include "objfitter.h"

//===========================================================================
//===========================================================================
//=========================== ObjectFitter ==================================
//===========================================================================
//===========================================================================

//============================= Matrix et Vector ============================
TMatrix<int_4>
ObjectFitter::FitResidus(TMatrix<int_4> const & mtx,GeneralFit& gfit,
                         double xorg,double yorg,double dx,double dy)
{return( ArrayFitter<int_4>::FitResidus(mtx,gfit,xorg,yorg,dx,dy) );}

TMatrix<r_4>
ObjectFitter::FitResidus(TMatrix<r_4> const & mtx,GeneralFit& gfit,
                         double xorg,double yorg,double dx,double dy)
{return( ArrayFitter<r_4>::FitResidus(mtx,gfit,xorg,yorg,dx,dy) );}

TMatrix<r_8>
ObjectFitter::FitResidus(TMatrix<r_8> const & mtx,GeneralFit& gfit,
                         double xorg,double yorg,double dx,double dy)
{return( ArrayFitter<r_8>::FitResidus(mtx,gfit,xorg,yorg,dx,dy) );}

TMatrix<int_4>
ObjectFitter::FitFunction(TMatrix<int_4> const & mtx,GeneralFit& gfit,
                          double xorg,double yorg,double dx,double dy)
{return( ArrayFitter<int_4>::FitFunction(mtx,gfit,xorg,yorg,dx,dy) );}

TMatrix<r_4>
ObjectFitter::FitFunction(TMatrix<r_4> const & mtx,GeneralFit& gfit,
                          double xorg,double yorg,double dx,double dy)
{return( ArrayFitter<r_4>::FitFunction(mtx,gfit,xorg,yorg,dx,dy) );}

TMatrix<r_8>
ObjectFitter::FitFunction(TMatrix<r_8> const & mtx,GeneralFit& gfit,
                          double xorg,double yorg,double dx,double dy)
{return( ArrayFitter<r_8>::FitFunction(mtx,gfit,xorg,yorg,dx,dy) );}

TVector<int_4>
ObjectFitter::FitResidus(TVector<int_4> const & vec,GeneralFit& gfit,
                         double xorg,double dx)
{return( ArrayFitter<int_4>::FitResidus(vec,gfit,xorg,dx) );}

TVector<r_4>
ObjectFitter::FitResidus(TVector<r_4> const & vec,GeneralFit& gfit,
                         double xorg,double dx)
{return( ArrayFitter<r_4>::FitResidus(vec,gfit,xorg,dx) );}

TVector<r_8>
ObjectFitter::FitResidus(TVector<r_8> const & vec,GeneralFit& gfit,
                         double xorg,double dx)
{return( ArrayFitter<r_8>::FitResidus(vec,gfit,xorg,dx) );}

TVector<int_4>
ObjectFitter::FitFunction(TVector<int_4> const & vec,GeneralFit& gfit,
                          double xorg,double dx)
{return( ArrayFitter<int_4>::FitFunction(vec,gfit,xorg,dx) );}

TVector<r_4>
ObjectFitter::FitFunction(TVector<r_4> const & vec,GeneralFit& gfit,
                          double xorg,double dx)
{return( ArrayFitter<r_4>::FitFunction(vec,gfit,xorg,dx) );}

TVector<r_8>
ObjectFitter::FitFunction(TVector<r_8> const & vec,GeneralFit& gfit,
                          double xorg,double dx)
{return( ArrayFitter<r_8>::FitFunction(vec,gfit,xorg,dx) );}

//============================= Matrix et Vector ============================
Image<uint_2> ObjectFitter::FitResidus(Image<uint_2> const & im,GeneralFit& gfit)
{ return( ImageFitter<uint_2>::FitResidus(im,gfit) ); }
Image<int_4> ObjectFitter::FitResidus(Image<int_4> const & im,GeneralFit& gfit)
{ return( ImageFitter<int_4>::FitResidus(im,gfit) ); }
Image<int_8> ObjectFitter::FitResidus(Image<int_8> const & im,GeneralFit& gfit)
{ return( ImageFitter<int_8>::FitResidus(im,gfit) ); }
Image<r_4> ObjectFitter::FitResidus(Image<r_4> const & im,GeneralFit& gfit)
{ return( ImageFitter<r_4>::FitResidus(im,gfit) ); }
Image<r_8> ObjectFitter::FitResidus(Image<r_8> const & im,GeneralFit& gfit)
{ return( ImageFitter<r_8>::FitResidus(im,gfit) ); }

Image<uint_2> ObjectFitter::FitFunction(Image<uint_2> const & im,GeneralFit& gfit)
{ return( ImageFitter<uint_2>::FitFunction(im,gfit) ); }
Image<int_4> ObjectFitter::FitFunction(Image<int_4> const & im,GeneralFit& gfit)
{ return( ImageFitter<int_4>::FitFunction(im,gfit) ); }
Image<int_8> ObjectFitter::FitFunction(Image<int_8> const & im,GeneralFit& gfit)
{ return( ImageFitter<int_8>::FitFunction(im,gfit) ); }
Image<r_4> ObjectFitter::FitFunction(Image<r_4> const & im,GeneralFit& gfit)
{ return( ImageFitter<r_4>::FitFunction(im,gfit) ); }
Image<r_8> ObjectFitter::FitFunction(Image<r_8> const & im,GeneralFit& gfit)
{ return( ImageFitter<r_8>::FitFunction(im,gfit) ); }

//=============================== Histo =====================================
/*! Retourne une classe contenant les residus du fit ``gfit''. */
Histo ObjectFitter::FitResidus(Histo const& hh, GeneralFit& gfit)
{
Histo h(hh);
if(h.NBins()<=0)
  throw(SzMismatchError("Histo::FitResidus: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("Histo::FitResidus: NULL pointer\n"));
TVector<r_8> par = gfit.GetParm();
for(int_4 i=0;i<h.NBins();i++) {
  r_8 x =  h.BinCenter(i);
  h(i) -=  f->Value(&x,par.Data());
}
return h;
}

/*! Retourne une classe contenant la fonction du fit ``gfit''. */
Histo ObjectFitter::FitFunction(Histo const& hh, GeneralFit& gfit)
{
Histo h(hh);
if(h.NBins()<=0)
  throw(SzMismatchError("Histo::FitFunction: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("Histo::FitFunction: NULL pointer\n"));
TVector<r_8> par = gfit.GetParm();
for(int_4 i=0;i<h.NBins();i++) {
  r_8 x =  h.BinCenter(i);
  h(i) =  f->Value(&x,par.Data());
}
return h;
}

/*!
  Fit de l'histogramme par ``gfit''.
  \param errtype,errscale,errmin : pour definir les erreurs
  \sa GeneralFitData::ComputeError(double,err,FitErrType,double,double,bool)
*/
int_4 ObjectFitter::Fit(Histo const& h, GeneralFit& gfit
      ,GeneralFitData::FitErrType errtype,double errscale,double errmin)
{
if(h.NBins()<=0) return -1000;

GeneralFitData mydata(1,h.NBins());

for(int_4 i=0;i<h.NBins();i++) {
  r_8 x =  h.BinCenter(i);
  r_8 f =  h(i);
  r_8 e = (h.HasErrors()) ? h.Error(i) : 1.;
  e = GeneralFitData::ComputeError(f,e,errtype,errscale,errmin);
  mydata.AddData1(x,f,e);
}

gfit.SetData(&mydata);

return gfit.Fit();
}

//============================== Histo2D ====================================
/*! Retourne une classe contenant les residus du fit ``gfit''. */
Histo2D ObjectFitter::FitResidus(Histo2D const& hh, GeneralFit& gfit)
{
Histo2D h(hh);
if(h.NBinX()<=0 || h.NBinY()<=0)
  throw(SzMismatchError("Histo2D::FitResidus: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("Histo2D::FitResidus: NULL pointer\n"));
TVector<r_8> par = gfit.GetParm();
for(int_4 i=0;i<h.NBinX();i++) for(int_4 j=0;j<h.NBinY();j++) {
  r_8 xc,yc;
  h.BinCenter(i,j,xc,yc);
  r_8 x[2] = {xc,yc};
  h(i,j) -= f->Value(x,par.Data());
}
return h;
}

/*! Retourne une classe contenant la fonction du fit ``gfit''. */
Histo2D ObjectFitter::FitFunction(Histo2D const& hh, GeneralFit& gfit)
{
Histo2D h(hh);
if(h.NBinX()<=0 || h.NBinY()<=0)
  throw(SzMismatchError("Histo2D::FitFunction: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("Histo2D::FitFunction: NULL pointer\n"));
TVector<r_8> par = gfit.GetParm();
for(int_4 i=0;i<h.NBinX();i++) for(int_4 j=0;j<h.NBinY();j++) {
  r_8 xc,yc;
  h.BinCenter(i,j,xc,yc);
  r_8 x[2] = {xc,yc};
  h(i,j) = f->Value(x,par.Data());
}
return h;
}

/*!
  Fit de l'histogramme par ``gfit''.
  \param errtype,errscale,errmin : pour definir les erreurs
  \sa GeneralFitData::ComputeError(double,err,FitErrType,double,double,bool)
*/
int_4 ObjectFitter::Fit(Histo2D const & h, GeneralFit& gfit
      ,GeneralFitData::FitErrType errtype,double errscale,double errmin)
{
if(h.NBinX()*h.NBinY()<=0) return -1000;

GeneralFitData mydata(2,h.NBinX()*h.NBinY());

for(int_4 i=0;i<h.NBinX();i++) for(int_4 j=0;j<h.NBinY();j++) {
  r_8 x,y; h.BinCenter(i,j,x,y);
  r_8 f = h(i,j);
  r_8 e = (h.HasErrors()) ? h.Error(i,j) : 1.;
  e = GeneralFitData::ComputeError(f,e,errtype,errscale,errmin);
  mydata.AddData2(x,y,f,e);
}

gfit.SetData(&mydata);

return gfit.Fit();
}

//===========================================================================
//===========================================================================
//========================== ArrayFitter<T> =================================
//===========================================================================
//===========================================================================

/*! Retourne une classe contenant les residus du fit ``gfit''. */
template <class T>
TMatrix<T> 
ArrayFitter<T>::FitResidus(TMatrix<T> const & mtx, GeneralFit& gfit,
                           double xorg,double yorg,double dx,double dy)
// Retourne une classe contenant les residus du fit ``gfit''.
// On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
// Les coordonnees de l'element (i,j) sont :
//    (i,j) -> x = xorg+j*dx , y = yorg+i*dy
{
if(mtx.NCols()<=0||mtx.NRows()<=0)
  throw(SzMismatchError("ArrayFitter::FitResidus(TMatrix<T>...) size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("ArrayFitter::FitResidus(TMatrix<T>...) GeneraFit==NULL\n"));
int npar =  gfit.GetNPar();
if(npar==0)
  throw(SzMismatchError("ArrayFitter::FitResidus(TMatrix<T>...) GeneraFit 0 parametre\n"));
double* par = new double[npar];
{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
TMatrix<T> m(mtx,false);
m.SetTemp(true);
for(uint_4 i=0;i<mtx.NRows();i++) for(uint_4 j=0;j<mtx.NCols();j++) {
  double x[2] = {xorg+j*dx,yorg+i*dy};
  m(i,j) -= f->Value(x,par);
}
delete [] par;
return m;
}


/*! Retourne une classe contenant la fonction du fit ``gfit''. */
template <class T>
TMatrix<T> 
ArrayFitter<T>::FitFunction(TMatrix<T> const & mtx, GeneralFit& gfit,
                            double xorg,double yorg,double dx,double dy)

// Retourne une classe contenant la fonction du fit ``gfit''.
// On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
// Les coordonnees de l'element (i,j) sont :
//    (i,j) -> x = xorg + j*dx , y = yorg + i*dy

{
if(mtx.NCols()<=0||mtx.NRows()<=0)
  throw(SzMismatchError("ArrayFitter::FitFunction(TMatrix<T>...) size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("ArrayFitter::FitFunction(TMatrix<T>...) GeneraFit==NULL\n"));
int npar =  gfit.GetNPar();
if(npar==0)
  throw(SzMismatchError("ArrayFitter::FitFunction(TMatrix<T>...) GeneraFit 0 parametre\n"));
double* par = new double[npar];
{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
TMatrix<T> m(mtx.NRows(), mtx.NCols());
m.SetTemp(true);
for(uint_4 i=0;i<mtx.NRows();i++) for(uint_4 j=0;j<mtx.NCols();j++) {
  double x[2] = {xorg+j*dx,yorg+i*dy};
  m(i,j) = f->Value(x,par);
}
delete [] par;
return m;

}

/*! Retourne une classe contenant les residus du fit ``gfit''. */
template <class T>
TVector<T> 
ArrayFitter<T>::FitResidus(TVector<T> const & vec, GeneralFit& gfit,
                           double xorg,double dx)
// Retourne une classe contenant la fonction fittee du fit ``gfit''.
// La coordonnee de l'element (i) est -> x = xorg + i*dx
{
if(vec.NElts()<=0)
  throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...)  size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit==NULL\n"));
int npar =  gfit.GetNPar();
if(npar==0)
  throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit 0 parametre\n"));
double* par = new double[npar];
{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
TVector<T> v(vec,false);
v.SetTemp(true);
for(uint_4 i=0;i<vec.NElts();i++) {
  double x = xorg+i*dx;
  v(i) -= f->Value(&x,par);
}
delete [] par;
return v;
}

/*! Retourne une classe contenant la fonction du fit ``gfit''. */
template <class T>
TVector<T> 
ArrayFitter<T>::FitFunction(TVector<T> const & vec, GeneralFit& gfit,
                            double xorg,double dx)
// Retourne une classe contenant la function fittee du fit ``gfit''.
// La coordonnee de l'element (i) est -> x = xorg + i*dx
{
if(vec.NElts()<=0)
  throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...)  size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit==NULL\n"));
int npar =  gfit.GetNPar();
if(npar==0)
  throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit 0 parametre\n"));
double* par = new double[npar];
{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
TVector<T> v(vec.NElts());
v.SetTemp(true);
for(uint_4 i=0;i<vec.NElts();i++) {
  double x = xorg+i*dx;
  v(i) = f->Value(&x,par);
}
delete [] par;
return v;
}

//===========================================================================
//===========================================================================
//========================== ImageFitter<T> =================================
//===========================================================================
//===========================================================================

/*! Retourne une classe contenant les residus du fit ``gfit''. */
template <class T>
Image<T> 
ImageFitter<T>::FitResidus(Image<T> const & ima, GeneralFit& gfit)
{
if(ima.XSize()<=0||ima.YSize()<=0)
  throw(SzMismatchError("ImageFitter::FitResidus(Image<T>...) size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("ImageFitter::FitResidus(Image<T>...) GeneraFit==NULL\n"));
int npar =  gfit.GetNPar();
if(npar==0)
  throw(SzMismatchError("ImageFitter::FitResidus(Image<T>...) GeneraFit 0 parametre\n"));
double* par = new double[npar];
{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
Image<T> im(ima,false);
im.SetTemp(true);
for(uint_4 i=0;i<ima.XSize();i++) for(uint_4 j=0;j<ima.YSize();j++) {
  double x[2] = {ima.XPos(i),ima.YPos(j)};
  im(i,j) -= f->Value(x,par);
}
delete [] par;
return im;
}

/*! Retourne une classe contenant la function fittee du fit ``gfit''. */
template <class T>
Image<T> 
ImageFitter<T>::FitFunction(Image<T> const & ima, GeneralFit& gfit)
{
if(ima.XSize()<=0||ima.YSize()<=0)
  throw(SzMismatchError("ImageFitter::FitFunction(Image<T>...) size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("ImageFitter::FitFunction(Image<T>...) GeneraFit==NULL\n"));
int npar =  gfit.GetNPar();
if(npar==0)
  throw(SzMismatchError("ImageFitter::FitFunction(Image<T>...) GeneraFit 0 parametre\n"));
double* par = new double[npar];
{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
Image<T> im(ima.XSize(),ima.YSize());
im.SetTemp(true);
for(uint_4 i=0;i<ima.XSize();i++) for(uint_4 j=0;j<ima.YSize();j++) {
  double x[2] = {ima.XPos(i),ima.YPos(j)};
  im(i,j) = f->Value(x,par);
}
delete [] par;
return im;
}

///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ArrayFitter<int_4>
#pragma define_template ArrayFitter<r_4>
#pragma define_template ArrayFitter<r_8>
#pragma define_template ArrayFitter< complex<r_4> >
#pragma define_template ArrayFitter< complex<r_8> >
#pragma define_template ImageFitter<uint_2>
#pragma define_template ImageFitter<int_4>
#pragma define_template ImageFitter<int_8>
#pragma define_template ImageFitter<r_4>
#pragma define_template ImageFitter<r_8>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ArrayFitter<int_4>;
template class ArrayFitter<r_4>;
template class ArrayFitter<r_8>;
template class ArrayFitter< complex<r_4> >;
template class ArrayFitter< complex<r_8> >;
template class ImageFitter<uint_2>;
template class ImageFitter<int_4>;
template class ImageFitter<int_8>;
template class ImageFitter<r_4>;
template class ImageFitter<r_8>;
#endif