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

Last change on this file since 1200 was 1200, checked in by ercodmgr, 25 years ago

bug dans fiwt (avec fenetre) cmv 23/9/00
File size: 64.4 KB
Line 
1//CMVBUG reste pb suivant sous Linux-egcs (OSF1-cxx c'est ok):
2//openppf cmvwppf.ppf
3//fitw h1g g2 p:4500,45,10,900
4//delobjs h1g
5//clic sur FIT -> OK
6//clic sur FIT 2sd fois -> core dump!!!
7//
8#include <stdio.h>
9#include <stdlib.h>
10#include <ctype.h>
11#include <iostream.h>
12#include <typeinfo>
13
14#include "strutil.h"
15#include "nbtri.h"
16#include "generalfit.h"
17#include "fct1dfit.h"
18#include "fct2dfit.h"
19#include "ntuple.h"
20#include "cimage.h"
21#include "histos.h"
22#include "histos2.h"
23#include "ntuple.h"
24#include "hisprof.h"
25
26#ifdef SANS_EVOLPLANCK
27#include "matrix.h"
28#include "cvector.h"
29#else
30#include "tmatrix.h"
31#include "tvector.h"
32#include "objfitter.h"
33#endif
34
35#include "piafitting.h"
36#include "nobjmgr.h"
37#include "pistdimgapp.h"
38
39////////////////////////////////////////////////////////////////////
40////---------- Classe de fenetre de fit interactive ------------////
41////////////////////////////////////////////////////////////////////
42class PIAFitterWind : public PIWindow {
43public :
44 PIAFitterWind(PIStdImgApp* par, PIAFitter* fiter,string func);
45 virtual ~PIAFitterWind();
46 virtual void Show();
47 virtual void Process(PIMessage msg, PIMsgHandler* sender, void* data=NULL);
48 void SetState(void);
49protected :
50 PIStdImgApp* mDap;
51 PIAFitter* mFitter;
52 string mFunc;
53 PILabel** lab; int nlab;
54 PIText** txt; int ntxt;
55 PIButton** but; int nbut;
56 PICheckBox** ckb; int nckb;
57 PIOptMenu** pom; int npom;
58 bool ReFillGData;
59 bool ReDecodeFunc;
60 double ErrValue;
61 int saveI1,saveI2, saveJ1, saveJ2;
62 double saveErrValue,saveXC,saveYC;
63};
64
65////////////////////////////////////////////////////////////////////////
66// --- Generation automatique du "usage" pour helpfit:
67// 1-/ copier le texte ci-dessous de piafitting.cc dans toto
68// 2-/ cat toto | sed 's?^//?usage+="\\n?' | sed 's?$? ";?'
69// NE PAS EFFACER LA SUITE: c'est elle qui doit etre modifiee,
70// le usage du helpfit sera regenere a l'aide des commandes ci-dessus.
71////////////////////////////////////////////////////////////////////////
72//| --------------- Fit Lineaire a 1 et 2 dimensions ---------------
73//| Syntaxe:
74//| fitlin nom pnn [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]
75//| avec:
76//| nom : cf commentaire ordre ''fit''
77//| pnn : fit polynome degre nn avec classe Poly (lineaire) 1D ou 2D
78//| o:aa,...,bb : cf commentaires ordre ''fit''
79//|
80////////////////////////////////////////////////////////////////////////
81//| --------------- Fit d'objets a 1 et 2 dimensions ---------------
82//| avec interface d'aide graphique
83//| Syntaxe:
84//| fitw nom func
85//| [p:p1,p2,...,pn s:s1,s2,...,sn m:m1,m2,...,mn M:M1,mM2,...,Mn f:f1,...,fn]
86//| [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]
87//| cf commentaire ordre ''fit''
88//|
89////////////////////////////////////////////////////////////////////////
90//| --------------- Fit d'objets a 1 et 2 dimensions ---------------
91//| Syntaxe:
92//| fit nom func
93//| [p:p1,p2,...,pn s:s1,s2,...,sn m:m1,m2,...,mn M:M1,mM2,...,Mn f:f1,...,fn]
94//| [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...]
95//|----- OBJET -----
96//| nom : nom de l'objet qui peut etre:
97//| fit-1D: Vector,Histo1D,HProf ou GeneraFitData(1D)
98//| fit-2D: Matrix,Histo2D,Image<T> ou GeneraFitData(2D)
99//|
100//|----- FUNCTION -----
101//| func : pnn : fit polynome degre nn avec GeneralFit (non-lineaire) 1D ou 2D
102//| : gnn : fit gaussienne (hauteur) + polynome de degre nn 1D
103//| : g : fit gaussienne (hauteur) 1D
104//| : enn : fit exponentielle + polynome de degre nn 1D
105//| : e : fit exponentielle 1D
106//| : Gnn : fit gaussienne (volume) + polynome de degre nn 1D
107//| : G : fit gaussienne (volume) 1D
108//| : : fit gaussienne+fond (volume) 2D
109//| : Gi : fit gaussienne+fond integree (volume) 2D
110//| : d : fit DL de gaussienne+fond (volume) 2D
111//| : di : fit DL de gaussienne+fond integree (volume) 2D
112//| : D : fit DL de gaussienne+fond avec coeff variable p6 (volume) 2D
113//| : Di : fit DL de gaussienne+fond integree avec coeff variable p6 (volume) 2D
114//| : M : fit Moffat+fond (expos=p6) (volume) 2D
115//| : Mi : fit Moffat+fond integree (expos=p6) (volume) 2D
116//| : Autre : fonction definie par l'utilisateur (cf commentaires ci apres)
117//|
118//|----- INIT PARAMETRES ET ETAT DU FIT -----
119//| p : p1,...,pn : valeur d'initialisation des parametres (def=0)
120//| s : s1,...,sn : valeur des steps de depart (def=1)
121//| m : m1,...,mn : valeur des minima (def=1)
122//| M : M1,...,Mn : valeur des maxima (def=-1) (max<=min : pas de limite)
123//| f : f1,...,fn : si >=1 parametre fixe sinon libre (def=0)
124//| - Remarque: si pi,si,mi ou Mi = '!' la valeur correspondante n'est pas changee
125//|
126//|----- OPTIONS -----
127//| o : options ''o:Eaa.b,eaa.b,f,r,caa.b,Xaa.b''
128//| F : initialisation a partir des resultats et de l'etat du fit precedent
129//| (option non prioritaire sur les definitions p:,s:,m:,M:,f:o:)
130//| a : initialisation automatique des parametres du fit (si disponible)
131//| (option non prioritaire sur les definitions p:,s:,m:,M:,f:o:
132//| et sur l'initialisation pa l'option ''F'')
133//| f : generation d'un Objet identique contenant la fonction fittee
134//| r : generation d'un Objet identique contenant les residus
135//| Xaa.b : aa.b valeur du DXi2 d'arret (def=1.e-3)
136//| Naa : aa nombre maximum d'iterations (def=100)
137//| la : niveau ''a'' de print: a=niveau de print Fit1/2D
138//| da : niveau ''a'' de debug: a=niveau de GeneralFit
139//| Ii1/i2 : numeros des bins X de l'histos utilises pour le fit [i1,i2]
140//|2D Jj1/j2 : numeros des bins Y de l'histos utilises pour le fit [j1,j2]
141//| - L'erreur est celle associee a l'objet (si elle existe),
142//| elle est mise a 1 sinon, sauf si E... ou e... est precise:
143//| Eaa.b : si |val|>=1 erreur = aa.b*sqrt(|val|)
144//| si |val|<1 erreur = aa.b
145//| si aa.b <=0 alors aa.b=1.0
146//| E seul est equivalent a E1.0
147//| eaa.b : erreur = aa.b
148//| si aa.b <=0 alors aa.b=1.0
149//| e seul est equivalent a e1.0
150//| xaa.b : demande de centrage: on fit x-aa.b au lieu de x)
151//| x : demande de centrage: on fit x-xc au lieu de x
152//| avec xc=abscisse du milieu de l'histogramme
153//| Actif pour exp+poly 1D, poly 1D
154//| Pour gauss+poly 1D, xc est le centre de la gaussienne.
155//|2D yaa.b et y : idem ''xaa.b et x'' mais pour y
156//|
157////////////////////////////////////////////////////////////////////////
158//| --------------- Def des fcts de fit par l'utilisateur ---------------
159//| Pour definir une fonction parametree utilisateur:
160//| Un exemple de fichier contenant des fonctions a fitter: PIext/userfitfunex.c
161//| > crfitfil nom_fichier nom_fonction nvar npar
162//| pour creer un squelette de fichier C contenant la fonction a fitter
163//| nom_fichier : nom du fichier ou est le code de la fonction
164//| c'est un fichier C de nom blabla.c
165//| nom_fonction : nom de la fonction ex: gauss2
166//| nvar : nombre de variables x[] (1D=1, 2D=2, etc...)
167//| npar : nombre de parametre p[]
168//| - La fonction est donc F(x[];p[]), soit par ex:
169//| double gauss2(double const* x,double const* p)
170//| - L'utilisateur a la possibilite de donner egalement la fonction
171//| retournant les derivees de la fonction par rapport aux parametres
172//| double gauss2_der(double const* x,double const* p,double* dp)
173//| ou dp[i] = dF(x[];p[])/dp[i]
174//| Par convention, cette fonction a le meme nom suivi de ''_der''
175//| > crfitfun nom_fichier nom_fonction nvar npar
176//| Permet de compiler le fichier ''nom_fichier'' et de linker
177//| la fonction ''nom_fonction''
178//| - Un meme fichier peut contenir plusieurs fonctions a fitter
179//|
180////////////////////////////////////////////////////////////////////////
181
182/* --Methode-- */
183PIAFitter::PIAFitter(PIACmd *piac, PIStdImgApp* app)
184: mApp(app), FWindFit(NULL)
185, mNObj(""), mObj(NULL), mGData(NULL)
186, mNPar(0), mNVar(0), mNBinX(0), mNBinY(0), mNData(0)
187, mPar(1), mStep(1), mMin(1), mMax(1), mFix(1)
188, mParSave(1), mStepSave(1), mMinSave(1), mMaxSave(1), mFixSave(1)
189, mFit(NULL),ApplyXYcentering(false)
190, mV(NULL), mH(NULL), mM(NULL), mH2(NULL), mIm(NULL), mG(NULL)
191, mFunction(NULL), mFName(""), mFunc(NULL)
192, mDlUFunc(NULL), mNameFitFunc(""), mFitFunc(NULL), mFitFuncDer(NULL)
193, mUFNVar(0), mUFNPar(0)
194{
195// Init et Set des options
196mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
197mParSave=0.; mStepSave=1.; mMinSave=1.; mMaxSave=-1.; mFixSave=0.;
198ResetOptions();
199mOptSave = mOpt;
200
201// enregistrement des ordres de fit
202string kw, usage, grp = "Fitting";
203
204kw = "helpfit";
205usage="";
206usage+="\n////////////////////////////////////////////////////////////////////// ";
207usage+="\n| --------------- Fit Lineaire a 1 et 2 dimensions --------------- ";
208usage+="\n| Syntaxe: ";
209usage+="\n| fitlin nom pnn [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...] ";
210usage+="\n| avec: ";
211usage+="\n| nom : cf commentaire ordre ''fit'' ";
212usage+="\n| pnn : fit polynome degre nn avec classe Poly (lineaire) 1D ou 2D ";
213usage+="\n| o:aa,...,bb : cf commentaires ordre ''fit'' ";
214usage+="\n| ";
215usage+="\n////////////////////////////////////////////////////////////////////// ";
216usage+="\n| --------------- Fit d'objets a 1 et 2 dimensions --------------- ";
217usage+="\n| avec interface d'aide graphique ";
218usage+="\n| Syntaxe: ";
219usage+="\n| fitw nom func ";
220usage+="\n| [p:p1,p2,...,pn s:s1,s2,...,sn m:m1,m2,...,mn M:M1,mM2,...,Mn f:f1,...,fn] ";
221usage+="\n| [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...] ";
222usage+="\n| cf commentaire ordre ''fit'' ";
223usage+="\n| ";
224usage+="\n////////////////////////////////////////////////////////////////////// ";
225usage+="\n| --------------- Fit d'objets a 1 et 2 dimensions --------------- ";
226usage+="\n| Syntaxe: ";
227usage+="\n| fit nom func ";
228usage+="\n| [p:p1,p2,...,pn s:s1,s2,...,sn m:m1,m2,...,mn M:M1,mM2,...,Mn f:f1,...,fn] ";
229usage+="\n| [o:.aa,bb,cc, o;dd,ee,ff o:gg,hh,jj,kk etc...] ";
230usage+="\n|----- OBJET ----- ";
231usage+="\n| nom : nom de l'objet qui peut etre: ";
232usage+="\n| fit-1D: Vector,Histo1D,HProf ou GeneraFitData(1D) ";
233usage+="\n| fit-2D: Matrix,Histo2D,Image<T> ou GeneraFitData(2D) ";
234usage+="\n| ";
235usage+="\n|----- FUNCTION ----- ";
236usage+="\n| func : pnn : fit polynome degre nn avec GeneralFit (non-lineaire) 1D ou 2D ";
237usage+="\n| : gnn : fit gaussienne (hauteur) + polynome de degre nn 1D ";
238usage+="\n| : g : fit gaussienne (hauteur) 1D ";
239usage+="\n| : enn : fit exponentielle + polynome de degre nn 1D ";
240usage+="\n| : e : fit exponentielle 1D ";
241usage+="\n| : Gnn : fit gaussienne (volume) + polynome de degre nn 1D ";
242usage+="\n| : G : fit gaussienne (volume) 1D ";
243usage+="\n| : : fit gaussienne+fond (volume) 2D ";
244usage+="\n| : Gi : fit gaussienne+fond integree (volume) 2D ";
245usage+="\n| : d : fit DL de gaussienne+fond (volume) 2D ";
246usage+="\n| : di : fit DL de gaussienne+fond integree (volume) 2D ";
247usage+="\n| : D : fit DL de gaussienne+fond avec coeff variable p6 (volume) 2D ";
248usage+="\n| : Di : fit DL de gaussienne+fond integree avec coeff variable p6 (volume) 2D ";
249usage+="\n| : M : fit Moffat+fond (expos=p6) (volume) 2D ";
250usage+="\n| : Mi : fit Moffat+fond integree (expos=p6) (volume) 2D ";
251usage+="\n| : Autre : fonction definie par l'utilisateur (cf commentaires ci apres) ";
252usage+="\n| ";
253usage+="\n|----- INIT PARAMETRES ET ETAT DU FIT ----- ";
254usage+="\n| p : p1,...,pn : valeur d'initialisation des parametres (def=0) ";
255usage+="\n| s : s1,...,sn : valeur des steps de depart (def=1) ";
256usage+="\n| m : m1,...,mn : valeur des minima (def=1) ";
257usage+="\n| M : M1,...,Mn : valeur des maxima (def=-1) (max<=min : pas de limite) ";
258usage+="\n| f : f1,...,fn : si >=1 parametre fixe sinon libre (def=0) ";
259usage+="\n| - Remarque: si pi,si,mi ou Mi = '!' la valeur correspondante n'est pas changee ";
260usage+="\n| ";
261usage+="\n|----- OPTIONS ----- ";
262usage+="\n| o : options ''o:Eaa.b,eaa.b,f,r,caa.b,Xaa.b'' ";
263usage+="\n| F : initialisation a partir des resultats et de l'etat du fit precedent ";
264usage+="\n| (option non prioritaire sur les definitions p:,s:,m:,M:,f:o:) ";
265usage+="\n| a : initialisation automatique des parametres du fit (si disponible) ";
266usage+="\n| (option non prioritaire sur les definitions p:,s:,m:,M:,f:o: ";
267usage+="\n| et sur l'initialisation pa l'option ''F'') ";
268usage+="\n| f : generation d'un Objet identique contenant la fonction fittee ";
269usage+="\n| r : generation d'un Objet identique contenant les residus ";
270usage+="\n| Xaa.b : aa.b valeur du DXi2 d'arret (def=1.e-3) ";
271usage+="\n| Naa : aa nombre maximum d'iterations (def=100) ";
272usage+="\n| la : niveau ''a'' de print: a=niveau de print Fit1/2D ";
273usage+="\n| da : niveau ''a'' de debug: a=niveau de GeneralFit ";
274usage+="\n| Ii1/i2 : numeros des bins X de l'histos utilises pour le fit [i1,i2] ";
275usage+="\n|2D Jj1/j2 : numeros des bins Y de l'histos utilises pour le fit [j1,j2] ";
276usage+="\n| - L'erreur est celle associee a l'objet (si elle existe), ";
277usage+="\n| elle est mise a 1 sinon, sauf si E... ou e... est precise: ";
278usage+="\n| Eaa.b : si |val|>=1 erreur = aa.b*sqrt(|val|) ";
279usage+="\n| si |val|<1 erreur = aa.b ";
280usage+="\n| si aa.b <=0 alors aa.b=1.0 ";
281usage+="\n| E seul est equivalent a E1.0 ";
282usage+="\n| eaa.b : erreur = aa.b ";
283usage+="\n| si aa.b <=0 alors aa.b=1.0 ";
284usage+="\n| e seul est equivalent a e1.0 ";
285usage+="\n| xaa.b : demande de centrage: on fit x-aa.b au lieu de x) ";
286usage+="\n| x : demande de centrage: on fit x-xc au lieu de x ";
287usage+="\n| avec xc=abscisse du milieu de l'histogramme ";
288usage+="\n| Actif pour exp+poly 1D, poly 1D ";
289usage+="\n| Pour gauss+poly 1D, xc est le centre de la gaussienne. ";
290usage+="\n|2D yaa.b et y : idem ''xaa.b et x'' mais pour y ";
291usage+="\n| ";
292usage+="\n////////////////////////////////////////////////////////////////////// ";
293usage+="\n| --------------- Def des fcts de fit par l'utilisateur --------------- ";
294usage+="\n| Pour definir une fonction parametree utilisateur: ";
295usage+="\n| Un exemple de fichier contenant des fonctions a fitter: PIext/userfitfunex.c ";
296usage+="\n| > crfitfil nom_fichier nom_fonction nvar npar ";
297usage+="\n| pour creer un squelette de fichier C contenant la fonction a fitter ";
298usage+="\n| nom_fichier : nom du fichier ou est le code de la fonction ";
299usage+="\n| c'est un fichier C de nom blabla.c ";
300usage+="\n| nom_fonction : nom de la fonction ex: gauss2 ";
301usage+="\n| nvar : nombre de variables x[] (1D=1, 2D=2, etc...) ";
302usage+="\n| npar : nombre de parametre p[] ";
303usage+="\n| - La fonction est donc F(x[];p[]), soit par ex: ";
304usage+="\n| double gauss2(double const* x,double const* p) ";
305usage+="\n| - L'utilisateur a la possibilite de donner egalement la fonction ";
306usage+="\n| retournant les derivees de la fonction par rapport aux parametres ";
307usage+="\n| double gauss2_der(double const* x,double const* p,double* dp) ";
308usage+="\n| ou dp[i] = dF(x[];p[])/dp[i] ";
309usage+="\n| Par convention, cette fonction a le meme nom suivi de ''_der'' ";
310usage+="\n| > crfitfun nom_fichier nom_fonction nvar npar ";
311usage+="\n| Permet de compiler le fichier ''nom_fichier'' et de linker ";
312usage+="\n| la fonction ''nom_fonction'' ";
313usage+="\n| - Un meme fichier peut contenir plusieurs fonctions a fitter ";
314usage+="\n| ";
315usage+="\n////////////////////////////////////////////////////////////////////// ";
316piac->RegisterHelp(kw,usage,grp);
317
318kw = "fit";
319usage = "Fitting function to DataObjects (Histo, Histo2D, Vector, ...)";
320usage += "\n Usage: fit nomobj func [Options]";
321usage += "\n [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]";
322usage += "\n Related commands: fitw fitlin crfitfun crfitfil";
323piac->RegisterCommand(kw,usage,this,grp);
324
325kw = "fitw";
326usage = "Fitting function to DataObjects with Interactive Window Help";
327usage += "\n Usage: fitw nomobj func [Options]";
328usage += "\n [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]";
329usage += "\n Related commands: fit fitlin crfitfun crfitfil";
330piac->RegisterCommand(kw,usage,this,grp);
331
332kw = "fitlin";
333usage = "Linear Fitting of Polynoms to DataObjects";
334usage += "\n Usage: fitlin nomobj func [o:... o:...]";
335usage += "\n Related commands: fit fitw crfitfun crfitfil";
336piac->RegisterCommand(kw,usage,this,grp);
337
338kw = "crfitfun";
339usage = "Creation et link de function utilisateur pour le fit";
340usage += "\n Usage: crfitfun file func nvar npar";
341usage += "\n Related commands: fit fitw fitlin crfitfil";
342piac->RegisterCommand(kw,usage,this,grp);
343
344kw = "crfitfil";
345usage = "Creation de fichier C pour function utilisateur du fit";
346usage += "\n Usage: crfitfil file func nvar npar";
347usage += "\n Related commands: fit fitw fitlin crfitfun";
348piac->RegisterCommand(kw,usage,this,grp);
349}
350
351/* --Methode-- */
352PIAFitter::~PIAFitter()
353{
354if(mFit != NULL) {delete mFit; mFit=NULL;}
355if(mFunction != NULL) {delete mFunction; mFunction=NULL;}
356if(mFunc != NULL) {delete mFunc; mFunc=NULL;}
357if(mGData != NULL) {delete mGData; mGData=NULL;}
358if(FWindFit != NULL) {delete FWindFit; FWindFit=NULL;}
359if (mDlUFunc != NULL) {delete mDlUFunc; mDlUFunc=NULL;}
360}
361
362/* --Methode-- */
363void PIAFitter::ResetDPointer(void)
364// Reset des pointeurs sur les donnees
365{
366mNObj = ""; mObj=NULL;
367mV=NULL; mH=NULL; mM=NULL; mH2=NULL; mIm=NULL; mG=NULL;
368}
369
370/* --Methode-- */
371void PIAFitter::ResetOptions(void)
372// Reset des options
373{
374mOpt.okres = mOpt.okfun = false;
375mOpt.polcx = mOpt.polcy = 0;
376mOpt.xc = mOpt.yc = 0.;
377mOpt.stc2 = 1.e-3;
378mOpt.nstep = 100;
379mOpt.err_e = mOpt.err_E = -1.;
380mOpt.lp = 1; mOpt.lpg = 0;
381mOpt.i1 = mOpt.j1 = mOpt.i2 = mOpt.j2 = -1;
382mOpt.fromlastfit = false;
383mOpt.autoinifit = false;
384}
385
386/* --Methode-- */
387void PIAFitter::DecodeOptions(string opt)
388// Decodage des options a partir des chaines de caracteres
389{
390ResetOptions();
391if(opt.length()<=0) return;
392
393size_t p,q; string dum;
394opt = "," + opt + ",";
395
396if(strstr(opt.c_str(),",F,")) // init options/param a partir du fit precedent
397 {mOpt = mOptSave; mOpt.fromlastfit = true;}
398if(strstr(opt.c_str(),",a,")) mOpt.autoinifit = true; // init auto par fit
399if(strstr(opt.c_str(),",r,")) mOpt.okres = true; // residus
400if(strstr(opt.c_str(),",f,")) mOpt.okfun = true; // fonction fittee
401if(strstr(opt.c_str(),",x")) { // Demande de centrage (fit X=x-xc)
402 mOpt.polcx = 2; // Le centrage est calcule automatiquement
403 p = opt.find(",x"); q = opt.find_first_of(',',p+1);
404 dum = opt.substr(p,q-p);
405 if(dum.length()>2) {
406 sscanf(dum.c_str(),",x%lf",&mOpt.xc);
407 mOpt.polcx = 1; // Le centrage est fixe par la valeur lue
408 }
409}
410if(strstr(opt.c_str(),",y")) { // Demande de centrage (fit Y=y-yc)
411 mOpt.polcy = 2; // Le centrage est calcule automatiquement
412 p = opt.find(",y"); q = opt.find_first_of(',',p+1);
413 dum = opt.substr(p,q-p);
414 if(dum.length()>2) {
415 sscanf(dum.c_str(),",y%lf",&mOpt.yc);
416 mOpt.polcy = 1; // Le centrage est fixe par la valeur lue
417 }
418}
419if(strstr(opt.c_str(),",E")) { // Erreurs imposees a "sqrt(val)" ou "aa.b*sqrt(val)"
420 p = opt.find(",E"); q = opt.find_first_of(',',p+1);
421 dum = opt.substr(p,q-p);
422 if(dum.length()>2) sscanf(dum.c_str(),",E%lf",&mOpt.err_E);
423 if(mOpt.err_E<=0.) mOpt.err_E = 1.; //mOpt.err_e=-1.;
424}
425if(strstr(opt.c_str(),",e")) { // Erreurs imposees a "1" ou "aa.b"
426 p = opt.find(",e"); q = opt.find_first_of(',',p+1);
427 dum = opt.substr(p,q-p);
428 if(dum.length()>2) sscanf(dum.c_str(),",e%lf",&mOpt.err_e);
429 if(mOpt.err_e<=0.) mOpt.err_e = 1.; //mOpt.err_E=-1.;
430}
431if(strstr(opt.c_str(),",X")) { // Valeur du StopChi2
432 p = opt.find(",X"); q = opt.find_first_of(',',p+1);
433 dum = opt.substr(p,q-p);
434 if(dum.length()>2) sscanf(dum.c_str(),",X%lf",&mOpt.stc2);
435}
436if(strstr(opt.c_str(),",N")) { // Nombre maximum d'iterations
437 p = opt.find(",N"); q = opt.find_first_of(',',p+1);
438 dum = opt.substr(p,q-p);
439 if(dum.length()>2) sscanf(dum.c_str(),",N%d",&mOpt.nstep);
440}
441if(strstr(opt.c_str(),",l")) { // niveau de print
442 p = opt.find(",l"); q = opt.find_first_of(',',p+1);
443 dum = opt.substr(p,q-p);
444 if(dum.length()>2) sscanf(dum.c_str(),",l%d",&mOpt.lp);
445}
446if(strstr(opt.c_str(),",d")) { // niveau de debug du fit
447 p = opt.find(",d"); q = opt.find_first_of(',',p+1);
448 dum = opt.substr(p,q-p);
449 if(dum.length()>2) sscanf(dum.c_str(),",d%d",&mOpt.lpg);
450}
451if(strstr(opt.c_str(),",I")) { // intervalle de fit selon X
452 p = opt.find(",I"); q = opt.find_first_of(',',p+1);
453 dum = opt.substr(p,q-p);
454 if(dum.length()>2) sscanf(dum.c_str(),",I%d/%d",&mOpt.i1,&mOpt.i2);
455}
456if(strstr(opt.c_str(),",J")) { // intervalle de fit selon Y
457 p = opt.find(",J"); q = opt.find_first_of(',',p+1);
458 dum = opt.substr(p,q-p);
459 if(dum.length()>2) sscanf(dum.c_str(),",J%d/%d",&mOpt.j1,&mOpt.j2);
460}
461
462return;
463}
464
465/* --Methode-- */
466void PIAFitter::DecodeObject(string obj)
467// Decodage de l'objet a fitter
468{
469NamedObjMgr omg;
470AnyDataObj* mobj = omg.GetObj(obj);
471if(mobj == NULL)
472 {cout<<"PIAFitter::DecodeObject Error , Pas d'objet de nom "<<obj<<endl;
473 return;}
474
475// Reset
476ResetDPointer();
477mNVar = mNBinX = mNBinY = mNData = 0;
478
479mObj = mobj; mNObj = obj;
480string ctyp = typeid(*mObj).name();
481
482// 1D
483if (typeid(*mObj) == typeid(Vector)) {
484 mNVar = 1;
485 mV = (Vector*) mObj; mNBinX = mV->NElts(); mNBinY = 1;
486 }
487else if ( (typeid(*mObj) == typeid(HProf)) || (typeid(*mObj) == typeid(Histo)) ) {
488 mNVar = 1;
489 mH = (Histo*) mObj; mNBinX = mH->NBins(); mNBinY = 1;
490 }
491// 2D
492else if (typeid(*mObj) == typeid(Matrix)) {
493 mNVar = 2;
494 mM = (Matrix*) mObj; mNBinX = mM->NCol(); mNBinY =mM->NRows();
495 }
496else if (typeid(*mObj) == typeid(Histo2D)) {
497 mNVar = 2;
498 mH2 = (Histo2D*) mObj; mNBinX = mH2->NBinX(); mNBinY = mH2->NBinY();
499 }
500else if (typeid(*mObj) == typeid(GeneralFitData)) {
501 mG = (GeneralFitData*) mObj; mNBinX = mG->NData(); mNBinY = 1;
502 if( mG->NVar()==1) mNVar = 1;
503 else if(mG->NVar()==2) mNVar = 2;
504 else
505 {cout<<"PIAFitter::DecodeObject Error: \n"
506 <<" GeneralFitData ne peut avoir que 1 ou 2 variables d'abscisse: "
507 <<((GeneralFitData*) mObj)->NVar()<<endl;
508 return;}
509 }
510#ifdef SANS_EVOLPLANCK
511else if (dynamic_cast<RzImage*>(mObj)) {
512 mNVar = 2;
513 mIm = (RzImage*) mObj; mNBinX = mIm->XSize(); mNBinY = mIm->YSize();
514 }
515#else
516//CMV_A_FAIRE else if (dynamic_cast<RzImage*>(mObj)) {
517//CMV_A_FAIRE mNVar = 2;
518//CMV_A_FAIRE mIm = (RzImage*) mObj; mNBinX = mIm->XSize(); mNBinY = mIm->YSize();
519//CMV_A_FAIRE }
520#endif
521else {
522 cout<<"PIAFitter::DecodeObject Error , Objet n'est pas un "
523 <<"Histo1D/HProf/Vector/Histo2D/Image/Matrix/GeneralFitData "<<ctyp<<endl;
524 return;
525 }
526
527mNData = mNBinX*mNBinY;
528if(mNData<=0) {
529 cout<<"L'objet a "<<mNBinX<<","<<mNBinY<<" bins ("<<mNData<<")"<<endl;
530 return;
531}
532
533return;
534}
535
536/* --Methode-- */
537void PIAFitter::CheckOptions(void)
538// Check consistence des choix des options
539// ATTENTION: cette methode a besoin que l'objet a fitter soit deja lu
540{
541if(mOpt.lp<0) mOpt.lp = 0;
542if(mOpt.lpg<0) mOpt.lpg = 0;
543mOpt.i1 = (mOpt.i1<0||mOpt.i1>=mNBinX)? 0: mOpt.i1;
544mOpt.i2 = (mOpt.i2<0||mOpt.i2>=mNBinX||mOpt.i2<=mOpt.i1)? mNBinX-1: mOpt.i2;
545if(mNVar>=2) {
546 mOpt.j1 = (mOpt.j1<0||mOpt.j1>=mNBinY)? 0: mOpt.j1;
547 mOpt.j2 = (mOpt.j2<0||mOpt.j2>=mNBinY||mOpt.j2<=mOpt.j1)? mNBinY-1: mOpt.j2;
548} else mOpt.j2 = mOpt.j1 = 0;
549if(mOpt.polcx==2) {
550 if(mV||mM) mOpt.xc = (mOpt.i2-mOpt.i1+1)/2.;
551 else if(mH) mOpt.xc = (mH->XMin()+mH->XMax())/2.;
552 else if(mH2) mOpt.xc = (mH2->XMin()+mH2->XMax())/2.;
553 else if(mG) {double mini,maxi; mG->GetMnMx(2,mini,maxi); mOpt.xc=(mini+maxi)/2.;}
554 else if(mIm) {mOpt.xc = mIm->XOrg() + mIm->XPixSize()*(mOpt.i2-mOpt.i1+1)/2.;}
555}
556if(mOpt.polcy==2 && mNVar>=2) {
557 if(mM) mOpt.yc = (mOpt.j2-mOpt.j1+1)/2.;
558 if(mH2) mOpt.yc = (mH2->YMin()+mH2->YMax())/2.;
559 if(mG) {double mini,maxi; mG->GetMnMx(12,mini,maxi); mOpt.yc=(mini+maxi)/2.;}
560 if(mIm) {mOpt.yc = mIm->YOrg() + mIm->YPixSize()*(mOpt.j2-mOpt.j1+1)/2.;}
561}
562if(mOpt.err_e>0.) mOpt.err_E=-1.;
563if(mOpt.err_E>0.) mOpt.err_e =-1.;
564if(mOpt.stc2<=0.) mOpt.stc2 = 1.e-3;
565if(mOpt.nstep<2) mOpt.nstep = 100;
566
567return;
568}
569
570/* --Methode-- */
571void PIAFitter::PrintOptions(void)
572// Print des options
573{
574cout<<"Fit["<<mNBinX<<","<<mNBinY<<"] ("<<mNData<<") dim="<<mNVar<<":"
575 <<" Int=["<<mOpt.i1<<","<<mOpt.i2<<"],["<<mOpt.j1<<","<<mOpt.j2<<"]"<<endl
576 <<" Cent="<<mOpt.polcx<<","<<mOpt.polcy<<","<<mOpt.xc<<"+x"<<","<<mOpt.yc<<"+y"
577 <<" TypE="<<mOpt.err_e<<","<<mOpt.err_E
578 <<" StpX2="<<mOpt.stc2<<" Nstep="<<mOpt.nstep<<endl
579 <<" Init.LFit="<<mOpt.fromlastfit<<" AutoIni.Fit="<<mOpt.autoinifit
580 <<" lp,lpg="<<mOpt.lp<<","<<mOpt.lpg<<endl;
581return;
582}
583
584/* --Methode-- */
585void PIAFitter::FillGData(void)
586// Fill generalfitdata pour le fit
587// ATTENTION: cette methode a besoin que les options soient decodees et checkees
588{
589if(mNData<=0) return;
590if(mGData!=NULL) {delete mGData; mGData=NULL;}
591mGData = new GeneralFitData(mNVar,mNData,0);
592
593for(int i=mOpt.i1;i<=mOpt.i2;i++) for(int j=mOpt.j1;j<=mOpt.j2;j++) {
594 double x,y=0.,f,e;
595
596 if(mV)
597 {x= (double) i; f=(*mV)(i); e=1.;}
598 else if(mH)
599 {x=mH->BinCenter(i); f=(*mH)(i); e=(mH->HasErrors())?mH->Error(i):1.;}
600 else if(mM)
601 {x=(double) i; y=(double) j; f=(*mM)(j,i); e=1.;}
602 else if(mH2)
603 {r_8 xf,yf; mH2->BinCenter(i,j,xf,yf); x=xf; y=yf;
604 f=(*mH2)(i,j); e=(mH2->HasErrors())?mH2->Error(i,j):1.;}
605 else if(mIm)
606 {x=mIm->XOrg()+(i+0.5)*mIm->XPixSize(); y=mIm->YOrg()+(j+0.5)*mIm->YPixSize();
607 f=mIm->DValue(i,j); e=1.;}
608 else if(mG&&mNVar==1) {x= mG->X(i); f=mG->Val(i); e=mG->EVal(i);}
609 else if(mG&&mNVar==2) {x= mG->X(i); y=mG->Y(i); f=mG->Val(i); e=mG->EVal(i);}
610 else x=y=f=e=0.;
611
612 // Gestion des erreurs a utiliser
613 if(mOpt.err_e>0.) e=mOpt.err_e;
614 else if(mOpt.err_E>0.) {e=(f<-1.||f>1.)?mOpt.err_E*sqrt(fabs(f)):mOpt.err_E;}
615
616 // Remplissage de generalfitdata
617 if(ApplyXYcentering) { // pour fitlin
618 if(mOpt.polcx>0) x-=mOpt.xc;
619 if(mNVar>=2 && mOpt.polcy>0) y-=mOpt.yc;
620 }
621 if(mNVar==1) mGData->AddData1(x,f,e);
622 else if(mNVar==2) mGData->AddData2(x,y,f,e);
623}
624
625if(mGData->NData()<=0) {
626 cout<<"Pas de donnees dans GeneralFitData: "<<mGData->NData()<<endl;
627 return;
628}
629
630if(mOpt.lp>2) {
631 mGData->PrintStatus();
632 mGData->PrintData(0);
633 mGData->PrintData(mGData->NData()-1);
634 double mini,maxi;
635 mGData->GetMnMx(2,mini,maxi);
636 cout<<"Limites ...x0 = ["<<mini<<","<<maxi<<"]"<<endl;
637 if(mNVar==2)
638 {mGData->GetMnMx(12,mini,maxi);
639 cout<<"Limites ...x1 = ["<<mini<<","<<maxi<<"]"<<endl;}
640 mGData->GetMnMx(0,mini,maxi);
641 cout<<"Limites ...y = ["<<mini<<","<<maxi<<"]"<<endl;
642 mGData->GetMnMx(1,mini,maxi);
643 cout<<"Limites ...ey = ["<<mini<<","<<maxi<<"]"<<endl;
644}
645
646return;
647}
648
649/* --Methode-- */
650void PIAFitter::DecodeFunction(string func)
651// Fonction a fitter
652// ATTENTION: cette methode a besoin que les donnees soient lues
653// (pour mopt.xc,yc et mNVar)
654{
655if(func.length()<=0)
656 {cout<<"PIAFitter::DecodeFunction Donnez un nom de fonction a fitter."<<endl;
657 return;}
658if(mFunction!=NULL) {delete mFunction; mFunction=NULL;}
659if(mFunc !=NULL) {delete mFunc; mFunc=NULL;}
660mNPar=0; mFName = "";
661
662if(func == mNameFitFunc) { // Fonction utilisateur
663 mFunc = new GeneralFunc(mUFNVar,mUFNPar,mFitFunc,mFitFuncDer);
664 mNPar = mUFNPar;
665 cout<<"Fonction utilisateur nvar="<<mUFNVar<<", npar="<<mNPar<<endl;
666} else if(func[0]=='p' && mNVar==1) { //polynome
667 int degre=0; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
668 cout<<"Fit polynome 1D de degre "<<degre<<endl;
669 Polyn1D* myf = new Polyn1D(degre,((mOpt.polcx>0)?mOpt.xc:0.));
670 mFunction = myf; mNPar = mFunction->NPar(); mFName = func;
671
672} else if(func[0]=='e' && mNVar==1) { //exponentielle
673 int degre=-1; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
674 cout<<"Fit d'exponentielle+polynome 1D de degre "<<degre<<endl;
675 Exp1DPol* myf;
676 if(degre>=0) myf = new Exp1DPol((unsigned int)degre,((mOpt.polcx>0)?mOpt.xc:0.));
677 else myf = new Exp1DPol(((mOpt.polcx>0)?mOpt.xc:0.));
678 mFunction = myf; mNPar = mFunction->NPar(); mFName = func;
679
680} else if(func[0]=='g' && mNVar==1) { //gaussienne en hauteur
681 int degre=-1; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
682 cout<<"Fit de Gaussienne_en_hauteur+polynome 1D de degre "<<degre<<endl;
683 Gauss1DPol* myf;
684 if(degre>=0) myf = new Gauss1DPol((unsigned int)degre,((mOpt.polcx)?true:false));
685 else { bool bfg = (mOpt.polcx)?true:false; myf = new Gauss1DPol(bfg); }
686 mFunction = myf; mNPar = mFunction->NPar(); mFName = func;
687
688} else if(func[0]=='G' && mNVar==1) { //gaussienne en volume
689 int degre=-1; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
690 cout<<"Fit de Gaussienne_en_volume+polynome 1D de degre "<<degre<<endl;
691 GaussN1DPol* myf;
692 if(degre>=0) myf = new GaussN1DPol((unsigned int)degre,((mOpt.polcx)?true:false));
693 else { bool bfg = (mOpt.polcx)?true:false; myf = new GaussN1DPol(bfg); }
694 mFunction = myf; mNPar = mFunction->NPar(); mFName = func;
695
696} else if(func[0]=='p' && mNVar==2) { //polynome 2D
697 int degre=0; if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
698 cout<<"Fit polynome 2D de degre "<<degre<<endl;
699 Polyn2D* myf =
700 new Polyn2D(degre,((mOpt.polcx>0)?mOpt.xc:0.),((mOpt.polcy>0)?mOpt.yc:0.));
701 mFunction = myf; mNPar = mFunction->NPar(); mFName = func;
702
703} else if(func[0]=='G' && mNVar==2) { //gaussienne+fond en volume
704 int integ=0; if(func.length()>1) if(func[1]=='i') integ=1;
705 cout<<"Fit de Gaussienne+Fond 2D integ="<<integ<<endl;
706 if(integ) {GauRhInt2D* myf = new GauRhInt2D; mFunction = myf;}
707 else {GauRho2D* myf = new GauRho2D; mFunction = myf;}
708 mNPar = mFunction->NPar(); mFName = func;
709
710} else if(func[0]=='d' && mNVar==2) { //DL gaussienne+fond en volume
711 int integ=0; if(func.length()>1) if(func[1]=='i') integ=1;
712 cout<<"Fit de DL de Gaussienne+Fond 2D integ="<<integ<<endl;
713 if(integ) {GdlRhInt2D* myf = new GdlRhInt2D; mFunction = myf;}
714 else {GdlRho2D* myf = new GdlRho2D; mFunction = myf;}
715 mNPar = mFunction->NPar(); mFName = func;
716
717} else if(func[0]=='D' && mNVar==2) { //DL gaussienne+fond avec coeff variable p6 en volume
718 int integ=0; if(func.length()>1) if(func[1]=='i') integ=1;
719 cout<<"Fit de DL de Gaussienne+Fond avec coeff variable (p6) 2D integ="<<integ<<endl;
720 if(integ) {Gdl1RhInt2D* myf = new Gdl1RhInt2D; mFunction = myf;}
721 else {Gdl1Rho2D* myf = new Gdl1Rho2D; mFunction = myf;}
722 mNPar = mFunction->NPar(); mFName = func;
723
724} else if(func[0]=='M' && mNVar==2) { //Moffat+fond (volume)
725 int integ=0; if(func.length()>1) if(func[1]=='i') integ=1;
726 cout<<"Fit de Moffat+Fond (expos=p6) 2D integ="<<integ<<endl;
727 if(integ) {MofRhInt2D* myf = new MofRhInt2D; mFunction = myf;}
728 else {MofRho2D* myf = new MofRho2D; mFunction = myf;}
729 mNPar = mFunction->NPar(); mFName = func;
730
731} else {
732 cout<<"Fonction "<<func<<" inconnue pour la dim "<<mNVar<<endl;
733 return;
734}
735
736return;
737}
738
739/* --Methode-- */
740void PIAFitter::ReSetParam(void)
741// Pour (re)dimensionner les tableaux des parametres
742// ATTENTION: cette methode a besoin que la fonction soit connue
743{
744if(mNPar==0) return;
745mPar.Realloc(mNPar); mPar=0.;
746mStep.Realloc(mNPar); mStep=1.;
747mMin.Realloc(mNPar); mMin=1.;
748mMax.Realloc(mNPar); mMax=-1.;
749mFix.Realloc(mNPar); mFix=0.;
750return;
751}
752
753/* --Methode-- */
754void PIAFitter::AutoIniFit(void)
755// Initialisation automatique des parametres du fit
756{
757if(mFunc) { // Fonction utilisateur
758 cout<<"PIAFitter::AutoIniFit : Fonction utilisateur, pas d'init auto"<<endl;
759} else if(mFName[0]=='p' && mNVar==1) { //poly 1D
760 cout<<"PIAFitter::AutoIniFit : Auto Init type polynome 1D"<<endl;
761 IniFitP1ou2D();
762
763} else if(mFName[0]=='g' && mNVar==1) { //gauss en haut 1D
764 cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_haut+pol 1D"<<endl;
765 IniFitGhP1D();
766
767} else if(mFName[0]=='G' && mNVar==1) { //gauss en vol 1D
768 cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_vol+pol 1D"<<endl;
769 IniFitGvP1D();
770
771} else if(mFName[0]=='p' && mNVar==2) { //poly 2D
772 cout<<"PIAFitter::AutoIniFit : Auto Init type polynome 2D"<<endl;
773 IniFitP1ou2D();
774
775} else if((mFName[0]=='G'||mFName[0]=='d') && mNVar==2) { //gauss+fond en vol
776 cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_vol+fond 2D"<<endl;
777 IniFitGv2D();
778
779} else if(mFName[0]=='D' && mNVar==2) { //DLgauss+fond en vol 2D
780 cout<<"PIAFitter::AutoIniFit : Auto Init type DLgauss_vol+fond 2D"<<endl;
781 IniFitGv2D();
782 Vector dum(mNPar);
783 dum = mPar; mPar(6)=0.5; mPar(7)=dum(6);
784 dum = mStep; mStep(6)=0.05; mStep(7)=dum(6);
785 dum = mMin; mMin(6)=0.; mMin(7)=dum(6);
786 dum = mMax; mMax(6)=10.; mMax(7)=dum(6);
787 dum = mFix; mFix(6)=1.1; mFix(7)=dum(6);
788
789} else if(mFName[0]=='M' && mNVar==2) { //Moffat+fond en vol 2D
790 cout<<"PIAFitter::AutoIniFit : Auto Init type moffat_vol+fond 2D"<<endl;
791 IniFitGv2D();
792 Vector dum(mNPar);
793 dum = mPar; mPar(6)=2.5; mPar(7)=dum(6);
794 dum = mStep; mStep(6)=0.05; mStep(7)=dum(6);
795 dum = mMin; mMin(6)=1.2; mMin(7)=dum(6);
796 dum = mMax; mMax(6)=10.; mMax(7)=dum(6);
797 dum = mFix; mFix(6)=1.1; mFix(7)=dum(6);
798
799} else {
800 cout<<"PIAFitter::AutoIniFit : Fonction "<<mFName
801 <<", \n pas d'init auto en dim "<<mNVar<<endl;
802 return;
803}
804
805return;
806}
807
808/* --Methode-- */
809void PIAFitter::InitParFromLastFit(void)
810// Initialisation du fit a partir des resultats du fit precedent
811// (Les arguments de la ligne de commande sont prioritaires)
812// ATTENTION: cette methode a besoin que la fonction soit connue
813{
814if(mNPar==0) return;
815int n = mParSave.NElts();
816if(n==0) return;
817
818for(int i=0;i<mNPar;i++) {
819 if(i>=n) break;
820 mPar(i) = mParSave(i);
821 mStep(i) = mStepSave(i);
822 mMin(i) = mMinSave(i);
823 mMax(i) = mMaxSave(i);
824 mFix(i) = mFixSave(i);
825}
826
827return;
828}
829
830/* --Methode-- */
831void PIAFitter::DecodeParam(string par,string step,string min,string max,string fix)
832// Decodage des parametres
833// ATTENTION: cette methode a besoin que la fonction soit connue
834// et que les options soient decodees
835{
836if(mNPar==0) return;
837ReSetParam();
838if(mOpt.autoinifit) AutoIniFit();
839if(mOpt.fromlastfit) InitParFromLastFit();
840Vector* v=NULL; string* s=NULL;
841for(int j=0;j<5;j++) {
842 if(j==0) {v=&mPar; s=&par;} else if(j==1) {v=&mStep; s=&step;}
843 else if(j==2) {v=&mMin; s=&min;} else if(j==3) {v=&mMax; s=&max;}
844 else if(j==4) {v=&mFix; s=&fix;}
845 if(s->length()>0) *s += ",";
846 for(int i=0;i<mNPar;i++) {
847 if(s->length()<=0) break;
848 if((*s)[0]!='!') sscanf(s->c_str(),"%lf",&(*v)(i));
849 size_t p = s->find_first_of(',') + 1;
850 if(p>=s->length()) *s = ""; else *s = s->substr(p);
851 }
852}
853return;
854}
855
856/* --Methode-- */
857int PIAFitter::DoFit(void)
858// Fit avec generalfit
859{
860if(mNPar==0) return -1000;
861if((mFunction==NULL && mFunc==NULL) || mGData==NULL) {
862 cout<<"PIAFitter::DoFit error: mFunction,mFunc="<<mFunction<<","<<mFunc
863 <<" mGData="<<mGData<<endl;
864 return -1001;
865}
866// classe GeneraFit
867if(mFit != NULL) {delete mFit; mFit= NULL;}
868if(mFunction) mFit = new GeneralFit(mFunction);
869 else if(mFunc) mFit = new GeneralFit(mFunc);
870if(!mFit) return(-1002);
871
872// Options et Set de GeneralFit
873mFit->SetDebug(mOpt.lpg);
874mFit->SetData(mGData);
875mFit->SetStopChi2(mOpt.stc2);
876mFit->SetMaxStep(mOpt.nstep);
877{for(int i=0;i<mNPar;i++) {
878 char str[10];
879 sprintf(str,"P%d",i);
880 mFit->SetParam(i,str,mPar(i),mStep(i),mMin(i),mMax(i));
881 if((int)(mFix(i)+0.001)>=1.) mFit->SetFix(i);
882}}
883if(mOpt.lp>1) mFit->PrintFit();
884
885// Fit
886#ifdef SANS_EVOLPLANCK
887mParSave = mPar; mStepSave = mStep; mFixSave = mFix;
888mMinSave = mMin; mMaxSave = mMax;
889#else
890mParSave.Clone(mPar); mStepSave.Clone(mStep); mFixSave.Clone(mFix);
891mMinSave.Clone(mMin); mMaxSave.Clone(mMax);
892#endif
893mOptSave = mOpt;
894double c2r = -1.;
895int rcfit = mFit->Fit();
896if(mOpt.lp>0) mFit->PrintFit();
897if(rcfit>0) {
898 c2r = mFit->GetChi2Red();
899 cout<<"C2r_Reduit = "<<c2r<<" nstep="<<mFit->GetNStep()<<" rc="<<rcfit<<endl;
900 for(int i=0;i<mNPar;i++) mParSave(i)=mFit->GetParm(i);
901} else {
902 cout<<"echec Fit, rc = "<<rcfit<<" nstep="<<mFit->GetNStep()<<endl;
903 mFit->PrintFitErr(rcfit);
904}
905
906return rcfit;
907}
908
909/* --Methode-- */
910void PIAFitter::FitFunRes(void)
911// Mise a disposition des resultats (fonction fitee et residus)
912{
913if(mFit && (mOpt.okres||mOpt.okfun)) {
914 NamedObjMgr omg;
915 string nomres = mNObj + "_res";
916 string nomfun = mNObj + "_fun";
917//####################
918#ifdef SANS_EVOLPLANCK
919//####################
920 if(mV) {
921 if(mOpt.okres) {Vector* ob = mV->FitResidus(*mFit); if(ob) omg.AddObj(ob,nomres);}
922 if(mOpt.okfun) {Vector* ob = mV->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
923 } else if(mH) {
924 if(mOpt.okres) {Histo* ob = mH->FitResidus(*mFit); if(ob) omg.AddObj(ob,nomres);}
925 if(mOpt.okfun) {Histo* ob = mH->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
926 } else if(mM) {
927 if(mOpt.okres) {Matrix* ob = mM->FitResidus(*mFit); if(ob) omg.AddObj(ob,nomres);}
928 if(mOpt.okfun) {Matrix* ob = mM->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
929 } else if(mH2) {
930 if(mOpt.okres)
931 {Histo2D* ob = mH2->FitResidus(*mFit); if(ob) omg.AddObj(ob,nomres);}
932 if(mOpt.okfun)
933 {Histo2D* ob = mH2->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
934 } else if(mIm) {
935 if(mOpt.okres) {
936 RzImage* ob = mIm->FitResidus(*mFit);
937 if(ob) {ImageR4* rob = new ImageR4(*ob); omg.AddObj(rob,nomres); delete ob;}
938 }
939 if(mOpt.okfun) {
940 RzImage* ob = mIm->FitFunction(*mFit);
941 if(ob) {ImageR4* rob = new ImageR4(*ob); omg.AddObj(rob,nomfun); delete ob;}
942 }
943 } else if(mG) {
944 if(mOpt.okres)
945 {GeneralFitData* ob = mG->FitResidus(*mFit); if(ob) omg.AddObj(ob,nomres);}
946 if(mOpt.okfun)
947 {GeneralFitData* ob = mG->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
948 }
949//####################
950#else
951//####################
952 if(mV) {
953 if(mOpt.okres) {
954 Vector* ob = new Vector(ObjectFitter::FitResidus(*mV, *mFit),true);
955 if(ob) omg.AddObj(ob,nomres);
956 }
957 if(mOpt.okfun) {
958 Vector* ob = new Vector(ObjectFitter::FitFunction(*mV, *mFit),true);
959 if(ob) omg.AddObj(ob,nomfun);
960 }
961 } else if(mH) {
962 if(mOpt.okres) {
963 Histo* ob = new Histo(ObjectFitter::FitResidus(*mH, *mFit));
964 if(ob) omg.AddObj(ob,nomres);
965 }
966 if(mOpt.okfun) {
967 Histo* ob = new Histo(ObjectFitter::FitFunction(*mH, *mFit));
968 if(ob) omg.AddObj(ob,nomfun);
969 }
970 } else if(mM) {
971 if(mOpt.okres) {
972 Matrix* ob = new Matrix(ObjectFitter::FitResidus(*mM, *mFit), true);
973 if(ob) omg.AddObj(ob,nomres);
974 }
975 if(mOpt.okfun) {
976 Matrix* ob = new Matrix(ObjectFitter::FitFunction(*mM, *mFit), true);
977 if(ob) omg.AddObj(ob,nomfun);
978 }
979 } else if(mH2) {
980 if(mOpt.okres) {
981 Histo2D* ob = new Histo2D(ObjectFitter::FitResidus(*mH2, *mFit));
982 if(ob) omg.AddObj(ob,nomres);
983 }
984 if(mOpt.okfun) {
985 Histo2D* ob = new Histo2D(ObjectFitter::FitFunction(*mH2, *mFit));
986 if(ob) omg.AddObj(ob,nomfun);
987 }
988//CMV_A_FAIRE } else if(mIm) {
989//CMV_A_FAIRE if(mOpt.okres) {
990//CMV_A_FAIRE RzImage* ob = mIm->FitResidus(*mFit);
991//CMV_A_FAIRE if(ob) {ImageR4* rob = new ImageR4(*ob); omg.AddObj(rob,nomres); delete ob;}
992//CMV _A_FAIRE }
993//CMV_A_FAIRE if(mOpt.okfun) {
994//CMV_A_FAIRE RzImage* ob = mIm->FitFunction(*mFit);
995//CMV_A_FAIRE if(ob) {ImageR4* rob = new ImageR4(*ob); omg.AddObj(rob,nomfun); delete ob;}
996//CMV_A_FAIRE }
997 } else if(mG) {
998 if(mOpt.okres) {
999 //GeneralFitData* ob = new GeneralFitData(mG->FitResidus(*mFit));
1000 GeneralFitData* ob = new GeneralFitData(ObjectFitter::FitResidus(*mG, *mFit));
1001 if(ob) omg.AddObj(ob,nomres);
1002 }
1003 if(mOpt.okfun) {
1004 //GeneralFitData* ob = new GeneralFitData(mG->FitFunction(*mFit));
1005 GeneralFitData* ob = new GeneralFitData(ObjectFitter::FitFunction(*mG, *mFit));
1006 if(ob) omg.AddObj(ob,nomfun);
1007 }
1008 }
1009//####################
1010#endif
1011//####################
1012}
1013return;
1014}
1015
1016/* --Methode-- */
1017void PIAFitter::LinFit(void)
1018// Fit lineaire de polynomes
1019{
1020if(mFName.length()<=0)
1021 {cout<<"PIAFitter::LinFit Donnez un nom de fonction a fitter."<<endl;
1022 return;}
1023if(!mGData)
1024 {cout<<"PIAFitter::LinFit pas de Data."<<endl; return;}
1025if(mFName[0]!='p')
1026 {cout<<"PIAFitter::LinFit fonction non prevue "<<mFName<<endl;
1027 return;}
1028int degre = 0;
1029if(mFName.length()>1) sscanf(mFName.c_str()+1,"%d",&degre);
1030
1031if(mNVar==1) { // Fit lineaire de polynome 1D
1032 cout<<"Fit (lineaire) 1D polynome de degre "<<degre<<endl;
1033 Poly p1(0);
1034 double c2rl = mGData->PolFit(0,p1,degre);
1035 cout<<"C2r_lineaire = "<<c2rl<<endl;
1036 if(mOpt.lp>0) cout<<p1<<endl;
1037} else if(mNVar==2) { // Fit lineaire de polynome 2D
1038 cout<<"Fit (lineaire) polynome 2D de degre "<<degre<<endl;
1039 Poly2 p2(0);
1040 double c2rl = mGData->PolFit(0,1,p2,degre);
1041 cout<<"C2r_lineaire = "<<c2rl<<endl;
1042 if(mOpt.lp>0) cout<<p2<<endl;
1043}
1044
1045return;
1046}
1047
1048/* --Methode-- */
1049void PIAFitter::LinkFitFunc(string const& fname, string const& funcname
1050 ,int nvar,int npar)
1051// Compilation et link dynamique de la fonction a fitter
1052{
1053if(nvar<1 || npar<1) return;
1054if(mDlUFunc != NULL) {delete mDlUFunc; mDlUFunc=NULL;}
1055mFitFunc = NULL; mFitFuncDer = NULL;
1056mUFNVar= mUFNPar = 0; mNameFitFunc = "";
1057
1058// On compile le fichier C, on fabrique un .so et on l'ouvre
1059mDlUFunc = PDynLinkMgr::BuildFromCFile(fname);
1060if(mDlUFunc == NULL)
1061 {cerr<<"PIAFitter::LinkFitFunc() Erreur building .so "<<endl; return;}
1062
1063// On recupere la fonction de fit
1064mFitFunc = (DlUserFitFunc) mDlUFunc->GetFunction(funcname);
1065if(mFitFunc == NULL) // Probleme, on n'a pas trouve etc ...
1066 {cerr<<"PIAFitter::LinkFitFunc() Erreur linking function"<<funcname<<endl;
1067 delete mDlUFunc; mDlUFunc = NULL; return;}
1068
1069// On recupere la fonction des derivees
1070string fnameder = funcname+"_der";
1071mFitFuncDer = (DlUserFitFuncDer) mDlUFunc->GetFunction(fnameder);
1072if(mFitFuncDer == NULL)
1073 {cout<<"PIAFitter::LinkFitFunc() No derive for"<<fnameder<<endl;}
1074
1075// Variables et pointeurs definissant le link a la fonction utilisateur
1076mNameFitFunc = funcname;
1077mUFNVar = nvar;
1078mUFNPar = npar;
1079cout<<"PIAFitter::LinkFitFunc() UserFitFunc="<<mNameFitFunc
1080 <<" NVar="<<mUFNVar<<" NPar="<<mUFNPar<<endl;
1081
1082return;
1083}
1084
1085/* --Methode-- */
1086void PIAFitter::PrepareFitFuncCFile(string const& fname, string const& funcname
1087 , int nvar, int npar)
1088// Pour ecrire un squelette des fichier C avec pour une fonction a fitter
1089{
1090if(fname.length()<=0 || funcname.length()<=0 || nvar<=0 || npar <=0) {
1091 cerr<<"PIAFitter::PrepareFitFuncCFile: error "
1092 <<" nvar="<<nvar<<" npar="<<npar<<endl
1093 <<" file_name="<<fname<<" func_name="<<funcname<<endl
1094 <<" ---- No File Created ! ----"<<endl;
1095 return;
1096}
1097
1098FILE* fip = fopen(fname.c_str(),"w");
1099if(fip==NULL)
1100 {cerr<<"PIAFitter::PrepareFitFuncCFile: error opening file for writing\n"
1101 <<fname<<endl; return;}
1102
1103fprintf(fip,"#include <stdio.h>\n#include <stdlib.h>\n#include <math.h>\n\n");
1104fprintf(fip,"/* - many functions can be written in this file: */\n");
1105fprintf(fip,"/* here is an exemple for creation of function %s */\n\n"
1106 ,funcname.c_str());
1107fprintf(fip,"double %s(double const* x,double const* p)\n",funcname.c_str());
1108fprintf(fip,"/* Function : F(x[0-%d];p[0-%d]) */\n",nvar-1,npar-1);
1109fprintf(fip,"/* variables x1=x[0],x2=x[1],...,x%d=x[%d] */\n",nvar,nvar-1);
1110fprintf(fip,"/* parameters p1=p[0],p2=p[1],...,p%d=p[%d] */\n",nvar,nvar-1);
1111fprintf(fip,"{\n val = F(x[];p[]);\n return(val);\n}\n\n");
1112string fder = funcname + "_der";
1113fprintf(fip,"/* That function is optional */\n");
1114fprintf(fip,"/*\n");
1115fprintf(fip,"double %s(double const* x,double const* p,double* dp)\n",fder.c_str());
1116fprintf(fip,"{\n");
1117for(int i=0;i<npar;i++) fprintf(fip," dp[%d] = dF(x[];p[])/dp[%d]\n",i,i);
1118fprintf(fip," val = F(x[];p[]);\n return(val);\n}\n\n");
1119fprintf(fip,"*/\n");
1120
1121fclose(fip);
1122}
1123
1124/* --Methode-- */
1125int PIAFitter::Execute(string& kw, vector<string>& tokens)
1126// Execution des commandes de fit : "fit" et "fitw"
1127{
1128////// Les commandes pour creer les fonctions de fit utilisateur
1129if(kw == "crfitfil") {
1130 if (tokens.size() < 4)
1131 {cout<<"Usage:crfitfil filename funcname nvar npar"<<endl;return(-1);}
1132 int nvar = atoi(tokens[2].c_str()), npar = atoi(tokens[3].c_str());
1133 PrepareFitFuncCFile(tokens[0],tokens[1],nvar,npar);
1134 return(0);
1135} else if (kw == "crfitfun") {
1136 if (tokens.size() < 4)
1137 {cout<<"Usage:crfitfun filename funcname nvar npar"<<endl;return(-1);}
1138 int nvar = atoi(tokens[2].c_str()), npar = atoi(tokens[3].c_str());
1139 LinkFitFunc(tokens[0],tokens[1],nvar,npar);
1140 return(0);
1141}
1142
1143////// Les commandes pour fitter
1144if(tokens.size() < 2) {
1145 cout
1146 <<"Usage:fit nomobj func [Options]\n"
1147 <<" ou fitw nomobj func [Options]\n"
1148 <<" avec Options: [p:p1,...,pn s:s1,...,sn m:m1,...,mn \n"
1149 <<" M:M1,...,Mn f:f1,...,fn o:.,.,. o:.,.,.]\n"
1150 <<" ou fitlin nomobj func [o:.,.,. o:.,.,.]\n";
1151 return(-1);
1152}
1153
1154// decodage des arguments de la commande
1155string mSp,mSs,mSm,mSM,mSf,mSo;
1156mSp=""; mSs=""; mSm=""; mSM=""; mSf=""; mSo="";
1157if(tokens.size()>2)
1158 for(int ip=2;ip<(int)tokens.size();ip++) {
1159 if(tokens[ip].length()<=2) continue;
1160 const char *c = tokens[ip].c_str();
1161 if(c[1]!=':') continue;
1162 if(c[0]=='p') mSp=c+2;
1163 else if(c[0]=='s') mSs=c+2;
1164 else if(c[0]=='m') mSm=c+2;
1165 else if(c[0]=='M') mSM=c+2;
1166 else if(c[0]=='f') mSf=c+2;
1167 else if(c[0]=='o') {mSo += ","; mSo += c+2;}
1168 }
1169
1170// Execution des commandes
1171DecodeOptions(mSo);
1172DecodeObject(tokens[0]);
1173CheckOptions();
1174if(mOpt.lp>0) PrintOptions();
1175if(kw == "fitlin") ApplyXYcentering=true; else ApplyXYcentering=false;
1176FillGData();
1177DecodeFunction(tokens[1]);
1178DecodeParam(mSp,mSs,mSm,mSM,mSf);
1179if(kw == "fit") {
1180 int rc = DoFit();
1181 if(rc>=0) FitFunRes();
1182} else if(kw == "fitlin") {
1183 LinFit();
1184} else if(kw == "fitw") {
1185 if(FWindFit!=NULL) {delete FWindFit; FWindFit=NULL;}
1186 FWindFit = new PIAFitterWind(mApp,this,tokens[1]);
1187 FWindFit->Show();
1188}
1189
1190return(0);
1191}
1192
1193/* --Methode-- */
1194void PIAFitter::IniFitP1ou2D(void)
1195{
1196if(mNPar==0) return;
1197mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
1198if(!mGData) return;
1199double mean,sigma;
1200int rc = mGData->GetMeanSigma(0,mean,sigma);
1201if(rc<=0) return;
1202mPar(0) = mean;
1203return;
1204}
1205
1206void PIAFitter::IniFitGhP1D(void)
1207{
1208if(mNPar==0) return;
1209mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
1210if(mNVar!=1 || mNPar<3) return;
1211double h,v,x0,y0,sx,sy,f;
1212int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
1213if(rc!=0) return;
1214mPar(0)=h; mStep(0)=h/1000.;
1215mPar(1)=x0; mStep(1)=sx/10.;
1216mPar(2)=sx; mStep(2)=sx/10.; mMin(2)=0.; mMax(2)=20.*sx;
1217if(mNPar>3) {mPar(3)=f; mStep(3)=f/1000.;} else mPar(0)=h+f;
1218}
1219
1220void PIAFitter::IniFitGvP1D(void)
1221{
1222if(mNPar==0) return;
1223mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
1224if(mNVar!=1 || mNPar<3) return;
1225double h,v,x0,y0,sx,sy,f;
1226int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
1227if(rc!=0) return;
1228mPar(0)=v; mStep(0)=v/1000.;
1229mPar(1)=x0; mStep(1)=sx/10.;
1230mPar(2)=sx; mStep(2)=sx/10.; mMin(2)=0.; mMax(2)=20.*sx;
1231if(mNPar>3) {mPar(3)=f; mStep(3)=f/1000.;} else mPar(0)=v+mNData*f;
1232}
1233
1234void PIAFitter::IniFitGv2D(void)
1235{
1236if(mNPar==0) return;
1237mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
1238if(mNVar!=2 || mNPar<6) return;
1239double h,v,x0,y0,sx,sy,f;
1240int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
1241if(rc!=0) return;
1242mPar(0)=v; mStep(0)=v/1000.;
1243mPar(1)=x0; mStep(1)=sx/10.;
1244mPar(2)=y0; mStep(2)=sy/10.;
1245mPar(3)=sx; mStep(3)=sx/10.; mMin(3)=0.; mMax(3)=20.*sx;
1246mPar(4)=sy; mStep(4)=sy/10.; mMin(4)=0.; mMax(4)=20.*sy;
1247mPar(5)=0.; mStep(5)=0.005; mMin(5)=-1.; mMax(5)=1.;
1248if(mNPar>6) {mPar(6)=f; mStep(6)=f/1000.;} else mPar(0)=v+mNData*f;
1249}
1250
1251int PIAFitter::IniFitGaus(GeneralFitData& g,double& h,double& v
1252 ,double& x0,double& y0,double& sx,double& sy,double& f)
1253// Pour initialiser une gaussienne 1D ou 2D
1254// return code =0 si OK
1255{
1256h=v=x0=y0=sx=sy=f=0.;
1257if(!&g) return(-1);
1258int nvar = g.NVar();
1259int ndata = g.NData();
1260if(nvar<1 || 2<nvar || ndata<=0) return(-2);
1261
1262// X et Y ordonnees croissant
1263double *X=NULL, *Y=NULL;
1264X = new double[ndata];
1265if(nvar==2) Y = new double[ndata];
1266{for(int i=0;i<ndata;i++) {X[i] = g.X(i); if(Y) Y[i] = g.Y(i);}}
1267qsort(X,(size_t) ndata,sizeof(double),qSort_Dble);
1268if(Y) qsort(Y,(size_t) ndata,sizeof(double),qSort_Dble);
1269
1270// determination des zones externes=[0,1/3] et [2/3,1]
1271// et des zones internes=[1/3,2/3]
1272int k1=ndata/3,k2=2*ndata/3;
1273double xmin,xmax, ymin=-1.,ymax=1.;
1274xmin = X[k1]; xmax = X[k2];
1275if(Y) {ymin = Y[k1]; ymax = Y[k2];}
1276
1277// recherche des parametres
1278double nb=0, nc=0;
1279{for(int i=0;i<ndata;i++) {
1280 double val,x,y=0.;
1281 x = g.X(i); if(Y) y = g.Y(i); val = g.Val(i);
1282 if(xmin<x && x<xmax && ymin<y && y<ymax) { // c'est le centre
1283 if(nc==0 || val>h) {h=val; x0=x; y0=y;}
1284 v += val; nc++;
1285 } else {f += val; nb++;} // c'est le bord
1286
1287}}
1288if(nb==0 || nc==0) {delete [] X; if(Y) delete [] Y; return(-3);}
1289
1290// Calcul des parametres
1291f /= (double) nb; h -= f; v -= f*nc;
1292sx = (xmax-xmin)/2.; if(Y) sy = (ymax-ymin)/2.;
1293
1294delete [] X; if(Y) delete [] Y;
1295return(0);
1296}
1297
1298
1299////////////////////////////////////////////////////////////////////
1300////---------- Classe de fenetre de fit interactive ------------////
1301////////////////////////////////////////////////////////////////////
1302
1303/* --Methode-- */
1304PIAFitterWind::PIAFitterWind(PIStdImgApp* par, PIAFitter* fiter,string func)
1305: PIWindow((PIMsgHandler*)par, "PIAFitter", PIWK_normal, 240, 240, 150, 150)
1306, mDap(par), mFitter(fiter), mFunc(func)
1307, ReFillGData(false), ReDecodeFunc(false)
1308, ErrValue(1.)
1309, saveI1(0), saveI2(0), saveJ1(0), saveJ2(0)
1310, saveErrValue(1.), saveXC(0.), saveYC(0.)
1311{
1312#ifdef SANS_EVOLPLANCK
1313COMPATIBLE_ASSERT(mFitter);
1314#else
1315ASSERT(mFitter);
1316#endif
1317
1318int npar = mFitter->mNPar;
1319
1320// Alloc de la taille
1321nlab = 13+npar; lab = new PILabel*[nlab]; {for(int i=0;i<nlab;i++) lab[i]=NULL;}
1322ntxt = 10+4*npar; txt = new PIText*[ntxt]; {for(int i=0;i<ntxt;i++) txt[i]=NULL;}
1323nbut = 6; but = new PIButton*[nbut]; {for(int i=0;i<nbut;i++) but[i]=NULL;}
1324nckb = 2+npar; ckb = new PICheckBox*[nckb]; {for(int i=0;i<nckb;i++) ckb[i]=NULL;}
1325npom = 3; pom = new PIOptMenu*[npom]; {for(int i=0;i<npom;i++) pom[i]=NULL;}
1326
1327// On definit la taille a partir de la taille par defaut des composantes
1328// bsx,bsy = environ 6 lettres
1329int bsx,bsy;
1330par->PrefCompSz(bsx, bsy);
1331int spx = (bsx>=10) ? bsx/10 : 1; // intervalle entre lettres X
1332int spy = (bsy>=6) ? bsy/6 : 1; // intervalle entre lettres Y
1333
1334int wszx = (5*bsx+5*spx)+spx;
1335int wszy = 8*(bsy+spy)+npar*(bsy+spy)+2*spy;
1336SetSize(wszx, wszy);
1337int cpx,cpy;
1338
1339// new ligne
1340cpx = spx; cpy = spy;
1341string dum = "Object: "+mFitter->mNObj+" , Fun: ";
1342if(mFitter->mFunc) dum+=mFitter->mNameFitFunc; else dum+=mFitter->mFName;
1343lab[0] = new PILabel(this,dum.c_str(),4*bsx,bsy,cpx,cpy);
1344lab[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1345
1346// new ligne
1347cpx = spx; cpy += bsy+spy;
1348but[0] = new PIButton(this,"Fr.Last",111,bsx,bsy,cpx,cpy);
1349but[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1350cpx += bsx+spx;
1351but[1] = new PIButton(this,"Default",222,bsx,bsy,cpx,cpy);
1352but[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1353cpx += bsx+spx;
1354ckb[0] = new PICheckBox(this,"Gen Fun",1001,bsx,bsy,cpx,cpy);
1355ckb[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1356cpx += bsx+2*spx;
1357ckb[1] = new PICheckBox(this,"Res",1002,bsx,bsy,cpx,cpy);
1358ckb[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1359
1360// new ligne
1361cpx = spx; cpy += bsy+spy;
1362pom[0] = new PIOptMenu(this,"Err def",1.5*bsx,bsy,cpx,cpy);
1363 pom[0]->AppendItem("Err def",2101);
1364 pom[0]->AppendItem("Err cste",2102);
1365 pom[0]->AppendItem("Err sqrt",2103);
1366 pom[0]->SetValue(2101);
1367pom[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1368cpx += int(1.5*bsx+spx);
1369txt[0] = new PIText(this,"",bsx,bsy,cpx,cpy);
1370txt[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1371cpx += bsx+spx;
1372lab[1] = new PILabel(this,"Prt,dbg",bsx,bsy,cpx,cpy);
1373lab[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1374cpx += bsx+spx;
1375txt[1] = new PIText(this,"",bsx,bsy,cpx,cpy);
1376txt[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1377
1378// new ligne
1379cpx = spx; cpy += bsy+spy;
1380lab[2] = new PILabel(this,"Stop Xi2",1.25*bsx,bsy,cpx,cpy);
1381lab[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1382cpx += int(1.25*bsx+spx);
1383txt[2] = new PIText(this,"",bsx,bsy,cpx,cpy);
1384txt[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1385cpx += bsx+spx;
1386lab[3] = new PILabel(this,"Iter",bsx, bsy, cpx, cpy);
1387lab[3]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1388cpx += bsx+spx;
1389txt[3] = new PIText(this,"",bsx,bsy,cpx,cpy);
1390txt[3]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1391
1392// new ligne
1393cpx = spx; cpy += bsy+spy;
1394lab[4] = new PILabel(this,"Range X",bsx,bsy,cpx,cpy);
1395lab[4]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1396cpx += bsx+spx;
1397txt[4] = new PIText(this,"",0.75*bsx,bsy,cpx,cpy);
1398txt[4]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1399cpx += int(0.75*bsx+spx);
1400txt[5] = new PIText(this,"",0.75*bsx,bsy,cpx,cpy);
1401txt[5]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1402cpx += int(0.75*bsx+spx);
1403lab[5] = new PILabel(this,"Y",bsx/2,bsy,cpx,cpy);
1404lab[5]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1405cpx += int(bsx/2.+spx);
1406txt[6] = new PIText(this,"",0.75*bsx,bsy,cpx,cpy);
1407txt[6]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1408cpx += int(0.75*bsx+spx);
1409txt[7] = new PIText(this,"",0.75*bsx,bsy,cpx,cpy);
1410txt[7]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1411
1412// new ligne
1413cpx = spx; cpy += bsy+spy;
1414lab[6] = new PILabel(this,"Cent",bsx/2,bsy,cpx,cpy);
1415lab[6]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1416cpx += int(bsx/2.+spx);
1417pom[1] = new PIOptMenu(this,"X.No",1.25*bsx,bsy,cpx,cpy);
1418 pom[1]->AppendItem("X.No",2001);
1419 pom[1]->AppendItem("X.Value",2002);
1420 pom[1]->AppendItem("X.Auto",2003);
1421 pom[1]->SetValue(2001);
1422pom[1]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1423cpx += int(1.25*bsx+spx);
1424txt[8] = new PIText(this,"",bsx,bsy,cpx,cpy);
1425txt[8]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1426cpx += bsx+spx;
1427pom[2] = new PIOptMenu(this,"Y.No",1.25*bsx,bsy,cpx,cpy);
1428 pom[2]->AppendItem("Y.No",2011);
1429 pom[2]->AppendItem("Y.Value",2012);
1430 pom[2]->AppendItem("Y.Auto",2013);
1431 pom[2]->SetValue(2011);
1432pom[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1433cpx += int(1.25*bsx+spx);
1434txt[9] = new PIText(this,"",bsx,bsy,cpx,cpy);
1435txt[9]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1436
1437// new ligne
1438cpx = spx; cpy += bsy+spy;
1439but[2] = new PIButton(this, "Dismiss",777,bsx,bsy,cpx,cpy);
1440but[2]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1441cpx += bsx+spx;
1442but[3] = new PIButton(this, "Check",666,bsx,bsy,cpx,cpy);
1443but[3]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1444cpx += bsx+spx;
1445but[4] = new PIButton(this, "Fit",333,bsx,bsy,cpx,cpy);
1446but[4]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1447cpx += bsx+spx;
1448but[5] = new PIButton(this, "Ini Auto",444,bsx,bsy,cpx,cpy);
1449but[5]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1450
1451// new ligne
1452cpx = spx; cpy += bsy+spy;
1453lab[7] = new PILabel(this,"Par",bsx/2,bsy,cpx,cpy);
1454lab[7]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1455cpx += int(bsx/2.+spx);
1456lab[8] = new PILabel(this,"Fx",bsx/3,bsy,cpx,cpy);
1457lab[8]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1458cpx += int(bsx/3.+spx);
1459lab[9] = new PILabel(this,"Init",bsx,bsy,cpx,cpy);
1460lab[9]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1461cpx += bsx+spx;
1462lab[10] = new PILabel(this,"Step",bsx,bsy,cpx,cpy);
1463lab[10]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1464cpx += bsx+spx;
1465lab[11] = new PILabel(this,"Min",bsx,bsy,cpx,cpy);
1466lab[11]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1467cpx += bsx+spx;
1468lab[12] = new PILabel(this,"Max",bsx,bsy,cpx,cpy);
1469lab[12]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1470
1471// new lines: Parametres (npar lignes)
1472if(npar>0) {
1473 for(int i=0;i<npar;i++) {
1474 cpx = spx; cpy += bsy+spy;
1475 char str[8]; sprintf(str,"P%d",i);
1476 lab[13+i] = new PILabel(this,str,bsx/2,bsy,cpx,cpy);
1477 lab[13+i]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1478 cpx += int(bsx/2.+spx);
1479 // Attention pas de message entre 1500 et 2000 (permet 500 parametres!)
1480 ckb[2+i] = new PICheckBox(this,"",1500+i,bsx/3,bsy,cpx,cpy+spy);
1481 ckb[2+i]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1482 cpx+=int(bsx/3.+spx);
1483 for(int j=0;j<4;j++) {
1484 txt[10+4*i+j] = new PIText(this,"",bsx,bsy,cpx,cpy);
1485 txt[10+4*i+j]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
1486 cpx += bsx+spx;
1487 }
1488 }
1489}
1490
1491return;
1492}
1493
1494/* --Methode-- */
1495PIAFitterWind::~PIAFitterWind()
1496{
1497int i;
1498for(i=0;i<nlab;i++) if(lab[i]) {delete lab[i]; lab[i]=NULL;}
1499for(i=0;i<ntxt;i++) if(txt[i]) {delete txt[i]; txt[i]=NULL;}
1500for(i=0;i<nbut;i++) if(but[i]) {delete but[i]; but[i]=NULL;}
1501for(i=0;i<nckb;i++) if(ckb[i]) {delete ckb[i]; ckb[i]=NULL;}
1502for(i=0;i<npom;i++) if(pom[i]) {delete pom[i]; pom[i]=NULL;}
1503if(lab) {delete lab; lab=NULL; nlab=0;}
1504if(txt) {delete txt; txt=NULL; ntxt=0;}
1505if(but) {delete but; but=NULL; nbut=0;}
1506if(ckb) {delete ckb; ckb=NULL; nckb=0;}
1507if(pom) {delete pom; pom=NULL; npom=0;}
1508}
1509
1510/* --Methode-- */
1511void PIAFitterWind::Show()
1512{
1513//mDap->SetBlocked();
1514SetState();
1515PIWindow::Show();
1516return;
1517}
1518
1519/* --Methode-- */
1520void PIAFitterWind::SetState(void)
1521// Ecriture et positionnement des ckb et txt
1522{
1523int npar = mFitter->mNPar;
1524
1525// Les check-boxes
1526ckb[0]->SetState(mFitter->mOpt.okfun);
1527ckb[1]->SetState(mFitter->mOpt.okres);
1528if(npar>0)
1529 for(int i=0;i<npar;i++) ckb[2+i]->SetState((mFitter->mFix(i)>=1.)?true:false);
1530
1531// Les champs textes
1532char str[128]; string dum;
1533sprintf(str,"%f",ErrValue); dum=str; txt[0]->SetText(dum);
1534sprintf(str,"%d,%d",mFitter->mOpt.lp,mFitter->mOpt.lpg); dum=str; txt[1]->SetText(dum);
1535sprintf(str,"%f",mFitter->mOpt.stc2); dum=str; txt[2]->SetText(dum);
1536sprintf(str,"%d",mFitter->mOpt.nstep); dum=str; txt[3]->SetText(dum);
1537sprintf(str,"%d",mFitter->mOpt.i1); dum=str; txt[4]->SetText(dum);
1538sprintf(str,"%d",mFitter->mOpt.i2); dum=str; txt[5]->SetText(dum);
1539sprintf(str,"%d",mFitter->mOpt.j1); dum=str; txt[6]->SetText(dum);
1540sprintf(str,"%d",mFitter->mOpt.j2); dum=str; txt[7]->SetText(dum);
1541sprintf(str,"%f",mFitter->mOpt.xc); dum=str; txt[8]->SetText(dum);
1542sprintf(str,"%f",mFitter->mOpt.yc); dum=str; txt[9]->SetText(dum);
1543if(npar>0) {
1544 for(int i=0;i<npar;i++) {
1545 sprintf(str,"%f",mFitter->mPar(i)); dum=str; txt[10+4*i+0]->SetText(dum);
1546 sprintf(str,"%f",mFitter->mStep(i)); dum=str; txt[10+4*i+1]->SetText(dum);
1547 sprintf(str,"%f",mFitter->mMin(i)); dum=str; txt[10+4*i+2]->SetText(dum);
1548 sprintf(str,"%f",mFitter->mMax(i)); dum=str; txt[10+4*i+3]->SetText(dum);
1549 }
1550}
1551
1552return;
1553}
1554
1555/* --Methode-- */
1556void PIAFitterWind::Process(PIMessage msg, PIMsgHandler* sender, void* data)
1557{
1558int lp = (mFitter->mOpt.lp>2)?1:0;
1559int npar = mFitter->mNPar;
1560
1561// Decodage des messages, init et dismiss
1562msg = UserMsg(msg);
1563if(msg ==777) {
1564 if(lp) cout<<"Dismiss"<<endl;
1565 //mDap->SetReady();
1566 Hide();
1567 return;
1568}
1569
1570// L'objet existe t-il encore?
1571NamedObjMgr omg;
1572if(omg.GetObj(mFitter->mNObj) == NULL) {
1573 cout<<"PIAFitterWind::Process Error , Pas d'objet de nom "
1574 <<mFitter->mNObj<<endl;
1575 mFitter->ResetDPointer(); return;
1576}
1577
1578// On recupere les champs textes
1579string dum; char str[128];
1580// Attention gestion du type d'erreur
1581dum=txt[0]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1582 sscanf(str,"%lf",&ErrValue); if(ErrValue<=0.) ErrValue=1.;
1583 if(mFitter->mOpt.err_e>0. ) mFitter->mOpt.err_e = ErrValue;
1584 if(mFitter->mOpt.err_E>0. ) mFitter->mOpt.err_E = ErrValue;
1585 if(saveErrValue!=ErrValue) {
1586 if(lp) cout<<"ErrValue changed: "<<ErrValue<<endl;
1587 ReFillGData = true;
1588 saveErrValue = ErrValue;
1589 }
1590dum=txt[1]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1591 sscanf(str,"%d,%d",&(mFitter->mOpt.lp),&(mFitter->mOpt.lpg));
1592dum=txt[2]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1593 sscanf(str,"%lf",&(mFitter->mOpt.stc2));
1594dum=txt[3]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1595 sscanf(str,"%d",&(mFitter->mOpt.nstep));
1596// Attention si les valeurs de mOpt.i1/.i2/.j1/.j2 ont ete changees
1597dum=txt[4]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1598 sscanf(str,"%d",&(mFitter->mOpt.i1));
1599dum=txt[5]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1600 sscanf(str,"%d",&(mFitter->mOpt.i2));
1601dum=txt[6]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1602 sscanf(str,"%d",&(mFitter->mOpt.j1));
1603dum=txt[7]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1604 sscanf(str,"%d",&(mFitter->mOpt.j2));
1605 if(saveI1!=mFitter->mOpt.i1 || saveI2!=mFitter->mOpt.i2 ||
1606 saveJ1!=mFitter->mOpt.j1 || saveJ2!=mFitter->mOpt.j2 ) {
1607 if(lp)
1608 cout<<"Range changed: I("<<saveI1<<","<<saveI2
1609 <<") -> I("<<mFitter->mOpt.i1<<","<<mFitter->mOpt.i2
1610 <<") ou J("<<saveJ1<<","<<saveJ2
1611 <<") -> J("<<mFitter->mOpt.j1<<","<<mFitter->mOpt.j2
1612 <<")"<<endl;
1613 ReFillGData = true;
1614 saveI1=mFitter->mOpt.i1; saveI2=mFitter->mOpt.i2;
1615 saveJ1=mFitter->mOpt.j1; saveJ2=mFitter->mOpt.j2;
1616 }
1617dum=txt[8]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1618 sscanf(str,"%lf",&(mFitter->mOpt.xc));
1619dum=txt[9]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1620 sscanf(str,"%lf",&(mFitter->mOpt.yc));
1621 if(saveXC!=mFitter->mOpt.xc || saveYC!=mFitter->mOpt.yc) {
1622 if(lp) cout<<"mOpt.xc,yc changed: "<<mFitter->mOpt.xc
1623 <<" , "<<mFitter->mOpt.xc<<endl;
1624 ReFillGData = true; ReDecodeFunc = true;
1625 saveXC = mFitter->mOpt.xc; saveYC = mFitter->mOpt.yc;
1626 }
1627
1628if(npar>0) {
1629 for(int i=0;i<npar;i++) {
1630 dum=txt[10+4*i+0]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1631 sscanf(str,"%lf",&(mFitter->mPar(i)));
1632 dum=txt[10+4*i+1]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1633 sscanf(str,"%lf",&(mFitter->mStep(i)));
1634 dum=txt[10+4*i+2]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1635 sscanf(str,"%lf",&(mFitter->mMin(i)));
1636 dum=txt[10+4*i+3]->GetText(); strcpy(str,dum.c_str()); strip(str,'B',' ');
1637 sscanf(str,"%lf",&(mFitter->mMax(i)));
1638 }
1639}
1640
1641// Decodage des messages, suite
1642if(msg ==666) {
1643 if(lp) cout<<"Check"<<endl;
1644 mFitter->CheckOptions();
1645 if(ReFillGData) mFitter->FillGData();
1646 ReFillGData = false;
1647}
1648else if(msg ==111) {
1649 if(lp) cout<<"Update from last"<<endl;
1650 mFitter->ResetOptions();
1651 mFitter->mOpt = mFitter->mOptSave;
1652 mFitter->ReSetParam();
1653 mFitter->InitParFromLastFit();
1654 ReFillGData = true;
1655}
1656else if(msg ==222) {
1657 if(lp) cout<<"Default options/par"<<endl;
1658 mFitter->ResetOptions();
1659 mFitter->ReSetParam();
1660 ReFillGData = true;
1661}
1662else if(msg ==333) {
1663 if(lp) cout<<"Do the fit"<<endl;
1664 mFitter->CheckOptions();
1665 if(ReFillGData) {
1666 if(lp) cout<<"ReFill Data"<<endl;
1667 mFitter->FillGData();
1668 }
1669 ReFillGData = false;
1670 if(ReDecodeFunc) mFitter->DecodeFunction(mFunc);
1671 ReDecodeFunc = false;
1672 int rc = mFitter->DoFit();
1673 if(rc>=0) mFitter->FitFunRes();
1674}
1675else if(msg ==444) {
1676 if(lp) cout<<"Automatic initialisation"<<endl;
1677 mFitter->AutoIniFit();
1678}
1679else if(msg==1001) {
1680 if(ckb[0]->GetState()) {
1681 if(lp) cout<<"Gen fitted function"<<endl;
1682 mFitter->mOpt.okfun = true;
1683 } else if(!ckb[0]->GetState()) {
1684 if(lp) cout<<"No Gen fitted function"<<endl;
1685 mFitter->mOpt.okfun = false;
1686 }
1687}
1688else if(msg==1002) {
1689 if(ckb[1]->GetState()) {
1690 if(lp) cout<<"Gen fitted residus"<<endl;
1691 mFitter->mOpt.okres = true;
1692 } else if(!ckb[1]->GetState()) {
1693 if(lp) cout<<"No Gen fitted residus"<<endl;
1694 mFitter->mOpt.okres = false;
1695 }
1696}
1697else if(2101<=msg && msg<=2103) {
1698 mFitter->mOpt.err_e = mFitter->mOpt.err_E = -1.;
1699 if(msg==2102) mFitter->mOpt.err_e = ErrValue;
1700 else if(msg==2103) mFitter->mOpt.err_E = ErrValue;
1701 if(lp) cout<<"Err_cste="<<mFitter->mOpt.err_e
1702 <<" Err_sqrt="<<mFitter->mOpt.err_E
1703 <<" ErrValue="<<ErrValue<<endl;
1704 ReFillGData = true;
1705}
1706else if(2001<=msg && msg<=2003) {
1707 mFitter->mOpt.polcx = msg-2001;
1708 if(lp) cout<<"Centrage X polcx="<<mFitter->mOpt.polcx
1709 <<" xc="<<mFitter->mOpt.xc<<endl;
1710 ReFillGData = true; ReDecodeFunc = true;
1711}
1712else if(2011<=msg && msg<=2013) {
1713 mFitter->mOpt.polcy = msg-2011;
1714 if(lp) cout<<"Centrage Y polcy="<<mFitter->mOpt.polcy
1715 <<" yc="<<mFitter->mOpt.yc<<endl;
1716 ReFillGData = true; ReDecodeFunc = true;
1717}
1718else if(msg>=1500 && msg<2000) {
1719 int ip = msg-1500;
1720 if(ckb[2+ip]->GetState()) {
1721 if(lp) cout<<"Fix Param "<<ip<<endl;
1722 mFitter->mFix(ip) = 1.001;
1723 } else if(!ckb[2+ip]->GetState()) {
1724 if(lp) cout<<"Not Fix Param "<<ip<<endl;
1725 mFitter->mFix(ip) = 0.;
1726 }
1727}
1728
1729mFitter->CheckOptions();
1730SetState();
1731if(mFitter->mOpt.lp>1) mFitter->PrintOptions();
1732
1733return;
1734}
Note: See TracBrowser for help on using the repository browser.