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