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