Context Navigation

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

Visit:

Last change on this file since 466 was 466, checked in by ercodmgr, 26 years ago

Changement nom /func en /autoc
Nouvelles commandes PAW-like n/plot,...,reset,h/plot/2d etc..
Creation de helpfit pour commenatires fit,fitlin,fitw,crfitfun etc..

cmv 14/10/99

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: