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