source: Sophya/trunk/SophyaLib/NTools/objfitter.cc@ 1204

Last change on this file since 1204 was 1204, checked in by ansari, 25 years ago

pour les adaptateurs de fit cmv 29/9/00

File size: 12.1 KB
Line 
1#include "objfitter.h"
2
3//===========================================================================
4//===========================================================================
5//=========================== ObjectFitter ==================================
6//===========================================================================
7//===========================================================================
8
9TMatrix<int_4>
10ObjectFitter::FitResidus(TMatrix<int_4> const & mtx, GeneralFit& gfit,
11 double xorg, double yorg, double dx, double dy)
12{
13 return( ArrayFitter<int_4>::FitResidus(mtx, gfit, xorg, yorg, dx, dy) );
14}
15
16TMatrix<r_4>
17ObjectFitter::FitResidus(TMatrix<r_4> const & mtx, GeneralFit& gfit,
18 double xorg, double yorg, double dx, double dy)
19{
20 return( ArrayFitter<r_4>::FitResidus(mtx, gfit, xorg, yorg, dx, dy) );
21}
22
23TMatrix<r_8>
24ObjectFitter::FitResidus(TMatrix<r_8> const & mtx, GeneralFit& gfit,
25 double xorg, double yorg, double dx, double dy)
26{
27 return( ArrayFitter<r_8>::FitResidus(mtx, gfit, xorg, yorg, dx, dy) );
28}
29
30TMatrix<int_4>
31ObjectFitter::FitFunction(TMatrix<int_4> const & mtx, GeneralFit& gfit,
32 double xorg, double yorg, double dx, double dy)
33{
34 return( ArrayFitter<int_4>::FitFunction(mtx, gfit, xorg, yorg, dx, dy) );
35}
36
37TMatrix<r_4>
38ObjectFitter::FitFunction(TMatrix<r_4> const & mtx, GeneralFit& gfit,
39 double xorg, double yorg, double dx, double dy)
40{
41 return( ArrayFitter<r_4>::FitFunction(mtx, gfit, xorg, yorg, dx, dy) );
42}
43
44TMatrix<r_8>
45ObjectFitter::FitFunction(TMatrix<r_8> const & mtx, GeneralFit& gfit,
46 double xorg, double yorg, double dx, double dy)
47{
48 return( ArrayFitter<r_8>::FitFunction(mtx, gfit, xorg, yorg, dx, dy) );
49}
50
51TVector<int_4>
52ObjectFitter::FitResidus(TVector<int_4> const & vec, GeneralFit& gfit,
53 double xorg, double dx)
54{
55 return( ArrayFitter<int_4>::FitResidus(vec, gfit, xorg, dx) );
56}
57
58TVector<r_4>
59ObjectFitter::FitResidus(TVector<r_4> const & vec, GeneralFit& gfit,
60 double xorg, double dx)
61{
62 return( ArrayFitter<r_4>::FitResidus(vec, gfit, xorg, dx) );
63}
64
65TVector<r_8>
66ObjectFitter::FitResidus(TVector<r_8> const & vec, GeneralFit& gfit,
67 double xorg, double dx)
68{
69 return( ArrayFitter<r_8>::FitResidus(vec, gfit, xorg, dx) );
70}
71
72TVector<int_4>
73ObjectFitter::FitFunction(TVector<int_4> const & vec, GeneralFit& gfit,
74 double xorg, double dx)
75{
76 return( ArrayFitter<int_4>::FitFunction(vec, gfit, xorg, dx) );
77}
78
79TVector<r_4>
80ObjectFitter::FitFunction(TVector<r_4> const & vec, GeneralFit& gfit,
81 double xorg, double dx)
82{
83 return( ArrayFitter<r_4>::FitFunction(vec, gfit, xorg, dx) );
84}
85
86TVector<r_8>
87ObjectFitter::FitFunction(TVector<r_8> const & vec, GeneralFit& gfit,
88 double xorg, double dx)
89{
90 return( ArrayFitter<r_8>::FitFunction(vec, gfit, xorg, dx) );
91}
92
93
94//=============================== Histo =====================================
95/*! Retourne une classe contenant les residus du fit ``gfit''. */
96Histo ObjectFitter::FitResidus(Histo const& hh, GeneralFit& gfit)
97{
98Histo h(hh);
99if(h.NBins()<=0)
100 throw(SzMismatchError("Histo::FitResidus: size mismatch\n"));
101GeneralFunction* f = gfit.GetFunction();
102if(f==NULL)
103 throw(NullPtrError("Histo::FitResidus: NULL pointer\n"));
104TVector<r_8> par = gfit.GetParm();
105for(int_4 i=0;i<h.NBins();i++) {
106 r_8 x = h.BinCenter(i);
107 h(i) -= f->Value(&x,par.Data());
108}
109return h;
110}
111
112/*! Retourne une classe contenant la fonction du fit ``gfit''. */
113Histo ObjectFitter::FitFunction(Histo const& hh, GeneralFit& gfit)
114{
115Histo h(hh);
116if(h.NBins()<=0)
117 throw(SzMismatchError("Histo::FitFunction: size mismatch\n"));
118GeneralFunction* f = gfit.GetFunction();
119if(f==NULL)
120 throw(NullPtrError("Histo::FitFunction: NULL pointer\n"));
121TVector<r_8> par = gfit.GetParm();
122for(int_4 i=0;i<h.NBins();i++) {
123 r_8 x = h.BinCenter(i);
124 h(i) = f->Value(&x,par.Data());
125}
126return h;
127}
128
129/*!
130 Fit de l'histogramme par ``gfit''.
131 \verbatim
132 typ_err = 0 :
133 - erreur attachee au bin si elle existe
134 - sinon 1
135 typ_err = 1 :
136 - erreur attachee au bin si elle existe
137 - sinon max( sqrt(abs(bin) ,1 )
138 typ_err = 2 :
139 - erreur forcee a 1
140 typ_err = 3 :
141 - erreur forcee a max( sqrt(abs(bin) ,1 )
142 typ_err = 4 :
143 - erreur forcee a 1, nulle si bin a zero.
144 typ_err = 5 :
145 - erreur forcee a max( sqrt(abs(bin) ,1 ),
146 nulle si bin a zero.
147 \endverbatim
148*/
149int_4 ObjectFitter::Fit(Histo const& h, GeneralFit& gfit,unsigned short typ_err)
150{
151if(h.NBins()<=0) return -1000;
152if(typ_err>5) typ_err=0;
153
154GeneralFitData mydata(1,h.NBins());
155
156for(int_4 i=0;i<h.NBins();i++) {
157 r_8 x = h.BinCenter(i);
158 r_8 f = h(i);
159 r_8 saf = sqrt(fabs( f)); if(saf<1.) saf=1.;
160 r_8 e=0.;
161 if(typ_err==0) {if(h.HasErrors()) e=h.Error(i); else e=1.;}
162 else if(typ_err==1) {if(h.HasErrors()) e=h.Error(i); else e=saf;}
163 else if(typ_err==2) e=1.;
164 else if(typ_err==3) e=saf;
165 else if(typ_err==4) e=(f==0.)?0.:1.;
166 else if(typ_err==5) e=(f==0.)?0.:saf;
167 mydata.AddData1(x,f,e);
168}
169
170gfit.SetData(&mydata);
171
172return gfit.Fit();
173}
174
175//============================== Histo2D ====================================
176/*! Retourne une classe contenant les residus du fit ``gfit''. */
177Histo2D ObjectFitter::FitResidus(Histo2D const& hh, GeneralFit& gfit)
178{
179Histo2D h(hh);
180if(h.NBinX()<=0 || h.NBinY()<=0)
181 throw(SzMismatchError("Histo2D::FitResidus: size mismatch\n"));
182GeneralFunction* f = gfit.GetFunction();
183if(f==NULL)
184 throw(NullPtrError("Histo2D::FitResidus: NULL pointer\n"));
185TVector<r_8> par = gfit.GetParm();
186for(int_4 i=0;i<h.NBinX();i++) for(int_4 j=0;j<h.NBinY();j++) {
187 r_8 xc,yc;
188 h.BinCenter(i,j,xc,yc);
189 r_8 x[2] = {xc,yc};
190 h(i,j) -= f->Value(x,par.Data());
191}
192return h;
193}
194
195/*! Retourne une classe contenant la fonction du fit ``gfit''. */
196Histo2D ObjectFitter::FitFunction(Histo2D const& hh, GeneralFit& gfit)
197{
198Histo2D h(hh);
199if(h.NBinX()<=0 || h.NBinY()<=0)
200 throw(SzMismatchError("Histo2D::FitFunction: size mismatch\n"));
201GeneralFunction* f = gfit.GetFunction();
202if(f==NULL)
203 throw(NullPtrError("Histo2D::FitFunction: NULL pointer\n"));
204TVector<r_8> par = gfit.GetParm();
205for(int_4 i=0;i<h.NBinX();i++) for(int_4 j=0;j<h.NBinY();j++) {
206 r_8 xc,yc;
207 h.BinCenter(i,j,xc,yc);
208 r_8 x[2] = {xc,yc};
209 h(i,j) = f->Value(x,par.Data());
210}
211return h;
212}
213
214/*!
215 Fit de l'histogramme par ``gfit''.
216 \verbatim
217 typ_err = 0 :
218 - erreur attachee au bin si elle existe
219 - sinon 1
220 typ_err = 1 :
221 - erreur attachee au bin si elle existe
222 - sinon max( sqrt(abs(bin) ,1 )
223 typ_err = 2 :
224 - erreur forcee a 1
225 typ_err = 3 :
226 - erreur forcee a max( sqrt(abs(bin) ,1 )
227 typ_err = 4 :
228 - erreur forcee a 1, nulle si bin a zero.
229 typ_err = 5 :
230 - erreur forcee a max( sqrt(abs(bin) ,1 ),
231 nulle si bin a zero.
232 \endverbatim
233*/
234int_4 ObjectFitter::Fit(Histo2D const & h, GeneralFit& gfit,unsigned short typ_err)
235{
236if(h.NBinX()*h.NBinY()<=0) return -1000;
237if(typ_err>5) typ_err=0;
238
239GeneralFitData mydata(2,h.NBinX()*h.NBinY());
240
241for(int_4 i=0;i<h.NBinX();i++) for(int_4 j=0;j<h.NBinY();j++) {
242 r_8 x,y;
243 h.BinCenter(i,j,x,y);
244 r_8 f = h(i,j);
245 r_8 saf = sqrt(fabs(f)); if(saf<1.) saf=1.;
246 r_8 e=0.;
247 if(typ_err==0) {if(h.HasErrors()) e=h.Error(i,j); else e=1.;}
248 else if(typ_err==1) {if(h.HasErrors()) e=h.Error(i,j); else e=saf;}
249 else if(typ_err==2) e=1.;
250 else if(typ_err==3) e=saf;
251 else if(typ_err==4) e=(f==0.)?0.:1.;
252 else if(typ_err==5) e=(f==0.)?0.:saf;
253 mydata.AddData2(x,y,f,e);
254}
255
256gfit.SetData(&mydata);
257
258return gfit.Fit();
259}
260
261//===========================================================================
262//===========================================================================
263//========================== ArrayFitter<T> =================================
264//===========================================================================
265//===========================================================================
266
267template <class T>
268TMatrix<T>
269ArrayFitter<T>::FitResidus(TMatrix<T> const & mtx, GeneralFit& gfit,
270 double xorg,double yorg,double dx,double dy)
271// Retourne une classe contenant les residus du fit ``gfit''.
272// On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
273// Les coordonnees de l'element (i,j) sont :
274// (i,j) -> x = xorg+j*dx , y = yorg+i*dy
275{
276if(mtx.NCols()<=0||mtx.NRows()<=0)
277 throw(SzMismatchError("ArrayFitter::FitResidus(TMatrix<T>...) size mismatch\n"));
278GeneralFunction* f = gfit.GetFunction();
279if(f==NULL)
280 throw(NullPtrError("ArrayFitter::FitResidus(TMatrix<T>...) GeneraFit==NULL\n"));
281int npar = gfit.GetNPar();
282if(npar==0)
283 throw(SzMismatchError("ArrayFitter::FitResidus(TMatrix<T>...) GeneraFit 0 parametre\n"));
284double* par = new double[npar];
285{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
286TMatrix<T> m(mtx);
287m.SetTemp(true);
288for(uint_4 i=0;i<mtx.NRows();i++) for(uint_4 j=0;j<mtx.NCols();j++) {
289 double x[2] = {xorg+j*dx,yorg+i*dy};
290 m(i,j) -= f->Value(x,par);
291}
292delete [] par;
293return m;
294}
295
296
297template <class T>
298TMatrix<T>
299ArrayFitter<T>::FitFunction(TMatrix<T> const & mtx, GeneralFit& gfit,
300 double xorg,double yorg,double dx,double dy)
301
302// Retourne une classe contenant la fonction du fit ``gfit''.
303// On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
304// Les coordonnees de l'element (i,j) sont :
305// (i,j) -> x = xorg + j*dx , y = yorg + i*dy
306
307{
308if(mtx.NCols()<=0||mtx.NRows()<=0)
309 throw(SzMismatchError("ArrayFitter::FitFunction(TMatrix<T>...) size mismatch\n"));
310GeneralFunction* f = gfit.GetFunction();
311if(f==NULL)
312 throw(NullPtrError("ArrayFitter::FitFunction(TMatrix<T>...) GeneraFit==NULL\n"));
313int npar = gfit.GetNPar();
314if(npar==0)
315 throw(SzMismatchError("ArrayFitter::FitFunction(TMatrix<T>...) GeneraFit 0 parametre\n"));
316double* par = new double[npar];
317{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
318TMatrix<T> m(mtx.NRows(), mtx.NCols());
319m.SetTemp(true);
320for(uint_4 i=0;i<mtx.NRows();i++) for(uint_4 j=0;j<mtx.NCols();j++) {
321 double x[2] = {xorg+j*dx,yorg+i*dy};
322 m(i,j) = f->Value(x,par);
323}
324delete [] par;
325return m;
326
327}
328
329template <class T>
330TVector<T>
331ArrayFitter<T>::FitResidus(TVector<T> const & vec, GeneralFit& gfit,
332 double xorg,double dx)
333// Retourne une classe contenant les residus du fit ``gfit''.
334// La coordonnee de l'element (i) est -> x = xorg + i*dx
335{
336if(vec.NElts()<=0)
337 throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) size mismatch\n"));
338GeneralFunction* f = gfit.GetFunction();
339if(f==NULL)
340 throw(NullPtrError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit==NULL\n"));
341int npar = gfit.GetNPar();
342if(npar==0)
343 throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit 0 parametre\n"));
344double* par = new double[npar];
345{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
346TVector<T> v(vec);
347v.SetTemp(true);
348for(uint_4 i=0;i<vec.NElts();i++) {
349 double x = xorg+i*dx;
350 v(i) -= f->Value(&x,par);
351}
352delete [] par;
353return v;
354}
355
356template <class T>
357TVector<T>
358ArrayFitter<T>::FitFunction(TVector<T> const & vec, GeneralFit& gfit,
359 double xorg,double dx)
360// Retourne une classe contenant les residus du fit ``gfit''.
361// La coordonnee de l'element (i) est -> x = xorg + i*dx
362{
363if(vec.NElts()<=0)
364 throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) size mismatch\n"));
365GeneralFunction* f = gfit.GetFunction();
366if(f==NULL)
367 throw(NullPtrError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit==NULL\n"));
368int npar = gfit.GetNPar();
369if(npar==0)
370 throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit 0 parametre\n"));
371double* par = new double[npar];
372{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
373TVector<T> v(vec.NElts());
374v.SetTemp(true);
375for(uint_4 i=0;i<vec.NElts();i++) {
376 double x = xorg+i*dx;
377 v(i) = f->Value(&x,par);
378}
379delete [] par;
380return v;
381}
382
383
384///////////////////////////////////////////////////////////////
385#ifdef __CXX_PRAGMA_TEMPLATES__
386#pragma define_template ArrayFitter<int_4>
387#pragma define_template ArrayFitter<r_4>
388#pragma define_template ArrayFitter<r_8>
389#pragma define_template ArrayFitter< complex<r_4> >
390#pragma define_template ArrayFitter< complex<r_8> >
391#endif
392
393#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
394template class ArrayFitter<int_4>;
395template class ArrayFitter<r_4>;
396template class ArrayFitter<r_8>;
397template class ArrayFitter< complex<r_4> >;
398template class ArrayFitter< complex<r_8> >;
399#endif
400
Note: See TracBrowser for help on using the repository browser.