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

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

Deplacement des methodes d'ajustement dans nouvelle classe Reza 6/8/99
File size: 22.8 KB
RevLine 
[361]1#include <stdio.h>
2#include <stdlib.h>
3#include <ctype.h>
4
5#include <iostream.h>
6
7#include <typeinfo>
8
9#include "piafitting.h"
10#include "strutil.h"
11
12#include "generalfit.h"
13
14#include "fct1dfit.h"
15#include "fct2dfit.h"
16
17#include "matrix.h"
18#include "cvector.h"
19#include "ntuple.h"
20#include "cimage.h"
21
22#include "histos.h"
23#include "histos2.h"
24#include "ntuple.h"
25#include "hisprof.h"
26
27#include "nobjmgr.h"
28#include "pistdimgapp.h"
29
30//// ---------- Classe de fenetre de fit interactive ------------
31class PIAFitterWind : public PIWindow {
32public :
33 PIAFitterWind(PIStdImgApp* par, PIAFitter* fiter);
34 virtual ~PIAFitterWind();
35 virtual void Show();
36 virtual void Process(PIMessage msg, PIMsgHandler* sender, void* data=NULL);
37 inline void SetObjName(string& nom) { oname = nom; }
38protected :
39 PIStdImgApp* dap;
40 PIAFitter* fitter;
41 PIOptMenu* pom[2];
42 PILabel* lab[2]; // Maximum 20 parametres
43 PIText* txt[2]; // " " "
44 PIButton* but[2];
45 PICheckBox* ckb[2];
46 string oname;
47};
48
49///////////////////// Fit 1D et 2D //////////////////////////
50
51PIAFitterWind* fwind = NULL;
52
53/* --Methode-- */
54PIAFitter::PIAFitter(PIACmd *piac, PIStdImgApp* app)
55{
56string kw, usage;
57
58kw = "fit";
59usage = "Fitting function to DataObjects (Histo, Histo2D, Vector, ...)";
60usage += "\n Usage: fit nomobj func [Options]";
61usage += "\n [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]";
62piac->RegisterCommand(kw, usage, this, "Fitting");
63
64kw = "fitw";
65usage = "Interactive Window Fit";
66usage += "\n Usage: fit nomobj ";
67piac->RegisterCommand(kw, usage, this, "Fitting");
68
69// Je cree une seule copie de PIAFitterWind
70// Attention, ca ne marchera pas si on detruit this
71if (fwind == NULL) fwind = new PIAFitterWind(app, this);
72}
73
74/* --Methode-- */
75PIAFitter::~PIAFitter()
76{
77}
78
79int PIAFitter::Execute(string& kw, vector<string>& tokens)
80{
81// Fit 1D sur objets 1D. Egalement Fit 2D sur objets 2D.
82if (kw == "fit") {
83 if (tokens.size() < 2) {
84 cout <<"Usage:fit nomobj func \n"
85 <<" [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]\n";
86 return(0);
87 }
88 string p=""; string s=""; string m=""; string M=""; string O="";
89 if (tokens.size()>2)
90 for(int ip=2;ip<tokens.size();ip++) {
91 if(tokens[ip].length()<=2) continue;
92 const char *c = tokens[ip].c_str();
93 if(c[1]!=':') continue;
94 if(c[0]=='p') p=c+2;
95 else if(c[0]=='s') s=c+2;
96 else if(c[0]=='m') m=c+2;
97 else if(c[0]=='M') M=c+2;
98 else if(c[0]=='o') {O += ","; O += c+2;}
99 }
100 Fit12D(tokens[0],tokens[1],p,s,m,M,O);
101}
102else if (kw == "fitw") {
103 cout << " *DBG* Affichage PIAFitterWind ** " << endl;
104 fwind->Show();
105 }
106
107return(0);
108}
109
110/* --Methode-- cmv 13/10/98 */
111void PIAFitter::Fit12D(string& nom, string& func,
112 string par,string step,string min,string max,
113 string opt)
114//| --------------- Fit d'objets a 1 et 2 dimensions ---------------
115//| nom : nom de l'objet qui peut etre:
116//| fit-1D: Vector,Histo1D,HProf ou GeneraFitData(1D)
117//| fit-2D: Matrix,Histo2D,Image<T> ou GeneraFitData(2D)
118//| func : pnn : fit polynome degre nn avec classe Poly (lineaire) 1D ou 2D
119//| : Pnn : fit polynome degre nn avec GeneralFit (non-lineaire) 1D ou 2D
120//| : gnn : fit gaussienne (hauteur) + polynome de degre nn 1D
121//| : g : fit gaussienne (hauteur) 1D
122//| : enn : fit exponentielle + polynome de degre nn 1D
123//| : e : fit exponentielle 1D
124//| : Gnn : fit gaussienne (volume) + polynome de degre nn 1D
125//| : G : fit gaussienne (volume) 1D
126//| : : fit gaussienne+fond (volume) 2D
127//| : Gi : fit gaussienne+fond integree (volume) 2D
128//| : d : fit DL de gaussienne+fond (volume) 2D
129//| : di : fit DL de gaussienne+fond integree (volume) 2D
130//| : D : fit DL de gaussienne+fond avec coeff variable p6 (volume) 2D
131//| : Di : fit DL de gaussienne+fond integree avec coeff variable p6 (volume) 2D
132//| : M : fit Moffat+fond (expos=p6) (volume) 2D
133//| : Mi : fit Moffat+fond integree (expos=p6) (volume) 2D
134//| par : p1,...,pn : valeur d'initialisation des parametres (def=0)
135//| step : s1,...,sn : valeur des steps de depart (def=1)
136//| min : m1,...,mn : valeur des minima (def=1)
137//| max : M1,...,Mn : valeur des maxima (def=-1) (max<=min : pas de limite)
138//| opt : options "Eaa.b,eaa.b,f,r,caa.b,Xaa.b"
139//| f : generation d'un Objet identique contenant la fonction fittee
140//| r : generation d'un Objet identique contenant les residus
141//| Xaa.b : aa.b valeur du DXi2 d'arret (def=1.e-3)
142//| Naa : aa nombre maximum d'iterations (def=100)
143//| la.b : niveau "a.b" de print: a=niveau de print Fit1/2D
144//| b=niveau de debug GeneralFit
145//| Ii1/i2 numeros des bins X de l'histos utilises pour le fit [i1,i2]
146//|2D Jj1/j2 numeros des bins Y de l'histos utilises pour le fit [j1,j2]
147//| - L'erreur est celle associee a l'objet (si elle existe),
148//| elle est mise a 1 sinon, sauf si E... ou e... est precise:
149//| Eaa.b : si |val|>=1 erreur = aa.b*sqrt(|val|)
150//| si |val|<1 erreur = aa.b
151//| si aa.b <=0 alors aa.b=1.0
152//| E seul est equivalent a E1.0
153//| eaa.b : erreur = aa.b
154//| si aa.b <=0 alors aa.b=1.0
155//| e seul est equivalent a e1.0
156//| xaa.b : demande de centrage: on fit x-aa.b au lieu de x)
157//| x : demande de centrage: on fit x-xc au lieu de x
158//| avec xc=abscisse du milieu de l'histogramme
159//| Actif pour exp+poly 1D, poly 1D
160//| pour gauss+poly 1D, xc est le centre de la gaussienne.
161//|2D yaa.b et y : idem "xaa.b et x" mais pour y
162{
163
164NamedObjMgr omg;
165AnyDataObj* obj = omg.GetObj(nom);
166if (obj == NULL) {
167 cout<<"PIAFitter::Fit12D() Error , Pas d'objet de nom "<<nom<<endl;
168 return;
169}
170if(func.length()<=0)
171 {cout<<"PIAFitter::Fit12D() Donnez un nom de fonction a fitter."<<endl;
172 return;}
173string ctyp = typeid(*obj).name();
174
175int ndim = 0, nbinx=0, nbiny=0, ndata = 0;
176Vector* v = NULL; Histo* h = NULL;
177Matrix* m = NULL; Histo2D* h2 = NULL; RzImage* im = NULL;
178GeneralFitData* g = NULL;
179
180 // 1D
181if (typeid(*obj) == typeid(Vector)) {
182 ndim = 1;
183 v = (Vector*) obj; nbinx = v->NElts(); nbiny = 1;
184 }
185else if ( (typeid(*obj) == typeid(HProf)) || (typeid(*obj) == typeid(Histo)) ) {
186 ndim = 1;
187 h = (Histo*) obj; nbinx = h->NBins(); nbiny = 1;
188 }
189else if (typeid(*obj) == typeid(Matrix)) {
190 ndim = 2;
191 m = (Matrix*) obj; nbinx = m->NCol(); nbiny = m->NRows();
192 }
193else if (typeid(*obj) == typeid(Histo2D)) {
194 ndim = 2;
195 h2 = (Histo2D*) obj; nbinx = h2->NBinX(); nbiny = h2->NBinY();
196 }
197else if (typeid(*obj) == typeid(GeneralFitData)) {
198 g = (GeneralFitData*) obj; nbinx = g->NData(); nbiny = 1;
199 if( g->NVar()==1) ndim = 1;
200 else if(g->NVar()==2) ndim = 2;
201 else {
202 cout<<"GeneralFitData ne peut avoir que 1 ou 2 variables d'abscisse: "
203 <<((GeneralFitData*) obj)->NVar()<<endl; return; }
204 }
205else if (dynamic_cast<RzImage*>(obj)) {
206 ndim = 2;
207 im = (RzImage*) obj; nbinx = im->XSize(); nbiny = im->YSize();
208 }
209else {
210 cout<<"PIAFitter::Fit12D() Error , Objet n'est pas un "
211 <<"Histo1D/HProf/Vector/Histo2D/Image/Matrix/GeneralFitData "<<ctyp<<endl;
212 return;
213 }
214
215ndata = nbinx*nbiny;
216if(ndata<=0)
217 {cout<<"L'objet a "<<nbinx<<","<<nbiny<<" bins ("<<ndata<<")"<<endl; return;}
218
219// Decodage des options et des parametres, mise en forme
220Vector Par(1); Vector Step(1); Vector Min(1); Vector Max(1); DFOPTIONS O;
221DecodeFitsOptions(par,step,min,max,opt,Par,Step,Min,Max,O);
222O.i1 = (O.i1<0||O.i1>=nbinx)? 0: O.i1;
223O.i2 = (O.i2<0||O.i2>=nbinx||O.i2<O.i1)? nbinx-1: O.i2;
224if(ndim>=2) {
225 O.j1 = (O.j1<0||O.j1>=nbiny)? 0: O.j1;
226 O.j2 = (O.j2<0||O.j2>=nbiny||O.j2<O.j1)? nbiny-1: O.j2;
227} else O.j2 = O.j1 = 0;
228if(O.polcx==2) {
229 if(v||m) O.xc = (O.i2-O.i1+1)/2.;
230 else if(h) O.xc = (h->XMin()+h->XMax())/2.;
231 else if(h2) O.xc = (h2->XMin()+h2->XMax())/2.;
232 else if(g) {double mini,maxi; g->GetMinMax(2,mini,maxi); O.xc=(mini+maxi)/2.;}
233 else if(im) {O.xc = im->XOrg() * im->XPxSize()*(O.i2-O.i1+1)/2.;}
234}
235if(O.polcy==2 && ndim>=2) {
236 if(m) O.yc = (O.j2-O.j1+1)/2.;
237 if(h2) O.yc = (h2->YMin()+h2->YMax())/2.;
238 if(g) {double mini,maxi; g->GetMinMax(12,mini,maxi); O.yc=(mini+maxi)/2.;}
239 if(im) {O.yc = im->YOrg() * im->YPxSize()*(O.j2-O.j1+1)/2.;}
240}
241if(O.lp>0)
242 cout<<"Fit["<<nbinx<<","<<nbiny<<"] ("<<ndata<<") dim="<<ndim<<":"
243 <<" Int=["<<O.i1<<","<<O.i2<<"],["<<O.j1<<","<<O.j2<<"]"<<endl
244 <<" Cent="<<O.polcx<<","<<O.polcy<<","<<O.xc<<"+x"<<","<<O.yc<<"+y"
245 <<" TypE="<<O.err_e<<","<<O.err_E
246 <<" StpX2="<<O.stc2<<" Nstep="<<O.nstep
247 <<" lp,lpg="<<O.lp<<","<<O.lpg<<endl;
248
249///////////////////////////////////
250// Remplissage de GeneralFitData //
251///////////////////////////////////
252GeneralFitData mydata(ndim,ndata,0);
253{for(int i=O.i1;i<=O.i2;i++) for(int j=O.j1;j<=O.j2;j++) {
254 double x,y,f,e;
255
256 if(v)
257 {x= (double) i; f=(*v)(i); e=1.;}
258 else if(h)
259 {x=h->BinCenter(i); f=(*h)(i); e=(h->HasErrors())?h->Error(i):1.;}
260 else if(m)
261 {x=(double) i; y=(double) j; f=(*m)(j,i); e=1.;}
262 else if(h2)
263 {float xf,yf; h2->BinCenter(i,j,xf,yf); x=(double)xf; y=(double)yf;
264 f=(*h2)(i,j); e=(h2->HasErrors())?h2->Error(i,j):1.;}
265 else if(im)
266 {x=im->XOrg()+(i+0.5)*im->XPxSize(); y=im->YOrg()+(j+0.5)*im->YPxSize();
267 f=im->DValue(i,j); e=1.;}
268 else if(g&&ndim==1) {x= g->X(i); f=g->Val(i); e=g->EVal(i);}
269 else if(g&&ndim==2) {x= g->X(i); y= g->Y(i); f=g->Val(i); e=g->EVal(i);}
270 else x=y=f=e=0.;
271
272 // Gestion des erreurs a utiliser
273 if(O.err_e>0.) e=O.err_e;
274 else if(O.err_E>0.) {e=(f<-1.||f>1.)?O.err_E*sqrt(fabs(f)):O.err_E;}
275
276 // Remplissage de generalfit
277 if(func[0]=='p') {x -= O.xc; if(ndim>=2) y -= O.yc;}
278 if(ndim==1) mydata.AddData1(x,f,e);
279 else if(ndim==2) mydata.AddData2(x,y,f,e);
280}}
281if(mydata.NData()<=0)
282 {cout<<"Pas de donnees dans GeneralFitData: "<<mydata.NData()<<endl;
283 return;}
284if(O.lpg>1) {
285 mydata.PrintStatus();
286 mydata.PrintData(0);
287 mydata.PrintData(mydata.NData()-1);
288}
289
290////////////////////////////////////////////
291// Identification de la fonction a fitter //
292////////////////////////////////////////////
293GeneralFunction* myfunc = NULL;
294if(func[0]=='p' && ndim==1) {
295 // Fit de polynome sans passer par les GeneralFit
296 int degre = 0;
297 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
298 cout<<"Fit (lineaire) 1D polynome de degre "<<degre<<endl;
299 Poly p1(0);
300 double c2rl = mydata.PolFit(0,p1,degre);
301 cout<<"C2r_lineaire = "<<c2rl<<endl;
302 if(O.lp>0) cout<<p1<<endl;
303 return;
304
305} else if(func[0]=='P' && ndim==1) {
306 // Fit de polynome
307 int degre = 0;
308 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
309 cout<<"Fit polynome 1D de degre "<<degre<<endl;
310 Polyn1D* myf = new Polyn1D(degre,O.xc);
311 myfunc = myf;
312
313} else if(func[0]=='e' && ndim==1) {
314 // Fit d'exponentielle
315 int degre =-1;
316 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
317 cout<<"Fit d'exponentielle+polynome 1D de degre "<<degre<<endl;
318 Exp1DPol* myf;
319 if(degre>=0) myf = new Exp1DPol((unsigned int)degre,O.xc);
320 else myf = new Exp1DPol(O.xc);
321 myfunc = myf;
322
323} else if(func[0]=='g' && ndim==1) {
324 // Fit de gaussienne en hauteur
325 int degre =-1;
326 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
327 cout<<"Fit de Gaussienne_en_hauteur+polynome 1D de degre "<<degre<<endl;
328 Gauss1DPol* myf;
329 if(degre>=0) myf = new Gauss1DPol((unsigned int)degre,((O.polcx)?true:false));
330 else { bool bfg = (O.polcx)?true:false; myf = new Gauss1DPol(bfg); }
331 myfunc = myf;
332
333} else if(func[0]=='G' && ndim==1) {
334 // Fit de gaussienne en volume
335 int degre =-1;
336 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
337 cout<<"Fit de Gaussienne_en_volume+polynome 1D de degre "<<degre<<endl;
338 GaussN1DPol* myf;
339 if(degre>=0) myf = new GaussN1DPol((unsigned int)degre,((O.polcx)?true:false));
340 else { bool bfg = (O.polcx)?true:false; myf = new GaussN1DPol(bfg); }
341 myfunc = myf;
342
343} else if(func[0]=='p' && ndim==2) {
344 // Fit de polynome 2D sans passer par les GeneralFit
345 int degre = 0;
346 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
347 cout<<"Fit (lineaire) polynome 2D de degre "<<degre<<endl;
348 Poly2 p2(0);
349 double c2rl = mydata.PolFit(0,1,p2,degre);
350 cout<<"C2r_lineaire = "<<c2rl<<endl;
351 if(O.lp>0) cout<<p2<<endl;
352 return;
353
354} else if(func[0]=='P' && ndim==2) {
355 // Fit de polynome 2D
356 int degre = 0;
357 if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
358 cout<<"Fit polynome 2D de degre "<<degre<<endl;
359 Polyn2D* myf = new Polyn2D(degre,O.xc,O.yc);
360 myfunc = myf;
361
362} else if(func[0]=='G' && ndim==2) {
363 // Fit de gaussienne+fond en volume
364 int integ = 0;
365 if(func.length()>1) if(func[1]=='i') integ=1;
366 cout<<"Fit de Gaussienne+Fond 2D integ="<<integ<<endl;
367 if(integ) {GauRhInt2D* myf = new GauRhInt2D; myfunc = myf;}
368 else {GauRho2D* myf = new GauRho2D; myfunc = myf;}
369
370} else if(func[0]=='d' && ndim==2) {
371 // Fit de DL gaussienne+fond en volume
372 int integ = 0;
373 if(func.length()>1) if(func[1]=='i') integ=1;
374 cout<<"Fit de DL de Gaussienne+Fond 2D integ="<<integ<<endl;
375 if(integ) {GdlRhInt2D* myf = new GdlRhInt2D; myfunc = myf;}
376 else {GdlRho2D* myf = new GdlRho2D; myfunc = myf;}
377
378} else if(func[0]=='D' && ndim==2) {
379 // Fit de DL gaussienne+fond avec coeff variable p6 en volume
380 int integ = 0;
381 if(func.length()>1) if(func[1]=='i') integ=1;
382 cout<<"Fit de DL de Gaussienne+Fond avec coeff variable (p6) 2D integ="<<integ<<endl;
383 if(integ) {Gdl1RhInt2D* myf = new Gdl1RhInt2D; myfunc = myf;}
384 else {Gdl1Rho2D* myf = new Gdl1Rho2D; myfunc = myf;}
385
386} else if(func[0]=='M' && ndim==2) {
387 // Fit de Moffat+fond (volume)
388 int integ = 0;
389 if(func.length()>1) if(func[1]=='i') integ=1;
390 cout<<"Fit de Moffat+Fond (expos=p6) 2D integ="<<integ<<endl;
391 if(integ) {MofRhInt2D* myf = new MofRhInt2D; myfunc = myf;}
392 else {MofRho2D* myf = new MofRho2D; myfunc = myf;}
393
394} else {
395 cout<<"Fonction "<<func<<" inconnue pour la dim "<<ndim<<endl;
396 return;
397}
398
399/////////////////////////
400// Fit avec generalfit //
401/////////////////////////
402if(myfunc->NPar()>Par.NElts())
403 {cout<<"Trop de parametres: "<<myfunc->NPar()<<">"<<Par.NElts()<<endl;
404 if(myfunc) delete myfunc; return;}
405GeneralFit myfit(myfunc);
406myfit.SetDebug(O.lpg);
407myfit.SetData(&mydata);
408myfit.SetStopChi2(O.stc2);
409myfit.SetMaxStep(O.nstep);
410{for(int i=0;i<myfunc->NPar();i++) {
411 char str[10];
412 sprintf(str,"P%d",i);
413 myfit.SetParam(i,str,Par(i),Step(i),Min(i),Max(i));
414}}
415if(O.lp>1) myfit.PrintFit();
416double c2r = -1.;
417int rcfit = (double) myfit.Fit();
418if(O.lp>0) myfit.PrintFit();
419if(rcfit>0) {
420 c2r = myfit.GetChi2Red();
421 cout<<"C2r_Reduit = "<<c2r<<" nstep="<<myfit.GetNStep()<<" rc="<<rcfit<<endl;
422 Vector ParFit(myfunc->NPar());
423 for(int i=0;i<myfunc->NPar();i++) ParFit(i)=myfit.GetParm(i);
424} else {
425 cout<<"echec Fit, rc = "<<rcfit<<" nstep="<<myfit.GetNStep()<<endl;
426 myfit.PrintFitErr(rcfit);
427}
428
429// Mise a disposition des resultats
430if(rcfit>=0 && myfunc && (O.okres>0||O.okfun>0)) {
431 string nomres = nom + "res";
432 string nomfun = nom + "fun";
433 if(v) {
434 if(O.okres) {Vector* ob = v->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
435 if(O.okfun) {Vector* ob = v->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
436 } else if(h) {
437 if(O.okres) {Histo* ob = h->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
438 if(O.okfun) {Histo* ob = h->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
439 } else if(m) {
440 if(O.okres) {Matrix* ob = m->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
441 if(O.okfun) {Matrix* ob = m->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
442 } else if(h2) {
443 if(O.okres) {Histo2D* ob = h2->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
444 if(O.okfun) {Histo2D* ob = h2->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
445 } else if(im) {
446 if(O.okres) {RzImage* ob = im->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
447 if(O.okfun) {RzImage* ob = im->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
448 } else if(g) {
449 if(O.okres) {GeneralFitData* ob = g->FitResidus(myfit); if(ob) omg.AddObj(ob,nomres);}
450 if(O.okfun) {GeneralFitData* ob = g->FitFunction(myfit); if(ob) omg.AddObj(ob,nomfun);}
451 }
452}
453
454// Nettoyage
455if(myfunc) delete myfunc;
456return;
457}
458
459
460/* --Function static propre aux routines de fit 1D et 2D-- cmv 13/10/98 */
461void PIAFitter::DecodeFitsOptions(string par,string step,string min,string max,string opt
462 ,Vector& Par,Vector& Step,Vector& Min,Vector& Max,DFOPTIONS& O)
463//| Pour decoder les "string" et remplir les vecteurs du fit (cf commentaires dans Fit1D)
464{
465// set des vecteurs et decodage des string correspondantes
466int NParMax = 100;
467Par.Realloc(NParMax); Step.Realloc(NParMax);
468Min.Realloc(NParMax); Max.Realloc(NParMax);
469{
470 Vector* v=NULL; string* s=NULL;
471 {for(int i=0;i<NParMax;i++) {Par(i)=0.; Step(i)=1.; Min(i)=1.; Max(i)=-1.;}}
472 for(int j=0;j<4;j++) {
473 if(j==0) {v=&Par; s=&par;}
474 else if(j==1) {v=&Step; s=&step;}
475 else if(j==2) {v=&Min; s=&min;}
476 else if(j==3) {v=&Max; s=&max;}
477 if(s->length()>0) *s += ",";
478 for(int i=0;i<NParMax;i++) {
479 if(s->length()<=0) break;
480 sscanf(s->c_str(),"%lf",&(*v)(i));
481 size_t p = s->find_first_of(',') + 1;
482 if(p>=s->length()) *s = ""; else *s = s->substr(p);
483 }
484 }
485}
486
487// Decodage de options de opt
488O.okres = O.okfun = 0;
489O.polcx = O.polcy = 0;
490O.xc = O.yc = 0.;
491O.stc2 = 1.e-3;
492O.nstep = 100;
493O.err_e = O.err_E = -1.;
494O.lp = 1; O.lpg = 0;
495O.i1 = O.j1 = O.i2 = O.j2 = -1;
496
497if(opt.length()<=0) return;
498opt = "," + opt + ",";
499
500if(strstr(opt.c_str(),",r,")) O.okres = 1; // residus
501if(strstr(opt.c_str(),",f,")) O.okfun = 1; // fonction fittee
502if(strstr(opt.c_str(),",x")) { // Demande de centrage (fit X=x-xc)
503 O.polcx = 2; // Le centrage est calcule automatiquement
504 size_t p = opt.find(",x");
505 size_t q = opt.find_first_of(',',p+1);
506 string dum = opt.substr(p,q-p);
507 if(dum.length()>2) {
508 sscanf(dum.c_str(),",x%lf",&O.xc);
509 O.polcx = 1; // Le centrage est fixe par la valeur lue
510 }
511}
512if(strstr(opt.c_str(),",y")) { // Demande de centrage (fit Y=y-yc)
513 O.polcy = 2; // Le centrage est calcule automatiquement
514 size_t p = opt.find(",y");
515 size_t q = opt.find_first_of(',',p+1);
516 string dum = opt.substr(p,q-p);
517 if(dum.length()>2) {
518 sscanf(dum.c_str(),",y%lf",&O.yc);
519 O.polcy = 1; // Le centrage est fixe par la valeur lue
520 }
521}
522if(strstr(opt.c_str(),",E")) { // Erreurs imposees a "sqrt(val)" ou "aa.b*sqrt(val)"
523 size_t p = opt.find(",E");
524 size_t q = opt.find_first_of(',',p+1);
525 string dum = opt.substr(p,q-p);
526 if(dum.length()>2) sscanf(dum.c_str(),",E%lf",&O.err_E);
527 if(O.err_E<=0.) O.err_E = 1.;
528 O.err_e=-1.;
529}
530if(strstr(opt.c_str(),",e")) { // Erreurs imposees a "1" ou "aa.b"
531 size_t p = opt.find(",e");
532 size_t q = opt.find_first_of(',',p+1);
533 string dum = opt.substr(p,q-p);
534 if(dum.length()>2) sscanf(dum.c_str(),",e%lf",&O.err_e);
535 if(O.err_e<=0.) O.err_e = 1.;
536 O.err_E=-1.;
537}
538if(strstr(opt.c_str(),",X")) { // Valeur du StopChi2
539 size_t p = opt.find(",X");
540 size_t q = opt.find_first_of(',',p+1);
541 string dum = opt.substr(p,q-p);
542 if(dum.length()>2) sscanf(dum.c_str(),",X%lf",&O.stc2);
543 if(O.stc2<=0.) O.stc2 = 1.e-3;
544}
545if(strstr(opt.c_str(),",N")) { // Nombre maximum d'iterations
546 size_t p = opt.find(",N");
547 size_t q = opt.find_first_of(',',p+1);
548 string dum = opt.substr(p,q-p);
549 if(dum.length()>2) sscanf(dum.c_str(),",N%d",&O.nstep);
550 if(O.nstep<2) O.nstep = 100;
551}
552if(strstr(opt.c_str(),",l")) { // niveau de print
553 size_t p = opt.find(",l");
554 size_t q = opt.find_first_of(',',p+1);
555 string dum = opt.substr(p,q-p);
556 float ab;
557 if(dum.length()>2) sscanf(dum.c_str(),",l%f",&ab);
558 if(ab<0) ab = 0.;
559 O.lp = (int) ab; O.lpg = int(10.*(ab-(float)O.lp+0.01));
560}
561if(strstr(opt.c_str(),",I")) { // intervalle de fit selon X
562 size_t p = opt.find(",I");
563 size_t q = opt.find_first_of(',',p+1);
564 string dum = opt.substr(p,q-p);
565 if(dum.length()>2) sscanf(dum.c_str(),",I%d/%d",&O.i1,&O.i2);
566}
567if(strstr(opt.c_str(),",J")) { // intervalle de fit selon Y
568 size_t p = opt.find(",J");
569 size_t q = opt.find_first_of(',',p+1);
570 string dum = opt.substr(p,q-p);
571 if(dum.length()>2) sscanf(dum.c_str(),",J%d/%d",&O.j1,&O.j2);
572}
573return;
574}
575
576
577// ---------------- Fenetre de fit interactive -------
578
579
580/* --Methode-- */
581PIAFitterWind::PIAFitterWind(PIStdImgApp* par, PIAFitter* fiter)
582 : PIWindow((PIMsgHandler*)par, "PIAFitter", PIWK_normal, 240, 240, 150, 150)
583{
584dap = par;
585fitter = fiter;
586int bsx, bsy, spx, spy;
587
588// On definit la taille a partir de la taille par defaut des composantes
589// PIApplicationPrefCompSize(bsx, bsy);
590par->PrefCompSz(bsx, bsy);
591spx = bsx/10;
592spy = bsy/4;
593
594int wszx = 5*spx+3*bsx;
595int wszy = 5*bsy+6*spy;
596SetSize(wszx, wszy);
597int cpx = spx;
598int cpy = spy;
599int csx = cpx;
600int csy = cpy;
601lab[0] = new PILabel(this, "Object", 1.5*bsx, bsy, cpx, cpy);
602cpy += spy+bsy;
603lab[1] = new PILabel(this, "Params", 1.5*bsx, bsy, cpx, cpy);
604cpx += bsx*1.5+spx;
605
606cpy = spy;
607txt[0] = new PIText(this, "", 3.5*bsx, bsy, cpx, cpy);
608cpy += spy+bsy;
609txt[1] = new PIText(this, "1.", 3.5*bsx, bsy, cpx, cpy);
610
611cpy += spy+bsy;
612ckb[0] = new PICheckBox(this,"Gen.Func", 1001, 2.5*bsx, bsy, cpx, cpy);
613cpx += bsx*1.5+spx;
614ckb[1] = new PICheckBox(this,"Gen.Resid", 1002, 2.5*bsx, bsy, cpx, cpy);
615
616cpy += spy+bsy;
617pom[0] = new PIOptMenu(this, "FitFunc", 2.5*bsx, bsy, cpx, cpy);
618pom[0]->AppendItem("Poly0", 100);
619pom[0]->AppendItem("Poly1", 101);
620pom[0]->AppendItem("Poly2", 102);
621pom[0]->AppendItem("Gaussienne", 103);
622pom[0]->AppendItem("Gauss+Poly0",104);
623pom[0]->AppendItem("Gauss+Poly1",105);
624pom[0]->AppendItem("Gauss+Poly2",106);
625
626cpx += bsx*1.5+spx;
627pom[1] = new PIOptMenu(this, "Bidon", 2.5*bsx, bsy, cpx, cpy);
628pom[1]->AppendItem("Un", 150);
629pom[1]->AppendItem("Deux",250);
630pom[1]->AppendItem("Trois",350);
631
632cpy += spy+bsy;
633cpx = 0.5*spx;
634but[0] = new PIButton(this, "DoFit", 555, bsx*2, bsy, cpx, cpy);
635cpx += 3*spx;
636but[1] = new PIButton(this, "Dismiss", 777, bsx*2, bsy, cpx, cpy);
637
638}
639/* --Methode-- */
640PIAFitterWind::~PIAFitterWind()
641{
642for(int i=0; i<2; i++) {
643 delete lab[i];
644 delete txt[i];
645 delete pom[i];
646 delete ckb[i];
647 delete but[i];
648 }
649}
650/* --Methode-- */
651void PIAFitterWind::Show()
652{
653// Si on veut initialiser des trucs au moment ou apparait la fenetre
654txt[0]->SetText(oname);
655PIWindow::Show();
656}
657
658/* --Methode-- */
659void PIAFitterWind::Process(PIMessage msg, PIMsgHandler* sender, void* data)
660{
661
662char *mf[7] = {"p0","p1","p2","g","g0","g1","g2"};
663msg = UserMsg(msg);
664if (msg == 777) { Hide(); return; } // On cache la fenetre
665
666if (msg == 555) { // On fait le fit
667 vector<string> args;
668 string cmd = oname;
669 string a;
670 args.push_back(oname);
671 int imf = pom[0]->GetValue()-100;
672 if ((imf < 0) || (imf > 6)) imf = 0;
673 a = mf[imf];
674 cmd += " " + a + " ";
675 args.push_back(a);
676 a = txt[1]->GetText(); // Definition des parametres
677 cmd += a;
678 args.push_back(a);
679 if (ckb[0]->GetState() || ckb[1]->GetState()) {
680 a = "o:";
681 if (ckb[0]->GetState()) a += "f";
682 if (ckb[1]->GetState()) a += "r";
683 cmd += " " + a;
684 args.push_back(a);
685 }
686 cout << "FitComm= fit " << cmd << endl;
687 string kw = "fit";
688 fitter->Execute(kw, args);
689 }
690
691}
Note: See TracBrowser for help on using the repository browser.