ROOT   Reference Guide
Searching...
No Matches
TFumili.cxx
Go to the documentation of this file.
1// @(#)root/fumili:$Id$
2// Author: Stanislav Nesterov 07/05/2003
3
4
5/** \class TFumili
6
7### FUMILI minimization package
8
9FUMILI is based on ideas, proposed by I.N. Silin [See NIM A440, 2000 (p431)].
10It was converted from FORTRAN to C by Sergey Yaschenko <s.yaschenko@fz-juelich.de>
11
12
13FUMILI is used to minimize Chi-square function or to search maximum of
14likelihood function.
15
16Experimentally measured values \f$F_i\f$ are fitted with theoretical
17functions \f$f_i({\vec x}_i,\vec\theta\,\,)\f$, where \f${\vec x}_i\f$ are
18coordinates, and \f$\vec\theta\f$ -- vector of parameters.
19
20For better convergence Chi-square function has to be the following form
21
22\f[
23{\chi^2\over2}={1\over2}\sum^n_{i=1}\left(f_i(\vec
24x_i,\vec\theta\,\,)-F_i\over\sigma_i\right)^2 \tag{1}
25\f]
26
27where \f$\sigma_i\f$ are errors of measured function.
28
29The minimum condition is
30
31\f[
32{\partial\chi^2\over\partial\theta_i}=\sum^n_{j=1}{1\over\sigma^2_j}\cdot
33{\partial f_j\over\partial\theta_i}\left[f_j(\vec
35\f]
36
37where m is the quantity of parameters.
38
39Expanding left part of (2) over parameter increments and
40retaining only linear terms one gets
41
42\f[
43\left(\partial\chi^2\over\theta_i\right)_{\vec\theta={\vec\theta}^0}
44+\sum_k\left(\partial^2\chi^2\over\partial\theta_i\partial\theta_k\right)_{
45\vec\theta={\vec\theta}^0}\cdot(\theta_k-\theta_k^0)
46= 0\tag{3}
47\f]
48
49Here \f${\vec\theta}_0\f$ is some initial value of parameters. In general case:
50
51\f[
52{\partial^2\chi^2\over\partial\theta_i\partial\theta_k}=
53\sum^n_{j=1}{1\over\sigma^2_j}{\partial f_j\over\theta_i}
54{\partial f_j\over\theta_k} +
55\sum^n_{j=1}{(f_j - F_j)\over\sigma^2_j}\cdot
56{\partial^2f_j\over\partial\theta_i\partial\theta_k}\tag{4}
57\f]
58
59In FUMILI algorithm for second derivatives of Chi-square approximate
60expression is used when last term in (4) is discarded. It is often
61done, not always wittingly, and sometimes causes troubles, for example,
62if user wants to limit parameters with positive values by writing down
63\f$\theta_i^2\f$ instead of \f$\theta_i\f$. FUMILI will fail if one tries
64minimize \f$\chi^2 = g^2(\vec\theta)\f$ where g is arbitrary function.
65
66Approximate value is:
67\f[{\partial^2\chi^2\over\partial\theta_i\partial\theta_k}\approx
68Z_{ik}=
69\sum^n_{j=1}{1\over\sigma^2_j}{\partial f_j\over\theta_i}
70{\partial f_j\over\theta_k}\tag{5}
71\f]
72
73Then the equations for parameter increments are
74\f[\left(\partial\chi^2\over\partial\theta_i\right)_{\vec\theta={\vec\theta}^0}
75+\sum_k Z_{ik}\cdot(\theta_k-\theta^0_k) = 0,
77\f]
78
79Remarkable feature of algorithm is the technique for step
80restriction. For an initial value of parameter \f${\vec\theta}^0\f$ a
81parallelepiped \f$P_0\f$ is built with the center at \f${\vec\theta}^0\f$ and
82axes parallel to coordinate axes \f$\theta_i\f$. The lengths of
83parallelepiped sides along i-th axis is \f$2b_i\f$, where \f$b_i\f$ is such a
84value that the functions \f$f_j(\vec\theta)\f$ are quasi-linear all over
85the parallelepiped.
86
87FUMILI takes into account simple linear inequalities in the form:
88\f[
89\theta_i^{\rm min}\le\theta_i\le\theta^{\rm max}_i\tag{7}
90\f]
91
92They form parallelepiped \f$P\f$ (\f$P_0\f$ may be deformed by \f$P\f$).
93Very similar step formulae are used in FUMILI for negative logarithm
94of the likelihood function with the same idea - linearization of
95function argument.
96
97*/
98
99
100#include "TFumili.h"
101
102#include <iostream>
103#include "TGraphAsymmErrors.h"
104#include "TF1.h"
105#include "TF2.h"
106#include "TF3.h"
107#include "TH1.h"
108#include "TMath.h"
109#include "TROOT.h"
110#include "TList.h"
111#include "TVirtualFitter.h"
112
113
114extern void H1FitChisquareFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
115extern void H1FitLikelihoodFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
116extern void GraphFitChisquareFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
117
118
120
122// Machine dependent values fiXME!!
123// But don't set min=max=0 if param is unlimited
124static const Double_t gMAXDOUBLE=1e300;
125static const Double_t gMINDOUBLE=-1e300;
126
127////////////////////////////////////////////////////////////////////////////////
128
130{//----------- FUMILI constructor ---------
131 // maxpar is the maximum number of parameters used with TFumili object
132 //
133 fMaxParam = TMath::Max(maxpar,25);
134 if (fMaxParam>200) fMaxParam=200;
135 BuildArrays();
136
137 fNumericDerivatives = true;
138 fLogLike = false;
141 fWARN = true;
142 fDEBUG = false;
143 fNlog = 0;
144 fSumLog = 0;
145 fNED1 = 0;
146 fNED2 = 0;
148 fEXDA = 0;
149 fFCN = 0;
150 fNfcn = 0;
151 fRP = 1.e-15; //precision
152 fS = 1e10;
153 fEPS =0.01;
154 fENDFLG = 0;
155 fNlimMul = 2;
156 fNmaxIter= 150;
157 fNstepDec= 3;
158 fLastFixed = -1;
159
160 fAKAPPA = 0.;
161 fGT = 0.;
162 for (int i = 0; i<5; ++i) fINDFLG[i] = 0;
163
164 SetName("Fumili");
165 gFumili = this;
167}
168
169////////////////////////////////////////////////////////////////////////////////
170///
171/// Allocates memory for internal arrays. Called by TFumili::TFumili
172///
173
176 fA = new Double_t[fMaxParam];
177 fPL0 = new Double_t[fMaxParam];
178 fPL = new Double_t[fMaxParam];
180 fDA = new Double_t[fMaxParam];
181 fAMX = new Double_t[fMaxParam];
182 fAMN = new Double_t[fMaxParam];
183 fR = new Double_t[fMaxParam];
184 fDF = new Double_t[fMaxParam];
185 fGr = new Double_t[fMaxParam];
187
188 // fX = new Double_t[10];
189
190 Int_t zSize = fMaxParam*(fMaxParam+1)/2;
191 fZ0 = new Double_t[zSize];
192 fZ = new Double_t[zSize];
193
194 for (Int_t i=0;i<fMaxParam;i++) {
195 fA[i] =0.;
196 fDF[i]=0.;
197 fAMN[i]=gMINDOUBLE;
198 fAMX[i]=gMAXDOUBLE;
199 fPL0[i]=.1;
200 fPL[i] =.1;
201 fParamError[i]=0.;
202 fANames[i]=Form("%d",i);
203 }
204}
205
206////////////////////////////////////////////////////////////////////////////////
207/// TFumili destructor
208
210 DeleteArrays();
211 if (gROOT && !gROOT->TestBit(TObject::kInvalidObject))
212 gROOT->GetListOfSpecials()->Remove(this);
213 if (gFumili == this) gFumili = 0;
214}
215
216////////////////////////////////////////////////////////////////////////////////
217/// return a chisquare equivalent
218
220{
221 Double_t amin = 0;
222 H1FitChisquareFumili(npar,params,amin,params,1);
223 return 2*amin;
224}
225
226////////////////////////////////////////////////////////////////////////////////
227///
228/// Resets all parameter names, values and errors to zero
229///
230/// Argument opt is ignored
231///
232/// NB: this procedure doesn't reset parameter limits
233
235{
237 fNfcn = 0;
238 for (Int_t i=0;i<fNpar;i++) {
239 fA[i] =0.;
240 fDF[i] =0.;
241 fPL0[i] =.1;
242 fPL[i] =.1;
243 fAMN[i] = gMINDOUBLE;
244 fAMX[i] = gMAXDOUBLE;
245 fParamError[i]=0.;
246 fANames[i]=Form("%d",i);
247 }
248}
249
250////////////////////////////////////////////////////////////////////////////////
251/// Deallocates memory. Called from destructor TFumili::~TFumili
252
254 delete[] fCmPar;
255 delete[] fANames;
256 delete[] fDF;
257 // delete[] fX;
258 delete[] fZ0;
259 delete[] fZ;
260 delete[] fGr;
261 delete[] fA;
262 delete[] fPL0;
263 delete[] fPL;
264 delete[] fDA;
265 delete[] fAMN;
266 delete[] fAMX;
267 delete[] fParamError;
268 delete[] fR;
269}
270
271////////////////////////////////////////////////////////////////////////////////
272///
273/// Calculates partial derivatives of theoretical function
274///
275/// Input:
276/// - fX - vector of data point
277///
278/// Output:
279/// - DF - array of derivatives
280///
281/// ARITHM.F: Converted from CERNLIB
282
284 Double_t ff,ai,hi,y,pi;
285 y = EvalTFN(df,fX);
286 for (Int_t i=0;i<fNpar;i++) {
287 df[i]=0;
288 if(fPL0[i]>0.) {
289 ai = fA[i]; // save current parameter value
290 hi = 0.01*fPL0[i]; // diff step
291 pi = fRP*TMath::Abs(ai);
292 if (hi<pi) hi = pi; // if diff step is less than precision
293 fA[i] = ai+hi;
294
295 if (fA[i]>fAMX[i]) { // if param is out of limits
296 fA[i] = ai-hi;
297 hi = -hi;
298 if (fA[i]<fAMN[i]) { // again out of bounds
299 fA[i] = fAMX[i]; // set param to high limit
300 hi = fAMX[i]-ai;
301 if (fAMN[i]-ai+hi<0) { // if hi < (ai-fAMN)
302 fA[i]=fAMN[i];
303 hi=fAMN[i]-ai;
304 }
305 }
306 }
307 ff = EvalTFN(df,fX);
308 df[i] = (ff-y)/hi;
309 fA[i] = ai;
310 }
311 }
312}
313
314
315////////////////////////////////////////////////////////////////////////////////
316/// Evaluate the minimisation function
317///
318/// Input parameters:
319/// - npar: number of currently variable parameters
320/// - par: array of (constant and variable) parameters
321/// - flag: Indicates what is to be calculated
323///
324/// Output parameters:
325/// - fval: The calculated function value.
326/// - grad: The vector of first derivatives.
327///
328/// The meaning of the parameters par is of course defined by the user,
329/// who uses the values of those parameters to calculate their function value.
330/// The starting values must be specified by the user.
331///
332/// Inside FCN user has to define Z-matrix by means TFumili::GetZ
333/// and TFumili::Derivatives,
334/// set theoretical function by means of TFumili::SetUserFunc,
335/// but first - pass number of parameters by TFumili::SetParNumber
336///
337/// Later values are determined by Fumili as it searches for the minimum
338/// or performs whatever analysis is requested by the user.
339///
340/// The default function calls the function specified in SetFCN
341
342Int_t TFumili::Eval(Int_t& npar, Double_t *grad, Double_t &fval, Double_t *par, Int_t flag)
343{
345 return npar;
346}
347
348
349////////////////////////////////////////////////////////////////////////////////
350/// Evaluate theoretical function
351/// - df: array of partial derivatives
352/// - X: vector of theoretical function argument
353
355{
356 // for the time being disable possibility to compute derivatives
357 //if(fTFN)
358 // return (*fTFN)(df,X,fA);
359 //else if(fTFNF1) {
360
361 TF1 *f1 = (TF1*)fUserFunc;
362 return f1->EvalPar(X,fA);
363 //}
364 //return 0.;
365}
366
367////////////////////////////////////////////////////////////////////////////////
368///
369/// Execute MINUIT commands. MINImize, SIMplex, MIGrad and FUMili all
370/// will call TFumili::Minimize method.
371///
372/// For full command list see
373/// MINUIT. Reference Manual. CERN Program Library Long Writeup D506.
374///
375/// Improvement and errors calculation are not yet implemented as well
376/// as Monte-Carlo seeking and minimization.
377/// Contour commands are also unsupported.
378///
379/// - command : command string
380/// - args : array of arguments
381/// - nargs : number of arguments
382
383Int_t TFumili::ExecuteCommand(const char *command, Double_t *args, Int_t nargs){
384 TString comand = command;
385 static TString clower = "abcdefghijklmnopqrstuvwxyz";
386 static TString cupper = "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
387 const Int_t nntot = 40;
388 const char *cname[nntot] = {
389 "MINImize ", // 0 checked
390 "SEEk ", // 1 none
391 "SIMplex ", // 2 checked same as 0
392 "MIGrad ", // 3 checked same as 0
393 "MINOs ", // 4 none
394 "SET xxx ", // 5 lot of stuff
395 "SHOw xxx ", // 6 -----------
396 "TOP of pag", // 7 .
397 "fiX ", // 8 .
398 "REStore ", // 9 .
399 "RELease ", // 10 .
400 "SCAn ", // 11 not yet implemented
401 "CONtour ", // 12 not yet implemented
402 "HESse ", // 13 not yet implemented
403 "SAVe ", // 14 obsolete
404 "IMProve ", // 15 not yet implemented
405 "CALl fcn ", // 16 .
406 "STAndard ", // 17 .
407 "END ", // 18 .
408 "EXIt ", // 19 .
409 "RETurn ", // 20 .
410 "CLEar ", // 21 .
411 "HELP ", // 22 not yet implemented
412 "MNContour ", // 23 not yet implemented
413 "STOp ", // 24 .
414 "JUMp ", // 25 not yet implemented
415 " ", //
416 " ", //
417 "FUMili ", // 28 checked same as 0
418 " ", //
419 " ", //
420 " ", //
421 " ", //
422 "COVARIANCE", // 33
423 "PRINTOUT ", // 34
425 "MATOUT ", // 36
426 "ERROR DEF ", // 37
427 "LIMITS ", // 38
428 "PUNCH "}; // 39
429
430
431 fCword = comand;
432 fCword.ToUpper();
433 if (nargs<=0) fCmPar[0] = 0;
434 Int_t i;
435 for(i=0;i<fMaxParam;i++) {
436 if(i<nargs) fCmPar[i] = args[i];
437 }
438 /*
439 fNmaxIter = int(fCmPar[0]);
440 if (fNmaxIter <= 0) {
441 fNmaxIter = fNpar*10 + 20 + fNpar*M*5;
442 }
443 fEPS = fCmPar[1];
444 */
445 //*-*- look for command in list CNAME . . . . . . . . . .
446 TString ctemp = fCword(0,3);
447 Int_t ind;
448 for (ind = 0; ind < nntot; ++ind) {
449 if (strncmp(ctemp.Data(),cname[ind],3) == 0) break;
450 }
451 if (ind==nntot) return -3; // Unknown command - input ignored
452 if (fCword(0,4) == "MINO") ind=3;
453 switch (ind) {
454 case 0: case 3: case 2: case 28:
455 // MINImize [maxcalls] [tolerance]
456 // also SIMplex, MIGrad and FUMili
457 if(nargs>=1)
458 fNmaxIter=TMath::Max(Int_t(fCmPar[0]),fNmaxIter); // fiXME!!
459 if(nargs==2)
460 fEPS=fCmPar[1];
461 return Minimize();
462 case 1:
463 // SEEk not implemented in this package
464 return -10;
465
466 case 4: // MINos errors analysis not implemented
467 return -10;
468
469 case 5: case 6: // SET xxx & SHOW xxx
470 return ExecuteSetCommand(nargs);
471
472 case 7: // Obsolete command
473 Printf("1");
474 return 0;
475 case 8: // fiX <parno> ....
476 if (nargs<1) return -1; // No parameters specified
477 for (i=0;i<nargs;i++) {
478 Int_t parnum = Int_t(fCmPar[i])-1;
479 FixParameter(parnum);
480 }
481 return 0;
482 case 9: // REStore <code>
483 if (nargs<1) return 0;
484 if(fCmPar[0]==0.)
485 for (i=0;i<fNpar;i++)
487 else
488 if(fCmPar[0]==1.) {
490 std::cout <<fLastFixed<<std::endl;
491 }
492 return 0;
493 case 10: // RELease <parno> ...
494 if (nargs<1) return -1; // No parameters specified
495 for (i=0;i<nargs;i++) {
496 Int_t parnum = Int_t(fCmPar[i])-1;
497 ReleaseParameter(parnum);
498 }
499 return 0;
500 case 11: // SCAn not implemented
501 return -10;
502 case 12: // CONt not implemented
503 return -10;
504
505 case 13: // HESSe not implemented
506 return -10;
507 case 14: // SAVe
508 Printf("SAVe command is obsolete");
509 return -10;
510 case 15: // IMProve not implemented
511 return -10;
512 case 16: // CALl fcn <iflag>
513 {if(nargs<1) return -1;
514 Int_t flag = Int_t(fCmPar[0]);
515 Double_t fval;
516 Eval(fNpar,fGr,fval,fA,flag);
517 return 0;}
518 case 17: // STAndard must call function STAND
519 return 0;
520 case 18: case 19:
521 case 20: case 24: {
522 Double_t fval;
523 Int_t flag = 3;
524 Eval(fNpar,fGr,fval,fA,flag);
525 return 0;
526 }
527 case 21:
528 Clear();
529 return 0;
530 case 22: //HELp not implemented
531 case 23: //MNContour not implemented
532 case 25: // JUMp not implemented
533 return -10;
534 case 26: case 27: case 29: case 30: case 31: case 32:
535 return 0; // blank commands
536 case 33: case 34: case 35: case 36: case 37: case 38:
537 case 39:
538 Printf("Obsolete command. Use corresponding SET command instead");
539 return -10;
540 default:
541 break;
542 }
543 return 0;
544}
545
546////////////////////////////////////////////////////////////////////////////////
547/// Called from TFumili::ExecuteCommand in case
548/// of "SET xxx" and "SHOW xxx".
549
551 static Int_t nntot = 30;
552 static const char *cname[30] = {
553 "FCN value ", // 0 .
554 "PARameters", // 1 .
555 "LIMits ", // 2 .
556 "COVariance", // 3 .
557 "CORrelatio", // 4 .
558 "PRInt levl", // 5 not implemented yet
560 "GRAdient ", // 7 .
561 "ERRor def ", // 8 not sure how to implement - by time being ignored
562 "INPut file", // 9 not implemented
563 "WIDth page", // 10 not implemented yet
564 "LINes page", // 11 not implemented yet
565 "NOWarnings", // 12 .
566 "WARnings ", // 13 .
567 "RANdom gen", // 14 not implemented
568 "TITle ", // 15 ignored
569 "STRategy ", // 16 ignored
570 "EIGenvalue", // 17 not implemented yet
571 "PAGe throw", // 18 ignored
572 "MINos errs", // 19 not implemented yet
573 "EPSmachine", // 20 .
574 "OUTputfile", // 21 not implemented
575 "BATch ", // 22 ignored
576 "INTeractiv", // 23 ignored
577 "VERsion ", // 24 .
578 "reserve ", // 25 .
579 "NODebug ", // 26 .
580 "DEBug ", // 27 .
581 "SHOw ", // 28 err
582 "SET "};// 29 err
583
584 TString cfname, cmode, ckind, cwarn, copt, ctemp, ctemp2;
585 Int_t i, ind;
586 Bool_t setCommand=kFALSE;
587 for (ind = 0; ind < nntot; ++ind) {
588 ctemp = cname[ind];
589 ckind = ctemp(0,3);
590 ctemp2 = fCword(4,6);
591 if (strstr(ctemp2.Data(),ckind.Data())) break;
592 }
593 ctemp2 = fCword(0,3);
594 if(ctemp2.Contains("SET")) setCommand=true;
595 if(ctemp2.Contains("HEL") || ctemp2.Contains("SHO")) setCommand=false;
596
597 if (ind>=nntot) return -3;
598
599 switch(ind) {
600 case 0: // SET FCN value illegal // SHOw only
601 if(!setCommand) Printf("FCN=%f",fS);
602 return 0;
603 case 1: // PARameter <parno> <value>
604 {
605 if (nargs<2 && setCommand) return -1;
606 Int_t parnum;
607 Double_t val;
608 if(setCommand) {
609 parnum = Int_t(fCmPar[0])-1;
610 val= fCmPar[1];
611 if(parnum<0 || parnum>=fNpar) return -2; //no such parameter
612 fA[parnum] = val;
613 } else {
614 if (nargs>0) {
615 parnum = Int_t(fCmPar[0])-1;
616 if(parnum<0 || parnum>=fNpar) return -2; //no such parameter
617 Printf("Parameter %s = %E",fANames[parnum].Data(),fA[parnum]);
618 } else
619 for (i=0;i<fNpar;i++)
620 Printf("Parameter %s = %E",fANames[i].Data(),fA[i]);
621
622 }
623 return 0;
624 }
625 case 2: // LIMits [parno] [ <lolim> <uplim> ]
626 {
627 Int_t parnum;
628 Double_t lolim,uplim;
629 if (nargs<1) {
630 for(i=0;i<fNpar;i++)
631 if(setCommand) {
632 fAMN[i] = gMINDOUBLE;
633 fAMX[i] = gMAXDOUBLE;
634 } else
635 Printf("Limits for param %s: Low=%E, High=%E",
636 fANames[i].Data(),fAMN[i],fAMX[i]);
637 } else {
638 parnum = Int_t(fCmPar[0])-1;
639 if(parnum<0 || parnum>=fNpar)return -1;
640 if(setCommand) {
641 if(nargs>2) {
642 lolim = fCmPar[1];
643 uplim = fCmPar[2];
644 if(uplim==lolim) return -1;
645 if(lolim>uplim) {
646 Double_t tmp = lolim;
647 lolim = uplim;
648 uplim = tmp;
649 }
650 } else {
651 lolim = gMINDOUBLE;
652 uplim = gMAXDOUBLE;
653 }
654 fAMN[parnum] = lolim;
655 fAMX[parnum] = uplim;
656 } else
657 Printf("Limits for param %s Low=%E, High=%E",
658 fANames[parnum].Data(),fAMN[parnum],fAMX[parnum]);
659 }
660 return 0;
661 }
662 case 3:
663 {
664 if(setCommand) return 0;
665 Printf("\nCovariant matrix ");
666 Int_t l = 0,nn=0,nnn=0;
667 for (i=0;i<fNpar;i++) if(fPL0[i]>0.) nn++;
668 for (i=0;i<nn;i++) {
669 for(;fPL0[nnn]<=0.;nnn++) { }
670 printf("%5s: ",fANames[nnn++].Data());
671 for (Int_t j=0;j<=i;j++)
672 printf("%11.2E",fZ[l++]);
673 std::cout<<std::endl;
674 }
675 std::cout<<std::endl;
676 return 0;
677 }
678 case 4:
679 if(setCommand) return 0;
680 Printf("\nGlobal correlation factors (maximum correlation of the parameter\n with arbitrary linear combination of other parameters)");
681 for(i=0;i<fNpar;i++) {
682 printf("%5s: ",fANames[i].Data());
683 printf("%11.3E\n",TMath::Sqrt(1-1/((fR[i]!=0.)?fR[i]:1.)) );
684 }
685 std::cout<<std::endl;
686 return 0;
687 case 5: // PRIntout not implemented
688 return -10;
690 if(!setCommand) return 0;
692 return 0;
694 if(!setCommand) return 0;
696 return 0;
697 case 8: // ERRordef - now ignored
698 return 0;
699 case 9: // INPut - not implemented
700 return -10;
701 case 10: // WIDthpage - not implemented
702 return -10;
703 case 11: // LINesperpage - not implemented
704 return -10;
705 case 12: //NOWarnings
706 if(!setCommand) return 0;
707 fWARN = false;
708 return 0;
709 case 13: // WARnings
710 if(!setCommand) return 0;
711 fWARN = true;
712 return 0;
713 case 14: // RANdomgenerator - not implemented
714 return -10;
715 case 15: // TITle - ignored
716 return 0;
717 case 16: // STRategy - ignored
718 return 0;
719 case 17: // EIGenvalues - not implemented
720 return -10;
721 case 18: // PAGethrow - ignored
722 return 0;
723 case 19: // MINos errors - not implemented
724 return -10;
725 case 20: //EPSmachine
726 if(!setCommand) {
727 Printf("Relative floating point precision RP=%E",fRP);
728 } else
729 if (nargs>0) {
730 Double_t pres=fCmPar[0];
731 if (pres<1e-5 && pres>1e-34) fRP=pres;
732 }
733 return 0;
734 case 21: // OUTputfile - not implemented
735 return -10;
736 case 22: // BATch - ignored
737 return 0;
738 case 23: // INTerative - ignored
739 return 0;
740 case 24: // VERsion
741 if(setCommand) return 0;
742 Printf("FUMILI-ROOT version 0.1");
743 return 0;
744 case 25: // reserved
745 return 0;
746 case 26: // NODebug
747 if(!setCommand) return 0;
748 fDEBUG = false;
749 return 0;
750 case 27: // DEBug
751 if(!setCommand) return 0;
752 fDEBUG = true;
753 return 0;
754 case 28:
755 case 29:
756 return -3;
757 default:
758 break;
759 }
760 return -3;
761}
762
763////////////////////////////////////////////////////////////////////////////////
764/// Fixes parameter number ipar
765
767 if(ipar>=0 && ipar<fNpar && fPL0[ipar]>0.) {
768 fPL0[ipar] = -fPL0[ipar];
769 fLastFixed = ipar;
770 }
771}
772
773////////////////////////////////////////////////////////////////////////////////
774/// Return a pointer to the covariance matrix
775
777{
778 return fZ;
779
780}
781
782////////////////////////////////////////////////////////////////////////////////
783/// Return element i,j from the covariance matrix
784
786{
787 if (!fZ) return 0;
788 if (i < 0 || i >= fNpar || j < 0 || j >= fNpar) {
789 Error("GetCovarianceMatrixElement","Illegal arguments i=%d, j=%d",i,j);
790 return 0;
791 }
792 return fZ[j+fNpar*i];
793}
794
795////////////////////////////////////////////////////////////////////////////////
796/// Return the total number of parameters (free + fixed)
797
799{
800 return fNpar;
801}
802
803////////////////////////////////////////////////////////////////////////////////
804/// Return the number of free parameters
805
807{
808 Int_t nfree = fNpar;
809 for (Int_t i=0;i<fNpar;i++) {
810 if (IsFixed(i)) nfree--;
811 }
812 return nfree;
813}
814
815////////////////////////////////////////////////////////////////////////////////
816/// Return error of parameter ipar
817
819{
820 if (ipar<0 || ipar>=fNpar) return 0;
821 else return fParamError[ipar];
822}
823
824////////////////////////////////////////////////////////////////////////////////
825/// Return current value of parameter ipar
826
828{
829 if (ipar<0 || ipar>=fNpar) return 0;
830 else return fA[ipar];
831}
832
833////////////////////////////////////////////////////////////////////////////////
834/// Get various ipar parameter attributes:
835///
836/// - cname: parameter name
837/// - value: parameter value
838/// - verr: parameter error
839/// - vlow: lower limit
840/// - vhigh: upper limit
841///
842/// WARNING! parname must be suitably dimensioned in the calling function.
843
844Int_t TFumili::GetParameter(Int_t ipar,char *cname,Double_t &value,Double_t &verr,Double_t &vlow, Double_t &vhigh) const
845{
846 if (ipar<0 || ipar>=fNpar) {
847 value = 0;
848 verr = 0;
849 vlow = 0;
850 vhigh = 0;
851 return -1;
852 }
853 strcpy(cname,fANames[ipar].Data());
854 value = fA[ipar];
855 verr = fParamError[ipar];
856 vlow = fAMN[ipar];
857 vhigh = fAMX[ipar];
858 return 0;
859}
860
861////////////////////////////////////////////////////////////////////////////////
862/// Return name of parameter ipar
863
864const char *TFumili::GetParName(Int_t ipar) const
865{
866 if (ipar < 0 || ipar > fNpar) return "";
867 return fANames[ipar].Data();
868}
869
870////////////////////////////////////////////////////////////////////////////////
871/// Return errors after MINOs
872/// not implemented
873
874Int_t TFumili::GetErrors(Int_t ipar,Double_t &eplus, Double_t &eminus, Double_t &eparab, Double_t &globcc) const
875{
876 eparab = 0;
877 globcc = 0;
878 if (ipar<0 || ipar>=fNpar) {
879 eplus = 0;
880 eminus = 0;
881 return -1;
882 }
883 eplus=fParamError[ipar];
884 eminus=-eplus;
885 return 0;
886}
887
888////////////////////////////////////////////////////////////////////////////////
889/// Return global fit parameters
890/// - amin : chisquare
891/// - edm : estimated distance to minimum
892/// - errdef
893/// - nvpar : number of variable parameters
894/// - nparx : total number of parameters
895
896Int_t TFumili::GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const
897{
898 amin = 2*fS;
899 edm = fGT; //
900 errdef = 0; // ??
901 nparx = fNpar;
902 nvpar = 0;
903 for(Int_t ii=0; ii<fNpar; ii++) {
904 if(fPL0[ii]>0.) nvpar++;
905 }
906 return 0;
907}
908
909////////////////////////////////////////////////////////////////////////////////
910/// Return Sum(log(i) i=0,n
911/// used by log-likelihood fits
912
914{
915 if (n < 0) return 0;
916 if (n > fNlog) {
917 if (fSumLog) delete [] fSumLog;
918 fNlog = 2*n+1000;
919 fSumLog = new Double_t[fNlog+1];
920 Double_t fobs = 0;
921 for (Int_t j=0;j<=fNlog;j++) {
922 if (j > 1) fobs += TMath::Log(j);
923 fSumLog[j] = fobs;
924 }
925 }
926 if (fSumLog) return fSumLog[n];
927 return 0;
928}
929
930////////////////////////////////////////////////////////////////////////////////
931/// Inverts packed diagonal matrix Z by square-root method.
932/// Matrix elements corresponding to
933/// fix parameters are removed.
934///
935/// - n: number of variable parameters
936
938{
939 static Double_t am = 3.4e138;
940 static Double_t rp = 5.0e-14;
941 Double_t ap, aps, c, d;
942 Double_t *r_1=fR;
943 Double_t *pl_1=fPL;
944 Double_t *z_1=fZ;
945 Int_t i, k, l, ii, ki, li, kk, ni, ll, nk, nl, ir, lk;
946 if (n < 1) {
947 return;
948 }
949 --pl_1;
950 --r_1;
951 --z_1;
952 aps = am / n;
953 aps = sqrt(aps);
954 ap = 1.0e0 / (aps * aps);
955 ir = 0;
956 for (i = 1; i <= n; ++i) {
957 L1:
958 ++ir;
959 if (pl_1[ir] <= 0.0e0) goto L1;
960 else goto L2;
961 L2:
962 ni = i * (i - 1) / 2;
963 ii = ni + i;
964 k = n + 1;
965 if (z_1[ii] <= rp * TMath::Abs(r_1[ir]) || z_1[ii] <= ap) {
966 goto L19;
967 }
968 z_1[ii] = 1.0e0 / sqrt(z_1[ii]);
969 nl = ii - 1;
970 L3:
971 if (nl - ni <= 0) goto L5;
972 else goto L4;
973 L4:
974 z_1[nl] *= z_1[ii];
975 if (TMath::Abs(z_1[nl]) >= aps) {
976 goto L16;
977 }
978 --nl;
979 goto L3;
980 L5:
981 if (i - n >= 0) goto L12;
982 else goto L6;
983 L6:
984 --k;
985 nk = k * (k - 1) / 2;
986 nl = nk;
987 kk = nk + i;
988 d = z_1[kk] * z_1[ii];
989 c = d * z_1[ii];
990 l = k;
991 L7:
992 ll = nk + l;
993 li = nl + i;
994 z_1[ll] -= z_1[li] * c;
995 --l;
996 nl -= l;
997 if (l - i <= 0) goto L9;
998 else goto L7;
999 L8:
1000 ll = nk + l;
1001 li = ni + l;
1002 z_1[ll] -= z_1[li] * d;
1003 L9:
1004 --l;
1005 if (l <= 0) goto L10;
1006 else goto L8;
1007 L10:
1008 z_1[kk] = -c;
1009 if (k - i - 1 <= 0) goto L11;
1010 else goto L6;
1011 L11:
1012 ;
1013 }
1014 L12:
1015 for (i = 1; i <= n; ++i) {
1016 for (k = i; k <= n; ++k) {
1017 nl = k * (k - 1) / 2;
1018 ki = nl + i;
1019 d = 0.0e0;
1020 for (l = k; l <= n; ++l) {
1021 li = nl + i;
1022 lk = nl + k;
1023 d += z_1[li] * z_1[lk];
1024 nl += l;
1025 }
1026 ki = k * (k - 1) / 2 + i;
1027 z_1[ki] = d;
1028 }
1029 }
1030 L15:
1031 return;
1032 L16:
1033 k = i + nl - ii;
1034 ir = 0;
1035 for (i = 1; i <= k; ++i) {
1036 L17:
1037 ++ir;
1038 if (pl_1[ir] <= 0.0e0) {
1039 goto L17;
1040 }
1041 }
1042 L19:
1043 pl_1[ir] = -2.0e0;
1044 r_1[ir] = 0.0e0;
1045 fINDFLG[0] = ir - 1;
1046 goto L15;
1047}
1048
1049////////////////////////////////////////////////////////////////////////////////
1050/// Return kTRUE if parameter ipar is fixed, kFALSE otherwise)
1051
1053{
1054 if(ipar < 0 || ipar >= fNpar) {
1055 Warning("IsFixed","Illegal parameter number :%d",ipar);
1056 return kFALSE;
1057 }
1058 if (fPL0[ipar] < 0) return kTRUE;
1059 else return kFALSE;
1060}
1061
1062////////////////////////////////////////////////////////////////////////////////
1063/// Main minimization procedure
1064///
1065/// This function is called after setting theoretical function
1066/// by means of TFumili::SetUserFunc and initializing parameters.
1067/// Optionally one can set FCN function (see TFumili::SetFCN and TFumili::Eval)
1068/// If FCN is undefined then user has to provide data arrays by calling
1069/// TFumili::SetData procedure.
1070///
1071/// TFumili::Minimize return following values:
1072/// - 0 - fit is converged
1073/// - -2 - function is not decreasing (or bad derivatives)
1074/// - -3 - error estimations are infinite
1075/// - -4 - maximum number of iterations is exceeded
1076
1078{
1079 Int_t i;
1080 // Flag3 - is fit is chi2 or likelihood? 0 - chi2, 1 - likelihood
1081 fINDFLG[2]=0;
1082 //
1083 // Are the parameters outside of the boundaries ?
1084 //
1085 Int_t parn;
1086
1087 if(fFCN) {
1088 Eval(parn,fGr,fS,fA,9); fNfcn++;
1089 }
1090 for( i = 0; i < fNpar; i++) {
1091 if(fA[i] > fAMX[i]) fA[i] = fAMX[i];
1092 if(fA[i] < fAMN[i]) fA[i] = fAMN[i];
1093 }
1094
1095 Int_t nn2, n, fixFLG, ifix1, fi, nn3, nn1, n0;
1096 Double_t t1;
1097 Double_t sp, t, olds=0;
1098 Double_t bi, aiMAX=0, amb;
1099 Double_t afix, sigi, akap;
1100 Double_t alambd, al, bm, abi, abm;
1101 Int_t l1, k, ifix;
1102
1103 nn2=0;
1104
1105 // Number of parameters;
1106 n=fNpar;
1107 fixFLG=0;
1108
1109 // Exit flag
1110 fENDFLG=0;
1111
1112 // Flag2
1113 fINDFLG[1] = 0;
1114 ifix1=-1;
1115 fi=0;
1116 nn3=0;
1117
1118 // Initialize param.step limits
1119 for( i=0; i < n; i++) {
1120 fR[i]=0.;
1121 if ( fEPS > 0.) fParamError[i] = 0.;
1122 fPL[i] = fPL0[i];
1123 }
1124
1125L3: // Start Iteration
1126
1127 nn1 = 1;
1128 t1 = 1.;
1129
1130L4: // New iteration
1131
1132 // fS - objective function value - zero first
1133 fS = 0.;
1134 // n0 - number of variable parameters in fit
1135 n0 = 0;
1136 for( i = 0; i < n; i++) {
1138 if (fPL0[i] > .0) {
1139 n0=n0+1;
1140 // new iteration - new parallelepiped
1141 if (fPL[i] > .0) fPL0[i]=fPL[i];
1142 }
1143 }
1144 Int_t nn0;
1145 // Calculate number of fZ-matrix elements as nn0=1+2+..+n0
1146 nn0 = n0*(n0+1)/2;
1147 // if (nn0 >= 1) ????
1148 // fZ-matrix is initialized
1149 for( i=0; i < nn0; i++) fZ[i]=0.;
1150
1151 // Flag1
1152 fINDFLG[0] = 0;
1153 Int_t ijkl=1;
1154
1155 // Calculate fS - objective function, fGr - gradients, fZ - fZ-matrix
1156 if(fFCN) {
1157 Eval(parn,fGr,fS,fA,2);
1158 fNfcn++;
1159 } else
1160 ijkl = SGZ();
1161 if(!ijkl) return 10;
1162 if (ijkl == -1) fINDFLG[0]=1;
1163
1164 // sp - scaled on fS machine precision
1165 sp=fRP*TMath::Abs(fS);
1166
1167 // save fZ-matrix
1168 for( i=0; i < nn0; i++) fZ0[i] = fZ[i];
1169 if (nn3 > 0) {
1170 if (nn1 <= fNstepDec) {
1171 t=2.*(fS-olds-fGT);
1172 if (fINDFLG[0] == 0) {
1173 if (TMath::Abs(fS-olds) <= sp && -fGT <= sp) goto L19;
1174 if( 0.59*t < -fGT) goto L19;
1175 t = -fGT/t;
1176 if (t < 0.25 ) t = 0.25;
1177 }
1178 else t = 0.25;
1179 fGT = fGT*t;
1180 t1 = t1*t;
1181 nn2=0;
1182 for( i = 0; i < n; i++) {
1183 if (fPL[i] > 0.) {
1184 fA[i]=fA[i]-fDA[i];
1185 fPL[i]=fPL[i]*t;
1186 fDA[i]=fDA[i]*t;
1187 fA[i]=fA[i]+fDA[i];
1188 }
1189 }
1190 nn1=nn1+1;
1191 goto L4;
1192 }
1193 }
1194
1195L19:
1196
1197 if(fINDFLG[0] != 0) {
1198 fENDFLG=-4;
1199 printf("trying to execute an illegal jump at L85\n");
1200 //goto L85;
1201 }
1202
1203
1204 Int_t k1, k2, i1, j, l;
1205 k1 = 1;
1206 k2 = 1;
1207 i1 = 1;
1208 // In this cycle we removed from fZ contributions from fixed parameters
1209 // We'll get fixed parameters after boundary check
1210 for( i = 0; i < n; i++) {
1211 if (fPL0[i] > .0) {
1212 // if parameter was fixed - release it
1213 if (fPL[i] == 0.) fPL[i]=fPL0[i];
1214 if (fPL[i] > .0) { // ??? it is already non-zero
1215 // if derivative is negative and we above maximum
1216 // or vice versa then fix parameter again and increment k1 by i1
1217 if ((fA[i] >= fAMX[i] && fGr[i] < 0.) ||
1218 (fA[i] <= fAMN[i] && fGr[i] > 0.)) {
1219 fPL[i] = 0.;
1220 k1 = k1 + i1; // i1 stands for fZ-matrix row-number multiplier
1221 /// - skip this row
1222 // in case we are fixing parameter number i
1223 } else {
1224 for( j=0; j <= i; j++) {// cycle on columns of fZ-matrix
1225 if (fPL0[j] > .0) {
1226 // if parameter is not fixed then fZ = fZ0
1227 // Now matrix fZ of other dimension
1228 if (fPL[j] > .0) {
1229 fZ[k2 -1] = fZ0[k1 -1];
1230 k2=k2+1;
1231 }
1232 k1=k1+1;
1233 }
1234 }
1235 }
1236 }
1237 else k1 = k1 + i1; // In case of negative fPL[i] - after mconvd
1238 i1=i1+1; // Next row of fZ0
1239 }
1240 }
1241
1242 // INVERT fZ-matrix (mconvd() procedure)
1243 i1 = 1;
1244 l = 1;
1245 for( i = 0; i < n; i++) {// extract diagonal elements to fR-vector
1246 if (fPL[i] > .0) {
1247 fR[i] = fZ[l - 1];
1248 i1 = i1+1;
1249 l = l + i1;
1250 }
1251 }
1252
1253 n0 = i1 - 1;
1254 InvertZ(n0);
1255
1256 // fZ matrix now is inverted
1257 if (fINDFLG[0] != 0) { // problems
1258 // some PLs now have negative values, try to reduce fZ-matrix again
1259 fINDFLG[0] = 0;
1260 fINDFLG[1] = 1; // errors can be infinite
1261 fixFLG = fixFLG + 1;
1262 fi = 0;
1263 goto L19;
1264 }
1265
1266 // ... CALCULATE THEORETICAL STEP TO MINIMUM
1267 i1 = 1;
1268 for( i = 0; i < n; i++) {
1269 fDA[i]=0.; // initial step is zero
1270 if (fPL[i] > .0) { // for non-fixed parameters
1271 l1=1;
1272 for( l = 0; l < n; l++) {
1273 if (fPL[l] > .0) {
1274 // Calculate offset of Z^-1(i1,l1) element in packed matrix
1275 // because we skip fixed param numbers we need also i,l
1276 if (i1 <= l1 ) k=l1*(l1-1)/2+i1;
1277 else k=i1*(i1-1)/2+l1;
1278 // dA_i = \sum (-Z^{-1}_{il}*grad(fS)_l)
1279 fDA[i]=fDA[i]-fGr[l]*fZ[k - 1];
1280 l1=l1+1;
1281 }
1282 }
1283 i1=i1+1;
1284 }
1285 }
1286 // ... CHECK FOR PARAMETERS ON BOUNDARY
1287
1288 afix=0.;
1289 ifix = -1;
1290 i1 = 1;
1291 l = i1;
1292 for( i = 0; i < n; i++)
1293 if (fPL[i] > .0) {
1294 sigi = TMath::Sqrt(TMath::Abs(fZ[l - 1])); // calculate \sqrt{Z^{-1}_{ii}}
1295 fR[i] = fR[i]*fZ[l - 1]; // Z_ii * Z^-1_ii
1296 if (fEPS > .0) fParamError[i]=sigi;
1297 if ((fA[i] >= fAMX[i] && fDA[i] > 0.) || (fA[i] <= fAMN[i]
1298 && fDA[i] < .0)) {
1299 // if parameter out of bounds and if step is making things worse
1300
1301 akap = TMath::Abs(fDA[i]/sigi);
1302 // let's found maximum of dA/sigi - the worst of parameter steps
1303 if (akap > afix) {
1304 afix=akap;
1305 ifix=i;
1306 ifix1=i;
1307 }
1308 }
1309 i1=i1+1;
1310 l=l+i1;
1311 }
1312 if (ifix != -1) {
1313 // so the worst parameter is found - fix it and exclude,
1314 // reduce fZ-matrix again
1315 fPL[ifix] = -1.;
1316 fixFLG = fixFLG + 1;
1317 fi = 0;
1318 //.. REPEAT CALCULATION OF THEORETICAL STEP AFTER fiXING EACH PARAMETER
1319 goto L19;
1320 }
1321
1322 //... CALCULATE STEP CORRECTION FACTOR
1323
1324 alambd = 1.;
1325 fAKAPPA = 0.;
1326 Int_t imax;
1327 imax = -1;
1328
1329
1330 for( i = 0; i < n; i++) {
1331 if (fPL[i] > .0) {
1332 bm = fAMX[i] - fA[i];
1333 abi = fA[i] + fPL[i]; // upper parameter limit
1334 abm = fAMX[i];
1335 if (fDA[i] <= .0) {
1336 bm = fA[i] - fAMN[i];
1337 abi = fA[i] - fPL[i]; // lower parameter limit
1338 abm = fAMN[i];
1339 }
1340 bi = fPL[i];
1341 // if parallelepiped boundary is crossing limits
1342 // then reduce it (deforming)
1343 if ( bi > bm) {
1344 bi = bm;
1345 abi = abm;
1346 }
1347 // if calculated step is out of bounds
1348 if ( TMath::Abs(fDA[i]) > bi) {
1349 // decrease step splitter alambdA if needed
1350 al = TMath::Abs(bi/fDA[i]);
1351 if (alambd > al) {
1352 imax=i;
1353 aiMAX=abi;
1354 alambd=al;
1355 }
1356 }
1357 // fAKAPPA - parameter will be <fEPS if fit is converged
1358 akap = TMath::Abs(fDA[i]/fParamError[i]);
1359 if (akap > fAKAPPA) fAKAPPA=akap;
1360 }
1361 }
1362 //... CALCULATE NEW CORRECTED STEP
1363 fGT = 0.;
1364 amb = 1.e18;
1365 // alambd - multiplier to split theoretical step dA
1366 if (alambd > .0) amb = 0.25/alambd;
1367 for( i = 0; i < n; i++) {
1368 if (fPL[i] > .0) {
1369 if (nn2 > fNlimMul ) {
1370 if (TMath::Abs(fDA[i]/fPL[i]) > amb ) {
1371 fPL[i] = 4.*fPL[i]; // increase parallelepiped
1372 t1=4.; // flag - that fPL was increased
1373 }
1374 }
1375 // cut step
1376 fDA[i] = fDA[i]*alambd;
1377 // expected function value change in next iteration
1378 fGT = fGT + fDA[i]*fGr[i];
1379 }
1380 }
1381
1382 //.. CHECK IF MINIMUM ATTAINED AND SET EXIT MODE
1383 // if expected fGT smaller than precision
1384 // and other stuff
1385 if (-fGT <= sp && t1 < 1. && alambd < 1.)fENDFLG = -1; // function is not decreasing
1386
1387 if (fENDFLG >= 0) {
1388 if (fAKAPPA < TMath::Abs(fEPS)) { // fit is converging
1389 if (fixFLG == 0)
1390 fENDFLG=1; // successful fit
1391 else {// we have fixed parameters
1392 if (fENDFLG == 0) {
1393 //... CHECK IF fiXING ON BOUND IS CORRECT
1394 fENDFLG = 1;
1395 fixFLG = 0;
1396 ifix1=-1;
1397 // release fixed parameters
1398 for( i = 0; i < fNpar; i++) fPL[i] = fPL0[i];
1399 fINDFLG[1] = 0;
1400 // and repeat iteration
1401 goto L19;
1402 } else {
1403 if( ifix1 >= 0) {
1404 fi = fi + 1;
1405 fENDFLG = 0;
1406 }
1407 }
1408 }
1409 } else { // fit does not converge
1410 if( fixFLG != 0) {
1411 if( fi > fixFLG ) {
1412 //... CHECK IF fiXING ON BOUND IS CORRECT
1413 fENDFLG = 1;
1414 fixFLG = 0;
1415 ifix1=-1;
1416 for( i = 0; i < fNpar; i++) fPL[i] = fPL0[i];
1417 fINDFLG[1] = 0;
1418 goto L19;
1419 } else {
1420 fi = fi + 1;
1421 fENDFLG = 0;
1422 }
1423 } else {
1424 fi = fi + 1;
1425 fENDFLG = 0;
1426 }
1427 }
1428 }
1429
1430// L85:
1431 // iteration number limit is exceeded
1432 if(fENDFLG == 0 && nn3 >= fNmaxIter) fENDFLG=-3;
1433
1434 // fit errors are infinite;
1435 if(fENDFLG > 0 && fINDFLG[1] > 0) fENDFLG=-2;
1436
1437 //MONITO (fS,fNpar,nn3,IT,fEPS,fGT,fAKAPPA,alambd);
1438 if (fENDFLG == 0) { // make step
1439 for ( i = 0; i < n; i++) fA[i] = fA[i] + fDA[i];
1440 if (imax >= 0) fA[imax] = aiMAX;
1441 olds=fS;
1442 nn2=nn2+1;
1443 nn3=nn3+1;
1444 } else {
1445 // fill covariant matrix VL
1446 // fill parameter error matrix up
1447 Int_t il;
1448 il = 0;
1449 for( Int_t ip = 0; ip < fNpar; ip++) {
1450 if( fPL0[ip] > .0) {
1451 for( Int_t jp = 0; jp <= ip; jp++) {
1452 if(fPL0[jp] > .0) {
1453 // VL[ind(ip,jp)] = fZ[il];
1454 il = il + 1;
1455 }
1456 }
1457 }
1458 }
1459 return fENDFLG - 1;
1460 }
1461 goto L3;
1462}
1463
1464////////////////////////////////////////////////////////////////////////////////
1465/// Prints fit results.
1466///
1467/// ikode is the type of printing parameters
1468/// p is function value
1469///
1470/// - ikode = 1 - print values, errors and limits
1471/// - ikode = 2 - print values, errors and steps
1472/// - ikode = 3 - print values, errors, steps and derivatives
1473/// - ikode = 4 - print only values and errors
1474
1476{
1477 TString exitStatus="";
1478 TString xsexpl="";
1479 TString colhdu[3],colhdl[3],cx2,cx3;
1480 switch (fENDFLG) {
1481 case 1:
1482 exitStatus="CONVERGED";
1483 break;
1484 case -1:
1485 exitStatus="CONST FCN";
1486 xsexpl="****\n* FUNCTION IS NOT DECREASING OR BAD DERIVATIVES\n****";
1487 break;
1488 case -2:
1489 exitStatus="ERRORS INF";
1490 xsexpl="****\n* ESTIMATED ERRORS ARE INfiNITE\n****";
1491 break;
1492 case -3:
1493 exitStatus="MAX ITER.";
1494 xsexpl="****\n* MAXIMUM NUMBER OF ITERATIONS IS EXCEEDED\n****";
1495 break;
1496 case -4:
1497 exitStatus="ZERO PROBAB";
1498 xsexpl="****\n* PROBABILITY OF LIKLIHOOD FUNCTION IS NEGATIVE OR ZERO\n****";
1499 break;
1500 default:
1501 exitStatus="UNDEfiNED";
1502 xsexpl="****\n* fiT IS IN PROGRESS\n****";
1503 break;
1504 }
1505 if (ikode == 1) {
1506 colhdu[0] = " ";
1507 colhdl[0] = " ERROR ";
1508 colhdu[1] = " PHYSICAL";
1509 colhdu[2] = " LIMITS ";
1510 colhdl[1] = " NEGATIVE ";
1511 colhdl[2] = " POSITIVE ";
1512 }
1513 if (ikode == 2) {
1514 colhdu[0] = " ";
1515 colhdl[0] = " ERROR ";
1516 colhdu[1] = " INTERNAL ";
1517 colhdl[1] = " STEP SIZE ";
1518 colhdu[2] = " INTERNAL ";
1519 colhdl[2] = " VALUE ";
1520 }
1521 if (ikode == 3) {
1522 colhdu[0] = " ";
1523 colhdl[0] = " ERROR ";
1524 colhdu[1] = " STEP ";
1525 colhdl[1] = " SIZE ";
1526 colhdu[2] = " fiRST ";
1527 colhdl[2] = " DERIVATIVE";
1528 }
1529 if (ikode == 4) {
1530 colhdu[0] = " PARABOLIC ";
1531 colhdl[0] = " ERROR ";
1532 colhdu[1] = " MINOS ";
1533 colhdu[2] = "ERRORS ";
1534 colhdl[1] = " NEGATIVE ";
1535 colhdl[2] = " POSITIVE ";
1536 }
1537 if(fENDFLG<1)Printf("%s", xsexpl.Data());
1538 Printf(" FCN=%g FROM FUMILI STATUS=%-10s %9d CALLS OF FCN",
1539 p,exitStatus.Data(),fNfcn);
1540 Printf(" EDM=%g ",-fGT);
1541 Printf(" EXT PARAMETER %-14s%-14s%-14s",
1542 (const char*)colhdu[0].Data()
1543 ,(const char*)colhdu[1].Data()
1544 ,(const char*)colhdu[2].Data());
1545 Printf(" NO. NAME VALUE %-14s%-14s%-14s",
1546 (const char*)colhdl[0].Data()
1547 ,(const char*)colhdl[1].Data()
1548 ,(const char*)colhdl[2].Data());
1549
1550 for (Int_t i=0;i<fNpar;i++) {
1551 if (ikode==3) {
1552 cx2 = Form("%14.5e",fDA[i]);
1553 cx3 = Form("%14.5e",fGr[i]);
1554
1555 }
1556 if (ikode==1) {
1557 cx2 = Form("%14.5e",fAMN[i]);
1558 cx3 = Form("%14.5e",fAMX[i]);
1559 }
1560 if (ikode==2) {
1561 cx2 = Form("%14.5e",fDA[i]);
1562 cx3 = Form("%14.5e",fA[i]);
1563 }
1564 if(ikode==4){
1565 cx2 = " *undefined* ";
1566 cx3 = " *undefined* ";
1567 }
1568 if(fPL0[i]<=0.) { cx2=" *fixed* ";cx3=""; }
1569 Printf("%4d %-11s%14.5e%14.5e%-14s%-14s",i+1
1570 ,fANames[i].Data(),fA[i],fParamError[i]
1571 ,cx2.Data(),cx3.Data());
1572 }
1573}
1574
1575////////////////////////////////////////////////////////////////////////////////
1576/// Releases parameter number ipar
1577
1579 if(ipar>=0 && ipar<fNpar && fPL0[ipar]<=0.) {
1580 fPL0[ipar] = -fPL0[ipar];
1581 if (fPL0[ipar] == 0. || fPL0[ipar]>=1.) fPL0[ipar]=.1;
1582 }
1583}
1584
1585////////////////////////////////////////////////////////////////////////////////
1586/// Sets pointer to data array provided by user.
1587/// Necessary if SetFCN is not called.
1588///
1589/// - numpoints: number of experimental points
1590/// - vecsize: size of data point vector + 2
1591/// (for N-dimensional fit vecsize=N+2)
1592/// - exdata: data array with following format
1593///
1594/// - exdata[0] = ExpValue_0 - experimental data value number 0
1595/// - exdata[1] = ExpSigma_0 - error of value number 0
1596/// - exdata[2] = X_0[0]
1597/// - exdata[3] = X_0[1]
1598///
1599/// - exdata[vecsize-1] = X_0[vecsize-3]
1600/// - exdata[vecsize] = ExpValue_1
1601/// - exdata[vecsize+1] = ExpSigma_1
1602/// - exdata[vecsize+2] = X_1[0]
1603///
1604/// - exdata[vecsize*(numpoints-1)] = ExpValue_(numpoints-1)
1605///
1606/// - exdata[vecsize*numpoints-1] = X_(numpoints-1)[vecsize-3]
1607
1608void TFumili::SetData(Double_t *exdata,Int_t numpoints,Int_t vecsize){
1609 if(exdata){
1610 fNED1 = numpoints;
1611 fNED2 = vecsize;
1612 fEXDA = exdata;
1613 }
1614}
1615
1616
1617////////////////////////////////////////////////////////////////////////////////
1618/// ret fit method (chisquare or log-likelihood)
1619
1621{
1622 if (!strcmp(name,"H1FitChisquare")) SetFCN(H1FitChisquareFumili);
1623 if (!strcmp(name,"H1FitLikelihood")) SetFCN(H1FitLikelihoodFumili);
1624 if (!strcmp(name,"GraphFitChisquare")) SetFCN(GraphFitChisquareFumili);
1625}
1626
1627////////////////////////////////////////////////////////////////////////////////
1628/// Sets for parameter number ipar initial parameter value,
1629/// name parname, initial error verr and limits vlow and vhigh
1630/// - If vlow = vhigh but not equil to zero, parameter will be fixed.
1631/// - If vlow = vhigh = 0, parameter is released and its limits are discarded
1632
1633Int_t TFumili::SetParameter(Int_t ipar,const char *parname,Double_t value,Double_t verr,Double_t vlow, Double_t vhigh) {
1634 if (ipar<0 || ipar>=fNpar) return -1;
1635 fANames[ipar] = parname;
1636 fA[ipar] = value;
1637 fParamError[ipar] = verr;
1638 if(vlow<vhigh) {
1639 fAMN[ipar] = vlow;
1640 fAMX[ipar] = vhigh;
1641 } else {
1642 if(vhigh<vlow) {
1643 fAMN[ipar] = vhigh;
1644 fAMX[ipar] = vlow;
1645 }
1646 if(vhigh==vlow) {
1647 if(vhigh==0.) {
1648 ReleaseParameter(ipar);
1649 fAMN[ipar] = gMINDOUBLE;
1650 fAMX[ipar] = gMAXDOUBLE;
1651 }
1652 if(vlow!=0) FixParameter(ipar);
1653 }
1654 }
1655 return 0;
1656}
1657
1658////////////////////////////////////////////////////////////////////////////////
1659/// Evaluates objective function ( chi-square ), gradients and
1660/// Z-matrix using data provided by user via TFumili::SetData
1661
1663{
1664 fS = 0.;
1665 Int_t i,j,l,k2=1,k1,ki=0;
1666 Double_t *x = new Double_t[fNED2];
1667 Double_t *df = new Double_t[fNpar];
1668 Int_t nx = fNED2-2;
1669 for (l=0;l<fNED1;l++) { // cycle on all exp. points
1670 k1 = k2;
1671 if (fLogLike) {
1673 nx = fNED2;
1674 k1 -= 2;
1675 }
1676
1677 for (i=0;i<nx;i++){
1678 ki += 1+i;
1679 x[i] = fEXDA[ki];
1680 }
1681 // Double_t y = ARITHM(df,x);
1682 Double_t y = EvalTFN(df,x);
1684 Double_t sig=1.;
1685 if(fLogLike) { // Likelihood method
1686 if(y>0.) {
1687 fS = fS - log(y);
1688 y = -y;
1689 sig= y;
1690 } else { //
1691 delete [] x;
1692 delete [] df;
1693 fS = 1e10;
1694 return -1; // indflg[0] = 1;
1695 }
1696 } else { // Chi2 method
1697 sig = fEXDA[k2]; // sigma of experimental point
1698 y = y - fEXDA[k1-1]; // f(x_i) - F_i
1699 fS = fS + (y*y/(sig*sig))*.5; // simple chi2/2
1700 }
1701 Int_t n = 0;
1702 for (i=0;i<fNpar;i++) {
1703 if (fPL0[i]>0){
1704 df[n] = df[i]/sig; // left only non-fixed param derivatives div by Sig
1705 fGr[i] += df[n]*(y/sig);
1706 n++;
1707 }
1708 }
1709 l = 0;
1710 for (i=0;i<n;i++)
1711 for (j=0;j<=i;j++)
1712 fZ[l++] += df[i]*df[j];
1713 k2 += fNED2;
1714 }
1715
1716 delete[] df;
1717 delete[] x;
1718 return 1;
1719}
1720
1721
1722////////////////////////////////////////////////////////////////////////////////
1723/// Minimization function for H1s using a Chisquare method.
1724/// Default method (function evaluated at center of bin)
1725/// for each point the cache contains the following info
1726/// - 1D : bc,e,xc (bin content, error, x of center of bin)
1727/// - 2D : bc,e,xc,yc
1728/// - 3D : bc,e,xc,yc,zc
1729
1731{
1732 Foption_t fitOption = GetFitOption();
1733 if (fitOption.Integral) {
1734 FitChisquareI(npar,gin,f,u,flag);
1735 return;
1736 }
1737 Double_t cu,eu,fu,fsum;
1738 Double_t x[3];
1739 Double_t *zik=0;
1740 Double_t *pl0=0;
1741
1742 TH1 *hfit = (TH1*)GetObjectFit();
1743 TF1 *f1 = (TF1*)GetUserFunc();
1744 Int_t nd = hfit->GetDimension();
1745 Int_t j;
1746
1747 npar = f1->GetNpar();
1748 SetParNumber(npar);
1749 if(flag == 9) return;
1750 zik = GetZ();
1751 pl0 = GetPL0();
1752
1753 Double_t *df = new Double_t[npar];
1754 f1->InitArgs(x,u);
1755 f = 0;
1756
1757 Int_t npfit = 0;
1758 Double_t *cache = fCache;
1759 for (Int_t i=0;i<fNpoints;i++) {
1760 if (nd > 2) x[2] = cache[4];
1761 if (nd > 1) x[1] = cache[3];
1762 x[0] = cache[2];
1763 cu = cache[0];
1765 fu = f1->EvalPar(x,u);
1766 if (TF1::RejectedPoint()) {cache += fPointSize; continue;}
1767 eu = cache[1];
1768 Derivatives(df,x);
1769 Int_t n = 0;
1770 fsum = (fu-cu)/eu;
1771 if (flag!=1) {
1772 for (j=0;j<npar;j++) {
1773 if (pl0[j]>0) {
1774 df[n] = df[j]/eu;
1775 // left only non-fixed param derivatives / by Sigma
1776 gin[j] += df[n]*fsum;
1777 n++;
1778 }
1779 }
1780 Int_t l = 0;
1781 for (j=0;j<n;j++)
1782 for (Int_t k=0;k<=j;k++)
1783 zik[l++] += df[j]*df[k];
1784 }
1785 f += .5*fsum*fsum;
1786 npfit++;
1787 cache += fPointSize;
1788 }
1789 f1->SetNumberFitPoints(npfit);
1790 delete [] df;
1791}
1792
1793////////////////////////////////////////////////////////////////////////////////
1794/// Minimization function for H1s using a Chisquare method.
1795/// The "I"ntegral method is used
1796/// for each point the cache contains the following info
1797/// - 1D : bc,e,xc,xw (bin content, error, x of center of bin, x bin width of bin)
1798/// - 2D : bc,e,xc,xw,yc,yw
1799/// - 3D : bc,e,xc,xw,yc,yw,zc,zw
1800
1802{
1803 Double_t cu,eu,fu,fsum;
1804 Double_t x[3];
1805 Double_t *zik=0;
1806 Double_t *pl0=0;
1807
1808 TH1 *hfit = (TH1*)GetObjectFit();
1809 TF1 *f1 = (TF1*)GetUserFunc();
1810 Int_t nd = hfit->GetDimension();
1811 Int_t j;
1812
1813 f1->InitArgs(x,u);
1814 npar = f1->GetNpar();
1815 SetParNumber(npar);
1816 if(flag == 9) return;
1817 zik = GetZ();
1818 pl0 = GetPL0();
1819
1820 Double_t *df=new Double_t[npar];
1821 f = 0;
1822
1823 Int_t npfit = 0;
1824 Double_t *cache = fCache;
1825 for (Int_t i=0;i<fNpoints;i++) {
1826 cu = cache[0];
1828 f1->SetParameters(u);
1829 if (nd < 2) {
1830 fu = f1->Integral(cache[2] - 0.5*cache[3],cache[2] + 0.5*cache[3])/cache[3];
1831 } else if (nd < 3) {
1832 fu = ((TF2*)f1)->Integral(cache[2] - 0.5*cache[3],cache[2] + 0.5*cache[3],cache[4] - 0.5*cache[5],cache[4] + 0.5*cache[5])/(cache[3]*cache[5]);
1833 } else {
1834 fu = ((TF3*)f1)->Integral(cache[2] - 0.5*cache[3],cache[2] + 0.5*cache[3],cache[4] - 0.5*cache[5],cache[4] + 0.5*cache[5],cache[6] - 0.5*cache[7],cache[6] + 0.5*cache[7])/(cache[3]*cache[5]*cache[7]);
1835 }
1836 if (TF1::RejectedPoint()) {cache += fPointSize; continue;}
1837 eu = cache[1];
1838 Derivatives(df,x);
1839 Int_t n = 0;
1840 fsum = (fu-cu)/eu;
1841 if (flag!=1) {
1842 for (j=0;j<npar;j++) {
1843 if (pl0[j]>0){
1844 df[n] = df[j]/eu;
1845 // left only non-fixed param derivatives / by Sigma
1846 gin[j] += df[n]*fsum;
1847 n++;
1848 }
1849 }
1850 Int_t l = 0;
1851 for (j=0;j<n;j++)
1852 for (Int_t k=0;k<=j;k++)
1853 zik[l++] += df[j]*df[k];
1854 }
1855 f += .5*fsum*fsum;
1856 npfit++;
1857 cache += fPointSize;
1858 }
1859 f1->SetNumberFitPoints(npfit);
1860 delete[] df;
1861}
1862
1863////////////////////////////////////////////////////////////////////////////////
1864/// Minimization function for H1s using a Likelihood method.
1865/// Basically, it forms the likelihood by determining the Poisson
1866/// probability that given a number of entries in a particular bin,
1867/// the fit would predict it's value. This is then done for each bin,
1868/// and the sum of the logs is taken as the likelihood.
1869///
1870/// Default method (function evaluated at center of bin)
1871/// for each point the cache contains the following info
1872/// - 1D : bc,e,xc (bin content, error, x of center of bin)
1873/// - 2D : bc,e,xc,yc
1874/// - 3D : bc,e,xc,yc,zc
1875
1877{
1878 Foption_t fitOption = GetFitOption();
1879 if (fitOption.Integral) {
1880 FitLikelihoodI(npar,gin,f,u,flag);
1881 return;
1882 }
1883 Double_t cu,fu,fobs,fsub;
1884 Double_t dersum[100];
1885 Double_t x[3];
1886 Int_t icu;
1887
1888 TH1 *hfit = (TH1*)GetObjectFit();
1889 TF1 *f1 = (TF1*)GetUserFunc();
1890 Int_t nd = hfit->GetDimension();
1891 Int_t j;
1892 Double_t *zik = GetZ();
1893 Double_t *pl0 = GetPL0();
1894
1895 npar = f1->GetNpar();
1896 SetParNumber(npar);
1897 if(flag == 9) return;
1898 Double_t *df=new Double_t[npar];
1899 if (flag == 2) for (j=0;j<npar;j++) dersum[j] = gin[j] = 0;
1900 f1->InitArgs(x,u);
1901 f = 0;
1902
1903 Int_t npfit = 0;
1904 Double_t *cache = fCache;
1905 for (Int_t i=0;i<fNpoints;i++) {
1906 if (nd > 2) x[2] = cache[4];
1907 if (nd > 1) x[1] = cache[3];
1908 x[0] = cache[2];
1909 cu = cache[0];
1911 fu = f1->EvalPar(x,u);
1912 if (TF1::RejectedPoint()) {cache += fPointSize; continue;}
1913 if (flag == 2) {
1914 for (j=0;j<npar;j++) {
1915 dersum[j] += 1; //should be the derivative
1916 //grad[j] += dersum[j]*(fu-cu)/eu; dersum[j] = 0;
1917 }
1918 }
1919 if (fu < 1.e-9) fu = 1.e-9;
1920 icu = Int_t(cu);
1921 fsub = -fu +icu*TMath::Log(fu);
1922 fobs = GetSumLog(icu);
1923 fsub -= fobs;
1924 Derivatives(df,x);
1925 int n=0;
1926 // Here we need gradients of Log likelihood function
1927 //
1928 for (j=0;j<npar;j++) {
1929 if (pl0[j]>0){
1930 df[n] = df[j]*(icu/fu-1);
1931 gin[j] -= df[n];
1932 n++;
1933 }
1934 }
1935 Int_t l = 0;
1936 // Z-matrix here - production of first derivatives
1937 // of log-likelihood function
1938 for (j=0;j<n;j++)
1939 for (Int_t k=0;k<=j;k++)
1940 zik[l++] += df[j]*df[k];
1941
1942 f -= fsub;
1943 npfit++;
1944 cache += fPointSize;
1945 }
1946 f *= 2;
1947 f1->SetNumberFitPoints(npfit);
1948 delete[] df;
1949}
1950
1951////////////////////////////////////////////////////////////////////////////////
1952/// Minimization function for H1s using a Likelihood method.
1953/// Basically, it forms the likelihood by determining the Poisson
1954/// probability that given a number of entries in a particular bin,
1955/// the fit would predict it's value. This is then done for each bin,
1956/// and the sum of the logs is taken as the likelihood.
1957///
1958/// The "I"ntegral method is used
1959/// for each point the cache contains the following info
1960/// - 1D : bc,e,xc,xw (bin content, error, x of center of bin, x bin width of bin)
1961/// - 2D : bc,e,xc,xw,yc,yw
1962/// - 3D : bc,e,xc,xw,yc,yw,zc,zw
1963
1965{
1966 Double_t cu,fu,fobs,fsub;
1967 Double_t dersum[100];
1968 Double_t x[3];
1969 Int_t icu;
1970
1971 TH1 *hfit = (TH1*)GetObjectFit();
1972 TF1 *f1 = (TF1*)GetUserFunc();
1973 Int_t nd = hfit->GetDimension();
1974 Int_t j;
1975 Double_t *zik = GetZ();
1976 Double_t *pl0 = GetPL0();
1977
1978 Double_t *df=new Double_t[npar];
1979
1980 npar = f1->GetNpar();
1981 SetParNumber(npar);
1982 if(flag == 9) {delete [] df; return;}
1983 if (flag == 2) for (j=0;j<npar;j++) dersum[j] = gin[j] = 0;
1984 f1->InitArgs(x,u);
1985 f = 0;
1986
1987 Int_t npfit = 0;
1988 Double_t *cache = fCache;
1989 for (Int_t i=0;i<fNpoints;i++) {
1990 if (nd > 2) x[2] = cache[4];
1991 if (nd > 1) x[1] = cache[3];
1992 x[0] = cache[2];
1993 cu = cache[0];
1995 if (nd < 2) {
1996 fu = f1->Integral(cache[2] - 0.5*cache[3],cache[2] + 0.5*cache[3])/cache[3];
1997 } else if (nd < 3) {
1998 fu = ((TF2*)f1)->Integral(cache[2] - 0.5*cache[3],cache[2] + 0.5*cache[3],cache[4] - 0.5*cache[5],cache[4] + 0.5*cache[5])/(cache[3]*cache[5]);
1999 } else {
2000 fu = ((TF3*)f1)->Integral(cache[2] - 0.5*cache[3],cache[2] + 0.5*cache[3],cache[4] - 0.5*cache[5],cache[4] + 0.5*cache[5],cache[6] - 0.5*cache[7],cache[6] + 0.5*cache[7])/(cache[3]*cache[5]*cache[7]);
2001 }
2002 if (TF1::RejectedPoint()) {cache += fPointSize; continue;}
2003 if (flag == 2) {
2004 for (j=0;j<npar;j++) {
2005 dersum[j] += 1; //should be the derivative
2006 //grad[j] += dersum[j]*(fu-cu)/eu; dersum[j] = 0;
2007 }
2008 }
2009 if (fu < 1.e-9) fu = 1.e-9;
2010 icu = Int_t(cu);
2011 fsub = -fu +icu*TMath::Log(fu);
2012 fobs = GetSumLog(icu);
2013 fsub -= fobs;
2014 Derivatives(df,x);
2015 int n=0;
2016 // Here we need gradients of Log likelihood function
2017 //
2018 for (j=0;j<npar;j++) {
2019 if (pl0[j]>0){
2020 df[n] = df[j]*(icu/fu-1);
2021 gin[j] -= df[n];
2022 n++;
2023 }
2024 }
2025 Int_t l = 0;
2026 // Z-matrix here - production of first derivatives
2027 // of log-likelihood function
2028 for (j=0;j<n;j++)
2029 for (Int_t k=0;k<=j;k++)
2030 zik[l++] += df[j]*df[k];
2031
2032 f -= fsub;
2033 npfit++;
2034 cache += fPointSize;
2035 }
2036 f *= 2;
2037 f1->SetNumberFitPoints(npfit);
2038 delete[] df;
2039}
2040
2041
2042//______________________________________________________________________________
2043//
2044// STATIC functions
2045//______________________________________________________________________________
2046
2047////////////////////////////////////////////////////////////////////////////////
2048/// Minimization function for H1s using a Chisquare method.
2049
2051{
2053 hFitter->FitChisquare(npar, gin, f, u, flag);
2054}
2055
2056////////////////////////////////////////////////////////////////////////////////
2057/// Minimization function for H1s using a Likelihood method.
2058/// Basically, it forms the likelihood by determining the Poisson
2059/// probability that given a number of entries in a particular bin,
2060/// the fit would predict it's value. This is then done for each bin,
2061/// and the sum of the logs is taken as the likelihood.
2062/// PDF: P=exp(-f(x_i))/[F_i]!*(f(x_i))^[F_i]
2063/// where F_i - experimental value, f(x_i) - expected theoretical value
2064/// [F_i] - integer part of F_i.
2065/// drawback is that if F_i>Int_t - GetSumLog will fail
2066/// for big F_i is faster to use Euler's Gamma-function
2067
2069{
2070
2072 hFitter->FitLikelihood(npar, gin, f, u, flag);
2073}
2074
2075////////////////////////////////////////////////////////////////////////////////
2076/// Minimization function for Graphs using a Chisquare method.
2077/// In case of a TGraphErrors object, ex, the error along x, is projected
2078/// along the y-direction by calculating the function at the points x-exlow and
2079/// x+exhigh.
2080///
2081/// The chisquare is computed as the sum of the quantity below at each point:
2082///
2083/// (y - f(x))**2
2084/// -----------------------------------
2085/// ey**2 + (0.5*(exl + exh)*f'(x))**2
2086///
2087/// where x and y are the point coordinates and f'(x) is the derivative of function f(x).
2088/// This method to approximate the uncertainty in y because of the errors in x, is called
2089/// "effective variance" method.
2090/// The improvement, compared to the previously used method (f(x+ exhigh) - f(x-exlow))/2
2091/// is of (error of x)**2 order.
2092///
2093/// NOTE:
2094///
2095/// 1. By using the "effective variance" method a simple linear regression
2096/// becomes a non-linear case , which takes several iterations
2097/// instead of 0 as in the linear case .
2098///
2099/// 2. The effective variance technique assumes that there is no correlation
2100/// between the x and y coordinate .
2101///
2102/// In case the function lies below (above) the data point, ey is ey_low (ey_high).
2103
2105 Double_t *u, Int_t flag)
2106{
2107 Double_t cu,eu,exl,exh,ey,eux,fu,fsum;
2108 Double_t x[1];
2109 Int_t i, bin, npfits=0;
2110
2112 TGraph *gr = (TGraph*)grFitter->GetObjectFit();
2113 TF1 *f1 = (TF1*)grFitter->GetUserFunc();
2114 Foption_t fitOption = grFitter->GetFitOption();
2115
2116 Int_t n = gr->GetN();
2117 Double_t *gx = gr->GetX();
2118 Double_t *gy = gr->GetY();
2119 npar = f1->GetNpar();
2120
2121 grFitter->SetParNumber(npar);
2122
2123 if(flag == 9) return;
2124 Double_t *zik = grFitter->GetZ();
2125 Double_t *pl0 = grFitter->GetPL0();
2126 Double_t *df = new Double_t[npar];
2127
2128
2129 f1->InitArgs(x,u);
2130 f = 0;
2131 for (bin=0;bin<n;bin++) {
2132 x[0] = gx[bin];
2133 if (!f1->IsInside(x)) continue;
2134 cu = gy[bin];
2136 fu = f1->EvalPar(x,u);
2137 if (TF1::RejectedPoint()) continue;
2138 npfits++;
2139 Double_t eusq=1.;
2140 if (fitOption.W1) {
2141 eu = 1.;
2142 } else {
2143 exh = gr->GetErrorXhigh(bin);
2144 exl = gr->GetErrorXlow(bin);
2145 ey = gr->GetErrorY(bin);
2146 if (exl < 0) exl = 0;
2147 if (exh < 0) exh = 0;
2148 if (ey < 0) ey = 0;
2149 if (exh > 0 && exl > 0) {
2150// "Effective variance" method for projecting errors
2151 eux = 0.5*(exl + exh)*f1->Derivative(x[0], u);
2152 } else
2153 eux = 0.;
2154 eu = ey*ey+eux*eux;
2155 if (eu <= 0) eu = 1;
2156 eusq = TMath::Sqrt(eu);
2157 }
2158 grFitter->Derivatives(df,x);
2159 n = 0;
2160 fsum = (fu-cu)/eusq;
2161 for (i=0;i<npar;i++) {
2162 if (pl0[i]>0){
2163 df[n] = df[i]/eusq;
2164 // left only non-fixed param derivatives / by Sigma
2165 gin[i] += df[n]*fsum;
2166 n++;
2167 }
2168 }
2169 Int_t l = 0;
2170 for (i=0;i<n;i++)
2171 for (Int_t j=0;j<=i;j++)
2172 zik[l++] += df[i]*df[j];
2173 f += .5*fsum*fsum;
2174
2175 }
2176 delete [] df;
2177 f1->SetNumberFitPoints(npfits);
2178}
2179
#define d(i)
Definition RSha256.hxx:102
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define e(i)
Definition RSha256.hxx:103
int Int_t
Definition RtypesCore.h:45
const Bool_t kFALSE
Definition RtypesCore.h:92
double Double_t
Definition RtypesCore.h:59
const Bool_t kTRUE
Definition RtypesCore.h:91
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
static const Double_t gMAXDOUBLE
Definition TFumili.cxx:124
void H1FitChisquareFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Minimization function for H1s using a Chisquare method.
Definition TFumili.cxx:2050
void GraphFitChisquareFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Minimization function for Graphs using a Chisquare method.
Definition TFumili.cxx:2104
TFumili * gFumili
Definition TFumili.cxx:121
void H1FitLikelihoodFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Minimization function for H1s using a Likelihood method.
Definition TFumili.cxx:2068
static const Double_t gMINDOUBLE
Definition TFumili.cxx:125
R__EXTERN TFumili * gFumili
Definition TFumili.h:117
char name[80]
Definition TGX11.cxx:110
#define hi
double sqrt(double)
double log(double)
#define gROOT
Definition TROOT.h:406
char * Form(const char *fmt,...)
void Printf(const char *fmt,...)
1-Dim function class
Definition TF1.h:213
static void RejectPoint(Bool_t reject=kTRUE)
Static function to set the global flag to reject points the fgRejectPoint global flag is tested by al...
Definition TF1.cxx:3658
virtual Int_t GetNpar() const
Definition TF1.h:481
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition TF1.cxx:2515
virtual Double_t Derivative(Double_t x, Double_t *params=0, Double_t epsilon=0.001) const
Returns the first derivative of the function at point x, computed by Richardson's extrapolation metho...
Definition TF1.cxx:1102
virtual void SetNumberFitPoints(Int_t npfits)
Definition TF1.h:624
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition TF1.cxx:1463
virtual void InitArgs(const Double_t *x, const Double_t *params)
Definition TF1.cxx:2467
static Bool_t RejectedPoint()
See TF1::RejectPoint above.
Definition TF1.cxx:3667
virtual void SetParameters(const Double_t *params)
Definition TF1.h:644
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition TF1.h:598
A 2-Dim function with parameters.
Definition TF2.h:29
A 3-Dim function with parameters.
Definition TF3.h:28
Bool_t fWARN
warnings
Definition TFumili.h:29
Int_t fNED12
fNED1+fNED2
Definition TFumili.h:18
void DeleteArrays()
Deallocates memory. Called from destructor TFumili::~TFumili.
Definition TFumili.cxx:253
Bool_t fNumericDerivatives
Definition TFumili.h:32
virtual Int_t GetNumberFreeParameters() const
Return the number of free parameters.
Definition TFumili.cxx:806
Int_t fNED2
K - Length of vector X plus 2 (for chi2)
Definition TFumili.h:17
virtual Double_t Chisquare(Int_t npar, Double_t *params) const
return a chisquare equivalent
Definition TFumili.cxx:219
Int_t fNpar
fNpar - number of parameters
Definition TFumili.h:19
virtual Double_t GetParError(Int_t ipar) const
Return error of parameter ipar.
Definition TFumili.cxx:818
virtual void PrintResults(Int_t k, Double_t p) const
Prints fit results.
Definition TFumili.cxx:1475
virtual ~TFumili()
TFumili destructor.
Definition TFumili.cxx:209
virtual Int_t GetErrors(Int_t ipar, Double_t &eplus, Double_t &eminus, Double_t &eparab, Double_t &globcc) const
Return errors after MINOs not implemented.
Definition TFumili.cxx:874
virtual void SetFitMethod(const char *name)
ret fit method (chisquare or log-likelihood)
Definition TFumili.cxx:1620
virtual void FitLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Minimization function for H1s using a Likelihood method.
Definition TFumili.cxx:1876
virtual Double_t GetSumLog(Int_t)
Return Sum(log(i) i=0,n used by log-likelihood fits.
Definition TFumili.cxx:913
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)
Execute MINUIT commands.
Definition TFumili.cxx:383
Double_t * fEXDA
[fNED12] experimental data poInt_ter
Definition TFumili.h:41
virtual Double_t GetCovarianceMatrixElement(Int_t i, Int_t j) const
Return element i,j from the covariance matrix.
Definition TFumili.cxx:785
Int_t SGZ()
Evaluates objective function ( chi-square ), gradients and Z-matrix using data provided by user via T...
Definition TFumili.cxx:1662
Int_t fNlog
Definition TFumili.h:14
virtual void FitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Minimization function for H1s using a Chisquare method.
Definition TFumili.cxx:1730
void Derivatives(Double_t *, Double_t *)
Calculates partial derivatives of theoretical function.
Definition TFumili.cxx:283
Double_t * fAMN
[fMaxParam] Minimum param value
Definition TFumili.h:51
TString * fANames
[fMaxParam] Parameter names
Definition TFumili.h:62
virtual void FixParameter(Int_t ipar)
Fixes parameter number ipar.
Definition TFumili.cxx:766
Double_t * GetPL0() const
Definition TFumili.h:95
Double_t * fPL
[fMaxParam] Limits for parameters step. If <0, then parameter is fixed
Definition TFumili.h:46
Int_t Eval(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *par, Int_t flag)
Evaluate the minimisation function.
Definition TFumili.cxx:342
Int_t fMaxParam
Definition TFumili.h:13
void SetParNumber(Int_t ParNum)
Definition TFumili.h:112
void SetData(Double_t *, Int_t, Int_t)
Sets pointer to data array provided by user.
Definition TFumili.cxx:1608
virtual Int_t GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const
Return global fit parameters.
Definition TFumili.cxx:896
Int_t fINDFLG[5]
internal flags;
Definition TFumili.h:25
Double_t EvalTFN(Double_t *, Double_t *)
Evaluate theoretical function.
Definition TFumili.cxx:354
Double_t * GetZ() const
Definition TFumili.h:102
Double_t * fParamError
[fMaxParam] Parameter errors
Definition TFumili.h:39
Int_t fENDFLG
End flag of fit.
Definition TFumili.h:24
Double_t * fR
[fMaxParam] Correlation factors
Definition TFumili.h:52
Double_t * fDA
[fMaxParam] Parameter step
Definition TFumili.h:49
Int_t fNstepDec
fNstepDec - maximum number of step decreasing counter
Definition TFumili.h:20
Double_t * fZ0
[fMaxParam2] Matrix of approximate second derivatives of objective function This matrix is diagonal a...
Definition TFumili.h:34
Double_t * fPL0
[fMaxParam] Step initial bounds
Definition TFumili.h:45
virtual Bool_t IsFixed(Int_t ipar) const
Return kTRUE if parameter ipar is fixed, kFALSE otherwise)
Definition TFumili.cxx:1052
virtual void ReleaseParameter(Int_t ipar)
Releases parameter number ipar.
Definition TFumili.cxx:1578
virtual Double_t * GetCovarianceMatrix() const
Return a pointer to the covariance matrix.
Definition TFumili.cxx:776
Double_t * fA
[fMaxParam] Fit parameter array
Definition TFumili.h:44
Double_t fAKAPPA
Definition TFumili.h:60
Int_t Minimize()
Main minimization procedure.
Definition TFumili.cxx:1077
Int_t fNmaxIter
fNmaxIter - maximum number of iterations
Definition TFumili.h:22
TFumili(Int_t maxpar=25)
Definition TFumili.cxx:129
Int_t ExecuteSetCommand(Int_t)
Called from TFumili::ExecuteCommand in case of "SET xxx" and "SHOW xxx".
Definition TFumili.cxx:550
Double_t fS
fS - objective function value (return)
Definition TFumili.h:57
Double_t fEPS
fEPS - required precision of parameters. If fEPS<0 then
Definition TFumili.h:58
virtual void FitChisquareI(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Minimization function for H1s using a Chisquare method.
Definition TFumili.cxx:1801
Int_t fNfcn
Number of FCN calls;.
Definition TFumili.h:15
Int_t fLastFixed
Last fixed parameter number.
Definition TFumili.h:23
void BuildArrays()
Allocates memory for internal arrays.
Definition TFumili.cxx:174
virtual Int_t GetNumberTotalParameters() const
Return the total number of parameters (free + fixed)
Definition TFumili.cxx:798
Double_t * fZ
[fMaxParam2] Invers fZ0 matrix - covariance matrix
Definition TFumili.h:37
Bool_t fLogLike
LogLikelihood flag.
Definition TFumili.h:31
Int_t fNED1
Number of experimental vectors X=(x1,x2,...xK)
Definition TFumili.h:16
Double_t * fGr
Definition TFumili.h:38
Double_t fGT
Expected function change in next iteration.
Definition TFumili.h:61
virtual Int_t SetParameter(Int_t ipar, const char *parname, Double_t value, Double_t verr, Double_t vlow, Double_t vhigh)
Sets for parameter number ipar initial parameter value, name parname, initial error verr and limits v...
Definition TFumili.cxx:1633
TString fCword
Command string.
Definition TFumili.h:63
virtual void FitLikelihoodI(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Minimization function for H1s using a Likelihood method.
Definition TFumili.cxx:1964
Double_t fRP
Precision of fit ( machine zero on CDC 6000) quite old yeh?
Definition TFumili.h:59
virtual Double_t GetParameter(Int_t ipar) const
Return current value of parameter ipar.
Definition TFumili.cxx:827
Double_t * fCmPar
[fMaxParam] parameters of commands
Definition TFumili.h:55
Double_t * fDF
[fMaxParam] First derivatives of theoretical function
Definition TFumili.h:54
virtual void Clear(Option_t *opt="")
Resets all parameter names, values and errors to zero.
Definition TFumili.cxx:234
virtual const char * GetParName(Int_t ipar) const
Return name of parameter ipar.
Definition TFumili.cxx:864
Double_t * fSumLog
[fNlog]
Definition TFumili.h:40
Double_t * fAMX
[fMaxParam] Maximum param value
Definition TFumili.h:50
Bool_t fDEBUG
debug info
Definition TFumili.h:30
Int_t fNlimMul
fNlimMul - after fNlimMul successful iterations permits four-fold increasing of fPL
Definition TFumili.h:21
Definition TFumili.h:28
void InvertZ(Int_t)
Inverts packed diagonal matrix Z by square-root method.
Definition TFumili.cxx:937
Double_t GetErrorXhigh(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetErrorXlow(Int_t bin) const
This function is called by GraphFitChisquare.
Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
Double_t * GetY() const
Definition TGraph.h:132
Int_t GetN() const
Definition TGraph.h:124
Double_t * GetX() const
Definition TGraph.h:131
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:58
virtual Int_t GetDimension() const
Definition TH1.h:282
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition TNamed.cxx:140
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:879
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:893
@ kInvalidObject
if object ctor succeeded but object should not be used
Definition TObject.h:68
Basic string class.
Definition TString.h:136
const char * Data() const
Definition TString.h:369
void ToUpper()
Change string to upper case.
Definition TString.cxx:1158
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:624
virtual TObject * GetObjectFit() const
TObject * fUserFunc
virtual Foption_t GetFitOption() const
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
To set the address of the minimization objective function called by the native compiler (see function...
Double_t * fCache
void(* fFCN)(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
static TVirtualFitter * GetFitter()
static: return the current Fitter
virtual TObject * GetUserFunc() const
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Double_t ey[n]
Definition legend1.C:17
TGraphErrors * gr
Definition legend1.C:25
TF1 * f1
Definition legend1.C:11
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:212
Double_t Log(Double_t x)
Definition TMath.h:760
Double_t Sqrt(Double_t x)
Definition TMath.h:691
Short_t Abs(Short_t d)
Definition TMathBase.h:120
int W1
Definition Foption.h:36
int Integral
Definition Foption.h:44
auto * l
Definition textangle.C:4
auto * t1
Definition textangle.C:20