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