1 | #include "objfitter.h"
|
---|
2 |
|
---|
3 | //===========================================================================
|
---|
4 | //===========================================================================
|
---|
5 | //=========================== ObjectFitter ==================================
|
---|
6 | //===========================================================================
|
---|
7 | //===========================================================================
|
---|
8 |
|
---|
9 | TMatrix<int_4>
|
---|
10 | ObjectFitter::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 |
|
---|
16 | TMatrix<r_4>
|
---|
17 | ObjectFitter::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 |
|
---|
23 | TMatrix<r_8>
|
---|
24 | ObjectFitter::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 |
|
---|
30 | TMatrix<int_4>
|
---|
31 | ObjectFitter::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 |
|
---|
37 | TMatrix<r_4>
|
---|
38 | ObjectFitter::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 |
|
---|
44 | TMatrix<r_8>
|
---|
45 | ObjectFitter::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 |
|
---|
51 | TVector<int_4>
|
---|
52 | ObjectFitter::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 |
|
---|
58 | TVector<r_4>
|
---|
59 | ObjectFitter::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 |
|
---|
65 | TVector<r_8>
|
---|
66 | ObjectFitter::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 |
|
---|
72 | TVector<int_4>
|
---|
73 | ObjectFitter::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 |
|
---|
79 | TVector<r_4>
|
---|
80 | ObjectFitter::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 |
|
---|
86 | TVector<r_8>
|
---|
87 | ObjectFitter::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''. */
|
---|
96 | Histo ObjectFitter::FitResidus(Histo const& hh, GeneralFit& gfit)
|
---|
97 | {
|
---|
98 | Histo h(hh);
|
---|
99 | if(h.NBins()<=0)
|
---|
100 | throw(SzMismatchError("Histo::FitResidus: size mismatch\n"));
|
---|
101 | GeneralFunction* f = gfit.GetFunction();
|
---|
102 | if(f==NULL)
|
---|
103 | throw(NullPtrError("Histo::FitResidus: NULL pointer\n"));
|
---|
104 | TVector<r_8> par = gfit.GetParm();
|
---|
105 | for(int_4 i=0;i<h.NBins();i++) {
|
---|
106 | r_8 x = h.BinCenter(i);
|
---|
107 | h(i) -= f->Value(&x,par.Data());
|
---|
108 | }
|
---|
109 | return h;
|
---|
110 | }
|
---|
111 |
|
---|
112 | /*! Retourne une classe contenant la fonction du fit ``gfit''. */
|
---|
113 | Histo ObjectFitter::FitFunction(Histo const& hh, GeneralFit& gfit)
|
---|
114 | {
|
---|
115 | Histo h(hh);
|
---|
116 | if(h.NBins()<=0)
|
---|
117 | throw(SzMismatchError("Histo::FitFunction: size mismatch\n"));
|
---|
118 | GeneralFunction* f = gfit.GetFunction();
|
---|
119 | if(f==NULL)
|
---|
120 | throw(NullPtrError("Histo::FitFunction: NULL pointer\n"));
|
---|
121 | TVector<r_8> par = gfit.GetParm();
|
---|
122 | for(int_4 i=0;i<h.NBins();i++) {
|
---|
123 | r_8 x = h.BinCenter(i);
|
---|
124 | h(i) = f->Value(&x,par.Data());
|
---|
125 | }
|
---|
126 | return 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 | */
|
---|
149 | int_4 ObjectFitter::Fit(Histo const& h, GeneralFit& gfit,unsigned short typ_err)
|
---|
150 | {
|
---|
151 | if(h.NBins()<=0) return -1000;
|
---|
152 | if(typ_err>5) typ_err=0;
|
---|
153 |
|
---|
154 | GeneralFitData mydata(1,h.NBins());
|
---|
155 |
|
---|
156 | for(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 |
|
---|
170 | gfit.SetData(&mydata);
|
---|
171 |
|
---|
172 | return gfit.Fit();
|
---|
173 | }
|
---|
174 |
|
---|
175 | //============================== Histo2D ====================================
|
---|
176 | /*! Retourne une classe contenant les residus du fit ``gfit''. */
|
---|
177 | Histo2D ObjectFitter::FitResidus(Histo2D const& hh, GeneralFit& gfit)
|
---|
178 | {
|
---|
179 | Histo2D h(hh);
|
---|
180 | if(h.NBinX()<=0 || h.NBinY()<=0)
|
---|
181 | throw(SzMismatchError("Histo2D::FitResidus: size mismatch\n"));
|
---|
182 | GeneralFunction* f = gfit.GetFunction();
|
---|
183 | if(f==NULL)
|
---|
184 | throw(NullPtrError("Histo2D::FitResidus: NULL pointer\n"));
|
---|
185 | TVector<r_8> par = gfit.GetParm();
|
---|
186 | for(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 | }
|
---|
192 | return h;
|
---|
193 | }
|
---|
194 |
|
---|
195 | /*! Retourne une classe contenant la fonction du fit ``gfit''. */
|
---|
196 | Histo2D ObjectFitter::FitFunction(Histo2D const& hh, GeneralFit& gfit)
|
---|
197 | {
|
---|
198 | Histo2D h(hh);
|
---|
199 | if(h.NBinX()<=0 || h.NBinY()<=0)
|
---|
200 | throw(SzMismatchError("Histo2D::FitFunction: size mismatch\n"));
|
---|
201 | GeneralFunction* f = gfit.GetFunction();
|
---|
202 | if(f==NULL)
|
---|
203 | throw(NullPtrError("Histo2D::FitFunction: NULL pointer\n"));
|
---|
204 | TVector<r_8> par = gfit.GetParm();
|
---|
205 | for(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 | }
|
---|
211 | return 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 | */
|
---|
234 | int_4 ObjectFitter::Fit(Histo2D const & h, GeneralFit& gfit,unsigned short typ_err)
|
---|
235 | {
|
---|
236 | if(h.NBinX()*h.NBinY()<=0) return -1000;
|
---|
237 | if(typ_err>5) typ_err=0;
|
---|
238 |
|
---|
239 | GeneralFitData mydata(2,h.NBinX()*h.NBinY());
|
---|
240 |
|
---|
241 | for(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 |
|
---|
256 | gfit.SetData(&mydata);
|
---|
257 |
|
---|
258 | return gfit.Fit();
|
---|
259 | }
|
---|
260 |
|
---|
261 | //===========================================================================
|
---|
262 | //===========================================================================
|
---|
263 | //========================== ArrayFitter<T> =================================
|
---|
264 | //===========================================================================
|
---|
265 | //===========================================================================
|
---|
266 |
|
---|
267 | template <class T>
|
---|
268 | TMatrix<T>
|
---|
269 | ArrayFitter<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 | {
|
---|
276 | if(mtx.NCols()<=0||mtx.NRows()<=0)
|
---|
277 | throw(SzMismatchError("ArrayFitter::FitResidus(TMatrix<T>...) size mismatch\n"));
|
---|
278 | GeneralFunction* f = gfit.GetFunction();
|
---|
279 | if(f==NULL)
|
---|
280 | throw(NullPtrError("ArrayFitter::FitResidus(TMatrix<T>...) GeneraFit==NULL\n"));
|
---|
281 | int npar = gfit.GetNPar();
|
---|
282 | if(npar==0)
|
---|
283 | throw(SzMismatchError("ArrayFitter::FitResidus(TMatrix<T>...) GeneraFit 0 parametre\n"));
|
---|
284 | double* par = new double[npar];
|
---|
285 | {for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
|
---|
286 | TMatrix<T> m(mtx);
|
---|
287 | m.SetTemp(true);
|
---|
288 | for(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 | }
|
---|
292 | delete [] par;
|
---|
293 | return m;
|
---|
294 | }
|
---|
295 |
|
---|
296 |
|
---|
297 | template <class T>
|
---|
298 | TMatrix<T>
|
---|
299 | ArrayFitter<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 | {
|
---|
308 | if(mtx.NCols()<=0||mtx.NRows()<=0)
|
---|
309 | throw(SzMismatchError("ArrayFitter::FitFunction(TMatrix<T>...) size mismatch\n"));
|
---|
310 | GeneralFunction* f = gfit.GetFunction();
|
---|
311 | if(f==NULL)
|
---|
312 | throw(NullPtrError("ArrayFitter::FitFunction(TMatrix<T>...) GeneraFit==NULL\n"));
|
---|
313 | int npar = gfit.GetNPar();
|
---|
314 | if(npar==0)
|
---|
315 | throw(SzMismatchError("ArrayFitter::FitFunction(TMatrix<T>...) GeneraFit 0 parametre\n"));
|
---|
316 | double* par = new double[npar];
|
---|
317 | {for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
|
---|
318 | TMatrix<T> m(mtx.NRows(), mtx.NCols());
|
---|
319 | m.SetTemp(true);
|
---|
320 | for(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 | }
|
---|
324 | delete [] par;
|
---|
325 | return m;
|
---|
326 |
|
---|
327 | }
|
---|
328 |
|
---|
329 | template <class T>
|
---|
330 | TVector<T>
|
---|
331 | ArrayFitter<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 | {
|
---|
336 | if(vec.NElts()<=0)
|
---|
337 | throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) size mismatch\n"));
|
---|
338 | GeneralFunction* f = gfit.GetFunction();
|
---|
339 | if(f==NULL)
|
---|
340 | throw(NullPtrError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit==NULL\n"));
|
---|
341 | int npar = gfit.GetNPar();
|
---|
342 | if(npar==0)
|
---|
343 | throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit 0 parametre\n"));
|
---|
344 | double* par = new double[npar];
|
---|
345 | {for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
|
---|
346 | TVector<T> v(vec);
|
---|
347 | v.SetTemp(true);
|
---|
348 | for(uint_4 i=0;i<vec.NElts();i++) {
|
---|
349 | double x = xorg+i*dx;
|
---|
350 | v(i) -= f->Value(&x,par);
|
---|
351 | }
|
---|
352 | delete [] par;
|
---|
353 | return v;
|
---|
354 | }
|
---|
355 |
|
---|
356 | template <class T>
|
---|
357 | TVector<T>
|
---|
358 | ArrayFitter<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 | {
|
---|
363 | if(vec.NElts()<=0)
|
---|
364 | throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) size mismatch\n"));
|
---|
365 | GeneralFunction* f = gfit.GetFunction();
|
---|
366 | if(f==NULL)
|
---|
367 | throw(NullPtrError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit==NULL\n"));
|
---|
368 | int npar = gfit.GetNPar();
|
---|
369 | if(npar==0)
|
---|
370 | throw(SzMismatchError("ArrayFitter::FitResidus(TVector<T>...) GeneraFit 0 parametre\n"));
|
---|
371 | double* par = new double[npar];
|
---|
372 | {for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
|
---|
373 | TVector<T> v(vec.NElts());
|
---|
374 | v.SetTemp(true);
|
---|
375 | for(uint_4 i=0;i<vec.NElts();i++) {
|
---|
376 | double x = xorg+i*dx;
|
---|
377 | v(i) = f->Value(&x,par);
|
---|
378 | }
|
---|
379 | delete [] par;
|
---|
380 | return 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)
|
---|
394 | template class ArrayFitter<int_4>;
|
---|
395 | template class ArrayFitter<r_4>;
|
---|
396 | template class ArrayFitter<r_8>;
|
---|
397 | template class ArrayFitter< complex<r_4> >;
|
---|
398 | template class ArrayFitter< complex<r_8> >;
|
---|
399 | #endif
|
---|
400 |
|
---|