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

#include "objfitter.h"

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

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<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<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<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) );
}


//=============================== 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''.
  \verbatim
    typ_err = 0 :
       - erreur attachee au bin si elle existe
       - sinon 1
    typ_err = 1 :
       - erreur attachee au bin si elle existe
       - sinon max( sqrt(abs(bin) ,1 )
    typ_err = 2 :
       - erreur forcee a 1
    typ_err = 3 :
       - erreur forcee a max( sqrt(abs(bin) ,1 )
    typ_err = 4 :
       - erreur forcee a 1, nulle si bin a zero.
    typ_err = 5 :
       - erreur forcee a max( sqrt(abs(bin) ,1 ),
         nulle si bin a zero.
  \endverbatim
*/
int_4 ObjectFitter::Fit(Histo const& h, GeneralFit& gfit,unsigned short typ_err)
{
if(h.NBins()<=0) return -1000;
if(typ_err>5) typ_err=0;

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 saf = sqrt(fabs( f)); if(saf<1.) saf=1.;
  r_8 e=0.;
  if(typ_err==0)      {if(h.HasErrors()) e=h.Error(i); else e=1.;}
  else if(typ_err==1) {if(h.HasErrors()) e=h.Error(i); else e=saf;}
  else if(typ_err==2) e=1.;
  else if(typ_err==3) e=saf;
  else if(typ_err==4) e=(f==0.)?0.:1.;
  else if(typ_err==5) e=(f==0.)?0.:saf;
  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''.
  \verbatim
    typ_err = 0 :
       - erreur attachee au bin si elle existe
       - sinon 1
    typ_err = 1 :
       - erreur attachee au bin si elle existe
       - sinon max( sqrt(abs(bin) ,1 )
    typ_err = 2 :
       - erreur forcee a 1
    typ_err = 3 :
       - erreur forcee a max( sqrt(abs(bin) ,1 )
    typ_err = 4 :
       - erreur forcee a 1, nulle si bin a zero.
    typ_err = 5 :
       - erreur forcee a max( sqrt(abs(bin) ,1 ),
         nulle si bin a zero.
  \endverbatim
*/
int_4 ObjectFitter::Fit(Histo2D const & h, GeneralFit& gfit,unsigned short typ_err)
{
if(h.NBinX()*h.NBinY()<=0) return -1000;
if(typ_err>5) typ_err=0;

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 saf = sqrt(fabs(f)); if(saf<1.) saf=1.;
  r_8 e=0.;
  if(typ_err==0)      {if(h.HasErrors()) e=h.Error(i,j); else e=1.;}
  else if(typ_err==1) {if(h.HasErrors()) e=h.Error(i,j); else e=saf;}
  else if(typ_err==2) e=1.;
  else if(typ_err==3) e=saf;
  else if(typ_err==4) e=(f==0.)?0.:1.;
  else if(typ_err==5) e=(f==0.)?0.:saf;
  mydata.AddData2(x,y,f,e);
}

gfit.SetData(&mydata);

return gfit.Fit();
}

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

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);
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;
}


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;

}

template <class T>
TVector<T> 
ArrayFitter<T>::FitResidus(TVector<T> const & vec, GeneralFit& gfit,
                         double xorg,double dx)
// Retourne une classe contenant les residus 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);
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;
}

template <class T>
TVector<T> 
ArrayFitter<T>::FitFunction(TVector<T> const & vec, GeneralFit& gfit,
                          double xorg,double dx)
// Retourne une classe contenant les residus 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;
}


///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template ArrayFitter<r_4>
#pragma define_template ArrayFitter<r_8>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class ArrayFitter<r_4>;
template class ArrayFitter<r_8>;
#endif