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