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

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

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