source: Sophya/trunk/SophyaPI/PIext/piafitting.cc@ 808

//CMVBUG reste pb suivant sous Linux-egcs (OSF1-cxx c'est ok):
//openppf cmvwppf.ppf
//fitw h1g g2 p:4500,45,10,900
//delobjs h1g
//clic sur FIT  -> OK
//clic sur FIT 2sd fois -> core dump!!!
//
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <iostream.h>
#include <typeinfo>

#include "strutil.h"
#include "nbtri.h"
#include "generalfit.h"
#include "fct1dfit.h"
#include "fct2dfit.h"
#include "ntuple.h"
#include "cimage.h"
#include "histos.h"
#include "histos2.h"
#include "ntuple.h"
#include "hisprof.h"

#ifdef SANS_EVOLPLANCK
#include "matrix.h"
#include "cvector.h"
#else
#include "tmatrix.h"
#include "tvector.h"
#include "objfitter.h"
#endif

#include "piafitting.h"
#include "nobjmgr.h"
#include "pistdimgapp.h"

////////////////////////////////////////////////////////////////////
////---------- Classe de fenetre de fit interactive ------------////
////////////////////////////////////////////////////////////////////
class PIAFitterWind :  public PIWindow {
public :
                PIAFitterWind(PIStdImgApp* par, PIAFitter* fiter);
  virtual       ~PIAFitterWind();
  virtual void  Show();
  virtual void  Process(PIMessage msg, PIMsgHandler* sender, void* data=NULL);
  void SetState(void);
protected :
  PIStdImgApp* mDap;
  PIAFitter*   mFitter;
  PILabel**    lab; int nlab;
  PIText**     txt; int ntxt;
  PIButton**   but; int nbut;
  PICheckBox** ckb; int nckb;
  PIOptMenu**  pom; int npom;
  bool ReFillGData;
  int saveI1,saveI2, saveJ1, saveJ2;
  double ErrValue;
};

////////////////////////////////////////////////////////////////////////
// --- Generation automatique du "usage" pour helpfit:
// 1-/ copier le texte ci-dessous de piafitting.cc dans toto
// 2-/ cat toto | sed 's?^//?usage+="\\n?' | sed 's?$?    ";?'
// NE PAS EFFACER LA SUITE: c'est elle qui doit etre modifiee,
// le usage du helpfit sera regenere a l'aide des commandes ci-dessus.
////////////////////////////////////////////////////////////////////////
//| --------------- Fit Lineaire a 1 et 2 dimensions ---------------
//| Syntaxe:
//| fitlin nom pnn [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]
//| avec:
//| nom  : cf commentaire ordre ''fit''
//| pnn  : fit polynome degre nn avec classe Poly (lineaire) 1D ou 2D
//| o:aa,...,bb : cf commentaires ordre ''fit''
//|
////////////////////////////////////////////////////////////////////////
//| --------------- Fit d'objets a 1 et 2 dimensions ---------------
//|                 avec interface d'aide graphique
//| Syntaxe:
//| fitw nom func
//|     [p:p1,p2,...,pn s:s1,s2,...,sn m:m1,m2,...,mn M:M1,mM2,...,Mn f:f1,...,fn]
//|     [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]
//| cf commentaire ordre ''fit''
//|
////////////////////////////////////////////////////////////////////////
//| --------------- Fit d'objets a 1 et 2 dimensions ---------------
//| Syntaxe:
//| fit nom func
//|     [p:p1,p2,...,pn s:s1,s2,...,sn m:m1,m2,...,mn M:M1,mM2,...,Mn f:f1,...,fn]
//|     [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]
//|----- OBJET -----
//| nom  : nom de l'objet qui peut etre:
//|        fit-1D:  Vector,Histo1D,HProf ou GeneraFitData(1D)
//|        fit-2D:  Matrix,Histo2D,Image<T> ou GeneraFitData(2D)
//|
//|----- FUNCTION -----
//| func : pnn : fit polynome degre nn avec GeneralFit (non-lineaire) 1D ou 2D
//|      : gnn : fit gaussienne (hauteur) + polynome de degre nn 1D
//|      : g   : fit gaussienne (hauteur) 1D
//|      : enn : fit exponentielle + polynome de degre nn 1D
//|      : e   : fit exponentielle 1D
//|      : Gnn : fit gaussienne (volume) + polynome de degre nn 1D
//|      : G   : fit gaussienne (volume) 1D
//|      :     : fit gaussienne+fond (volume) 2D
//|      : Gi  : fit gaussienne+fond integree (volume) 2D
//|      : d   : fit DL de gaussienne+fond (volume) 2D
//|      : di  : fit DL de gaussienne+fond integree (volume) 2D
//|      : D   : fit DL de gaussienne+fond avec coeff variable p6 (volume) 2D
//|      : Di  : fit DL de gaussienne+fond integree avec coeff variable p6 (volume) 2D
//|      : M   : fit Moffat+fond (expos=p6) (volume) 2D
//|      : Mi  : fit Moffat+fond integree (expos=p6) (volume) 2D
//|      : Autre : fonction definie par l'utilisateur (cf commentaires ci apres)
//|
//|----- INIT PARAMETRES ET ETAT DU FIT -----
//| p    : p1,...,pn : valeur d'initialisation des parametres (def=0)
//| s    : s1,...,sn : valeur des steps de depart (def=1)
//| m    : m1,...,mn : valeur des minima (def=1)
//| M    : M1,...,Mn : valeur des maxima (def=-1) (max<=min : pas de limite)
//| f    : f1,...,fn : si >=1 parametre fixe sinon libre (def=0)
//| - Remarque: si pi,si,mi ou Mi = '!' la valeur correspondante n'est pas changee
//|
//|----- OPTIONS -----
//| o    : options ''o:Eaa.b,eaa.b,f,r,caa.b,Xaa.b''
//|      F : initialisation a partir des resultats et de l'etat du fit precedent
//|          (option non prioritaire sur les definitions p:,s:,m:,M:,f:o:)
//|      a : initialisation automatique des parametres du fit (si disponible)
//|          (option non prioritaire sur les definitions p:,s:,m:,M:,f:o:
//|           et sur l'initialisation pa l'option ''F'')
//|      f : generation d'un Objet identique contenant la fonction fittee
//|      r : generation d'un Objet identique contenant les residus
//|      Xaa.b : aa.b valeur du DXi2 d'arret (def=1.e-3)
//|      Naa : aa nombre maximum d'iterations (def=100)
//|      la : niveau ''a'' de print: a=niveau de print Fit1/2D
//|      da : niveau ''a'' de debug: a=niveau de GeneralFit
//|      Ii1/i2 : numeros des bins X de l'histos utilises pour le fit [i1,i2]
//|2D    Jj1/j2 : numeros des bins Y de l'histos utilises pour le fit [j1,j2]
//|    - L'erreur est celle associee a l'objet (si elle existe),
//|      elle est mise a 1 sinon, sauf si E... ou e... est precise:
//|      Eaa.b : si |val|>=1 erreur = aa.b*sqrt(|val|)
//|              si |val|<1  erreur = aa.b
//|              si aa.b <=0 alors aa.b=1.0
//|              E seul est equivalent a E1.0
//|      eaa.b : erreur = aa.b
//|              si aa.b <=0 alors aa.b=1.0
//|              e seul est equivalent a e1.0
//|      xaa.b : demande de centrage: on fit x-aa.b au lieu de x)
//|      x : demande de centrage: on fit x-xc au lieu de x
//|          avec xc=abscisse du milieu de l'histogramme
//|          Actif pour  exp+poly 1D, poly 1D
//|          Pour gauss+poly 1D, xc est le centre de la gaussienne.
//|2D    yaa.b et y : idem ''xaa.b et x'' mais pour y
//|
////////////////////////////////////////////////////////////////////////
//| --------------- Def des fcts de fit par l'utilisateur ---------------
//| Pour definir une fonction parametree utilisateur:
//| Un exemple de fichier contenant des fonctions a fitter: PIext/userfitfunex.c
//| > crfitfil nom_fichier nom_fonction nvar npar
//|     pour creer un squelette de fichier C contenant la fonction a fitter
//|            nom_fichier : nom du fichier ou est le code de la fonction
//|                          c'est un fichier C de nom blabla.c
//|            nom_fonction : nom de la fonction ex: gauss2
//|            nvar : nombre de variables x[] (1D=1, 2D=2, etc...)
//|            npar : nombre de parametre p[]
//|   - La fonction est donc F(x[];p[]), soit par ex:
//|            double gauss2(double const* x,double const* p)
//|   - L'utilisateur a la possibilite de donner egalement la fonction
//|     retournant les derivees de la fonction par rapport aux parametres
//|            double gauss2_der(double const* x,double const* p,double* dp)
//|     ou dp[i] = dF(x[];p[])/dp[i]
//|     Par convention, cette fonction a le meme nom suivi de ''_der''
//| > crfitfun nom_fichier nom_fonction nvar npar
//|     Permet de compiler le fichier ''nom_fichier'' et de linker
//|     la fonction ''nom_fonction''
//|   - Un meme fichier peut contenir plusieurs fonctions a fitter
//|
////////////////////////////////////////////////////////////////////////

/* --Methode-- */
PIAFitter::PIAFitter(PIACmd *piac, PIStdImgApp* app)
: mApp(app), FWindFit(NULL)
, mNObj(""), mObj(NULL), mGData(NULL)
, mNPar(0), mNVar(0),  mNBinX(0), mNBinY(0), mNData(0)
, mPar(1), mStep(1), mMin(1), mMax(1), mFix(1)
, mParSave(1), mStepSave(1), mMinSave(1), mMaxSave(1), mFixSave(1)
, mFit(NULL)
, mV(NULL), mH(NULL), mM(NULL), mH2(NULL), mIm(NULL), mG(NULL)
, mFunction(NULL), mFName(""), mFunc(NULL)
, mDlUFunc(NULL), mNameFitFunc(""), mFitFunc(NULL), mFitFuncDer(NULL)
, mUFNVar(0), mUFNPar(0)
{
// Init et Set des options
mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
mParSave=0.; mStepSave=1.; mMinSave=1.; mMaxSave=-1.; mFixSave=0.;
ResetOptions();
mOptSave = mOpt;

// enregistrement des ordres de fit
string kw, usage, grp = "Fitting";

kw = "helpfit";
usage="";
usage+="\n//////////////////////////////////////////////////////////////////////    ";
usage+="\n| --------------- Fit Lineaire a 1 et 2 dimensions ---------------    ";
usage+="\n| Syntaxe:    ";
usage+="\n| fitlin nom pnn [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]    ";
usage+="\n| avec:    ";
usage+="\n| nom  : cf commentaire ordre ''fit''    ";
usage+="\n| pnn  : fit polynome degre nn avec classe Poly (lineaire) 1D ou 2D    ";
usage+="\n| o:aa,...,bb : cf commentaires ordre ''fit''    ";
usage+="\n|    ";
usage+="\n//////////////////////////////////////////////////////////////////////    ";
usage+="\n| --------------- Fit d'objets a 1 et 2 dimensions ---------------    ";
usage+="\n|                 avec interface d'aide graphique    ";
usage+="\n| Syntaxe:    ";
usage+="\n| fitw nom func    ";
usage+="\n|     [p:p1,p2,...,pn s:s1,s2,...,sn m:m1,m2,...,mn M:M1,mM2,...,Mn f:f1,...,fn]    ";
usage+="\n|     [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]    ";
usage+="\n| cf commentaire ordre ''fit''    ";
usage+="\n|    ";
usage+="\n//////////////////////////////////////////////////////////////////////    ";
usage+="\n| --------------- Fit d'objets a 1 et 2 dimensions ---------------    ";
usage+="\n| Syntaxe:    ";
usage+="\n| fit nom func    ";
usage+="\n|     [p:p1,p2,...,pn s:s1,s2,...,sn m:m1,m2,...,mn M:M1,mM2,...,Mn f:f1,...,fn]    ";
usage+="\n|     [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]    ";
usage+="\n|----- OBJET -----    ";
usage+="\n| nom  : nom de l'objet qui peut etre:    ";
usage+="\n|        fit-1D:  Vector,Histo1D,HProf ou GeneraFitData(1D)    ";
usage+="\n|        fit-2D:  Matrix,Histo2D,Image<T> ou GeneraFitData(2D)    ";
usage+="\n|    ";
usage+="\n|----- FUNCTION -----    ";
usage+="\n| func : pnn : fit polynome degre nn avec GeneralFit (non-lineaire) 1D ou 2D    ";
usage+="\n|      : gnn : fit gaussienne (hauteur) + polynome de degre nn 1D    ";
usage+="\n|      : g   : fit gaussienne (hauteur) 1D    ";
usage+="\n|      : enn : fit exponentielle + polynome de degre nn 1D    ";
usage+="\n|      : e   : fit exponentielle 1D    ";
usage+="\n|      : Gnn : fit gaussienne (volume) + polynome de degre nn 1D    ";
usage+="\n|      : G   : fit gaussienne (volume) 1D    ";
usage+="\n|      :     : fit gaussienne+fond (volume) 2D    ";
usage+="\n|      : Gi  : fit gaussienne+fond integree (volume) 2D    ";
usage+="\n|      : d   : fit DL de gaussienne+fond (volume) 2D    ";
usage+="\n|      : di  : fit DL de gaussienne+fond integree (volume) 2D    ";
usage+="\n|      : D   : fit DL de gaussienne+fond avec coeff variable p6 (volume) 2D    ";
usage+="\n|      : Di  : fit DL de gaussienne+fond integree avec coeff variable p6 (volume) 2D    ";
usage+="\n|      : M   : fit Moffat+fond (expos=p6) (volume) 2D    ";
usage+="\n|      : Mi  : fit Moffat+fond integree (expos=p6) (volume) 2D    ";
usage+="\n|      : Autre : fonction definie par l'utilisateur (cf commentaires ci apres)    ";
usage+="\n|    ";
usage+="\n|----- INIT PARAMETRES ET ETAT DU FIT -----    ";
usage+="\n| p    : p1,...,pn : valeur d'initialisation des parametres (def=0)    ";
usage+="\n| s    : s1,...,sn : valeur des steps de depart (def=1)    ";
usage+="\n| m    : m1,...,mn : valeur des minima (def=1)    ";
usage+="\n| M    : M1,...,Mn : valeur des maxima (def=-1) (max<=min : pas de limite)    ";
usage+="\n| f    : f1,...,fn : si >=1 parametre fixe sinon libre (def=0)    ";
usage+="\n| - Remarque: si pi,si,mi ou Mi = '!' la valeur correspondante n'est pas changee    ";
usage+="\n|    ";
usage+="\n|----- OPTIONS -----    ";
usage+="\n| o    : options ''o:Eaa.b,eaa.b,f,r,caa.b,Xaa.b''    ";
usage+="\n|      F : initialisation a partir des resultats et de l'etat du fit precedent    ";
usage+="\n|          (option non prioritaire sur les definitions p:,s:,m:,M:,f:o:)    ";
usage+="\n|      a : initialisation automatique des parametres du fit (si disponible)    ";
usage+="\n|          (option non prioritaire sur les definitions p:,s:,m:,M:,f:o:    ";
usage+="\n|           et sur l'initialisation pa l'option ''F'')    ";
usage+="\n|      f : generation d'un Objet identique contenant la fonction fittee    ";
usage+="\n|      r : generation d'un Objet identique contenant les residus    ";
usage+="\n|      Xaa.b : aa.b valeur du DXi2 d'arret (def=1.e-3)    ";
usage+="\n|      Naa : aa nombre maximum d'iterations (def=100)    ";
usage+="\n|      la : niveau ''a'' de print: a=niveau de print Fit1/2D    ";
usage+="\n|      da : niveau ''a'' de debug: a=niveau de GeneralFit    ";
usage+="\n|      Ii1/i2 : numeros des bins X de l'histos utilises pour le fit [i1,i2]    ";
usage+="\n|2D    Jj1/j2 : numeros des bins Y de l'histos utilises pour le fit [j1,j2]    ";
usage+="\n|    - L'erreur est celle associee a l'objet (si elle existe),    ";
usage+="\n|      elle est mise a 1 sinon, sauf si E... ou e... est precise:    ";
usage+="\n|      Eaa.b : si |val|>=1 erreur = aa.b*sqrt(|val|)    ";
usage+="\n|              si |val|<1  erreur = aa.b    ";
usage+="\n|              si aa.b <=0 alors aa.b=1.0    ";
usage+="\n|              E seul est equivalent a E1.0    ";
usage+="\n|      eaa.b : erreur = aa.b    ";
usage+="\n|              si aa.b <=0 alors aa.b=1.0    ";
usage+="\n|              e seul est equivalent a e1.0    ";
usage+="\n|      xaa.b : demande de centrage: on fit x-aa.b au lieu de x)    ";
usage+="\n|      x : demande de centrage: on fit x-xc au lieu de x    ";
usage+="\n|          avec xc=abscisse du milieu de l'histogramme    ";
usage+="\n|          Actif pour  exp+poly 1D, poly 1D    ";
usage+="\n|          Pour gauss+poly 1D, xc est le centre de la gaussienne.    ";
usage+="\n|2D    yaa.b et y : idem ''xaa.b et x'' mais pour y    ";
usage+="\n|    ";
usage+="\n//////////////////////////////////////////////////////////////////////    ";
usage+="\n| --------------- Def des fcts de fit par l'utilisateur ---------------    ";
usage+="\n| Pour definir une fonction parametree utilisateur:    ";
usage+="\n| Un exemple de fichier contenant des fonctions a fitter: PIext/userfitfunex.c    ";
usage+="\n| > crfitfil nom_fichier nom_fonction nvar npar    ";
usage+="\n|     pour creer un squelette de fichier C contenant la fonction a fitter    ";
usage+="\n|            nom_fichier : nom du fichier ou est le code de la fonction    ";
usage+="\n|                          c'est un fichier C de nom blabla.c    ";
usage+="\n|            nom_fonction : nom de la fonction ex: gauss2    ";
usage+="\n|            nvar : nombre de variables x[] (1D=1, 2D=2, etc...)    ";
usage+="\n|            npar : nombre de parametre p[]    ";
usage+="\n|   - La fonction est donc F(x[];p[]), soit par ex:    ";
usage+="\n|            double gauss2(double const* x,double const* p)    ";
usage+="\n|   - L'utilisateur a la possibilite de donner egalement la fonction    ";
usage+="\n|     retournant les derivees de la fonction par rapport aux parametres    ";
usage+="\n|            double gauss2_der(double const* x,double const* p,double* dp)    ";
usage+="\n|     ou dp[i] = dF(x[];p[])/dp[i]    ";
usage+="\n|     Par convention, cette fonction a le meme nom suivi de ''_der''    ";
usage+="\n| > crfitfun nom_fichier nom_fonction nvar npar    ";
usage+="\n|     Permet de compiler le fichier ''nom_fichier'' et de linker    ";
usage+="\n|     la fonction ''nom_fonction''    ";
usage+="\n|   - Un meme fichier peut contenir plusieurs fonctions a fitter    ";
usage+="\n|    ";
usage+="\n//////////////////////////////////////////////////////////////////////    ";
piac->RegisterHelp(kw,usage,grp);

kw = "fit";
usage = "Fitting function to DataObjects (Histo, Histo2D, Vector, ...)";
usage += "\n Usage: fit nomobj func [Options]";
usage += "\n [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]";
usage += "\n  Related commands: fitw fitlin crfitfun crfitfil";  
piac->RegisterCommand(kw,usage,this,grp);

kw = "fitw";
usage = "Fitting function to DataObjects with Interactive Window Help";
usage += "\n Usage: fitw nomobj func [Options]";
usage += "\n [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]";
usage += "\n  Related commands: fit fitlin crfitfun crfitfil";  
piac->RegisterCommand(kw,usage,this,grp);

kw = "fitlin";
usage = "Linear Fitting of Polynoms to DataObjects";
usage += "\n Usage: fitlin nomobj func [o:... o:...]";
usage += "\n  Related commands: fit fitw crfitfun crfitfil";  
piac->RegisterCommand(kw,usage,this,grp);

kw = "crfitfun";
usage = "Creation et link de function utilisateur pour le fit";
usage += "\n Usage: crfitfun file func nvar npar";
usage += "\n  Related commands: fit fitw fitlin crfitfil";  
piac->RegisterCommand(kw,usage,this,grp);

kw = "crfitfil";
usage = "Creation de fichier C pour function utilisateur du fit";
usage += "\n Usage: crfitfil file func nvar npar";
usage += "\n  Related commands: fit fitw fitlin crfitfun";  
piac->RegisterCommand(kw,usage,this,grp);
}

/* --Methode-- */
PIAFitter::~PIAFitter()
{
if(mFit != NULL)      {delete mFit;      mFit=NULL;}
if(mFunction != NULL) {delete mFunction; mFunction=NULL;}
if(mFunc != NULL)     {delete mFunc;     mFunc=NULL;}
if(mGData != NULL)    {delete mGData;    mGData=NULL;}
if(FWindFit != NULL)  {delete FWindFit;  FWindFit=NULL;}
if (mDlUFunc != NULL) {delete mDlUFunc;  mDlUFunc=NULL;}
}

/* --Methode-- */
void PIAFitter::ResetDPointer(void)
// Reset des pointeurs sur les donnees
{
mNObj = ""; mObj=NULL;
mV=NULL; mH=NULL; mM=NULL; mH2=NULL; mIm=NULL; mG=NULL;
}

/* --Methode-- */
void PIAFitter::ResetOptions(void)
// Reset des options
{
mOpt.okres = mOpt.okfun = false;
mOpt.polcx = mOpt.polcy = 0;
mOpt.xc = mOpt.yc = 0.;
mOpt.stc2 = 1.e-3;
mOpt.nstep = 100;
mOpt.err_e = mOpt.err_E = -1.;
mOpt.lp = 1; mOpt.lpg = 0;
mOpt.i1 = mOpt.j1 = mOpt.i2 = mOpt.j2 = -1;
mOpt.fromlastfit = false;
mOpt.autoinifit = false;
}

/* --Methode-- */
void PIAFitter::DecodeOptions(string opt)
// Decodage des options a partir des chaines de caracteres
{
ResetOptions();
if(opt.length()<=0) return;

size_t p,q; string dum;
opt = "," + opt + ",";

if(strstr(opt.c_str(),",F,"))  // init options/param a partir du fit precedent
  {mOpt = mOptSave; mOpt.fromlastfit = true;}
if(strstr(opt.c_str(),",a,")) mOpt.autoinifit = true; // init auto par fit
if(strstr(opt.c_str(),",r,")) mOpt.okres = true;  // residus
if(strstr(opt.c_str(),",f,")) mOpt.okfun = true;  // fonction fittee
if(strstr(opt.c_str(),",x")) { // Demande de centrage (fit X=x-xc)
  mOpt.polcx = 2; // Le centrage est calcule automatiquement
  p = opt.find(",x"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) {
    sscanf(dum.c_str(),",x%lf",&mOpt.xc);
    mOpt.polcx = 1; // Le centrage est fixe par la valeur lue
  }
}
if(strstr(opt.c_str(),",y")) { // Demande de centrage (fit Y=y-yc)
  mOpt.polcy = 2; // Le centrage est calcule automatiquement
  p = opt.find(",y"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) {
    sscanf(dum.c_str(),",y%lf",&mOpt.yc);
    mOpt.polcy = 1; // Le centrage est fixe par la valeur lue
  }
}
if(strstr(opt.c_str(),",E")) { // Erreurs imposees a "sqrt(val)" ou "aa.b*sqrt(val)"
  p = opt.find(",E"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",E%lf",&mOpt.err_E);
  if(mOpt.err_E<=0.) mOpt.err_E = 1.;  //mOpt.err_e=-1.;
}
if(strstr(opt.c_str(),",e")) { // Erreurs imposees a "1" ou "aa.b"
  p = opt.find(",e"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",e%lf",&mOpt.err_e);
  if(mOpt.err_e<=0.) mOpt.err_e = 1.;  //mOpt.err_E=-1.;
}
if(strstr(opt.c_str(),",X")) { // Valeur du StopChi2
  p = opt.find(",X"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",X%lf",&mOpt.stc2);
}
if(strstr(opt.c_str(),",N")) { // Nombre maximum d'iterations
  p = opt.find(",N"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",N%d",&mOpt.nstep);
}
if(strstr(opt.c_str(),",l")) { // niveau de print
  p = opt.find(",l"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",l%d",&mOpt.lp);
}
if(strstr(opt.c_str(),",d")) { // niveau de debug du fit
  p = opt.find(",d"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",d%d",&mOpt.lpg);
}
if(strstr(opt.c_str(),",I")) { // intervalle de fit selon X
  p = opt.find(",I"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",I%d/%d",&mOpt.i1,&mOpt.i2);
}
if(strstr(opt.c_str(),",J")) { // intervalle de fit selon Y
  p = opt.find(",J"); q = opt.find_first_of(',',p+1);
  dum = opt.substr(p,q-p);
  if(dum.length()>2) sscanf(dum.c_str(),",J%d/%d",&mOpt.j1,&mOpt.j2);
}

return;
}

/* --Methode-- */
void PIAFitter::DecodeObject(string obj)
// Decodage de l'objet a fitter
{
NamedObjMgr omg;
AnyDataObj* mobj = omg.GetObj(obj);
if(mobj == NULL)
  {cout<<"PIAFitter::DecodeObject Error , Pas d'objet de nom "<<obj<<endl;
   return;}

// Reset
ResetDPointer();
mNVar = mNBinX = mNBinY = mNData = 0;

mObj = mobj; mNObj = obj;
string ctyp = typeid(*mObj).name();

// 1D
if (typeid(*mObj) == typeid(Vector)) {
  mNVar = 1;
  mV = (Vector*) mObj; mNBinX = mV->NElts(); mNBinY = 1;
  }
else if ( (typeid(*mObj) == typeid(HProf)) || (typeid(*mObj) == typeid(Histo)) ) {
  mNVar = 1;
  mH = (Histo*)  mObj; mNBinX = mH->NBins(); mNBinY = 1;
  }
// 2D
else if (typeid(*mObj) == typeid(Matrix)) {
  mNVar = 2;
  mM = (Matrix*) mObj; mNBinX = mM->NCol(); mNBinY =mM->NRows();
  }
else if (typeid(*mObj) == typeid(Histo2D)) {
  mNVar = 2;
  mH2 = (Histo2D*) mObj; mNBinX = mH2->NBinX(); mNBinY = mH2->NBinY();
  }
else if (typeid(*mObj) == typeid(GeneralFitData)) {
  mG = (GeneralFitData*) mObj; mNBinX = mG->NData(); mNBinY = 1;
  if(     mG->NVar()==1) mNVar = 1;
  else if(mG->NVar()==2) mNVar = 2;
  else
     {cout<<"PIAFitter::DecodeObject Error: \n"
         <<" GeneralFitData ne peut avoir que 1 ou 2 variables d'abscisse: "
         <<((GeneralFitData*) mObj)->NVar()<<endl;
     return;}
  }
#ifdef SANS_EVOLPLANCK
else if (dynamic_cast<RzImage*>(mObj)) {
  mNVar = 2;
  mIm = (RzImage*) mObj; mNBinX = mIm->XSize(); mNBinY = mIm->YSize();
  }
#else
//CMV_A_FAIRE else if (dynamic_cast<RzImage*>(mObj)) {
//CMV_A_FAIRE   mNVar = 2;
//CMV_A_FAIRE   mIm = (RzImage*) mObj; mNBinX = mIm->XSize(); mNBinY = mIm->YSize();
//CMV_A_FAIRE   }
#endif
else {
  cout<<"PIAFitter::DecodeObject Error , Objet n'est pas un "
      <<"Histo1D/HProf/Vector/Histo2D/Image/Matrix/GeneralFitData "<<ctyp<<endl;
  return;
  }

mNData = mNBinX*mNBinY;
if(mNData<=0) {
  cout<<"L'objet a "<<mNBinX<<","<<mNBinY<<" bins ("<<mNData<<")"<<endl;
  return;
}

return;
}

/* --Methode-- */
void PIAFitter::CheckOptions(void)
// Check consistence des choix des options
// ATTENTION: cette methode a besoin que l'objet a fitter soit deja lu
{
if(mOpt.lp<0) mOpt.lp = 0;
if(mOpt.lpg<0) mOpt.lpg = 0;
mOpt.i1 = (mOpt.i1<0||mOpt.i1>=mNBinX)? 0: mOpt.i1;
mOpt.i2 = (mOpt.i2<0||mOpt.i2>=mNBinX||mOpt.i2<=mOpt.i1)? mNBinX-1: mOpt.i2;
if(mNVar>=2) {
  mOpt.j1 = (mOpt.j1<0||mOpt.j1>=mNBinY)? 0: mOpt.j1;
  mOpt.j2 = (mOpt.j2<0||mOpt.j2>=mNBinY||mOpt.j2<=mOpt.j1)? mNBinY-1: mOpt.j2;
} else mOpt.j2 = mOpt.j1 = 0;
if(mOpt.polcx==2) {
  if(mV||mM)   mOpt.xc = (mOpt.i2-mOpt.i1+1)/2.;
  else if(mH)  mOpt.xc = (mH->XMin()+mH->XMax())/2.;
  else if(mH2) mOpt.xc = (mH2->XMin()+mH2->XMax())/2.;
  else if(mG)  {double mini,maxi; mG->GetMinMax(2,mini,maxi); mOpt.xc=(mini+maxi)/2.;}
  else if(mIm) {mOpt.xc = mIm->XOrg() * mIm->XPxSize()*(mOpt.i2-mOpt.i1+1)/2.;}
}
if(mOpt.polcy==2 && mNVar>=2) {
  if(mM)  mOpt.yc = (mOpt.j2-mOpt.j1+1)/2.;
  if(mH2) mOpt.yc = (mH2->YMin()+mH2->YMax())/2.;
  if(mG)  {double mini,maxi; mG->GetMinMax(12,mini,maxi); mOpt.yc=(mini+maxi)/2.;}
  if(mIm) {mOpt.yc = mIm->YOrg() * mIm->YPxSize()*(mOpt.j2-mOpt.j1+1)/2.;}
}
if(mOpt.err_e>0.) mOpt.err_E=-1.;
if(mOpt.err_E>0.) mOpt.err_e =-1.;
if(mOpt.stc2<=0.) mOpt.stc2 = 1.e-3;
if(mOpt.nstep<2) mOpt.nstep = 100;

return;
}

/* --Methode-- */
void PIAFitter::PrintOptions(void)
// Print des options
{
cout<<"Fit["<<mNBinX<<","<<mNBinY<<"] ("<<mNData<<") dim="<<mNVar<<":"
    <<" Int=["<<mOpt.i1<<","<<mOpt.i2<<"],["<<mOpt.j1<<","<<mOpt.j2<<"]"<<endl
    <<" Cent="<<mOpt.polcx<<","<<mOpt.polcy<<","<<mOpt.xc<<"+x"<<","<<mOpt.yc<<"+y"
    <<" TypE="<<mOpt.err_e<<","<<mOpt.err_E
    <<" StpX2="<<mOpt.stc2<<" Nstep="<<mOpt.nstep<<endl
    <<" Init.LFit="<<mOpt.fromlastfit<<" AutoIni.Fit="<<mOpt.autoinifit
    <<" lp,lpg="<<mOpt.lp<<","<<mOpt.lpg<<endl;
return;
}

/* --Methode-- */
void PIAFitter::FillGData(void)
// Fill generalfitdata pour le fit
// ATTENTION: cette methode a besoin que les options soient decodees et checkees
{
if(mNData<=0) return;
if(mGData!=NULL) {delete mGData; mGData=NULL;}
mGData = new GeneralFitData(mNVar,mNData,0);

for(int i=mOpt.i1;i<=mOpt.i2;i++) for(int j=mOpt.j1;j<=mOpt.j2;j++) {
  double x,y=0.,f,e;

  if(mV)
    {x= (double) i; f=(*mV)(i); e=1.;}
  else if(mH)
    {x=mH->BinCenter(i); f=(*mH)(i); e=(mH->HasErrors())?mH->Error(i):1.;}
  else if(mM)
    {x=(double) i; y=(double) j; f=(*mM)(j,i); e=1.;}
  else if(mH2)
    {float xf,yf; mH2->BinCenter(i,j,xf,yf); x=(double)xf; y=(double)yf;
     f=(*mH2)(i,j); e=(mH2->HasErrors())?mH2->Error(i,j):1.;}
  else if(mIm)
    {x=mIm->XOrg()+(i+0.5)*mIm->XPxSize(); y=mIm->YOrg()+(j+0.5)*mIm->YPxSize();
     f=mIm->DValue(i,j); e=1.;}
  else if(mG&&mNVar==1) {x= mG->X(i); f=mG->Val(i); e=mG->EVal(i);}
  else if(mG&&mNVar==2) {x= mG->X(i); y=mG->Y(i); f=mG->Val(i); e=mG->EVal(i);}
  else x=y=f=e=0.;

  // Gestion des erreurs a utiliser
  if(mOpt.err_e>0.) e=mOpt.err_e;
  else if(mOpt.err_E>0.) {e=(f<-1.||f>1.)?mOpt.err_E*sqrt(fabs(f)):mOpt.err_E;}

  // Remplissage de generalfitdata
  if(mFName[0]=='p') {x -= mOpt.xc; if(mNVar>=2) y -= mOpt.yc;}
  if(mNVar==1)      mGData->AddData1(x,f,e);
  else if(mNVar==2) mGData->AddData2(x,y,f,e);
}

if(mGData->NData()<=0) {
  cout<<"Pas de donnees dans GeneralFitData: "<<mGData->NData()<<endl;
  return;
}

if(mOpt.lpg>1) {
  mGData->PrintStatus();
  mGData->PrintData(0);
  mGData->PrintData(mGData->NData()-1);
}

return;
}

/* --Methode-- */
void PIAFitter::DecodeFunction(string func)
// Fonction a fitter
// ATTENTION: cette methode a besoin que les donnees soient lues
{
if(func.length()<=0)
  {cout<<"PIAFitter::DecodeFunction Donnez un nom de fonction a fitter."<<endl;
   return;}
if(mFunction!=NULL) {delete mFunction; mFunction=NULL;}
if(mFunc    !=NULL) {delete mFunc;     mFunc=NULL;}
mNPar=0; mFName = "";

if(func == mNameFitFunc) { // Fonction utilisateur
  mFunc = new GeneralFunc(mUFNVar,mUFNPar,mFitFunc,mFitFuncDer);
  mNPar = mUFNPar;
  cout<<"Fonction utilisateur nvar="<<mUFNVar<<", npar="<<mNPar<<endl;
} else if(func[0]=='p' && mNVar==1) { //polynome
  int degre=0; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
  cout<<"Fit polynome 1D de degre "<<degre<<endl;
  Polyn1D* myf = new Polyn1D(degre,mOpt.xc);
  mFunction = myf; mNPar = mFunction->NPar(); mFName = func;

} else if(func[0]=='e' && mNVar==1) { //exponentielle
  int degre=-1; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
  cout<<"Fit d'exponentielle+polynome 1D de degre "<<degre<<endl;
  Exp1DPol* myf;
  if(degre>=0) myf = new Exp1DPol((unsigned int)degre,mOpt.xc);
       else    myf = new Exp1DPol(mOpt.xc);
  mFunction = myf; mNPar = mFunction->NPar(); mFName = func;

} else if(func[0]=='g' && mNVar==1) { //gaussienne en hauteur
  int degre=-1; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
  cout<<"Fit de Gaussienne_en_hauteur+polynome 1D de degre "<<degre<<endl;
  Gauss1DPol* myf;
  if(degre>=0) myf = new Gauss1DPol((unsigned int)degre,((mOpt.polcx)?true:false));
  else { bool bfg = (mOpt.polcx)?true:false;   myf = new Gauss1DPol(bfg); }
  mFunction = myf; mNPar = mFunction->NPar(); mFName = func;

} else if(func[0]=='G' && mNVar==1) { //gaussienne en volume
  int degre=-1; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
  cout<<"Fit de Gaussienne_en_volume+polynome 1D de degre "<<degre<<endl;
  GaussN1DPol* myf;
  if(degre>=0) myf = new GaussN1DPol((unsigned int)degre,((mOpt.polcx)?true:false));
  else  { bool bfg = (mOpt.polcx)?true:false;   myf = new GaussN1DPol(bfg); }
  mFunction = myf; mNPar = mFunction->NPar(); mFName = func;

} else if(func[0]=='p' && mNVar==2) { //polynome 2D
  int degre=0; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
  cout<<"Fit polynome 2D de degre "<<degre<<endl;
  Polyn2D* myf = new Polyn2D(degre,mOpt.xc,mOpt.yc);
  mFunction = myf; mNPar = mFunction->NPar(); mFName = func;

} else if(func[0]=='G' && mNVar==2) { //gaussienne+fond en volume
  int integ=0; if(func.length()>1) if(func[1]=='i') integ=1;
  cout<<"Fit de Gaussienne+Fond 2D integ="<<integ<<endl;
  if(integ) {GauRhInt2D* myf = new GauRhInt2D; mFunction = myf;}
    else    {GauRho2D*   myf = new GauRho2D;   mFunction = myf;}
  mNPar = mFunction->NPar(); mFName = func;

} else if(func[0]=='d' && mNVar==2) { //DL gaussienne+fond en volume
  int integ=0; if(func.length()>1) if(func[1]=='i') integ=1;
  cout<<"Fit de DL de Gaussienne+Fond 2D integ="<<integ<<endl;
  if(integ) {GdlRhInt2D* myf = new GdlRhInt2D; mFunction = myf;}
    else    {GdlRho2D*   myf = new GdlRho2D;   mFunction = myf;}
  mNPar = mFunction->NPar(); mFName = func;

} else if(func[0]=='D' && mNVar==2) { //DL gaussienne+fond avec coeff variable p6 en volume
  int integ=0; if(func.length()>1) if(func[1]=='i') integ=1;
  cout<<"Fit de DL de Gaussienne+Fond avec coeff variable (p6) 2D integ="<<integ<<endl;
  if(integ) {Gdl1RhInt2D* myf = new Gdl1RhInt2D; mFunction = myf;}
    else    {Gdl1Rho2D*   myf = new Gdl1Rho2D;   mFunction = myf;}
  mNPar = mFunction->NPar(); mFName = func;

} else if(func[0]=='M' && mNVar==2) { //Moffat+fond (volume)
  int integ=0; if(func.length()>1) if(func[1]=='i') integ=1;
  cout<<"Fit de Moffat+Fond (expos=p6) 2D integ="<<integ<<endl;
  if(integ) {MofRhInt2D* myf = new MofRhInt2D; mFunction = myf;}
    else    {MofRho2D*   myf = new MofRho2D;   mFunction = myf;}
   mNPar = mFunction->NPar(); mFName = func;

} else {
  cout<<"Fonction "<<func<<" inconnue pour la dim "<<mNVar<<endl;
  return;
}

return;
}

/* --Methode-- */
void PIAFitter::ReSetParam(void)
// Pour (re)dimensionner les tableaux des parametres
// ATTENTION: cette methode a besoin que la fonction soit connue
{
if(mNPar==0) return;
mPar.Realloc(mNPar);  mPar=0.;
mStep.Realloc(mNPar); mStep=1.;
mMin.Realloc(mNPar);  mMin=1.;
mMax.Realloc(mNPar);  mMax=-1.;
mFix.Realloc(mNPar);  mFix=0.;
return;
}

/* --Methode-- */
void PIAFitter::AutoIniFit(void)
// Initialisation automatique des parametres du fit
{
if(mFunc) { // Fonction utilisateur
  cout<<"PIAFitter::AutoIniFit : Fonction utilisateur, pas d'init auto"<<endl;
} else if(mFName[0]=='p' && mNVar==1) { //poly 1D
  cout<<"PIAFitter::AutoIniFit : Auto Init type polynome 1D"<<endl;
  IniFitP1ou2D();

} else if(mFName[0]=='g' && mNVar==1) { //gauss en haut 1D
  cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_haut+pol 1D"<<endl;
  IniFitGhP1D();

} else if(mFName[0]=='G' && mNVar==1) { //gauss en vol 1D
  cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_vol+pol 1D"<<endl;
  IniFitGvP1D();

} else if(mFName[0]=='p' && mNVar==2) { //poly 2D
  cout<<"PIAFitter::AutoIniFit : Auto Init type polynome 2D"<<endl;
  IniFitP1ou2D();

} else if((mFName[0]=='G'||mFName[0]=='d') && mNVar==2) { //gauss+fond en vol
  cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_vol+fond 2D"<<endl;
  IniFitGv2D();

} else if(mFName[0]=='D' && mNVar==2) { //DLgauss+fond en vol 2D
  cout<<"PIAFitter::AutoIniFit : Auto Init type DLgauss_vol+fond 2D"<<endl;
  IniFitGv2D();
  Vector dum(mNPar);
  dum = mPar;  mPar(6)=0.5;   mPar(7)=dum(6);
  dum = mStep; mStep(6)=0.05; mStep(7)=dum(6);
  dum = mMin;  mMin(6)=0.;    mMin(7)=dum(6);
  dum = mMax;  mMax(6)=10.;   mMax(7)=dum(6);
  dum = mFix;  mFix(6)=1.1;   mFix(7)=dum(6);

} else if(mFName[0]=='M' && mNVar==2) { //Moffat+fond en vol 2D
  cout<<"PIAFitter::AutoIniFit : Auto Init type moffat_vol+fond 2D"<<endl;
  IniFitGv2D();
  Vector dum(mNPar);
  dum = mPar;  mPar(6)=2.5;   mPar(7)=dum(6);
  dum = mStep; mStep(6)=0.05; mStep(7)=dum(6);
  dum = mMin;  mMin(6)=1.2;   mMin(7)=dum(6);
  dum = mMax;  mMax(6)=10.;   mMax(7)=dum(6);
  dum = mFix;  mFix(6)=1.1;   mFix(7)=dum(6);

} else {
  cout<<"PIAFitter::AutoIniFit : Fonction "<<mFName
      <<", \n  pas d'init auto en dim "<<mNVar<<endl;
  return;
}

return;
}

/* --Methode-- */
void PIAFitter::InitParFromLastFit(void)
// Initialisation du fit a partir des resultats du fit precedent
// (Les arguments de la ligne de commande sont prioritaires)
// ATTENTION: cette methode a besoin que la fonction soit connue
{
if(mNPar==0) return;
int n = mParSave.NElts();
if(n==0) return;

for(int i=0;i<mNPar;i++) {
  if(i>=n) break;
  mPar(i)  = mParSave(i);
  mStep(i) = mStepSave(i);
  mMin(i)  = mMinSave(i);
  mMax(i)  = mMaxSave(i);
  mFix(i)  = mFixSave(i);
}

return;
}

/* --Methode-- */
void PIAFitter::DecodeParam(string par,string step,string min,string max,string fix)
// Decodage des parametres
// ATTENTION: cette methode a besoin que la fonction soit connue
//                            et que les options soient decodees
{
if(mNPar==0) return;
ReSetParam();
if(mOpt.autoinifit)  AutoIniFit();
if(mOpt.fromlastfit) InitParFromLastFit();
Vector* v=NULL; string* s=NULL;
for(int j=0;j<5;j++) {
  if(j==0) {v=&mPar; s=&par;} else if(j==1) {v=&mStep; s=&step;}
  else if(j==2) {v=&mMin; s=&min;} else if(j==3) {v=&mMax; s=&max;}
  else if(j==4) {v=&mFix; s=&fix;}
  if(s->length()>0) *s += ",";
  for(int i=0;i<mNPar;i++) {
    if(s->length()<=0) break;
    if((*s)[0]!='!') sscanf(s->c_str(),"%lf",&(*v)(i));
    size_t p = s->find_first_of(',') + 1;
    if(p>=s->length()) *s = ""; else *s = s->substr(p);
  }
}
return;
}

/* --Methode-- */
int PIAFitter::DoFit(void)
// Fit avec generalfit
{
if(mNPar==0) return -1000;
if((mFunction==NULL && mFunc==NULL) || mGData==NULL) {
  cout<<"PIAFitter::DoFit error: mFunction,mFunc="<<mFunction<<","<<mFunc
      <<" mGData="<<mGData<<endl;
  return -1001;
}
// classe GeneraFit
if(mFit != NULL) {delete mFit; mFit= NULL;}
if(mFunction)    mFit = new GeneralFit(mFunction);
  else if(mFunc) mFit = new GeneralFit(mFunc);
if(!mFit) return(-1002);

// Options et Set de GeneralFit
mFit->SetDebug(mOpt.lpg);
mFit->SetData(mGData);
mFit->SetStopChi2(mOpt.stc2);
mFit->SetMaxStep(mOpt.nstep);
{for(int i=0;i<mNPar;i++) {
  char str[10];
  sprintf(str,"P%d",i);
  mFit->SetParam(i,str,mPar(i),mStep(i),mMin(i),mMax(i));
  if((int)(mFix(i)+0.001)>=1.) mFit->SetFix(i);
}}
if(mOpt.lp>1) mFit->PrintFit();

// Fit
mParSave = mPar;
mStepSave = mStep; mMinSave = mMin; mMaxSave = mMax; mFixSave = mFix;
mOptSave = mOpt;
double c2r = -1.;
int rcfit = mFit->Fit();
if(mOpt.lp>0) mFit->PrintFit();
if(rcfit>0) {
  c2r = mFit->GetChi2Red();
  cout<<"C2r_Reduit = "<<c2r<<" nstep="<<mFit->GetNStep()<<" rc="<<rcfit<<endl;
  for(int i=0;i<mNPar;i++) mParSave(i)=mFit->GetParm(i);
} else {
  cout<<"echec Fit, rc = "<<rcfit<<"  nstep="<<mFit->GetNStep()<<endl;
  mFit->PrintFitErr(rcfit);
}

return rcfit;
}

/* --Methode-- */
void PIAFitter::FitFunRes(void)
// Mise a disposition des resultats (fonction fitee et residus)
{
if(mFit && (mOpt.okres||mOpt.okfun)) {
  NamedObjMgr omg;
  string nomres = mNObj + "_res";
  string nomfun = mNObj + "_fun";
//####################
#ifdef SANS_EVOLPLANCK
//####################
  if(mV) {
    if(mOpt.okres) {Vector* ob = mV->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
    if(mOpt.okfun) {Vector* ob = mV->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
  } else if(mH) {
    if(mOpt.okres) {Histo* ob = mH->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
    if(mOpt.okfun) {Histo* ob = mH->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
  } else if(mM) {
    if(mOpt.okres) {Matrix* ob = mM->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
    if(mOpt.okfun) {Matrix* ob = mM->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
  } else if(mH2) {
    if(mOpt.okres)
      {Histo2D* ob = mH2->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
    if(mOpt.okfun)
      {Histo2D* ob = mH2->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
  } else if(mIm) {
    if(mOpt.okres) {
      RzImage* ob = mIm->FitResidus(*mFit);  
      if(ob) {ImageR4* rob = new ImageR4(*ob); omg.AddObj(rob,nomres); delete ob;}
    }
    if(mOpt.okfun) {
      RzImage* ob = mIm->FitFunction(*mFit);
      if(ob) {ImageR4* rob = new ImageR4(*ob); omg.AddObj(rob,nomfun); delete ob;}
    }
  } else if(mG) {
    if(mOpt.okres)
      {GeneralFitData* ob = mG->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
    if(mOpt.okfun)
      {GeneralFitData* ob = mG->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
  }
//####################
#else
//####################
  if(mV) {
    if(mOpt.okres) {
      Vector* ob = new Vector(ObjectFitter::FitResidus(*mV, *mFit),true);  
      if(ob) omg.AddObj(ob,nomres);
    }
    if(mOpt.okfun) {
      Vector* ob = new Vector(ObjectFitter::FitFunction(*mV, *mFit),true); 
      if(ob) omg.AddObj(ob,nomfun);
    }
  } else if(mH) {
    if(mOpt.okres) {Histo* ob =  new Histo(mH->FitResidus(*mFit));  if(ob) omg.AddObj(ob,nomres);}
    if(mOpt.okfun) {Histo* ob =  new Histo(mH->FitFunction(*mFit)); if(ob) omg.AddObj(ob,nomfun);}
  } else if(mM) {
    if(mOpt.okres) {
      Matrix* ob = new Matrix(ObjectFitter::FitResidus(*mM, *mFit), true);  
      if(ob) omg.AddObj(ob,nomres);
    }
    if(mOpt.okfun) {
      Matrix* ob = new Matrix(ObjectFitter::FitFunction(*mM, *mFit), true); 
      if(ob) omg.AddObj(ob,nomfun);
    }
  } else if(mH2) {
    if(mOpt.okres)
      {Histo2D* ob = new Histo2D(mH2->FitResidus(*mFit));  if(ob) omg.AddObj(ob,nomres);}
    if(mOpt.okfun)
      {Histo2D* ob = new Histo2D(mH2->FitFunction(*mFit)); if(ob) omg.AddObj(ob,nomfun);}
//CMV_A_FAIRE  } else if(mIm) {
//CMV_A_FAIRE    if(mOpt.okres) {
//CMV_A_FAIRE      RzImage* ob = mIm->FitResidus(*mFit);  
//CMV_A_FAIRE      if(ob) {ImageR4* rob = new ImageR4(*ob); omg.AddObj(rob,nomres); delete ob;}
//CMV _A_FAIRE    }
//CMV_A_FAIRE    if(mOpt.okfun) {
//CMV_A_FAIRE      RzImage* ob = mIm->FitFunction(*mFit);
//CMV_A_FAIRE      if(ob) {ImageR4* rob = new ImageR4(*ob); omg.AddObj(rob,nomfun); delete ob;}
//CMV_A_FAIRE    }
  } else if(mG) {
    if(mOpt.okres)
      {GeneralFitData* ob = new GeneralFitData(mG->FitResidus(*mFit));  if(ob) omg.AddObj(ob,nomres);}
    if(mOpt.okfun)
      {GeneralFitData* ob = new GeneralFitData(mG->FitFunction(*mFit)); if(ob) omg.AddObj(ob,nomfun);}
  }
//####################
#endif
//####################
}
return;
}

/* --Methode-- */
void PIAFitter::LinFit(void)
// Fit lineaire de polynomes
{
if(mFName.length()<=0)
  {cout<<"PIAFitter::LinFit Donnez un nom de fonction a fitter."<<endl;
   return;}
if(!mGData)
  {cout<<"PIAFitter::LinFit pas de Data."<<endl; return;}
if(mFName[0]!='p') 
  {cout<<"PIAFitter::LinFit fonction non prevue "<<mFName<<endl;
   return;}
int degre = 0;
if(mFName.length()>1) sscanf(mFName.c_str()+1,"%d",&degre);

if(mNVar==1) {  // Fit lineaire de polynome 1D
  cout<<"Fit (lineaire) 1D polynome de degre "<<degre<<endl;
  Poly p1(0);
  double c2rl = mGData->PolFit(0,p1,degre);
  cout<<"C2r_lineaire = "<<c2rl<<endl;
  if(mOpt.lp>0) cout<<p1<<endl;
} else if(mNVar==2) {  // Fit lineaire de polynome 2D
  cout<<"Fit (lineaire) polynome 2D de degre "<<degre<<endl;
  Poly2 p2(0);
  double c2rl = mGData->PolFit(0,1,p2,degre);
  cout<<"C2r_lineaire = "<<c2rl<<endl;
  if(mOpt.lp>0) cout<<p2<<endl;
}

return;
}

/* --Methode-- */
void PIAFitter::LinkFitFunc(string const& fname, string const& funcname
                           ,int nvar,int npar)
// Compilation et link dynamique de la fonction a fitter
{
if(nvar<1 || npar<1) return;
if(mDlUFunc != NULL) {delete mDlUFunc;  mDlUFunc=NULL;}
mFitFunc = NULL; mFitFuncDer = NULL; 
mUFNVar= mUFNPar = 0; mNameFitFunc = "";

// On compile le fichier C, on fabrique un .so et on l'ouvre
mDlUFunc = PDynLinkMgr::BuildFromCFile(fname);
if(mDlUFunc == NULL)
  {cerr<<"PIAFitter::LinkFitFunc() Erreur building .so "<<endl; return;}

// On recupere la fonction de fit
mFitFunc = (DlUserFitFunc) mDlUFunc->GetFunction(funcname);
if(mFitFunc == NULL)  // Probleme, on n'a pas trouve etc ... 
  {cerr<<"PIAFitter::LinkFitFunc() Erreur linking function"<<funcname<<endl;
   delete mDlUFunc; mDlUFunc = NULL; return;}

// On recupere la fonction des derivees
string fnameder = funcname+"_der";
mFitFuncDer = (DlUserFitFuncDer) mDlUFunc->GetFunction(fnameder);
if(mFitFuncDer == NULL)
  {cout<<"PIAFitter::LinkFitFunc() No derive for"<<fnameder<<endl;}

// Variables et pointeurs definissant le link a la fonction utilisateur
mNameFitFunc = funcname;
mUFNVar = nvar;
mUFNPar = npar;
cout<<"PIAFitter::LinkFitFunc() UserFitFunc="<<mNameFitFunc
    <<" NVar="<<mUFNVar<<" NPar="<<mUFNPar<<endl;

return;
}

/* --Methode-- */
void PIAFitter::PrepareFitFuncCFile(string const& fname, string const& funcname
                                   , int nvar, int npar)
// Pour ecrire un squelette des fichier C avec pour une fonction a fitter
{
if(fname.length()<=0 || funcname.length()<=0 || nvar<=0 || npar <=0) {
  cerr<<"PIAFitter::PrepareFitFuncCFile: error "
      <<" nvar="<<nvar<<" npar="<<npar<<endl
      <<" file_name="<<fname<<" func_name="<<funcname<<endl
      <<" ---- No File Created ! ----"<<endl;
  return;
}

FILE* fip = fopen(fname.c_str(),"w");
if(fip==NULL)
  {cerr<<"PIAFitter::PrepareFitFuncCFile: error opening file for writing\n"
       <<fname<<endl; return;}

fprintf(fip,"#include <stdio.h>\n#include <stdlib.h>\n#include <math.h>\n\n");
fprintf(fip,"/* - many functions can be written in this file: */\n");
fprintf(fip,"/*   here is an exemple for creation of function %s */\n\n"
       ,funcname.c_str());
fprintf(fip,"double %s(double const* x,double const* p)\n",funcname.c_str());
fprintf(fip,"/* Function : F(x[0-%d];p[0-%d]) */\n",nvar-1,npar-1);
fprintf(fip,"/* variables  x1=x[0],x2=x[1],...,x%d=x[%d] */\n",nvar,nvar-1);
fprintf(fip,"/* parameters p1=p[0],p2=p[1],...,p%d=p[%d] */\n",nvar,nvar-1);
fprintf(fip,"{\n  val = F(x[];p[]);\n  return(val);\n}\n\n");
string fder = funcname + "_der";
fprintf(fip,"/* That function is optional */\n");
fprintf(fip,"/*\n");
fprintf(fip,"double %s(double const* x,double const* p,double* dp)\n",fder.c_str());
fprintf(fip,"{\n");
for(int i=0;i<npar;i++) fprintf(fip,"  dp[%d] = dF(x[];p[])/dp[%d]\n",i,i);
fprintf(fip,"  val = F(x[];p[]);\n  return(val);\n}\n\n");
fprintf(fip,"*/\n");

fclose(fip);
}

/* --Methode-- */
int PIAFitter::Execute(string& kw, vector<string>& tokens)
// Execution des commandes de fit : "fit" et "fitw"
{
////// Les commandes pour creer les fonctions de fit utilisateur
if(kw == "crfitfil") {
  if (tokens.size() < 4)
    {cout<<"Usage:crfitfil filename funcname nvar npar"<<endl;return(-1);}
  int nvar = atoi(tokens[2].c_str()), npar = atoi(tokens[3].c_str());
  PrepareFitFuncCFile(tokens[0],tokens[1],nvar,npar);
  return(0);
} else if (kw == "crfitfun") {
  if (tokens.size() < 4)
    {cout<<"Usage:crfitfun filename funcname nvar npar"<<endl;return(-1);}
  int nvar = atoi(tokens[2].c_str()), npar = atoi(tokens[3].c_str());
  LinkFitFunc(tokens[0],tokens[1],nvar,npar);
  return(0);
}

////// Les commandes pour fitter
if(tokens.size() < 2) {
  cout
  <<"Usage:fit  nomobj func [Options]\n"
  <<" ou   fitw nomobj func [Options]\n"
  <<" avec Options: [p:p1,...,pn s:s1,...,sn m:m1,...,mn \n"
  <<"                M:M1,...,Mn f:f1,...,fn o:.,.,. o:.,.,.]\n"
  <<" ou   fitlin nomobj func [o:.,.,. o:.,.,.]\n";
  return(-1);
}

// decodage des arguments de la commande
string mSp,mSs,mSm,mSM,mSf,mSo;
mSp=""; mSs=""; mSm=""; mSM=""; mSf=""; mSo="";
if(tokens.size()>2)
  for(int ip=2;ip<(int)tokens.size();ip++) {
    if(tokens[ip].length()<=2) continue;
    const char *c = tokens[ip].c_str();
    if(c[1]!=':') continue;
    if(c[0]=='p')      mSp=c+2;
    else if(c[0]=='s') mSs=c+2;
    else if(c[0]=='m') mSm=c+2;
    else if(c[0]=='M') mSM=c+2;
    else if(c[0]=='f') mSf=c+2;
    else if(c[0]=='o') {mSo += ","; mSo += c+2;}
  }

// Execution des commandes
DecodeOptions(mSo);
DecodeObject(tokens[0]);
CheckOptions();
if(mOpt.lp>0) PrintOptions();
FillGData();
DecodeFunction(tokens[1]);
DecodeParam(mSp,mSs,mSm,mSM,mSf);
if (kw == "fit") {
  int rc = DoFit();
  if(rc>=0) FitFunRes();
} else if (kw == "fitlin") {
  LinFit();
} else if (kw == "fitw") {
  if(FWindFit!=NULL) {delete FWindFit; FWindFit=NULL;}
  FWindFit = new PIAFitterWind(mApp,this);
  FWindFit->Show();
}

return(0);
}

/* --Methode-- */
void PIAFitter::IniFitP1ou2D(void)
{
if(mNPar==0) return;
mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
if(!mGData) return;
double mean,sigma;
int rc = mGData->GetMeanSigma(0,mean,sigma);
if(rc<=0) return;
mPar(0) = mean;
return;
}

void PIAFitter::IniFitGhP1D(void)
{
if(mNPar==0) return;
mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
if(mNVar!=1 || mNPar<3) return;
double h,v,x0,y0,sx,sy,f;
int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
if(rc!=0) return;
mPar(0)=h;  mStep(0)=h/1000.;
mPar(1)=x0; mStep(1)=sx/10.;
mPar(2)=sx; mStep(2)=sx/10.; mMin(2)=0.; mMax(2)=20.*sx; 
if(mNPar>3) {mPar(3)=f; mStep(3)=f/1000.;} else mPar(0)=h+f;
}

void PIAFitter::IniFitGvP1D(void)
{
if(mNPar==0) return;
mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
if(mNVar!=1 || mNPar<3) return;
double h,v,x0,y0,sx,sy,f;
int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
if(rc!=0) return;
mPar(0)=v;  mStep(0)=v/1000.;
mPar(1)=x0; mStep(1)=sx/10.;
mPar(2)=sx; mStep(2)=sx/10.; mMin(2)=0.; mMax(2)=20.*sx; 
if(mNPar>3) {mPar(3)=f; mStep(3)=f/1000.;} else mPar(0)=v+mNData*f;
}

void PIAFitter::IniFitGv2D(void)
{
if(mNPar==0) return;
mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
if(mNVar!=2 || mNPar<6) return;
double h,v,x0,y0,sx,sy,f;
int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
if(rc!=0) return;
mPar(0)=v;  mStep(0)=v/1000.;
mPar(1)=x0; mStep(1)=sx/10.;
mPar(2)=y0; mStep(2)=sy/10.;
mPar(3)=sx; mStep(3)=sx/10.; mMin(3)=0.;  mMax(3)=20.*sx; 
mPar(4)=sy; mStep(4)=sy/10.; mMin(4)=0.;  mMax(4)=20.*sy; 
mPar(5)=0.; mStep(5)=0.005;  mMin(5)=-1.; mMax(5)=1.; 
if(mNPar>6) {mPar(6)=f; mStep(6)=f/1000.;} else mPar(0)=v+mNData*f;
}

int PIAFitter::IniFitGaus(GeneralFitData& g,double& h,double& v
    ,double& x0,double& y0,double& sx,double& sy,double& f)
// Pour initialiser une gaussienne 1D ou 2D
// return code =0 si OK
{
h=v=x0=y0=sx=sy=f=0.;
if(!&g) return(-1);
int nvar = g.NVar();
int ndata = g.NData();
if(nvar<1 || 2<nvar || ndata<=0) return(-2);

// X et Y ordonnees croissant
double *X=NULL, *Y=NULL;
X = new double[ndata];
if(nvar==2) Y = new double[ndata];
{for(int i=0;i<ndata;i++) {X[i] = g.X(i); if(Y) Y[i] = g.Y(i);}}
qsort(X,(size_t) ndata,sizeof(double),qSort_Dble);
if(Y) qsort(Y,(size_t) ndata,sizeof(double),qSort_Dble);

// determination des zones externes=[0,1/3] et [2/3,1]
//            et des zones internes=[1/3,2/3]
int k1=ndata/3,k2=2*ndata/3;
double xmin,xmax, ymin=-1.,ymax=1.;
xmin = X[k1]; xmax = X[k2];
if(Y) {ymin = Y[k1]; ymax = Y[k2];}

// recherche des parametres
double nb=0, nc=0;
{for(int i=0;i<ndata;i++) {
   double val,x,y=0.;
   x = g.X(i); if(Y) y = g.Y(i); val = g.Val(i);
   if(xmin<x && x<xmax && ymin<y && y<ymax) { // c'est le centre
     if(nc==0 || val>h) {h=val; x0=x; y0=y;}
     v += val; nc++;
   } else {f += val; nb++;} // c'est le bord

}}
if(nb==0 || nc==0) {delete [] X; if(Y) delete [] Y; return(-3);}

// Calcul des parametres
f /= (double) nb; h -= f; v -= f*nc;
sx = (xmax-xmin)/2.; if(Y) sy = (ymax-ymin)/2.; 

delete [] X; if(Y) delete [] Y;
return(0);
}


////////////////////////////////////////////////////////////////////
////---------- Classe de fenetre de fit interactive ------------////
////////////////////////////////////////////////////////////////////

/* --Methode-- */
PIAFitterWind::PIAFitterWind(PIStdImgApp* par, PIAFitter* fiter)
: PIWindow((PIMsgHandler*)par, "PIAFitter", PIWK_normal, 240, 240, 150, 150)
, mDap(par), mFitter(fiter), ReFillGData(false)
, saveI1(0), saveI2(0), saveJ1(0), saveJ2(0), ErrValue(1.)
{
#ifdef SANS_EVOLPLANCK
COMPATIBLE_ASSERT(mFitter);
#else
ASSERT(mFitter);
#endif

int npar = mFitter->mNPar;

// Alloc de la taille
nlab = 13+npar;   lab = new PILabel*[nlab];    {for(int i=0;i<nlab;i++) lab[i]=NULL;}
ntxt = 10+4*npar; txt = new PIText*[ntxt];     {for(int i=0;i<ntxt;i++) txt[i]=NULL;}
nbut = 6;         but = new PIButton*[nbut];   {for(int i=0;i<nbut;i++) but[i]=NULL;}
nckb = 2+npar;    ckb = new PICheckBox*[nckb]; {for(int i=0;i<nckb;i++) ckb[i]=NULL;}
npom = 3;         pom = new PIOptMenu*[npom];  {for(int i=0;i<npom;i++) pom[i]=NULL;}

// On definit la taille a partir de la taille par defaut des composantes
// bsx,bsy = environ 6 lettres
int bsx,bsy;
par->PrefCompSz(bsx, bsy);
int spx = (bsx>=10) ? bsx/10 : 1;  // intervalle entre lettres X
int spy = (bsy>=6) ? bsy/6 : 1;    // intervalle entre lettres Y

int wszx = (5*bsx+5*spx)+spx;
int wszy = 8*(bsy+spy)+npar*(bsy+spy)+2*spy;
SetSize(wszx, wszy);
int cpx,cpy;

// new ligne
cpx = spx; cpy = spy;
string dum = "Object: "+mFitter->mNObj+" , Fun: ";
if(mFitter->mFunc) dum+=mFitter->mNameFitFunc; else dum+=mFitter->mFName;
lab[0] = new PILabel(this,dum.c_str(),4*bsx,bsy,cpx,cpy);
lab[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);

// new ligne
cpx = spx; cpy += bsy+spy;
but[0] = new PIButton(this,"Fr.Last",111,bsx,bsy,cpx,cpy);
but[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
but[1] = new PIButton(this,"Default",222,bsx,bsy,cpx,cpy);
but[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
ckb[0] = new PICheckBox(this,"Gen Fun",1001,bsx,bsy,cpx,cpy);
ckb[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+2*spx;
ckb[1] = new PICheckBox(this,"Res",1002,bsx,bsy,cpx,cpy);
ckb[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);

// new ligne
cpx = spx; cpy += bsy+spy;
pom[0] = new PIOptMenu(this,"Err def",1.5*bsx,bsy,cpx,cpy);
  pom[0]->AppendItem("Err def",2101);
  pom[0]->AppendItem("Err cste",2102);
  pom[0]->AppendItem("Err sqrt",2103);
  pom[0]->SetValue(2101);
pom[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(1.5*bsx+spx);
txt[0] = new PIText(this,"",bsx,bsy,cpx,cpy);
txt[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
lab[1] = new PILabel(this,"Prt,dbg",bsx,bsy,cpx,cpy);
lab[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
txt[1] = new PIText(this,"",bsx,bsy,cpx,cpy);
txt[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);

// new ligne
cpx = spx; cpy += bsy+spy;
lab[2] = new PILabel(this,"Stop Xi2",1.25*bsx,bsy,cpx,cpy);
lab[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(1.25*bsx+spx);
txt[2] = new PIText(this,"",bsx,bsy,cpx,cpy);
txt[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
lab[3] = new PILabel(this,"Iter",bsx, bsy, cpx, cpy);
lab[3]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
txt[3] = new PIText(this,"",bsx,bsy,cpx,cpy);
txt[3]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);

// new ligne
cpx = spx; cpy += bsy+spy;
lab[4] = new PILabel(this,"Range X",bsx,bsy,cpx,cpy);
lab[4]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
txt[4] = new PIText(this,"",0.75*bsx,bsy,cpx,cpy);
txt[4]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(0.75*bsx+spx);
txt[5] = new PIText(this,"",0.75*bsx,bsy,cpx,cpy);
txt[5]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(0.75*bsx+spx);
lab[5] = new PILabel(this,"Y",bsx/2,bsy,cpx,cpy);
lab[5]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(bsx/2.+spx);
txt[6] = new PIText(this,"",0.75*bsx,bsy,cpx,cpy);
txt[6]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(0.75*bsx+spx);
txt[7] = new PIText(this,"",0.75*bsx,bsy,cpx,cpy);
txt[7]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);

// new ligne
cpx = spx; cpy += bsy+spy;
lab[6] = new PILabel(this,"Cent",bsx/2,bsy,cpx,cpy);
lab[6]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(bsx/2.+spx);
pom[1] = new PIOptMenu(this,"X.No",1.25*bsx,bsy,cpx,cpy);
  pom[1]->AppendItem("X.No",2001);
  pom[1]->AppendItem("X.Value",2002);
  pom[1]->AppendItem("X.Auto",2003);
  pom[1]->SetValue(2001);
pom[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(1.25*bsx+spx);
txt[8] = new PIText(this,"",bsx,bsy,cpx,cpy);
txt[8]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
pom[2] = new PIOptMenu(this,"Y.No",1.25*bsx,bsy,cpx,cpy);
  pom[2]->AppendItem("Y.No",2011);
  pom[2]->AppendItem("Y.Value",2012);
  pom[2]->AppendItem("Y.Auto",2013);
  pom[2]->SetValue(2011);
pom[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(1.25*bsx+spx);
txt[9] = new PIText(this,"",bsx,bsy,cpx,cpy);
txt[9]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);

// new ligne
cpx = spx; cpy += bsy+spy;
but[2] = new PIButton(this, "Dismiss",777,bsx,bsy,cpx,cpy);
but[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
but[3] = new PIButton(this, "Check",666,bsx,bsy,cpx,cpy);
but[3]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
but[4] = new PIButton(this, "Fit",333,bsx,bsy,cpx,cpy);
but[4]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
but[5] = new PIButton(this, "Ini Auto",444,bsx,bsy,cpx,cpy);
but[5]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);

// new ligne
cpx = spx; cpy += bsy+spy;
lab[7] = new PILabel(this,"Par",bsx/2,bsy,cpx,cpy);
lab[7]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(bsx/2.+spx);
lab[8] = new PILabel(this,"Fx",bsx/3,bsy,cpx,cpy);
lab[8]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += int(bsx/3.+spx);
lab[9] = new PILabel(this,"Init",bsx,bsy,cpx,cpy);
lab[9]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
lab[10] = new PILabel(this,"Step",bsx,bsy,cpx,cpy);
lab[10]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
lab[11] = new PILabel(this,"Min",bsx,bsy,cpx,cpy);
lab[11]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
cpx += bsx+spx;
lab[12] = new PILabel(this,"Max",bsx,bsy,cpx,cpy);
lab[12]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);

// new lines: Parametres (npar lignes)
if(npar>0) {
  for(int i=0;i<npar;i++) {
    cpx = spx; cpy += bsy+spy;
    char str[8]; sprintf(str,"P%d",i);
    lab[13+i] = new PILabel(this,str,bsx/2,bsy,cpx,cpy);
    lab[13+i]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
    cpx += int(bsx/2.+spx);
    // Attention pas de message entre 1500 et 2000 (permet 500 parametres!)
    ckb[2+i] = new PICheckBox(this,"",1500+i,bsx/3,bsy,cpx,cpy+spy);
    ckb[2+i]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
    cpx+=int(bsx/3.+spx);
    for(int j=0;j<4;j++) {
      txt[10+4*i+j] = new PIText(this,"",bsx,bsy,cpx,cpy);
      txt[10+4*i+j]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
      cpx += bsx+spx;
    }
  }
}

return;
}

/* --Methode-- */
PIAFitterWind::~PIAFitterWind()
{
int i;
for(i=0;i<nlab;i++) if(lab[i]) {delete lab[i]; lab[i]=NULL;}
for(i=0;i<ntxt;i++) if(txt[i]) {delete txt[i]; txt[i]=NULL;}
for(i=0;i<nbut;i++) if(but[i]) {delete but[i]; but[i]=NULL;}
for(i=0;i<nckb;i++) if(ckb[i]) {delete ckb[i]; ckb[i]=NULL;}
for(i=0;i<npom;i++) if(pom[i]) {delete pom[i]; pom[i]=NULL;}
if(lab) {delete lab; lab=NULL; nlab=0;}
if(txt) {delete txt; txt=NULL; ntxt=0;}
if(but) {delete but; but=NULL; nbut=0;}
if(ckb) {delete ckb; ckb=NULL; nckb=0;}
if(pom) {delete pom; pom=NULL; npom=0;}
}

/* --Methode-- */
void PIAFitterWind::Show()
{
//mDap->SetBlocked();
SetState();
PIWindow::Show();
return;
}

/* --Methode-- */
void PIAFitterWind::SetState(void)
// Ecriture et positionnement des ckb et txt
{
int npar = mFitter->mNPar;

// Les check-boxes
ckb[0]->SetState(mFitter->mOpt.okfun);
ckb[1]->SetState(mFitter->mOpt.okres);
if(npar>0)
  for(int i=0;i<npar;i++) ckb[2+i]->SetState((mFitter->mFix(i)>=1.)?true:false);

// Les champs textes
char str[128]; string dum;
sprintf(str,"%f",ErrValue); dum=str; txt[0]->SetText(dum);
sprintf(str,"%d,%d",mFitter->mOpt.lp,mFitter->mOpt.lpg); dum=str; txt[1]->SetText(dum);
sprintf(str,"%f",mFitter->mOpt.stc2);  dum=str; txt[2]->SetText(dum);
sprintf(str,"%d",mFitter->mOpt.nstep); dum=str; txt[3]->SetText(dum);
sprintf(str,"%d",mFitter->mOpt.i1);    dum=str; txt[4]->SetText(dum);
sprintf(str,"%d",mFitter->mOpt.i2);    dum=str; txt[5]->SetText(dum);
sprintf(str,"%d",mFitter->mOpt.j1);    dum=str; txt[6]->SetText(dum);
sprintf(str,"%d",mFitter->mOpt.j2);    dum=str; txt[7]->SetText(dum);
sprintf(str,"%f",mFitter->mOpt.xc);    dum=str; txt[8]->SetText(dum);
sprintf(str,"%f",mFitter->mOpt.yc);    dum=str; txt[9]->SetText(dum);
if(npar>0) {
  for(int i=0;i<npar;i++) {
    sprintf(str,"%f",mFitter->mPar(i));  dum=str; txt[10+4*i+0]->SetText(dum);
    sprintf(str,"%f",mFitter->mStep(i)); dum=str; txt[10+4*i+1]->SetText(dum);
    sprintf(str,"%f",mFitter->mMin(i));  dum=str; txt[10+4*i+2]->SetText(dum);
    sprintf(str,"%f",mFitter->mMax(i));  dum=str; txt[10+4*i+3]->SetText(dum);
  }
}

return;
}

/* --Methode-- */
void PIAFitterWind::Process(PIMessage msg, PIMsgHandler* sender, void* data)
{
int lp = (mFitter->mOpt.lp>2)?1:0;;
int npar = mFitter->mNPar;

// Decodage des messages, init et dismiss
msg = UserMsg(msg);
if(msg ==777) {
  if(lp) cout<<"Dismiss"<<endl;
  //mDap->SetReady();
  Hide();
  return;
}

// L'objet existe t-il encore?
NamedObjMgr omg;
if(omg.GetObj(mFitter->mNObj) == NULL) {
  cout<<"PIAFitterWind::Process Error , Pas d'objet de nom "
      <<mFitter->mNObj<<endl;
  mFitter->ResetDPointer(); return;
}

// On recupere les champs textes
string dum; char str[128];
// Attention gestion du type d'erreur
dum=txt[0]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%lf",&ErrValue); if(ErrValue<=0.) ErrValue=1.;
  if(mFitter->mOpt.err_e>0. ) mFitter->mOpt.err_e = ErrValue;
  if(mFitter->mOpt.err_E>0. ) mFitter->mOpt.err_E = ErrValue;
dum=txt[1]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%d,%d",&(mFitter->mOpt.lp),&(mFitter->mOpt.lpg));
dum=txt[2]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%lf",&(mFitter->mOpt.stc2));
dum=txt[3]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%d",&(mFitter->mOpt.nstep));
// Attention si les valeurs de mOpt.i1/.i2/.j1/.j2 ont ete changees
dum=txt[4]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%d",&(mFitter->mOpt.i1));
dum=txt[5]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%d",&(mFitter->mOpt.i2));
dum=txt[6]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%d",&(mFitter->mOpt.j1));
dum=txt[7]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%d",&(mFitter->mOpt.j2));
  if(saveI1!=mFitter->mOpt.i1 || saveI2!=mFitter->mOpt.i2 ||
     saveJ1!=mFitter->mOpt.j1 || saveJ2!=mFitter->mOpt.j2 ) {
    if(lp)
      cout<<"Range changed: I("<<saveI1<<","<<saveI2
          <<") -> I("<<mFitter->mOpt.i1<<","<<mFitter->mOpt.i2
          <<")  ou   J("<<saveJ1<<","<<saveJ2
          <<") -> J("<<mFitter->mOpt.j1<<","<<mFitter->mOpt.j2
          <<")"<<endl;
    ReFillGData = true;
    saveI1=mFitter->mOpt.i1; saveI2=mFitter->mOpt.i2;
    saveJ1=mFitter->mOpt.j1; saveJ2=mFitter->mOpt.j2;
  }
dum=txt[8]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%lf",&(mFitter->mOpt.xc));
dum=txt[9]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
  sscanf(str,"%lf",&(mFitter->mOpt.yc));

if(npar>0) {
  for(int i=0;i<npar;i++) {
    dum=txt[10+4*i+0]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
      sscanf(str,"%lf",&(mFitter->mPar(i)));
    dum=txt[10+4*i+1]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
      sscanf(str,"%lf",&(mFitter->mStep(i)));
    dum=txt[10+4*i+2]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
      sscanf(str,"%lf",&(mFitter->mMin(i)));
    dum=txt[10+4*i+3]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
      sscanf(str,"%lf",&(mFitter->mMax(i)));
  }
}

// Decodage des messages, suite
if(msg ==666) {
  if(lp) cout<<"Check"<<endl;
  mFitter->CheckOptions();
  if(ReFillGData) mFitter->FillGData();
  ReFillGData = false;
  if(mFitter->mFit) mFitter->mFit->PrintStatus();
}
else if(msg ==111) {
  if(lp) cout<<"Update from last"<<endl;
  mFitter->ResetOptions();
  mFitter->mOpt = mFitter->mOptSave;
  mFitter->ReSetParam();
  mFitter->InitParFromLastFit();
  ReFillGData = true;
}
else if(msg ==222) {
  if(lp) cout<<"Default options/par"<<endl;
  mFitter->ResetOptions();
  mFitter->ReSetParam();
  ReFillGData = true;
}
else if(msg ==333) {
  if(lp) cout<<"Do the fit"<<endl;
  mFitter->CheckOptions();
  if(ReFillGData) mFitter->FillGData();
  ReFillGData = false;
  int rc = mFitter->DoFit();
  if(rc>=0) mFitter->FitFunRes();
}
else if(msg ==444) {
  if(lp) cout<<"Automatic initialisation"<<endl;
  mFitter->AutoIniFit();
}
else if(msg==1001) {
  if(ckb[0]->GetState()) {
    if(lp) cout<<"Gen fitted function"<<endl;
    mFitter->mOpt.okfun = true;
  } else if(!ckb[0]->GetState()) {
    if(lp) cout<<"No Gen fitted function"<<endl;
    mFitter->mOpt.okfun = false;
  }
}
else if(msg==1002) {
  if(ckb[1]->GetState()) {
    if(lp) cout<<"Gen fitted residus"<<endl;
    mFitter->mOpt.okres = true;
  } else if(!ckb[1]->GetState()) {
    if(lp) cout<<"No Gen fitted residus"<<endl;
    mFitter->mOpt.okres = false;
  }
}
else if(2101<=msg && msg<=2103) {
  mFitter->mOpt.err_e = mFitter->mOpt.err_E = -1.;
  if(msg==2102)      mFitter->mOpt.err_e = ErrValue;
  else if(msg==2103) mFitter->mOpt.err_E = ErrValue;
  if(lp) cout<<"Err_cste="<<mFitter->mOpt.err_e
             <<" Err_sqrt="<<mFitter->mOpt.err_E
             <<" ErrValue="<<ErrValue<<endl;
  ReFillGData = true;
}
else if(2001<=msg && msg<=2003) {
  mFitter->mOpt.polcx = msg-2001;
  if(lp) cout<<"Centrage X polcx="<<mFitter->mOpt.polcx
             <<" xc="<<mFitter->mOpt.xc<<endl;
  ReFillGData = true;
}
else if(2011<=msg && msg<=2013) {
  mFitter->mOpt.polcy = msg-2011;
  if(lp) cout<<"Centrage Y polcy="<<mFitter->mOpt.polcy
             <<" yc="<<mFitter->mOpt.yc<<endl;
  ReFillGData = true;
}
else if(msg>=1500 && msg<2000) {
  int ip = msg-1500;
  if(ckb[2+ip]->GetState()) {
    if(lp) cout<<"Fix Param "<<ip<<endl;
    mFitter->mFix(ip) = 1.001;
  } else if(!ckb[2+ip]->GetState()) {
    if(lp) cout<<"Not Fix Param "<<ip<<endl;
    mFitter->mFix(ip) = 0.;
  }
}

mFitter->CheckOptions();
SetState();
if(mFitter->mOpt.lp>1) mFitter->PrintOptions();

return;
}