ROOT   6.10/09 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
16 likelihood function.
17
18 <p>Experimentally measured values $F_i$ are fitted with theoretical
19 functions $f_i({\vec x}_i,\vec\theta\,\,)$, where ${\vec x}_i$ are
20 coordinates, 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 26 x_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 34 x_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
39 retaining 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
48 case:
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
58 expression is used when last term in (4) is discarded. It is often
59 done, not always wittingly, and sometimes causes troubles, for example,
60 if 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
62 minimize $\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 66 Z_{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
78 restriction. For an initial value of parameter ${\vec\theta}^0$ a
79 parallelepiped $P_0$ is built with the center at ${\vec\theta}^0$ and
80 axes parallel to coordinate axes $\theta_i$. The lengths of
81 parallelepiped sides along i-th axis is $2b_i$, where $b_i$ is such a
82 value that the functions $f_j(\vec\theta)$ are quasi-linear all over
83 the 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$).
91 Very similar step formulae are used in FUMILI for negative logarithm
92 of the likelihood function with the same idea - linearization of
93 function 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
113 extern void H1FitChisquareFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
114 extern void H1FitLikelihoodFumili(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
115 extern 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
123 static const Double_t gMAXDOUBLE=1e300;
124 static 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;
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 ////////////////////////////////////////////////////////////////////////////////
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 {
241  fNpar = fMaxParam;
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
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
353 {
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
394 Int_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
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++)
497  ReleaseParameter(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
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;
705  if(!setCommand) return 0;
707  return 0;
709  if(!setCommand) return 0;
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
860 Int_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
880 const 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
890 Int_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
912 Int_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
1158 L3: // Start Iteration
1159
1160  nn1 = 1;
1161  t1 = 1.;
1162
1163 L4: // 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++) {
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
1228 L19:
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
1645 void 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
1657 void TFumili::SetFitMethod(const char *name)
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
1672 Int_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);
1723  if(fNumericDerivatives) Derivatives(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 /// ======================================================
2091
2093 {
2094  TFumili *hFitter = (TFumili*)TVirtualFitter::GetFitter();
2095  hFitter->FitChisquare(npar, gin, f, u, flag);
2096 }
2097
2098 ////////////////////////////////////////////////////////////////////////////////
2099 /// -*-*-*-*Minimization function for H1s using a Likelihood method*-*-*-*-*-*
2100 /// =======================================================
2101 /// Basically, it forms the likelihood by determining the Poisson
2102 /// probability that given a number of entries in a particular bin,
2103 /// the fit would predict it's value. This is then done for each bin,
2104 /// and the sum of the logs is taken as the likelihood.
2105 /// PDF: P=exp(-f(x_i))/[F_i]!*(f(x_i))^[F_i]
2106 /// where F_i - experimental value, f(x_i) - expected theoretical value
2107 /// [F_i] - integer part of F_i.
2108 /// drawback is that if F_i>Int_t - GetSumLog will fail
2109 /// for big F_i is faster to use Euler's Gamma-function
2110
2112 {
2113
2114  TFumili *hFitter = (TFumili*)TVirtualFitter::GetFitter();
2115  hFitter->FitLikelihood(npar, gin, f, u, flag);
2116 }
2117
2118
2119
2120 ////////////////////////////////////////////////////////////////////////////////
2121 ///*-*-*-*-*-*Minimization function for Graphs using a Chisquare method*-*-*-*-*
2122 ///*-* =========================================================
2123 ///
2124 /// In case of a TGraphErrors object, ex, the error along x, is projected
2125 /// along the y-direction by calculating the function at the points x-exlow and
2126 /// x+exhigh.
2127 ///
2128 /// The chisquare is computed as the sum of the quantity below at each point:
2129 ///
2130 /// (y - f(x))**2
2131 /// -----------------------------------
2132 /// ey**2 + (0.5*(exl + exh)*f'(x))**2
2133 ///
2134 /// where x and y are the point coordinates and f'(x) is the derivative of function f(x).
2135 /// This method to approximate the uncertainty in y because of the errors in x, is called
2136 /// "effective variance" method.
2137 /// The improvement, compared to the previously used method (f(x+ exhigh) - f(x-exlow))/2
2138 /// is of (error of x)**2 order.
2139 /// NOTE:
2140 /// 1) By using the "effective variance" method a simple linear regression
2141 /// becomes a non-linear case , which takes several iterations
2142 /// instead of 0 as in the linear case .
2143 ///
2144 /// 2) The effective variance technique assumes that there is no correlation
2145 /// between the x and y coordinate .
2146 ///
2147 /// In case the function lies below (above) the data point, ey is ey_low (ey_high).
2148
2150  Double_t *u, Int_t flag)
2151 {
2152  Double_t cu,eu,exl,exh,ey,eux,fu,fsum;
2153  Double_t x[1];
2154  Int_t i, bin, npfits=0;
2155
2156  TFumili *grFitter = (TFumili*)TVirtualFitter::GetFitter();
2157  TGraph *gr = (TGraph*)grFitter->GetObjectFit();
2158  TF1 *f1 = (TF1*)grFitter->GetUserFunc();
2159  Foption_t fitOption = grFitter->GetFitOption();
2160
2161  Int_t n = gr->GetN();
2162  Double_t *gx = gr->GetX();
2163  Double_t *gy = gr->GetY();
2164  npar = f1->GetNpar();
2165
2166  grFitter->SetParNumber(npar);
2167
2168  if(flag == 9) return;
2169  Double_t *zik = grFitter->GetZ();
2170  Double_t *pl0 = grFitter->GetPL0();
2171  Double_t *df = new Double_t[npar];
2172
2173
2174  f1->InitArgs(x,u);
2175  f = 0;
2176  for (bin=0;bin<n;bin++) {
2177  x[0] = gx[bin];
2178  if (!f1->IsInside(x)) continue;
2179  cu = gy[bin];
2181  fu = f1->EvalPar(x,u);
2182  if (TF1::RejectedPoint()) continue;
2183  npfits++;
2184  Double_t eusq=1.;
2185  if (fitOption.W1) {
2186  eu = 1.;
2187  } else {
2188  exh = gr->GetErrorXhigh(bin);
2189  exl = gr->GetErrorXlow(bin);
2190  ey = gr->GetErrorY(bin);
2191  if (exl < 0) exl = 0;
2192  if (exh < 0) exh = 0;
2193  if (ey < 0) ey = 0;
2194  if (exh > 0 && exl > 0) {
2195 // "Effective variance" method for projecting errors
2196  eux = 0.5*(exl + exh)*f1->Derivative(x[0], u);
2197  } else
2198  eux = 0.;
2199  eu = ey*ey+eux*eux;
2200  if (eu <= 0) eu = 1;
2201  eusq = TMath::Sqrt(eu);
2202  }
2203  grFitter->Derivatives(df,x);
2204  n = 0;
2205  fsum = (fu-cu)/eusq;
2206  for (i=0;i<npar;i++) {
2207  if (pl0[i]>0){
2208  df[n] = df[i]/eusq;
2209  // left only non-fixed param derivatives / by Sigma
2210  gin[i] += df[n]*fsum;
2211  n++;
2212  }
2213  }
2214  Int_t l = 0;
2215  for (i=0;i<n;i++)
2216  for (Int_t j=0;j<=i;j++)
2217  zik[l++] += df[i]*df[j];
2218  f += .5*fsum*fsum;
2219
2220  }
2221  delete [] df;
2222  f1->SetNumberFitPoints(npfits);
2223 }
2224
const int nx
Definition: kalman.C:16
double par[1]
Definition: unuranDistr.cxx:38
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
Double_t * fCache
Double_t * fDF
Definition: TFumili.h:62
virtual Double_t Chisquare(Int_t npar, Double_t *params) const
return a chisquare equivalent
Definition: TFumili.cxx:222
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:588
virtual Double_t * GetCovarianceMatrix() const
return a pointer to the covariance matrix
Definition: TFumili.cxx:791
virtual Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1406
virtual Double_t GetCovarianceMatrixElement(Int_t i, Int_t j) const
return element i,j from the covariance matrix
Definition: TFumili.cxx:800
Double_t Log(Double_t x)
Definition: TMath.h:649
Double_t * fR
Definition: TFumili.h:60
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
const double pi
const char Option_t
Definition: RtypesCore.h:62
Double_t fRP
Definition: TFumili.h:67
virtual void FixParameter(Int_t ipar)
Fixes parameter number ipar.
Definition: TFumili.cxx:781
Int_t fNlog
Definition: TFumili.h:22
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:131
TFumili(Int_t maxpar=25)
Definition: TFumili.cxx:129
virtual Double_t GetParameter(Int_t ipar) const
return current value of parameter ipar
Definition: TFumili.cxx:843
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1112
#define gROOT
Definition: TROOT.h:375
virtual Int_t GetNumberTotalParameters() const
return the total number of parameters (free + fixed)
Definition: TFumili.cxx:814
Basic string class.
Definition: TString.h:129
virtual void SetNumberFitPoints(Int_t npfits)
Definition: TF1.h:568
Double_t * fDA
Definition: TFumili.h:57
Double_t * fParamError
Definition: TFumili.h:47
virtual void PrintResults(Int_t k, Double_t p) const
Prints fit results.
Definition: TFumili.cxx:1509
int Int_t
Definition: RtypesCore.h:41
Double_t * GetPL0() const
Definition: TFumili.h:103
bool Bool_t
Definition: RtypesCore.h:59
Double_t * fPL
Definition: TFumili.h:54
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
IntegralOneDim or analytical integral.
Definition: TF1.cxx:2309
TLatex * t1
Definition: textangle.C:20
Double_t * fEXDA
Definition: TFumili.h:49
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
Int_t Minimize()
Definition: TFumili.cxx:1089
Bool_t fNumericDerivatives
Definition: TFumili.h:40
Int_t fNED1
Definition: TFumili.h:24
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&#39;s extrapolation metho...
Definition: TF1.cxx:867
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
virtual Int_t GetDimension() const
Definition: TH1.h:263
double sqrt(double)
void InvertZ(Int_t)
Inverts packed diagonal matrix Z by square-root method.
Definition: TFumili.cxx:958
Double_t x[n]
Definition: legend1.C:17
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:2092
virtual const char * GetParName(Int_t ipar) const
return name of parameter ipar
Definition: TFumili.cxx:880
void BuildArrays()
Allocates memory for internal arrays.
Definition: TFumili.cxx:174
TString * fANames
Definition: TFumili.h:70
Int_t fNpar
Definition: TFumili.h:27
static const Double_t gMINDOUBLE
Definition: TFumili.cxx:124
std::vector< std::vector< double > > Data
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
you should not use this method at all Int_t Int_t Double_t bm
Definition: TRolke.cxx:630
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
TString fCword
Definition: TFumili.h:71
Bool_t fDEBUG
Definition: TFumili.h:38
Bool_t fLogLike
Definition: TFumili.h:39
Double_t * fZ0
Definition: TFumili.h:42
virtual Double_t GetParError(Int_t ipar) const
return error of parameter ipar
Definition: TFumili.cxx:834
virtual Foption_t GetFitOption() const
Double_t * fCmPar
Definition: TFumili.h:63
TObject * fUserFunc
Int_t fNED2
Definition: TFumili.h:25
virtual Double_t GetErrorXlow(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1426
Int_t ExecuteSetCommand(Int_t)
Called from TFumili::ExecuteCommand in case of "SET xxx" and "SHOW xxx".
Definition: TFumili.cxx:565
if object ctor succeeded but object should not be used
Definition: TObject.h:65
A 3-Dim function with parameters.
Definition: TF3.h:28
static const Double_t gMAXDOUBLE
Definition: TFumili.cxx:123
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*-*-*-*-*-* Basically, it forms the likelihood by determining the Poisson probability that given a number of entries in a particular bin, the fit would predict it&#39;s value.
Definition: TFumili.cxx:2111
void DeleteArrays()
Deallocates memory.
Definition: TFumili.cxx:261
Double_t fAKAPPA
Definition: TFumili.h:68
virtual Bool_t IsFixed(Int_t ipar) const
return kTRUE if parameter ipar is fixed, kFALSE othersise)
Definition: TFumili.cxx:1076
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:873
char * Form(const char *fmt,...)
static TVirtualFitter * GetFitter()
static: return the current Fitter
Int_t GetN() const
Definition: TGraph.h:122
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:3459
TLine * l
Definition: textangle.C:4
Double_t fS
Definition: TFumili.h:65
virtual Int_t GetNumberFreeParameters() const
return the number of free parameters
Definition: TFumili.cxx:822
int Integral
Definition: Foption.h:44
Int_t fMaxParam
Definition: TFumili.h:21
A 2-Dim function with parameters.
Definition: TF2.h:29
virtual Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs)
Execute MINUIT commands.
Definition: TFumili.cxx:394
Double_t * fGr
Definition: TFumili.h:46
virtual void SetFitMethod(const char *name)
ret fit method (chisquare or loglikelihood)
Definition: TFumili.cxx:1657
TGraphErrors * gr
Definition: legend1.C:25
int W1
Definition: Foption.h:36
#define Printf
Definition: TGeoToOCC.h:18
const Bool_t kFALSE
Definition: RtypesCore.h:92
Double_t * GetX() const
Definition: TGraph.h:129
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 The "I"ntegral method is used for each point t...
Definition: TFumili.cxx:1842
Int_t fNED12
Definition: TFumili.h:26
Bool_t fWARN
Definition: TFumili.h:37
Int_t fLastFixed
Definition: TFumili.h:31
Double_t * GetZ() const
Definition: TFumili.h:110
Double_t * fAMX
Definition: TFumili.h:58
#define ClassImp(name)
Definition: Rtypes.h:336
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*-*-*-*-*-* Basically, it forms the likelihood...
Definition: TFumili.cxx:2005
double f(double x)
Int_t fENDFLG
Definition: TFumili.h:32
double Double_t
Definition: RtypesCore.h:55
Double_t * fSumLog
Definition: TFumili.h:48
Int_t fNstepDec
Definition: TFumili.h:28
virtual ~TFumili()
TFumili destructor.
Definition: TFumili.cxx:212
void Derivatives(Double_t *, Double_t *)
Calculates partial derivatives of theoretical function.
Definition: TFumili.cxx:293
Double_t y[n]
Definition: legend1.C:17
Double_t ey[n]
Definition: legend1.C:17
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:572
The TH1 histogram class.
Definition: TH1.h:56
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
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...
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 Default method (function evaluated at center o...
Definition: TFumili.cxx:1770
static const double eu
Definition: Vavilov.cxx:44
virtual void InitArgs(const Double_t *x, const Double_t *params)
Definition: TF1.cxx:2261
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition: TF1.h:542
Double_t * fAMN
Definition: TFumili.h:59
virtual Int_t GetNpar() const
Definition: TF1.h:435
Double_t * GetY() const
Definition: TGraph.h:130
virtual void ReleaseParameter(Int_t ipar)
Releases parameter number ipar.
Definition: TFumili.cxx:1613
static Bool_t RejectedPoint()
See TF1::RejectPoint above.
Definition: TF1.cxx:3468
virtual Double_t GetSumLog(Int_t)
return Sum(log(i) i=0,n used by log likelihood fits
Definition: TFumili.cxx:931
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
1-Dim function class
Definition: TF1.h:150
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
TF1 * f1
Definition: legend1.C:11
void SetData(Double_t *, Int_t, Int_t)
Sets pointer to data array provided by user.
Definition: TFumili.cxx:1645
void(* fFCN)(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
Double_t * fZ
Definition: TFumili.h:45
Double_t * fPL0
Definition: TFumili.h:53
virtual void Clear(Option_t *opt="")
Resets all parameter names, values and errors to zero.
Definition: TFumili.cxx:239
Int_t fNfcn
Definition: TFumili.h:23
Double_t fEPS
Definition: TFumili.h:66
float type_of_call hi(const int &, const int &)
void SetParNumber(Int_t ParNum)
Definition: TFumili.h:120
Double_t * fA
Definition: TFumili.h:52
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*-*-*-*-*-* Basically, it forms the likelihood...
Definition: TFumili.cxx:1917
const int nn
Double_t Sqrt(Double_t x)
Definition: TMath.h:591
Definition: TFumili.h:36
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1226
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:2149
virtual TObject * GetUserFunc() const
virtual TObject * GetObjectFit() const
const Bool_t kTRUE
Definition: RtypesCore.h:91
Int_t fINDFLG[5]
Definition: TFumili.h:33
const Int_t n
Definition: legend1.C:16
virtual Double_t GetErrorXhigh(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1416
Double_t fGT
Definition: TFumili.h:69
Int_t SGZ()
Evaluates objective function ( chi-square ), gradients and Z-matrix using data provided by user via T...
Definition: TFumili.cxx:1702
double log(double)
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:859
TFumili * gFumili
Definition: TFumili.cxx:120
Int_t fNmaxIter
Definition: TFumili.h:30
Int_t fNlimMul
Definition: TFumili.h:29
const char * Data() const
Definition: TString.h:347