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

Last change on this file since 4003 was 3572, checked in by cmv, 17 years ago

char* -> const char* pour regler les problemes de deprecated string const... + comparaison unsigned signed + suppression EVOL_PLANCK rz+cmv 07/02/2009

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