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