Logo ROOT  
Reference Guide
TMinuit.cxx
Go to the documentation of this file.
1 // @(#)root/minuit:$Id$
2 // Author: Rene Brun, Frederick James 12/08/95
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /**
13 
14  \defgroup MinuitOld TMinuit
15  \ingroup Math
16 
17  The Minuit Minimization package.
18  Direct C++ implementation of the Minuit minimization package.
19  This package was originally written in Fortran by Fred James
20  and part of PACKLIB (patch D506)
21  It has been converted to a C++ class, TMinuit, by R.Brun.
22 */
23 
24 /** \class TMinuit
25 
26 Implementation in C++ of the Minuit package written by Fred James.
27 This is a straightforward conversion of the original Fortran version.
28 
29 The main changes are:
30 
31  - The variables in the various Minuit labelled common blocks
32  have been changed to the TMinuit class data members.
33 
34  - The internal arrays with a maximum dimension depending on the
35  maximum number of parameters are now data members arrays with
36  a dynamic dimension such that one can fit very large problems
37  by simply initialising the TMinuit constructor with the maximum
38  number of parameters.
39 
40  - The include file Minuit.h has been commented as much as possible
41  using existing comments in the code or the printed documentation
42 
43  - The original Minuit subroutines are now member functions.
44 
45  - Constructors and destructor have been added.
46 
47  - Instead of passing the FCN function in the argument
48  list, the addresses of this function is stored as pointer
49  in the data members of the class. This is by far more elegant
50  and flexible in an interactive environment.
51  The member function SetFCN can be used to define this pointer.
52 
53  - The ROOT static function Printf is provided to replace all
54  format statements and to print on currently defined output file.
55  - The functions SetObjectFit(TObject * obj)/GetObjectFit() can be
56  used inside the FCN function to set/get a referenced object
57  instead of using global variables.
58 
59 
60 ## Basic concepts of MINUIT
61 
62 The [MINUIT](https://root.cern.ch/sites/d35c7d8c.web.cern.ch/files/minuit.pdf)
63 package acts on a multiparameter Fortran function to which one
64 must give the generic name <TT>FCN</TT>. In the ROOT implementation,
65 the function <TT>FCN</TT> is defined via the MINUIT SetFCN member function
66 when an Histogram.Fit command is invoked.
67 The value of <TT>FCN</TT> will in general depend on one
68 or more variable parameters.
69 
70 To take a simple example, in case of ROOT histograms (classes TH1C,TH1S,TH1F,TH1D)
71 the Fit function defines the Minuit fitting function as being H1FitChisquare
72 or H1FitLikelihood depending on the options selected.
73 H1FitChisquare
74 calculates the chisquare between the user fitting function (gaussian, polynomial,
75 user defined,etc) and the data for given values of the parameters.
76 It is the task of MINUIT to find those values of the parameters
77 which give the lowest value of chisquare.
78 
79 ### Basic concepts - The transformation for parameters with limits.
80 
81 For variable parameters with limits, MINUIT uses the following
82 transformation:
83 
84 \f[
85 P_{\mathrm{int}} = \arcsin
86  \left( 2\: \frac{P_{\mathrm{ext}}-a}{b-a} - 1 \right)
87 P_{\mathrm{ext}} = a + \frac{b - a}{2} \left( \sin P_{\mathrm{int}} + 1 \right)
88 \f]
89 
90 so that the internal value \f$P_{\mathrm{int}}\f$ can take on any value, while
91 the external value \f$P_{\mathrm{ext}}\f$ can take on values only between the lower
92 limit \f$a\f$ and the upper limit \f$b\f$.
93 Since the transformation is necessarily non-linear, it would transform a
94 nice linear problem into a nasty non-linear one, which is the reason why
95 limits should be avoided if not necessary.
96 In addition, the transformation
97 does require some computer time, so it slows down the computation a little
98 bit, and more importantly, it introduces additional numerical inaccuracy into
99 the problem in addition to what is introduced in the numerical calculation
100 of the FCN value.
101 The effects of non-linearity and numerical roundoff both
102 become more important as the external value gets closer to one of the limits
103 (expressed as the distance to nearest limit divided by distance between limits).
104 The user must therefore be aware of the fact that, for example,
105 if he puts limits of \f$(0,10^{10})\f$ on a parameter, then the values \f$0.0\f$
106 and \f$1.0\f$ will be indistinguishable to the accuracy of most machines.
107 
108 The transformation also affects the parameter error matrix, of course,
109 so Minuit does a transformation of the error matrix (and the
110 ``parabolic'' parameter errors) when there are parameter limits.
111 Users should however realize that the transformation is only a linear
112 approximation, and that it cannot give a meaningful result if one or more
113 parameters is very close to a limit, where
114 \f$\partial P_{\mathrm{ext}} / \partial P_{\mathrm{int}} \approx 0\f$.
115 Therefore, it is recommended that:
116 
117  1. Limits on variable parameters should be used only when needed in order
118 to prevent the parameter from taking on unphysical values.
119 
120  2. When a satisfactory minimum has been found using limits, the limits
121 should then be removed if possible, in order to perform or re-perform the
122 error analysis without limits.
123 
124 
125 ### How to get the right answer from MINUIT.
126 
127 MINUIT offers the user a choice of several minimization algorithms. The
128 MIGRAD algorithm is in general the best minimizer for
129 nearly all functions. It is a variable-metric method with inexact line
130 search, a stable metric updating scheme, and checks for
131 positive-definiteness. Its main weakness is that it depends heavily on
132 knowledge of the first derivatives, and fails miserably if they are very
133 inaccurate.
134 
135 If parameter limits are needed, in spite of the side effects, then the
136 user should be aware of the following techniques to alleviate problems
137 caused by limits:
138 
139 #### Getting the right minimum with limits.
140 
141 If MIGRAD converges normally to a point where no parameter is near one of
142 its limits, then the existence of limits has probably not prevented MINUIT
143 from finding the right minimum. On the other hand, if one or more
144 parameters is near its limit at the minimum, this may be because the true
145 minimum is indeed at a limit, or it may be because the minimizer has
146 become ``blocked'' at a limit. This may normally happen only if the
147 parameter is so close to a limit (internal value at an odd multiple of
148 \f$\pm \frac{\pi}{2}\f$
149 that MINUIT prints a warning to this effect when it prints
150 the parameter values.
151 
152 The minimizer can become blocked at a limit, because at a limit
153 the derivative seen by the minimizer
154 \f$\partial F / \partial P_{\mathrm{int}}\f$
155 is zero no matter what the real derivative
156 \f$\partial F / \partial P_{\mathrm{ext}}\f$ is.
157 
158 \f[
159 \frac{\partial F}{\partial P_{\mathrm{int}}} =
160 \frac{\partial F}{\partial P_{\mathrm{ext}}}
161 \frac{\partial P_{\mathrm{ext}}}{\partial P_{\mathrm{int}}} =
162 \frac{\partial F}{\partial P_{\mathrm{ext}}} = 0
163 \f]
164 
165 #### Getting the right parameter errors with limits.
166 
167 In the best case, where the minimum is far from any limits, MINUIT will
168 correctly transform the error matrix, and the parameter errors it reports
169 should be accurate and very close to those you would have got without
170 limits. In other cases (which should be more common, since otherwise you
171 wouldn't need limits), the very meaning of parameter errors becomes
172 problematic. Mathematically, since the limit is an absolute constraint on
173 the parameter, a parameter at its limit has no error, at least in one
174 direction. The error matrix, which can assign only symmetric errors, then
175 becomes essentially meaningless.
176 
177 ### Interpretation of Parameter Errors:
178 
179 There are two kinds of problems that can arise: the reliability of
180 MINUIT's error estimates, and their statistical interpretation, assuming
181 they are accurate.
182 
183 ### Statistical interpretation:
184 
185 For discussion of basic concepts, such as the meaning of the elements of
186 the error matrix, or setting of exact confidence levels see:
187 
188  1. F.James.
189  Determining the statistical Significance of experimental Results.
190  Technical Report DD/81/02 and CERN Report 81-03, CERN, 1981.
191 
192  2. W.T.Eadie, D.Drijard, F.James, M.Roos, and B.Sadoulet.
193  Statistical Methods in Experimental Physics.
194  North-Holland, 1971.
195 
196 ### Reliability of MINUIT error estimates.
197 
198 MINUIT always carries around its own current estimates of the parameter
199 errors, which it will print out on request, no matter how accurate they
200 are at any given point in the execution. For example, at initialization,
201 these estimates are just the starting step sizes as specified by the user.
202 After a HESSE step, the errors are usually quite accurate,
203 unless there has been a problem. MINUIT, when it prints out error values,
204 also gives some indication of how reliable it thinks they are. For
205 example, those marked <TT>CURRENT GUESS ERROR</TT> are only working values
206 not to be believed, and <TT>APPROXIMATE ERROR</TT> means that they have
207 been calculated but there is reason to believe that they may not be
208 accurate.
209 
210 If no mitigating adjective is given, then at least MINUIT believes the
211 errors are accurate, although there is always a small chance that MINUIT
212 has been fooled. Some visible signs that MINUIT may have been fooled are:
213 
214 
215  1. Warning messages produced during the minimization or error analysis.
216 
217  2. Failure to find new minimum.
218 
219  3. Value of <TT>EDM</TT> too big (estimated Distance to Minimum).
220 
221  4. Correlation coefficients exactly equal to zero, unless some parameters
222  are known to be uncorrelated with the others.
223 
224  5. Correlation coefficients very close to one (greater than 0.99). This
225  indicates both an exceptionally difficult problem, and one which has been
226  badly parameterised so that individual errors are not very meaningful
227  because they are so highly correlated.
228 
229  6. Parameter at limit. This condition, signalled by a MINUIT warning
230  message, may make both the function minimum and parameter errors
231  unreliable. See the discussion above ``Getting the right parameter errors
232  with limits''.
233 
234 
235 The best way to be absolutely sure of the errors, is to use
236 ``independent'' calculations and compare them, or compare the calculated
237 errors with a picture of the function. Theoretically, the covariance
238 matrix for a ``physical'' function must be positive-definite at the
239 minimum, although it may not be so for all points far away from the
240 minimum, even for a well-determined physical problem. Therefore, if MIGRAD
241 reports that it has found a non-positive-definite covariance matrix, this
242 may be a sign of one or more of the following:
243 
244 ##### A non-physical region:
245 
246 On its way to the minimum, MIGRAD may have traversed a region which has
247 unphysical behaviour, which is of course not a serious problem as long as
248 it recovers and leaves such a region.
249 
250 ##### An underdetermined problem:
251 
252 If the matrix is not positive-definite even at the minimum, this may mean
253 that the solution is not well-defined, for example that there are more
254 unknowns than there are data points, or that the parameterisation of the
255 fit contains a linear dependence. If this is the case, then MINUIT (or any
256 other program) cannot solve your problem uniquely, and the error matrix
257 will necessarily be largely meaningless, so the user must remove the
258 under-determinedness by reformulating the parameterisation. MINUIT cannot
259 do this itself.
260 
261 ##### Numerical inaccuracies:
262 
263 It is possible that the apparent lack of positive-definiteness is in fact
264 only due to excessive roundoff errors in numerical calculations in the
265 user function or not enough precision. This is unlikely in general, but
266 becomes more likely if the number of free parameters is very large, or if
267 
268 the parameters are badly scaled (not all of the same order of magnitude),
269 and correlations are also large. In any case, whether the
270 non-positive-definiteness is real or only numerical is largely irrelevant,
271 since in both cases the error matrix will be unreliable and the minimum
272 suspicious.
273 
274 ##### An ill-posed problem:
275 
276 For questions of parameter dependence, see the discussion above on
277 positive-definiteness.
278 
279 Possible other mathematical problems are the following:
280 
281 ##### Excessive numerical roundoff:
282 
283 Be especially careful of exponential and factorial functions which get big
284 very quickly and lose accuracy.
285 
286 ##### Starting too far from the solution:
287 
288 The function may have unphysical local minima, especially at infinity in
289 some variables.
290 
291 ##### Minuit parameter errors in the presence of limits
292 This concerns the way Minuit reports the symmetric (or parabolic) errors
293 on parameters. It does not apply to the errors reported from Minos, which
294 are in general asymmetric.
295 
296 The symmetric errors reported by Minuit are always calculated from
297 the covariance matrix, assuming that this matrix has been calculated,
298 usually as the result of a Migrad minimization or a direct
299 calculation by Hesse which inverts the second derivative matrix.
300 
301 When there are no limits on the parameter in question, the error reported
302 by Minuit should therefore be exactly equal to the square root of the
303 corresponding diagonal element of the error matrix reported by Minuit.
304 
305 However, when there are limits on the parameter, there is a transformation
306 between the internal parameter values seen by Minuit (which are unbounded)
307 and the external parameter values seen by the user in FCN (which remain
308 inside the desired limits). Therefore the internal error matrix kept by
309 Minuit must be transformed to an external error matrix for the user.
310 This is done by multiplying the (I,J)th element by DEXDIN(I)*DEXDIN(J),
311 where DEXDIN is the derivative of the external value with respect to the
312 internal value at the minimum. This is a linearisation of the
313 transformation, and is the only way to produce an error matrix in external
314 coordinates meaningful to the user. But when reporting the individual
315 parabolic errors for limited parameters, Minuit can do a little better, so
316 it does. In this case, Minuit actually transforms the ends of the
317 internal "error bar" to external coordinates and reports the length of
318 this transformed interval. Strictly speaking, it is now asymmetric, but
319 since the origin of the asymmetry is only an artificial transformation it
320 does not make much sense, so the transformed errors are symmetrized.
321 
322 The result of all the above is that for parameters with limits, the error
323 reported by Minuit is not exactly equal to the square root of the diagonal
324 element of the error matrix. The difference is a measure of how much the
325 limits deform the problem. If possible, it is suggested not to use limits
326 on parameters, and the problem goes away. If for some reason limits are
327 necessary, and you are sensitive to the difference between the two ways of
328 calculating the errors, it is suggested to use Minos errors which take
329 into account the non-linearities much more precisely.
330 
331 @ingroup MinuitOld
332 */
333 
334 #include <stdlib.h>
335 #include <stdio.h>
336 
337 #include "TROOT.h"
338 #include "TList.h"
339 #include "TMinuit.h"
340 #include "TMath.h"
341 #include "TError.h"
342 #include "TPluginManager.h"
343 #include "TClass.h"
344 
345 #include <atomic>
346 
348 
349 static const char charal[29] = " .ABCDEFGHIJKLMNOPQRSTUVWXYZ";
350 
352 
353 ////////////////////////////////////////////////////////////////////////////////
354 /// Minuit normal constructor
355 ///
356 
357 TMinuit::TMinuit(): TNamed("MINUIT","The Minimization package")
358 {
359  if (TMinuit::Class()->IsCallingNew() != TClass::kRealNew) {
360  //preset all pointers to null
361  fCpnam = 0;
362  fU = 0;
363  fAlim = 0;
364  fBlim = 0;
365  fPstar = 0;
366  fGin = 0;
367  fNvarl = 0;
368  fNiofex = 0;
369 
370  fNexofi = 0;
371  fIpfix = 0;
372  fErp = 0;
373  fErn = 0;
374  fWerr = 0;
375  fGlobcc = 0;
376  fX = 0;
377  fXt = 0;
378  fDirin = 0;
379  fXs = 0;
380  fXts = 0;
381  fDirins = 0;
382  fGrd = 0;
383  fG2 = 0;
384  fGstep = 0;
385  fDgrd = 0;
386  fGrds = 0;
387  fG2s = 0;
388  fGsteps = 0;
389  fPstst = 0;
390  fPbar = 0;
391  fPrho = 0;
392  fWord7 = 0;
393  fVhmat = 0;
394  fVthmat = 0;
395  fP = 0;
396  fXpt = 0;
397  fYpt = 0;
398  fChpt = 0;
399  fCONTgcc = 0;
400  fCONTw = 0;
401  fFIXPyy = 0;
402  fGRADgf = 0;
403  fHESSyy = 0;
404  fIMPRdsav = 0;
405  fIMPRy = 0;
406  fMATUvline = 0;
407  fMIGRflnu = 0;
408  fMIGRstep = 0;
409  fMIGRgs = 0;
410  fMIGRvg = 0;
411  fMIGRxxs = 0;
412  fMNOTxdev = 0;
413  fMNOTw = 0;
414  fMNOTgcc = 0;
415  fPSDFs = 0;
416  fSEEKxmid = 0;
417  fSEEKxbest = 0;
418  fSIMPy = 0;
419  fVERTq = 0;
420  fVERTs = 0;
421  fVERTpp = 0;
422  fCOMDplist = 0;
423  fPARSplist = 0;
424 
425  fUp = 0;
426  fEpsi = 0;
427  fApsi = 0;
428  fXmidcr = 0;
429  fYmidcr = 0;
430  fXdircr = 0;
431  fYdircr = 0;
432 
433  fStatus = 0;
434  fEmpty = 0;
435  fObjectFit = 0;
436  fMethodCall = 0;
437  fPlot = 0;
439 
440  } else {
441  BuildArrays(25);
442 
443  fUp = 0;
444  fEpsi = 0;
445  fApsi = 0;
446  fXmidcr = 0;
447  fYmidcr = 0;
448  fXdircr = 0;
449  fYdircr = 0;
450 
451  fStatus = 0;
452  fEmpty = 0;
453  fObjectFit = 0;
454  fMethodCall = 0;
455  fPlot = 0;
458  mninit(5,6,7);
459  }
460 
461  fFCN = 0;
462  {
464  gROOT->GetListOfSpecials()->Add(this);
465  }
466  gMinuit = this;
467 }
468 
469 ////////////////////////////////////////////////////////////////////////////////
470 /// Minuit normal constructor
471 ///
472 /// maxpar is the maximum number of parameters used with this TMinuit object.
473 
474 TMinuit::TMinuit(Int_t maxpar): TNamed("MINUIT","The Minimization package")
475 {
476  fFCN = 0;
477 
478  BuildArrays(maxpar);
479 
480  fStatus = 0;
481  fEmpty = 0;
482  fObjectFit = 0;
483  fMethodCall = 0;
484  fPlot = 0;
487 
488  mninit(5,6,7);
489  {
491  gROOT->GetListOfSpecials()->Add(this);
492  }
493  gMinuit = this;
494 }
495 
496 ////////////////////////////////////////////////////////////////////////////////
497 /// Private TMinuit copy ctor. TMinuit can not be copied.
498 
499 TMinuit::TMinuit(const TMinuit &minuit) : TNamed(minuit)
500 {
501  Error("TMinuit", "can not copy construct TMinuit");
502 }
503 
504 ////////////////////////////////////////////////////////////////////////////////
505 /// Minuit default destructor
506 
508 {
509  DeleteArrays();
510  delete fPlot;
511  delete fMethodCall;
512  {
514  if (gROOT != 0 && gROOT->GetListOfSpecials() != 0) gROOT->GetListOfSpecials()->Remove(this);
515  }
516  if (gMinuit == this) gMinuit = nullptr;
517 }
518 
519 ////////////////////////////////////////////////////////////////////////////////
520 /// Create internal Minuit arrays for the maxpar parameters
521 
523 {
524  fMaxpar = 25;
525  if (maxpar >= fMaxpar) fMaxpar = maxpar+1;
526  fMaxpar1= fMaxpar*(fMaxpar+1);
527  fMaxpar2= 2*fMaxpar;
528  fMaxpar5= fMaxpar1/2;
529  fMaxcpt = 101;
530  fCpnam = new TString[fMaxpar2];
531  fU = new Double_t[fMaxpar2];
532  fAlim = new Double_t[fMaxpar2];
533  fBlim = new Double_t[fMaxpar2];
534  fPstar = new Double_t[fMaxpar2];
535  fGin = new Double_t[fMaxpar2];
536  fNvarl = new Int_t[fMaxpar2];
537  fNiofex = new Int_t[fMaxpar2];
538 
539  fNexofi = new Int_t[fMaxpar];
540  fIpfix = new Int_t[fMaxpar];
541  fErp = new Double_t[fMaxpar];
542  fErn = new Double_t[fMaxpar];
543  fWerr = new Double_t[fMaxpar];
544  fGlobcc = new Double_t[fMaxpar];
545  fX = new Double_t[fMaxpar];
546  fXt = new Double_t[fMaxpar];
547  fDirin = new Double_t[fMaxpar];
548  fXs = new Double_t[fMaxpar];
549  fXts = new Double_t[fMaxpar];
550  fDirins = new Double_t[fMaxpar];
551  fGrd = new Double_t[fMaxpar];
552  fG2 = new Double_t[fMaxpar];
553  fGstep = new Double_t[fMaxpar];
554  fDgrd = new Double_t[fMaxpar];
555  fGrds = new Double_t[fMaxpar];
556  fG2s = new Double_t[fMaxpar];
557  fGsteps = new Double_t[fMaxpar];
558  fPstst = new Double_t[fMaxpar];
559  fPbar = new Double_t[fMaxpar];
560  fPrho = new Double_t[fMaxpar];
561  fWord7 = new Double_t[fMaxpar];
562  fVhmat = new Double_t[fMaxpar5];
563  fVthmat = new Double_t[fMaxpar5];
564  fP = new Double_t[fMaxpar1];
565  fXpt = new Double_t[fMaxcpt];
566  fYpt = new Double_t[fMaxcpt];
567  fChpt = new char[fMaxcpt+1];
568  // initialisation of dynamic arrays used internally in some functions
569  // these arrays had a fix dimension in Minuit
570  fCONTgcc = new Double_t[fMaxpar];
571  fCONTw = new Double_t[fMaxpar];
572  fFIXPyy = new Double_t[fMaxpar];
573  fGRADgf = new Double_t[fMaxpar];
574  fHESSyy = new Double_t[fMaxpar];
575  fIMPRdsav = new Double_t[fMaxpar];
576  fIMPRy = new Double_t[fMaxpar];
577  fMATUvline = new Double_t[fMaxpar];
578  fMIGRflnu = new Double_t[fMaxpar];
579  fMIGRstep = new Double_t[fMaxpar];
580  fMIGRgs = new Double_t[fMaxpar];
581  fMIGRvg = new Double_t[fMaxpar];
582  fMIGRxxs = new Double_t[fMaxpar];
583  fMNOTxdev = new Double_t[fMaxpar];
584  fMNOTw = new Double_t[fMaxpar];
585  fMNOTgcc = new Double_t[fMaxpar];
586  fPSDFs = new Double_t[fMaxpar];
587  fSEEKxmid = new Double_t[fMaxpar];
588  fSEEKxbest = new Double_t[fMaxpar];
589  fSIMPy = new Double_t[fMaxpar];
590  fVERTq = new Double_t[fMaxpar];
591  fVERTs = new Double_t[fMaxpar];
592  fVERTpp = new Double_t[fMaxpar];
593  fCOMDplist = new Double_t[fMaxpar];
594  fPARSplist = new Double_t[fMaxpar];
595 
596  for (int i = 0; i < fMaxpar; i++) {
597  fErp[i] = 0;
598  fErn[i] = 0;
599  }
600 }
601 
602 ////////////////////////////////////////////////////////////////////////////////
603 /// Make a clone of an object using the Streamer facility.
604 /// Function pointer is copied to Clone
605 
606 TObject *TMinuit::Clone(const char *newname) const
607 {
608  TMinuit *named = (TMinuit*)TNamed::Clone(newname);
609  named->fFCN=fFCN;
610  return named;
611 }
612 
613 ////////////////////////////////////////////////////////////////////////////////
614 /// Execute a Minuit command
615 ///
616 /// Equivalent to MNEXCM except that the command is given as a character string.
617 /// See TMinuit::mnhelp for the full list of available commands
618 /// See also the
619 /// [complete documentation of all the available commands](https://root.cern.ch/sites/d35c7d8c.web.cern.ch/files/minuit.pdf)
620 ///
621 /// Returns the status of the execution:
622 /// - 0: command executed normally
623 /// - 1: command is blank, ignored
624 /// - 2: command line unreadable, ignored
625 /// - 3: unknown command, ignored
626 /// - 4: abnormal termination (e.g., MIGRAD not converged)
627 /// - 5: command is a request to read PARAMETER definitions
628 /// - 6: 'SET INPUT' command
629 /// - 7: 'SET TITLE' command
630 /// - 8: 'SET COVAR' command
631 /// - 9: reserved
632 /// - 10: END command
633 /// - 11: EXIT or STOP command
634 /// - 12: RETURN command
635 
636 Int_t TMinuit::Command(const char *command)
637 {
638  Int_t status = 0;
639  mncomd(command,status);
640  return status;
641 }
642 
643 ////////////////////////////////////////////////////////////////////////////////
644 /// Creates a TGraph object describing the n-sigma contour of a
645 /// TMinuit fit. The contour of the parameters pa1 and pa2 is calculated
646 /// using npoints (>=4) points. The TMinuit status will be
647 /// - 0 on success and
648 /// - -1 if errors in the calling sequence (pa1, pa2 not variable)
649 /// - 1 if less than four points can be found
650 /// - 2 if npoints<4
651 /// - n>3 if only n points can be found (n < npoints)
652 /// The status can be obtained via TMinuit::GetStatus().
653 ///
654 /// To get the n-sigma contour the ERRDEF parameter in Minuit has to set
655 /// to n^2. The fcn function has to be set before the routine is called.
656 ///
657 /// The TGraph object is created via the interpreter. The user must cast it
658 /// to a TGraph*. Note that the TGraph is created with npoints+1 in order to
659 /// close the contour (setting last point equal to first point).
660 ///
661 /// You can find an example in $ROOTSYS/tutorials/fit/fitcont.C
662 
664 {
665  if (npoints<4) {
666  // we need at least 4 points
667  fStatus= 2;
668  return (TObject *)0;
669  }
670  Int_t npfound;
671  Double_t *xcoor = new Double_t[npoints+1];
672  Double_t *ycoor = new Double_t[npoints+1];
673  mncont(pa1,pa2,npoints,xcoor,ycoor,npfound);
674  if (npfound<4) {
675  // mncont did go wrong
676  Warning("Contour","Cannot find more than 4 points, no TGraph returned");
677  fStatus= (npfound==0 ? 1 : npfound);
678  delete [] xcoor;
679  delete [] ycoor;
680  return (TObject *)0;
681  }
682  if (npfound!=npoints) {
683  // mncont did go wrong
684  Warning("Contour","Returning a TGraph with %d points only",npfound);
685  npoints = npfound;
686  }
687  fStatus=0;
688  // create graph via the PluginManager
689  xcoor[npoints] = xcoor[0]; // add first point at end to get closed polyline
690  ycoor[npoints] = ycoor[0];
691  TObject *gr = 0;
692  TPluginHandler *h;
693  if ((h = gROOT->GetPluginManager()->FindHandler("TMinuitGraph"))) {
694  if (h->LoadPlugin() != -1)
695  gr = (TObject*)h->ExecPlugin(3,npoints+1,xcoor,ycoor);
696  }
697  delete [] xcoor;
698  delete [] ycoor;
699  return gr;
700 }
701 
702 ////////////////////////////////////////////////////////////////////////////////
703 /// Define a parameter
704 
705 Int_t TMinuit::DefineParameter( Int_t parNo, const char *name, Double_t initVal, Double_t initErr, Double_t lowerLimit, Double_t upperLimit )
706 {
707  Int_t err;
708 
709  TString sname = name;
710  mnparm( parNo, sname, initVal, initErr, lowerLimit, upperLimit, err);
711 
712  return err;
713 }
714 
715 ////////////////////////////////////////////////////////////////////////////////
716 /// Delete internal Minuit arrays
717 
719 {
720  if (fEmpty) return;
721  delete [] fCpnam;
722  delete [] fU;
723  delete [] fAlim;
724  delete [] fBlim;
725  delete [] fErp;
726  delete [] fErn;
727  delete [] fWerr;
728  delete [] fGlobcc;
729  delete [] fNvarl;
730  delete [] fNiofex;
731  delete [] fNexofi;
732  delete [] fX;
733  delete [] fXt;
734  delete [] fDirin;
735  delete [] fXs;
736  delete [] fXts;
737  delete [] fDirins;
738  delete [] fGrd;
739  delete [] fG2;
740  delete [] fGstep;
741  delete [] fGin;
742  delete [] fDgrd;
743  delete [] fGrds;
744  delete [] fG2s;
745  delete [] fGsteps;
746  delete [] fIpfix;
747  delete [] fVhmat;
748  delete [] fVthmat;
749  delete [] fP;
750  delete [] fPstar;
751  delete [] fPstst;
752  delete [] fPbar;
753  delete [] fPrho;
754  delete [] fWord7;
755  delete [] fXpt;
756  delete [] fYpt;
757  delete [] fChpt;
758 
759  delete [] fCONTgcc;
760  delete [] fCONTw;
761  delete [] fFIXPyy;
762  delete [] fGRADgf;
763  delete [] fHESSyy;
764  delete [] fIMPRdsav;
765  delete [] fIMPRy;
766  delete [] fMATUvline;
767  delete [] fMIGRflnu;
768  delete [] fMIGRstep;
769  delete [] fMIGRgs;
770  delete [] fMIGRvg;
771  delete [] fMIGRxxs;
772  delete [] fMNOTxdev;
773  delete [] fMNOTw;
774  delete [] fMNOTgcc;
775  delete [] fPSDFs;
776  delete [] fSEEKxmid;
777  delete [] fSEEKxbest;
778  delete [] fSIMPy;
779  delete [] fVERTq;
780  delete [] fVERTs;
781  delete [] fVERTpp;
782  delete [] fCOMDplist;
783  delete [] fPARSplist;
784 
785  fEmpty = 1;
786 }
787 
788 ////////////////////////////////////////////////////////////////////////////////
789 /// Evaluate the minimisation function
790 /// Input parameters:
791 /// - npar: number of currently variable parameters
792 /// - par: array of (constant and variable) parameters
793 /// - flag: Indicates what is to be calculated (see example below)
794 /// - grad: array of gradients
795 /// Output parameters:
796 /// - fval: The calculated function value.
797 /// - grad: The (optional) vector of first derivatives).
798 ///
799 /// The meaning of the parameters par is of course defined by the user,
800 /// who uses the values of those parameters to calculate their function value.
801 /// The starting values must be specified by the user.
802 /// Later values are determined by Minuit as it searches for the minimum
803 /// or performs whatever analysis is requested by the user.
804 ///
805 /// Note that this virtual function may be redefined in a class derived from TMinuit.
806 /// The default function calls the function specified in SetFCN
807 ///
808 /// Example of Minimisation function:
809 
810 Int_t TMinuit::Eval(Int_t npar, Double_t *grad, Double_t &fval, Double_t *par, Int_t flag)
811 {
812 /*
813  if (flag == 1) {
814  read input data,
815  calculate any necessary constants, etc.
816  }
817  if (flag == 2) {
818  calculate GRAD, the first derivatives of FVAL
819  (this is optional)
820  }
821  Always calculate the value of the function, FVAL,
822  which is usually a chisquare or log likelihood.
823  if (iflag == 3) {
824  will come here only after the fit is finished.
825  Perform any final calculations, output fitted data, etc.
826  }
827 */
828 // See concrete examples in TH1::H1FitChisquare, H1FitLikelihood
829 
830  if (fFCN) (*fFCN)(npar,grad,fval,par,flag);
831  return 0;
832 }
833 
834 ////////////////////////////////////////////////////////////////////////////////
835 /// fix a parameter
836 
838 {
839  Int_t err;
840  Double_t tmp[1];
841  tmp[0] = parNo+1; //set internal Minuit numbering
842 
843  mnexcm( "FIX", tmp, 1, err );
844 
845  return err;
846 }
847 
848 ////////////////////////////////////////////////////////////////////////////////
849 /// return parameter value and error
850 
851 Int_t TMinuit::GetParameter( Int_t parNo, Double_t &currentValue, Double_t &currentError ) const
852 {
853  Int_t err;
854  TString name; // ignored
855  Double_t bnd1, bnd2; // ignored
856 
857  mnpout( parNo, name, currentValue, currentError, bnd1, bnd2, err );
858 
859  return err;
860 }
861 
862 ////////////////////////////////////////////////////////////////////////////////
863 /// returns the number of currently fixed parameters
864 
866 {
867  return fNpfix;
868 }
869 
870 ////////////////////////////////////////////////////////////////////////////////
871 /// returns the number of currently free parameters
872 
874 {
875  return fNpar;
876 }
877 
878 ////////////////////////////////////////////////////////////////////////////////
879 /// returns the total number of parameters that have been defined
880 /// as fixed or free. The constant parameters are not counted.
881 
883 {
884  return fNpar + fNpfix;
885 }
886 
887 ////////////////////////////////////////////////////////////////////////////////
888 /// invokes the MIGRAD minimizer
889 
891 {
892  Int_t err;
893  Double_t tmp[1];
894  tmp[0] = 0;
895 
896  mnexcm( "MIGRAD", tmp, 0, err );
897 
898  return err;
899 }
900 
901 ////////////////////////////////////////////////////////////////////////////////
902 /// release a parameter
903 
905 {
906  Int_t err;
907  Double_t tmp[1];
908  tmp[0] = parNo+1; //set internal Minuit numbering
909 
910  mnexcm( "RELEASE", tmp, 1, err );
911 
912  return err;
913 }
914 
915 ////////////////////////////////////////////////////////////////////////////////
916 /// To get the n-sigma contour the error def parameter "up" has to set to n^2.
917 
919 {
920  Int_t err;
921 
922  mnexcm( "SET ERRDEF", &up, 1, err );
923 
924  return err;
925 }
926 
927 ////////////////////////////////////////////////////////////////////////////////
928 /// To set the address of the minimization function
929 
930 void TMinuit::SetFCN(void (*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
931 {
932  fFCN = fcn;
933 }
934 
935 ////////////////////////////////////////////////////////////////////////////////
936 /// Static function called when SetFCN is called in interactive mode
937 
938 void InteractiveFCNm(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
939 {
941  if (!m) return;
942 
943  Long_t args[5];
944  args[0] = (Long_t)&npar;
945  args[1] = (Long_t)gin;
946  args[2] = (Long_t)&f;
947  args[3] = (Long_t)u;
948  args[4] = (Long_t)flag;
949  m->SetParamPtrs(args);
950  Double_t result;
951  m->Execute(result);
952 }
953 
954 ////////////////////////////////////////////////////////////////////////////////
955 /// set Minuit print level.
956 ///
957 /// printlevel:
958 /// - = -1 quiet (also suppress all warnings)
959 /// - = 0 normal
960 /// - = 1 verbose
961 
963 {
964  Int_t err;
965  Double_t tmp[1];
966  tmp[0] = printLevel;
967 
968  mnexcm( "SET PRINT", tmp, 1, err );
969 
970  if (printLevel <=-1) mnexcm("SET NOWarnings",tmp,0,err);
971 
972  return err;
973 }
974 
975 ////////////////////////////////////////////////////////////////////////////////
976 /// Initialize AMIN
977 ///
978 /// Called from many places. Initializes the value of AMIN by
979 /// calling the user function. Prints out the function value and
980 /// parameter values if Print Flag value is high enough.
981 
983 {
984  /* Local variables */
985  Double_t fnew;
986  Int_t nparx;
987 
988  nparx = fNpar;
989  if (fISW[4] >= 1) {
990  Printf(" FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.");
991  }
992  mnexin(fX);
993  Eval(nparx, fGin, fnew, fU, 4); ++fNfcn;
994  fAmin = fnew;
995  fEDM = fBigedm;
996 }
997 
998 ////////////////////////////////////////////////////////////////////////////////
999 /// Compute reasonable histogram intervals
1000 ///
1001 /// Function TO DETERMINE REASONABLE HISTOGRAM INTERVALS
1002 /// GIVEN ABSOLUTE UPPER AND LOWER BOUNDS A1 AND A2
1003 /// AND DESIRED MAXIMUM NUMBER OF BINS NAA
1004 /// PROGRAM MAKES REASONABLE BINNING FROM BL TO BH OF WIDTH BWID
1005 /// F. JAMES, AUGUST, 1974 , stolen for Minuit, 1988
1006 
1007 void TMinuit::mnbins(Double_t a1, Double_t a2, Int_t naa, Double_t &bl, Double_t &bh, Int_t &nb, Double_t &bwid)
1008 {
1009  /* Local variables */
1010  Double_t awid,ah, al, sigfig, sigrnd, alb;
1011  Int_t kwid, lwid, na=0, log_;
1012 
1013  al = TMath::Min(a1,a2);
1014  ah = TMath::Max(a1,a2);
1015  if (al == ah) ah = al + 1;
1016 
1017 // IF NAA .EQ. -1 , PROGRAM USES BWID INPUT FROM CALLING ROUTINE
1018  if (naa == -1) goto L150;
1019 L10:
1020  na = naa - 1;
1021  if (na < 1) na = 1;
1022 
1023 // GET NOMINAL BIN WIDTH IN EXPON FORM
1024 L20:
1025  awid = (ah-al) / Double_t(na);
1026  log_ = Int_t(TMath::Log10(awid));
1027  if (awid <= 1) --log_;
1028  sigfig = awid*TMath::Power(10, -log_);
1029 // ROUND MANTISSA UP TO 2, 2.5, 5, OR 10
1030  if (sigfig > 2) goto L40;
1031  sigrnd = 2;
1032  goto L100;
1033 L40:
1034  if (sigfig > 2.5) goto L50;
1035  sigrnd = 2.5;
1036  goto L100;
1037 L50:
1038  if (sigfig > 5) goto L60;
1039  sigrnd = 5;
1040  goto L100;
1041 L60:
1042  sigrnd = 1;
1043  ++log_;
1044 L100:
1045  bwid = sigrnd*TMath::Power(10, log_);
1046  goto L200;
1047 // GET NEW BOUNDS FROM NEW WIDTH BWID
1048 L150:
1049  if (bwid <= 0) goto L10;
1050 L200:
1051  alb = al / bwid;
1052  lwid = Int_t(alb);
1053  if (alb < 0) --lwid;
1054  bl = bwid*Double_t(lwid);
1055  alb = ah / bwid + 1;
1056  kwid = Int_t(alb);
1057  if (alb < 0) --kwid;
1058  bh = bwid*Double_t(kwid);
1059  nb = kwid - lwid;
1060  if (naa > 5) goto L240;
1061  if (naa == -1) return;
1062 // REQUEST FOR ONE BIN IS DIFFICULT CASE
1063  if (naa > 1 || nb == 1) return;
1064  bwid *= 2;
1065  nb = 1;
1066  return;
1067 L240:
1068  if (nb << 1 != naa) return;
1069  ++na;
1070  goto L20;
1071 }
1072 
1073 ////////////////////////////////////////////////////////////////////////////////
1074 /// Transform FCN to find further minima
1075 ///
1076 /// Called only from MNIMPR. Transforms the function FCN
1077 /// by dividing out the quadratic part in order to find further
1078 /// minima. Calculates `ycalf = (f-fmin)/(x-xmin)*v*(x-xmin)`
1079 
1080 void TMinuit::mncalf(Double_t *pvec, Double_t &ycalf)
1081 {
1082  /* Local variables */
1083  Int_t ndex, i, j, m, n, nparx;
1084  Double_t denom, f;
1085 
1086  nparx = fNpar;
1087  mninex(&pvec[0]);
1088  Eval(nparx, fGin, f, fU, 4); ++fNfcn;
1089  for (i = 1; i <= fNpar; ++i) {
1090  fGrd[i-1] = 0;
1091  for (j = 1; j <= fNpar; ++j) {
1092  m = TMath::Max(i,j);
1093  n = TMath::Min(i,j);
1094  ndex = m*(m-1) / 2 + n;
1095  fGrd[i-1] += fVthmat[ndex-1]*(fXt[j-1] - pvec[j-1]);
1096  }
1097  }
1098  denom = 0;
1099  for (i = 1; i <= fNpar; ++i) {denom += fGrd[i-1]*(fXt[i-1] - pvec[i-1]); }
1100  if (denom <= 0) {
1101  fDcovar = 1;
1102  fISW[1] = 0;
1103  denom = 1;
1104  }
1105  ycalf = (f - fApsi) / denom;
1106 }
1107 
1108 ////////////////////////////////////////////////////////////////////////////////
1109 /// Resets the parameter list to UNDEFINED
1110 ///
1111 /// Called from MINUIT and by option from MNEXCM
1112 
1114 {
1115  Int_t i;
1116 
1117  fNpfix = 0;
1118  fNu = 0;
1119  fNpar = 0;
1120  fNfcn = 0;
1121  fNwrmes[0] = 0;
1122  fNwrmes[1] = 0;
1123  for (i = 1; i <= fMaxext; ++i) {
1124  fU[i-1] = 0;
1125  fCpnam[i-1] = fCundef;
1126  fNvarl[i-1] = -1;
1127  fNiofex[i-1] = 0;
1128  }
1129  mnrset(1);
1130  fCfrom = "CLEAR ";
1131  fNfcnfr = fNfcn;
1132  fCstatu = "UNDEFINED ";
1133  fLnolim = kTRUE;
1134  fLphead = kTRUE;
1135 }
1136 
1137 ////////////////////////////////////////////////////////////////////////////////
1138 /// Print function contours in two variables, on line printer
1139 ///
1140 /// input arguments: parx, pary, devs, ngrid
1141 
1142 void TMinuit::mncntr(Int_t ike1, Int_t ike2, Int_t &ierrf)
1143 {
1144  static const char *const clabel = "0123456789ABCDEFGHIJ";
1145 
1146  /* Local variables */
1147  Double_t d__1, d__2;
1148  Double_t fcna[115], fcnb[115], contur[20];
1149  Double_t ylabel, fmn, fmx, xlo, ylo, xup, yup;
1150  Double_t devs, xsav, ysav, bwidx, bwidy, unext, ff, xb4;
1151  Int_t i, ngrid, ixmid, nparx, ix, nx, ny, ki1, ki2, ixzero, iy, ics;
1152  TString chmid, chln, chzero;
1153 
1154  Int_t ke1 = ike1+1;
1155  Int_t ke2 = ike2+1;
1156  if (ke1 <= 0 || ke2 <= 0) goto L1350;
1157  if (ke1 > fNu || ke2 > fNu) goto L1350;
1158  ki1 = fNiofex[ke1-1];
1159  ki2 = fNiofex[ke2-1];
1160  if (ki1 <= 0 || ki2 <= 0) goto L1350;
1161  if (ki1 == ki2) goto L1350;
1162 
1163  if (fISW[1] < 1) {
1164  mnhess();
1165  mnwerr();
1166  }
1167  nparx = fNpar;
1168  xsav = fU[ke1-1];
1169  ysav = fU[ke2-1];
1170  devs = fWord7[2];
1171  if (devs <= 0) devs = 2;
1172  xlo = fU[ke1-1] - devs*fWerr[ki1-1];
1173  xup = fU[ke1-1] + devs*fWerr[ki1-1];
1174  ylo = fU[ke2-1] - devs*fWerr[ki2-1];
1175  yup = fU[ke2-1] + devs*fWerr[ki2-1];
1176  ngrid = Int_t(fWord7[3]);
1177  if (ngrid <= 0) {
1178  ngrid = 25;
1179 // Computing MIN
1180  nx = TMath::Min(fNpagwd - 15,ngrid);
1181 // Computing MIN
1182  ny = TMath::Min(fNpagln - 7,ngrid);
1183  } else {
1184  nx = ngrid;
1185  ny = ngrid;
1186  }
1187  if (nx < 11) nx = 11;
1188  if (ny < 11) ny = 11;
1189  if (nx >= 115) nx = 114;
1190 
1191 // ask if parameter outside limits
1192  if (fNvarl[ke1-1] > 1) {
1193  if (xlo < fAlim[ke1-1]) xlo = fAlim[ke1-1];
1194  if (xup > fBlim[ke1-1]) xup = fBlim[ke1-1];
1195  }
1196  if (fNvarl[ke2-1] > 1) {
1197  if (ylo < fAlim[ke2-1]) ylo = fAlim[ke2-1];
1198  if (yup > fBlim[ke2-1]) yup = fBlim[ke2-1];
1199  }
1200  bwidx = (xup - xlo) / Double_t(nx);
1201  bwidy = (yup - ylo) / Double_t(ny);
1202  ixmid = Int_t(((xsav - xlo)*Double_t(nx) / (xup - xlo)) + 1);
1203  if (ixmid < 1) ixmid = 1;
1204  if (fAmin == fUndefi) mnamin();
1205 
1206  for (i = 1; i <= 20; ++i) { contur[i-1] = fAmin + fUp*(i-1)*(i-1); }
1207  contur[0] += fUp*.01;
1208 // fill FCNB to prepare first row, and find column zero/
1209  fU[ke2-1] = yup;
1210  ixzero = 0;
1211  xb4 = 1;
1212 //TH
1213  chmid.Resize(nx+1);
1214  chzero.Resize(nx+1);
1215  chln.Resize(nx+1);
1216  for (ix = 1; ix <= nx + 1; ++ix) {
1217  fU[ke1-1] = xlo + Double_t(ix-1)*bwidx;
1218  Eval(nparx, fGin, ff, fU, 4);
1219  fcnb[ix-1] = ff;
1220  if (xb4 < 0 && fU[ke1-1] > 0) ixzero = ix - 1;
1221  xb4 = fU[ke1-1];
1222  chmid[ix-1] = '*';
1223  chzero[ix-1] = '-';
1224  }
1225  Printf(" Y-AXIS: PARAMETER %3d: %s",ke2,(const char*)fCpnam[ke2-1]);
1226  if (ixzero > 0) {
1227  chzero[ixzero-1] = '+';
1228  chln = " ";
1229  Printf(" X=0");
1230  }
1231 // loop over rows
1232  for (iy = 1; iy <= ny; ++iy) {
1233  unext = fU[ke2-1] - bwidy;
1234 // prepare this line background pattern for contour
1235  chln = " ";
1236 // TH
1237  chln.Resize(nx+1);
1238  chln[ixmid-1] = '*';
1239  if (ixzero != 0) chln[ixzero-1] = ':';
1240  if (fU[ke2-1] > ysav && unext < ysav) chln = chmid;
1241  if (fU[ke2-1] > 0 && unext < 0) chln = chzero;
1242  fU[ke2-1] = unext;
1243  ylabel = fU[ke2-1] + bwidy*.5;
1244 // move FCNB to FCNA and fill FCNB with next row
1245  for (ix = 1; ix <= nx + 1; ++ix) {
1246  fcna[ix-1] = fcnb[ix-1];
1247  fU[ke1-1] = xlo + Double_t(ix-1)*bwidx;
1248  Eval(nparx, fGin, ff, fU, 4);
1249  fcnb[ix-1] = ff;
1250  }
1251 // look for contours crossing the FCNxy squares
1252  for (ix = 1; ix <= nx; ++ix) {
1253  d__1 = TMath::Max(fcna[ix-1],fcnb[ix-1]),
1254  d__2 = TMath::Max(fcna[ix],fcnb[ix]);
1255  fmx = TMath::Max(d__1,d__2);
1256  d__1 = TMath::Min(fcna[ix-1],fcnb[ix-1]),
1257  d__2 = TMath::Min(fcna[ix],fcnb[ix]);
1258  fmn = TMath::Min(d__1,d__2);
1259  for (ics = 1; ics <= 20; ++ics) {
1260  if (contur[ics-1] > fmn) goto L240;
1261  }
1262  continue;
1263 L240:
1264  if (contur[ics-1] < fmx) chln[ix-1] = clabel[ics-1];
1265  }
1266 // print a row of the contour plot
1267  Printf(" %12.4g %s",ylabel,(const char*)chln);
1268  }
1269 // contours printed, label x-axis
1270  chln = " ";
1271  chln(0,1) = 'I';
1272  chln(ixmid-1,1) = 'I';
1273  chln(nx-1,1) = 'I';
1274  Printf(" %s",(const char*)chln);
1275 
1276 // the hardest of all: print x-axis scale!
1277  chln = " ";
1278  if (nx <= 26) {
1279  Printf(" %12.4g%s%12.4g",xlo,(const char*)chln,xup);
1280  Printf(" %s%12.4g",(const char*)chln,xsav);
1281  } else {
1282  Printf(" %12.4g%s%12.4g%s%12.4g",xlo,(const char*)chln,xsav,(const char*)chln,xup);
1283  }
1284  Printf(" X-AXIS: PARAMETER %3d %s ONE COLUMN=%12.4g"
1285  ,ke1,(const char*)fCpnam[ke1-1],bwidx);
1286  Printf(" FUNCTION VALUES: F(I)=%12.4g +%12.4g *I**2",fAmin,fUp);
1287 // finished. reset input values
1288  fU[ke1-1] = xsav;
1289  fU[ke2-1] = ysav;
1290  ierrf = 0;
1291  return;
1292 L1350:
1293  Printf(" INVALID PARAMETER NUMBER(S) REQUESTED. IGNORED.");
1294  ierrf = 1;
1295 }
1296 
1297 ////////////////////////////////////////////////////////////////////////////////
1298 /// Reads a command string and executes
1299 ///
1300 /// Called by user. 'Reads' a command string and executes.
1301 /// Equivalent to MNEXCM except that the command is given as a
1302 /// character string.
1303 ///
1304 /// ICONDN =
1305 /// - 0: command executed normally
1306 /// - 1: command is blank, ignored
1307 /// - 2: command line unreadable, ignored
1308 /// - 3: unknown command, ignored
1309 /// - 4: abnormal termination (e.g., MIGRAD not converged)
1310 /// - 5: command is a request to read PARAMETER definitions
1311 /// - 6: 'SET INPUT' command
1312 /// - 7: 'SET TITLE' command
1313 /// - 8: 'SET COVAR' command
1314 /// - 9: reserved
1315 /// - 10: END command
1316 /// - 11: EXIT or STOP command
1317 /// - 12: RETURN command
1318 ///
1319 
1320 void TMinuit::mncomd(const char *crdbin, Int_t &icondn)
1321 {
1322  /* Local variables */
1323  Int_t ierr, ipos, i, llist, lenbuf, lnc;
1324  Bool_t leader;
1325  TString comand, crdbuf, ctemp;
1326 
1327  crdbuf = crdbin;
1328  crdbuf.ToUpper();
1329  lenbuf = crdbuf.Length();
1330  icondn = 0;
1331 // record not case-sensitive, get upper case, strip leading blanks
1332  leader = kTRUE;
1333  ipos = 1;
1334  for (i = 1; i <= TMath::Min(20,lenbuf); ++i) {
1335  if (crdbuf[i-1] == '\'') break;
1336  if (crdbuf[i-1] == ' ') {
1337  if (leader) ++ipos;
1338  continue;
1339  }
1340  leader = kFALSE;
1341  }
1342 
1343 // blank or null command
1344  if (ipos > lenbuf) {
1345  Printf(" BLANK COMMAND IGNORED.");
1346  icondn = 1;
1347  return;
1348  }
1349 // preemptive commands
1350 // if command is 'PARAMETER'
1351  if (crdbuf(ipos-1,3) == "PAR") {
1352  icondn = 5;
1353  fLphead = kTRUE;
1354  return;
1355  }
1356 // if command is 'SET INPUT'
1357  if (crdbuf(ipos-1,3) == "SET INP") {
1358  icondn = 6;
1359  fLphead = kTRUE;
1360  return;
1361  }
1362 // if command is 'SET TITLE'
1363  if (crdbuf(ipos-1,7) == "SET TIT") {
1364  icondn = 7;
1365  fLphead = kTRUE;
1366  return;
1367  }
1368 // if command is 'SET COVARIANCE'
1369  if (crdbuf(ipos-1,7) == "SET COV") {
1370  icondn = 8;
1371  fLphead = kTRUE;
1372  return;
1373  }
1374 // crack the command
1375  ctemp = crdbuf(ipos-1,lenbuf-ipos+1);
1376  mncrck(ctemp, 20, comand, lnc, fMaxpar, fCOMDplist, llist, ierr, fIsyswr);
1377  if (ierr > 0) {
1378  Printf(" COMMAND CANNOT BE INTERPRETED");
1379  icondn = 2;
1380  return;
1381  }
1382 
1383  mnexcm(comand.Data(), fCOMDplist, llist, ierr);
1384  icondn = ierr;
1385 }
1386 
1387 ////////////////////////////////////////////////////////////////////////////////
1388 /// Find points along a contour where FCN is minimum
1389 ///
1390 /// Find NPTU points along a contour where the function
1391 ///
1392 /// FMIN (X(KE1),X(KE2)) = AMIN+UP
1393 ///
1394 /// where FMIN is the minimum of FCN with respect to all
1395 /// the other NPAR-2 variable parameters (if any).
1396 ///
1397 /// IERRF on return will be equal to the number of points found:
1398 /// - NPTU if normal termination with NPTU points found
1399 /// - -1 if errors in the calling sequence (KE1, KE2 not variable)
1400 /// - 0 if less than four points can be found (using MNMNOT)
1401 /// - n>3 if only n points can be found (n < NPTU)
1402 ///
1403 /// input arguments: parx, pary, devs, ngrid
1404 
1405 void TMinuit::mncont(Int_t ike1, Int_t ike2, Int_t nptu, Double_t *xptu, Double_t *yptu, Int_t &ierrf)
1406 {
1407  /* System generated locals */
1408  Int_t i__1;
1409 
1410  /* Local variables */
1411  Double_t d__1, d__2;
1412  Double_t dist, xdir, ydir, aopt, u1min, u2min;
1413  Double_t abest, scalx, scaly;
1414  Double_t a1, a2, val2mi, val2pl, dc, sclfac, bigdis, sigsav;
1415  Int_t nall, iold, line, mpar, ierr, inew, move, next, i, j, nfcol, iercr;
1416  Int_t idist=0, npcol, kints, i2, i1, lr, nfcnco=0, ki1, ki2, ki3, ke3;
1417  Int_t nowpts, istrav, nfmxin, isw2, isw4;
1418  Bool_t ldebug;
1419 
1420  /* Function Body */
1421  Int_t ke1 = ike1+1;
1422  Int_t ke2 = ike2+1;
1423  ldebug = fIdbg[6] >= 1;
1424  if (ke1 <= 0 || ke2 <= 0) goto L1350;
1425  if (ke1 > fNu || ke2 > fNu) goto L1350;
1426  ki1 = fNiofex[ke1-1];
1427  ki2 = fNiofex[ke2-1];
1428  if (ki1 <= 0 || ki2 <= 0) goto L1350;
1429  if (ki1 == ki2) goto L1350;
1430  if (nptu < 4) goto L1400;
1431 
1432  nfcnco = fNfcn;
1433  fNfcnmx = (nptu + 5)*100*(fNpar + 1);
1434 // The minimum
1435  mncuve();
1436  u1min = fU[ke1-1];
1437  u2min = fU[ke2-1];
1438  ierrf = 0;
1439  fCfrom = "MNContour ";
1440  fNfcnfr = nfcnco;
1441  if (fISW[4] >= 0) {
1442  Printf(" START MNCONTOUR CALCULATION OF %4d POINTS ON CONTOUR.",nptu);
1443  if (fNpar > 2) {
1444  if (fNpar == 3) {
1445  ki3 = 6 - ki1 - ki2;
1446  ke3 = fNexofi[ki3-1];
1447  Printf(" EACH POINT IS A MINIMUM WITH RESPECT TO PARAMETER %3d %s",ke3,(const char*)fCpnam[ke3-1]);
1448  } else {
1449  Printf(" EACH POINT IS A MINIMUM WITH RESPECT TO THE OTHER %3d VARIABLE PARAMETERS.",fNpar - 2);
1450  }
1451  }
1452  }
1453 
1454 // Find the first four points using MNMNOT
1455 // first two points
1456  mnmnot(ke1, ke2, val2pl, val2mi);
1457  if (fErn[ki1-1] == fUndefi) {
1458  xptu[0] = fAlim[ke1-1];
1459  mnwarn("W", "MNContour ", "Contour squeezed by parameter limits.");
1460  } else {
1461  if (fErn[ki1-1] >= 0) goto L1500;
1462  xptu[0] = u1min + fErn[ki1-1];
1463  }
1464  yptu[0] = val2mi;
1465 
1466  if (fErp[ki1-1] == fUndefi) {
1467  xptu[2] = fBlim[ke1-1];
1468  mnwarn("W", "MNContour ", "Contour squeezed by parameter limits.");
1469  } else {
1470  if (fErp[ki1-1] <= 0) goto L1500;
1471  xptu[2] = u1min + fErp[ki1-1];
1472  }
1473  yptu[2] = val2pl;
1474  scalx = 1 / (xptu[2] - xptu[0]);
1475 // next two points
1476  mnmnot(ke2, ke1, val2pl, val2mi);
1477  if (fErn[ki2-1] == fUndefi) {
1478  yptu[1] = fAlim[ke2-1];
1479  mnwarn("W", "MNContour ", "Contour squeezed by parameter limits.");
1480  } else {
1481  if (fErn[ki2-1] >= 0) goto L1500;
1482  yptu[1] = u2min + fErn[ki2-1];
1483  }
1484  xptu[1] = val2mi;
1485  if (fErp[ki2-1] == fUndefi) {
1486  yptu[3] = fBlim[ke2-1];
1487  mnwarn("W", "MNContour ", "Contour squeezed by parameter limits.");
1488  } else {
1489  if (fErp[ki2-1] <= 0) goto L1500;
1490  yptu[3] = u2min + fErp[ki2-1];
1491  }
1492  xptu[3] = val2pl;
1493  scaly = 1 / (yptu[3] - yptu[1]);
1494  nowpts = 4;
1495  next = 5;
1496  if (ldebug) {
1497  Printf(" Plot of four points found by MINOS");
1498  fXpt[0] = u1min;
1499  fYpt[0] = u2min;
1500  fChpt[0] = ' ';
1501 // Computing MIN
1502  nall = TMath::Min(nowpts + 1,101);
1503  for (i = 2; i <= nall; ++i) {
1504  fXpt[i-1] = xptu[i-2];
1505  fYpt[i-1] = yptu[i-2];
1506  }
1507  sprintf(fChpt,"%s"," ABCD");
1508  mnplot(fXpt, fYpt, fChpt, nall, fNpagwd, fNpagln);
1509  }
1510 
1511 // save some values before fixing
1512  isw2 = fISW[1];
1513  isw4 = fISW[3];
1514  sigsav = fEDM;
1515  istrav = fIstrat;
1516  dc = fDcovar;
1517  fApsi = fEpsi*.5;
1518  abest = fAmin;
1519  mpar = fNpar;
1520  nfmxin = fNfcnmx;
1521  for (i = 1; i <= mpar; ++i) { fXt[i-1] = fX[i-1]; }
1522  i__1 = mpar*(mpar + 1) / 2;
1523  for (j = 1; j <= i__1; ++j) { fVthmat[j-1] = fVhmat[j-1]; }
1524  for (i = 1; i <= mpar; ++i) {
1525  fCONTgcc[i-1] = fGlobcc[i-1];
1526  fCONTw[i-1] = fWerr[i-1];
1527  }
1528 // fix the two parameters in question
1529  kints = fNiofex[ke1-1];
1530  mnfixp(kints-1, ierr);
1531  kints = fNiofex[ke2-1];
1532  mnfixp(kints-1, ierr);
1533 // Fill in the rest of the points
1534  for (inew = next; inew <= nptu; ++inew) {
1535 // find the two neighbouring points with largest separation
1536  bigdis = 0;
1537  for (iold = 1; iold <= inew - 1; ++iold) {
1538  i2 = iold + 1;
1539  if (i2 == inew) i2 = 1;
1540  d__1 = scalx*(xptu[iold-1] - xptu[i2-1]);
1541  d__2 = scaly*(yptu[iold-1] - yptu[i2-1]);
1542  dist = d__1*d__1 + d__2*d__2;
1543  if (dist > bigdis) {
1544  bigdis = dist;
1545  idist = iold;
1546  }
1547  }
1548  i1 = idist;
1549  i2 = i1 + 1;
1550  if (i2 == inew) i2 = 1;
1551 // next point goes between I1 and I2
1552  a1 = .5;
1553  a2 = .5;
1554 L300:
1555  fXmidcr = a1*xptu[i1-1] + a2*xptu[i2-1];
1556  fYmidcr = a1*yptu[i1-1] + a2*yptu[i2-1];
1557  xdir = yptu[i2-1] - yptu[i1-1];
1558  ydir = xptu[i1-1] - xptu[i2-1];
1559  sclfac = TMath::Max(TMath::Abs(xdir*scalx),TMath::Abs(ydir*scaly));
1560  fXdircr = xdir / sclfac;
1561  fYdircr = ydir / sclfac;
1562  fKe1cr = ke1;
1563  fKe2cr = ke2;
1564 // Find the contour crossing point along DIR
1565  fAmin = abest;
1566  mncros(aopt, iercr);
1567  if (iercr > 1) {
1568 // If cannot find mid-point, try closer to point 1
1569  if (a1 > .5) {
1570  if (fISW[4] >= 0) {
1571  Printf(" MNCONT CANNOT FIND NEXT POINT ON CONTOUR. ONLY %3d POINTS FOUND.",nowpts);
1572  }
1573  goto L950;
1574  }
1575  mnwarn("W", "MNContour ", "Cannot find midpoint, try closer.");
1576  a1 = .75;
1577  a2 = .25;
1578  goto L300;
1579  }
1580 // Contour has been located, insert new point in list
1581  for (move = nowpts; move >= i1 + 1; --move) {
1582  xptu[move] = xptu[move-1];
1583  yptu[move] = yptu[move-1];
1584  }
1585  ++nowpts;
1586  xptu[i1] = fXmidcr + fXdircr*aopt;
1587  yptu[i1] = fYmidcr + fYdircr*aopt;
1588  }
1589 L950:
1590 
1591  ierrf = nowpts;
1592  fCstatu = "SUCCESSFUL";
1593  if (nowpts < nptu) fCstatu = "INCOMPLETE";
1594 
1595 // make a lineprinter plot of the contour
1596  if (fISW[4] >= 0) {
1597  fXpt[0] = u1min;
1598  fYpt[0] = u2min;
1599  fChpt[0] = ' ';
1600  nall = TMath::Min(nowpts + 1,101);
1601  for (i = 2; i <= nall; ++i) {
1602  fXpt[i-1] = xptu[i-2];
1603  fYpt[i-1] = yptu[i-2];
1604  fChpt[i-1] = 'X';
1605  }
1606  fChpt[nall] = 0;
1607  Printf(" Y-AXIS: PARAMETER %3d %s",ke2,(const char*)fCpnam[ke2-1]);
1608 
1609  mnplot(fXpt, fYpt, fChpt, nall, fNpagwd, fNpagln);
1610 
1611  Printf(" X-AXIS: PARAMETER %3d %s",ke1,(const char*)fCpnam[ke1-1]);
1612  }
1613 // print out the coordinates around the contour
1614  if (fISW[4] >= 1) {
1615  npcol = (nowpts + 1) / 2;
1616  nfcol = nowpts / 2;
1617  Printf("%5d POINTS ON CONTOUR. FMIN=%13.5e ERRDEF=%11.3g",nowpts,abest,fUp);
1618  Printf(" %s%s%s%s",(const char*)fCpnam[ke1-1],
1619  (const char*)fCpnam[ke2-1],
1620  (const char*)fCpnam[ke1-1],
1621  (const char*)fCpnam[ke2-1]);
1622  for (line = 1; line <= nfcol; ++line) {
1623  lr = line + npcol;
1624  Printf(" %5d%13.5e%13.5e %5d%13.5e%13.5e",line,xptu[line-1],yptu[line-1],lr,xptu[lr-1],yptu[lr-1]);
1625  }
1626  if (nfcol < npcol) {
1627  Printf(" %5d%13.5e%13.5e",npcol,xptu[npcol-1],yptu[npcol-1]);
1628  }
1629  }
1630 // contour finished. reset v
1631  fItaur = 1;
1632  mnfree(1);
1633  mnfree(1);
1634  i__1 = mpar*(mpar + 1) / 2;
1635  for (j = 1; j <= i__1; ++j) { fVhmat[j-1] = fVthmat[j-1]; }
1636  for (i = 1; i <= mpar; ++i) {
1637  fGlobcc[i-1] = fCONTgcc[i-1];
1638  fWerr[i-1] = fCONTw[i-1];
1639  fX[i-1] = fXt[i-1];
1640  }
1641  mninex(fX);
1642  fEDM = sigsav;
1643  fAmin = abest;
1644  fISW[1] = isw2;
1645  fISW[3] = isw4;
1646  fDcovar = dc;
1647  fItaur = 0;
1648  fNfcnmx = nfmxin;
1649  fIstrat = istrav;
1650  fU[ke1-1] = u1min;
1651  fU[ke2-1] = u2min;
1652  goto L2000;
1653 // Error returns
1654 L1350:
1655  Printf(" INVALID PARAMETER NUMBERS.");
1656  goto L1450;
1657 L1400:
1658  Printf(" LESS THAN FOUR POINTS REQUESTED.");
1659 L1450:
1660  ierrf = -1;
1661  fCstatu = "USER ERROR";
1662  goto L2000;
1663 L1500:
1664  Printf(" MNCONT UNABLE TO FIND FOUR POINTS.");
1665  fU[ke1-1] = u1min;
1666  fU[ke2-1] = u2min;
1667  ierrf = 0;
1668  fCstatu = "FAILED";
1669 L2000:
1670  fCfrom = "MNContour ";
1671  fNfcnfr = nfcnco;
1672 }
1673 
1674 ////////////////////////////////////////////////////////////////////////////////
1675 /// Cracks the free-format input
1676 ///
1677 /// Cracks the free-format input, expecting zero or more
1678 /// alphanumeric fields (which it joins into COMAND(1:LNC))
1679 /// followed by one or more numeric fields separated by
1680 /// blanks and/or one comma. The numeric fields are put into
1681 /// the LLIST (but at most MXP) elements of PLIST.
1682 ///
1683 /// IERR :
1684 /// - = 0 if no errors,
1685 /// - = 1 if error(s).
1686 
1687 void TMinuit::mncrck(TString cardbuf, Int_t maxcwd, TString &comand, Int_t &lnc,
1688  Int_t mxp, Double_t *plist, Int_t &llist, Int_t &ierr, Int_t)
1689 {
1690  /* Initialized data */
1691 
1692  char *cnull = 0;
1693  const char *cnumer = "123456789-.0+";
1694 
1695  /* Local variables */
1696  Int_t ifld, iend, lend, left, nreq, ipos, kcmnd, nextb, ic, ibegin, ltoadd;
1697  Int_t ielmnt, lelmnt[25], nelmnt;
1698  TString ctemp;
1699  char *celmnt[25];
1700  char command[25];
1701 
1702  /* Function Body */
1703  char *crdbuf = (char*)cardbuf.Data();
1704  lend = cardbuf.Length();
1705  ielmnt = 0;
1706  nextb = 1;
1707  ierr = 0;
1708 // loop over words CELMNT
1709 L10:
1710  for (ipos = nextb; ipos <= lend; ++ipos) {
1711  ibegin = ipos;
1712  if (crdbuf[ipos-1] == ' ') continue;
1713  if (crdbuf[ipos-1] == ',') goto L250;
1714  goto L150;
1715  }
1716  goto L300;
1717 L150:
1718 // found beginning of word, look for end
1719  for (ipos = ibegin + 1; ipos <= lend; ++ipos) {
1720  if (crdbuf[ipos-1] == ' ') goto L250;
1721  if (crdbuf[ipos-1] == ',') goto L250;
1722  }
1723  ipos = lend + 1;
1724 L250:
1725  iend = ipos - 1;
1726  ++ielmnt;
1727  if (iend >= ibegin) celmnt[ielmnt-1] = &crdbuf[ibegin-1];
1728  else celmnt[ielmnt-1] = cnull;
1729  lelmnt[ielmnt-1] = iend - ibegin + 1;
1730  if (lelmnt[ielmnt-1] > 19) {
1731  Printf(" MINUIT WARNING: INPUT DATA WORD TOO LONG.");
1732  ctemp = cardbuf(ibegin-1,iend-ibegin+1);
1733  Printf(" ORIGINAL:%s",ctemp.Data());
1734  Printf(" TRUNCATED TO:%s",celmnt[ielmnt-1]);
1735  lelmnt[ielmnt-1] = 19;
1736  }
1737  if (ipos >= lend) goto L300;
1738  if (ielmnt >= 25) goto L300;
1739 // look for comma or beginning of next word
1740  for (ipos = iend + 1; ipos <= lend; ++ipos) {
1741  if (crdbuf[ipos-1] == ' ') continue;
1742  nextb = ipos;
1743  if (crdbuf[ipos-1] == ',') nextb = ipos + 1;
1744  goto L10;
1745  }
1746 // All elements found, join the alphabetic ones to
1747 // form a command
1748 L300:
1749  nelmnt = ielmnt;
1750  command[0] = ' '; command[1] = 0;
1751  lnc = 1;
1752  plist[0] = 0;
1753  llist = 0;
1754  if (ielmnt == 0) goto L900;
1755  kcmnd = 0;
1756  for (ielmnt = 1; ielmnt <= nelmnt; ++ielmnt) {
1757  if ( celmnt[ielmnt-1] == cnull) goto L450;
1758  for (ic = 1; ic <= 13; ++ic) {
1759  if (*celmnt[ielmnt-1] == cnumer[ic-1]) goto L450;
1760  }
1761  if (kcmnd >= maxcwd) continue;
1762  left = maxcwd - kcmnd;
1763  ltoadd = lelmnt[ielmnt-1];
1764  if (ltoadd > left) ltoadd = left;
1765  strncpy(&command[kcmnd],celmnt[ielmnt-1],ltoadd);
1766  kcmnd += ltoadd;
1767  if (kcmnd == maxcwd) continue;
1768  command[kcmnd] = ' ';
1769  ++kcmnd;
1770  command[kcmnd] = 0;
1771  }
1772  lnc = kcmnd;
1773  goto L900;
1774 L450:
1775  lnc = kcmnd;
1776 // we have come to a numeric field
1777  llist = 0;
1778  for (ifld = ielmnt; ifld <= nelmnt; ++ifld) {
1779  ++llist;
1780  if (llist > mxp) {
1781  nreq = nelmnt - ielmnt + 1;
1782  Printf(" MINUIT WARNING IN MNCRCK: ");
1783  Printf(" COMMAND HAS INPUT %5d NUMERIC FIELDS, BUT MINUIT CAN ACCEPT ONLY%3d",nreq,mxp);
1784  goto L900;
1785  }
1786  if (celmnt[ifld-1] == cnull) plist[llist-1] = 0;
1787  else {
1788  sscanf(celmnt[ifld-1],"%lf",&plist[llist-1]);
1789  }
1790  }
1791 // end loop over numeric fields
1792 L900:
1793  if (lnc <= 0) lnc = 1;
1794  comand = command;
1795 }
1796 
1797 ////////////////////////////////////////////////////////////////////////////////
1798 /// Find point where MNEVAL=AMIN+UP
1799 ///
1800 /// Find point where MNEVAL=AMIN+UP, along the line through
1801 /// XMIDCR,YMIDCR with direction XDIRCR,YDIRCR, where X and Y
1802 /// are parameters KE1CR and KE2CR. If KE2CR=0 (from MINOS),
1803 /// only KE1CR is varied. From MNCONT, both are varied.
1804 /// Crossing point is at
1805 ///
1806 /// (U(KE1),U(KE2)) = (XMID,YMID) + AOPT*(XDIR,YDIR)
1807 
1808 void TMinuit::mncros(Double_t &aopt, Int_t &iercr)
1809 {
1810  /* Local variables */
1811  Double_t alsb[3], flsb[3], bmin, bmax, zmid, sdev, zdir, zlim;
1812  Double_t coeff[3], aleft, aulim, fdist, adist, aminsv;
1813  Double_t anext, fnext, slope, s1, s2, x1, x2, ecarmn, ecarmx;
1814  Double_t determ, rt, smalla, aright, aim, tla, tlf, dfda,ecart;
1815  Int_t iout=0, i, ileft, ierev, maxlk, ibest, ik, it;
1816  Int_t noless, iworst=0, iright, itoohi, kex, ipt;
1817  Bool_t ldebug;
1818  const char *chsign;
1819  x2 = 0;
1820 
1821  ldebug = fIdbg[6] >= 1;
1822  aminsv = fAmin;
1823 // convergence when F is within TLF of AIM and next prediction
1824 // of AOPT is within TLA of previous value of AOPT
1825  aim = fAmin + fUp;
1826  tlf = fUp*.01;
1827  tla = .01;
1828  fXpt[0] = 0;
1829  fYpt[0] = aim;
1830  fChpt[0] = ' ';
1831  ipt = 1;
1832  if (fKe2cr == 0) {
1833  fXpt[1] = -1;
1834  fYpt[1] = fAmin;
1835  fChpt[1] = '.';
1836  ipt = 2;
1837  }
1838 // find the largest allowed A
1839  aulim = 100;
1840  for (ik = 1; ik <= 2; ++ik) {
1841  if (ik == 1) {
1842  kex = fKe1cr;
1843  zmid = fXmidcr;
1844  zdir = fXdircr;
1845  } else {
1846  if (fKe2cr == 0) continue;
1847  kex = fKe2cr;
1848  zmid = fYmidcr;
1849  zdir = fYdircr;
1850  }
1851  if (fNvarl[kex-1] <= 1) continue;
1852  if (zdir == 0) continue;
1853  zlim = fAlim[kex-1];
1854  if (zdir > 0) zlim = fBlim[kex-1];
1855  aulim = TMath::Min(aulim,(zlim - zmid) / zdir);
1856  }
1857 // LSB = Line Search Buffer
1858 // first point
1859  anext = 0;
1860  aopt = anext;
1861  fLimset = kFALSE;
1862  if (aulim < aopt + tla) fLimset = kTRUE;
1863  mneval(anext, fnext, ierev);
1864 // debug printout:
1865  if (ldebug) {
1866  Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
1867  }
1868  if (ierev > 0) goto L900;
1869  if (fLimset && fnext <= aim) goto L930;
1870  ++ipt;
1871  fXpt[ipt-1] = anext;
1872  fYpt[ipt-1] = fnext;
1873  fChpt[ipt-1] = charal[ipt-1];
1874  alsb[0] = anext;
1875  flsb[0] = fnext;
1876  fnext = TMath::Max(fnext,aminsv + fUp*.1);
1877  aopt = TMath::Sqrt(fUp / (fnext - aminsv)) - 1;
1878  if (TMath::Abs(fnext - aim) < tlf) goto L800;
1879 
1880  if (aopt < -.5)aopt = -.5;
1881  if (aopt > 1) aopt = 1;
1882  fLimset = kFALSE;
1883  if (aopt > aulim) {
1884  aopt = aulim;
1885  fLimset = kTRUE;
1886  }
1887  mneval(aopt, fnext, ierev);
1888 // debug printout:
1889  if (ldebug) {
1890  Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
1891  }
1892  if (ierev > 0) goto L900;
1893  if (fLimset && fnext <= aim) goto L930;
1894  alsb[1] = aopt;
1895  ++ipt;
1896  fXpt[ipt-1] = alsb[1];
1897  fYpt[ipt-1] = fnext;
1898  fChpt[ipt-1] = charal[ipt-1];
1899  flsb[1] = fnext;
1900  dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
1901 // DFDA must be positive on the contour
1902  if (dfda > 0) goto L460;
1903 L300:
1904  mnwarn("D", "MNCROS ", "Looking for slope of the right sign");
1905  maxlk = 15 - ipt;
1906  for (it = 1; it <= maxlk; ++it) {
1907  alsb[0] = alsb[1];
1908  flsb[0] = flsb[1];
1909  aopt = alsb[0] + Double_t(it)*.2;
1910  fLimset = kFALSE;
1911  if (aopt > aulim) {
1912  aopt = aulim;
1913  fLimset = kTRUE;
1914  }
1915  mneval(aopt, fnext, ierev);
1916 // debug printout:
1917  if (ldebug) {
1918  Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
1919  }
1920  if (ierev > 0) goto L900;
1921  if (fLimset && fnext <= aim) goto L930;
1922  alsb[1] = aopt;
1923  ++ipt;
1924  fXpt[ipt-1] = alsb[1];
1925  fYpt[ipt-1] = fnext;
1926  fChpt[ipt-1] = charal[ipt-1];
1927  flsb[1] = fnext;
1928  dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
1929  if (dfda > 0) goto L450;
1930  }
1931  mnwarn("W", "MNCROS ", "Cannot find slope of the right sign");
1932  goto L950;
1933 L450:
1934 // we have two points with the right slope
1935 L460:
1936  aopt = alsb[1] + (aim - flsb[1]) / dfda;
1937  fdist = TMath::Min(TMath::Abs(aim - flsb[0]),TMath::Abs(aim - flsb[1]));
1938  adist = TMath::Min(TMath::Abs(aopt - alsb[0]),TMath::Abs(aopt - alsb[1]));
1939  tla = .01;
1940  if (TMath::Abs(aopt) > 1) tla = TMath::Abs(aopt)*.01;
1941  if (adist < tla && fdist < tlf) goto L800;
1942  if (ipt >= 15) goto L950;
1943  bmin = TMath::Min(alsb[0],alsb[1]) - 1;
1944  if (aopt < bmin) aopt = bmin;
1945  bmax = TMath::Max(alsb[0],alsb[1]) + 1;
1946  if (aopt > bmax) aopt = bmax;
1947 // Try a third point
1948  fLimset = kFALSE;
1949  if (aopt > aulim) {
1950  aopt = aulim;
1951  fLimset = kTRUE;
1952  }
1953  mneval(aopt, fnext, ierev);
1954 // debug printout:
1955  if (ldebug) {
1956  Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
1957  }
1958  if (ierev > 0) goto L900;
1959  if (fLimset && fnext <= aim) goto L930;
1960  alsb[2] = aopt;
1961  ++ipt;
1962  fXpt[ipt-1] = alsb[2];
1963  fYpt[ipt-1] = fnext;
1964  fChpt[ipt-1] = charal[ipt-1];
1965  flsb[2] = fnext;
1966 // now we have three points, ask how many <AIM
1967  ecarmn = TMath::Abs(fnext-aim);
1968  ibest = 3;
1969  ecarmx = 0;
1970  noless = 0;
1971  for (i = 1; i <= 3; ++i) {
1972  ecart = TMath::Abs(flsb[i-1] - aim);
1973  if (ecart > ecarmx) { ecarmx = ecart; iworst = i; }
1974  if (ecart < ecarmn) { ecarmn = ecart; ibest = i; }
1975  if (flsb[i-1] < aim) ++noless;
1976  }
1977 // if at least one on each side of AIM, fit a parabola
1978  if (noless == 1 || noless == 2) goto L500;
1979 // if all three are above AIM, third must be closest to AIM
1980  if (noless == 0 && ibest != 3) goto L950;
1981 // if all three below, and third is not best, then slope
1982 // has again gone negative, look for positive slope.
1983  if (noless == 3 && ibest != 3) {
1984  alsb[1] = alsb[2];
1985  flsb[1] = flsb[2];
1986  goto L300;
1987  }
1988 // in other cases, new straight line thru last two points
1989  alsb[iworst-1] = alsb[2];
1990  flsb[iworst-1] = flsb[2];
1991  dfda = (flsb[1] - flsb[0]) / (alsb[1] - alsb[0]);
1992  goto L460;
1993 // parabola fit
1994 L500:
1995  mnpfit(alsb, flsb, 3, coeff, sdev);
1996  if (coeff[2] <= 0) {
1997  mnwarn("D", "MNCROS ", "Curvature is negative near contour line.");
1998  }
1999  determ = coeff[1]*coeff[1] - coeff[2]*4*(coeff[0] - aim);
2000  if (determ <= 0) {
2001  mnwarn("D", "MNCROS ", "Problem 2, impossible determinant");
2002  goto L950;
2003  }
2004 // Find which root is the right one
2005  rt = TMath::Sqrt(determ);
2006  x1 = (-coeff[1] + rt) / (coeff[2]*2);
2007  x2 = (-coeff[1] - rt) / (coeff[2]*2);
2008  s1 = coeff[1] + x1*2*coeff[2];
2009  s2 = coeff[1] + x2*2*coeff[2];
2010  if (s1*s2 > 0) {
2011  Printf(" MNCONTour problem 1");
2012  }
2013  aopt = x1;
2014  slope = s1;
2015  if (s2 > 0) {
2016  aopt = x2;
2017  slope = s2;
2018  }
2019 // ask if converged
2020  tla = .01;
2021  if (TMath::Abs(aopt) > 1) tla = TMath::Abs(aopt)*.01;
2022  if (TMath::Abs(aopt - alsb[ibest-1]) < tla && TMath::Abs(flsb[ibest-1] - aim) < tlf) {
2023  goto L800;
2024  }
2025  if (ipt >= 15) goto L950;
2026 
2027 // see if proposed point is in acceptable zone between L and R
2028 // first find ILEFT, IRIGHT, IOUT and IBEST
2029  ileft = 0;
2030  iright = 0;
2031  ibest = 1;
2032  ecarmx = 0;
2033  ecarmn = TMath::Abs(aim - flsb[0]);
2034  for (i = 1; i <= 3; ++i) {
2035  ecart = TMath::Abs(flsb[i-1] - aim);
2036  if (ecart < ecarmn) { ecarmn = ecart; ibest = i; }
2037  if (ecart > ecarmx) { ecarmx = ecart; }
2038  if (flsb[i-1] > aim) {
2039  if (iright == 0) iright = i;
2040  else if (flsb[i-1] > flsb[iright-1]) iout = i;
2041  else { iout = iright; iright = i; }
2042  }
2043  else if (ileft == 0) ileft = i;
2044  else if (flsb[i-1] < flsb[ileft-1]) iout = i;
2045  else { iout = ileft; ileft = i; }
2046  }
2047 // avoid keeping a very bad point next time around
2048  if (ecarmx > TMath::Abs(flsb[iout-1] - aim)*10) {
2049  aopt = aopt*.5 + (alsb[iright-1] + alsb[ileft-1])*.25;
2050  }
2051 // knowing ILEFT and IRIGHT, get acceptable window
2052  smalla = tla*.1;
2053  if (slope*smalla > tlf) smalla = tlf / slope;
2054  aleft = alsb[ileft-1] + smalla;
2055  aright = alsb[iright-1] - smalla;
2056 // move proposed point AOPT into window if necessary
2057  if (aopt < aleft) aopt = aleft;
2058  if (aopt > aright) aopt = aright;
2059  if (aleft > aright) aopt = (aleft + aright)*.5;
2060 
2061 // see if proposed point outside limits (should be impossible!)
2062  fLimset = kFALSE;
2063  if (aopt > aulim) {
2064  aopt = aulim;
2065  fLimset = kTRUE;
2066  }
2067 // Evaluate function at new point AOPT
2068  mneval(aopt, fnext, ierev);
2069 // debug printout:
2070  if (ldebug) {
2071  Printf(" MNCROS: calls=%8d AIM=%10.5f F,A=%10.5f%10.5f",fNfcn,aim,fnext,aopt);
2072  }
2073  if (ierev > 0) goto L900;
2074  if (fLimset && fnext <= aim) goto L930;
2075  ++ipt;
2076  fXpt[ipt-1] = aopt;
2077  fYpt[ipt-1] = fnext;
2078  fChpt[ipt-1] = charal[ipt-1];
2079 // Replace odd point by new one
2080  alsb[iout-1] = aopt;
2081  flsb[iout-1] = fnext;
2082 // the new point may not be the best, but it is the only one
2083 // which could be good enough to pass convergence criteria
2084  ibest = iout;
2085  goto L500;
2086 
2087 // Contour has been located, return point to MNCONT OR MINOS
2088 L800:
2089  iercr = 0;
2090  goto L1000;
2091 // error in the minimization
2092 L900:
2093  if (ierev == 1) goto L940;
2094  goto L950;
2095 // parameter up against limit
2096 L930:
2097  iercr = 1;
2098  goto L1000;
2099 // too many calls to FCN
2100 L940:
2101  iercr = 2;
2102  goto L1000;
2103 // cannot find next point
2104 L950:
2105  iercr = 3;
2106 // in any case
2107 L1000:
2108  if (ldebug) {
2109  itoohi = 0;
2110  for (i = 1; i <= ipt; ++i) {
2111  if (fYpt[i-1] > aim + fUp) {
2112  fYpt[i-1] = aim + fUp;
2113  fChpt[i-1] = '+';
2114  itoohi = 1;
2115  }
2116  }
2117  fChpt[ipt] = 0;
2118  chsign = "POSI";
2119  if (fXdircr < 0) chsign = "NEGA";
2120  if (fKe2cr == 0) {
2121  Printf(" %sTIVE MINOS ERROR, PARAMETER %3d",chsign,fKe1cr);
2122  }
2123  if (itoohi == 1) {
2124  Printf("POINTS LABELLED '+' WERE TOO HIGH TO PLOT.");
2125  }
2126  if (iercr == 1) {
2127  Printf("RIGHTMOST POINT IS UP AGAINST LIMIT.");
2128  }
2129  mnplot(fXpt, fYpt, fChpt, ipt, fNpagwd, fNpagln);
2130  }
2131 }
2132 
2133 ////////////////////////////////////////////////////////////////////////////////
2134 /// Makes sure that the current point is a local minimum
2135 ///
2136 /// Makes sure that the current point is a local
2137 /// minimum and that the error matrix exists,
2138 /// or at least something good enough for MINOS and MNCONT
2139 
2141 {
2142  /* Local variables */
2143  Double_t dxdi, wint;
2144  Int_t ndex, iext, i, j;
2145 
2146  if (fISW[3] < 1) {
2147  Printf(" FUNCTION MUST BE MINIMIZED BEFORE CALLING %s",(const char*)fCfrom);
2148  fApsi = fEpsi;
2149  mnmigr();
2150  }
2151  if (fISW[1] < 3) {
2152  mnhess();
2153  if (fISW[1] < 1) {
2154  mnwarn("W", fCfrom, "NO ERROR MATRIX. WILL IMPROVISE.");
2155  for (i = 1; i <= fNpar; ++i) {
2156  ndex = i*(i-1) / 2;
2157  for (j = 1; j <= i-1; ++j) {
2158  ++ndex;
2159  fVhmat[ndex-1] = 0;
2160  }
2161  ++ndex;
2162  if (fG2[i-1] <= 0) {
2163  wint = fWerr[i-1];
2164  iext = fNexofi[i-1];
2165  if (fNvarl[iext-1] > 1) {
2166  mndxdi(fX[i-1], i-1, dxdi);
2167  if (TMath::Abs(dxdi) < .001) wint = .01;
2168  else wint /= TMath::Abs(dxdi);
2169  }
2170  fG2[i-1] = fUp / (wint*wint);
2171  }
2172  fVhmat[ndex-1] = 2 / fG2[i-1];
2173  }
2174  fISW[1] = 1;
2175  fDcovar = 1;
2176  } else mnwerr();
2177  }
2178 }
2179 
2180 ////////////////////////////////////////////////////////////////////////////////
2181 /// Calculates the first derivatives of FCN (GRD)
2182 ///
2183 /// Calculates the first derivatives of FCN (GRD),
2184 /// either by finite differences or by transforming the user-
2185 /// supplied derivatives to internal coordinates,
2186 /// according to whether fISW[2] is zero or one.
2187 
2189 {
2190  /* Local variables */
2191  Double_t step, dfmin, stepb4, dd, df, fs1;
2192  Double_t tlrstp, tlrgrd, epspri, optstp, stpmax, stpmin, fs2, grbfor=0, d1d2, xtf;
2193  Int_t icyc, ncyc, iint, iext, i, nparx;
2194  Bool_t ldebug;
2195 
2196  nparx = fNpar;
2197  ldebug = fIdbg[2] >= 1;
2198  if (fAmin == fUndefi) mnamin();
2199  if (fISW[2] == 1) goto L100;
2200 
2201  if (ldebug) {
2202 // make sure starting at the right place
2203  mninex(fX);
2204  nparx = fNpar;
2205  Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
2206  if (fs1 != fAmin) {
2207  df = fAmin - fs1;
2208  mnwarn("D", "MNDERI", TString::Format("function value differs from AMIN by %12.3g",df));
2209  fAmin = fs1;
2210  }
2211  Printf(" FIRST DERIVATIVE DEBUG PRINTOUT. MNDERI");
2212  Printf(" PAR DERIV STEP MINSTEP OPTSTEP D1-D2 2ND DRV");
2213  }
2214  dfmin = fEpsma2*8*(TMath::Abs(fAmin) + fUp);
2215  if (fIstrat <= 0) {
2216  ncyc = 2;
2217  tlrstp = .5;
2218  tlrgrd = .1;
2219  } else if (fIstrat == 1) {
2220  ncyc = 3;
2221  tlrstp = .3;
2222  tlrgrd = .05;
2223  } else {
2224  ncyc = 5;
2225  tlrstp = .1;
2226  tlrgrd = .02;
2227  }
2228 // loop over variable parameters
2229  for (i = 1; i <= fNpar; ++i) {
2230  epspri = fEpsma2 + TMath::Abs(fGrd[i-1]*fEpsma2);
2231 // two-point derivatives always assumed necessary
2232 // maximum number of cycles over step size depends on strategy
2233  xtf = fX[i-1];
2234  stepb4 = 0;
2235 // loop as little as possible here!/
2236  for (icyc = 1; icyc <= ncyc; ++icyc) {
2237 // theoretically best step
2238  optstp = TMath::Sqrt(dfmin / (TMath::Abs(fG2[i-1]) + epspri));
2239 // step cannot decrease by more than a factor of ten
2240  step = TMath::Max(optstp,TMath::Abs(fGstep[i-1]*.1));
2241 // but if parameter has limits, max step size = 0.5
2242  if (fGstep[i-1] < 0 && step > .5) step = .5;
2243 // and not more than ten times the previous step
2244  stpmax = TMath::Abs(fGstep[i-1])*10;
2245  if (step > stpmax) step = stpmax;
2246 // minimum step size allowed by machine precision
2247  stpmin = TMath::Abs(fEpsma2*fX[i-1])*8;
2248  if (step < stpmin) step = stpmin;
2249 // end of iterations if step change less than factor 2
2250  if (TMath::Abs((step - stepb4) / step) < tlrstp) goto L50;
2251 // take step positive
2252  stepb4 = step;
2253  if (fGstep[i-1] > 0) fGstep[i-1] = TMath::Abs(step);
2254  else fGstep[i-1] = -TMath::Abs(step);
2255  stepb4 = step;
2256  fX[i-1] = xtf + step;
2257  mninex(fX);
2258  Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
2259 // take step negative
2260  fX[i-1] = xtf - step;
2261  mninex(fX);
2262  Eval(nparx, fGin, fs2, fU, 4); ++fNfcn;
2263  grbfor = fGrd[i-1];
2264  fGrd[i-1] = (fs1 - fs2) / (step*2);
2265  fG2[i-1] = (fs1 + fs2 - fAmin*2) / (step*step);
2266  fX[i-1] = xtf;
2267  if (ldebug) {
2268  d1d2 = (fs1 + fs2 - fAmin*2) / step;
2269  Printf("%4d%11.3g%11.3g%10.2g%10.2g%10.2g%10.2g",i,fGrd[i-1],step,stpmin,optstp,d1d2,fG2[i-1]);
2270  }
2271 // see if another iteration is necessary
2272  if (TMath::Abs(grbfor - fGrd[i-1]) / (TMath::Abs(fGrd[i-1]) + dfmin/step) < tlrgrd)
2273  goto L50;
2274  }
2275 // end of ICYC loop. too many iterations
2276  if (ncyc == 1) goto L50;
2277  mnwarn("D", "MNDERI", TString::Format("First derivative not converged. %g%g",fGrd[i-1],grbfor));
2278 L50:
2279  ;
2280  }
2281  mninex(fX);
2282  return;
2283 // derivatives calc by fcn
2284 L100:
2285  for (iint = 1; iint <= fNpar; ++iint) {
2286  iext = fNexofi[iint-1];
2287  if (fNvarl[iext-1] <= 1) {
2288  fGrd[iint-1] = fGin[iext-1];
2289  } else {
2290  dd = (fBlim[iext-1] - fAlim[iext-1])*.5*TMath::Cos(fX[iint-1]);
2291  fGrd[iint-1] = fGin[iext-1]*dd;
2292  }
2293  }
2294 }
2295 
2296 ////////////////////////////////////////////////////////////////////////////////
2297 /// Calculates the transformation factor between ext/internal values
2298 ///
2299 /// calculates the transformation factor between external and
2300 /// internal parameter values. this factor is one for
2301 /// parameters which are not limited. called from MNEMAT.
2302 
2303 void TMinuit::mndxdi(Double_t pint, Int_t ipar, Double_t &dxdi)
2304 {
2305  Int_t i = fNexofi[ipar];
2306  dxdi = 1;
2307  if (fNvarl[i-1] > 1) {
2308  dxdi = TMath::Abs((fBlim[i-1] - fAlim[i-1])*TMath::Cos(pint))*.5;
2309  }
2310 }
2311 
2312 ////////////////////////////////////////////////////////////////////////////////
2313 /// Compute matrix eigen values
2314 
2315 void TMinuit::mneig(Double_t *a, Int_t ndima, Int_t n, Int_t mits, Double_t *work, Double_t precis, Int_t &ifault)
2316 {
2317  /* System generated locals */
2318  Int_t a_offset;
2319  Double_t d__1;
2320 
2321  /* Local variables */
2322  Double_t b, c, f, h, r, s, hh, gl, pr, pt;
2323  Int_t i, j, k, l, m=0, i0, i1, j1, m1, n1;
2324 
2325 // PRECIS is the machine precision EPSMAC
2326  /* Parameter adjustments */
2327  a_offset = ndima + 1;
2328  a -= a_offset;
2329  --work;
2330 
2331  /* Function Body */
2332  ifault = 1;
2333 
2334  i = n;
2335  for (i1 = 2; i1 <= n; ++i1) {
2336  l = i-2;
2337  f = a[i + (i-1)*ndima];
2338  gl = 0;
2339 
2340  if (l < 1) goto L25;
2341 
2342  for (k = 1; k <= l; ++k) {
2343  d__1 = a[i + k*ndima];
2344  gl += d__1*d__1;
2345  }
2346 L25:
2347  h = gl + f*f;
2348 
2349  if (gl > 1e-35) goto L30;
2350 
2351  work[i] = 0;
2352  work[n + i] = f;
2353  goto L65;
2354 L30:
2355  ++l;
2356  gl = TMath::Sqrt(h);
2357  if (f >= 0) gl = -gl;
2358  work[n + i] = gl;
2359  h -= f*gl;
2360  a[i + (i-1)*ndima] = f - gl;
2361  f = 0;
2362  for (j = 1; j <= l; ++j) {
2363  a[j + i*ndima] = a[i + j*ndima] / h;
2364  gl = 0;
2365  for (k = 1; k <= j; ++k) { gl += a[j + k*ndima]*a[i + k*ndima]; }
2366  if (j >= l) goto L47;
2367  j1 = j + 1;
2368  for (k = j1; k <= l; ++k) { gl += a[k + j*ndima]*a[i + k*ndima]; }
2369 L47:
2370  work[n + j] = gl / h;
2371  f += gl*a[j + i*ndima];
2372  }
2373  hh = f / (h + h);
2374  for (j = 1; j <= l; ++j) {
2375  f = a[i + j*ndima];
2376  gl = work[n + j] - hh*f;
2377  work[n + j] = gl;
2378  for (k = 1; k <= j; ++k) {
2379  a[j + k*ndima] = a[j + k*ndima] - f*work[n + k] - gl*a[i + k*ndima];
2380  }
2381  }
2382  work[i] = h;
2383 L65:
2384  --i;
2385  }
2386  work[1] = 0;
2387  work[n + 1] = 0;
2388  for (i = 1; i <= n; ++i) {
2389  l = i-1;
2390  if (work[i] == 0 || l == 0) goto L100;
2391 
2392  for (j = 1; j <= l; ++j) {
2393  gl = 0;
2394  for (k = 1; k <= l; ++k) { gl += a[i + k*ndima]*a[k + j*ndima]; }
2395  for (k = 1; k <= l; ++k) { a[k + j*ndima] -= gl*a[k + i*ndima]; }
2396  }
2397 L100:
2398  work[i] = a[i + i*ndima];
2399  a[i + i*ndima] = 1;
2400  if (l == 0) continue;
2401 
2402  for (j = 1; j <= l; ++j) {
2403  a[i + j*ndima] = 0;
2404  a[j + i*ndima] = 0;
2405  }
2406  }
2407 
2408  n1 = n - 1;
2409  for (i = 2; i <= n; ++i) {
2410  i0 = n + i-1;
2411  work[i0] = work[i0 + 1];
2412  }
2413  work[n + n] = 0;
2414  b = 0;
2415  f = 0;
2416  for (l = 1; l <= n; ++l) {
2417  j = 0;
2418  h = precis*(TMath::Abs(work[l]) + TMath::Abs(work[n + l]));
2419  if (b < h) b = h;
2420  for (m1 = l; m1 <= n; ++m1) {
2421  m = m1;
2422  if (TMath::Abs(work[n + m]) <= b) goto L150;
2423  }
2424 
2425 L150:
2426  if (m == l) goto L205;
2427 
2428 L160:
2429  if (j == mits) return;
2430  ++j;
2431  pt = (work[l + 1] - work[l]) / (work[n + l]*2);
2432  r = TMath::Sqrt(pt*pt + 1);
2433  pr = pt + r;
2434  if (pt < 0) pr = pt - r;
2435 
2436  h = work[l] - work[n + l] / pr;
2437  for (i = l; i <= n; ++i) { work[i] -= h; }
2438  f += h;
2439  pt = work[m];
2440  c = 1;
2441  s = 0;
2442  m1 = m - 1;
2443  i = m;
2444  for (i1 = l; i1 <= m1; ++i1) {
2445  j = i;
2446  --i;
2447  gl = c*work[n + i];
2448  h = c*pt;
2449  if (TMath::Abs(pt) >= TMath::Abs(work[n + i])) goto L180;
2450 
2451  c = pt / work[n + i];
2452  r = TMath::Sqrt(c*c + 1);
2453  work[n + j] = s*work[n + i]*r;
2454  s = 1 / r;
2455  c /= r;
2456  goto L190;
2457 L180:
2458  c = work[n + i] / pt;
2459  r = TMath::Sqrt(c*c + 1);
2460  work[n + j] = s*pt*r;
2461  s = c / r;
2462  c = 1 / r;
2463 L190:
2464  pt = c*work[i] - s*gl;
2465  work[j] = h + s*(c*gl + s*work[i]);
2466  for (k = 1; k <= n; ++k) {
2467  h = a[k + j*ndima];
2468  a[k + j*ndima] = s*a[k + i*ndima] + c*h;
2469  a[k + i*ndima] = c*a[k + i*ndima] - s*h;
2470  }
2471  }
2472  work[n + l] = s*pt;
2473  work[l] = c*pt;
2474 
2475  if (TMath::Abs(work[n + l]) > b) goto L160;
2476 
2477 L205:
2478  work[l] += f;
2479  }
2480  for (i = 1; i <= n1; ++i) {
2481  k = i;
2482  pt = work[i];
2483  i1 = i + 1;
2484  for (j = i1; j <= n; ++j) {
2485  if (work[j] >= pt) continue;
2486  k = j;
2487  pt = work[j];
2488  }
2489 
2490  if (k == i) continue;
2491 
2492  work[k] = work[i];
2493  work[i] = pt;
2494  for (j = 1; j <= n; ++j) {
2495  pt = a[j + i*ndima];
2496  a[j + i*ndima] = a[j + k*ndima];
2497  a[j + k*ndima] = pt;
2498  }
2499  }
2500  ifault = 0;
2501 }
2502 
2503 ////////////////////////////////////////////////////////////////////////////////
2504 /// Calculates the external error matrix from the internal matrix
2505 ///
2506 /// Note that if the matrix is declared like Double_t matrix[5][5]
2507 /// in the calling program, one has to call mnemat with, eg
2508 ///
2509 /// gMinuit->mnemat(&matrix[0][0],5);
2510 
2511 void TMinuit::mnemat(Double_t *emat, Int_t ndim)
2512 {
2513  /* System generated locals */
2514  Int_t emat_dim1, emat_offset;
2515 
2516  /* Local variables */
2517  Double_t dxdi, dxdj;
2518  Int_t i, j, k, npard, k2, kk, iz, nperln, kga, kgb;
2519  TString ctemp;
2520 
2521  /* Parameter adjustments */
2522  emat_dim1 = ndim;
2523  emat_offset = emat_dim1 + 1;
2524  emat -= emat_offset;
2525 
2526  /* Function Body */
2527  if (fISW[1] < 1) return;
2528  if (fISW[4] >= 2) {
2529  Printf(" EXTERNAL ERROR MATRIX. NDIM=%4d NPAR=%3d ERR DEF=%g",ndim,fNpar,fUp);
2530  }
2531 // size of matrix to be printed
2532  npard = fNpar;
2533  if (ndim < fNpar) {
2534  npard = ndim;
2535  if (fISW[4] >= 0) {
2536  Printf(" USER-DIMENSIONED ARRAY EMAT NOT BIG ENOUGH. REDUCED MATRIX CALCULATED.");
2537  }
2538  }
2539 // NPERLN is the number of elements that fit on one line
2540 
2541  nperln = (fNpagwd - 5) / 10;
2542  nperln = TMath::Min(nperln,13);
2543  if (fISW[4] >= 1 && npard > nperln) {
2544  Printf(" ELEMENTS ABOVE DIAGONAL ARE NOT PRINTED.");
2545  }
2546 // I counts the rows of the matrix
2547  for (i = 1; i <= npard; ++i) {
2548  mndxdi(fX[i-1], i-1, dxdi);
2549  kga = i*(i-1) / 2;
2550  for (j = 1; j <= i; ++j) {
2551  mndxdi(fX[j-1], j-1, dxdj);
2552  kgb = kga + j;
2553  emat[i + j*emat_dim1] = dxdi*fVhmat[kgb-1]*dxdj*fUp;
2554  emat[j + i*emat_dim1] = emat[i + j*emat_dim1];
2555  }
2556  }
2557 // IZ is number of columns to be printed in row I
2558  if (fISW[4] >= 2) {
2559  for (i = 1; i <= npard; ++i) {
2560  iz = npard;
2561  if (npard >= nperln) iz = i;
2562  ctemp = " ";
2563  for (k = 1; nperln < 0 ? k >= iz : k <= iz; k += nperln) {
2564  k2 = k + nperln - 1;
2565  if (k2 > iz) k2 = iz;
2566  for (kk = k; kk <= k2; ++kk) {
2567  ctemp += TString::Format("%10.3e ",emat[i + kk*emat_dim1]);
2568  }
2569  Printf("%s",(const char*)ctemp);
2570  }
2571  }
2572  }
2573 }
2574 
2575 ////////////////////////////////////////////////////////////////////////////////
2576 /// Utility routine to get MINOS errors
2577 ///
2578 /// Called by user.
2579 ///
2580 /// NUMBER is the parameter number
2581 ///
2582 /// values returned by MNERRS:
2583 /// - EPLUS, EMINUS are MINOS errors of parameter NUMBER,
2584 /// - EPARAB is 'parabolic' error (from error matrix).
2585 /// (Errors not calculated are set = 0)
2586 /// - GCC is global correlation coefficient from error matrix
2587 
2588 void TMinuit::mnerrs(Int_t number, Double_t &eplus, Double_t &eminus, Double_t &eparab, Double_t &gcc)
2589 {
2590  Double_t dxdi;
2591  Int_t ndiag, iin, iex;
2592 
2593  iex = number+1;
2594 
2595  if (iex > fNu || iex <= 0) goto L900;
2596  iin = fNiofex[iex-1];
2597  if (iin <= 0) goto L900;
2598 
2599 // IEX is external number, IIN is internal number
2600  eplus = fErp[iin-1];
2601  if (eplus == fUndefi) eplus = 0;
2602  eminus = fErn[iin-1];
2603  if (eminus == fUndefi) eminus = 0;
2604  mndxdi(fX[iin-1], iin-1, dxdi);
2605  ndiag = iin*(iin + 1) / 2;
2606  eparab = TMath::Abs(dxdi*TMath::Sqrt(TMath::Abs(fUp*fVhmat[ndiag- 1])));
2607 // global correlation coefficient
2608  gcc = 0;
2609  if (fISW[1] < 2) return;
2610  gcc = fGlobcc[iin-1];
2611  return;
2612 // ERROR. parameter number not valid
2613 L900:
2614  eplus = 0;
2615  eminus = 0;
2616  eparab = 0;
2617  gcc = 0;
2618 }
2619 
2620 ////////////////////////////////////////////////////////////////////////////////
2621 /// Evaluates the function being analysed by MNCROS
2622 ///
2623 /// Evaluates the function being analysed by MNCROS, which is
2624 /// generally the minimum of FCN with respect to all remaining
2625 /// variable parameters. The class data members contains the
2626 /// data necessary to know the values of U(KE1CR) and U(KE2CR)
2627 /// to be used, namely U(KE1CR) = XMIDCR + ANEXT*XDIRCR
2628 /// and (if KE2CR .NE. 0) U(KE2CR) = YMIDCR + ANEXT*YDIRCR
2629 
2630 void TMinuit::mneval(Double_t anext, Double_t &fnext, Int_t &ierev)
2631 {
2632  Int_t nparx;
2633 
2634  fU[fKe1cr-1] = fXmidcr + anext*fXdircr;
2635  if (fKe2cr != 0) fU[fKe2cr-1] = fYmidcr + anext*fYdircr;
2636  mninex(fX);
2637  nparx = fNpar;
2638  Eval(nparx, fGin, fnext, fU, 4); ++fNfcn;
2639  ierev = 0;
2640  if (fNpar > 0) {
2641  fItaur = 1;
2642  fAmin = fnext;
2643  fISW[0] = 0;
2644  mnmigr();
2645  fItaur = 0;
2646  fnext = fAmin;
2647  if (fISW[0] >= 1) ierev = 1;
2648  if (fISW[3] < 1) ierev = 2;
2649  }
2650 }
2651 
2652 ////////////////////////////////////////////////////////////////////////////////
2653 /// Interprets a command and takes appropriate action
2654 ///
2655 /// either directly by skipping to the corresponding code in
2656 /// MNEXCM, or by setting up a call to a function
2657 ///
2658 /// recognized MINUIT commands:
2659 /// obsolete commands:
2660 /// IERFLG is now (94.5) defined the same as ICONDN in MNCOMD =
2661 /// - 0: command executed normally
2662 /// - 1: command is blank, ignored
2663 /// - 2: command line unreadable, ignored
2664 /// - 3: unknown command, ignored
2665 /// - 4: abnormal termination (e.g., MIGRAD not converged)
2666 /// - 9: reserved
2667 /// - 10: END command
2668 /// - 11: EXIT or STOP command
2669 /// - 12: RETURN command
2670 ///
2671 /// see also
2672 /// [the possible list of all Minuit commands](https://root.cern.ch/sites/d35c7d8c.web.cern.ch/files/minuit.pdf).
2673 
2674 void TMinuit::mnexcm(const char *command, Double_t *plist, Int_t llist, Int_t &ierflg)
2675 {
2676  /* Initialized data */
2677 
2678  TString comand = command;
2679  static const char *const cname[40] = {
2680  "MINImize ",
2681  "SEEk ",
2682  "SIMplex ",
2683  "MIGrad ",
2684  "MINOs ",
2685  "SET xxx ",
2686  "SHOw xxx ",
2687  "TOP of pag",
2688  "FIX ",
2689  "REStore ",
2690  "RELease ",
2691  "SCAn ",
2692  "CONtour ",
2693  "HESse ",
2694  "SAVe ",
2695  "IMProve ",
2696  "CALl fcn ",
2697  "STAndard ",
2698  "END ",
2699  "EXIt ",
2700  "RETurn ",
2701  "CLEar ",
2702  "HELP ",
2703  "MNContour ",
2704  "STOp ",
2705  "JUMp ",
2706  " ",
2707  " ",
2708  " ",
2709  " ",
2710  " ",
2711  " ",
2712  " ",
2713  "COVARIANCE",
2714  "PRINTOUT ",
2715  "GRADIENT ",
2716  "MATOUT ",
2717  "ERROR DEF ",
2718  "LIMITS ",
2719  "PUNCH "};
2720 
2721  Int_t nntot = 40;
2722 
2723  /* Local variables */
2724  Double_t step, xptu[101], yptu[101], f, rno;
2725  Int_t icol, kcol, ierr, iint, iext, lnow, nptu, i, iflag, ierrf;
2726  Int_t ilist, nparx, izero, nf, lk, it, iw, inonde, nsuper;
2727  Int_t it2, ke1, ke2, nowprt, kll, krl;
2728  TString chwhy, c26, cvblnk, cneway, comd;
2729  TString ctemp;
2730  Bool_t lfreed, ltofix, lfixed;
2731 
2732 // alphabetical order of command names!
2733 
2734  /* Function Body */
2735 
2736  lk = comand.Length();
2737  if (lk > 20) lk = 20;
2738  fCword = comand;
2739  fCword.ToUpper();
2740 // Copy the first MAXP arguments into WORD7, making
2741 // sure that WORD7(1)=0 if LLIST=0
2742  for (iw = 1; iw <= fMaxpar; ++iw) {
2743  fWord7[iw-1] = 0;
2744  if (iw <= llist) fWord7[iw-1] = plist[iw-1];
2745  }
2746  ++fIcomnd;
2747  fNfcnlc = fNfcn;
2748  if (fCword(0,7) != "SET PRI" || fWord7[0] >= 0) {
2749  if (fISW[4] >= 0) {
2750  lnow = llist;
2751  if (lnow > 4) lnow = 4;
2752  Printf(" **********");
2753  ctemp.Form(" **%5d **%s",fIcomnd,(const char*)fCword);
2754  for (i = 1; i <= lnow; ++i) {
2755  ctemp += TString::Format("%12.4g",plist[i-1]);
2756  }
2757  Printf("%s",(const char*)ctemp);
2758  inonde = 0;
2759  if (llist > lnow) {
2760  kll = llist;
2761  if (llist > fMaxpar) {
2762  inonde = 1;
2763  kll = fMaxpar;
2764  }
2765  Printf(" ***********");
2766  for (i = lnow + 1; i <= kll; ++i) {
2767  Printf("%12.4g",plist[i-1]);
2768  }
2769  }
2770  Printf(" **********");
2771  if (inonde > 0) {
2772  Printf(" ERROR: ABOVE CALL TO MNEXCM TRIED TO PASS MORE THAN %d PARAMETERS.", fMaxpar);
2773  }
2774  }
2775  }
2776  fNfcnmx = Int_t(fWord7[0]);
2777  if (fNfcnmx <= 0) {
2778  fNfcnmx = fNpar*100 + 200 + fNpar*fNpar*5;
2779  }
2780  fEpsi = fWord7[1];
2781  if (fEpsi <= 0) {
2782  fEpsi = fUp*.1;
2783  }
2784  fLnewmn = kFALSE;
2785  fLphead = kTRUE;
2786  fISW[0] = 0;
2787  ierflg = 0;
2788 // look for command in list CNAME
2789  ctemp = fCword(0,3);
2790  for (i = 1; i <= nntot; ++i) {
2791  if (strncmp(ctemp.Data(),cname[i-1],3) == 0) goto L90;
2792  }
2793  Printf("UNKNOWN COMMAND IGNORED:%s", comand.Data());
2794  ierflg = 3;
2795  return;
2796 // normal case: recognized MINUIT command
2797 L90:
2798  if (fCword(0,4) == "MINO") i = 5;
2799  if (i != 6 && i != 7 && i != 8 && i != 23) {
2800  fCfrom = cname[i-1];
2801  fNfcnfr = fNfcn;
2802  }
2803 // 1 2 3 4 5 6 7 8 9 10
2804  switch (i) {
2805  case 1: goto L400;
2806  case 2: goto L200;
2807  case 3: goto L300;
2808  case 4: goto L400;
2809  case 5: goto L500;
2810  case 6: goto L700;
2811  case 7: goto L700;
2812  case 8: goto L800;
2813  case 9: goto L900;
2814  case 10: goto L1000;
2815  case 11: goto L1100;
2816  case 12: goto L1200;
2817  case 13: goto L1300;
2818  case 14: goto L1400;
2819  case 15: goto L1500;
2820  case 16: goto L1600;
2821  case 17: goto L1700;
2822  case 18: goto L1800;
2823  case 19: goto L1900;
2824  case 20: goto L1900;
2825  case 21: goto L1900;
2826  case 22: goto L2200;
2827  case 23: goto L2300;
2828  case 24: goto L2400;
2829  case 25: goto L1900;
2830  case 26: goto L2600;
2831  case 27: goto L3300;
2832  case 28: goto L3300;
2833  case 29: goto L3300;
2834  case 30: goto L3300;
2835  case 31: goto L3300;
2836  case 32: goto L3300;
2837  case 33: goto L3300;
2838  case 34: goto L3400;
2839  case 35: goto L3500;
2840  case 36: goto L3600;
2841  case 37: goto L3700;
2842  case 38: goto L3800;
2843  case 39: goto L3900;
2844  case 40: goto L4000;
2845  }
2846 // seek
2847 L200:
2848  mnseek();
2849  return;
2850 // simplex
2851 L300:
2852  mnsimp();
2853  if (fISW[3] < 1) ierflg = 4;
2854  return;
2855 // migrad, minimize
2856 L400:
2857  nf = fNfcn;
2858  fApsi = fEpsi;
2859  mnmigr();
2860  mnwerr();
2861  if (fISW[3] >= 1) return;
2862  ierflg = 4;
2863  if (fISW[0] == 1) return;
2864  if (fCword(0,3) == "MIG") return;
2865 
2866  fNfcnmx = fNfcnmx + nf - fNfcn;
2867  nf = fNfcn;
2868  mnsimp();
2869  if (fISW[0] == 1) return;
2870  fNfcnmx = fNfcnmx + nf - fNfcn;
2871  mnmigr();
2872  if (fISW[3] >= 1) ierflg = 0;
2873  mnwerr();
2874  return;
2875 // minos
2876 L500:
2877  nsuper = fNfcn + ((fNpar + 1) << 1)*fNfcnmx;
2878 // possible loop over new minima
2879  fEpsi = fUp*.1;
2880 L510:
2881  fCfrom = cname[i-1]; // ensure that mncuve complains about MINOS not MIGRAD
2882  mncuve();
2883  mnmnos();
2884  if (! fLnewmn) return;
2885  mnrset(0);
2886  mnmigr();
2887  mnwerr();
2888  if (fNfcn < nsuper) goto L510;
2889  Printf(" TOO MANY FUNCTION CALLS. MINOS GIVES UP");
2890  ierflg = 4;
2891  return;
2892 // set, show
2893 L700:
2894  mnset();
2895  return;
2896 // top of page
2897 
2898 L800:
2899  Printf("1");
2900  return;
2901 // fix
2902 L900:
2903  ltofix = kTRUE;
2904 // (also release)
2905 L901:
2906  lfreed = kFALSE;
2907  lfixed = kFALSE;
2908  if (llist == 0) {
2909  Printf("%s: NO PARAMETERS REQUESTED ",(const char*)fCword);
2910  return;
2911  }
2912  for (ilist = 1; ilist <= llist; ++ilist) {
2913  iext = Int_t(plist[ilist-1]);
2914  chwhy = " IS UNDEFINED.";
2915  if (iext <= 0) goto L930;
2916  if (iext > fNu) goto L930;
2917  if (fNvarl[iext-1] < 0) goto L930;
2918  chwhy = " IS CONSTANT. ";
2919  if (fNvarl[iext-1] == 0) goto L930;
2920  iint = fNiofex[iext-1];
2921  if (ltofix) {
2922  chwhy = " ALREADY FIXED.";
2923  if (iint == 0) goto L930;
2924  mnfixp(iint-1, ierr);
2925  if (ierr == 0) lfixed = kTRUE;
2926  else ierflg = 4;
2927  } else {
2928  chwhy = " ALREADY VARIABLE.";
2929  if (iint > 0) goto L930;
2930  krl = -abs(iext);
2931  mnfree(krl);
2932  lfreed = kTRUE;
2933  }
2934  continue;
2935 L930:
2936  if (fISW[4] >= 0) Printf(" PARAMETER %4d %s IGNORED.",iext,(const char*)chwhy);
2937  }
2938  if (lfreed || lfixed) mnrset(0);
2939  if (lfreed) {
2940  fISW[1] = 0;
2941  fDcovar = 1;
2942  fEDM = fBigedm;
2943  fISW[3] = 0;
2944  }
2945  mnwerr();
2946  if (fISW[4] > 1) mnprin(5, fAmin);
2947  return;
2948 // restore
2949 L1000:
2950  it = Int_t(fWord7[0]);
2951  if (it > 1 || it < 0) goto L1005;
2952  lfreed = fNpfix > 0;
2953  mnfree(it);
2954  if (lfreed) {
2955  mnrset(0);
2956  fISW[1] = 0;
2957  fDcovar = 1;
2958  fEDM = fBigedm;
2959  }
2960  return;
2961 L1005:
2962  Printf(" IGNORED. UNKNOWN ARGUMENT:%4d",it);
2963  ierflg = 3;
2964  return;
2965 // release
2966 L1100:
2967  ltofix = kFALSE;
2968  goto L901;
2969 // scan
2970 L1200:
2971  iext = Int_t(fWord7[0]);
2972  if (iext <= 0) goto L1210;
2973  it2 = 0;
2974  if (iext <= fNu) it2 = fNiofex[iext-1];
2975  if (it2 <= 0) goto L1250;
2976 
2977 L1210:
2978  mnscan();
2979  return;
2980 L1250:
2981  Printf(" PARAMETER %4d NOT VARIABLE.",iext);
2982  ierflg = 3;
2983  return;
2984 // contour
2985 L1300:
2986  ke1 = Int_t(fWord7[0]);
2987  ke2 = Int_t(fWord7[1]);
2988  if (ke1 == 0) {
2989  if (fNpar == 2) {
2990  ke1 = fNexofi[0];
2991  ke2 = fNexofi[1];
2992  } else {
2993  Printf("%s: NO PARAMETERS REQUESTED ",(const char*)fCword);
2994  ierflg = 3;
2995  return;
2996  }
2997  }
2998  fNfcnmx = 1000;
2999  mncntr(ke1-1, ke2-1, ierrf);
3000  if (ierrf > 0) ierflg = 3;
3001  return;
3002 // hesse
3003 L1400:
3004  mnhess();
3005  mnwerr();
3006  if (fISW[4] >= 0) mnprin(2, fAmin);
3007  if (fISW[4] >= 1) mnmatu(1);
3008  return;
3009 // save
3010 L1500:
3011  mnsave();
3012  return;
3013 // improve
3014 L1600:
3015  mncuve();
3016  mnimpr();
3017  if (fLnewmn) goto L400;
3018  ierflg = 4;
3019  return;
3020 // call fcn
3021 L1700:
3022  iflag = Int_t(fWord7[0]);
3023  nparx = fNpar;
3024  f = fUndefi;
3025  Eval(nparx, fGin, f, fU, iflag); ++fNfcn;
3026  nowprt = 0;
3027  if (f != fUndefi) {
3028  if (fAmin == fUndefi) {
3029  fAmin = f;
3030  nowprt = 1;
3031  } else if (f < fAmin) {
3032  fAmin = f;
3033  nowprt = 1;
3034  }
3035  if (fISW[4] >= 0 && iflag <= 5 && nowprt == 1) {
3036  mnprin(5, fAmin);
3037  }
3038  if (iflag == 3) fFval3 = f;
3039  }
3040  if (iflag > 5) mnrset(1);
3041  return;
3042 // standard
3043 L1800:
3044 // stand();
3045  return;
3046 // return, stop, end, exit
3047 L1900:
3048  it = Int_t(fWord7[0]);
3049  if (fFval3 != fAmin && it == 0) {
3050  iflag = 3;
3051  if (fISW[4] >= 0) Printf(" CALL TO USER FUNCTION WITH IFLAG = 3");
3052  nparx = fNpar;
3053  Eval(nparx, fGin, f, fU, iflag); ++fNfcn;
3054  }
3055  ierflg = 11;
3056  if (fCword(0,3) == "END") ierflg = 10;
3057  if (fCword(0,3) == "RET") ierflg = 12;
3058  return;
3059 // clear
3060 L2200:
3061  mncler();
3062  if (fISW[4] >= 1) {
3063  Printf(" MINUIT MEMORY CLEARED. NO PARAMETERS NOW DEFINED.");
3064  }
3065  return;
3066 // help
3067 L2300:
3068  kcol = 0;
3069  for (icol = 5; icol <= lk; ++icol) {
3070  if (fCword[icol-1] == ' ') continue;
3071  kcol = icol;
3072  goto L2320;
3073  }
3074 L2320:
3075  if (kcol == 0) comd = "* ";
3076  else comd = fCword(kcol-1,lk-kcol+1);
3077  mnhelp(comd);
3078  return;
3079 // MNContour
3080 L2400:
3081  fEpsi = fUp*.05;
3082  ke1 = Int_t(fWord7[0]);
3083  ke2 = Int_t(fWord7[1]);
3084  if (ke1 == 0 && fNpar == 2) {
3085  ke1 = fNexofi[0];
3086  ke2 = fNexofi[1];
3087  }
3088  nptu = Int_t(fWord7[2]);
3089  if (nptu <= 0) nptu = 20;
3090  if (nptu > 101) nptu = 101;
3091  fNfcnmx = (nptu + 5)*100*(fNpar + 1);
3092  mncont(ke1-1, ke2-1, nptu, xptu, yptu, ierrf);
3093  if (ierrf < nptu) ierflg = 4;
3094  if (ierrf == -1) ierflg = 3;
3095  return;
3096 // jump
3097 L2600:
3098  step = fWord7[0];
3099  if (step <= 0) step = 2;
3100  rno = 0;
3101  izero = 0;
3102  for (i = 1; i <= fNpar; ++i) {
3103  mnrn15(rno, izero);
3104  rno = rno*2 - 1;
3105  fX[i-1] += rno*step*fWerr[i-1];
3106  }
3107  mninex(fX);
3108  mnamin();
3109  mnrset(0);
3110  return;
3111 // blank line
3112 L3300:
3113  Printf(" BLANK COMMAND IGNORED.");
3114  ierflg = 1;
3115  return;
3116 // obsolete commands
3117 // covariance
3118 L3400:
3119  Printf(" THE *COVARIANCE* COMMAND IS OSBSOLETE. THE COVARIANCE MATRIX IS NOW SAVED IN A DIFFERENT FORMAT WITH THE *SAVE* COMMAND AND READ IN WITH:*SET COVARIANCE*");
3120  ierflg = 3;
3121  return;
3122 // printout
3123 L3500:
3124  cneway = "SET PRInt ";
3125  goto L3100;
3126 // gradient
3127 L3600:
3128  cneway = "SET GRAd ";
3129  goto L3100;
3130 // matout
3131 L3700:
3132  cneway = "SHOW COVar";
3133  goto L3100;
3134 // error def
3135 L3800:
3136  cneway = "SET ERRdef";
3137  goto L3100;
3138 // limits
3139 L3900:
3140  cneway = "SET LIMits";
3141  goto L3100;
3142 // punch
3143 L4000:
3144  cneway = "SAVE ";
3145 // come from obsolete commands
3146 L3100:
3147  Printf(" OBSOLETE COMMAND:%s PLEASE USE:%s",(const char*)fCword
3148  ,(const char*)cneway);
3149  fCword = cneway;
3150  if (fCword == "SAVE ") goto L1500;
3151  goto L700;
3152 //
3153 }
3154 
3155 ////////////////////////////////////////////////////////////////////////////////
3156 /// Transforms the external parameter values U to internal values
3157 ///
3158 /// Transforms the external parameter values U to internal
3159 /// values in the dense array PINT.
3160 
3162 {
3163  Double_t pinti;
3164  Int_t iint, iext;
3165 
3166  fLimset = kFALSE;
3167  for (iint = 1; iint <= fNpar; ++iint) {
3168  iext = fNexofi[iint-1];
3169  mnpint(fU[iext-1], iext-1, pinti);
3170  pint[iint-1] = pinti;
3171  }
3172 }
3173 
3174 ////////////////////////////////////////////////////////////////////////////////
3175 /// Removes parameter IINT from the internal parameter list
3176 ///
3177 /// and arranges the rest of the list to fill the hole.
3178 
3179 void TMinuit::mnfixp(Int_t iint1, Int_t &ierr)
3180 {
3181  /* Local variables */
3182  Double_t yyover;
3183  Int_t kold, nold, ndex, knew, iext, i, j, m, n, lc, ik;
3184 
3185 // first see if it can be done
3186  ierr = 0;
3187  Int_t iint = iint1+1;
3188  if (iint > fNpar || iint <= 0) {
3189  ierr = 1;
3190  Printf(" MINUIT ERROR. ARGUMENT TO MNFIXP=%4d",iint);
3191  return;
3192  }
3193  iext = fNexofi[iint-1];
3194  if (fNpfix >= fMaxpar) {
3195  ierr = 1;
3196  Printf(" MINUIT CANNOT FIX PARAMETER %4d MAXIMUM NUMBER THAT CAN BE FIXED IS %d",iext,fMaxpar);
3197  return;
3198  }
3199 // reduce number of variable parameters by one
3200 
3201  fNiofex[iext-1] = 0;
3202  nold = fNpar;
3203  --fNpar;
3204 // save values in case parameter is later restored
3205 
3206  ++fNpfix;
3207  fIpfix[fNpfix-1] = iext;
3208  lc = iint;
3209  fXs[fNpfix-1] = fX[lc-1];
3210  fXts[fNpfix-1] = fXt[lc-1];
3211  fDirins[fNpfix-1] = fWerr[lc-1];
3212  fGrds[fNpfix-1] = fGrd[lc-1];
3213  fG2s[fNpfix-1] = fG2[lc-1];
3214  fGsteps[fNpfix-1] = fGstep[lc-1];
3215 // shift values for other parameters to fill hole
3216  for (ik = iext + 1; ik <= fNu; ++ik) {
3217  if (fNiofex[ik-1] > 0) {
3218  lc = fNiofex[ik-1] - 1;
3219  fNiofex[ik-1] = lc;
3220  fNexofi[lc-1] = ik;
3221  fX[lc-1] = fX[lc];
3222  fXt[lc-1] = fXt[lc];
3223  fDirin[lc-1] = fDirin[lc];
3224  fWerr[lc-1] = fWerr[lc];
3225  fGrd[lc-1] = fGrd[lc];
3226  fG2[lc-1] = fG2[lc];
3227  fGstep[lc-1] = fGstep[lc];
3228  }
3229  }
3230  if (fISW[1] <= 0) return;
3231 // remove one row and one column from variance matrix
3232  if (fNpar <= 0) return;
3233  for (i = 1; i <= nold; ++i) {
3234  m = TMath::Max(i,iint);
3235  n = TMath::Min(i,iint);
3236  ndex = m*(m-1) / 2 + n;
3237  fFIXPyy[i-1] = fVhmat[ndex-1];
3238  }
3239  yyover = 1 / fFIXPyy[iint-1];
3240  knew = 0;
3241  kold = 0;
3242  for (i = 1; i <= nold; ++i) {
3243  for (j = 1; j <= i; ++j) {
3244  ++kold;
3245  if (j == iint || i == iint) continue;
3246  ++knew;
3247  fVhmat[knew-1] = fVhmat[kold-1] - fFIXPyy[j-1]*fFIXPyy[i-1]*yyover;
3248  }
3249  }
3250 }
3251 
3252 ////////////////////////////////////////////////////////////////////////////////
3253 /// Restores one or more fixed parameter(s) to variable status
3254 ///
3255 /// Restores one or more fixed parameter(s) to variable status
3256 /// by inserting it into the internal parameter list at the
3257 /// appropriate place.
3258 ///
3259 /// - K = 0 means restore all parameters
3260 /// - K = 1 means restore the last parameter fixed
3261 /// - K = -I means restore external parameter I (if possible)
3262 /// - IQ = fix-location where internal parameters were stored
3263 /// - IR = external number of parameter being restored
3264 /// - IS = internal number of parameter being restored
3265 
3267 {
3268  /* Local variables */
3269  Double_t grdv, xv, dirinv, g2v, gstepv, xtv;
3270  Int_t i, ipsav, ka, lc, ik, iq, ir, is;
3271 
3272  if (k > 1) {
3273  Printf(" CALL TO MNFREE IGNORED. ARGUMENT GREATER THAN ONE");
3274  }
3275  if (fNpfix < 1) {
3276  Printf(" CALL TO MNFREE IGNORED. THERE ARE NO FIXED PARAMETERS");
3277  }
3278  if (k == 1 || k == 0) goto L40;
3279 
3280 // release parameter with specified external number
3281  ka = abs(k);
3282  if (fNiofex[ka-1] == 0) goto L15;
3283  Printf(" IGNORED. PARAMETER SPECIFIED IS ALREADY VARIABLE.");
3284  return;
3285 L15:
3286  if (fNpfix < 1) goto L21;
3287  for (ik = 1; ik <= fNpfix; ++ik) { if (fIpfix[ik-1] == ka) goto L24; }
3288 L21:
3289  Printf(" PARAMETER %4d NOT FIXED. CANNOT BE RELEASED.",ka);
3290  return;
3291 L24:
3292  if (ik == fNpfix) goto L40;
3293 
3294 // move specified parameter to end of list
3295  ipsav = ka;
3296  xv = fXs[ik-1];
3297  xtv = fXts[ik-1];
3298  dirinv = fDirins[ik-1];
3299  grdv = fGrds[ik-1];
3300  g2v = fG2s[ik-1];
3301  gstepv = fGsteps[ik-1];
3302  for (i = ik + 1; i <= fNpfix; ++i) {
3303  fIpfix[i-2] = fIpfix[i-1];
3304  fXs[i-2] = fXs[i-1];
3305  fXts[i-2] = fXts[i-1];
3306  fDirins[i-2] = fDirins[i-1];
3307  fGrds[i-2] = fGrds[i-1];
3308  fG2s[i-2] = fG2s[i-1];
3309  fGsteps[i-2] = fGsteps[i-1];
3310  }
3311  fIpfix[fNpfix-1] = ipsav;
3312  fXs[fNpfix-1] = xv;
3313  fXts[fNpfix-1] = xtv;
3314  fDirins[fNpfix-1] = dirinv;
3315  fGrds[fNpfix-1] = grdv;
3316  fG2s[fNpfix-1] = g2v;
3317  fGsteps[fNpfix-1] = gstepv;
3318 // restore last parameter in fixed list -- IPFIX(NPFIX)
3319 L40:
3320  if (fNpfix < 1) goto L300;
3321  ir = fIpfix[fNpfix-1];
3322  is = 0;
3323  for (ik = fNu; ik >= ir; --ik) {
3324  if (fNiofex[ik-1] > 0) {
3325  lc = fNiofex[ik-1] + 1;
3326  is = lc - 1;
3327  fNiofex[ik-1] = lc;
3328  fNexofi[lc-1] = ik;
3329  fX[lc-1] = fX[lc-2];
3330  fXt[lc-1] = fXt[lc-2];
3331  fDirin[lc-1] = fDirin[lc-2];
3332  fWerr[lc-1] = fWerr[lc-2];
3333  fGrd[lc-1] = fGrd[lc-2];
3334  fG2[lc-1] = fG2[lc-2];
3335  fGstep[lc-1] = fGstep[lc-2];
3336  }
3337  }
3338  ++fNpar;
3339  if (is == 0) is = fNpar;
3340  fNiofex[ir-1] = is;
3341  fNexofi[is-1] = ir;
3342  iq = fNpfix;
3343  fX[is-1] = fXs[iq-1];
3344  fXt[is-1] = fXts[iq-1];
3345  fDirin[is-1] = fDirins[iq-1];
3346  fWerr[is-1] = fDirins[iq-1];
3347  fGrd[is-1] = fGrds[iq-1];
3348  fG2[is-1] = fG2s[iq-1];
3349  fGstep[is-1] = fGsteps[iq-1];
3350  --fNpfix;
3351  fISW[1] = 0;
3352  fDcovar = 1;
3353  if (fISW[4] - fItaur >= 1) {
3354  Printf(" PARAMETER %4d %s RESTORED TO VARIABLE.",ir,
3355  (const char*)fCpnam[ir-1]);
3356  }
3357  if (k == 0) goto L40;
3358 L300:
3359 // if different from internal, external values are taken
3360  mnexin(fX);
3361 }
3362 
3363 ////////////////////////////////////////////////////////////////////////////////
3364 /// Interprets the SET GRAD command
3365 ///
3366 /// - Called from MNSET
3367 /// - Interprets the SET GRAD command, which informs MINUIT whether
3368 /// - the first derivatives of FCN will be calculated by the user
3369 /// - inside FCN. It can check the user derivative calculation
3370 /// - by comparing it with a finite difference approximation.
3371 
3373 {
3374  /* Local variables */
3375  Double_t fzero, err;
3376  Int_t i, nparx, lc, istsav;
3377  Bool_t lnone;
3378 
3379  fISW[2] = 1;
3380  nparx = fNpar;
3381  if (fWord7[0] > 0) goto L2000;
3382 
3383 // get user-calculated first derivatives from FCN
3384  for (i = 1; i <= fNu; ++i) { fGin[i-1] = fUndefi; }
3385  mninex(fX);
3386  Eval(nparx, fGin, fzero, fU, 2); ++fNfcn;
3387  mnderi();
3388  for (i = 1; i <= fNpar; ++i) { fGRADgf[i-1] = fGrd[i-1]; }
3389 // get MINUIT-calculated first derivatives
3390  fISW[2] = 0;
3391  istsav = fIstrat;
3392  fIstrat = 2;
3393  mnhes1();
3394  fIstrat = istsav;
3395  Printf(" CHECK OF GRADIENT CALCULATION IN FCN");
3396  Printf(" PARAMETER G(IN FCN) G(MINUIT) DG(MINUIT) AGREEMENT");
3397  fISW[2] = 1;
3398  lnone = kFALSE;
3399  for (lc = 1; lc <= fNpar; ++lc) {
3400  i = fNexofi[lc-1];
3401  const char *cwd = "GOOD";
3402  err = fDgrd[lc-1];
3403  if (TMath::Abs(fGRADgf[lc-1] - fGrd[lc-1]) > err) {
3404  cwd = " BAD";
3405  fISW[2] = 0;
3406  }
3407  if (fGin[i-1] == fUndefi) {
3408  cwd = "NONE";
3409  lnone = kTRUE;
3410  fGRADgf[lc-1] = 0;
3411  fISW[2] = 0;
3412  }
3413  Printf(" %5d %10s%12.4e%12.4e%12.4e %s",i
3414  ,(const char*)fCpnam[i-1]
3415  ,fGRADgf[lc-1],fGrd[lc-1],err,cwd);
3416  }
3417  if (lnone) {
3418  Printf(" AGREEMENT=NONE MEANS FCN DID NOT CALCULATE THE DERIVATIVE");
3419  }
3420  if (fISW[2] == 0) {
3421  Printf(" MINUIT DOES NOT ACCEPT DERIVATIVE CALCULATIONS BY FCN");
3422  Printf(" TO FORCE ACCEPTANCE, ENTER *SET GRAD 1*");
3423  }
3424 
3425 L2000:
3426  return;
3427 }
3428 
3429 ////////////////////////////////////////////////////////////////////////////////
3430 /// interface to Minuit help
3431 
3432 void TMinuit::mnhelp(const char *command)
3433 {
3434  TString comd = command;
3435  mnhelp(comd);
3436 }
3437 
3438 ////////////////////////////////////////////////////////////////////////////////
3439 /// HELP routine for MINUIT interactive commands
3440 ///
3441 /// - COMD ='*' or "" prints a global help for all commands
3442 /// - COMD =Command_name: print detailed help for one command.
3443 /// Note that at least 3 characters must be given for the command
3444 /// name.
3445 ///
3446 /// Author: Rene Brun
3447 /// comments extracted from the MINUIT documentation file.
3448 
3450 {
3451 //______________________________________________________________________________
3452 //
3453 // Global HELP: Summary of all commands
3454 //
3455  comd.ToUpper();
3456  if( comd.Length() == 0 || comd[0] == '*' || comd[0] == '?' || comd[0] == 0 || comd=="HELP" ) {
3457  Printf(" ==>List of MINUIT Interactive commands:");
3458  Printf(" CLEar Reset all parameter names and values undefined");
3459  Printf(" CONtour Make contour map of the user function");
3460  Printf(" EXIT Exit from Interactive Minuit");
3461  Printf(" FIX Cause parameter(s) to remain constant");
3462  Printf(" HESse Calculate the Hessian or error matrix.");
3463  Printf(" IMPROVE Search for a new minimum around current minimum");
3464  Printf(" MIGrad Minimize by the method of Migrad");
3465  Printf(" MINImize MIGRAD + SIMPLEX method if Migrad fails");
3466  Printf(" MINOs Exact (non-linear) parameter error analysis");
3467  Printf(" MNContour Calculate one MINOS function contour");
3468  Printf(" PARameter Define or redefine new parameters and values");
3469  Printf(" RELease Make previously FIXed parameters variable again");
3470  Printf(" REStore Release last parameter fixed");
3471  Printf(" SAVe Save current parameter values on a file");
3472  Printf(" SCAn Scan the user function by varying parameters");
3473  Printf(" SEEk Minimize by the method of Monte Carlo");
3474  Printf(" SET Set various MINUIT constants or conditions");
3475  Printf(" SHOw Show values of current constants or conditions");
3476  Printf(" SIMplex Minimize by the method of Simplex");
3477  goto L99;
3478  }
3479 
3480 //______________________________________________________________________________
3481 //
3482 // Command CLEAR
3483 //
3484  if( !strncmp(comd.Data(),"CLE",3) ) {
3485  Printf(" ***>CLEAR");
3486  Printf(" Resets all parameter names and values to undefined.");
3487  Printf(" Must normally be followed by a PARameters command or ");
3488  Printf(" equivalent, in order to define parameter values.");
3489  goto L99;
3490  }
3491 //______________________________________________________________________________
3492 //
3493 // Command CONTOUR
3494 //
3495  if( !strncmp(comd.Data(),"CON",3) ) {
3496  Printf(" ***>CONTOUR <par1> <par2> [devs] [ngrid]");
3497  Printf(" Instructs Minuit to trace contour lines of the user function");
3498  Printf(" with respect to the two parameters whose external numbers");
3499  Printf(" are <par1> and <par2>.");
3500  Printf(" Other variable parameters of the function, if any, will have");
3501  Printf(" their values fixed at the current values during the contour");
3502  Printf(" tracing. The optional parameter [devs] (default value 2.)");
3503  Printf(" gives the number of standard deviations in each parameter");
3504  Printf(" which should lie entirely within the plotting area.");
3505  Printf(" Optional parameter [ngrid] (default value 25 unless page");
3506  Printf(" size is too small) determines the resolution of the plot,");
3507  Printf(" i.e. the number of rows and columns of the grid at which the");
3508  Printf(" function will be evaluated. [See also MNContour.]");
3509  goto L99;
3510  }
3511 //______________________________________________________________________________
3512 //
3513 // Command END
3514 //
3515  if( !strncmp(comd.Data(),"END",3) ) {
3516  Printf(" ***>END");
3517  Printf(" Signals the end of a data block (i.e., the end of a fit),");
3518  Printf(" and implies that execution should continue, because another");
3519  Printf(" Data Block follows. A Data Block is a set of Minuit data");
3520  Printf(" consisting of");
3521  Printf(" (1) A Title,");
3522  Printf(" (2) One or more Parameter Definitions,");
3523  Printf(" (3) A blank line, and");
3524  Printf(" (4) A set of Minuit Commands.");
3525  Printf(" The END command is used when more than one Data Block is to");
3526  Printf(" be used with the same FCN function. It first causes Minuit");
3527  Printf(" to issue a CALL FCN with IFLAG=3, in order to allow FCN to");
3528  Printf(" perform any calculations associated with the final fitted");
3529  Printf(" parameter values, unless a CALL FCN 3 command has already");
3530  Printf(" been executed at the current FCN value.");
3531  goto L99;
3532  }
3533 //______________________________________________________________________________
3534 //
3535 // Command EXIT
3536 //
3537  if( !strncmp(comd.Data(),"EXI",3) ) {
3538  Printf(" ***>EXIT");
3539  Printf(" Signals the end of execution.");
3540  Printf(" The EXIT command first causes Minuit to issue a CALL FCN");
3541  Printf(" with IFLAG=3, to allow FCN to perform any calculations");
3542  Printf(" associated with the final fitted parameter values, unless a");
3543  Printf(" CALL FCN 3 command has already been executed.");
3544  goto L99;
3545  }
3546 //______________________________________________________________________________
3547 //
3548 // Command FIX
3549 //
3550  if( !strncmp(comd.Data(),"FIX",3) ) {
3551  Printf(" ***>FIX} <parno> [parno] ... [parno]");
3552  Printf(" Causes parameter(s) <parno> to be removed from the list of");
3553  Printf(" variable parameters, and their value(s) will remain constant");
3554  Printf(" during subsequent minimizations, etc., until another command");
3555  Printf(" changes their value(s) or status.");
3556  goto L99;
3557  }
3558 //______________________________________________________________________________
3559 //
3560 // Command HESSE
3561 //
3562  if( !strncmp(comd.Data(),"HES",3) ) {
3563  Printf(" ***>HESse [maxcalls]");
3564  Printf(" Calculate, by finite differences, the Hessian or error matrix.");
3565  Printf(" That is, it calculates the full matrix of second derivatives");
3566  Printf(" of the function with respect to the currently variable");
3567  Printf(" parameters, and inverts it, printing out the resulting error");
3568  Printf(" matrix. The optional argument [maxcalls] specifies the");
3569  Printf(" (approximate) maximum number of function calls after which");
3570  Printf(" the calculation will be stopped.");
3571  goto L99;
3572  }
3573 //______________________________________________________________________________
3574 //
3575 // Command IMPROVE
3576 //
3577  if( !strncmp(comd.Data(),"IMP",3) ) {
3578  Printf(" ***>IMPROVE [maxcalls]");
3579  Printf(" If a previous minimization has converged, and the current");
3580  Printf(" values of the parameters therefore correspond to a local");
3581  Printf(" minimum of the function, this command requests a search for");
3582  Printf(" additional distinct local minima.");
3583  Printf(" The optional argument [maxcalls] specifies the (approximate");
3584  Printf(" maximum number of function calls after which the calculation");
3585  Printf(" will be stopped.");
3586  goto L99;
3587  }
3588 //______________________________________________________________________________
3589 //
3590 // Command MIGRAD
3591 //
3592  if( !strncmp(comd.Data(),"MIG",3) ) {
3593  Printf(" ***>MIGrad [maxcalls] [tolerance]");
3594  Printf(" Causes minimization of the function by the method of Migrad,");
3595  Printf(" the most efficient and complete single method, recommended");
3596  Printf(" for general functions (see also MINImize).");
3597  Printf(" The minimization produces as a by-product the error matrix");
3598  Printf(" of the parameters, which is usually reliable unless warning");
3599  Printf(" messages are produced.");
3600  Printf(" The optional argument [maxcalls] specifies the (approximate)");
3601  Printf(" maximum number of function calls after which the calculation");
3602  Printf(" will be stopped even if it has not yet converged.");
3603  Printf(" The optional argument [tolerance] specifies required tolerance");
3604  Printf(" on the function value at the minimum.");
3605  Printf(" The default tolerance is 0.1, and the minimization will stop");
3606  Printf(" when the estimated vertical distance to the minimum (EDM) is");
3607  Printf(" less than 0.001*[tolerance]*UP (see [SET ERRordef]).");
3608  goto L99;
3609  }
3610 //______________________________________________________________________________
3611 //
3612 // Command MINIMIZE
3613 //
3614  if( !strncmp(comd.Data(),"MINI",4) ) {
3615  Printf(" ***>MINImize [maxcalls] [tolerance]");
3616  Printf(" Causes minimization of the function by the method of Migrad,");
3617  Printf(" as does the MIGrad command, but switches to the SIMplex method");
3618  Printf(" if Migrad fails to converge. Arguments are as for MIGrad.");
3619  Printf(" Note that command requires four characters to be unambiguous.");
3620  goto L99;
3621  }
3622 //______________________________________________________________________________
3623 //
3624 // Command MINOS
3625 //
3626  if( !strncmp(comd.Data(),"MIN0",4) ) {
3627  Printf(" ***>MINOs [maxcalls] [parno] [parno] ...");
3628  Printf(" Causes a Minos error analysis to be performed on the parameters");
3629  Printf(" whose numbers [parno] are specified. If none are specified,");
3630  Printf(" Minos errors are calculated for all variable parameters.");
3631  Printf(" Minos errors may be expensive to calculate, but are very");
3632  Printf(" reliable since they take account of non-linearities in the");
3633  Printf(" problem as well as parameter correlations, and are in general");
3634  Printf(" asymmetric.");
3635  Printf(" The optional argument [maxcalls] specifies the (approximate)");
3636  Printf(" maximum number of function calls per parameter requested,");
3637  Printf(" after which the calculation will stop for that parameter.");
3638  goto L99;
3639  }
3640 //______________________________________________________________________________
3641 //
3642 // Command MNCONTOUR
3643 //
3644  if( !strncmp(comd.Data(),"MNC",3) ) {
3645  Printf(" ***>MNContour <par1> <par2> [npts]");
3646  Printf(" Calculates one function contour of FCN with respect to");
3647  Printf(" parameters par1 and par2, with FCN minimized always with");
3648  Printf(" respect to all other NPAR-2 variable parameters (if any).");
3649  Printf(" Minuit will try to find npts points on the contour (default 20)");
3650  Printf(" If only two parameters are variable at the time, it is not");
3651  Printf(" necessary to specify their numbers. To calculate more than");
3652  Printf(" one contour, it is necessary to SET ERRordef to the appropriate");
3653  Printf(" value and issue the MNContour command for each contour.");
3654  goto L99;
3655  }
3656 //______________________________________________________________________________
3657 //
3658 // Command PARAMETER
3659 //
3660  if( !strncmp(comd.Data(),"PAR",3) ) {
3661  Printf(" ***>PARameters");
3662  Printf(" followed by one or more parameter definitions.");
3663  Printf(" Parameter definitions are of the form:");
3664  Printf(" <number> ''name'' <value> <step> [lolim] [uplim] ");
3665  Printf(" for example:");
3666  Printf(" 3 ''K width'' 1.2 0.1");
3667  Printf(" the last definition is followed by a blank line or a zero.");
3668  goto L99;
3669  }
3670 //______________________________________________________________________________
3671 //
3672 // Command RELEASE
3673 //
3674  if( !strncmp(comd.Data(),"REL",3) ) {
3675  Printf(" ***>RELease <parno> [parno] ... [parno]");
3676  Printf(" If <parno> is the number of a previously variable parameter");
3677  Printf(" which has been fixed by a command: FIX <parno>, then that");
3678  Printf(" parameter will return to variable status. Otherwise a warning");
3679  Printf(" message is printed and the command is ignored.");
3680  Printf(" Note that this command operates only on parameters which were");
3681  Printf(" at one time variable and have been FIXed. It cannot make");
3682  Printf(" constant parameters variable; that must be done by redefining");
3683  Printf(" the parameter with a PARameters command.");
3684  goto L99;
3685  }
3686 //______________________________________________________________________________
3687 //
3688 // Command RESTORE
3689 //
3690  if( !strncmp(comd.Data(),"RES",3) ) {
3691  Printf(" ***>REStore [code]");
3692  Printf(" If no [code] is specified, this command restores all previously");
3693  Printf(" FIXed parameters to variable status. If [code]=1, then only");
3694  Printf(" the last parameter FIXed is restored to variable status.");
3695  Printf(" If code is neither zero nor one, the command is ignored.");
3696  goto L99;
3697  }
3698 //______________________________________________________________________________
3699 //
3700 // Command RETURN
3701 //
3702  if( !strncmp(comd.Data(),"RET",3) ) {
3703  Printf(" ***>RETURN");
3704  Printf(" Signals the end of a data block, and instructs Minuit to return");
3705  Printf(" to the program which called it. The RETurn command first");
3706  Printf(" causes Minuit to CALL FCN with IFLAG=3, in order to allow FCN");
3707  Printf(" to perform any calculations associated with the final fitted");
3708  Printf(" parameter values, unless a CALL FCN 3 command has already been");
3709  Printf(" executed at the current FCN value.");
3710  goto L99;
3711  }
3712 //______________________________________________________________________________
3713 //
3714 // Command SAVE
3715 //
3716  if( !strncmp(comd.Data(),"SAV",3) ) {
3717  Printf(" ***>SAVe");
3718  Printf(" Causes the current parameter values to be saved on a file in");
3719  Printf(" such a format that they can be read in again as Minuit");
3720  Printf(" parameter definitions. If the covariance matrix exists, it is");
3721  Printf(" also output in such a format. The unit number is by default 7,");
3722  Printf(" or that specified by the user in their call to MINTIO or");
3723  Printf(" MNINIT. The user is responsible for opening the file previous");
3724  Printf(" to issuing the [SAVe] command (except where this can be done");
3725  Printf(" interactively).");
3726  goto L99;
3727  }
3728 //______________________________________________________________________________
3729 //
3730 // Command SCAN
3731 //
3732  if( !strncmp(comd.Data(),"SCA",3) ) {
3733  Printf(" ***>SCAn [parno] [numpts] [from] [to]");
3734  Printf(" Scans the value of the user function by varying parameter");
3735  Printf(" number [parno], leaving all other parameters fixed at the");
3736  Printf(" current value. If [parno] is not specified, all variable");
3737  Printf(" parameters are scanned in sequence.");
3738  Printf(" The number of points [numpts] in the scan is 40 by default,");
3739  Printf(" and cannot exceed 100. The range of the scan is by default");
3740  Printf(" 2 standard deviations on each side of the current best value,");
3741  Printf(" but can be specified as from [from] to [to].");
3742  Printf(" After each scan, if a new minimum is found, the best parameter");
3743  Printf(" values are retained as start values for future scans or");
3744  Printf(" minimizations. The curve resulting from each scan is plotted");
3745  Printf(" on the output unit in order to show the approximate behaviour");
3746  Printf(" of the function.");
3747  Printf(" This command is not intended for minimization, but is sometimes");
3748  Printf(" useful for debugging the user function or finding a");
3749  Printf(" reasonable starting point.");
3750  goto L99;
3751  }
3752 //______________________________________________________________________________
3753 //
3754 // Command SEEK
3755 //
3756  if( !strncmp(comd.Data(),"SEE",3) ) {
3757  Printf(" ***>SEEk [maxcalls] [devs]");
3758  Printf(" Causes a Monte Carlo minimization of the function, by choosing");
3759  Printf(" random values of the variable parameters, chosen uniformly");
3760  Printf(" over a hypercube centered at the current best value.");
3761  Printf(" The region size is by default 3 standard deviations on each");
3762  Printf(" side, but can be changed by specifying the value of [devs].");
3763  goto L99;
3764  }
3765 //______________________________________________________________________________
3766 //
3767 // Command SET
3768 //
3769  if( !strncmp(comd.Data(),"SET",3) ) {
3770  Printf(" ***>SET <option_name>");
3771  Printf(" SET BATch");
3772  Printf(" Informs Minuit that it is running in batch mode.");
3773 
3774  Printf(" ");
3775  Printf(" SET EPSmachine <accuracy>");
3776  Printf(" Informs Minuit that the relative floating point arithmetic");
3777  Printf(" precision is <accuracy>. Minuit determines the nominal");
3778  Printf(" precision itself, but the SET EPSmachine command can be");
3779  Printf(" used to override Minuit own determination, when the user");
3780  Printf(" knows that the FCN function value is not calculated to");
3781  Printf(" the nominal machine accuracy. Typical values of <accuracy>");
3782  Printf(" are between 10**-5 and 10**-14.");
3783 
3784  Printf(" ");
3785  Printf(" SET ERRordef <up>");
3786  Printf(" Sets the value of UP (default value= 1.), defining parameter");
3787  Printf(" errors. Minuit defines parameter errors as the change");
3788  Printf(" in parameter value required to change the function value");
3789  Printf(" by UP. Normally, for chisquared fits UP=1, and for negative");
3790  Printf(" log likelihood, UP=0.5.");
3791 
3792  Printf(" ");
3793  Printf(" SET GRAdient [force]");
3794  Printf(" Informs Minuit that the user function is prepared to");
3795  Printf(" calculate its own first derivatives and return their values");
3796  Printf(" in the array GRAD when IFLAG=2 (see specs of FCN).");
3797  Printf(" If [force] is not specified, Minuit will calculate");
3798  Printf(" the FCN derivatives by finite differences at the current");
3799  Printf(" point and compare with the user calculation at that point,");
3800  Printf(" accepting the user values only if they agree.");
3801  Printf(" If [force]=1, Minuit does not do its own derivative");
3802  Printf(" calculation, and uses the derivatives calculated in FCN.");
3803 
3804  Printf(" ");
3805  Printf(" SET INPut [unitno] [filename]");
3806  Printf(" Causes Minuit, in data-driven mode only, to read subsequent");
3807  Printf(" commands (or parameter definitions) from a different input");
3808  Printf(" file. If no [unitno] is specified, reading reverts to the");
3809  Printf(" previous input file, assuming that there was one.");
3810  Printf(" If [unitno] is specified, and that unit has not been opened,");
3811  Printf(" then Minuit attempts to open the file [filename]} if a");
3812  Printf(" name is specified. If running in interactive mode and");
3813  Printf(" [filename] is not specified and [unitno] is not opened,");
3814  Printf(" Minuit prompts the user to enter a file name.");
3815  Printf(" If the word REWIND is added to the command (note:no blanks");
3816  Printf(" between INPUT and REWIND), the file is rewound before");
3817  Printf(" reading. Note that this command is implemented in standard");
3818  Printf(" Fortran 77 and the results may depend on the system;");
3819  Printf(" for example, if a filename is given under VM/CMS, it must");
3820  Printf(" be preceded by a slash.");
3821 
3822  Printf(" ");
3823  Printf(" SET INTeractive");
3824  Printf(" Informs Minuit that it is running interactively.");
3825 
3826  Printf(" ");
3827  Printf(" SET LIMits [parno] [lolim] [uplim]");
3828  Printf(" Allows the user to change the limits on one or all");
3829  Printf(" parameters. If no arguments are specified, all limits are");
3830  Printf(" removed from all parameters. If [parno] alone is specified,");
3831  Printf(" limits are removed from parameter [parno].");
3832  Printf(" If all arguments are specified, then parameter [parno] will");
3833  Printf(" be bounded between [lolim] and [uplim].");
3834  Printf(" Limits can be specified in either order, Minuit will take");
3835  Printf(" the smaller as [lolim] and the larger as [uplim].");
3836  Printf(" However, if [lolim] is equal to [uplim], an error condition");
3837  Printf(" results.");
3838 
3839  Printf(" ");
3840  Printf(" SET LINesperpage");
3841  Printf(" Sets the number of lines for one page of output.");
3842  Printf(" Default value is 24 for interactive mode");
3843 
3844  Printf(" ");
3845  Printf(" SET NOGradient");
3846  Printf(" The inverse of SET GRAdient, instructs Minuit not to");
3847  Printf(" use the first derivatives calculated by the user in FCN.");
3848 
3849  Printf(" ");
3850  Printf(" SET NOWarnings");
3851  Printf(" Supresses Minuit warning messages.");
3852 
3853  Printf(" ");
3854  Printf(" SET OUTputfile <unitno>");
3855  Printf(" Instructs Minuit to write further output to unit <unitno>.");
3856 
3857  Printf(" ");
3858  Printf(" SET PAGethrow <integer>");
3859  Printf(" Sets the carriage control character for ``new page'' to");
3860  Printf(" <integer>. Thus the value 1 produces a new page, and 0");
3861  Printf(" produces a blank line, on some devices (see TOPofpage)");
3862 
3863 
3864  Printf(" ");
3865  Printf(" SET PARameter <parno> <value>");
3866  Printf(" Sets the value of parameter <parno> to <value>.");
3867  Printf(" The parameter in question may be variable, fixed, or");
3868  Printf(" constant, but must be defined.");
3869 
3870  Printf(" ");
3871  Printf(" SET PRIntout <level>");
3872  Printf(" Sets the print level, determining how much output will be");
3873  Printf(" produced. Allowed values and their meanings are displayed");
3874  Printf(" after a SHOw PRInt command, and are currently <level>=:");
3875  Printf(" [-1] no output except from SHOW commands");
3876  Printf(" [0] minimum output");
3877  Printf(" [1] default value, normal output");
3878  Printf(" [2] additional output giving intermediate results.");
3879  Printf(" [3] maximum output, showing progress of minimizations.");
3880  Printf(" Note: See also the SET WARnings command.");
3881 
3882  Printf(" ");
3883  Printf(" SET RANdomgenerator <seed>");
3884  Printf(" Sets the seed of the random number generator used in SEEk.");
3885  Printf(" This can be any integer between 10000 and 900000000, for");
3886  Printf(" example one which was output from a SHOw RANdom command of");
3887  Printf(" a previous run.");
3888 
3889  Printf(" ");
3890  Printf(" SET STRategy <level>");
3891  Printf(" Sets the strategy to be used in calculating first and second");
3892  Printf(" derivatives and in certain minimization methods.");
3893  Printf(" In general, low values of <level> mean fewer function calls");
3894  Printf(" and high values mean more reliable minimization.");
3895  Printf(" Currently allowed values are 0, 1 (default), and 2.");
3896 
3897  Printf(" ");
3898  Printf(" SET TITle");
3899  Printf(" Informs Minuit that the next input line is to be considered");
3900  Printf(" the (new) title for this task or sub-task. This is for");
3901  Printf(" the convenience of the user in reading their output.");
3902 
3903  Printf(" ");
3904  Printf(" SET WARnings");
3905  Printf(" Instructs Minuit to output warning messages when suspicious");
3906  Printf(" conditions arise which may indicate unreliable results.");
3907  Printf(" This is the default.");
3908 
3909  Printf(" ");
3910  Printf(" SET WIDthpage");
3911  Printf(" Informs Minuit of the output page width.");
3912  Printf(" Default values are 80 for interactive jobs");
3913  goto L99;
3914  }
3915 //______________________________________________________________________________
3916 //
3917 // Command SHOW
3918 //
3919  if( !strncmp(comd.Data(),"SHO",3) ) {
3920  Printf(" ***>SHOw <option_name>");
3921  Printf(" All SET XXXX commands have a corresponding SHOw XXXX command.");
3922  Printf(" In addition, the SHOw commands listed starting here have no");
3923  Printf(" corresponding SET command for obvious reasons.");
3924 
3925  Printf(" ");
3926  Printf(" SHOw CORrelations");
3927  Printf(" Calculates and prints the parameter correlations from the");
3928  Printf(" error matrix.");
3929 
3930  Printf(" ");
3931  Printf(" SHOw COVariance");
3932  Printf(" Prints the (external) covariance (error) matrix.");
3933 
3934  Printf(" ");
3935  Printf(" SHOw EIGenvalues");
3936  Printf(" Calculates and prints the eigenvalues of the covariance");
3937  Printf(" matrix.");
3938 
3939  Printf(" ");
3940  Printf(" SHOw FCNvalue");
3941  Printf(" Prints the current value of FCN.");
3942  goto L99;
3943  }
3944 //______________________________________________________________________________
3945 //
3946 // Command SIMPLEX
3947 //
3948  if( !strncmp(comd.Data(),"SIM",3) ) {
3949  Printf(" ***>SIMplex [maxcalls] [tolerance]");
3950  Printf(" Performs a function minimization using the simplex method of");
3951  Printf(" Nelder and Mead. Minimization terminates either when the");
3952  Printf(" function has been called (approximately) [maxcalls] times,");
3953  Printf(" or when the estimated vertical distance to minimum (EDM) is");
3954  Printf(" less than [tolerance].");
3955  Printf(" The default value of [tolerance] is 0.1*UP(see SET ERRordef).");
3956  goto L99;
3957  }
3958 //______________________________________________________________________________
3959 //
3960 // Command STANDARD
3961 //
3962  if( !strncmp(comd.Data(),"STA",3) ) {
3963  Printf(" ***>STAndard");
3964  goto L99;
3965  }
3966 //______________________________________________________________________________
3967 //
3968 // Command STOP
3969 //
3970  if( !strncmp(comd.Data(),"STO",3) ) {
3971  Printf(" ***>STOP");
3972  Printf(" Same as EXIT.");
3973  goto L99;
3974  }
3975 //______________________________________________________________________________
3976 //
3977 // Command TOPOFPAGE
3978 //
3979  if( !strncmp(comd.Data(),"TOP",3) ) {
3980  Printf(" ***>TOPofpage");
3981  Printf(" Causes Minuit to write the character specified in a");
3982  Printf(" SET PAGethrow command (default = 1) to column 1 of the output");
3983  Printf(" file, which may or may not position your output medium to");
3984  Printf(" the top of a page depending on the device and system.");
3985  goto L99;
3986  }
3987 //______________________________________________________________________________
3988  Printf(" Unknown MINUIT command. Type HELP for list of commands.");
3989 
3990 L99:
3991  return;
3992 }
3993 
3994 ////////////////////////////////////////////////////////////////////////////////
3995 /// Calculates the full second-derivative matrix of FCN
3996 ///
3997 /// by taking finite differences. When calculating diagonal
3998 /// elements, it may iterate so that step size is nearly that
3999 /// which gives function change= UP/10. The first derivatives
4000 /// of course come as a free side effect, but with a smaller
4001 /// step size in order to obtain a known accuracy.
4002 
4004 {
4005  /* Local variables */
4006  Double_t dmin_, dxdi, elem, wint, tlrg2, d, dlast, ztemp, g2bfor;
4007  Double_t df, aimsag, fs1, tlrstp, fs2, stpinm, g2i, sag=0, xtf, xti, xtj;
4008  Int_t icyc, ncyc, ndex, idrv, iext, npar2, i, j, ifail, npard, nparx, id, multpy;
4009  Bool_t ldebug;
4010 
4011  ldebug = fIdbg[3] >= 1;
4012  if (fAmin == fUndefi) {
4013  mnamin();
4014  }
4015  if (fIstrat <= 0) {
4016  ncyc = 3;
4017  tlrstp = .5;
4018  tlrg2 = .1;
4019  } else if (fIstrat == 1) {
4020  ncyc = 5;
4021  tlrstp = .3;
4022  tlrg2 = .05;
4023  } else {
4024  ncyc = 7;
4025  tlrstp = .1;
4026  tlrg2 = .02;
4027  }
4028  if (fISW[4] >= 2 || ldebug) {
4029  Printf(" START COVARIANCE MATRIX CALCULATION.");
4030  }
4031  fCfrom = "HESSE ";
4032  fNfcnfr = fNfcn;
4033  fCstatu = "OK ";
4034  npard = fNpar;
4035 // make sure starting at the right place
4036  mninex(fX);
4037  nparx = fNpar;
4038  Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
4039  if (fs1 != fAmin) {
4040  df = fAmin - fs1;
4041  mnwarn("D", "MNHESS", TString::Format("function value differs from AMIN by %g",df));
4042  }
4043  fAmin = fs1;
4044  if (ldebug) {
4045  Printf(" PAR D GSTEP D G2 GRD SAG ");
4046  }
4047 // diagonal elements .
4048 
4049 // fISW[1] = 1 if approx, 2 if not posdef, 3 if ok
4050 // AIMSAG is the sagitta we are aiming for in second deriv calc.
4051 
4052  aimsag = TMath::Sqrt(fEpsma2)*(TMath::Abs(fAmin) + fUp);
4053 // Zero the second derivative matrix
4054  npar2 = fNpar*(fNpar + 1) / 2;
4055  for (i = 1; i <= npar2; ++i) { fVhmat[i-1] = 0; }
4056 
4057 // Loop over variable parameters for second derivatives
4058  idrv = 2;
4059  for (id = 1; id <= npard; ++id) {
4060  i = id + fNpar - npard;
4061  iext = fNexofi[i-1];
4062  if (fG2[i-1] == 0) {
4063  mnwarn("W", "HESSE", Form("Second derivative enters zero, param %d",iext));
4064  wint = fWerr[i-1];
4065  if (fNvarl[iext-1] > 1) {
4066  mndxdi(fX[i-1], i-1, dxdi);
4067  if (TMath::Abs(dxdi) < .001) wint = .01;
4068  else wint /= TMath::Abs(dxdi);
4069  }
4070  fG2[i-1] = fUp / (wint*wint);
4071  }
4072  xtf = fX[i-1];
4073  dmin_ = fEpsma2*8*TMath::Abs(xtf);
4074 
4075 // find step which gives sagitta = AIMSAG
4076  d = TMath::Abs(fGstep[i-1]);
4077  int skip50 = 0;
4078  for (icyc = 1; icyc <= ncyc; ++icyc) {
4079 // loop here only if SAG=0
4080  for (multpy = 1; multpy <= 5; ++multpy) {
4081 // take two steps
4082  fX[i-1] = xtf + d;
4083  mninex(fX);
4084  nparx = fNpar;
4085  Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
4086  fX[i-1] = xtf - d;
4087  mninex(fX);
4088  Eval(nparx, fGin, fs2, fU, 4); ++fNfcn;
4089  fX[i-1] = xtf;
4090  sag = (fs1 + fs2 - fAmin*2)*.5;
4091  if (sag != 0) goto L30;
4092  if (fGstep[i-1] < 0) {
4093  if (d >= .5) goto L26;
4094  d *= 10;
4095  if (d > .5) d = .51;
4096  continue;
4097  }
4098  d *= 10;
4099  }
4100 L26:
4101  mnwarn("W", "HESSE", TString::Format("Second derivative zero for parameter%d",iext));
4102  goto L390;
4103 // SAG is not zero
4104 L30:
4105  g2bfor = fG2[i-1];
4106  fG2[i-1] = sag*2 / (d*d);
4107  fGrd[i-1] = (fs1 - fs2) / (d*2);
4108  if (ldebug) {
4109  Printf("%4d%2d%12.5g%12.5g%12.5g%12.5g%12.5g",i,idrv,fGstep[i-1],d,fG2[i-1],fGrd[i-1],sag);
4110  }
4111  if (fGstep[i-1] > 0) fGstep[i-1] = TMath::Abs(d);
4112  else fGstep[i-1] = -TMath::Abs(d);
4113  fDirin[i-1] = d;
4114  fHESSyy[i-1]= fs1;
4115  dlast = d;
4116  d = TMath::Sqrt(aimsag*2 / TMath::Abs(fG2[i-1]));
4117 // if parameter has limits, max int step size = 0.5
4118  stpinm = .5;
4119  if (fGstep[i-1] < 0) d = TMath::Min(d,stpinm);
4120  if (d < dmin_) d = dmin_;
4121 // see if converged
4122  if (TMath::Abs((d - dlast) / d) < tlrstp ||
4123  TMath::Abs((fG2[i-1] - g2bfor) / fG2[i-1]) < tlrg2) {
4124  skip50 = 1;
4125  break;
4126  }
4127  d = TMath::Min(d,dlast*102);
4128  d = TMath::Max(d,dlast*.1);
4129  }
4130 // end of step size loop
4131  if (!skip50)
4132  mnwarn("D", "MNHESS", TString::Format("Second Deriv. SAG,AIM= %d%g%g",iext,sag,aimsag));
4133 
4134  ndex = i*(i + 1) / 2;
4135  fVhmat[ndex-1] = fG2[i-1];
4136  }
4137 // end of diagonal second derivative loop
4138  mninex(fX);
4139 // refine the first derivatives
4140  if (fIstrat > 0) mnhes1();
4141  fISW[1] = 3;
4142  fDcovar = 0;
4143 // off-diagonal elements
4144 
4145  if (fNpar == 1) goto L214;
4146  for (i = 1; i <= fNpar; ++i) {
4147  for (j = 1; j <= i-1; ++j) {
4148  xti = fX[i-1];
4149  xtj = fX[j-1];
4150  fX[i-1] = xti + fDirin[i-1];
4151  fX[j-1] = xtj + fDirin[j-1];
4152  mninex(fX);
4153  Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
4154  fX[i-1] = xti;
4155  fX[j-1] = xtj;
4156  elem = (fs1 + fAmin - fHESSyy[i-1] - fHESSyy[j-1]) / (
4157  fDirin[i-1]*fDirin[j-1]);
4158  ndex = i*(i-1) / 2 + j;
4159  fVhmat[ndex-1] = elem;
4160  }
4161  }
4162 L214:
4163  mninex(fX);
4164 // verify matrix positive-definite
4165  mnpsdf();
4166  for (i = 1; i <= fNpar; ++i) {
4167  for (j = 1; j <= i; ++j) {
4168  ndex = i*(i-1) / 2 + j;
4169  fP[i + j*fMaxpar - fMaxpar-1] = fVhmat[ndex-1];
4170  fP[j + i*fMaxpar - fMaxpar-1] = fP[i + j*fMaxpar - fMaxpar-1];
4171  }
4172  }
4173  mnvert(fP, fMaxint, fMaxint, fNpar, ifail);
4174  if (ifail > 0) {
4175  mnwarn("W", "HESSE", "Matrix inversion fails.");
4176  goto L390;
4177  }
4178 // calculate e d m
4179  fEDM = 0;
4180 
4181  for (i = 1; i <= fNpar; ++i) {
4182 // off-diagonal elements
4183  ndex = i*(i-1) / 2;
4184  for (j = 1; j <= i-1; ++j) {
4185  ++ndex;
4186  ztemp = fP[i + j*fMaxpar - fMaxpar-1]*2;
4187  fEDM += fGrd[i-1]*ztemp*fGrd[j-1];
4188  fVhmat[ndex-1] = ztemp;
4189  }
4190 // diagonal elements
4191  ++ndex;
4192  fVhmat[ndex-1] = fP[i + i*fMaxpar - fMaxpar-1]*2;
4193  fEDM += fP[i + i*fMaxpar - fMaxpar-1]*(fGrd[i-1]*fGrd[i-1]);
4194  }
4195  if (fISW[4] >= 1 && fISW[1] == 3 && fItaur == 0) {
4196  Printf(" COVARIANCE MATRIX CALCULATED SUCCESSFULLY");
4197  }
4198  goto L900;
4199 // failure to invert 2nd deriv matrix
4200 L390:
4201  fISW[1] = 1;
4202  fDcovar = 1;
4203  fCstatu = "FAILED ";
4204  if (fISW[4] >= 0) {
4205  Printf(" MNHESS FAILS AND WILL RETURN DIAGONAL MATRIX. ");
4206  }
4207  for (i = 1; i <= fNpar; ++i) {
4208  ndex = i*(i-1) / 2;
4209  for (j = 1; j <= i-1; ++j) {
4210  ++ndex;
4211  fVhmat[ndex-1] = 0;
4212  }
4213  ++ndex;
4214  g2i = fG2[i-1];
4215  if (g2i <= 0) g2i = 1;
4216  fVhmat[ndex-1] = 2 / g2i;
4217  }
4218 L900:
4219  return;
4220 }
4221 
4222 ////////////////////////////////////////////////////////////////////////////////
4223 /// Calculate first derivatives (GRD) and uncertainties (DGRD)
4224 ///
4225 /// and appropriate step sizes GSTEP
4226 /// Called from MNHESS and MNGRAD
4227 
4229 {
4230  /* Local variables */
4231  Double_t dmin_, d, dfmin, dgmin=0, change, chgold, grdold=0, epspri;
4232  Double_t fs1, optstp, fs2, grdnew=0, sag, xtf;
4233  Int_t icyc, ncyc=0, idrv, i, nparx;
4234  Bool_t ldebug;
4235 
4236  ldebug = fIdbg[5] >= 1;
4237  if (fIstrat <= 0) ncyc = 1;
4238  if (fIstrat == 1) ncyc = 2;
4239  if (fIstrat > 1) ncyc = 6;
4240  idrv = 1;
4241  nparx = fNpar;
4242  dfmin = fEpsma2*4*(TMath::Abs(fAmin) + fUp);
4243 // main loop over parameters
4244  for (i = 1; i <= fNpar; ++i) {
4245  xtf = fX[i-1];
4246  dmin_ = fEpsma2*4*TMath::Abs(xtf);
4247  epspri = fEpsma2 + TMath::Abs(fGrd[i-1]*fEpsma2);
4248  optstp = TMath::Sqrt(dfmin / (TMath::Abs(fG2[i-1]) + epspri));
4249  d = TMath::Abs(fGstep[i-1])*.2;
4250  if (d > optstp) d = optstp;
4251  if (d < dmin_) d = dmin_;
4252  chgold = 1e4;
4253 // iterate reducing step size
4254  for (icyc = 1; icyc <= ncyc; ++icyc) {
4255  fX[i-1] = xtf + d;
4256  mninex(fX);
4257  Eval(nparx, fGin, fs1, fU, 4); ++fNfcn;
4258  fX[i-1] = xtf - d;
4259  mninex(fX);
4260  Eval(nparx, fGin, fs2, fU, 4); ++fNfcn;
4261  fX[i-1] = xtf;
4262 // check if step sizes appropriate
4263  sag = (fs1 + fs2 - fAmin*2)*.5;
4264  grdold = fGrd[i-1];
4265  grdnew = (fs1 - fs2) / (d*2);
4266  dgmin = fEpsmac*(TMath::Abs(fs1) + TMath::Abs(fs2)) / d;
4267  if (ldebug) {
4268  Printf("%4d%2d%12.5g%12.5g%12.5g%12.5g%12.5g",i,idrv,fGstep[i-1],d,fG2[i-1],grdnew,sag);
4269  }
4270  if (grdnew == 0) goto L60;
4271  change = TMath::Abs((grdold - grdnew) / grdnew);
4272  if (change > chgold && icyc > 1) goto L60;
4273  chgold = change;
4274  fGrd[i-1] = grdnew;
4275  if (fGstep[i-1] > 0) fGstep[i-1] = TMath::Abs(d);
4276  else fGstep[i-1] = -TMath::Abs(d);
4277 // decrease step until first derivative changes by <5%
4278  if (change < .05) goto L60;
4279  if (TMath::Abs(grdold - grdnew) < dgmin) goto L60;
4280  if (d < dmin_) {
4281  mnwarn("D", "MNHES1", "Step size too small for 1st drv.");
4282  goto L60;
4283  }
4284  d *= .2;
4285  }
4286 // loop satisfied = too many iter
4287  mnwarn("D", "MNHES1", TString::Format("Too many iterations on D1.%g%g",grdold,grdnew));
4288 L60:
4289  fDgrd[i-1] = TMath::Max(dgmin,TMath::Abs(grdold - grdnew));
4290  }
4291 // end of first deriv. loop
4292  mninex(fX);
4293 }
4294 
4295 ////////////////////////////////////////////////////////////////////////////////
4296 /// Attempts to improve on a good local minimum
4297 ///
4298 /// Attempts to improve on a good local minimum by finding a
4299 /// better one. The quadratic part of FCN is removed by MNCALF
4300 /// and this transformed function is minimised using the simplex
4301 /// method from several random starting points.
4302 ///
4303 /// ref. -- Goldstein and Price, Math.Comp. 25, 569 (1971)
4304 
4306 {
4307  /* Initialized data */
4308 
4309  Double_t rnum = 0;
4310 
4311  /* Local variables */
4312  Double_t amax, ycalf, ystar, ystst;
4313  Double_t pb, ep, wg, xi, sigsav, reg, sig2;
4314  Int_t npfn, ndex, loop=0, i, j, ifail, iseed=0;
4315  Int_t jhold, nloop, nparx, nparp1, jh, jl, iswtr;
4316 
4317  if (fNpar <= 0) return;
4318  if (fAmin == fUndefi) mnamin();
4319  fCstatu = "UNCHANGED ";
4320  fItaur = 1;
4321  fEpsi = fUp*.1;
4322  npfn = fNfcn;
4323  nloop = Int_t(fWord7[1]);
4324  if (nloop <= 0) nloop = fNpar + 4;
4325  nparx = fNpar;
4326  nparp1 = fNpar + 1;
4327  wg = 1 / Double_t(fNpar);
4328  sigsav = fEDM;
4329  fApsi = fAmin;
4330  iswtr = fISW[4] - 2*fItaur;
4331  for (i = 1; i <= fNpar; ++i) {
4332  fXt[i-1] = fX[i-1];
4333  fIMPRdsav[i-1] = fWerr[i-1];
4334  for (j = 1; j <= i; ++j) {
4335  ndex = i*(i-1) / 2 + j;
4336  fP[i + j*fMaxpar - fMaxpar-1] = fVhmat[ndex-1];
4337  fP[j + i*fMaxpar - fMaxpar-1] = fP[i + j*fMaxpar - fMaxpar-1];
4338  }
4339  }
4340  mnvert(fP, fMaxint, fMaxint, fNpar, ifail);
4341  if (ifail >= 1) goto L280;
4342 // Save inverted matrix in VT
4343  for (i = 1; i <= fNpar; ++i) {
4344  ndex = i*(i-1) / 2;
4345  for (j = 1; j <= i; ++j) {
4346  ++ndex;
4347  fVthmat[ndex-1] = fP[i + j*fMaxpar - fMaxpar-1];
4348  }
4349  }
4350  loop = 0;
4351 
4352 L20:
4353  for (i = 1; i <= fNpar; ++i) {
4354  fDirin[i-1] = fIMPRdsav[i-1]*2;
4355  mnrn15(rnum, iseed);
4356  fX[i-1] = fXt[i-1] + fDirin[i-1]*2*(rnum - .5);
4357  }
4358  ++loop;
4359  reg = 2;
4360  if (fISW[4] >= 0) {
4361  Printf("START ATTEMPT NO.%2d TO FIND NEW MINIMUM",loop);
4362  }
4363 L30:
4364  mncalf(fX, ycalf);
4365  fAmin = ycalf;
4366 // set up random simplex
4367  jl = nparp1;
4368  jh = nparp1;
4369  fIMPRy[nparp1-1] = fAmin;
4370  amax = fAmin;
4371  for (i = 1; i <= fNpar; ++i) {
4372  xi = fX[i-1];
4373  mnrn15(rnum, iseed);
4374  fX[i-1] = xi - fDirin[i-1]*(rnum - .5);
4375  mncalf(fX, ycalf);
4376  fIMPRy[i-1] = ycalf;
4377  if (fIMPRy[i-1] < fAmin) {
4378  fAmin = fIMPRy[i-1];
4379  jl = i;
4380  } else if (fIMPRy[i-1] > amax) {
4381  amax = fIMPRy[i-1];
4382  jh = i;
4383  }
4384  for (j = 1; j <= fNpar; ++j) { fP[j + i*fMaxpar - fMaxpar-1] = fX[j-1]; }
4385  fP[i + nparp1*fMaxpar - fMaxpar-1] = xi;
4386  fX[i-1] = xi;
4387  }
4388 
4389  fEDM = fAmin;
4390  sig2 = fEDM;
4391 // start main loop
4392 L50:
4393  if (fAmin < 0) goto L95;
4394  if (fISW[1] <= 2) goto L280;
4395  ep = fAmin*.1;
4396  if (sig2 < ep && fEDM < ep) goto L100;
4397  sig2 = fEDM;
4398  if (fNfcn - npfn > fNfcnmx) goto L300;
4399 // calculate new point * by reflection
4400  for (i = 1; i <= fNpar; ++i) {
4401  pb = 0;
4402  for (j = 1; j <= nparp1; ++j) { pb += wg*fP[i + j*fMaxpar - fMaxpar-1]; }
4403  fPbar[i-1] = pb - wg*fP[i + jh*fMaxpar - fMaxpar-1];
4404  fPstar[i-1] = fPbar[i-1]*2 - fP[i + jh*fMaxpar - fMaxpar-1]*1;
4405  }
4406  mncalf(fPstar, ycalf);
4407  ystar = ycalf;
4408  if (ystar >= fAmin) goto L70;
4409 // point * better than jl, calculate new point **
4410  for (i = 1; i <= fNpar; ++i) {
4411  fPstst[i-1] = fPstar[i-1]*2 + fPbar[i- 1]*-1;
4412  }
4413  mncalf(fPstst, ycalf);
4414  ystst = ycalf;
4415  if (ystst < fIMPRy[jl-1]) goto L67;
4416  mnrazz(ystar, fPstar, fIMPRy, jh, jl);
4417  goto L50;
4418 L67:
4419  mnrazz(ystst, fPstst, fIMPRy, jh, jl);
4420  goto L50;
4421 // point * is not as good as jl
4422 L70:
4423  if (ystar >= fIMPRy[jh-1]) goto L73;
4424  jhold = jh;
4425  mnrazz(ystar, fPstar, fIMPRy, jh, jl);
4426  if (jhold != jh) goto L50;
4427 // calculate new point **
4428 L73:
4429  for (i = 1; i <= fNpar; ++i) {
4430  fPstst[i-1] = fP[i + jh*fMaxpar - fMaxpar-1]*.5 + fPbar[i-1]*.5;
4431  }
4432  mncalf(fPstst, ycalf);
4433  ystst = ycalf;
4434  if (ystst > fIMPRy[jh-1]) goto L30;
4435 // point ** is better than jh
4436  if (ystst < fAmin) goto L67;
4437  mnrazz(ystst, fPstst, fIMPRy, jh, jl);
4438  goto L50;
4439 // end main loop
4440 L95:
4441  if (fISW[4] >= 0) {
4442  Printf(" AN IMPROVEMENT ON THE PREVIOUS MINIMUM HAS BEEN FOUND");
4443  }
4444  reg = .1;
4445 // ask if point is new
4446 L100:
4447  mninex(fX);
4448  Eval(nparx, fGin, fAmin, fU, 4); ++fNfcn;
4449  for (i = 1; i <= fNpar; ++i) {
4450  fDirin[i-1] = reg*fIMPRdsav[i-1];
4451  if (TMath::Abs(fX[i-1] - fXt[i-1]) > fDirin[i-1]) goto L150;
4452  }
4453  goto L230;
4454 L150:
4455  fNfcnmx = fNfcnmx + npfn - fNfcn;
4456  npfn = fNfcn;
4457  mnsimp();
4458  if (fAmin >= fApsi) goto L325;
4459  for (i = 1; i <= fNpar; ++i) {
4460  fDirin[i-1] = fIMPRdsav[i-1]*.1;
4461  if (TMath::Abs(fX[i-1] - fXt[i-1]) > fDirin[i-1]) goto L250;
4462  }
4463 L230:
4464  if (fAmin < fApsi) goto L350;
4465  goto L325;
4466 /* truly new minimum */
4467 L250:
4468  fLnewmn = kTRUE;
4469  if (fISW[1] >= 1) {
4470  fISW[1] = 1;
4471  fDcovar = TMath::Max(fDcovar,.5);
4472  } else fDcovar = 1;
4473  fItaur = 0;
4474  fNfcnmx = fNfcnmx + npfn - fNfcn;
4475  fCstatu = "NEW MINIMU";
4476  if (fISW[4] >= 0) {
4477  Printf(" IMPROVE HAS FOUND A TRULY NEW MINIMUM");
4478  Printf(" *************************************");
4479  }
4480  return;
4481 // return to previous region
4482 L280:
4483  if (fISW[4] > 0) {
4484  Printf(" COVARIANCE MATRIX WAS NOT POSITIVE-DEFINITE");
4485  }
4486  goto L325;
4487 L300:
4488  fISW[0] = 1;
4489 L325:
4490  for (i = 1; i <= fNpar; ++i) {
4491  fDirin[i-1] = fIMPRdsav[i-1]*.01;
4492  fX[i-1] = fXt[i-1];
4493  }
4494  fAmin = fApsi;
4495  fEDM = sigsav;
4496 L350:
4497  mninex(fX);
4498  if (fISW[4] > 0) {
4499  Printf(" IMPROVE HAS RETURNED TO REGION OF ORIGINAL MINIMUM");
4500  }
4501  fCstatu = "UNCHANGED ";
4502  mnrset(0);
4503  if (fISW[1] < 2) goto L380;
4504  if (loop < nloop && fISW[0] < 1) goto L20;
4505 L380:
4506  if (iswtr >= 0) mnprin(5, fAmin);
4507  fItaur = 0;
4508 }
4509 
4510 ////////////////////////////////////////////////////////////////////////////////
4511 /// Transforms from internal coordinates (PINT) to external (U)
4512 ///
4513 /// The minimising routines which work in
4514 /// internal coordinates call this routine before calling FCN.
4515 
4517 {
4518  Int_t i, j;
4519 
4520  for (j = 0; j < fNpar; ++j) {
4521  i = fNexofi[j]-1;
4522  if (fNvarl[i] == 1) {
4523  fU[i] = pint[j];
4524  } else {
4525  fU[i] = fAlim[i] + (TMath::Sin(pint[j]) + 1)*.5*(fBlim[i] - fAlim[i]);
4526  }
4527  }
4528 }
4529 
4530 ////////////////////////////////////////////////////////////////////////////////
4531 /// Main initialization member function for MINUIT
4532 ///
4533 /// It initializes some constants
4534 /// (including the logical I/O unit nos.),
4535 
4537 {
4538  /* Local variables */
4539  volatile Double_t epsp1;
4540  Double_t piby2, epstry, epsbak, distnn;
4541  Int_t i, idb;
4542 
4543 // I/O unit numbers
4544  fIsysrd = i1;
4545  fIsyswr = i2;
4546  fIstkwr[0] = fIsyswr;
4547  fNstkwr = 1;
4548  fIsyssa = i3;
4549  fNstkrd = 0;
4550 // version identifier
4551  fCvrsn = "95.03++ ";
4552 // some CONSTANT
4553  fMaxint = fMaxpar;
4554  fMaxext = 2*fMaxpar;
4555  fUndefi = -54321;
4556  fBigedm = 123456;
4557  fCundef = ")UNDEFINED";
4558  fCovmes[0] = "NO ERROR MATRIX ";
4559  fCovmes[1] = "ERR MATRIX APPROXIMATE";
4560  fCovmes[2] = "ERR MATRIX NOT POS-DEF";
4561  fCovmes[3] = "ERROR MATRIX ACCURATE ";
4562 // some starting values
4563  fNblock = 0;
4564  fIcomnd = 0;
4565  fCtitl = fCundef;
4566  fCfrom = "INPUT ";
4567  fNfcn = 0;
4568  fNfcnfr = fNfcn;
4569  fCstatu = "INITIALIZE";
4570  fISW[2] = 0;
4571  fISW[3] = 0;
4572  fISW[4] = 1;
4573 // fISW[5]=0 for batch jobs, =1 for interactive jobs
4574 // =-1 for originally interactive temporarily batch
4575 
4576  fISW[5] = 0;
4577 // if (intrac(&dummy)) fISW[5] = 1;
4578 // DEBUG options set to default values
4579  for (idb = 0; idb <= 10; ++idb) { fIdbg[idb] = 0; }
4580  fLrepor = kFALSE;
4581  fLwarn = kTRUE;
4582  fLimset = kFALSE;
4583  fLnewmn = kFALSE;
4584  fIstrat = 1;
4585  fItaur = 0;
4586 // default page dimensions and 'new page' carriage control integer
4587  fNpagwd = 120;
4588  fNpagln = 56;
4589  fNewpag = 1;
4590  if (fISW[5] > 0) {
4591  fNpagwd = 80;
4592  fNpagln = 30;
4593  fNewpag = 0;
4594  }
4595  fUp = 1;
4596  fUpdflt = fUp;
4597 // determine machine accuracy epsmac
4598  epstry = .5;
4599  for (i = 1; i <= 100; ++i) {
4600  epstry *= .5;
4601  epsp1 = epstry + 1;
4602  mntiny(epsp1, epsbak);
4603  if (epsbak < epstry) goto L35;
4604  }
4605  epstry = 1e-7;
4606  fEpsmac = epstry*4;
4607  Printf(" MNINIT UNABLE TO DETERMINE ARITHMETIC PRECISION. WILL ASSUME:%g",fEpsmac);
4608 L35:
4609  fEpsmac = epstry*8;
4611 // the vlims are a non-negligible distance from pi/2
4612 // used by MNPINT to set variables "near" the physical limits
4613  piby2 = TMath::ATan(1)*2;
4614  distnn = TMath::Sqrt(fEpsma2)*8;
4615  fVlimhi = piby2 - distnn;
4616  fVlimlo = -piby2 + distnn;
4617  mncler();
4618 // Printf(" MINUIT RELEASE %s INITIALIZED. DIMENSIONS 100/50 EPSMAC=%g",(const char*)fCvrsn,fEpsmac);
4619 }
4620 
4621 ////////////////////////////////////////////////////////////////////////////////
4622 /// Interprets the SET LIM command, to reset the parameter limits
4623 ///
4624 /// Called from MNSET
4625 
4627 {
4628  /* Local variables */
4629  Double_t dxdi, snew;
4630  Int_t kint, i2, newcod, ifx=0, inu;
4631 
4632  fCfrom = "SET LIM ";
4633  fNfcnfr = fNfcn;
4634  fCstatu = "NO CHANGE ";
4635  i2 = Int_t(fWord7[0]);
4636  if (i2 > fMaxext || i2 < 0) goto L900;
4637  if (i2 > 0) goto L30;
4638 // set limits on all parameters
4639  newcod = 4;
4640  if (fWord7[1] == fWord7[2]) newcod = 1;
4641  for (inu = 1; inu <= fNu; ++inu) {
4642  if (fNvarl[inu-1] <= 0) continue;
4643  if (fNvarl[inu-1] == 1 && newcod == 1) continue;
4644  kint = fNiofex[inu-1];
4645 // see if parameter has been fixed
4646  if (kint <= 0) {
4647  if (fISW[4] >= 0) {
4648  Printf(" LIMITS NOT CHANGED FOR FIXED PARAMETER:%4d",inu);
4649  }
4650  continue;
4651  }
4652  if (newcod == 1) {
4653 // remove limits from parameter
4654  if (fISW[4] > 0) {
4655  Printf(" LIMITS REMOVED FROM PARAMETER :%3d",inu);
4656  }
4657  fCstatu = "NEW LIMITS";
4658  mndxdi(fX[kint-1], kint-1, dxdi);
4659  snew = fGstep[kint-1]*dxdi;
4660  fGstep[kint-1] = TMath::Abs(snew);
4661  fNvarl[inu-1] = 1;
4662  } else {
4663 // put limits on parameter
4664  fAlim[inu-1] = TMath::Min(fWord7[1],fWord7[2]);
4665  fBlim[inu-1] = TMath::Max(fWord7[1],fWord7[2]);
4666  if (fISW[4] > 0) {
4667  Printf(" PARAMETER %3d LIMITS SET TO %15.5g%15.5g",inu,fAlim[inu-1],fBlim[inu-1]);
4668  }
4669  fNvarl[inu-1] = 4;
4670  fCstatu = "NEW LIMITS";
4671  fGstep[kint-1] = -.1;
4672  }
4673  }
4674  goto L900;
4675 // set limits on one parameter
4676 L30:
4677  if (fNvarl[i2-1] <= 0) {
4678  Printf(" PARAMETER %3d IS NOT VARIABLE.", i2);
4679  goto L900;
4680  }
4681  kint = fNiofex[i2-1];
4682 // see if parameter was fixed
4683  if (kint == 0) {
4684  Printf(" REQUEST TO CHANGE LIMITS ON FIXED PARAMETER:%3d",i2);
4685  for (ifx = 1; ifx <= fNpfix; ++ifx) {
4686  if (i2 == fIpfix[ifx-1]) goto L92;
4687  }
4688  Printf(" MINUIT BUG IN MNLIMS. SEE F. JAMES");
4689 L92:
4690  ;
4691  }
4692  if (fWord7[1] != fWord7[2]) goto L235;
4693 // remove limits
4694  if (fNvarl[i2-1] != 1) {
4695  if (fISW[4] > 0) {
4696  Printf(" LIMITS REMOVED FROM PARAMETER %2d",i2);
4697  }
4698  fCstatu = "NEW LIMITS";
4699  if (kint <= 0) {
4700  fGsteps[ifx-1] = TMath::Abs(fGsteps[ifx-1]);
4701  } else {
4702  mndxdi(fX[kint-1], kint-1, dxdi);
4703  if (TMath::Abs(dxdi) < .01) dxdi = .01;
4704  fGstep[kint-1] = TMath::Abs(fGstep[kint-1]*dxdi);
4705  fGrd[kint-1] *= dxdi;
4706  }
4707  fNvarl[i2-1] = 1;
4708  } else {
4709  Printf(" NO LIMITS SPECIFIED. PARAMETER %3d IS ALREADY UNLIMITED. NO CHANGE.",i2);
4710  }
4711  goto L900;
4712 // put on limits
4713 L235:
4714  fAlim[i2-1] = TMath::Min(fWord7[1],fWord7[2]);
4715  fBlim[i2-1] = TMath::Max(fWord7[1],fWord7[2]);
4716  fNvarl[i2-1] = 4;
4717  if (fISW[4] > 0) {
4718  Printf(" PARAMETER %3d LIMITS SET TO %15.5g%15.5g",i2,fAlim[i2-1],fBlim[i2-1]);
4719  }
4720  fCstatu = "NEW LIMITS";
4721  if (kint <= 0) fGsteps[ifx-1] = -.1;
4722  else fGstep[kint-1] = -.1;
4723 
4724 L900:
4725  if (fCstatu != "NO CHANGE ") {
4726  mnexin(fX);
4727  mnrset(1);
4728  }
4729 }
4730 
4731 ////////////////////////////////////////////////////////////////////////////////
4732 /// Perform a line search from position START
4733 ///
4734 /// along direction STEP, where the length of vector STEP
4735 /// gives the expected position of minimum.
4736 /// - FSTART is value of function at START
4737 /// - SLOPE (if non-zero) is df/dx along STEP at START
4738 /// - TOLER is initial tolerance of minimum in direction STEP
4739 ///
4740 /// SLAMBG and ALPHA control the maximum individual steps allowed.
4741 /// The first step is always =1. The max length of second step is SLAMBG.
4742 /// The max size of subsequent steps is the maximum previous successful
4743 /// step multiplied by ALPHA + the size of most recent successful step,
4744 /// but cannot be smaller than SLAMBG.
4745 
4746 void TMinuit::mnline(Double_t *start, Double_t fstart, Double_t *step, Double_t slope, Double_t toler)
4747 {
4748  /* Local variables */
4749  Double_t xpq[12], ypq[12], slam, sdev, coeff[3], denom, flast;
4750  Double_t fvals[3], xvals[3], f1, fvmin, xvmin, ratio, f2, f3 = 0., fvmax;
4751  Double_t toler8, toler9, overal, undral, slamin, slamax, slopem;
4752  Int_t i, nparx=0, nvmax=0, nxypt, kk, ipt;
4753  Bool_t ldebug;
4754  TString cmess;
4755  char chpq[13];
4756  int l65, l70, l80;
4757 
4758  /* Function Body */
4759  l65 = 0; l70 = 0; l80 = 0;
4760  ldebug = fIdbg[1] >= 1;
4761 // starting values for overall limits on total step SLAM
4762  overal = 1e3;
4763  undral = -100;
4764 // debug check if start is ok
4765  if (ldebug) {
4766  mninex(&start[0]);
4767  Eval(nparx, fGin, f1, fU, 4); ++fNfcn;
4768  if (f1 != fstart) {
4769  Printf(" MNLINE start point not consistent, F values, parameters=");
4770  for (kk = 1; kk <= fNpar; ++kk) {
4771  Printf(" %14.5e",fX[kk-1]);
4772  }
4773  }
4774  }
4775 // set up linear search along STEP
4776  fvmin = fstart;
4777  xvmin = 0;
4778  nxypt = 1;
4779  chpq[0] = charal[0];
4780  xpq[0] = 0;
4781  ypq[0] = fstart;
4782 // SLAMIN = smallest possible value of ABS(SLAM)
4783  slamin = 0;
4784  for (i = 1; i <= fNpar; ++i) {
4785  if (step[i-1] != 0) {
4786  ratio = TMath::Abs(start[i-1] / step[i-1]);
4787  if (slamin == 0) slamin = ratio;
4788  if (ratio < slamin) slamin = ratio;
4789  }
4790  fX[i-1] = start[i-1] + step[i-1];
4791  }
4792  if (slamin == 0) slamin = fEpsmac;
4793  slamin *= fEpsma2;
4794  nparx = fNpar;
4795 
4796  mninex(fX);
4797  Eval(nparx, fGin, f1, fU, 4); ++fNfcn;
4798  ++nxypt;
4799  chpq[nxypt-1] = charal[nxypt-1];
4800  xpq[nxypt-1] = 1;
4801  ypq[nxypt-1] = f1;
4802  if (f1 < fstart) {
4803  fvmin = f1;
4804  xvmin = 1;
4805  }
4806 // quadr interp using slope GDEL and two points
4807  slam = 1;
4808  toler8 = toler;
4809  slamax = 5;
4810  flast = f1;
4811 // can iterate on two-points (cut) if no imprvmnt
4812 
4813  do {
4814  denom = (flast - fstart - slope*slam)*2 / (slam*slam);
4815  slam = 1;
4816  if (denom != 0) slam = -slope / denom;
4817  if (slam < 0) slam = slamax;
4818  if (slam > slamax) slam = slamax;
4819  if (slam < toler8) slam = toler8;
4820  if (slam < slamin) {
4821  l80 = 1;
4822  break;
4823  }
4824  if (TMath::Abs(slam - 1) < toler8 && f1 < fstart) {
4825  l70 = 1;
4826  break;
4827  }
4828  if (TMath::Abs(slam - 1) < toler8) slam = toler8 + 1;
4829  if (nxypt >= 12) {
4830  l65 = 1;
4831  break;
4832  }
4833  for (i = 1; i <= fNpar; ++i) { fX[i-1] = start[i-1] + slam*step[i-1]; }
4834  mninex(fX);
4835  nparx = fNpar;
4836  Eval(nparx, fGin, f2, fU, 4); ++fNfcn;
4837  ++nxypt;
4838  chpq[nxypt-1] = charal[nxypt-1];
4839  xpq[nxypt-1] = slam;
4840  ypq[nxypt-1] = f2;
4841  if (f2 < fvmin) {
4842  fvmin = f2;
4843  xvmin = slam;
4844  }
4845  if (fstart == fvmin) {
4846  flast = f2;
4847  toler8 = toler*slam;
4848  overal = slam - toler8;
4849  slamax = overal;
4850  }
4851  } while (fstart == fvmin);
4852 
4853  if (!l65 && !l70 && !l80) {
4854 // quadr interp using 3 points
4855  xvals[0] = xpq[0];
4856  fvals[0] = ypq[0];
4857  xvals[1] = xpq[nxypt-2];
4858  fvals[1] = ypq[nxypt-2];
4859  xvals[2] = xpq[nxypt-1];
4860  fvals[2] = ypq[nxypt-1];
4861 // begin iteration, calculate desired step
4862  do {
4863  slamax = TMath::Max(slamax,TMath::Abs(xvmin)*2);
4864  mnpfit(xvals, fvals, 3, coeff, sdev);
4865  if (coeff[2] <= 0) {
4866  slopem = coeff[2]*2*xvmin + coeff[1];
4867  if (slopem <= 0) slam = xvmin + slamax;
4868  else slam = xvmin - slamax;
4869  } else {
4870  slam = -coeff[1] / (coeff[2]*2);
4871  if (slam > xvmin + slamax) slam = xvmin + slamax;
4872  if (slam < xvmin - slamax) slam = xvmin - slamax;
4873  }
4874  if (slam > 0) {
4875  if (slam > overal)
4876  slam = overal;
4877  else if (slam < undral)
4878  slam = undral;
4879  }
4880 
4881 // come here if step was cut below
4882  do {
4883  toler9 = TMath::Max(toler8,TMath::Abs(toler8*slam));
4884  for (ipt = 1; ipt <= 3; ++ipt) {
4885  if (TMath::Abs(slam - xvals[ipt-1]) < toler9) {
4886  l70 = 1;
4887  break;
4888  }
4889  }
4890  if (l70) break;
4891 // take the step
4892  if (nxypt >= 12) {
4893  l65 = 1;
4894  break;
4895  }
4896  for (i = 1; i <= fNpar; ++i) { fX[i-1] = start[i-1] + slam*step[i-1]; }
4897  mninex(fX);
4898  Eval(nparx, fGin, f3, fU, 4); ++fNfcn;
4899  ++nxypt;
4900  chpq[nxypt-1] = charal[nxypt-1];
4901  xpq[nxypt-1] = slam;
4902  ypq[nxypt-1] = f3;
4903 // find worst previous point out of three
4904  fvmax = fvals[0];
4905  nvmax = 1;
4906  if (fvals[1] > fvmax) {
4907  fvmax = fvals[1];
4908  nvmax = 2;
4909  }
4910  if (fvals[2] > fvmax) {
4911  fvmax = fvals[2];
4912  nvmax = 3;
4913  }
4914 // if latest point worse than all three previous, cut step
4915  if (f3 >= fvmax) {
4916  if (nxypt >= 12) {
4917  l65 = 1;
4918  break;
4919  }
4920  if (slam > xvmin) overal = TMath::Min(overal,slam - toler8);
4921  if (slam < xvmin) undral = TMath::Max(undral,slam + toler8);
4922  slam = (slam + xvmin)*.5;
4923  }
4924  } while (f3 >= fvmax);
4925 
4926 // prepare another iteration, replace worst previous point
4927  if (l65 || l70) break;
4928 
4929  xvals[nvmax-1] = slam;
4930  fvals[nvmax-1] = f3;
4931  if (f3 < fvmin) {
4932  fvmin = f3;
4933  xvmin = slam;
4934  } else {
4935  if (slam > xvmin) overal = TMath::Min(overal,slam - toler8);
4936  if (slam < xvmin) undral = TMath::Max(undral,slam + toler8);
4937  }
4938  } while (nxypt < 12);
4939  }
4940 
4941 // end of iteration
4942 // stop because too many iterations
4943  if (!l70 && !l80) {
4944  cmess = " LINE SEARCH HAS EXHAUSTED THE LIMIT OF FUNCTION CALLS ";
4945  if (ldebug) {
4946  Printf(" MNLINE DEBUG: steps=");
4947  for (kk = 1; kk <= fNpar; ++kk) {
4948  Printf(" %12.4g",step[kk-1]);
4949  }
4950  }
4951  }
4952 // stop because within tolerance
4953  if (l70) cmess = " LINE SEARCH HAS ATTAINED TOLERANCE ";
4954  if (l80) cmess = " STEP SIZE AT ARITHMETICALLY ALLOWED MINIMUM";
4955 
4956  fAmin = fvmin;
4957  for (i = 1; i <= fNpar; ++i) {
4958  fDirin[i-1] = step[i-1]*xvmin;
4959  fX[i-1] = start[i-1] + fDirin[i-1];
4960  }
4961  mninex(fX);
4962  if (xvmin < 0) {
4963  mnwarn("D", "MNLINE", " LINE MINIMUM IN BACKWARDS DIRECTION");
4964  }
4965  if (fvmin == fstart) {
4966  mnwarn("D", "MNLINE", " LINE SEARCH FINDS NO IMPROVEMENT ");
4967  }
4968  if (ldebug) {
4969  Printf(" AFTER %3d POINTS,%s",nxypt,(const char*)cmess);
4970  mnplot(xpq, ypq, chpq, nxypt, fNpagwd, fNpagln);
4971  }
4972 }
4973 
4974 ////////////////////////////////////////////////////////////////////////////////
4975 /// Prints the covariance matrix v when KODE=1
4976 ///
4977 /// always prints the global correlations, and
4978 /// calculates and prints the individual correlation coefficients
4979 
4981 {
4982  /* Local variables */
4983  Int_t ndex, i, j, m, n, ncoef, nparm, id, it, ix;
4984  Int_t nsofar, ndi, ndj, iso, isw2, isw5;
4985  TString ctemp;
4986 
4987  isw2 = fISW[1];
4988  if (isw2 < 1) {
4989  Printf("%s",(const char*)fCovmes[isw2]);
4990  return;
4991  }
4992  if (fNpar == 0) {
4993  Printf(" MNMATU: NPAR=0");
4994  return;
4995  }
4996 // external error matrix
4997  if (kode == 1) {
4998  isw5 = fISW[4];
4999  fISW[4] = 2;
5000  mnemat(fP, fMaxint);
5001  if (isw2 < 3) {
5002  Printf("%s",(const char*)fCovmes[isw2]);
5003  }
5004  fISW[4] = isw5;
5005  }
5006 // correlation coeffs
5007  if (fNpar <= 1) return;
5008  mnwerr();
5009 // NCOEF is number of coeff. that fit on one line, not to exceed 20
5010  ncoef = (fNpagwd - 19) / 6;
5011  ncoef = TMath::Min(ncoef,20);
5012  nparm = TMath::Min(fNpar,ncoef);
5013  Printf(" PARAMETER CORRELATION COEFFICIENTS ");
5014  ctemp = " NO. GLOBAL";
5015  for (id = 1; id <= nparm; ++id) {
5016  ctemp += TString::Format(" %6d",fNexofi[id-1]);
5017  }
5018  Printf("%s",(const char*)ctemp);
5019  for (i = 1; i <= fNpar; ++i) {
5020  ix = fNexofi[i-1];
5021  ndi = i*(i + 1) / 2;
5022  for (j = 1; j <= fNpar; ++j) {
5023  m = TMath::Max(i,j);
5024  n = TMath::Min(i,j);
5025  ndex = m*(m-1) / 2 + n;
5026  ndj = j*(j + 1) / 2;
5027  fMATUvline[j-1] = fVhmat[ndex-1] / TMath::Sqrt(TMath::Abs(fVhmat[ndi-1]*fVhmat[ndj-1]));
5028  }
5029  nparm = TMath::Min(fNpar,ncoef);
5030  ctemp.Form(" %3d %7.5f ",ix,fGlobcc[i-1]);
5031  for (it = 1; it <= nparm; ++it) {
5032  ctemp += TString::Format(" %6.3f",fMATUvline[it-1]);
5033  }
5034  Printf("%s",(const char*)ctemp);
5035  if (i <= nparm) continue;
5036  ctemp = " ";
5037  for (iso = 1; iso <= 10; ++iso) {
5038  nsofar = nparm;
5039  nparm = TMath::Min(fNpar,nsofar + ncoef);
5040  for (it = nsofar + 1; it <= nparm; ++it) {
5041  ctemp = ctemp + TString::Format(" %6.3f",fMATUvline[it-1]);
5042  }
5043  Printf("%s",(const char*)ctemp);
5044  if (i <= nparm) break;
5045  }
5046  }
5047  if (isw2 < 3) {
5048  Printf(" %s",(const char*)fCovmes[isw2]);
5049  }
5050 }
5051 
5052 ////////////////////////////////////////////////////////////////////////////////
5053 /// Performs a local function minimization
5054 ///
5055 /// Performs a local function minimization using basically the
5056 /// method of Davidon-Fletcher-Powell as modified by Fletcher
5057 ///
5058 /// ref. -- Fletcher, Comp.J. 13,317 (1970) "switching method"
5059 
5061 {
5062  /* Local variables */
5063  Double_t gdel, gami, vlen, dsum, gssq, vsum, d;
5064  Double_t fzero, fs, ri, delgam, rhotol;
5065  Double_t gdgssq, gvg, vgi;
5066  Int_t npfn, ndex, iext, i, j, m, n, npsdf, nparx;
5067  Int_t iswtr, lined2, kk, nfcnmg, nrstrt,iter;
5068  Bool_t ldebug;
5069  Double_t toler = 0.05;
5070 
5071  if (fNpar <= 0) return;
5072  if (fAmin == fUndefi) mnamin();
5073  ldebug = kFALSE; if ( fIdbg[4] >= 1) ldebug = kTRUE;
5074  fCfrom = "MIGRAD ";
5075  fNfcnfr = fNfcn;
5076  nfcnmg = fNfcn;
5077  fCstatu = "INITIATE ";
5078  iswtr = fISW[4] - 2*fItaur;
5079  npfn = fNfcn;
5080  nparx = fNpar;
5081  vlen = (Double_t) (fNpar*(fNpar + 1) / 2);
5082  nrstrt = 0;
5083  npsdf = 0;
5084  lined2 = 0;
5085  fISW[3] = -1;
5086  rhotol = fApsi*.001;
5087  if (iswtr >= 1) {
5088  Printf(" START MIGRAD MINIMIZATION. STRATEGY %2d. CONVERGENCE WHEN EDM .LT.%9.2e",fIstrat,rhotol);
5089  }
5090 // initialization strategy
5091  if (fIstrat < 2 || fISW[1] >= 3) goto L2;
5092 // come (back) here to restart completely
5093 L1:
5094  if (nrstrt > fIstrat) {
5095  fCstatu = "FAILED ";
5096  fISW[3] = -1;
5097  goto L230;
5098  }
5099 // get full covariance and gradient
5100  mnhess();
5101  mnwerr();
5102  npsdf = 0;
5103  if (fISW[1] >= 1) goto L10;
5104 // get gradient at start point
5105 L2:
5106  mninex(fX);
5107  if (fISW[2] == 1) {
5108  Eval(nparx, fGin, fzero, fU, 2); ++fNfcn;
5109  }
5110  mnderi();
5111  if (fISW[1] >= 1) goto L10;
5112 // sometimes start with diagonal matrix
5113  for (i = 1; i <= fNpar; ++i) {
5114  fMIGRxxs[i-1] = fX[i-1];
5115  fMIGRstep[i-1] = 0;
5116  }
5117 // do line search if second derivative negative
5118  ++lined2;
5119  if (lined2 < (fIstrat + 1)*fNpar) {
5120  for (i = 1; i <= fNpar; ++i) {
5121  if (fG2[i-1] > 0) continue;
5122  if (fGrd[i-1] > 0) fMIGRstep[i-1] = -TMath::Abs(fGstep[i-1]);
5123  else fMIGRstep[i-1] = TMath::Abs(fGstep[i-1]);
5124  gdel = fMIGRstep[i-1]*fGrd[i-1];
5125  fs = fAmin;
5126  mnline(fMIGRxxs, fs, fMIGRstep, gdel, toler);
5127  mnwarn("D", "MNMIGR", "Negative G2 line search");
5128  iext = fNexofi[i-1];
5129  if (ldebug) {
5130  Printf(" Negative G2 line search, param %3d %13.3g%13.3g",iext,fs,fAmin);
5131  }
5132  goto L2;
5133  }
5134  }
5135 // make diagonal error matrix
5136  for (i = 1; i <= fNpar; ++i) {
5137  ndex = i*(i-1) / 2;
5138  for (j = 1; j <= i-1; ++j) {
5139  ++ndex;
5140  fVhmat[ndex-1] = 0;
5141  }
5142  ++ndex;
5143  if (fG2[i-1] <= 0) fG2[i-1] = 1;
5144  fVhmat[ndex-1] = 2 / fG2[i-1];
5145  }
5146  fDcovar = 1;
5147  if (ldebug) {
5148  Printf(" DEBUG MNMIGR, STARTING MATRIX DIAGONAL, VHMAT=");
5149  for (kk = 1; kk <= Int_t(vlen); ++kk) {
5150  Printf(" %10.2g",fVhmat[kk-1]);
5151  }
5152  }
5153 // ready to start first iteration
5154 L10:
5155  ++nrstrt;
5156  if (nrstrt > fIstrat + 1) {
5157  fCstatu = "FAILED ";
5158  goto L230;
5159  }
5160  fs = fAmin;
5161 // get EDM and set up loop
5162  fEDM = 0;
5163  for (i = 1; i <= fNpar; ++i) {
5164  fMIGRgs[i-1] = fGrd[i-1];
5165  fMIGRxxs[i-1] = fX[i-1];
5166  ndex = i*(i-1) / 2;
5167  for (j = 1; j <= i-1; ++j) {
5168  ++ndex;
5169  fEDM += fMIGRgs[i-1]*fVhmat[ndex-1]*fMIGRgs[j-1];
5170  }
5171  ++ndex;
5172  fEDM += fMIGRgs[i-1]*fMIGRgs[i-1]*.5*fVhmat[ndex-1];
5173  }
5174  fEDM = fEDM*.5*(fDcovar*3 + 1);
5175  if (fEDM < 0) {
5176  mnwarn("W", "MIGRAD", "STARTING MATRIX NOT POS-DEFINITE.");
5177  fISW[1] = 0;
5178  fDcovar = 1;
5179  goto L2;
5180  }
5181  if (fISW[1] == 0) fEDM = fBigedm;
5182  iter = 0;
5183  mninex(fX);
5184  mnwerr();
5185  if (iswtr >= 1) mnprin(3, fAmin);
5186  if (iswtr >= 2) mnmatu(0);
5187 // start main loop
5188 L24:
5189  if (fNfcn - npfn >= fNfcnmx) goto L190;
5190  gdel = 0;
5191  gssq = 0;
5192  for (i = 1; i <= fNpar; ++i) {
5193  ri = 0;
5194  gssq += fMIGRgs[i-1]*fMIGRgs[i-1];
5195  for (j = 1; j <= fNpar; ++j) {
5196  m = TMath::Max(i,j);
5197  n = TMath::Min(i,j);
5198  ndex = m*(m-1) / 2 + n;
5199  ri += fVhmat[ndex-1]*fMIGRgs[j-1];
5200  }
5201  fMIGRstep[i-1] = ri*-.5;
5202  gdel += fMIGRstep[i-1]*fMIGRgs[i-1];
5203  }
5204  if (gssq == 0) {
5205  mnwarn("D", "MIGRAD", " FIRST DERIVATIVES OF FCN ARE ALL ZERO");
5206  goto L300;
5207  }
5208 // if gdel positive, V not posdef
5209  if (gdel >= 0) {
5210  mnwarn("D", "MIGRAD", " NEWTON STEP NOT DESCENT.");
5211  if (npsdf == 1) goto L1;
5212  mnpsdf();
5213  npsdf = 1;
5214  goto L24;
5215  }
5216 // do line search
5217  mnline(fMIGRxxs, fs, fMIGRstep, gdel, toler);
5218  if (fAmin == fs) goto L200;
5219  fCfrom = "MIGRAD ";
5220  fNfcnfr = nfcnmg;
5221  fCstatu = "PROGRESS ";
5222 // get gradient at new point
5223  mninex(fX);
5224  if (fISW[2] == 1) {
5225  Eval(nparx, fGin, fzero, fU, 2); ++fNfcn;
5226  }
5227  mnderi();
5228 // calculate new EDM
5229  npsdf = 0;
5230 L81:
5231  fEDM = 0;
5232  gvg = 0;
5233  delgam = 0;
5234  gdgssq = 0;
5235  for (i = 1; i <= fNpar; ++i) {
5236  ri = 0;
5237  vgi = 0;
5238  for (j = 1; j <= fNpar; ++j) {
5239  m = TMath::Max(i,j);
5240  n = TMath::Min(i,j);
5241  ndex = m*(m-1) / 2 + n;
5242  vgi += fVhmat[ndex-1]*(fGrd[j-1] - fMIGRgs[j-1]);
5243  ri += fVhmat[ndex-1]*fGrd[j-1];
5244  }
5245  fMIGRvg[i-1] = vgi*.5;
5246  gami = fGrd[i-1] - fMIGRgs[i-1];
5247  gdgssq += gami*gami;
5248  gvg += gami*fMIGRvg[i-1];
5249  delgam += fDirin[i-1]*gami;
5250  fEDM += fGrd[i-1]*ri*.5;
5251  }
5252  fEDM = fEDM*.5*(fDcovar*3 + 1);
5253 // if EDM negative, not positive-definite
5254  if (fEDM < 0 || gvg <= 0) {
5255  mnwarn("D", "MIGRAD", "NOT POS-DEF. EDM OR GVG NEGATIVE.");
5256  fCstatu = "NOT POSDEF";
5257  if (npsdf == 1) goto L230;
5258  mnpsdf();
5259  npsdf = 1;
5260  goto L81;
5261  }
5262 // print information about this iteration
5263  ++iter;
5264  if (iswtr >= 3 || (iswtr == 2 && iter % 10 == 1)) {
5265  mnwerr();
5266  mnprin(3, fAmin);
5267  }
5268  if (gdgssq == 0) {
5269  mnwarn("D", "MIGRAD", "NO CHANGE IN FIRST DERIVATIVES OVER LAST STEP");
5270  }
5271  if (delgam < 0) {
5272  mnwarn("D", "MIGRAD", "FIRST DERIVATIVES INCREASING ALONG SEARCH LINE");
5273  }
5274 // update covariance matrix
5275  fCstatu = "IMPROVEMNT";
5276  if (ldebug) {
5277  Printf(" VHMAT 1 =");
5278  for (kk = 1; kk <= 10; ++kk) {
5279  Printf(" %10.2g",fVhmat[kk-1]);
5280  }
5281  }
5282  dsum = 0;
5283  vsum = 0;
5284  for (i = 1; i <= fNpar; ++i) {
5285  for (j = 1; j <= i; ++j) {
5286  if(delgam == 0 || gvg == 0) d = 0;
5287  else d = fDirin[i-1]*fDirin[j-1] / delgam - fMIGRvg[i-1]*fMIGRvg[j-1] / gvg;
5288  dsum += TMath::Abs(d);
5289  ndex = i*(i-1) / 2 + j;
5290  fVhmat[ndex-1] += d*2;
5291  vsum += TMath::Abs(fVhmat[ndex-1]);
5292  }
5293  }
5294 // smooth local fluctuations by averaging DCOVAR
5295  fDcovar = (fDcovar + dsum / vsum)*.5;
5296  if (iswtr >= 3 || ldebug) {
5297  Printf(" RELATIVE CHANGE IN COV. MATRIX=%5.1f per cent",fDcovar*100);
5298  }
5299  if (ldebug) {
5300  Printf(" VHMAT 2 =");
5301  for (kk = 1; kk <= 10; ++kk) {
5302  Printf(" %10.3g",fVhmat[kk-1]);
5303  }
5304  }
5305  if (delgam <= gvg) goto L135;
5306  for (i = 1; i <= fNpar; ++i) {
5307  fMIGRflnu[i-1] = fDirin[i-1] / delgam - fMIGRvg[i-1] / gvg;
5308  }
5309  for (i = 1; i <= fNpar; ++i) {
5310  for (j = 1; j <= i; ++j) {
5311  ndex = i*(i-1) / 2 + j;
5312  fVhmat[ndex-1] += gvg*2*fMIGRflnu[i-1]*fMIGRflnu[j-1];
5313  }
5314  }
5315 L135:
5316 // and see if converged
5317  if (fEDM < rhotol*.1) goto L300;
5318 // if not, prepare next iteration
5319  for (i = 1; i <= fNpar; ++i) {
5320  fMIGRxxs[i-1] = fX[i-1];
5321  fMIGRgs[i-1] = fGrd[i-1];
5322  }
5323  fs = fAmin;
5324  if (fISW[1] == 0 && fDcovar < .5) fISW[1] = 1;
5325  if (fISW[1] == 3 && fDcovar > .1) fISW[1] = 1;
5326  if (fISW[1] == 1 && fDcovar < .05) fISW[1] = 3;
5327  goto L24;
5328 // end main loop
5329 // call limit in MNMIGR
5330 L190:
5331  fISW[0] = 1;
5332  if (fISW[4] >= 0) {
5333  Printf(" CALL LIMIT EXCEEDED IN MIGRAD.");
5334  }
5335  fCstatu = "CALL LIMIT";
5336  goto L230;
5337 // fails to improve
5338 L200:
5339  if (iswtr >= 1) {
5340  Printf(" MIGRAD FAILS TO FIND IMPROVEMENT");
5341  }
5342  for (i = 1; i <= fNpar; ++i) { fX[i-1] = fMIGRxxs[i-1]; }
5343  if (fEDM < rhotol) goto L300;
5344  if (fEDM < TMath::Abs(fEpsma2*fAmin)) {
5345  if (iswtr >= 0) {
5346  Printf(" MACHINE ACCURACY LIMITS FURTHER IMPROVEMENT.");
5347  }
5348  goto L300;
5349  }
5350  if (fIstrat < 1) {
5351  if (fISW[4] >= 0) {
5352  Printf(" MIGRAD FAILS WITH STRATEGY=0. WILL TRY WITH STRATEGY=1.");
5353  }
5354  fIstrat = 1;
5355  }
5356  goto L1;
5357 // fails to converge
5358 L230:
5359  if (iswtr >= 0) {
5360  Printf(" MIGRAD TERMINATED WITHOUT CONVERGENCE.");
5361  }
5362  if (fISW[1] == 3) fISW[1] = 1;
5363  fISW[3] = -1;
5364  goto L400;
5365 // apparent convergence
5366 L300:
5367  if (iswtr >= 0) {
5368  Printf(" MIGRAD MINIMIZATION HAS CONVERGED.");
5369  }
5370  if (fItaur == 0) {
5371  if (fIstrat >= 2 || (fIstrat == 1 && fISW[1] < 3)) {
5372  if (fISW[4] >= 0) {
5373  Printf(" MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.");
5374  }
5375  mnhess();
5376  mnwerr();
5377  npsdf = 0;
5378  if (fEDM > rhotol) goto L10;
5379  }
5380  }
5381  fCstatu = "CONVERGED ";
5382  fISW[3] = 1;
5383 // come here in any case
5384 L400:
5385  fCfrom = "MIGRAD ";
5386  fNfcnfr = nfcnmg;
5387  mninex(fX);
5388  mnwerr();
5389  if (iswtr >= 0) mnprin(3, fAmin);
5390  if (iswtr >= 1) mnmatu(1);
5391 }
5392 
5393 ////////////////////////////////////////////////////////////////////////////////
5394 /// Performs a MINOS error analysis
5395 ///
5396 /// Performs a MINOS error analysis on those parameters for
5397 /// which it is requested on the MINOS command by calling
5398 /// MNMNOT for each parameter requested.
5399 
5401 {
5402  /* Local variables */
5403  Double_t val2mi, val2pl;
5404  Int_t nbad, ilax, ilax2, ngood, nfcnmi, iin, knt;
5405 
5406  if (fNpar <= 0) goto L700;
5407  ngood = 0;
5408  nbad = 0;
5409  nfcnmi = fNfcn;
5410 // loop over parameters requested
5411  for (knt = 1; knt <= fNpar; ++knt) {
5412  if (Int_t(fWord7[1]) == 0) {
5413  ilax = fNexofi[knt-1];
5414  } else {
5415  if (knt >= 7) break;
5416  ilax = Int_t(fWord7[knt]);
5417  if (ilax == 0) break;
5418  if (ilax > 0 && ilax <= fNu) {
5419  if (fNiofex[ilax-1] > 0) goto L565;
5420  }
5421  Printf(" PARAMETER NUMBER %3d NOT A VARIABLE. IGNORED.",ilax);
5422  continue;
5423  }
5424 L565:
5425 // calculate one pair of M E s
5426  ilax2 = 0;
5427  mnmnot(ilax, ilax2, val2pl, val2mi);
5428  if (fLnewmn) goto L650;
5429 // update NGOOD and NBAD
5430  iin = fNiofex[ilax-1];
5431  if (fErp[iin-1] > 0) ++ngood;
5432  else ++nbad;
5433  if (fErn[iin-1] < 0) ++ngood;
5434  else ++nbad;
5435  }
5436 // end of loop
5437 // printout final values
5438  fCfrom = "MINOS ";
5439  fNfcnfr = nfcnmi;
5440  fCstatu = "UNCHANGED ";
5441  if (ngood == 0 && nbad == 0) goto L700;
5442  if (ngood > 0 && nbad == 0) fCstatu = "SUCCESSFUL";
5443  if (ngood == 0 && nbad > 0) fCstatu = "FAILURE ";
5444  if (ngood > 0 && nbad > 0) fCstatu = "PROBLEMS ";
5445  if (fISW[4] >= 0) mnprin(4, fAmin);
5446  if (fISW[4] >= 2) mnmatu(0);
5447  return;
5448 // new minimum found
5449 L650:
5450  fCfrom = "MINOS ";
5451  fNfcnfr = nfcnmi;
5452  fCstatu = "NEW MINIMU";
5453  if (fISW[4] >= 0) mnprin(4, fAmin);
5454  Printf(" NEW MINIMUM FOUND. GO BACK TO MINIMIZATION STEP.");
5455  Printf(" =================================================");
5456  Printf(" V");
5457