forked from GooFit/GooFit
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathTMinuit.cc
7798 lines (7414 loc) · 258 KB
/
TMinuit.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
// This file is derived from TMinuit.cxx from ROOT v5.30/06, by Rene Brun and
// Fons Rademakers. It is distributed under the GNU Lesser General Public License,
// version 2.1. The copyright notice in that file reads:
/*************************************************************************
> * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
> * All rights reserved. *
> * *
> * For the licensing terms see $ROOTSYS/LICENSE. *
> * For the list of contributors see $ROOTSYS/README/CREDITS. *
> *************************************************************************/
//______________________________________________________________________________
//*-*-*-*-*-*-*-*-*-*-*-*The Minimization package*-*--*-*-*-*-*-*-*-*-*-*-*
//*-* ======================== *
//*-* *
//*-* This package was originally written in Fortran by Fred James *
//*-* and part of PACKLIB (patch D506) *
//*-* *
//*-* It has been converted to a C++ class by R.Brun *
//*-* The current implementation in C++ is a straightforward conversion *
//*-* of the original Fortran version: The main changes are: *
//*-* *
//*-* - The variables in the various Minuit labelled common blocks *
//*-* have been changed to the TMinuit class data members. *
//*-* - The internal arrays with a maximum dimension depending on the *
//*-* maximum number of parameters are now data members arrays with *
//*-* a dynamic dimension such that one can fit very large problems *
//*-* by simply initialising the TMinuit constructor with the maximum *
//*-* number of parameters. *
//*-* - The include file Minuit.h has been commented as much as possible*
//*-* using existing comments in the code or the printed documentation*
//*-* - The original Minuit subroutines are now member functions. *
//*-* - Constructors and destructor have been added. *
//*-* - Instead of passing the FCN function in the argument *
//*-* list, the addresses of this function is stored as pointer *
//*-* in the data members of the class. This is by far more elegant *
//*-* and flexible in an interactive environment. *
//*-* The member function SetFCN can be used to define this pointer. *
//*-* - The ROOT static function Printf is provided to replace all *
//*-* format statements and to print on currently defined output file.*
//*-* - The functions SetObjectFit(TObject *obj)/GetObjectFit() can be *
//*-* used inside the FCN function to set/get a referenced object *
//*-* instead of using global variables. *
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//BEGIN_HTML <!--
/* -->
<P>
<H2><A NAME=H2Basic-concepts-of-MINUIT.html>Basic concepts of MINUIT</A></H2>
<P>
The <A HREF="http://wwwinfo.cern.ch/asdoc/minuit/minmain.html"> MINUIT</A> package acts on a multiparameter Fortran function to which one
must give the generic name <TT>FCN</TT>. In the ROOT implementation,
the function <TT>FCN</TT> is defined via the MINUIT SetFCN member function
when an Histogram.Fit command is invoked.
The value of <TT>FCN</TT> will in general depend on one
or more variable parameters.
<P>
To take a simple example, in case of ROOT histograms (classes TH1C,TH1S,TH1F,TH1D)
the Fit function defines the Minuit fitting function as being H1FitChisquare
or H1FitLikelihood depending on the options selected.
H1FitChisquare
calculates the chisquare between the user fitting function (gaussian, polynomial,
user defined,etc) and the data for given values of the parameters.
It is the task of MINUIT to find those values of the parameters
which give the lowest value of chisquare.
<P>
<H3>Basic concepts - The transformation for parameters with limits.</H3>
<P>
For variable parameters with limits, MINUIT uses the following
transformation:
<P>
<PRE>
P = arcsin(2((P -a)/(b- a))-1) P = a+((b- a)/(2))(sinP +1)
int ext ext int
</PRE>
<P>
so that the internal value P can take on any value, while the external
int
value P can take on values only between the lower limit a and the
ext
upper limit b. Since the transformation is necessarily non-linear, it
would transform a nice linear problem into a nasty non-linear one, which
is the reason why limits should be avoided if not necessary. In addition,
the transformation does require some computer time, so it slows down the
computation a little bit, and more importantly, it introduces additional
numerical inaccuracy into the problem in addition to what is introduced in
the numerical calculation of the <TT>FCN</TT> value. The effects of
non-linearity and numerical roundoff both become more important as the
external value gets closer to one of the limits (expressed as the distance
to nearest limit divided by distance between limits). The user must
therefore be aware of the fact that, for example, if he puts limits of
(0,10^10 ) on a parameter, then the values 0.0 and 1. 0 will be
indistinguishable to the accuracy of most machines.
<P>
The transformation also affects the parameter error matrix, of course, so
MINUIT does a transformation of the error matrix (and the ``parabolic''
parameter errors) when there are parameter limits. Users should however
realize that the transformation is only a linear approximation, and that
it cannot give a meaningful result if one or more parameters is very close
to a limit, where partial Pext /partial Pint #0. Therefore, it is
recommended that:
<P>
<OL>
<LI>Limits on variable parameters should be used only when needed in order
to prevent the parameter from taking on unphysical values.
<LI>When a satisfactory minimum has been found using limits, the limits
should then be removed if possible, in order to perform or re-perform the
error analysis without limits.
</OL>
<P>
<H3>How to get the right answer from MINUIT.</H3>
<P>
MINUIT offers the user a choice of several minimization algorithms. The
MIGRAD algorithm is in general the best minimizer for
nearly all functions. It is a variable-metric method with inexact line
search, a stable metric updating scheme, and checks for
positive-definiteness. Its main weakness is that it depends heavily on
knowledge of the first derivatives, and fails miserably if they are very
inaccurate.
<P>
If parameter limits are needed, in spite of the side effects, then the
user should be aware of the following techniques to alleviate problems
caused by limits:
<P>
<H4>Getting the right minimum with limits.</H4>
<P>
If MIGRAD converges normally to a point where no parameter is near one of
its limits, then the existence of limits has probably not prevented MINUIT
from finding the right minimum. On the other hand, if one or more
parameters is near its limit at the minimum, this may be because the true
minimum is indeed at a limit, or it may be because the minimizer has
become ``blocked'' at a limit. This may normally happen only if the
parameter is so close to a limit (internal value at an odd multiple of #((pi)/(2))
that MINUIT prints a warning to this effect when it prints
the parameter values.
The minimizer can become blocked at a limit, because at a limit the
derivative seen by the minimizer partial F/partial Pint is zero no matter
what the real derivative partial F/partial Pext is.
<P>
<P>
<PRE>
((partial F)/(partial P ))= ((partial F)/(partial P ))((partial P )/(partial P )) =((partial F)/(partial P ))= 0
int ext ext int ext
</PRE>
<P>
<P>
<H4>Getting the right parameter errors with limits.</H4>
<P>
In the best case, where the minimum is far from any limits, MINUIT will
correctly transform the error matrix, and the parameter errors it reports
should be accurate and very close to those you would have got without
limits. In other cases (which should be more common, since otherwise you
wouldn't need limits), the very meaning of parameter errors becomes
problematic. Mathematically, since the limit is an absolute constraint on
the parameter, a parameter at its limit has no error, at least in one
direction. The error matrix, which can assign only symmetric errors, then
becomes essentially meaningless.
<P>
<H3>Interpretation of Parameter Errors:</H3>
<P>
There are two kinds of problems that can arise: the reliability of
MINUIT's error estimates, and their statistical interpretation, assuming
they are accurate.
<P>
<H3>Statistical interpretation:</H3>
<P>
For discussuion of basic concepts, such as the meaning of the elements of
the error matrix, or setting of exact confidence levels see:
<ol>
<li>F.James.
Determining the statistical Significance of experimental Results.
Technical Report DD/81/02 and CERN Report 81-03, CERN, 1981.</li>
<li>W.T.Eadie, D.Drijard, F.James, M.Roos, and B.Sadoulet.
Statistical Methods in Experimental Physics.
North-Holland, 1971.</li>
</ol>
<P>
<H3>Reliability of MINUIT error estimates.</H3>
<P>
MINUIT always carries around its own current estimates of the parameter
errors, which it will print out on request, no matter how accurate they
are at any given point in the execution. For example, at initialization,
these estimates are just the starting step sizes as specified by the user.
After a HESSE step, the errors are usually quite accurate,
unless there has been a problem. MINUIT, when it prints out error values,
also gives some indication of how reliable it thinks they are. For
example, those marked <TT>CURRENT GUESS ERROR</TT> are only working values
not to be believed, and <TT>APPROXIMATE ERROR</TT> means that they have
been calculated but there is reason to believe that they may not be
accurate.
<P>
If no mitigating adjective is given, then at least MINUIT believes the
errors are accurate, although there is always a small chance that MINUIT
has been fooled. Some visible signs that MINUIT may have been fooled are:
<P>
<OL>
<LI>Warning messages produced during the minimization or error analysis.
<LI>Failure to find new minimum.
<LI>Value of <TT>EDM</TT> too big (estimated Distance to Minimum).
<LI>Correlation coefficients exactly equal to zero, unless some parameters
are known to be uncorrelated with the others.
<LI>Correlation coefficients very close to one (greater than 0.99). This
indicates both an exceptionally difficult problem, and one which has been
badly parameterized so that individual errors are not very meaningful
because they are so highly correlated.
<LI>Parameter at limit. This condition, signalled by a MINUIT warning
message, may make both the function minimum and parameter errors
unreliable. See the discussion above ``Getting the right parameter errors
with limits''.
</OL>
<P>
The best way to be absolutely sure of the errors, is to use
``independent'' calculations and compare them, or compare the calculated
errors with a picture of the function. Theoretically, the covariance
matrix for a ``physical'' function must be positive-definite at the
minimum, although it may not be so for all points far away from the
minimum, even for a well-determined physical problem. Therefore, if MIGRAD
reports that it has found a non-positive-definite covariance matrix, this
may be a sign of one or more of the following:
<P>
<H5>A non-physical region:</H5>
<P>
On its way to the minimum, MIGRAD may have traversed a region which has
unphysical behaviour, which is of course not a serious problem as long as
it recovers and leaves such a region.
<P>
<H5>An underdetermined problem:</H5>
<P>
If the matrix is not positive-definite even at the minimum, this may mean
that the solution is not well-defined, for example that there are more
unknowns than there are data points, or that the parameterization of the
fit contains a linear dependence. If this is the case, then MINUIT (or any
other program) cannot solve your problem uniquely, and the error matrix
will necessarily be largely meaningless, so the user must remove the
underdeterminedness by reformulating the parameterization. MINUIT cannot
do this itself.
<P>
<H5>Numerical inaccuracies:</H5>
<P>
It is possible that the apparent lack of positive-definiteness is in fact
only due to excessive roundoff errors in numerical calculations in the
user function or not enough precision. This is unlikely in general, but
becomes more likely if the number of free parameters is very large, or if
the parameters are badly scaled (not all of the same order of magnitude),
and correlations are also large. In any case, whether the
non-positive-definiteness is real or only numerical is largely irrelevant,
since in both cases the error matrix will be unreliable and the minimum
suspicious.
<P>
<H5>An ill-posed problem:</H5>
<P>
For questions of parameter dependence, see the discussion above on
positive-definiteness.
<P>
Possible other mathematical problems are the following:
<P>
<H5>Excessive numerical roundoff:</H5>
<P>
Be especially careful of exponential and factorial functions which get big
very quickly and lose accuracy.
<P>
<H5>Starting too far from the solution:</H5>
<P>
The function may have unphysical local minima, especially at infinity in
some variables.
<H5>Minuit parameter errors in the presence of limits</H5>
This concerns the way Minuit reports the symmetric (or parabolic) errors
on parameters. It does not apply to the errors reported from Minos, which
are in general asymmetric.
<P>
The symmetric errors reported by Minuit are always calculated from
the covariance matrix, assuming that this matrix has been calculated,
usually as the result of a Migrad minimization or a direct
calculation by Hesse which inverts the second derivative matrix.
<P>
When there are no limits on the parameter in question, the error reported
by Minuit should therefore be exactly equal to the square root of the
corresponding diagonal element of the error matrix reported by Minuit.
<P>
However, when there are limits on the parameter, there is a transformation
between the internal parameter values seen by Minuit (which are unbounded)
and the external parameter values seen by the user in FCN (which remain
inside the desired limits). Therefore the internal error matrix kept by
Minuit must be transformed to an external error matrix for the user.
This is done by multiplying the (I,J)th element by DEXDIN(I)*DEXDIN(J),
where DEXDIN is the derivative of the external value with respect to the
internal value at the minimum. This is a linearization of the
transformation, and is the only way to produce an error matrix in external
coordinates meaningful to the user. But when reporting the individual
parabolic errors for limited parameters, Minuit can do a little better, so
it does. In this case, Minuit actually transforms the ends of the
internal "error bar" to external coordinates and reports the length of
this transformed interval. Strictly speaking, it is now asymmetric, but
since the origin of the asymmetry is only an artificial transformation it
does not make much sense, so the transformed errors are symmetrized.
<P>
The result of all the above is that for parameters with limits, the error
reported by Minuit is not exactly equal to the square root of the diagonal
element of the error matrix. The difference is a measure of how much the
limits deform the problem. If possible, it is suggested not to use limits
on parameters, and the problem goes away. If for some reason limits are
necessary, and you are sensitive to the difference between the two ways of
calculating the errors, it is suggested to use Minos errors which take
into account the non-linearities much more precisely.
<!--*/
// -->END_HTML
#include <stdlib.h>
#include <stdio.h>
#include <cstring>
#include <cmath>
#include "TMinuit.hh"
TMinuit* TMinuit::gMinuit = 0;
const char charal[29] = " .ABCDEFGHIJKLMNOPQRSTUVWXYZ";
// Stupid little helper method because std::string doesn't have basic functionality...
void toUpper (string& str) {
for (unsigned int i = 0; i < str.size(); ++i) {
str[i] = toupper(str[i]);
}
}
bool contains (string& toLookIn, string toLookFor) {
return (toLookIn.find(toLookFor) != string::npos);
}
#include "stdarg.h"
// Why oh why can't C++ do strings?
string local_Form (const char* va_fmt, ...) {
static char strbuffer[10000];
va_list ap;
va_start(ap, va_fmt);
vsnprintf(strbuffer, 10000, va_fmt, ap); // Use vsn, not sn, when we are passing a vararg list and not a strict vararg.
va_end(ap);
string ret(strbuffer);
return ret;
}
//______________________________________________________________________________
//______________________________________________________________________________
TMinuit::TMinuit (int maxpar) {
//*-*-*-*-*-*-*-*-*-*-*Minuit normal constructor*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ========================
//
// maxpar is the maximum number of parameters used with this TMinuit object.
fFCN = 0;
BuildArrays(maxpar);
fStatus = 0;
fEmpty = 0;
//fObjectFit = 0;
//fMethodCall = 0;
//fPlot = 0;
//fGraphicsMode = true;
SetMaxIterations();
mninit(5,6,7);
gMinuit = this;
}
//______________________________________________________________________________
TMinuit::~TMinuit()
{
//*-*-*-*-*-*-*-*-*-*-*Minuit default destructor*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* =========================
DeleteArrays();
//delete fPlot;
//delete fMethodCall;
//gROOT->GetListOfSpecials()->Remove(this);
if (gMinuit == this) gMinuit = 0;
}
//______________________________________________________________________________
void TMinuit::BuildArrays(int maxpar)
{
//*-*-*-*-*-*-*Create internal Minuit arrays for the maxpar parameters*-*-*-*
//*-* =======================================================
fMaxpar = 25;
if (maxpar >= fMaxpar) fMaxpar = maxpar+1;
fMaxpar1= fMaxpar*(fMaxpar+1);
fMaxpar2= 2*fMaxpar;
fMaxpar5= fMaxpar1/2;
fMaxcpt = 101;
fCpnam = new string[fMaxpar2];
fU = new double[fMaxpar2];
fAlim = new double[fMaxpar2];
fBlim = new double[fMaxpar2];
fPstar = new double[fMaxpar2];
fGin = new double[fMaxpar2];
fNvarl = new int[fMaxpar2];
fNiofex = new int[fMaxpar2];
fNexofi = new int[fMaxpar];
fIpfix = new int[fMaxpar];
fErp = new double[fMaxpar];
fErn = new double[fMaxpar];
fWerr = new double[fMaxpar];
fGlobcc = new double[fMaxpar];
fX = new double[fMaxpar];
fXt = new double[fMaxpar];
fDirin = new double[fMaxpar];
fXs = new double[fMaxpar];
fXts = new double[fMaxpar];
fDirins = new double[fMaxpar];
fGrd = new double[fMaxpar];
fG2 = new double[fMaxpar];
fGstep = new double[fMaxpar];
fDgrd = new double[fMaxpar];
fGrds = new double[fMaxpar];
fG2s = new double[fMaxpar];
fGsteps = new double[fMaxpar];
fPstst = new double[fMaxpar];
fPbar = new double[fMaxpar];
fPrho = new double[fMaxpar];
fWord7 = new double[fMaxpar];
fVhmat = new double[fMaxpar5];
fVthmat = new double[fMaxpar5];
fP = new double[fMaxpar1];
fXpt = new double[fMaxcpt];
fYpt = new double[fMaxcpt];
fChpt = new char[fMaxcpt+1];
// initialisation of dynamic arrays used internally in some functions
// these arrays had a fix dimension in Minuit
fCONTgcc = new double[fMaxpar];
fCONTw = new double[fMaxpar];
fFIXPyy = new double[fMaxpar];
fGRADgf = new double[fMaxpar];
fHESSyy = new double[fMaxpar];
fIMPRdsav = new double[fMaxpar];
fIMPRy = new double[fMaxpar];
fMATUvline = new double[fMaxpar];
fMIGRflnu = new double[fMaxpar];
fMIGRstep = new double[fMaxpar];
fMIGRgs = new double[fMaxpar];
fMIGRvg = new double[fMaxpar];
fMIGRxxs = new double[fMaxpar];
fMNOTxdev = new double[fMaxpar];
fMNOTw = new double[fMaxpar];
fMNOTgcc = new double[fMaxpar];
fPSDFs = new double[fMaxpar];
fSEEKxmid = new double[fMaxpar];
fSEEKxbest = new double[fMaxpar];
fSIMPy = new double[fMaxpar];
fVERTq = new double[fMaxpar];
fVERTs = new double[fMaxpar];
fVERTpp = new double[fMaxpar];
fCOMDplist = new double[fMaxpar];
fPARSplist = new double[fMaxpar];
for (int i = 0; i < fMaxpar; i++) {
fErp[i] = 0;
fErn[i] = 0;
}
}
//______________________________________________________________________________
int TMinuit::Command(const char *command)
{
// execute a Minuit command
// Equivalent to MNEXCM except that the command is given as a
// character string.
// See TMinuit::mnhelp for the full list of available commands
// See also http://wwwasdoc.web.cern.ch/wwwasdoc/minuit/node18.html for
// a complete documentation of all the available commands
//
// Returns the status of the execution:
// = 0: command executed normally
// 1: command is blank, ignored
// 2: command line unreadable, ignored
// 3: unknown command, ignored
// 4: abnormal termination (e.g., MIGRAD not converged)
// 5: command is a request to read PARAMETER definitions
// 6: 'SET INPUT' command
// 7: 'SET TITLE' command
// 8: 'SET COVAR' command
// 9: reserved
// 10: END command
// 11: EXIT or STOP command
// 12: RETURN command
//
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
int status = 0;
mncomd(command,status);
return status;
}
//______________________________________________________________________________
#ifdef NEEDS_CAREFUL_CONVERSION
TObject *TMinuit::Contour(int npoints, int pa1, int pa2)
{
// Creates a TGraph object describing the n-sigma contour of a
// TMinuit fit. The contour of the parameters pa1 and pa2 is calculated
// unsing npoints (>=4) points. The TMinuit status will be
// 0 on success and
// -1 if errors in the calling sequence (pa1, pa2 not variable)
// 1 if less than four points can be found
// 2 if npoints<4
// n>3 if only n points can be found (n < npoints)
// The status can be obtained via TMinuit::GetStatus().
//
// To get the n-sigma contour the ERRDEF parameter in Minuit has to set
// to n^2. The fcn function has to be set before the routine is called.
//
// The TGraph object is created via the interpreter. The user must cast it
// to a TGraph*. Note that the TGraph is created with npoints+1 in order to
// close the contour (setting last point equal to first point).
//
// You can find an example in $ROOTSYS/tutorials/fit/fitcont.C
if (npoints<4) {
// we need at least 4 points
fStatus= 2;
return (TObject *)0;
}
int npfound;
double *xcoor = new double[npoints+1];
double *ycoor = new double[npoints+1];
mncont(pa1,pa2,npoints,xcoor,ycoor,npfound);
if (npfound<4) {
// mncont did go wrong
Warning("Contour","Cannot find more than 4 points, no TGraph returned");
fStatus= (npfound==0 ? 1 : npfound);
delete [] xcoor;
delete [] ycoor;
return (TObject *)0;
}
if (npfound!=npoints) {
// mncont did go wrong
Warning("Contour","Returning a TGraph with %d points only",npfound);
npoints = npfound;
}
fStatus=0;
// create graph via the PluginManager
xcoor[npoints] = xcoor[0]; // add first point at end to get closed polyline
ycoor[npoints] = ycoor[0];
TObject *gr = 0;
TPluginHandler *h;
if ((h = gROOT->GetPluginManager()->FindHandler("TMinuitGraph"))) {
if (h->LoadPlugin() != -1)
gr = (TObject*)h->ExecPlugin(3,npoints+1,xcoor,ycoor);
}
delete [] xcoor;
delete [] ycoor;
return gr;
}
#endif
//______________________________________________________________________________
int TMinuit::DefineParameter( int parNo, const char *name, double initVal, double initErr, double lowerLimit, double upperLimit )
{
// Define a parameter
int err;
string sname = name;
mnparm( parNo, sname, initVal, initErr, lowerLimit, upperLimit, err);
return err;
}
//______________________________________________________________________________
void TMinuit::DeleteArrays()
{
//*-*-*-*-*-*-*-*-*-*-*-*Delete internal Minuit arrays*-*-*-*-*-*-*-*-*
//*-* =============================
if (fEmpty) return;
delete [] fCpnam;
delete [] fU;
delete [] fAlim;
delete [] fBlim;
delete [] fErp;
delete [] fErn;
delete [] fWerr;
delete [] fGlobcc;
delete [] fNvarl;
delete [] fNiofex;
delete [] fNexofi;
delete [] fX;
delete [] fXt;
delete [] fDirin;
delete [] fXs;
delete [] fXts;
delete [] fDirins;
delete [] fGrd;
delete [] fG2;
delete [] fGstep;
delete [] fGin;
delete [] fDgrd;
delete [] fGrds;
delete [] fG2s;
delete [] fGsteps;
delete [] fIpfix;
delete [] fVhmat;
delete [] fVthmat;
delete [] fP;
delete [] fPstar;
delete [] fPstst;
delete [] fPbar;
delete [] fPrho;
delete [] fWord7;
delete [] fXpt;
delete [] fYpt;
delete [] fChpt;
delete [] fCONTgcc;
delete [] fCONTw;
delete [] fFIXPyy;
delete [] fGRADgf;
delete [] fHESSyy;
delete [] fIMPRdsav;
delete [] fIMPRy;
delete [] fMATUvline;
delete [] fMIGRflnu;
delete [] fMIGRstep;
delete [] fMIGRgs;
delete [] fMIGRvg;
delete [] fMIGRxxs;
delete [] fMNOTxdev;
delete [] fMNOTw;
delete [] fMNOTgcc;
delete [] fPSDFs;
delete [] fSEEKxmid;
delete [] fSEEKxbest;
delete [] fSIMPy;
delete [] fVERTq;
delete [] fVERTs;
delete [] fVERTpp;
delete [] fCOMDplist;
delete [] fPARSplist;
fEmpty = 1;
}
//______________________________________________________________________________
int TMinuit::Eval(int npar, double *grad, double &fval, double *par, int flag)
{
// Evaluate the minimisation function
// Input parameters:
// npar: number of currently variable parameters
// par: array of (constant and variable) parameters
// flag: Indicates what is to be calculated (see example below)
// grad: array of gradients
// Output parameters:
// fval: The calculated function value.
// grad: The (optional) vector of first derivatives).
//
// The meaning of the parameters par is of course defined by the user,
// who uses the values of those parameters to calculate his function value.
// The starting values must be specified by the user.
// Later values are determined by Minuit as it searches for the minimum
// or performs whatever analysis is requested by the user.
//
// Note that this virtual function may be redefined in a class derived from TMinuit.
// The default function calls the function specified in SetFCN
//
// Example of Minimisation function:
/*
if (flag == 1) {
read input data,
calculate any necessary constants, etc.
}
if (flag == 2) {
calculate GRAD, the first derivatives of FVAL
(this is optional)
}
Always calculate the value of the function, FVAL,
which is usually a chisquare or log likelihood.
if (iflag == 3) {
will come here only after the fit is finished.
Perform any final calculations, output fitted data, etc.
}
*/
// See concrete examples in TH1::H1FitChisquare, H1FitLikelihood
if (fFCN) (*fFCN)(npar,grad,fval,par,flag);
return 0;
}
//______________________________________________________________________________
int TMinuit::FixParameter( int parNo)
{
// fix a parameter
int err;
double tmp[1];
tmp[0] = parNo+1; //set internal Minuit numbering
mnexcm( "FIX", tmp, 1, err );
return err;
}
//______________________________________________________________________________
int TMinuit::GetParameter( int parNo, double ¤tValue, double ¤tError ) const
{
// return parameter value and error
int err;
string name; // ignored
double bnd1, bnd2; // ignored
mnpout( parNo, name, currentValue, currentError, bnd1, bnd2, err );
return err;
}
//______________________________________________________________________________
int TMinuit::GetNumFixedPars() const
{
// returns the number of currently fixed parameters
return fNpfix;
}
//______________________________________________________________________________
int TMinuit::GetNumFreePars() const
{
// returns the number of currently free parameters
return fNpar;
}
//______________________________________________________________________________
int TMinuit::GetNumPars() const
{
// returns the total number of parameters that have been defined.
// (fixed and free)
return fNpar + fNpfix;
}
//______________________________________________________________________________
int TMinuit::Migrad()
{
// invokes the MIGRAD minimizer
int err;
double tmp[1];
tmp[0] = 0;
mnexcm( "MIGRAD", tmp, 0, err );
return err;
}
//______________________________________________________________________________
int TMinuit::Release( int parNo)
{
// release a parameter
int err;
double tmp[1];
tmp[0] = parNo+1; //set internal Minuit numbering
mnexcm( "RELEASE", tmp, 1, err );
return err;
}
//______________________________________________________________________________
int TMinuit::SetErrorDef( double up )
{
// To get the n-sigma contour the error def parameter "up" has to set to n^2.
int err;
mnexcm( "SET ERRDEF", &up, 1, err );
return err;
}
//______________________________________________________________________________
void TMinuit::SetFCN(void (*fcn)(int &, double *, double &f, double *, int))
{
//*-*-*-*-*-*-*To set the address of the minimization function*-*-*-*-*-*-*-*
//*-* ===============================================
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
fFCN = fcn;
}
//______________________________________________________________________________
int TMinuit::SetPrintLevel( int printLevel )
{
//set Minuit print level
// printlevel = -1 quiet (also suppresse all warnings)
// = 0 normal
// = 1 verbose
int err;
double tmp[1];
tmp[0] = printLevel;
mnexcm( "SET PRINT", tmp, 1, err );
if (printLevel <=-1) mnexcm("SET NOWarnings",tmp,0,err);
return err;
}
//______________________________________________________________________________
void TMinuit::mnamin()
{
//*-*-*-*-*-*-*-*-*-*-*-*-*Initialize AMIN*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
//*-* ===============
//*-*C Called from many places. Initializes the value of AMIN by
//*-*C calling the user function. Prints out the function value and
//*-*C parameter values if Print Flag value is high enough.
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
/* Local variables */
double fnew;
int nparx;
nparx = fNpar;
if (fISW[4] >= 1) {
printf(" FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.\n");
}
mnexin(fX);
Eval(nparx, fGin, fnew, fU, 4); ++fNfcn;
fAmin = fnew;
fEDM = fBigedm;
} /* mnamin_ */
//______________________________________________________________________________
void TMinuit::mnbins(double a1, double a2, int naa, double &bl, double &bh, int &nb, double &bwid)
{
//*-*-*-*-*-*-*-*-*-*-*Compute reasonable histogram intervals*-*-*-*-*-*-*-*-*
//*-* ======================================
//*-* Function TO DETERMINE REASONABLE HISTOGRAM INTERVALS
//*-* GIVEN ABSOLUTE UPPER AND LOWER BOUNDS A1 AND A2
//*-* AND DESIRED MAXIMUM NUMBER OF BINS NAA
//*-* PROGRAM MAKES REASONABLE BINNING FROM BL TO BH OF WIDTH BWID
//*-* F. JAMES, AUGUST, 1974 , stolen for Minuit, 1988
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
/* Local variables */
double awid,ah, al, sigfig, sigrnd, alb;
int kwid, lwid, na=0, log_;
al = min(a1,a2);
ah = max(a1,a2);
if (al == ah) ah = al + 1;
//*-*- IF NAA .EQ. -1 , PROGRAM USES BWID INPUT FROM CALLING ROUTINE
if (naa == -1) goto L150;
L10:
na = naa - 1;
if (na < 1) na = 1;
//*-*- GET NOMINAL BIN WIDTH IN EXPON FORM
L20:
awid = (ah-al) / double(na);
log_ = int(log10(awid));
if (awid <= 1) --log_;
sigfig = awid*pow(10.0, -log_);
//*-*- ROUND MANTISSA UP TO 2, 2.5, 5, OR 10
if (sigfig > 2) goto L40;
sigrnd = 2;
goto L100;
L40:
if (sigfig > 2.5) goto L50;
sigrnd = 2.5;
goto L100;
L50:
if (sigfig > 5) goto L60;
sigrnd = 5;
goto L100;
L60:
sigrnd = 1;
++log_;
L100:
bwid = sigrnd*pow(10.0, log_);
goto L200;
//*-*- GET NEW BOUNDS FROM NEW WIDTH BWID
L150:
if (bwid <= 0) goto L10;
L200:
alb = al / bwid;
lwid = int(alb);
if (alb < 0) --lwid;
bl = bwid*double(lwid);
alb = ah / bwid + 1;
kwid = int(alb);
if (alb < 0) --kwid;
bh = bwid*double(kwid);
nb = kwid - lwid;
if (naa > 5) goto L240;
if (naa == -1) return;
//*-*- REQUEST FOR ONE BIN IS DIFFICULT CASE
if (naa > 1 || nb == 1) return;
bwid *= 2;
nb = 1;
return;
L240:
if (nb << 1 != naa) return;
++na;
goto L20;
} /* mnbins_ */
//______________________________________________________________________________
void TMinuit::mncalf(double *pvec, double &ycalf)
{
//*-*-*-*-*-*-*-*-*-*Transform FCN to find further minima*-*-*-*-*-*-*-*-*-*
//*-* ====================================
//*-* Called only from MNIMPR. Transforms the function FCN
//*-* by dividing out the quadratic part in order to find further
//*-* minima. Calculates ycalf = (f-fmin)/(x-xmin)*v*(x-xmin)
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
/* Local variables */
int ndex, i, j, m, n, nparx;
double denom, f;
nparx = fNpar;
mninex(&pvec[0]);
Eval(nparx, fGin, f, fU, 4); ++fNfcn;
for (i = 1; i <= fNpar; ++i) {
fGrd[i-1] = 0;
for (j = 1; j <= fNpar; ++j) {
m = max(i,j);
n = min(i,j);
ndex = m*(m-1) / 2 + n;
fGrd[i-1] += fVthmat[ndex-1]*(fXt[j-1] - pvec[j-1]);
}
}
denom = 0;
for (i = 1; i <= fNpar; ++i) {denom += fGrd[i-1]*(fXt[i-1] - pvec[i-1]); }
if (denom <= 0) {
fDcovar = 1;
fISW[1] = 0;
denom = 1;
}
ycalf = (f - fApsi) / denom;
} /* mncalf_ */
//______________________________________________________________________________
void TMinuit::mncler()
{
//*-*-*-*-*-*-*-*-*-*-*Resets the parameter list to UNDEFINED*-*-*-*-*-*-*-*
//*-* ======================================
//*-* Called from MINUIT and by option from MNEXCM
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
int i;
fNpfix = 0;
fNu = 0;
fNpar = 0;
fNfcn = 0;
fNwrmes[0] = 0;
fNwrmes[1] = 0;
for (i = 1; i <= fMaxext; ++i) {
fU[i-1] = 0;
fCpnam[i-1] = fCundef;
fNvarl[i-1] = -1;
fNiofex[i-1] = 0;
}
mnrset(1);
fCfrom = "CLEAR ";
fNfcnfr = fNfcn;
fCstatu = "UNDEFINED ";
fLnolim = true;
fLphead = true;
} /* mncler_ */
//______________________________________________________________________________
void TMinuit::mncntr(int ike1, int ike2, int &ierrf)
{
//*-*-*-*-*Print function contours in two variables, on line printer*-*-*-*-*
//*-* =========================================================
//*-*
//*-* input arguments: parx, pary, devs, ngrid
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
static string clabel = "0123456789ABCDEFGHIJ";
/* Local variables */
double d__1, d__2;
double fcna[115], fcnb[115], contur[20];
double ylabel, fmn, fmx, xlo, ylo, xup, yup;
double devs, xsav, ysav, bwidx, bwidy, unext, ff, xb4;
int i, ngrid, ixmid, nparx, ix, nx, ny, ki1, ki2, ixzero, iy, ics;
string chmid, chln, chzero;
int ke1 = ike1+1;